diff --git a/src/formats/mol2format.cpp b/src/formats/mol2format.cpp index 034759d26..56adb551c 100644 --- a/src/formats/mol2format.cpp +++ b/src/formats/mol2format.cpp @@ -15,6 +15,8 @@ GNU General Public License for more details. #include #include +#include +#include using namespace std; namespace OpenBabel @@ -100,6 +102,17 @@ namespace OpenBabel return(true); } + static unsigned int TotalNumberOfBonds(OBAtom* atom) + { + return atom->GetImplicitHydrogen() + atom->GetValence(); + } + static bool IsOxygenOrSulfur(OBAtom *atom) + { + switch (atom->GetAtomicNum()) { + case 8: case 16: return true; + default: return false; + } + } ///////////////////////////////////////////////////////////////// bool MOL2Format::ReadMolecule(OBBase* pOb, OBConversion* pConv) @@ -216,6 +229,7 @@ namespace OpenBabel char temp_type[BUFF_SIZE], resname[BUFF_SIZE], atmid[BUFF_SIZE]; int elemno, resnum = -1; int isotope = 0; + bool has_explicit_hydrogen = false; ttab.SetFromType("SYB"); for (i = 0;i < natoms;i++) @@ -302,6 +316,8 @@ namespace OpenBabel atom.SetAtomicNum(elemno); if (isotope) atom.SetIsotope(isotope); + else if (elemno == 1) + has_explicit_hydrogen = true; ttab.SetToType("INT"); ttab.Translate(str1,str); atom.SetType(str1); @@ -353,7 +369,8 @@ namespace OpenBabel break; } - int start,end,order; + int start, end; + bool needs_kekulization = false; for (i = 0; i < nbonds; i++) { if (!ifs.getline(buffer,BUFF_SIZE)) @@ -361,170 +378,124 @@ namespace OpenBabel sscanf(buffer,"%*d %d %d %1024s",&start,&end,temp_type); str = temp_type; - order = 1; - if (str == "ar" || str == "AR" || str == "Ar") - order = 5; + unsigned int flags = 0; + int order; + if (str == "ar" || str == "AR" || str == "Ar") { + order = 1; + flags = OB_AROMATIC_BOND; + needs_kekulization = true; + } else if (str == "AM" || str == "am" || str == "Am") order = 1; else order = atoi(str.c_str()); - mol.AddBond(start,end,order); + mol.AddBond(start, end, order, flags); } - // Make a pass to ensure that there are no double bonds - // between atoms which are also involved in aromatic bonds - // as that may ill-condition kekulization (fixes potential - // issues with molecules like CEWYIM30 (MMFF94 validation suite) - // Patch by Paolo Tosco 2012-06-07 - int idx1, idx2; - bool idx1arom, idx2arom; - FOR_BONDS_OF_MOL(bond, mol) { - if (bond->GetBO() != 2) - continue; - idx1 = bond->GetBeginAtom()->GetIdx(); - idx2 = bond->GetEndAtom()->GetIdx(); - idx1arom = idx2arom = false; - FOR_BONDS_OF_MOL(bond2, mol) { - if (&*bond == &*bond2) - continue; - if ((bond2->GetBeginAtom()->GetIdx() == idx1 || bond2->GetEndAtom()->GetIdx() == idx1) - && bond2->GetBO() == 5) - idx1arom = true; - else if ((bond2->GetBeginAtom()->GetIdx() == idx2 || bond2->GetEndAtom()->GetIdx() == idx2) - && bond2->GetBO() == 5) - idx2arom = true; - if (idx1arom && idx2arom) { - bond->SetBO(1); + // TODO: Add a test case for the statement below of Paolo Tosco + // - I am currently assuming that is not a problem for the + // the current kekulization code, but it needs to be + // checked + // "Make a pass to ensure that there are no double bonds + // between atoms which are also involved in aromatic bonds + // as that may ill-condition kekulization (fixes potential + // issues with molecules like CEWYIM30 (MMFF94 validation suite)" + + mol.SetAromaticPerceived(); // don't trigger reperception + + if (has_explicit_hydrogen) { + FOR_ATOMS_OF_MOL(atom, mol) { + unsigned int total_valence = TotalNumberOfBonds(&*atom); + switch (atom->GetAtomicNum()) { + case 8: + if (total_valence != 1) continue; + if (strcmp(atom->GetType(), "O2") != 0) continue; // TODO: the O.co2 type is lost by this point + { + OBAtomBondIter bit(&*atom); + if (!bit->IsAromatic() && bit->GetBondOrder() == 1) + atom->SetFormalCharge(-1); // set -1 charge on dangling O.co2 + } + break; + case 17: // Cl + if (total_valence == 0) + atom->SetFormalCharge(-1); break; } } } - // Now that bonds are added, make a pass to "de-aromatize" carboxylates - // and (di)thiocarboxylates - // Fixes PR#3092368 - OBAtom *carboxylCarbon, *oxysulf; - FOR_BONDS_OF_MOL(bond, mol) - { - if (bond->GetBO() != 5) - continue; - - if (bond->GetBeginAtom()->IsCarboxylOxygen() || IsThiocarboxylSulfur(bond->GetBeginAtom())) { - carboxylCarbon = bond->GetEndAtom(); - oxysulf = bond->GetBeginAtom(); + // Kekulization is neccessary if an aromatic bond is present + if (needs_kekulization) { + // "de-aromatize" carboxylates and (di)thiocarboxylates + // The typical case (in our test suite anyway) is a carboxylate binding to + // a metal ion. The two O's have charges of -0.5 and the bond orders are aromatic + FOR_ATOMS_OF_MOL(atom, mol) { + OBAtom* oxygenOrSulfur = &*atom; + // Look first for a terminal O/S + if (!IsOxygenOrSulfur(oxygenOrSulfur) || TotalNumberOfBonds(oxygenOrSulfur) != 1) continue; + OBAtomBondIter bitA(oxygenOrSulfur); + OBBond *bondA = &*bitA; + if (!bondA->IsAromatic()) continue; + // Look for the carbon + OBAtom *carbon = bondA->GetNbrAtom(oxygenOrSulfur); + if (carbon->GetAtomicNum() != 6) continue; + // Look for the other oxygen or sulfur + OBAtom* otherOxygenOrSulfur = (OBAtom*)0; + OBBond* bondB = (OBBond*)0; + FOR_BONDS_OF_ATOM(bitB, carbon) { + if (&*bitB == bondA || !bitB->IsAromatic()) continue; + OBAtom* nbr = bitB->GetNbrAtom(carbon); + if (IsOxygenOrSulfur(nbr) && TotalNumberOfBonds(nbr) == 1) { + otherOxygenOrSulfur = nbr; + bondB = &*bitB; + } } - else if (bond->GetEndAtom()->IsCarboxylOxygen() || IsThiocarboxylSulfur(bond->GetEndAtom())) { - carboxylCarbon = bond->GetBeginAtom(); - oxysulf = bond->GetEndAtom(); - } else // not a carboxylate - continue; + if (!otherOxygenOrSulfur) continue; - if (carboxylCarbon->HasDoubleBond()) { // we've already picked a double bond - bond->SetBO(1); // this should be a single bond, not "aromatic" - continue; - } + // Now set as C(=O)O + bondA->UnsetAromatic(); + oxygenOrSulfur->SetFormalCharge(-1); - // We need to choose a double bond - if (oxysulf->ExplicitHydrogenCount() == 1 || oxysulf->GetFormalCharge() == -1) { // single only - bond->SetBO(1); - continue; - } else - bond->SetBO(2); // we have to pick one, let's use this one + bondB->UnsetAromatic(); + bondB->SetBondOrder(2); } - // Make a pass to fix aromatic bond orders and formal charges - // involving nitrogen and oxygen atoms - before this patch - // the aromaticity of a molecule as simple as pyridinium - // cation could not be correctly perceived - // Patch by Paolo Tosco 2012-06-07 - OBAtom *carbon, *partner, *boundToNitrogen; - OBBitVec bv; - - bv.SetBitOn(nbonds); - bv.Clear(); - FOR_BONDS_OF_MOL(bond, mol) - { - if (bv[bond->GetIdx()] || (bond->GetBO() != 5)) - continue; - - // only bother for 6 membered rings (e.g., pyridinium) - // 5-membered rings like pyrrole, imidazole, or triazole are OK with nH - OBRing *ring = bond->FindSmallestRing(); - if ( !ring || ring->Size() != 6 ) - continue; - - if ((bond->GetBeginAtom()->IsCarbon() && bond->GetEndAtom()->IsNitrogen()) - || (bond->GetBeginAtom()->IsNitrogen() && bond->GetEndAtom()->IsCarbon())) { - carbon = (bond->GetBeginAtom()->IsCarbon() ? bond->GetBeginAtom() : bond->GetEndAtom()); - int min_n_h_bonded = 100; - int min_idx = mol.NumAtoms() + 1; - FOR_BONDS_OF_ATOM(bond2, carbon) { - if (bond2->GetBO() != 5) - continue; - partner = (bond2->GetBeginAtom() == carbon ? bond2->GetEndAtom() : bond2->GetBeginAtom()); - if (!ring->IsMember(partner)) - continue; // not in the same 6-membered ring - - if (partner->IsNitrogen() && partner->GetValence() == 3 && partner->GetFormalCharge() == 0) { - int n_h_bonded = 0; - FOR_BONDS_OF_ATOM(bond3, partner) { - boundToNitrogen = (bond3->GetBeginAtom() == partner ? bond3->GetEndAtom() : bond3->GetBeginAtom()); - if (boundToNitrogen->IsHydrogen()) - n_h_bonded++; - } - if (n_h_bonded < min_n_h_bonded || (n_h_bonded == min_n_h_bonded && partner->GetIdx() < min_idx)) { - min_n_h_bonded = n_h_bonded; - min_idx = partner->GetIdx(); - } - } + // First of all, set the atoms at the ends of the aromatic bonds to also + // be aromatic. This information is required for OBKekulize. + FOR_BONDS_OF_MOL(bond, mol) { + if (bond->IsAromatic()) { + bond->GetBeginAtom()->SetAromatic(); + bond->GetEndAtom()->SetAromatic(); } - FOR_BONDS_OF_ATOM(bond2, carbon) { - if (bond2->GetBO() != 5) - continue; - partner = (bond2->GetBeginAtom() == carbon ? bond2->GetEndAtom() : bond2->GetBeginAtom()); - if (partner->IsNitrogen() && partner->GetValence() == 3 && partner->GetFormalCharge() == 0) { - int n_ar_bond = 0; - FOR_BONDS_OF_ATOM(bond3, partner) { - boundToNitrogen = (bond3->GetBeginAtom() == partner ? bond3->GetEndAtom() : bond3->GetBeginAtom()); - if (boundToNitrogen->IsOxygen() && boundToNitrogen->GetValence() == 1) { - n_ar_bond = -1; - break; - } - if (bond3->GetBO() == 5) - ++n_ar_bond; - } - if (n_ar_bond == -1) - continue; - if (partner->GetIdx() == min_idx) { - partner->SetFormalCharge(1); - if (n_ar_bond == 1) { - bond2->SetBO(2); - } - } - else if (n_ar_bond == 1) { - bond2->SetBO(1); - } - } - bv.SetBitOn(bond2->GetIdx()); - } - } else if ((bond->GetBeginAtom()->IsCarbon() && bond->GetEndAtom()->IsOxygen()) - || (bond->GetBeginAtom()->IsOxygen() && bond->GetEndAtom()->IsCarbon())) { - OBAtom *atom1, *atom2; - atom1 = bond->GetBeginAtom(); - atom2 = bond->GetEndAtom(); - // set formal charges for pyrilium - // (i.e., this bond is a 6-membered ring, aromatic, and C-O) - if (atom1->IsOxygen() && atom1->IsInRingSize(6)) - atom1->SetFormalCharge(1); - else if (atom2->IsOxygen() && atom2->IsInRingSize(6)) - atom2->SetFormalCharge(1); + } + bool ok = OBKekulize(&mol); + if (!ok) { + stringstream errorMsg; + errorMsg << "Failed to kekulize aromatic bonds in MOL2 file"; + std::string title = mol.GetTitle(); + if (!title.empty()) + errorMsg << " (title is " << title << ")"; + errorMsg << endl; + obErrorLog.ThrowError(__FUNCTION__, errorMsg.str(), obWarning); + // return false; Should we return false for a kekulization failure? } } + + mol.EndModify(); + // Suggestion by Liu Zhiguo 2008-01-26 // Mol2 files define atom types -- there is no need to re-perceive mol.SetAtomTypesPerceived(); - mol.EndModify(); + + if (!has_explicit_hydrogen) { + // Guess how many hydrogens are present on each atom based on typical valencies + // TODO: implement the MOL2 valence model (if it exists) + FOR_ATOMS_OF_MOL(matom, mol) { + if (matom->GetImplicitHydrogen() == 0) + OBAtomAssignTypicalImplicitHydrogens(&*matom); + } + } //must add generic data after end modify - otherwise it will be blown away if (comment) diff --git a/src/obfunctions.cpp b/src/obfunctions.cpp index db2c7cc7e..9bb211eeb 100644 --- a/src/obfunctions.cpp +++ b/src/obfunctions.cpp @@ -111,7 +111,8 @@ namespace OpenBabel case 0: switch (bosum) { case 0: case 1: case 2: case 3: return 3; - case 4: case 5: return 5; + case 4: return 4; // don't round up to 5 for nitrogen + case 5: return 5; } break; case 1: diff --git a/test/files/gzip.in b/test/files/gzip.in index de3532f6c..bf5d0f9bc 100644 --- a/test/files/gzip.in +++ b/test/files/gzip.in @@ -19,10 +19,10 @@ CCC(C(C(=O)O)O)(O)C OC1C(OP(=O)(O)O)C(O)C(C(C1O)O)O O=CN(c1ccc(cc1)C(=O)NC(C(=O)O)CCC(=O)O)CC1CNc2c(N1)c(=O)nc([nH]2)N #c4.mol2.gz -O=C/C=C/c1ccc(o1)[N](=O)O -COC(=O)/C=C/c1ccc(o1)[N](=O)O -C[NH+](Cc1ccccc1)C.Cl -CCCCCCCC[C@H]1O[C@@H]1CCCCCCCC(=O)O.CCCCCCCCCC[C@H]1O[C@@H]1CCCCCC(=O)O.[Mg] +O=C/C=C/c1ccc(o1)[N](=O)[O-] +COC(=O)/C=C/c1ccc(o1)[N](=O)[O-] +C[NH+](Cc1ccccc1)C.[Cl-] +CCCCCCCC[C@H]1O[C@@H]1CCCCCCCC(=O)[O-].CCCCCCCCCC[C@H]1O[C@@H]1CCCCCC(=O)[O-].[Mg] #c5.smi.gz COc1ccc2c(c1)c1c([nH]2)ccc2c1c[n+](CCN1CCC(CC1)C1CCN(CC1)CC[n+]1ccc3c(c1)c1c(cc3)[nH]c3c1cc(OC)cc3)cc2 O=C(N(C)C)Nc1ccc(c(c1)Cl)Cl diff --git a/test/gziptest.cpp b/test/gziptest.cpp index 61cea1dd6..c1b4a9a51 100644 --- a/test/gziptest.cpp +++ b/test/gziptest.cpp @@ -45,12 +45,22 @@ void checkResults(const string& file, const vector& correctResults) while(conv.Read(&mol)) { OBConversion strconv; - OBMol mol2; - string gzippedstr = conv.WriteString(&mol); - strconv.SetInAndOutFormats(conv.GetInFormat(), conv.FindFormat("CAN"), conv.GetInGzipped()); - strconv.SetOptions("n",OBConversion::OUTOPTIONS); - strconv.ReadString(&mol2, gzippedstr); - string cansmi = strconv.WriteString(&mol2, true); + string cansmi; + if (conv.GetInFormat() == conv.FindFormat("mol2")) { + // Don't test roundtripping for mol2 as the atom types are changed by the reader (TODO: Fix this) + strconv.SetOutFormat("can"); + strconv.SetOptions("n", conv.OUTOPTIONS); + cansmi = strconv.WriteString(&mol, true); + } + else { + // Testing roundtrip through gz with Read/WriteString + OBMol mol2; + string gzippedstr = conv.WriteString(&mol); + strconv.SetInAndOutFormats(conv.GetInFormat(), conv.FindFormat("CAN"), conv.GetInGzipped()); + strconv.SetOptions("n", OBConversion::OUTOPTIONS); + strconv.ReadString(&mol2, gzippedstr); + cansmi = strconv.WriteString(&mol2, true); + } results.push_back(cansmi); } @@ -120,6 +130,11 @@ int gziptest(int argc, char* argv[]) string filepath; vector correctResults; //read from file + OBConversion conv; + conv.SetInAndOutFormats("smi", "can"); + conv.SetOptions("n", conv.OUTOPTIONS); + OBMol mol; + while(getline(ifs, line)) { if(line.length() == 0) @@ -148,7 +163,9 @@ int gziptest(int argc, char* argv[]) } else { - correctResults.push_back(line); + conv.ReadString(&mol, line); + std::string can = conv.WriteString(&mol, true); + correctResults.push_back(can); } } checkResults(filepath, correctResults);