C++ OBMol::AddHydrogens方法代码示例

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

if(argc<2){
cout << "Usage: ProgrameName InputFileName";
return 1;
}

ifstream ifs(argv[1]);
if(!ifs){
cout << "Cannot open input file";
return 1;
}

OpenBabel::OBConversion conv(&ifs, &cout);
if(conv.SetInAndOutFormats("mol","sdf")){
  OpenBabel::OBMol mol;
  if(conv.Read(&mol)){
    cout << "Mol.wt: "<<mol.GetMolWt()<<endl;//works even with implicit hydrogen
    mol.AddHydrogens(false,false); //ensure that hydrogens are all explicit
    cout << "No. of atoms: " << mol.NumAtoms()<<endl;
    cout << "No. of hvy atoms: " << mol.NumHvyAtoms() << endl;
    cout << "No. of bonds: " << mol.NumBonds() << endl;
    double exactMass = 0.0;
    FOR_ATOMS_OF_MOL(a, mol){
      cout << "iterating..."<<endl;
      exactMass += a->GetExactMass();
    }

void Molecule::addHydrogens(Atom *a,
                            const QList<unsigned long> &atomIds,
                            const QList<unsigned long> &bondIds)
{
    if (atomIds.size() != bondIds.size()) {
        qDebug() << "Error, addHydrogens called with atom & bond id lists of different size!";
    }

    // Construct an OBMol, call AddHydrogens and translate the changes
    OpenBabel::OBMol obmol = OBMol();
    if (a) {
        OpenBabel::OBAtom *obatom = obmol.GetAtom(a->index()+1);
        // Set implicit valence for unusual elements not handled by OpenBabel
        // PR#2803076
        switch (obatom->GetAtomicNum()) {
        case 3:
        case 11:
        case 19:
        case 37:
        case 55:
        case 85:
        case 87:
            obatom->SetImplicitValence(1);
            obatom->SetHyb(1);
            obmol.SetImplicitValencePerceived();
            break;

        case 4:
        case 12:
        case 20:
        case 38:
        case 56:
        case 88:
            obatom->SetImplicitValence(2);
            obatom->SetHyb(2);
            obmol.SetImplicitValencePerceived();
            break;

        case 84: // Po
            obatom->SetImplicitValence(2);
            obatom->SetHyb(3);
            obmol.SetImplicitValencePerceived();
            break;

        default: // do nothing
            break;
        }
        obmol.AddHydrogens(obatom);
    }
    else
        obmol.AddHydrogens();
    // All new atoms in the OBMol must be the additional hydrogens
    unsigned int numberAtoms = numAtoms();
    int j = 0;
    for (unsigned int i = numberAtoms+1; i <= obmol.NumAtoms(); ++i, ++j) {
        if (obmol.GetAtom(i)->IsHydrogen()) {
            OpenBabel::OBAtom *obatom = obmol.GetAtom(i);
            Atom *atom;
            if (atomIds.isEmpty())
                atom = addAtom();
            else if (j < atomIds.size())
                atom = addAtom(atomIds.at(j));
            else {
                qDebug() << "Error - not enough unique ids in addHydrogens.";
                break;
            }
            atom->setOBAtom(obatom);
            // Get the neighbor atom
            OpenBabel::OBBondIterator iter;
            OpenBabel::OBAtom *next = obatom->BeginNbrAtom(iter);
            Bond *bond;
            if (bondIds.isEmpty())
                bond = addBond();
            else // Already confirmed by atom ids
                bond = addBond(bondIds.at(j));
            bond->setEnd(Molecule::atom(atom->index()));
            bond->setBegin(Molecule::atom(next->GetIdx()-1));
        }
    }
    for (unsigned int i = 1; i <= numberAtoms; ++i) {
        // Warning -- OB atom index off-by-one here
        atom(i-1)->setPartialCharge(obmol.GetAtom(i)->GetPartialCharge());
    }
}

//.........这里部分代码省略.........
			if (!rebuild_mol && keep3DcoordsProtonate)
			{
				for(OpenBabel::OBAtomIterator a_it=obmol->BeginAtoms(); a_it!=obmol->EndAtoms(); ++a_it)
				{
					coord_backup.insert(make_pair((*a_it)->GetIdx(), (*a_it)->GetVector()));
				}
			}
			
			
			// Create 3D coordinates for individual fragments using template library
			OpenBabel::OBBuilder builder;
			builder.Build(*obmol);
			obmol->SetDimension(3);
			
			
			// Set up a constraints object for OpenBabel minimization
			// It is used to set positional constraints on non-hydrogen atoms when only hydrogens should be minimized
			OpenBabel::OBFFConstraints constraints;
			
			
			// Reset original coordinates and set positional constraints on non-hydrogen atoms
			if (!rebuild_mol && keep3DcoordsProtonate)
			{
				for(OpenBabel::OBAtomIterator a_it=obmol->BeginAtoms(); a_it!=obmol->EndAtoms(); ++a_it)
				{
					(*a_it)->SetVector(coord_backup[(*a_it)->GetIdx()]);
					
					constraints.AddAtomConstraint((*a_it)->GetIdx());
				}
			}
			
			
			// Now add hydrogens
			obmol->AddHydrogens(false, true, pH);
			
			
			// Setup MMFF94 force-field
			setup_ok = pFF->Setup(*obmol, constraints);
			
			
			if (setup_ok)
			{
				// Optimize bond-lengths and torsions
				pFF->SteepestDescent(70, 1.0e-4); // 250, 1e-04
				
				// Rotate fragments around rotatable bonds
				pFF->WeightedRotorSearch(10, 5); // 200, 25
				
				// Optimize bond-lengths and torsions again
				pFF->ConjugateGradients(70, 1.0e-4); // 250, 1e-06
				pFF->UpdateCoordinates(*obmol);
			}
			else
			{
				Log.level(20) << "[Warning:] Openbabel force-field setup failed for molecule " << i << ". Will only use template-based approach without force-field optimization for this molecule!" << endl;
			}
		}
		catch (...)
		{
			Log.level(20) << "Error while trying to generate 3D coordinates for molecule " << i<< " ";
			if (mol->getName() != "")
			{
				Log.level(20) << "'" << mol->getName() << "'";
			}
			Log.level(20) << ". Will write unmodified input coordinates to output-file." << endl;