Skip to content

Commit

Permalink
Fixup PDB reading with non-standard MD files (no element column)
Browse files Browse the repository at this point in the history 
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 <[email protected]>
  • Loading branch information
@ghutchis
ghutchis committed Nov 13, 2023
1 parent 2bf4c24 commit 1ebf20d
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 20 deletions.
15 changes: 12 additions & 3 deletions avogadro/io/pdbformat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@

#include <cctype>
#include <istream>
#include <iostream>
#include <string>

using Avogadro::Core::Array;
Expand Down Expand Up @@ -139,7 +140,7 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol)
auto altLoc = lexicalCast<string>(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);
Expand All @@ -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")) {
Expand Down
55 changes: 38 additions & 17 deletions avogadro/qtplugins/select/select.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,11 @@
#include <avogadro/qtgui/rwlayermanager.h>
#include <avogadro/qtgui/rwmolecule.h>

#include <QAction>
#include <QtCore/QDebug>
#include <QtCore/QRegularExpression>
#include <QtCore/QRegularExpressionMatch>
#include <QtGui/QKeySequence>
#include <QAction>
#include <QtWidgets/QInputDialog>

#include <QtCore/QStringList>
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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);
}
}
}
Expand All @@ -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
Expand All @@ -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);
}
}
}
Expand All @@ -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?
Expand All @@ -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);
}
}
}
Expand Down Expand Up @@ -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);
}
}
}
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}
}

Expand Down

0 comments on commit 1ebf20d

Please sign in to comment.