C++ PDBFile::close方法代码示例

int main(int argc, char* argv[])
{
	FDPB fdpb;

	// instantiate CommandlineParser object
	CommandlineParser parpars("CalculateSolvationFreeEnergy", "calculate solvation free energy of a protein using AMBER ", VERSION, String(__DATE__), "ForceFields");
	parpars.registerParameter("pdb",  "input pdb file ", INFILE,  true);
	parpars.registerParameter("epsilon_medium", "dielectric constant in medium", DOUBLE, false,  fdpb.options.getReal(FDPB::Option::SOLVENT_DC));
	parpars.registerParameter("epsilon_vacuum", "dielectric constant in vacuum", DOUBLE, false, 1);

	// the manual
	String man = String("This tool computes the solvation free energy of a pdb file using the Jackson-Sternberg approach (bonded energy using a force field and a non bonded energy (electrostatics only) by solving the Poisson-Boltzmann equation. Parameters are the dielectric constants for the medium (-epsilon_medium) and the vacuum (-epsilon_vacuum).");

	parpars.setToolManual(man);

	parpars.setSupportedFormats("pdb", "pdb");

	// parse the command line
	parpars.parse(argc, argv);


	PDBFile pdb;
	pdb.open(parpars.get("pdb"), std::ios::in);

	if (!pdb)
	{
		// if file does not exist: complain and abort
		Log.error() << "error opening " << parpars.get("pdb") << " for input." << std::endl;
		exit(2);
	}

	System sys;
	pdb >> sys;
	pdb.close();

	// normalize the names and build all bonds
	FragmentDB db("");
	sys.apply(db.normalize_names);
	sys.apply(db.build_bonds);

	// TODO: Ask ResidueChecker if everything is ok! see tool CalculateEnergy

	// create an AMBER force field without non-bonded interactions
	AmberFF FF(sys);

	// calculate the total energy
	float total_energy = FF.updateEnergy();
	//Log << FF.getResults() << std::endl;
	//Log << "   total energy using force field evaluation: " << total_energy << " kJ/mol" << std::endl;

	//Log << "removing non bonded energy terms ..." << std::endl;
	FF.removeComponent("Amber NonBonded");

	// calculate the internal energy (neglecting internal VdW interactions)
	float internal_energy = FF.updateEnergy();
	//Log << FF.getResults() << std::endl;
	Log << "  internal energy: " << internal_energy << " kJ/mol" << std::endl;

	// assign atom radii
	AssignRadiusProcessor radius_processor("radii/PARSE.siz");
	sys.apply(radius_processor);

	// calculate the electrostatic part of the solvation energy	
	//FDPB fdpb;

	float dielectric_const = fdpb.options.getReal(FDPB::Option::SOLVENT_DC);
	if (parpars.has("epsilon_medium"))
		dielectric_const = parpars.get("epsilon_medium").toFloat();
	fdpb.options[FDPB::Option::SOLVENT_DC] = dielectric_const;
	Log << "... using dielectric constant in medium: " << fdpb.options[FDPB::Option::SOLVENT_DC].toFloat() << std::endl;

	fdpb.setup(sys);
	fdpb.solve();

	float solvent_energy = fdpb.getEnergy();

	dielectric_const = 1.0;
	if (parpars.has("epsilon_vacuum"))
		dielectric_const = parpars.get("epsilon_vacuum").toFloat();
	fdpb.options[FDPB::Option::SOLVENT_DC] = dielectric_const;
	Log << "... using dielectric constant in vacuum: " << fdpb.options[FDPB::Option::SOLVENT_DC].toFloat() << std::endl;

	fdpb.setup(sys);
	fdpb.solve();

	float vacuum_energy = fdpb.getEnergy();
	Log << "\n  electrostatic solvation free energy: "	<< solvent_energy - vacuum_energy << std::endl;

	Log << "\n  total energy using a combination of force field and FDPB evaluation: "
		  << internal_energy - vacuum_energy + solvent_energy << " kJ/mol" << std::endl;

	return 0;
}

