Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixup PDB reading with non-standard MD files (no element column) #1450

Merged
merged 2 commits into from
Nov 13, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
21 changes: 15 additions & 6 deletions avogadro/io/pdbformat.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <avogadro/core/vector.h>

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

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 Expand Up @@ -215,9 +224,9 @@ bool PdbFormat::read(std::istream& in, Core::Molecule& mol)
else {
int b = lexicalCast<int>(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
}

Expand Down
55 changes: 38 additions & 17 deletions avogadro/qtplugins/select/select.cpp
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