Skip to content

Commit

Permalink
Merge pull request OpenChemistry#1848 from ghutchis/property-radical-…
Browse files Browse the repository at this point in the history 
…formula

Read radicals from SDF / Molfile
  • Loading branch information
@ghutchis
ghutchis authored Dec 7, 2024
2 parents a5bfe59 + f21926a commit 72fb7c9
Show file tree
Hide file tree
Showing 2 changed files with 99 additions and 8 deletions.
100 changes: 92 additions & 8 deletions avogadro/io/mdlformat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,9 @@ namespace Avogadro::Io {

using chargePair = std::pair<size_t, int>;
namespace {

// Helper function to handle partial charge property blocks
// e.g. PUBCHEM_MMFF94_PARTIAL_CHARGES
void handlePartialCharges(Core::Molecule& mol, std::string data)
{
// the string starts with the number of charges
Expand All @@ -59,6 +62,7 @@ void handlePartialCharges(Core::Molecule& mol, std::string data)
bool MdlFormat::read(std::istream& in, Core::Molecule& mol)
{
string buffer;
int spinMultiplicity = 1; // check for RAD lines

// The first line is the molecule name.
getline(in, buffer);
Expand Down Expand Up @@ -124,6 +128,9 @@ bool MdlFormat::read(std::istream& in, Core::Molecule& mol)
Atom newAtom = mol.addAtom(atomicNum);
newAtom.setPosition3d(pos);
// In case there's no CHG property
if (charge == 4) // doublet radical
spinMultiplicity += 1;

charge = (charge > 4) ? ((charge <= 7) ? 4 - charge : 0)
: ((charge < 4) ? charge : 0);
if (charge)
Expand Down Expand Up @@ -191,6 +198,49 @@ bool MdlFormat::read(std::istream& in, Core::Molecule& mol)
if (charge)
chargeList.emplace_back(index, charge);
}
} else if (prefix == "M RAD") {
// radical center
spinMultiplicity = 1; // reset and count
size_t entryCount(lexicalCast<int>(buffer.substr(6, 3), ok));
for (size_t i = 0; i < entryCount; i++) {
size_t index(lexicalCast<size_t>(buffer.substr(10 + 8 * i, 3), ok) - 1);
if (!ok) {
appendError("Error parsing radical atom index:" +
buffer.substr(10 + 8 * i, 3));
return false;
}
auto radical(lexicalCast<int>(buffer.substr(14 + 8 * i, 3), ok));
if (!ok) {
appendError("Error parsing radical type:" +
buffer.substr(14 + 8 * i, 3));
return false;
}
// we don't set radical centers, just count them
// for the total spin multiplicity
if (radical == 2) // doublet
spinMultiplicity += 1;
else if (radical == 3) // triplet
spinMultiplicity += 2;
}
} else if (prefix == "M ISO") {
// isotope
size_t entryCount(lexicalCast<int>(buffer.substr(6, 3), ok));
for (size_t i = 0; i < entryCount; i++) {
size_t index(lexicalCast<size_t>(buffer.substr(10 + 8 * i, 3), ok) - 1);
if (!ok) {
appendError("Error parsing isotope atom index:" +
buffer.substr(10 + 8 * i, 3));
return false;
}
auto isotope(lexicalCast<int>(buffer.substr(14 + 8 * i, 3), ok));
if (!ok) {
appendError("Error parsing isotope type:" +
buffer.substr(14 + 8 * i, 3));
return false;
}
// TODO: Implement isotope setting
// mol.atom(index).setIsotope(isotope);
}
}
}

Expand All @@ -206,6 +256,10 @@ bool MdlFormat::read(std::istream& in, Core::Molecule& mol)
mol.setFormalCharge(index, charge);
}

// Set the total spin multiplicity
if (spinMultiplicity > 1)
mol.setData("totalSpinMultiplicity", spinMultiplicity);

// Check that all atoms were handled.
if (mol.atomCount() != static_cast<size_t>(numAtoms) ||
mol.bondCount() != static_cast<size_t>(numBonds)) {
Expand Down Expand Up @@ -253,6 +307,7 @@ bool MdlFormat::read(std::istream& in, Core::Molecule& mol)
bool MdlFormat::readV3000(std::istream& in, Core::Molecule& mol)
{
string buffer;
int spinMultiplicity = 1; // check for RAD lines
// we should have M V30 BEGIN CTAB
getline(in, buffer);
if (trimmed(buffer) != "M V30 BEGIN CTAB") {
Expand Down Expand Up @@ -320,15 +375,40 @@ bool MdlFormat::readV3000(std::istream& in, Core::Molecule& mol)
// check for formal charge in the atom block
// CHG=1 for example
if (atomData.size() > 8) {
string chargeData = atomData[8];
if (startsWith(chargeData, "CHG=")) {
int charge = lexicalCast<int>(chargeData.substr(4), ok);
if (!ok) {
appendError("Failed to parse atom charge: " + chargeData);
return false;
// loop through the key=value pairs
for (size_t j = 8; j < atomData.size(); ++j) {
string key = atomData[j];
if (startsWith(key, "CHG=")) {
int charge = lexicalCast<int>(key.substr(4), ok);
if (!ok) {
appendError("Failed to parse atom charge: " + key);
return false;
}
newAtom.setFormalCharge(charge);
} else if (startsWith(key, "RAD=")) {
// radical center
int radical = lexicalCast<int>(key.substr(4), ok);
if (!ok) {
appendError("Failed to parse radical type: " + key);
return false;
}
// we don't set radical centers, just count them
// for the total spin multiplicity
if (radical == 2) // doublet
spinMultiplicity += 1;
else if (radical == 3) // triplet
spinMultiplicity += 2;
} else if (startsWith(key, "ISO=")) {
// isotope
int isotope = lexicalCast<int>(key.substr(4), ok);
if (!ok) {
appendError("Failed to parse isotope type: " + key);
return false;
}
// TODO: handle isotopes
// mol.atom(i).setIsotope(isotope);
}
newAtom.setFormalCharge(charge);
}
} // end of key-value loop
}
} // end of atom block
getline(in, buffer);
Expand Down Expand Up @@ -391,6 +471,10 @@ bool MdlFormat::readV3000(std::istream& in, Core::Molecule& mol)
}
}

// if we have a spin multiplicity, set it
if (spinMultiplicity > 1)
mol.setData("totalSpinMultiplicity", spinMultiplicity);

return true;
}

Expand Down
7 changes: 7 additions & 0 deletions avogadro/qtplugins/molecularproperties/molecularmodel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -194,6 +194,13 @@ QString formatFormula(Molecule* molecule)
else if (charge > 0)
formula += QString("<sup>+%1</sup>").arg(charge);

// add doublet or triplet for spin multiplicity as radical dot
int spinMultiplicity = molecule->totalSpinMultiplicity();
if (spinMultiplicity == 2)
formula += "<sup>•</sup>";
else if (spinMultiplicity == 3)
formula += "<sup>••</sup>";

return formula;
}

Expand Down

0 comments on commit 72fb7c9

Please sign in to comment.