From 1ebf20d6dcdd84611fd1f0f807e43468ac336682 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 13 Nov 2023 10:11:01 -0500 Subject: [PATCH 1/2] Fixup PDB reading with non-standard MD files (no element column) Also adjust "select water" to better handle bonds to extra sites (e.g., H-O-H plus extra atom for water model) Signed-off-by: Geoff Hutchison --- avogadro/io/pdbformat.cpp | 15 ++++++-- avogadro/qtplugins/select/select.cpp | 55 +++++++++++++++++++--------- 2 files changed, 50 insertions(+), 20 deletions(-) diff --git a/avogadro/io/pdbformat.cpp b/avogadro/io/pdbformat.cpp index d3f7e6aac..6f3120159 100644 --- a/avogadro/io/pdbformat.cpp +++ b/avogadro/io/pdbformat.cpp @@ -15,6 +15,7 @@ #include #include +#include #include using Avogadro::Core::Array; @@ -139,7 +140,7 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol) auto altLoc = lexicalCast(buffer.substr(16, 1), ok); string element; // Element symbol, right justified - unsigned char atomicNum; + unsigned char atomicNum = 255; if (buffer.size() >= 78) { element = buffer.substr(76, 2); element = trimmed(element); @@ -151,12 +152,20 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol) atomicNum = Elements::atomicNumberFromSymbol(element); if (atomicNum == 255) appendError("Invalid element"); - } else { + } + + if (atomicNum == 255) { // non-standard or old-school PDB file - try to parse the atom name element = trimmed(atomName); + // remove any trailing digits + while (element.size() && std::isdigit(element.back())) + element.pop_back(); + atomicNum = Elements::atomicNumberFromSymbol(element); - if (atomicNum == 255) + if (atomicNum == 255) { appendError("Invalid element"); + continue; // skip this invalid record + } } if (altLoc.compare("") && altLoc.compare("A")) { diff --git a/avogadro/qtplugins/select/select.cpp b/avogadro/qtplugins/select/select.cpp index eba27a9c3..d4987c1fe 100644 --- a/avogadro/qtplugins/select/select.cpp +++ b/avogadro/qtplugins/select/select.cpp @@ -11,11 +11,11 @@ #include #include +#include #include #include #include #include -#include #include #include @@ -181,7 +181,8 @@ void Select::selectElement(int element) for (Index i = 0; i < m_molecule->atomCount(); ++i) { if (m_molecule->atomicNumber(i) == element) { - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } else m_molecule->undoMolecule()->setAtomSelected(i, false, undoText); } @@ -223,22 +224,33 @@ void Select::selectWater() } else if (atomicNumber == 1) { // check if it's attached to a water oxygen auto bonds = m_molecule->bonds(i); - if (bonds.size() != 1 || !isWaterOxygen(bonds[0].getOtherAtom(i).index())) - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), undoText); - else - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + bool isWater = false; + // check the bonds for a water oxygen + for (auto& bond : bonds) { + if (isWaterOxygen(bond.getOtherAtom(i).index())) { + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); + isWater = true; + break; + } + } + if (!isWater) + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), + undoText); continue; } - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), + undoText); } // also select water residues (which may be isolated "O" atoms) for (const auto& residue : m_molecule->residues()) { if (residue.residueName() == "HOH") { for (auto atom : residue.residueAtoms()) { Index i = atom.index(); - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } } } @@ -258,7 +270,8 @@ void Select::selectBackboneAtoms() auto name = residue.getAtomName(atom); if (name == "CA" || name == "C" || name == "N" || name == "O") { Index i = atom.index(); - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } // also select hydrogens connected to the backbone atoms @@ -270,7 +283,8 @@ void Select::selectBackboneAtoms() if (otherName == "CA" || otherName == "C" || otherName == "N" || otherName == "O") { Index i = atom.index(); - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } } } @@ -292,7 +306,8 @@ void Select::selectSidechainAtoms() auto name = residue.getAtomName(atom); if (name != "CA" && name != "C" && name != "N" && name != "O") { Index i = atom.index(); - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } // or is it a hydrogen connected to a backbone atom? @@ -305,7 +320,8 @@ void Select::selectSidechainAtoms() if (otherName == "CA" || otherName == "C" || otherName == "N" || otherName == "O") { Index i = atom.index(); - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), + undoText); } } } @@ -357,7 +373,8 @@ void Select::enlargeSelection() Vector3 displacement = m_molecule->atomPosition3d(i) - center; Real distance = displacement.squaredNorm(); if (distance < maxDistance) { - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } } } @@ -391,9 +408,11 @@ void Select::shrinkSelection() Vector3 displacement = m_molecule->atomPosition3d(i) - center; Real distance = displacement.squaredNorm(); if (distance < maxDistance) - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); else - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(false, i), + undoText); } m_molecule->emitChanged(Molecule::Atoms); @@ -425,13 +444,15 @@ void Select::selectAtomIndex() int last = range.back().toInt(&ok2); if (ok1 && ok2) { for (int i = start; i <= last; ++i) - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } } } else { int i = item.toInt(&ok); if (ok) - m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), undoText); + m_molecule->undoMolecule()->setAtomSelected(i, evalSelect(true, i), + undoText); } } From 8b9aef5fdfb118d26c7f737b84d4275ffd424cd2 Mon Sep 17 00:00:00 2001 From: Geoff Hutchison Date: Mon, 13 Nov 2023 14:22:16 -0500 Subject: [PATCH 2/2] Fix formatting Signed-off-by: Geoff Hutchison --- avogadro/io/pdbformat.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/avogadro/io/pdbformat.cpp b/avogadro/io/pdbformat.cpp index 6f3120159..190ac45a1 100644 --- a/avogadro/io/pdbformat.cpp +++ b/avogadro/io/pdbformat.cpp @@ -14,8 +14,8 @@ #include #include -#include #include +#include #include using Avogadro::Core::Array; @@ -152,7 +152,7 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol) atomicNum = Elements::atomicNumberFromSymbol(element); if (atomicNum == 255) appendError("Invalid element"); - } + } if (atomicNum == 255) { // non-standard or old-school PDB file - try to parse the atom name @@ -224,9 +224,9 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol) else { int b = lexicalCast(buffer.substr(bCoords[i], 5), ok) - 1; if (!ok) { - appendError("Failed to parse bond connection b" + std::to_string(i) + - " " + buffer.substr(bCoords[i], 5)); - //return false; + appendError("Failed to parse bond connection b" + + std::to_string(i) + " " + buffer.substr(bCoords[i], 5)); + // return false; continue; // skip this invalid record }