mirror of
https://github.com/rdkit/rdkit.git
synced 2026-06-04 21:54:27 +08:00
8800e1556b
* Do not reset the ringInfo information when adding bonds This call was inconsistent (for example, the version of addBond() in ROMol did not do it) and is unnecessary since the standard assumption is molecules need to be re-sanitized after adding atoms and bonds * response to review clear the property cache on atoms after adding a bond. * add a property cache update to the reaction runner * add something to the release notes
1388 lines
49 KiB
C++
1388 lines
49 KiB
C++
//
|
|
// Copyright (C) 2003-2025 Greg Landrum and other RDKit contributors
|
|
//
|
|
// @@ All Rights Reserved @@
|
|
// This file is part of the RDKit.
|
|
// The contents are covered by the terms of the BSD license
|
|
// which is included in the file license.txt, found at the root
|
|
// of the RDKit source tree.
|
|
//
|
|
#include "RDKitBase.h"
|
|
#include <list>
|
|
#include "QueryAtom.h"
|
|
#include "QueryOps.h"
|
|
#include "MonomerInfo.h"
|
|
#include "Chirality.h"
|
|
#include <Geometry/Transform3D.h>
|
|
#include <Geometry/point.h>
|
|
#include <boost/algorithm/string/classification.hpp>
|
|
#include <boost/dynamic_bitset.hpp>
|
|
#include <boost/range/iterator_range.hpp>
|
|
|
|
constexpr double sq_dist_zero_tol = 1.e-4;
|
|
|
|
namespace RDKit {
|
|
|
|
// Local utility functionality:
|
|
namespace {
|
|
Atom *getAtomNeighborNot(ROMol *mol, const Atom *atom, const Atom *other) {
|
|
PRECONDITION(mol, "bad molecule");
|
|
PRECONDITION(atom, "bad atom");
|
|
PRECONDITION(atom->getDegree() > 1, "bad degree");
|
|
PRECONDITION(other, "bad atom");
|
|
Atom *res = nullptr;
|
|
|
|
ROMol::ADJ_ITER nbrIdx, endNbrs;
|
|
boost::tie(nbrIdx, endNbrs) = mol->getAtomNeighbors(atom);
|
|
while (nbrIdx != endNbrs) {
|
|
if (*nbrIdx != other->getIdx()) {
|
|
res = mol->getAtomWithIdx(*nbrIdx);
|
|
break;
|
|
}
|
|
++nbrIdx;
|
|
}
|
|
|
|
POSTCONDITION(res, "no neighbor found");
|
|
return res;
|
|
}
|
|
|
|
void AssignHsResidueInfo(RWMol &mol) {
|
|
int max_serial = 0;
|
|
unsigned int stopIdx = mol.getNumAtoms();
|
|
for (unsigned int aidx = 0; aidx < stopIdx; ++aidx) {
|
|
auto *info =
|
|
(AtomPDBResidueInfo *)(mol.getAtomWithIdx(aidx)->getMonomerInfo());
|
|
if (info && info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE &&
|
|
info->getSerialNumber() > max_serial) {
|
|
max_serial = info->getSerialNumber();
|
|
}
|
|
}
|
|
|
|
AtomPDBResidueInfo *current_info = nullptr;
|
|
int current_h_id = 0;
|
|
for (unsigned int aidx = 0; aidx < stopIdx; ++aidx) {
|
|
Atom *newAt = mol.getAtomWithIdx(aidx);
|
|
auto *info = (AtomPDBResidueInfo *)(newAt->getMonomerInfo());
|
|
if (info && info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE) {
|
|
ROMol::ADJ_ITER begin, end;
|
|
boost::tie(begin, end) = mol.getAtomNeighbors(newAt);
|
|
while (begin != end) {
|
|
if (mol.getAtomWithIdx(*begin)->getAtomicNum() == 1) {
|
|
// Make all Hs unique - increment id even for existing
|
|
++current_h_id;
|
|
// skip if hydrogen already has PDB info
|
|
auto *h_info = (AtomPDBResidueInfo *)mol.getAtomWithIdx(*begin)
|
|
->getMonomerInfo();
|
|
if (h_info &&
|
|
h_info->getMonomerType() == AtomMonomerInfo::PDBRESIDUE) {
|
|
continue;
|
|
}
|
|
// the hydrogens have unique names on residue basis (H1, H2, ...)
|
|
if (!current_info ||
|
|
current_info->getResidueNumber() != info->getResidueNumber() ||
|
|
current_info->getChainId() != info->getChainId()) {
|
|
current_h_id = 1;
|
|
current_info = info;
|
|
}
|
|
std::string h_label = std::to_string(current_h_id);
|
|
if (h_label.length() > 3) {
|
|
h_label = h_label.substr(h_label.length() - 3, 3);
|
|
}
|
|
while (h_label.length() < 3) {
|
|
h_label = h_label + " ";
|
|
}
|
|
h_label = "H" + h_label;
|
|
// wrap around id to '3H12'
|
|
h_label = h_label.substr(3, 1) + h_label.substr(0, 3);
|
|
AtomPDBResidueInfo *newInfo = new AtomPDBResidueInfo(
|
|
h_label, max_serial, "", info->getResidueName(),
|
|
info->getResidueNumber(), info->getChainId(), "", 1.0, 0.0,
|
|
info->getIsHeteroAtom());
|
|
mol.getAtomWithIdx(*begin)->setMonomerInfo(newInfo);
|
|
|
|
++max_serial;
|
|
}
|
|
++begin;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
std::map<unsigned int, std::vector<unsigned int>> getIsoMap(const ROMol &mol) {
|
|
std::map<unsigned int, std::vector<unsigned int>> isoMap;
|
|
for (auto atom : mol.atoms()) {
|
|
if (atom->hasProp(common_properties::_isotopicHs)) {
|
|
atom->clearProp(common_properties::_isotopicHs);
|
|
}
|
|
}
|
|
for (auto bond : mol.bonds()) {
|
|
auto ba = bond->getBeginAtom();
|
|
auto ea = bond->getEndAtom();
|
|
int ha = -1;
|
|
unsigned int iso;
|
|
if (ba->getAtomicNum() == 1 && ba->getIsotope() &&
|
|
ea->getAtomicNum() != 1) {
|
|
ha = ea->getIdx();
|
|
iso = ba->getIsotope();
|
|
} else if (ea->getAtomicNum() == 1 && ea->getIsotope() &&
|
|
ba->getAtomicNum() != 1) {
|
|
ha = ba->getIdx();
|
|
iso = ea->getIsotope();
|
|
}
|
|
if (ha == -1) {
|
|
continue;
|
|
}
|
|
auto &v = isoMap[ha];
|
|
v.push_back(iso);
|
|
}
|
|
return isoMap;
|
|
}
|
|
|
|
bool may_need_extra_H(const ROMol &mol, const Atom *atom) {
|
|
unsigned single_bonds = 0;
|
|
unsigned aromatic_bonds = 0;
|
|
for (auto bond : mol.atomBonds(atom)) {
|
|
if (bond->getBondType() == Bond::SINGLE) {
|
|
++single_bonds;
|
|
} else if (bond->getBondType() == Bond::AROMATIC) {
|
|
++aromatic_bonds;
|
|
} else {
|
|
return false;
|
|
}
|
|
}
|
|
return single_bonds == 1 && aromatic_bonds == 2 &&
|
|
atom->getTotalValence() == 3;
|
|
}
|
|
|
|
} // end of unnamed namespace
|
|
|
|
namespace MolOps {
|
|
|
|
namespace {
|
|
RDGeom::Point3D pickBisector(const RDGeom::Point3D &nbr1Vect,
|
|
const RDGeom::Point3D &nbr2Vect,
|
|
const RDGeom::Point3D &nbr3Vect) {
|
|
auto dirVect = nbr2Vect + nbr3Vect;
|
|
if (dirVect.lengthSq() < sq_dist_zero_tol) {
|
|
// nbr2Vect and nbr3Vect are anti-parallel (was #3854)
|
|
dirVect = nbr2Vect;
|
|
std::swap(dirVect.x, dirVect.y);
|
|
dirVect.x *= -1;
|
|
}
|
|
if (dirVect.dotProduct(nbr1Vect) < 0) {
|
|
dirVect *= -1;
|
|
}
|
|
|
|
return dirVect;
|
|
}
|
|
} // namespace
|
|
|
|
void setTerminalAtomCoords(ROMol &mol, unsigned int idx,
|
|
unsigned int otherIdx) {
|
|
// we will loop over all the coordinates
|
|
PRECONDITION(otherIdx != idx, "degenerate atoms");
|
|
Atom *atom = mol.getAtomWithIdx(idx);
|
|
PRECONDITION(mol.getAtomDegree(atom) == 1, "bad atom degree");
|
|
const Bond *bond = mol.getBondBetweenAtoms(otherIdx, idx);
|
|
PRECONDITION(bond, "no bond between atoms");
|
|
|
|
const Atom *otherAtom = mol.getAtomWithIdx(otherIdx);
|
|
double bondLength =
|
|
PeriodicTable::getTable()->getRb0(1) +
|
|
PeriodicTable::getTable()->getRb0(otherAtom->getAtomicNum());
|
|
|
|
RDGeom::Point3D dirVect(0, 0, 0);
|
|
|
|
RDGeom::Point3D perpVect, rotnAxis, nbrPerp;
|
|
RDGeom::Point3D nbr1Vect, nbr2Vect, nbr3Vect;
|
|
RDGeom::Transform3D tform;
|
|
RDGeom::Point3D otherPos, atomPos;
|
|
|
|
const Atom *nbr1 = nullptr, *nbr2 = nullptr, *nbr3 = nullptr;
|
|
const Bond *nbrBond;
|
|
ROMol::ADJ_ITER nbrIdx, endNbrs;
|
|
|
|
switch (otherAtom->getDegree()) {
|
|
case 1:
|
|
// --------------------------------------------------------------------------
|
|
// No other atoms present:
|
|
// --------------------------------------------------------------------------
|
|
// loop over the conformations and set the coordinates
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers();
|
|
cfi++) {
|
|
if ((*cfi)->is3D()) {
|
|
dirVect.z = 1;
|
|
} else {
|
|
dirVect.x = 1;
|
|
}
|
|
otherPos = (*cfi)->getAtomPos(otherIdx);
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
}
|
|
break;
|
|
|
|
case 2:
|
|
// --------------------------------------------------------------------------
|
|
// One other neighbor:
|
|
// --------------------------------------------------------------------------
|
|
nbr1 = getAtomNeighborNot(&mol, otherAtom, atom);
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers();
|
|
++cfi) {
|
|
otherPos = (*cfi)->getAtomPos(otherIdx);
|
|
RDGeom::Point3D nbr1Pos = (*cfi)->getAtomPos(nbr1->getIdx());
|
|
// get a normalized vector pointing away from the neighbor:
|
|
nbr1Vect = nbr1Pos - otherPos;
|
|
if (nbr1Vect.lengthSq() < sq_dist_zero_tol) {
|
|
// no difference, which likely indicates that we have redundant atoms.
|
|
// just put it on top of the heavy atom. This was #678
|
|
(*cfi)->setAtomPos(idx, otherPos);
|
|
continue;
|
|
}
|
|
nbr1Vect.normalize();
|
|
nbr1Vect *= -1;
|
|
|
|
// ok, nbr1Vect points away from the other atom, figure out where
|
|
// this H goes:
|
|
switch (otherAtom->getHybridization()) {
|
|
case Atom::SP3:
|
|
// get a perpendicular to nbr1Vect:
|
|
if ((*cfi)->is3D()) {
|
|
perpVect = nbr1Vect.getPerpendicular();
|
|
} else {
|
|
perpVect.z = 1.0;
|
|
}
|
|
// and move off it:
|
|
tform.SetRotation((180 - 109.471) * M_PI / 180., perpVect);
|
|
dirVect = tform * nbr1Vect;
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
case Atom::SP2:
|
|
// default 3D position is to just take an arbitrary perpendicular
|
|
// for 2D we take the normal to the xy plane
|
|
if ((*cfi)->is3D()) {
|
|
perpVect = nbr1Vect.getPerpendicular();
|
|
} else {
|
|
perpVect.z = 1.0;
|
|
}
|
|
if (nbr1->getDegree() > 1) {
|
|
// can we use the neighboring atom to establish a perpendicular?
|
|
nbrBond = mol.getBondBetweenAtoms(otherIdx, nbr1->getIdx());
|
|
if (nbrBond->getIsAromatic() ||
|
|
nbrBond->getBondType() == Bond::DOUBLE ||
|
|
nbrBond->getIsConjugated()) {
|
|
nbr2 = getAtomNeighborNot(&mol, nbr1, otherAtom);
|
|
nbr2Vect =
|
|
nbr1Pos.directionVector((*cfi)->getAtomPos(nbr2->getIdx()));
|
|
auto crossProd = nbr2Vect.crossProduct(nbr1Vect);
|
|
|
|
// if nbr1 and nbr2 are aligned, the perpendicular will be null,
|
|
// and we'll just keep the default calculated above. Otherwise
|
|
// we use the cross product
|
|
if (crossProd.lengthSq() >= sq_dist_zero_tol) {
|
|
perpVect = crossProd;
|
|
}
|
|
}
|
|
}
|
|
perpVect.normalize();
|
|
// rotate the nbr1Vect 60 degrees about perpVect and we're done:
|
|
tform.SetRotation(60. * M_PI / 180., perpVect);
|
|
dirVect = tform * nbr1Vect;
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
case Atom::SP:
|
|
// just lay the H along the vector:
|
|
dirVect = nbr1Vect;
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
default:
|
|
// FIX: handle other hybridizations
|
|
// for now, just lay the H along the vector:
|
|
dirVect = nbr1Vect;
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
}
|
|
}
|
|
break;
|
|
case 3:
|
|
// --------------------------------------------------------------------------
|
|
// Two other neighbors:
|
|
// --------------------------------------------------------------------------
|
|
boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(otherAtom);
|
|
while (nbrIdx != endNbrs) {
|
|
if (*nbrIdx != idx) {
|
|
if (!nbr1) {
|
|
nbr1 = mol.getAtomWithIdx(*nbrIdx);
|
|
} else {
|
|
nbr2 = mol.getAtomWithIdx(*nbrIdx);
|
|
}
|
|
}
|
|
++nbrIdx;
|
|
}
|
|
TEST_ASSERT(nbr1);
|
|
TEST_ASSERT(nbr2);
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers();
|
|
++cfi) {
|
|
// start along the average of the two vectors:
|
|
otherPos = (*cfi)->getAtomPos(otherIdx);
|
|
nbr1Vect = otherPos - (*cfi)->getAtomPos(nbr1->getIdx());
|
|
nbr2Vect = otherPos - (*cfi)->getAtomPos(nbr2->getIdx());
|
|
if (nbr1Vect.lengthSq() < sq_dist_zero_tol ||
|
|
nbr2Vect.lengthSq() < sq_dist_zero_tol) {
|
|
// no difference, which likely indicates that we have redundant atoms.
|
|
// just put it on top of the heavy atom. This was #678
|
|
(*cfi)->setAtomPos(idx, otherPos);
|
|
continue;
|
|
}
|
|
nbr1Vect.normalize();
|
|
nbr2Vect.normalize();
|
|
dirVect = nbr1Vect + nbr2Vect;
|
|
|
|
if (dirVect.lengthSq() < sq_dist_zero_tol) {
|
|
// nbr1Vect and nbr2Vect are non-null, but they may
|
|
// still cancel each other out
|
|
continue;
|
|
}
|
|
dirVect.normalize();
|
|
if ((*cfi)->is3D()) {
|
|
switch (otherAtom->getHybridization()) {
|
|
case Atom::SP3:
|
|
// get the perpendicular to the neighbors:
|
|
nbrPerp = nbr1Vect.crossProduct(nbr2Vect);
|
|
// and the perpendicular to that:
|
|
rotnAxis = nbrPerp.crossProduct(dirVect);
|
|
// and then rotate about that:
|
|
rotnAxis.normalize();
|
|
tform.SetRotation((109.471 / 2) * M_PI / 180., rotnAxis);
|
|
dirVect = tform * dirVect;
|
|
atomPos =
|
|
otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
case Atom::SP2:
|
|
// don't need to do anything here, the H atom goes right on the
|
|
// direction vector
|
|
atomPos =
|
|
otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
default:
|
|
// FIX: handle other hybridizations
|
|
// for now, just lay the H along the neighbor vector;
|
|
atomPos =
|
|
otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
break;
|
|
}
|
|
} else {
|
|
// don't need to do anything here, the H atom goes right on the
|
|
// direction vector
|
|
atomPos = otherPos + dirVect;
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
}
|
|
}
|
|
break;
|
|
case 4:
|
|
// --------------------------------------------------------------------------
|
|
// Three other neighbors:
|
|
// --------------------------------------------------------------------------
|
|
boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(otherAtom);
|
|
|
|
// We're using chiral tag for checking chirality, so we just take the
|
|
// initial order
|
|
while (nbrIdx != endNbrs) {
|
|
if (*nbrIdx != idx) {
|
|
if (!nbr1) {
|
|
nbr1 = mol.getAtomWithIdx(*nbrIdx);
|
|
} else if (!nbr2) {
|
|
nbr2 = mol.getAtomWithIdx(*nbrIdx);
|
|
} else {
|
|
nbr3 = mol.getAtomWithIdx(*nbrIdx);
|
|
}
|
|
}
|
|
++nbrIdx;
|
|
}
|
|
|
|
TEST_ASSERT(nbr1);
|
|
TEST_ASSERT(nbr2);
|
|
TEST_ASSERT(nbr3);
|
|
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers();
|
|
++cfi) {
|
|
otherPos = (*cfi)->getAtomPos(otherIdx);
|
|
nbr1Vect = otherPos - (*cfi)->getAtomPos(nbr1->getIdx());
|
|
nbr2Vect = otherPos - (*cfi)->getAtomPos(nbr2->getIdx());
|
|
nbr3Vect = otherPos - (*cfi)->getAtomPos(nbr3->getIdx());
|
|
if (nbr1Vect.lengthSq() < sq_dist_zero_tol ||
|
|
nbr2Vect.lengthSq() < sq_dist_zero_tol ||
|
|
nbr3Vect.lengthSq() < sq_dist_zero_tol) {
|
|
// no difference, which likely indicates that we have redundant atoms.
|
|
// just put it on top of the heavy atom. This was #678
|
|
(*cfi)->setAtomPos(idx, otherPos);
|
|
continue;
|
|
}
|
|
nbr1Vect.normalize();
|
|
nbr2Vect.normalize();
|
|
nbr3Vect.normalize();
|
|
|
|
// if three neighboring atoms are more or less planar, this
|
|
// is going to be in a quasi-random (but almost definitely bad)
|
|
// direction...
|
|
// correct for this (issue 2951221):
|
|
if ((*cfi)->is3D()) {
|
|
if (fabs(nbr3Vect.dotProduct(nbr1Vect.crossProduct(nbr2Vect))) <
|
|
0.1) {
|
|
// compute the normal:
|
|
dirVect = nbr1Vect.crossProduct(nbr2Vect);
|
|
|
|
// Each of the nbr vectors is non-null, but there might be pairs
|
|
// that cancel each other out. Try to find a direction from atoms
|
|
// that do not overlap.
|
|
if (dirVect.lengthSq() < sq_dist_zero_tol) {
|
|
// This definition of dirVect reverses the parity around otherIdx
|
|
// the change of sign restores it
|
|
dirVect = nbr1Vect.crossProduct(nbr3Vect) * -1;
|
|
}
|
|
if (dirVect.lengthSq() < sq_dist_zero_tol) {
|
|
dirVect = nbr2Vect.crossProduct(nbr3Vect);
|
|
}
|
|
// We couldn't find a good direction
|
|
if (dirVect.lengthSq() < sq_dist_zero_tol) {
|
|
continue;
|
|
}
|
|
|
|
std::string cipCode;
|
|
if (otherAtom->getPropIfPresent(common_properties::_CIPCode,
|
|
cipCode)) {
|
|
// the heavy atom is a chiral center, make sure
|
|
// that we went go the right direction to preserve
|
|
// its chirality. We use the chiral volume for this:
|
|
RDGeom::Point3D v1 = dirVect - nbr3Vect;
|
|
RDGeom::Point3D v2 = nbr1Vect - nbr3Vect;
|
|
RDGeom::Point3D v3 = nbr2Vect - nbr3Vect;
|
|
double vol = v1.dotProduct(v2.crossProduct(v3));
|
|
|
|
if ((otherAtom->getChiralTag() ==
|
|
Atom::ChiralType::CHI_TETRAHEDRAL_CCW &&
|
|
vol < 0) ||
|
|
(otherAtom->getChiralTag() ==
|
|
Atom::ChiralType::CHI_TETRAHEDRAL_CW &&
|
|
vol > 0)) {
|
|
dirVect *= -1;
|
|
}
|
|
}
|
|
} else {
|
|
dirVect = nbr1Vect + nbr2Vect + nbr3Vect;
|
|
}
|
|
} else {
|
|
// we're in flatland
|
|
|
|
// github #3879 and #908: find the two neighbors with the largest
|
|
// outer angle between them and then place the H to bisect that angle
|
|
// This is recommendation ST-1.1.4 from the 2006 IUPAC "Graphical
|
|
// representation of stereochemical configuration" guideline
|
|
auto angle12 = nbr1Vect.angleTo(nbr2Vect);
|
|
auto angle13 = nbr1Vect.angleTo(nbr3Vect);
|
|
auto angle23 = nbr2Vect.angleTo(nbr3Vect);
|
|
auto accum1 = angle12 + angle13;
|
|
auto accum2 = angle12 + angle23;
|
|
auto accum3 = angle13 + angle23;
|
|
if (accum1 <= accum2 && accum1 <= accum3) {
|
|
dirVect = pickBisector(nbr1Vect, nbr2Vect, nbr3Vect);
|
|
} else if (accum2 <= accum1 && accum2 <= accum3) {
|
|
dirVect = pickBisector(nbr2Vect, nbr1Vect, nbr3Vect);
|
|
} else {
|
|
dirVect = pickBisector(nbr3Vect, nbr1Vect, nbr2Vect);
|
|
}
|
|
}
|
|
|
|
dirVect.normalize();
|
|
atomPos = otherPos + dirVect * ((*cfi)->is3D() ? bondLength : 1.0);
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
}
|
|
break;
|
|
default:
|
|
// --------------------------------------------------------------------------
|
|
// FIX: figure out what to do here
|
|
// --------------------------------------------------------------------------
|
|
atomPos = otherPos + dirVect * bondLength;
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers();
|
|
++cfi) {
|
|
(*cfi)->setAtomPos(idx, atomPos);
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
namespace {
|
|
bool isQueryAtom(const RWMol &mol, const Atom &atom) {
|
|
if (atom.hasQuery()) {
|
|
return true;
|
|
}
|
|
for (const auto bnd : mol.atomBonds(&atom)) {
|
|
if (bnd->hasQuery()) {
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
} // namespace
|
|
void addHs(RWMol &mol, const AddHsParameters ¶ms,
|
|
const UINT_VECT *onlyOnAtoms) {
|
|
// when we hit each atom, clear its computed properties
|
|
// NOTE: it is essential that we not clear the ring info in the
|
|
// molecule's computed properties. We don't want to have to
|
|
// regenerate that. This caused Issue210 and Issue212:
|
|
mol.clearComputedProps(false);
|
|
|
|
// precompute the number of hydrogens we are going to add so that we can
|
|
// pre-allocate the necessary space on the conformations of the molecule
|
|
// for their coordinates
|
|
unsigned int numAddHyds = 0;
|
|
boost::dynamic_bitset<> onAtoms(mol.getNumAtoms());
|
|
if (onlyOnAtoms) {
|
|
for (auto atIdx : *onlyOnAtoms) {
|
|
onAtoms.set(atIdx);
|
|
}
|
|
} else {
|
|
onAtoms.set();
|
|
}
|
|
std::vector<unsigned int> numExplicitHs(mol.getNumAtoms(), 0);
|
|
std::vector<unsigned int> numImplicitHs(mol.getNumAtoms(), 0);
|
|
for (auto at : mol.atoms()) {
|
|
numExplicitHs[at->getIdx()] = at->getNumExplicitHs();
|
|
numImplicitHs[at->getIdx()] = at->getNumImplicitHs();
|
|
if (onAtoms[at->getIdx()]) {
|
|
if (params.skipQueries && isQueryAtom(mol, *at)) {
|
|
onAtoms.set(at->getIdx(), 0);
|
|
continue;
|
|
}
|
|
numAddHyds += at->getNumExplicitHs();
|
|
if (!params.explicitOnly) {
|
|
numAddHyds += at->getNumImplicitHs();
|
|
}
|
|
}
|
|
}
|
|
unsigned int nSize = mol.getNumAtoms() + numAddHyds;
|
|
|
|
// loop over the conformations of the molecule and allocate new space
|
|
// for the H locations (need to do this even if we aren't adding coords so
|
|
// that the conformers have the correct number of atoms).
|
|
for (auto cfi = mol.beginConformers(); cfi != mol.endConformers(); ++cfi) {
|
|
(*cfi)->reserve(nSize);
|
|
}
|
|
|
|
unsigned int stopIdx = mol.getNumAtoms();
|
|
for (unsigned int aidx = 0; aidx < stopIdx; ++aidx) {
|
|
if (!onAtoms[aidx]) {
|
|
continue;
|
|
}
|
|
|
|
Atom *newAt = mol.getAtomWithIdx(aidx);
|
|
|
|
std::vector<unsigned int> isoHs;
|
|
if (newAt->getPropIfPresent(common_properties::_isotopicHs, isoHs)) {
|
|
newAt->clearProp(common_properties::_isotopicHs);
|
|
}
|
|
std::vector<unsigned int>::const_iterator isoH = isoHs.begin();
|
|
unsigned int newIdx;
|
|
newAt->clearComputedProps();
|
|
// always convert explicit Hs
|
|
unsigned int onumexpl = numExplicitHs[aidx];
|
|
for (unsigned int i = 0; i < onumexpl; i++) {
|
|
newIdx = mol.addAtom(new Atom(1), false, true);
|
|
mol.addBond(aidx, newIdx, Bond::SINGLE);
|
|
auto hAtom = mol.getAtomWithIdx(newIdx);
|
|
hAtom->updatePropertyCache();
|
|
if (params.addCoords) {
|
|
setTerminalAtomCoords(mol, newIdx, aidx);
|
|
}
|
|
if (isoH != isoHs.end()) {
|
|
hAtom->setIsotope(*isoH);
|
|
++isoH;
|
|
}
|
|
}
|
|
// clear the local property
|
|
newAt->setNumExplicitHs(0);
|
|
|
|
if (!params.explicitOnly) {
|
|
// take care of implicits
|
|
for (unsigned int i = 0; i < numImplicitHs[aidx]; i++) {
|
|
newIdx = mol.addAtom(new Atom(1), false, true);
|
|
mol.addBond(aidx, newIdx, Bond::SINGLE);
|
|
// set the isImplicit label so that we can strip these back
|
|
// off later if need be.
|
|
auto hAtom = mol.getAtomWithIdx(newIdx);
|
|
hAtom->setProp(common_properties::isImplicit, 1);
|
|
hAtom->updatePropertyCache();
|
|
if (params.addCoords) {
|
|
setTerminalAtomCoords(mol, newIdx, aidx);
|
|
}
|
|
if (isoH != isoHs.end()) {
|
|
hAtom->setIsotope(*isoH);
|
|
++isoH;
|
|
}
|
|
}
|
|
}
|
|
// update the atom's derived properties (valence count, etc.)
|
|
// no sense in being strict here (was github #2782)
|
|
newAt->updatePropertyCache(false);
|
|
if (isoH != isoHs.end()) {
|
|
BOOST_LOG(rdWarningLog) << "extra H isotope information found on atom "
|
|
<< newAt->getIdx() << std::endl;
|
|
}
|
|
}
|
|
// take care of AtomPDBResidueInfo for Hs if root atom has it
|
|
if (params.addResidueInfo) {
|
|
AssignHsResidueInfo(mol);
|
|
}
|
|
}
|
|
|
|
namespace {
|
|
// returns whether or not an adjustment was made, in case we want that info
|
|
bool adjustStereoAtomsIfRequired(RWMol &mol, const Atom *atom,
|
|
const Atom *heavyAtom) {
|
|
PRECONDITION(atom != nullptr, "bad atom");
|
|
PRECONDITION(heavyAtom != nullptr, "bad heavy atom");
|
|
// nothing we can do if the degree is only 2 (and we should have covered
|
|
// that earlier anyway)
|
|
if (heavyAtom->getDegree() == 2) {
|
|
return false;
|
|
}
|
|
const auto &cbnd =
|
|
mol.getBondBetweenAtoms(atom->getIdx(), heavyAtom->getIdx());
|
|
if (!cbnd) {
|
|
return false;
|
|
}
|
|
for (const auto &nbri :
|
|
boost::make_iterator_range(mol.getAtomBonds(heavyAtom))) {
|
|
Bond *bnd = mol[nbri];
|
|
if (bnd->getBondType() == Bond::DOUBLE &&
|
|
bnd->getStereo() > Bond::STEREOANY) {
|
|
auto sAtomIt = std::find(bnd->getStereoAtoms().begin(),
|
|
bnd->getStereoAtoms().end(), atom->getIdx());
|
|
if (sAtomIt != bnd->getStereoAtoms().end()) {
|
|
// sAtomIt points to the position of this atom's index in the list.
|
|
// find the index of another atom attached to the heavy atom and
|
|
// use it to update sAtomIt
|
|
unsigned int dblNbrIdx = bnd->getOtherAtomIdx(heavyAtom->getIdx());
|
|
for (const auto &nbri :
|
|
boost::make_iterator_range(mol.getAtomNeighbors(heavyAtom))) {
|
|
const auto &nbr = mol[nbri];
|
|
if (nbr->getIdx() == dblNbrIdx || nbr->getIdx() == atom->getIdx()) {
|
|
continue;
|
|
}
|
|
*sAtomIt = nbr->getIdx();
|
|
bool madeAdjustment = true;
|
|
switch (bnd->getStereo()) {
|
|
case Bond::STEREOCIS:
|
|
bnd->setStereo(Bond::STEREOTRANS);
|
|
break;
|
|
case Bond::STEREOTRANS:
|
|
bnd->setStereo(Bond::STEREOCIS);
|
|
break;
|
|
default:
|
|
// I think we shouldn't need to do anything with E and Z...
|
|
madeAdjustment = false;
|
|
break;
|
|
}
|
|
return madeAdjustment;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
void molRemoveH(RWMol &mol, unsigned int idx, bool updateExplicitCount) {
|
|
auto atom = mol.getAtomWithIdx(idx);
|
|
PRECONDITION(atom->getAtomicNum() == 1, "idx corresponds to a non-Hydrogen");
|
|
for (const auto bond : mol.atomBonds(atom)) {
|
|
Atom *heavyAtom = bond->getOtherAtom(atom);
|
|
int heavyAtomNum = heavyAtom->getAtomicNum();
|
|
|
|
// we'll update the neighbor's explicit H count if we were told to
|
|
// *or* if the neighbor is chiral, in which case the H is needed
|
|
// in order to complete the coordination
|
|
// *or* if the neighbor has the noImplicit flag set:
|
|
if (updateExplicitCount || heavyAtom->getNoImplicit() ||
|
|
heavyAtom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
|
|
heavyAtom->setNumExplicitHs(heavyAtom->getNumExplicitHs() + 1);
|
|
} else {
|
|
// this is a special case related to Issue 228 and the
|
|
// "disappearing Hydrogen" problem discussed in MolOps::adjustHs
|
|
//
|
|
// If we remove a hydrogen from an aromatic N or P, or if
|
|
// the heavy atom it is connected to is not in its default
|
|
// valence state, we need to be *sure* to increment the
|
|
// explicit count, even if the H itself isn't marked as explicit
|
|
const INT_VECT &defaultVs =
|
|
PeriodicTable::getTable()->getValenceList(heavyAtomNum);
|
|
if (((heavyAtomNum == 7 || heavyAtomNum == 15 ||
|
|
may_need_extra_H(mol, heavyAtom)) &&
|
|
heavyAtom->getIsAromatic()) ||
|
|
(std::find(defaultVs.begin() + 1, defaultVs.end(),
|
|
heavyAtom->getTotalValence()) != defaultVs.end())) {
|
|
heavyAtom->setNumExplicitHs(heavyAtom->getNumExplicitHs() + 1);
|
|
}
|
|
}
|
|
|
|
// One other consequence of removing the H from the graph is
|
|
// that we may change the ordering of the bonds about a
|
|
// chiral center. This may change the chiral label at that
|
|
// atom. We deal with that by explicitly checking here:
|
|
if (heavyAtom->getChiralTag() != Atom::CHI_UNSPECIFIED) {
|
|
INT_LIST neighborIndices;
|
|
for (const auto &nbnd : mol.atomBonds(heavyAtom)) {
|
|
if (nbnd->getIdx() != bond->getIdx()) {
|
|
neighborIndices.push_back(nbnd->getIdx());
|
|
}
|
|
}
|
|
neighborIndices.push_back(bond->getIdx());
|
|
|
|
int nSwaps = heavyAtom->getPerturbationOrder(neighborIndices);
|
|
// std::cerr << "H: "<<atom->getIdx()<<" hvy:
|
|
// "<<heavyAtom->getIdx()<<" swaps: " << nSwaps<<std::endl;
|
|
if (nSwaps % 2) {
|
|
heavyAtom->invertChirality();
|
|
}
|
|
}
|
|
|
|
// If we are removing a H atom that defines bond stereo (e.g. imines),
|
|
// Then also remove the bond stereo information, as it is no longer valid.
|
|
if (heavyAtom->getDegree() == 2) {
|
|
for (auto &nbnd : mol.atomBonds(heavyAtom)) {
|
|
if (nbnd != bond) {
|
|
if (nbnd->getStereo() > Bond::STEREOANY) {
|
|
nbnd->setStereo(Bond::STEREONONE);
|
|
nbnd->getStereoAtoms().clear();
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
// if it's a wavy bond, then we need to
|
|
// mark the beginning atom with the _UnknownStereo tag.
|
|
// so that we know later that something was affecting its
|
|
// stereochem
|
|
if (bond->getBondDir() == Bond::UNKNOWN &&
|
|
bond->getBeginAtomIdx() == heavyAtom->getIdx()) {
|
|
heavyAtom->setProp(common_properties::_UnknownStereo, 1);
|
|
} else if (bond->getBondDir() == Bond::ENDDOWNRIGHT ||
|
|
bond->getBondDir() == Bond::ENDUPRIGHT) {
|
|
// if the direction is set on this bond and the atom it's connected to
|
|
// has no other single bonds with directions set, then we need to set
|
|
// direction on one of the other neighbors in order to avoid double
|
|
// bond stereochemistry possibly being lost. This was github #754
|
|
bool foundADir = false;
|
|
Bond *oBond = nullptr;
|
|
for (const auto &nbri :
|
|
boost::make_iterator_range(mol.getAtomBonds(heavyAtom))) {
|
|
Bond *nbnd = mol[nbri];
|
|
if (nbnd->getIdx() != bond->getIdx() &&
|
|
nbnd->getBondType() == Bond::SINGLE) {
|
|
if (nbnd->getBondDir() == Bond::NONE) {
|
|
oBond = nbnd;
|
|
} else {
|
|
foundADir = true;
|
|
}
|
|
}
|
|
}
|
|
if (!foundADir && oBond != nullptr) {
|
|
bool flipIt = (oBond->getBeginAtom() == heavyAtom) &&
|
|
(bond->getBeginAtom() == heavyAtom);
|
|
if (flipIt) {
|
|
oBond->setBondDir(bond->getBondDir() == Bond::ENDDOWNRIGHT
|
|
? Bond::ENDUPRIGHT
|
|
: Bond::ENDDOWNRIGHT);
|
|
} else {
|
|
oBond->setBondDir(bond->getBondDir());
|
|
}
|
|
}
|
|
// if this atom is one of the stereoatoms for a double bond we need
|
|
// to switch the stereo atom on this end to be the other neighbor
|
|
// This was part of github #1810
|
|
adjustStereoAtomsIfRequired(mol, atom, heavyAtom);
|
|
} else {
|
|
// if this atom is one of the stereoatoms for a double bond we need
|
|
// to switch the stereo atom on this end to be the other neighbor
|
|
// This was part of github #1810
|
|
adjustStereoAtomsIfRequired(mol, atom, heavyAtom);
|
|
}
|
|
|
|
// remove the bond from any SGroups that might include it.
|
|
for (auto &sg : getSubstanceGroups(mol)) {
|
|
sg.removeBondWithIdx(bond->getIdx());
|
|
}
|
|
}
|
|
|
|
// Finally, remove the atom from any SGroups that might include it, so that
|
|
// the SGroups don't get removed in removeAtom(). Since we allow removing
|
|
// SGroup SAP lvidx H atoms, we need to check for those and update them.
|
|
for (auto &sg : getSubstanceGroups(mol)) {
|
|
sg.removeAtomWithIdx(idx);
|
|
sg.removeParentAtomWithIdx(idx);
|
|
|
|
for (auto &sap : sg.getAttachPoints()) {
|
|
if (sap.lvIdx == static_cast<int>(idx)) {
|
|
sap.lvIdx = -1;
|
|
}
|
|
}
|
|
}
|
|
// computed properties will be cleared after all hydrogens are removed
|
|
bool clearProps = false;
|
|
mol.removeAtom(atom, clearProps);
|
|
}
|
|
|
|
bool shouldRemoveH(const RWMol &mol, const Atom *atom,
|
|
const RemoveHsParameters &ps) {
|
|
if (atom->getAtomicNum() != 1) {
|
|
return false;
|
|
}
|
|
if (!ps.removeWithQuery && atom->hasQuery()) {
|
|
return false;
|
|
}
|
|
if (!ps.removeDegreeZero && !atom->getDegree()) {
|
|
if (ps.showWarnings) {
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "WARNING: not removing hydrogen atom without neighbors"
|
|
<< std::endl;
|
|
}
|
|
return false;
|
|
}
|
|
if (!ps.removeHigherDegrees && atom->getDegree() > 1) {
|
|
return false;
|
|
}
|
|
if (!ps.removeIsotopes && !ps.removeAndTrackIsotopes && atom->getIsotope()) {
|
|
return false;
|
|
}
|
|
if (!ps.removeNonimplicit && !atom->hasProp(common_properties::isImplicit)) {
|
|
return false;
|
|
}
|
|
if (!ps.removeMapped && atom->getAtomMapNum()) {
|
|
return false;
|
|
}
|
|
|
|
if (ps.removeInSGroups) {
|
|
// If removing H in SGroups, do not remove H atoms in special
|
|
// roles in the SGroup
|
|
for (const auto &sg : getSubstanceGroups(mol)) {
|
|
// The H atom is one of the "caps" of the SGroup. Technically,
|
|
// it's not part of the group, but it defines its boundaries.
|
|
for (const auto &bond_idx : sg.getBonds()) {
|
|
if (sg.getBondType(bond_idx) == SubstanceGroup::BondType::XBOND) {
|
|
auto bond = mol.getBondWithIdx(bond_idx);
|
|
if (bond->getBeginAtom() == atom || bond->getEndAtom() == atom) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
|
|
for (const auto &sap : sg.getAttachPoints()) {
|
|
// The H atoms is an attach point. This would be weird, but is possible.
|
|
// (if it is a 'leaving atom' we don't care, though)
|
|
if (sap.aIdx == atom->getIdx()) {
|
|
return false;
|
|
}
|
|
}
|
|
|
|
for (const auto &cs : sg.getCStates()) {
|
|
// The bond to the H atom defines a CState
|
|
auto bond = mol.getBondWithIdx(cs.bondIdx);
|
|
if (bond->getBeginAtom() == atom || bond->getEndAtom() == atom) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
} else {
|
|
for (const auto &sg : getSubstanceGroups(mol)) {
|
|
if (sg.includesAtom(atom->getIdx())) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
if (!ps.removeHydrides && atom->getFormalCharge() == -1) {
|
|
return false;
|
|
}
|
|
bool removeIt = true;
|
|
if (atom->getDegree() &&
|
|
(!ps.removeDummyNeighbors || !ps.removeDefiningBondStereo ||
|
|
!ps.removeOnlyHNeighbors || !ps.removeNontetrahedralNeighbors ||
|
|
!ps.removeWithWedgedBond)) {
|
|
bool onlyHNeighbors = true;
|
|
for (const auto nbr : mol.atomNeighbors(atom)) {
|
|
// is it a dummy?
|
|
if (!ps.removeDummyNeighbors && nbr->getAtomicNum() < 1) {
|
|
if (ps.showWarnings) {
|
|
BOOST_LOG(rdWarningLog) << "WARNING: not removing hydrogen atom "
|
|
"with dummy atom neighbors"
|
|
<< std::endl;
|
|
}
|
|
return false;
|
|
}
|
|
// does it have non-tetrahedral stereo:
|
|
if (!ps.removeNontetrahedralNeighbors &&
|
|
Chirality::hasNonTetrahedralStereo(nbr)) {
|
|
if (ps.showWarnings) {
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "WARNING: not removing hydrogen atom "
|
|
"with neighbor that has non-tetrahedral stereochemistry"
|
|
<< std::endl;
|
|
}
|
|
return false;
|
|
}
|
|
if (!ps.removeOnlyHNeighbors && nbr->getAtomicNum() != 1) {
|
|
onlyHNeighbors = false;
|
|
}
|
|
if (!ps.removeWithWedgedBond) {
|
|
const auto bnd = mol.getBondBetweenAtoms(atom->getIdx(), nbr->getIdx());
|
|
if (bnd->getBondDir() == Bond::BEGINDASH ||
|
|
bnd->getBondDir() == Bond::BEGINWEDGE) {
|
|
if (ps.showWarnings) {
|
|
BOOST_LOG(rdWarningLog) << "WARNING: not removing hydrogen atom "
|
|
"with wedged bond"
|
|
<< std::endl;
|
|
}
|
|
return false;
|
|
}
|
|
}
|
|
// Check to see if the neighbor has a double bond and we're the only
|
|
// neighbor at this end. This was part of github #1810
|
|
if (!ps.removeDefiningBondStereo && nbr->getDegree() == 2) {
|
|
for (const auto bnd : mol.atomBonds(nbr)) {
|
|
if (bnd->getBondType() == Bond::DOUBLE &&
|
|
(bnd->getStereo() > Bond::STEREOANY ||
|
|
mol.getBondBetweenAtoms(atom->getIdx(), nbr->getIdx())
|
|
->getBondDir() > Bond::NONE)) {
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
if (removeIt && (!ps.removeOnlyHNeighbors && onlyHNeighbors)) {
|
|
return false;
|
|
}
|
|
}
|
|
return removeIt;
|
|
}
|
|
|
|
// Do not remove H atoms that are part of SGroups that only contain H atoms.
|
|
void filter_sgroup_emptying_hydrogens(const ROMol &mol,
|
|
boost::dynamic_bitset<> &atomsToRemove) {
|
|
for (const auto &sg : getSubstanceGroups(mol)) {
|
|
const auto &atoms = sg.getAtoms();
|
|
const auto &patoms = sg.getParentAtoms();
|
|
|
|
// If the SGroup already didn't have atoms, we don't care about it
|
|
if (atoms.empty() && patoms.empty()) {
|
|
continue;
|
|
}
|
|
|
|
auto would_remove_atom = [&atomsToRemove](const auto idx) {
|
|
return atomsToRemove[idx];
|
|
};
|
|
|
|
auto no_atoms = atoms.empty() ||
|
|
std::all_of(atoms.begin(), atoms.end(), would_remove_atom);
|
|
if (no_atoms) {
|
|
auto no_patoms =
|
|
patoms.empty() ||
|
|
std::all_of(patoms.begin(), patoms.end(), would_remove_atom);
|
|
if (no_patoms) {
|
|
for (auto atom : atoms) {
|
|
atomsToRemove.set(atom, false);
|
|
}
|
|
for (auto patom : patoms) {
|
|
atomsToRemove.set(patom, false);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
} // end of anonymous namespace
|
|
|
|
void removeHs(RWMol &mol, const RemoveHsParameters &ps, bool sanitize) {
|
|
if (ps.removeAndTrackIsotopes) {
|
|
// if there are any non-isotopic Hs remove them first
|
|
// to make sure chirality is preserved
|
|
bool needRemoveHs = false;
|
|
for (auto atom : mol.atoms()) {
|
|
if (atom->getAtomicNum() == 1 && atom->getIsotope() == 0) {
|
|
needRemoveHs = true;
|
|
break;
|
|
}
|
|
}
|
|
if (needRemoveHs) {
|
|
RemoveHsParameters psCopy(ps);
|
|
psCopy.removeAndTrackIsotopes = false;
|
|
psCopy.removeIsotopes = false;
|
|
removeHs(mol, psCopy, false);
|
|
}
|
|
}
|
|
for (auto atom : mol.atoms()) {
|
|
atom->updatePropertyCache(false);
|
|
}
|
|
if (ps.removeAndTrackIsotopes) {
|
|
for (const auto &pair : getIsoMap(mol)) {
|
|
mol.getAtomWithIdx(pair.first)
|
|
->setProp(common_properties::_isotopicHs, pair.second);
|
|
}
|
|
}
|
|
boost::dynamic_bitset<> atomsToRemove{mol.getNumAtoms(), 0};
|
|
|
|
for (auto atom : mol.atoms()) {
|
|
if (shouldRemoveH(mol, atom, ps)) {
|
|
atomsToRemove.set(atom->getIdx());
|
|
}
|
|
} // end of the loop over atoms
|
|
|
|
// Once we know which H atoms would be removed, filter out those that
|
|
// would cause any SGroups to become empty
|
|
if (ps.removeInSGroups) {
|
|
filter_sgroup_emptying_hydrogens(mol, atomsToRemove);
|
|
}
|
|
|
|
// now that we know which atoms need to be removed, go ahead and remove them
|
|
// NOTE: there's too much complexity around stereochemistry here
|
|
// to be able to safely use batch editing.
|
|
for (int idx = mol.getNumAtoms() - 1; idx >= 0; --idx) {
|
|
if (atomsToRemove[idx]) {
|
|
molRemoveH(mol, idx, ps.updateExplicitCount);
|
|
}
|
|
}
|
|
mol.clearComputedProps(true);
|
|
//
|
|
// If we didn't only remove implicit Hs, which are guaranteed to
|
|
// be the highest numbered atoms, we may have altered atom indices.
|
|
// This can screw up derived properties (such as ring members), so
|
|
// do some checks:
|
|
//
|
|
if (!atomsToRemove.empty() && ps.removeNonimplicit && sanitize) {
|
|
sanitizeMol(mol);
|
|
}
|
|
};
|
|
ROMol *removeHs(const ROMol &mol, const RemoveHsParameters &ps, bool sanitize) {
|
|
auto *res = new RWMol(mol);
|
|
try {
|
|
removeHs(*res, ps, sanitize);
|
|
} catch (const MolSanitizeException &) {
|
|
delete res;
|
|
throw;
|
|
}
|
|
return static_cast<ROMol *>(res);
|
|
}
|
|
void removeHs(RWMol &mol, bool implicitOnly, bool updateExplicitCount,
|
|
bool sanitize) {
|
|
RemoveHsParameters ps;
|
|
ps.removeNonimplicit = !implicitOnly;
|
|
ps.updateExplicitCount = updateExplicitCount;
|
|
removeHs(mol, ps, sanitize);
|
|
};
|
|
ROMol *removeHs(const ROMol &mol, bool implicitOnly, bool updateExplicitCount,
|
|
bool sanitize) {
|
|
auto *res = new RWMol(mol);
|
|
RemoveHsParameters ps;
|
|
ps.removeNonimplicit = !implicitOnly;
|
|
ps.updateExplicitCount = updateExplicitCount;
|
|
try {
|
|
removeHs(*res, ps, sanitize);
|
|
} catch (const MolSanitizeException &) {
|
|
delete res;
|
|
throw;
|
|
}
|
|
return static_cast<ROMol *>(res);
|
|
}
|
|
|
|
void removeAllHs(RWMol &mol, bool sanitize) {
|
|
RemoveHsParameters ps;
|
|
ps.removeDegreeZero = true;
|
|
ps.removeHigherDegrees = true;
|
|
ps.removeOnlyHNeighbors = true;
|
|
ps.removeIsotopes = true;
|
|
ps.removeDummyNeighbors = true;
|
|
ps.removeDefiningBondStereo = true;
|
|
ps.removeWithWedgedBond = true;
|
|
ps.removeWithQuery = true;
|
|
ps.removeNonimplicit = true;
|
|
ps.removeInSGroups = true;
|
|
ps.showWarnings = false;
|
|
ps.removeHydrides = true;
|
|
ps.removeNontetrahedralNeighbors = true;
|
|
removeHs(mol, ps, sanitize);
|
|
};
|
|
ROMol *removeAllHs(const ROMol &mol, bool sanitize) {
|
|
auto *res = new RWMol(mol);
|
|
try {
|
|
removeAllHs(*res, sanitize);
|
|
} catch (const MolSanitizeException &) {
|
|
delete res;
|
|
throw;
|
|
}
|
|
return static_cast<ROMol *>(res);
|
|
}
|
|
|
|
namespace {
|
|
enum class HydrogenType {
|
|
NotAHydrogen,
|
|
UnMergableQueryHydrogen,
|
|
QueryHydrogen
|
|
};
|
|
|
|
template <class Q>
|
|
std::pair<bool, bool> queryHasHs(Q queryAtom, bool inor = false) {
|
|
for (auto childit = queryAtom->beginChildren();
|
|
childit != queryAtom->endChildren(); ++childit) {
|
|
QueryAtom::QUERYATOM_QUERY::CHILD_TYPE query = *childit;
|
|
if (query->getDescription() == "AtomOr") {
|
|
return queryHasHs(query, true);
|
|
} else if (query->getDescription() == "AtomAtomicNum") {
|
|
if (static_cast<ATOM_EQUALS_QUERY *>(query.get())->getVal() == 1 &&
|
|
!query->getNegation()) {
|
|
return std::make_pair(true, inor);
|
|
}
|
|
} else if (query->getDescription() == "AtomType") {
|
|
auto val = static_cast<ATOM_EQUALS_QUERY *>(query.get())->getVal();
|
|
// 1001 == aromtic hydrogen (not a thing, really)
|
|
// 1 == aliphatic hydrogen
|
|
if ((val == 1001 || val == 1) && !query->getNegation()) {
|
|
return std::make_pair(true, inor);
|
|
}
|
|
}
|
|
}
|
|
return std::make_pair(false, inor);
|
|
;
|
|
}
|
|
|
|
HydrogenType isQueryH(const Atom *atom) {
|
|
PRECONDITION(atom, "bogus atom");
|
|
if (atom->getAtomicNum() == 1) {
|
|
// the simple case: the atom is flagged as being an H and
|
|
// has no query
|
|
if (!atom->hasQuery() ||
|
|
(!atom->getQuery()->getNegation() &&
|
|
atom->getQuery()->getDescription() == "AtomAtomicNum")) {
|
|
return HydrogenType::QueryHydrogen;
|
|
}
|
|
}
|
|
|
|
if (!(atom->getDegree() <= 1)) {
|
|
// bonded and unbonded H atoms will continue rest will be returned
|
|
return HydrogenType::NotAHydrogen;
|
|
}
|
|
|
|
if (atom->hasQuery() && atom->getQuery()->getNegation()) {
|
|
// we will not merge negated queries
|
|
return HydrogenType::NotAHydrogen;
|
|
}
|
|
|
|
if (atom->hasQuery()) {
|
|
std::pair<bool, bool> res = std::make_pair(false, false);
|
|
if (atom->getQuery()->getDescription() == "AtomOr") {
|
|
res = queryHasHs(atom->getQuery(), true);
|
|
} else if (atom->getQuery()->getDescription() == "AtomAnd") {
|
|
res = queryHasHs(atom->getQuery(), false);
|
|
}
|
|
if (res.first) { // hasH
|
|
if (res.second) { // inOr
|
|
BOOST_LOG(rdWarningLog)
|
|
<< "WARNING: merging explicit H queries involved "
|
|
"in ORs is not supported. This query will not "
|
|
"be merged"
|
|
<< std::endl;
|
|
return HydrogenType::UnMergableQueryHydrogen;
|
|
} else {
|
|
return HydrogenType::QueryHydrogen;
|
|
}
|
|
}
|
|
}
|
|
return HydrogenType::NotAHydrogen;
|
|
}
|
|
} // namespace
|
|
|
|
//
|
|
// This routine removes explicit hydrogens (and bonds to them) from
|
|
// the molecular graph and adds them as queries to the heavy atoms
|
|
// to which they are bound. If the heavy atoms (or atom queries)
|
|
// already have hydrogen-count queries, they will be updated.
|
|
//
|
|
// NOTE:
|
|
// - Hydrogens which aren't connected to a heavy atom will not be
|
|
// removed. This prevents molecules like "[H][H]" from having
|
|
// all atoms removed.
|
|
//
|
|
// - By default all hydrogens are removed, however if
|
|
// merge_unmapped_only is true, any hydrogen participating
|
|
// in an atom map will be retained
|
|
void mergeQueryHs(RWMol &mol, bool mergeUnmappedOnly, bool mergeIsotopes) {
|
|
std::vector<unsigned int> atomsToRemove;
|
|
|
|
boost::dynamic_bitset<> hatoms(mol.getNumAtoms());
|
|
for (unsigned int i = 0; i < mol.getNumAtoms(); ++i) {
|
|
hatoms[i] = isQueryH(mol.getAtomWithIdx(i)) == HydrogenType::QueryHydrogen;
|
|
}
|
|
unsigned int currIdx = 0, stopIdx = mol.getNumAtoms();
|
|
while (currIdx < stopIdx) {
|
|
Atom *atom = mol.getAtomWithIdx(currIdx);
|
|
if (!hatoms[currIdx]) {
|
|
unsigned int numHsToRemove = 0;
|
|
ROMol::ADJ_ITER begin, end;
|
|
boost::tie(begin, end) = mol.getAtomNeighbors(atom);
|
|
|
|
while (begin != end) {
|
|
if (hatoms[*begin]) {
|
|
Atom &bgn = *mol.getAtomWithIdx(*begin);
|
|
bool checkUnmapped =
|
|
!mergeUnmappedOnly ||
|
|
!bgn.hasProp(common_properties::molAtomMapNumber);
|
|
bool checkIsotope = mergeIsotopes || bgn.getIsotope() == 0;
|
|
if (checkUnmapped && checkIsotope) {
|
|
atomsToRemove.push_back(rdcast<unsigned int>(*begin));
|
|
++numHsToRemove;
|
|
}
|
|
}
|
|
++begin;
|
|
}
|
|
if (numHsToRemove) {
|
|
//
|
|
// We have H neighbors:
|
|
// Add the appropriate queries to compensate for their removal.
|
|
//
|
|
// Examples:
|
|
// C[H] -> [C;!H0]
|
|
// C([H])[H] -> [C;!H0;!H1]
|
|
//
|
|
// It would be more efficient to do this using range queries like:
|
|
// C([H])[H] -> [C;H{2-}]
|
|
// but that would produce non-standard SMARTS without the user
|
|
// having started with a non-standard SMARTS.
|
|
//
|
|
if (!atom->hasQuery()) {
|
|
// it wasn't a query atom, we need to replace it so that we can add
|
|
// a query:
|
|
ATOM_EQUALS_QUERY *tmp = makeAtomNumQuery(atom->getAtomicNum());
|
|
auto *newAt = new QueryAtom;
|
|
newAt->setQuery(tmp);
|
|
newAt->updateProps(*atom);
|
|
mol.replaceAtom(atom->getIdx(), newAt);
|
|
delete newAt;
|
|
atom = mol.getAtomWithIdx(currIdx);
|
|
}
|
|
for (unsigned int i = 0; i < numHsToRemove; ++i) {
|
|
ATOM_EQUALS_QUERY *tmp = makeAtomHCountQuery(i);
|
|
tmp->setNegation(true);
|
|
atom->expandQuery(tmp);
|
|
}
|
|
} // end of numHsToRemove test
|
|
|
|
// recurse if needed (was github isusue 544)
|
|
if (atom->hasQuery()) {
|
|
if (atom->getQuery()->getDescription() == "RecursiveStructure") {
|
|
auto *rsq = dynamic_cast<RecursiveStructureQuery *>(atom->getQuery());
|
|
CHECK_INVARIANT(rsq, "could not convert recursive structure query");
|
|
RWMol *rqm = new RWMol(*rsq->getQueryMol());
|
|
mergeQueryHs(*rqm, mergeUnmappedOnly, mergeIsotopes);
|
|
rsq->setQueryMol(rqm);
|
|
}
|
|
|
|
// FIX: shouldn't be repeating this code here
|
|
std::list<QueryAtom::QUERYATOM_QUERY::CHILD_TYPE> childStack(
|
|
atom->getQuery()->beginChildren(), atom->getQuery()->endChildren());
|
|
while (childStack.size()) {
|
|
QueryAtom::QUERYATOM_QUERY::CHILD_TYPE qry = childStack.front();
|
|
childStack.pop_front();
|
|
if (qry->getDescription() == "RecursiveStructure") {
|
|
auto *rsq = dynamic_cast<RecursiveStructureQuery *>(qry.get());
|
|
CHECK_INVARIANT(rsq, "could not convert recursive structure query");
|
|
RWMol *rqm = new RWMol(*rsq->getQueryMol());
|
|
mergeQueryHs(*rqm, mergeUnmappedOnly, mergeIsotopes);
|
|
rsq->setQueryMol(rqm);
|
|
} else if (qry->beginChildren() != qry->endChildren()) {
|
|
childStack.insert(childStack.end(), qry->beginChildren(),
|
|
qry->endChildren());
|
|
}
|
|
}
|
|
} // end of recursion loop
|
|
}
|
|
++currIdx;
|
|
}
|
|
mol.beginBatchEdit();
|
|
for (auto aidx : atomsToRemove) {
|
|
mol.removeAtom(aidx);
|
|
}
|
|
mol.commitBatchEdit();
|
|
};
|
|
ROMol *mergeQueryHs(const ROMol &mol, bool mergeUnmappedOnly,
|
|
bool mergeIsotopes) {
|
|
auto *res = new RWMol(mol);
|
|
mergeQueryHs(*res, mergeUnmappedOnly, mergeIsotopes);
|
|
return static_cast<ROMol *>(res);
|
|
};
|
|
|
|
bool needsHs(const ROMol &mol) {
|
|
for (const auto atom : mol.atoms()) {
|
|
bool includeNeighbors = false;
|
|
if (atom->getTotalNumHs(includeNeighbors)) {
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
std::pair<bool, bool> hasQueryHs(const ROMol &mol) {
|
|
bool queryHs = false;
|
|
// We don't care about announcing ORs or other items during isQueryH
|
|
RDLog::LogStateSetter blocker;
|
|
|
|
for (const auto atom : mol.atoms()) {
|
|
switch (isQueryH(atom)) {
|
|
case HydrogenType::UnMergableQueryHydrogen:
|
|
return std::make_pair(true, true);
|
|
case HydrogenType::QueryHydrogen:
|
|
queryHs = true;
|
|
break;
|
|
default: // HydrogenType::NotAHydrogen:
|
|
break;
|
|
}
|
|
if (atom->hasQuery()) {
|
|
if (atom->getQuery()->getDescription() == "RecursiveStructure") {
|
|
auto *rsq = dynamic_cast<RecursiveStructureQuery *>(atom->getQuery());
|
|
CHECK_INVARIANT(rsq, "could not convert recursive structure query");
|
|
auto res = hasQueryHs(*rsq->getQueryMol());
|
|
if (res.second) { // unmergableH implies queryH
|
|
return res;
|
|
}
|
|
queryHs |= res.first;
|
|
}
|
|
|
|
// FIX: shouldn't be repeating this code here -- yet again!
|
|
std::list<QueryAtom::QUERYATOM_QUERY::CHILD_TYPE> childStack(
|
|
atom->getQuery()->beginChildren(), atom->getQuery()->endChildren());
|
|
while (!childStack.empty()) {
|
|
QueryAtom::QUERYATOM_QUERY::CHILD_TYPE qry = childStack.front();
|
|
childStack.pop_front();
|
|
if (qry->getDescription() == "RecursiveStructure") {
|
|
auto *rsq = dynamic_cast<RecursiveStructureQuery *>(qry.get());
|
|
CHECK_INVARIANT(rsq, "could not convert recursive structure query");
|
|
auto res = hasQueryHs(*rsq->getQueryMol());
|
|
if (res.second) {
|
|
return res;
|
|
}
|
|
queryHs |= res.first;
|
|
} else {
|
|
childStack.insert(childStack.end(), qry->beginChildren(),
|
|
qry->endChildren());
|
|
}
|
|
}
|
|
}
|
|
} // end of recursion loop
|
|
|
|
return std::make_pair(queryHs, false);
|
|
}
|
|
|
|
} // namespace MolOps
|
|
} // namespace RDKit
|