int main(int argc, char** argv)
{

    if (argc != 3)
    {
        Log << "Usage:" << argv[0] << " <PDB infile> <PDB outfile> [<amber parameter file>]" << endl;

        return 1;
    }

    System system;

    PDBFile f;
    f.open(argv[1]);
    if (f.bad())
    {
        Log.error() << "cannot read PDB file " << argv[1] << endl;
        return 2;
    }
    f >> system;
    f.close();

    FragmentDB db("");

    // ResidueChecker check(db);
    // system.apply(check);


    db.normalize_names.setNamingStandard("Amber");
    system.apply(db.normalize_names);

    system.apply(db.build_bonds);

    Size cyx_counter = 0;
    Size hip_counter = 0;

    ResidueIterator it = system.beginResidue();
    for (; +it; ++it)
    {
        if (it->getFullName() == "CYS-S")
        {
            it->setName("CYX");
            cyx_counter++;
        }
        if (it->getFullName() == "HIS")
        {
            it->setName("HIP");
            hip_counter++;
        }
    }

    if (cyx_counter > 0)
    {
        Log.info() << "Renamed " << cyx_counter << " residues from CYS-S to CYX"
                   << endl;
    }

    if (hip_counter > 0)
    {
        Log.info() << "Renamed " << hip_counter << " residues from HIS to HIP"
                   << endl;
    }

    PDBFile g;
    g.open(argv[2], ios::out);
    if (g.bad())
    {
        Log.error() << "cannot write PDB file " << argv[2] << endl;
        return 2;
    }


    g << system;
    g.close();

    Log.info()
            << endl
            << "Conversion to AMBER naming scheme done. Please note that you might"
            << endl
            << "have to edit the resulting file by hand (if there are HIS or CYS"
            << endl
            << "residues in the original file, e. g.)"
            << endl
            << endl
            << "Good luck."
            << endl;

}

int main(int argc, char* argv[])
{
	// instantiate CommandlineParser object
	CommandlineParser parpars("CalculateEnergy", "calculate free energy of a protein ", VERSION, String(__DATE__), "ForceFields");
	parpars.registerMandatoryInputFile("pdb", "input pdb file ");

	// TODO: offer upload of a distinguished fragDB file?

	// choice of force field
	parpars.registerOptionalStringParameter("force_field", "force field (AMBER, MMFF94)", "AMBER");
	list<String> force_fields;
	force_fields.push_back("AMBER");
	force_fields.push_back("MMFF94");
	// TODO: shall we offer CHARM as well?
	parpars.setParameterRestrictions("force_field", force_fields);

	// TODO: shall we offer a force field parameter file upload?

	// TODO: check the naming!
	parpars.registerOptionalDoubleParameter("non_bond_cutoff",  "cutoff radius in calculations of nonbonded interactions", 20.0);
	parpars.registerOptionalDoubleParameter("elec_stat_cuton", "electrostatic cuton", 13.0);
	parpars.registerOptionalDoubleParameter("elec_stat_cutoff", "electrostatic cutoff", 15.0);
	parpars.registerFlag("dist_dep_dielec", "apply distance dependent dielectric constant", false);
	// NOTE: assign is the default
	//parpars.registerFlag("assign_typenames","automatically assign type names to the system", false);
	//parpars.registerFlag("assign_types","automatically assign types to the system if missing", false);
	//parpars.registerFlag("assign_charges", "automatically assign charges to the system if missing", false);
	// TODO: if we only allow PDBFile to upload then from where do we get stuff to overwrite?? 
	parpars.registerFlag("overwrite_types", "overwrite even non-empty type names", false);
	parpars.registerFlag("overwrite_charges","overwrite even non-zero charges", false);

	// the manual
	String man = String("This tool computes the free energy of a pdb file using a specified force field (-force_field) and force field parameters (-non_bond_cutoff, -elec_stat_cuton ... ).");

	parpars.setToolManual(man);

	parpars.setSupportedFormats("pdb", "pdb");

	// parse the command line
	parpars.parse(argc, argv);


	PDBFile pdb;
	pdb.open(parpars.get("pdb"), std::ios::in);

	if (!pdb)
	{
		// if file does not exist: complain and abort
		Log.error() << "error opening " << parpars.get("pdb") << " for input." << std::endl;
		exit(2);
	}

	System sys;
	pdb >> sys;
	pdb.close();

	// normalize the names and build all bonds
	FragmentDB db("");
	sys.apply(db.normalize_names);
	sys.apply(db.build_bonds);

	// check the structure
	Log.info() << " checking residues..." << std::endl;
	ResidueChecker rc(db);
	sys.apply(rc);

	if (!rc.getStatus())
	{
		Log.error() << "There are errors in the given structure. Use the ResidueChecker tool for further investigation." << std::endl;
		exit(2);
	}

	// create a force field	
	AmberFF* amber_force_field = NULL;
	MMFF94*  mmff_force_field = NULL;
	ForceField* force_field = NULL;

	if (parpars.has("force_field"))
	{
		String penalty_table = parpars.get("force_field");
		if (penalty_table == "AMBER")
		{
			amber_force_field = new AmberFF();
			force_field = amber_force_field;
			Log << " using the amber force field" << std::endl;
		}
		else if (penalty_table == "MMFF94")
		{
			mmff_force_field = new MMFF94();
			force_field = mmff_force_field;
			Log << " using the MMFF94 force field" << std::endl;
		}
		else
		{
			Log.error() << "Unknown force field " << parpars.get("force_field") << " Abort." << std::endl;
			exit(2);
		}
	}

	if (!amber_force_field && !mmff_force_field)
//.........这里部分代码省略.........