pub_soft/rdkit
3
0
Fork 0
mirror of https://github.com/rdkit/rdkit.git synced 2026-06-04 21:54:27 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
Release_2025_09
rdkit/Code/GraphMol/MolInteractionFields/MIFDescriptors.cpp
Ricardo Rodriguez e3f2e4651e Refactor iostreams includes (#8846) 
* refactor iostreams includes

* restore ostream to MonomerInfo.cpp
2025-10-31 13:42:16 +01:00

1983 lines
72 KiB
C++
Raw Permalink Blame History

//
// Copyright (c) 2014-2024, Novartis Institutes for BioMedical Research 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 "MIFDescriptors.h"
#include <Geometry/point.h>
#include <Geometry/UniformRealValueGrid3D.h>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/Substruct/SubstructMatch.h>
#include <RDGeneral/FileParseException.h>
#include <RDGeneral/BadFileException.h>
#include <GraphMol/ForceFieldHelpers/MMFF/AtomTyper.h>
#include <ForceField/MMFF/Nonbonded.h>
#include <ForceField/UFF/Nonbonded.h>
#include <ForceField/UFF/Params.h>
#include <vector>
#include <fstream>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923
#endif
namespace {
constexpr double M_55D = 0.95993108859688126730;
// angle C-O-H , C-N-H
constexpr double M_70_5D = 1.23045712265600235173;
constexpr double M_110D = 1.91986217719376253461;
constexpr double CUTOFF = 0.001;
constexpr double MIN_CUTOFF_VAL = CUTOFF * CUTOFF;
} // namespace
namespace RDMIF {
std::unique_ptr<RDGeom::UniformRealValueGrid3D> constructGrid(
const RDKit::ROMol &mol, int confId, double margin, double spacing) {
PRECONDITION(mol.getNumConformers(), "No conformers available for molecule.");
const std::vector<RDGeom::Point3D> &ptVect =
mol.getConformer(confId).getPositions();
RDGeom::Point3D originPt(0.0, 0.0, 0.0);
RDGeom::Point3D marginPt(margin, margin, margin);
const RDGeom::Point3D &firstPoint =
(!ptVect.empty() ? ptVect.front() : originPt);
auto minPt = firstPoint;
auto maxPt = minPt;
for (const auto &pt : ptVect) {
for (auto i = 0u; i < pt.dimension(); ++i) {
minPt[i] = std::min(minPt[i], pt[i]);
maxPt[i] = std::max(maxPt[i], pt[i]);
}
}
minPt -= marginPt;
maxPt += marginPt;
auto res = std::make_unique<RDGeom::UniformRealValueGrid3D>(
maxPt.x - minPt.x, maxPt.y - minPt.y, maxPt.z - minPt.z, spacing, &minPt);
return res;
}
VdWaals::VdWaals(const RDKit::ROMol &mol, int confId, double cutoff) {
d_cutoff = std::max(cutoff * cutoff, MIN_CUTOFF_VAL);
d_nAtoms = mol.getNumAtoms();
d_pos.reserve(3 * d_nAtoms);
d_R_star_ij.reserve(d_nAtoms);
d_wellDepth.reserve(d_nAtoms);
// this will throw a ConformerException if confId does not exist
const RDKit::Conformer &conf = mol.getConformer(confId);
d_mol.reset(new RDKit::ROMol(mol, false, conf.getId()));
}
void VdWaals::fillVectors() {
const auto &conf = d_mol->getConformer();
for (unsigned int i = 0; i < d_nAtoms; i++) {
const RDGeom::Point3D &pt = conf.getAtomPos(i);
d_pos.push_back(pt.x);
d_pos.push_back(pt.y);
d_pos.push_back(pt.z);
fillVdwParamVectors(i);
}
}
MMFFVdWaals::MMFFVdWaals(const RDKit::ROMol &mol, int confId,
unsigned int probeAtomType, bool scaling,
double cutoff)
: VdWaals::VdWaals(mol, confId, cutoff), d_scaling(scaling) {
d_props.reset(new RDKit::MMFF::MMFFMolProperties(*d_mol));
if (!d_props->isValid()) {
throw ValueErrorException(
"No MMFF atom types available for at least one atom in molecule.");
}
d_mmffVdW = RDKit::MMFF::DefaultParameters::getMMFFVdW();
d_probeParams = (*d_mmffVdW)(probeAtomType);
fillVectors();
}
void MMFFVdWaals::fillVdwParamVectors(unsigned int atomIdx) {
PRECONDITION(atomIdx < d_mol->getNumAtoms(), "atomIdx out of bounds");
const auto iAtomType = d_props->getMMFFAtomType(atomIdx);
auto params = (*d_mmffVdW)(iAtomType);
auto rStarIJ = ForceFields::MMFF::Utils::calcUnscaledVdWMinimum(
d_mmffVdW, params, d_probeParams);
d_R_star_ij.push_back(rStarIJ);
d_wellDepth.push_back(ForceFields::MMFF::Utils::calcUnscaledVdWWellDepth(
rStarIJ, params, d_probeParams));
// scaling for taking undirected H-Bonds into account
if (d_scaling) {
ForceFields::MMFF::Utils::scaleVdWParams(d_R_star_ij[atomIdx],
d_wellDepth[atomIdx], d_mmffVdW,
params, d_probeParams);
}
}
UFFVdWaals::UFFVdWaals(const RDKit::ROMol &mol, int confId,
const std::string &probeAtomType, double cutoff)
: VdWaals::VdWaals(mol, confId, cutoff) {
d_uffParamColl = ForceFields::UFF::ParamCollection::getParams();
d_probeParams = (*d_uffParamColl)(probeAtomType);
const auto [params, haveParams] = RDKit::UFF::getAtomTypes(mol);
if (!haveParams) {
throw ValueErrorException(
"No UFF atom types available for at least one atom in molecule.");
}
d_params = std::move(params);
fillVectors();
}
void UFFVdWaals::fillVdwParamVectors(unsigned int atomIdx) {
PRECONDITION(atomIdx < d_mol->getNumAtoms(), "atomIdx out of bounds");
d_R_star_ij.push_back(d_probeParams->x1 * d_params[atomIdx]->x1);
d_wellDepth.push_back(ForceFields::UFF::Utils::calcNonbondedDepth(
d_probeParams, d_params[atomIdx]));
}
double VdWaals::operator()(double x, double y, double z, double thres) const {
double res = 0.0;
for (unsigned int i = 0, j = 0; i < d_nAtoms; ++i) {
auto temp = x - d_pos[j++];
auto dist2 = temp * temp;
temp = y - d_pos[j++];
dist2 += temp * temp;
temp = z - d_pos[j++];
dist2 += temp * temp;
if (dist2 < thres) {
dist2 = std::max(dist2, d_cutoff);
res += calcEnergy(dist2, d_R_star_ij[i], d_wellDepth[i]);
}
}
return res;
}
double UFFVdWaals::calcEnergy(double dist2, double x_ij,
double wellDepth) const {
double r6 = x_ij / dist2;
r6 *= r6 * r6;
double r12 = r6 * r6;
return wellDepth * (r12 - 2.0 * r6);
}
double MMFFVdWaals::calcEnergy(double dist2, double R_star_ij,
double wellDepth) const {
return ForceFields::MMFF::Utils::calcVdWEnergy(sqrt(dist2), R_star_ij,
wellDepth);
}
namespace CoulombDetail {
constexpr double prefactor =
1 / (4.0 * 3.141592 * 8.854188) * 1.602 * 1.602 * 6.02214129 * 10000;
}
Coulomb::Coulomb(const std::vector<double> &charges,
const std::vector<RDGeom::Point3D> &positions,
double probeCharge, bool absVal, double alpha, double cutoff)
: d_nAtoms(charges.size()),
d_absVal(absVal),
d_cutoff(cutoff * cutoff),
d_probe(probeCharge),
d_alpha(alpha),
d_charges(charges) {
PRECONDITION(d_charges.size() == positions.size(),
"Lengths of positions and charges vectors do not match.");
d_pos.reserve(3 * d_nAtoms);
for (const auto &position : positions) {
d_pos.push_back(position.x);
d_pos.push_back(position.y);
d_pos.push_back(position.z);
}
if (fabs(d_alpha) < MIN_CUTOFF_VAL) {
d_softcore = false;
if (d_cutoff < MIN_CUTOFF_VAL) {
d_cutoff = CUTOFF;
}
} else {
d_softcore = true;
}
}
Coulomb::Coulomb(const RDKit::ROMol &mol, int confId, double probeCharge,
bool absVal, const std::string &prop, double alpha,
double cutoff)
: d_nAtoms(mol.getNumAtoms()),
d_absVal(absVal),
d_cutoff(cutoff * cutoff),
d_probe(probeCharge),
d_alpha(alpha) {
d_charges.reserve(d_nAtoms);
d_pos.reserve(3 * d_nAtoms);
RDKit::Conformer conf = mol.getConformer(confId);
for (unsigned int i = 0; i < d_nAtoms; ++i) {
d_charges.push_back(mol.getAtomWithIdx(i)->getProp<double>(prop));
const RDGeom::Point3D &pt = conf.getAtomPos(i);
d_pos.push_back(pt.x);
d_pos.push_back(pt.y);
d_pos.push_back(pt.z);
}
if (fabs(d_alpha) < MIN_CUTOFF_VAL) {
d_softcore = false;
if (d_cutoff < MIN_CUTOFF_VAL) {
d_cutoff = CUTOFF;
}
} else {
d_softcore = true;
}
}
double Coulomb::operator()(double x, double y, double z, double thres) const {
double res = 0.0, dist2, temp;
if (d_softcore) {
for (unsigned int i = 0, j = 0; i < d_nAtoms; i++) {
temp = x - d_pos[j++];
dist2 = temp * temp;
temp = y - d_pos[j++];
dist2 += temp * temp;
temp = z - d_pos[j++];
dist2 += temp * temp;
if (dist2 < thres) {
res += d_charges[i] * (1.0 / sqrt(d_alpha + dist2));
}
}
} else {
for (unsigned int i = 0, j = 0; i < d_nAtoms; i++) {
temp = x - d_pos[j++];
dist2 = temp * temp;
temp = y - d_pos[j++];
dist2 += temp * temp;
temp = z - d_pos[j++];
dist2 += temp * temp;
if (dist2 < thres) {
dist2 = std::max(dist2, d_cutoff);
res += d_charges[i] * (1.0 / sqrt(dist2));
}
}
}
res *= CoulombDetail::prefactor * d_probe;
if (d_absVal) {
res = -fabs(
res); // takes the negative absolute value of the interaction energy
}
return res;
}
CoulombDielectric::CoulombDielectric(
const std::vector<double> &charges,
const std::vector<RDGeom::Point3D> &positions, double probeCharge,
bool absVal, double alpha, double cutoff, double epsilon, double xi)
: d_nAtoms(charges.size()),
d_absVal(absVal),
d_cutoff(cutoff * cutoff),
d_probe(probeCharge),
d_epsilon(epsilon),
d_xi(xi),
d_alpha(alpha),
d_charges(charges) {
PRECONDITION(d_charges.size() == positions.size(),
"Lengths of position and charge vectors do not match.");
d_dielectric = (d_xi - d_epsilon) / (d_xi + d_epsilon);
std::vector<unsigned int> neighbors(positions.size(), 0);
d_dists.resize(d_nAtoms);
d_sp.reserve(d_nAtoms);
d_pos.reserve(3 * d_nAtoms);
for (unsigned int i = 0; i < positions.size(); ++i) {
d_pos.push_back(positions[i].x);
d_pos.push_back(positions[i].y);
d_pos.push_back(positions[i].z);
for (unsigned int j = i + 1; j < positions.size(); ++j) {
double dis = (positions[j] - positions[i]).length();
if (dis < 4.0) {
++neighbors[i];
++neighbors[j];
}
}
}
for (unsigned int neighbor : neighbors) {
switch (neighbor) {
case 0:
case 1:
case 2:
case 3:
case 4:
case 5:
case 6:
d_sp.push_back(0.0);
break;
case 7:
d_sp.push_back(.4);
break;
case 8:
d_sp.push_back(.9);
break;
case 9:
d_sp.push_back(1.4);
break;
case 10:
d_sp.push_back(1.9);
break;
case 11:
d_sp.push_back(2.6);
break;
default:
d_sp.push_back(4.0);
}
}
if (fabs(d_alpha) < MIN_CUTOFF_VAL) {
d_softcore = false;
if (d_cutoff < MIN_CUTOFF_VAL * MIN_CUTOFF_VAL) {
d_cutoff = CUTOFF * CUTOFF;
}
} else {
d_softcore = true;
}
}
CoulombDielectric::CoulombDielectric(const RDKit::ROMol &mol, int confId,
double probeCharge, bool absVal,
const std::string &prop, double alpha,
double cutoff, double epsilon, double xi)
: d_nAtoms(mol.getNumAtoms()),
d_absVal(absVal),
d_cutoff(cutoff * cutoff),
d_probe(probeCharge),
d_epsilon(epsilon),
d_xi(xi),
d_alpha(alpha) {
PRECONDITION(mol.getNumConformers() > 0, "No Conformers for Molecule");
d_charges.reserve(d_nAtoms);
d_sp.reserve(d_nAtoms);
d_dists.resize(d_nAtoms);
d_pos.reserve(3 * d_nAtoms);
RDKit::Conformer conf = mol.getConformer(confId);
for (unsigned int i = 0; i < d_nAtoms; ++i) {
d_charges.push_back(mol.getAtomWithIdx(i)->getProp<double>(prop));
const RDGeom::Point3D &pt = conf.getAtomPos(i);
d_pos.push_back(pt.x);
d_pos.push_back(pt.y);
d_pos.push_back(pt.z);
}
d_dielectric = (d_xi - d_epsilon) / (d_xi + d_epsilon);
std::vector<unsigned int> neighbors(d_charges.size(), 0);
for (unsigned int i = 0; i < d_charges.size() - 1; ++i) {
for (unsigned int j = i + 1; j < d_charges.size(); ++j) {
double temp = d_pos[j * 3] - d_pos[i * 3];
double dist2 = temp * temp;
temp = d_pos[j * 3 + 1] - d_pos[i * 3 + 1];
dist2 += temp * temp;
temp = d_pos[j * 3 + 2] - d_pos[i * 3 + 2];
dist2 += temp * temp;
if (dist2 < 16.0) {
neighbors[i]++;
neighbors[j]++;
}
}
}
for (unsigned int neighbor : neighbors) {
switch (neighbor) {
case 0:
case 1:
case 2:
case 3:
case 4:
case 5:
case 6:
d_sp.push_back(0.0);
break;
case 7:
d_sp.push_back(.4);
break;
case 8:
d_sp.push_back(.9);
break;
case 9:
d_sp.push_back(1.4);
break;
case 10:
d_sp.push_back(1.9);
break;
case 11:
d_sp.push_back(2.6);
break;
default:
d_sp.push_back(4.0);
}
}
if (fabs(d_alpha) < MIN_CUTOFF_VAL) {
d_softcore = false;
if (d_cutoff < MIN_CUTOFF_VAL * MIN_CUTOFF_VAL) {
d_cutoff = CUTOFF * CUTOFF;
}
} else {
d_softcore = true;
}
}
double CoulombDielectric::operator()(double x, double y, double z,
double thres) const {
int neigh = 0;
double res = 0.0, sq = 0.0;
for (unsigned int i = 0, j = 0; i < d_nAtoms; ++i, j += 3) {
double temp = x - d_pos[j];
double dist2 = temp * temp;
temp = y - d_pos[j + 1];
dist2 += temp * temp;
temp = z - d_pos[j + 2];
dist2 += temp * temp;
d_dists[i] = dist2;
if (dist2 < 16.0) {
neigh += 1;
}
}
switch (neigh) {
case 0:
case 1:
case 2:
case 3:
case 4:
case 5:
case 6:
sq = 0.0;
break;
case 7:
sq = .4;
break;
case 8:
sq = .9;
break;
case 9:
sq = 1.4;
break;
case 10:
sq = 1.9;
break;
case 11:
sq = 2.6;
break;
default:
sq = 4.0;
}
if (d_softcore) {
for (unsigned int i = 0; i < d_nAtoms; i++) {
if (d_dists[i] < thres) {
res += d_charges[i] *
(1 / sqrt(d_alpha + d_dists[i]) +
d_dielectric / sqrt(d_alpha + d_dists[i] + 4.0 * d_sp[i] * sq));
}
}
} else {
for (unsigned int i = 0; i < d_nAtoms; i++) {
if (d_dists[i] < thres) {
double dist2 = std::max(d_dists[i], d_cutoff);
res += d_charges[i] * (1 / sqrt(dist2) +
d_dielectric / sqrt(dist2 + 4.0 * d_sp[i] * sq));
}
}
}
res *= CoulombDetail::prefactor * (1 / d_xi) * d_probe;
if (d_absVal) {
res = -fabs(
res); // takes the negative absolute value of the interaction energy
}
return res;
}
namespace HBondDetail {
const double em[2][2] = {{-8.368, -11.715}, {-11.715, -16.736}};
const double rm[2][2] = {{3.2, 3.0}, {3.0, 2.8}};
const double K1 = 562.25380293;
const double K2 = 0.11697778;
const double bondlength[2] = {-0.972, -1.019};
double cos_2(double t, double, double) {
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
return temp;
} else {
return 0.0;
}
};
double cos_2_0(double t, double t_0 = 0.0, double t_i = 1.0) {
double temp;
if (t_i < M_55D) {
temp = cos(t_0);
temp *= temp;
return temp;
} else if (t_i < M_PI_2) {
temp = cos(t);
temp *= temp;
return temp;
} else {
return 0.0;
}
};
double cos_2_rot(double t, double, double) {
t -= M_70_5D;
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
return temp;
} else {
return 0.0;
}
};
double cos_4(double t, double, double) {
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
temp *= temp;
return temp;
} else {
return 0.0;
}
};
double cos_4_rot(double t, double, double) {
t -= M_70_5D;
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
temp *= temp;
return temp;
} else {
return 0.0;
}
};
double cos_6(double t, double, double) {
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
temp *= temp * temp;
return temp;
} else {
return 0.0;
}
};
double cos_6_rot(double t, double, double) {
t -= M_70_5D;
if (t < M_PI_2) {
double temp = cos(t);
temp *= temp;
temp *= temp * temp;
return temp;
} else {
return 0.0;
}
};
double cos_acc(double, double t_0, double t_i) {
double temp;
if (t_i < M_PI_2) {
temp = cos(t_0) * (0.9 + 0.1 * sin(2 * t_i));
return temp;
} else if (t_i < M_110D) {
temp = cos(t_i);
temp *= temp;
temp = K2 - temp;
temp *= temp * temp;
temp *= cos(t_0) * K1;
return temp;
} else {
return 0.0;
}
};
double no_dep(double, double, double) { return 1.0; };
} // namespace HBondDetail
HBond::HBond(const RDKit::ROMol &mol, int confId,
const std::string &probeAtomType, bool fixed, double cutoff)
: d_cutoff(cutoff * cutoff) {
if (d_cutoff < (MIN_CUTOFF_VAL * MIN_CUTOFF_VAL)) {
d_cutoff = CUTOFF * CUTOFF;
}
if (probeAtomType == "O") {
d_DAprop = 'D';
d_probetype = O;
} else if (probeAtomType == "OH") {
d_DAprop = 'A';
d_probetype = O;
} else if (probeAtomType == "N") {
d_DAprop = 'D';
d_probetype = N;
} else if (probeAtomType == "NH") {
d_DAprop = 'A';
d_probetype = N;
} else {
const std::string msg = "Probe atom type not supported: " + probeAtomType;
BOOST_LOG(rdErrorLog) << msg << std::endl;
throw ValueErrorException(msg);
}
d_nInteract = mol.getNumAtoms(); // number of atoms = highest possible number
// of interactions
std::vector<unsigned int> specialAtoms;
findSpecials(mol, confId, fixed, specialAtoms);
if (d_DAprop == 'A') {
if (fixed) {
findAcceptors(mol, confId, specialAtoms);
} else {
findAcceptorsUnfixed(mol, confId, specialAtoms);
}
} else if (d_DAprop == 'D') {
if (fixed) {
findDonors(mol, confId, specialAtoms);
} else {
findDonorsUnfixed(mol, confId, specialAtoms);
}
} else { // this should never be the case
BOOST_LOG(rdErrorLog) << "HBond: unknown target property d_DAprop: "
<< d_DAprop << std::endl;
}
d_nInteract = d_targettypes.size(); // updated to number of interactions
d_eneContrib.resize(d_nInteract, 0);
d_vectTargetProbe.resize(d_nInteract * 3, 0);
POSTCONDITION(
d_nInteract * 3 == d_pos.size(),
"Error in constructing H-Bond descriptor: Vector length mismatch (target atom types).");
POSTCONDITION(
d_nInteract * 3 == d_direction.size(),
"Error in constructing H-Bond descriptor: Vector length mismatch (bond directions).");
POSTCONDITION(
d_nInteract == d_function.size(),
"Error in constructing H-Bond descriptor: Vector length mismatch (angular functions).");
POSTCONDITION(
d_nInteract * 3 == d_plane.size(),
"Error in constructing H-Bond descriptor: Vector length mismatch (lone pair planes).");
POSTCONDITION(
d_nInteract == d_lengths.size(),
"Error in constructing H-Bond descriptor: Vector length mismatch (bondlengths).");
}
unsigned int HBond::findSpecials(const RDKit::ROMol &mol, int confId,
bool fixed,
std::vector<unsigned int> &specials) {
using namespace HBondDetail;
RDKit::MatchVectType matches;
RDKit::ROMol::ADJ_ITER nbrIdx, endNbrs;
RDGeom::Point3D pos, dir, hbonddir, plane, bondDirection[12];
unsigned int nbrs;
unsigned int match = 0, nMatches = 0;
const RDKit::Conformer &conf = mol.getConformer(confId);
// RDKit::RWMol thr = *RDKit::SmilesToMol("C[C@H]([C@@H](C(=O))N)O");
// //threonine, serine itself is a substructure of this
static const auto ser = RDKit::v2::SmilesParse::MolFromSmiles(
"[CH2]([C@@H](C(=O))N)O"); // serine
// RDKit::RWMol his = *RDKit::SmilesToMol("Cc1cnc[nH]1"); //imidazole residue,
// is correctly taken into account in 'normal' treatment
match = RDKit::SubstructMatch(mol, *ser, matches);
nMatches += match;
for (auto &matche : matches) {
const auto *atom = mol.getAtomWithIdx(matche.second);
if (atom->getAtomicNum() == 8) {
boost::tie(nbrIdx, endNbrs) = mol.getAtomNeighbors(atom);
pos = conf.getAtomPos(matche.second);
if (d_DAprop == 'A') {
nbrs = 0;
while (nbrIdx != endNbrs) { // loop over atoms
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() == 1) {
hbonddir = pos - conf.getAtomPos(*nbrIdx);
hbonddir.normalize();
} else {
bondDirection[nbrs] = pos - conf.getAtomPos(*nbrIdx);
bondDirection[nbrs].normalize();
++nbrs;
}
++nbrIdx;
}
if (nbrs == 2) {
dir = bondDirection[0] + hbonddir;
plane =
bondDirection[0].crossProduct(hbonddir); // X-O-Y plane vector
plane = plane.crossProduct(dir); // lp plane vector
if (fixed) {
addVectElements(O, &cos_2_0, pos, dir, plane);
} else {
addVectElements(O, &cos_4_rot, pos, bondDirection[0]);
}
specials.push_back(matche.second);
}
} else if (d_DAprop == 'D') {
if (fixed) {
while (nbrIdx != endNbrs) {
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // hydrogen neighbor necessary for H bond donation
dir = conf.getAtomPos(*nbrIdx) - pos;
addVectElements(O, &cos_6, pos, dir);
}
++nbrIdx;
}
} else {
while (nbrIdx != endNbrs) {
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // search for non-hydrogen bond direction
dir = pos - conf.getAtomPos(*nbrIdx);
addVectElements(O, &cos_6_rot, pos, dir);
}
++nbrIdx;
}
}
specials.push_back(matche.second);
} else { // this should never be the case
BOOST_LOG(rdErrorLog)
<< "HBond::operator(): unknown target property d_DAprop: "
<< d_DAprop << std::endl;
}
}
}
return nMatches;
}
/* General structure of findAcceptors, findAcceptorsUnfixed, findDonors,
* findDonorsUnfixed functions: loop over all atoms:
* - check whether atom was already treated in specialAtoms
* - switch ( atomicNum ): find atoms which are able to donate/accept
* hydrogen bonds if able:
* - switch ( number of neighbors ): find the right geometry
* - check for charge or aromatic neighborhood
* - calculate adequate hydrogen bond direction vector,
* lone pair plane vector
* - choose correct angular function
* - add atomtype, angular function, position of atom,
* hydrogen bond direction vector and lone pair plane vector to vectors for
* calculating the interaction returns number of added interactions
*/
unsigned int HBond::findAcceptors(const RDKit::ROMol &mol, int confId,
const std::vector<unsigned int> &specials) {
using namespace HBondDetail;
const RDKit::Conformer &conf =
mol.getConformer(confId); // get conformer of molecule
RDKit::ROMol::ADJ_ITER nbrIdx, endNbrs;
RDKit::ROMol::ADJ_ITER secnbrIdx, secendNbrs;
RDGeom::Point3D pos, dir, plane, bondDirection[12];
unsigned int nbrs; // no of neigbors
unsigned int interact = 0; // no of interactions
bool aromaticnbr; // aromatic neighbor atom?
for (unsigned int i = 0; i < d_nInteract; i++) { // loop over all atoms
if (std::find(specials.begin(), specials.end(), i) !=
specials.end()) { // check whether atom was already treated specially
interact++;
continue;
}
const auto *atom = mol.getAtomWithIdx(i); // get ptr to atom
switch (atom->getAtomicNum()) { // find atoms able to donate hydrogen bonds
case 7: // Nitrogen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get position of atom
nbrs = 0;
while (nbrIdx != endNbrs) { // loop over all neigbors
bondDirection[nbrs] =
pos - conf.getAtomPos(*nbrIdx); // store bond vectors in an array
bondDirection[nbrs].normalize();
++nbrs; // count neigbors
++nbrIdx;
}
switch (nbrs) { // number of neigbors
case 1: // sp, eg. nitriles
if (atom->getFormalCharge() <=
0) { // no positively charged nitrogens
addVectElements(
N, &cos_2, pos,
bondDirection[0]); // no differentiation between in-plane and
// out-of-plane angle, plane not needed
interact++;
}
break;
case 2: // eg. imines, heterocycles
if (atom->getFormalCharge() <=
0) { // no positively charged nitrogen
dir = bondDirection[0] +
bondDirection[1]; // get hydrogen bond direction
plane = bondDirection[0].crossProduct(
bondDirection[1]); // normal vector of lone pair plane =
// plane of three atoms
addVectElements(N, &cos_2, pos, dir, plane);
interact++;
}
break;
case 3: // amines, iminium ions
if (atom->getFormalCharge() <= 0 &&
!(RDKit::MolOps::atomHasConjugatedBond(
atom))) { // no positively charged nitrogen, no conjugated
// nitrogen (amide bonds!)
dir = bondDirection[0] + bondDirection[1] +
bondDirection[2]; // get hydrogen bond direction
addVectElements(N, &cos_2, pos,
dir); // no differentiation between in-plane and
// out-of-plane angle, plane not needed
interact++;
}
break;
case 0: // unbound nitrogen atoms, no hydrogen bonding
case 4: // ammonium ions, no hydrogen bonding
break;
default: // more than four bonds: not possible with nitrogen
BOOST_LOG(rdErrorLog)
<< "HBond: Nitrogen atom bound to more than 4 neighbor atoms: Atom: "
<< i << std::endl;
}
break;
case 8: // Oxygen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get position of atom
nbrs = 0;
aromaticnbr = false;
while (nbrIdx != endNbrs) { // loop over neighbors
bondDirection[nbrs] =
pos - conf.getAtomPos(*nbrIdx); // store bond vectors in an array
bondDirection[nbrs].normalize(); // normalization
if (mol.getAtomWithIdx(*nbrIdx)
->getIsAromatic()) { // check whether neighbor is aromatic
// (phenolic oxygens)
aromaticnbr = true;
}
++nbrs; // count neighbors
++nbrIdx;
}
switch (nbrs) { // no of neighbors
case 1: // carbonyl, carboxyl C=O, X=O (X=S,P,...), anions
// (alcoholates, carboxylates)
--nbrIdx;
boost::tie(secnbrIdx, secendNbrs) = mol.getAtomNeighbors(
mol.getAtomWithIdx(*nbrIdx)); // get neighbors of neighbor atom
while (secnbrIdx !=
secendNbrs) { // loop over neighbors of neighbor atom
if ((*secnbrIdx) !=
i) { // second neighbor should not be the oxygen itself
// bond direction of neighbor atom to second neighbor atom
dir = conf.getAtomPos(*secnbrIdx) - conf.getAtomPos(*nbrIdx);
break; // we only need one bond vector
}
++secnbrIdx;
}
plane = bondDirection[0].crossProduct(dir); // lp plane vector
if (atom->getFormalCharge() == 0) { // carbonyl, carboxyl C=O, X=O
addVectElements(O, &cos_acc, pos, bondDirection[0], plane);
interact++;
} else if (atom->getFormalCharge() ==
-1) { // anion, eg. phenolate, carboxylic acid
if (RDKit::MolOps::atomHasConjugatedBond(atom)) {
// charged oxygen in conjungated system, eg phenolates or
// carboxylates
addVectElements(O, &cos_acc, pos, bondDirection[0], plane);
} else {
// non-conjugated anion, eg sp3-alcoholate
addVectElements(O, &cos_2, pos, bondDirection[0], plane);
}
interact++;
}
break;
case 2: // alcohols, ethers, carboxyl OH
dir = bondDirection[0] + bondDirection[1]; // get bond direction
plane = bondDirection[0].crossProduct(
bondDirection[1]); // X-O-Y plane vector
plane = plane.crossProduct(dir); // lp plane vector
if (aromaticnbr) { // if oxygen bound to aromatic system, eg.
// phenol
addVectElements(O, &cos_2, pos, dir, plane);
} else { // all other
addVectElements(O, &cos_2_0, pos, dir, plane);
}
interact++;
break;
case 0: // oxygen atoms, no hydrogen bonding
case 3: // only with positively charged atom possible, no hydrogen
// bonding
break;
default: // more than 3 neighbors: not possible with oxygen
BOOST_LOG(rdErrorLog)
<< "HBond: Oxygen atom bound to more than 3 neighbor atoms: Atom: "
<< i << std::endl;
}
break;
// Halogens
case 9: // F
case 17: // Cl
case 35: // Br
case 53: // I
pos = conf.getAtomPos(i); // get position of atom
dir = RDGeom::Point3D(1.0, 1.0, 1.0); // no directionality needed
addVectElements(
N, &no_dep, pos,
dir); // type of halogens ~ nitrogen; no lp plane needed
interact++;
break;
default:
break;
}
}
return interact;
}
unsigned int HBond::findAcceptorsUnfixed(
const RDKit::ROMol &mol, int confId,
const std::vector<unsigned int> &specials) {
using namespace HBondDetail;
const RDKit::Conformer &conf =
mol.getConformer(confId); // get conformer of molecule
RDKit::ROMol::ADJ_ITER nbrIdx, endNbrs;
RDKit::ROMol::ADJ_ITER secnbrIdx, secendNbrs;
RDGeom::Point3D pos, dir, plane, bondDirection[12];
unsigned int nbrs,
nonhnbrs; // no of neighbors, no of nonhydrogen - neighbors
unsigned int interact = 0; // no of interactions
bool aromaticnbr; // aromatic neighbor atom?
for (unsigned int i = 0; i < d_nInteract; i++) { // loop over all atoms
if (std::find(specials.begin(), specials.end(), i) !=
specials.end()) { // check whether atom was already treated specially
interact++;
continue;
}
const auto *atom = mol.getAtomWithIdx(i); // get pointer to atom
switch (atom->getAtomicNum()) { // find atoms able to accept hydrogen bonds
// (O, N, halogens)
case 7: // Nitrogen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get position of atom
nbrs = 0;
nonhnbrs = 0;
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // hydrogen bond directions are not included
bondDirection[nonhnbrs] =
pos -
conf.getAtomPos(*nbrIdx); // store bond direction in an array
bondDirection[nonhnbrs].normalize();
++nonhnbrs; // count non hydrogen neighbors
}
++nbrs; // count neighbors
++nbrIdx;
}
switch (nbrs) { // number of neigbors
case 1: // sp, eg. nitriles, no difference to fixed bonds
if (atom->getFormalCharge() <=
0) { // no positively charged nitrogens
addVectElements(
N, &cos_2, pos,
bondDirection[0]); // no differentiation between in-plane and
// out-of-plane angle, plane not needed
interact++;
}
break;
case 2: // eg. imines, heterocycles
if (atom->getFormalCharge() <=
0) { // no positively charged nitrogen
if (nonhnbrs == 2) { // secondary imines, heterocycles
dir = bondDirection[0] +
bondDirection[1]; // get hydrogen bond direction
plane = bondDirection[0].crossProduct(
bondDirection[1]); // normal vector of lone pair plane =
// plane of three atoms
addVectElements(N, &cos_2, pos, dir, plane);
} else if (nonhnbrs == 1) { // primary imine, hydrogen is allowed
// to swap places
nbrIdx -= 2;
nbrs = nonhnbrs;
while (nbrIdx != endNbrs) {
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // only bond directions to hydrogens are included
bondDirection[nbrs] = pos - conf.getAtomPos(*nbrIdx);
bondDirection[nbrs].normalize();
++nbrs;
}
++nbrIdx;
}
plane = bondDirection[0].crossProduct(
bondDirection[1]); // normal vector of lone pair plane =
// plane of three atoms
addVectElements(N, &cos_acc, pos, bondDirection[0], plane);
} else { //[N-]H2 rotating
addVectElements(N, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
}
interact++;
}
break;
case 3: // amines, iminium ions
if (atom->getFormalCharge() <=
0) { // no iminium ions, no positively charged nitrogen
if (nonhnbrs == 0) { // ammonia
addVectElements(N, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
} else if (nonhnbrs == 1) { // primary amines, rotation
addVectElements(
N, &cos_2_rot, pos,
bondDirection[0]); // no differentiation between in-plane
// and out-of-plane angle, plane not
// needed
} else { // secondary amines (not flexible) and tertiary amines,
// same as fixed hydrogen bonds
if (atom->getFormalCharge() <= 0 &&
!(RDKit::MolOps::atomHasConjugatedBond(
atom))) { // positively charged nitrogen, no conjugated
// nitrogen (amide bonds!)
nbrIdx -= 3;
nbrs = nonhnbrs;
while (nbrIdx != endNbrs) {
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // only bond directions to hydrogens are included
bondDirection[nbrs] = pos - conf.getAtomPos(*nbrIdx);
bondDirection[nbrs].normalize();
++nbrs;
}
++nbrIdx;
}
dir = bondDirection[0] + bondDirection[1] +
bondDirection[2]; // hydrogen bond direction
addVectElements(
N, &cos_2, pos,
dir); // no differentiation between in-plane and
// out-of-plane angle, plane not needed
}
}
interact++;
}
break;
case 0: // unbound nitrogen atoms, no hydrogen bonding
case 4: // ammonium ions, no hydrogen bonding
break;
default: // more than four bonds: not possible with nitrogen
BOOST_LOG(rdErrorLog)
<< "HBond: Nitrogen atom bound to more than 4 neighbor atoms: Atom: "
<< i << std::endl;
}
break;
case 8: // Oxygen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get atom position
nbrs = 0;
nonhnbrs = 0;
aromaticnbr = false;
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // bond directions to hydrogen are not included
bondDirection[nonhnbrs] = pos - conf.getAtomPos(*nbrIdx);
bondDirection[nonhnbrs].normalize();
++nonhnbrs; // count non hydrogen neighbors
}
if (mol.getAtomWithIdx(*nbrIdx)
->getIsAromatic()) { // check whether aromatic neighbor
aromaticnbr = true;
}
++nbrs; // count neighbors
++nbrIdx;
}
switch (nbrs) { // no of neighbors
case 1: // carbonyl, carboxyl C=O, X=O (X=S,P,...), anions
// (alcoholates, carboxylates)
--nbrIdx;
boost::tie(secnbrIdx, secendNbrs) = mol.getAtomNeighbors(
mol.getAtomWithIdx(*nbrIdx)); // get neighbors of neighbor atom
while (secnbrIdx !=
secendNbrs) { // loop over neighbors of neighbor atom
if ((*secnbrIdx) !=
i) { // second neighbor should not be the oxygen itself
// bond direction of neighbor atom to second neighbor atom
dir = conf.getAtomPos(*secnbrIdx) - conf.getAtomPos(*nbrIdx);
break; // we only need one vector
}
++secnbrIdx;
}
plane = bondDirection[0].crossProduct(dir); // lp plane vector
if (atom->getFormalCharge() == 0) { // carbonyl, carboxyl C=O, X=O
addVectElements(O, &cos_acc, pos, bondDirection[0], plane);
interact++;
} else if (atom->getFormalCharge() ==
-1) { // anion, eg. phenolate, carboxylic acid
if (RDKit::MolOps::atomHasConjugatedBond(atom)) {
// charged oxygen in conjungated system, eg phenolates or
// carboxylates
addVectElements(O, &cos_acc, pos, bondDirection[0], plane);
} else { // all other
addVectElements(O, &cos_2, pos, bondDirection[0], plane);
}
interact++;
}
break;
case 2: // alcohols, ethers, carboxyl OH
if (nonhnbrs == 0) { // water
addVectElements(O, &no_dep, pos, RDGeom::Point3D(0.0, 0.0, 0.0));
} else if (nonhnbrs == 1) { // hydroxy groups
if (aromaticnbr) { // phenol
nbrIdx -= 2;
nbrs = nonhnbrs;
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // only bond directions to hydrogens are included
bondDirection[nbrs] =
pos - conf.getAtomPos(
*nbrIdx); // get O-H bond direction vector
bondDirection[nbrs].normalize();
++nbrs;
}
++nbrIdx;
}
plane = bondDirection[0].crossProduct(
bondDirection[1]); // X-O-Y plane vector
addVectElements(O, &cos_acc, pos, bondDirection[0], plane);
} else {
addVectElements(
O, &cos_2_rot, pos,
bondDirection[0]); // no plane information needed
}
} else { // ethers, same as fixed hydrogen bonds
dir = bondDirection[0] +
bondDirection[1]; // get hydrogen bond direction
plane = bondDirection[0].crossProduct(
bondDirection[1]); // X-O-Y plane vector
plane = plane.crossProduct(dir); // lp plane vector
addVectElements(O, &cos_2_0, pos, dir, plane);
}
interact++;
break;
case 0: // oxygen atoms, no hydrogen bonding
case 3: // only with positively charged atom possible, no hydrogen
// bonding
break;
default: // more than 3 neighbors: not possible with oxygen
BOOST_LOG(rdErrorLog)
<< "HBond: Oxygen atom bound to more than 3 neighbor atoms: Atom: "
<< i << std::endl;
}
break;
// Halogens
case 9: // F
case 17: // Cl
case 35: // Br
case 53: // I
pos = conf.getAtomPos(i); // get atoms position
dir = RDGeom::Point3D(1.0, 1.0, 1.0); // no directionality needed
addVectElements(
N, &no_dep, pos,
dir); // type of halogens ~ nitrogen; no lp plane needed
interact++;
break;
default:
break;
}
}
return interact;
}
unsigned int HBond::findDonors(const RDKit::ROMol &mol, int confId,
const std::vector<unsigned int> &specials) {
using namespace HBondDetail;
const RDKit::Conformer &conf =
mol.getConformer(confId); // get conformer of molecule
RDKit::ROMol::ADJ_ITER nbrIdx, endNbrs;
RDGeom::Point3D pos, dir;
unsigned int interact = 0;
for (unsigned int i = 0; i < d_nInteract; i++) { // loop over all atoms
if (std::find(specials.begin(), specials.end(), i) !=
specials.end()) { // check whether atom was already treated specially
interact++;
continue;
}
const auto *atom = mol.getAtomWithIdx(i); // get ptr to atom
switch (
atom->getAtomicNum()) { // find atoms able to donate hydrogen bonds (O,
// N, of course only with attached hydrogen)
case 7: // Nitrogen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neigbors
pos = conf.getAtomPos(i); // get position
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // hydrogen neighbor necessary for H bond donation
dir = conf.getAtomPos(*nbrIdx) -
pos; // bond direction vector, IMPORTANT: no normalization,
// because operator() needs length of vector in case of
// donors
addVectElements(N, &cos_2, pos, dir);
interact++;
}
++nbrIdx;
}
break;
case 8: // Oxygen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neigbors
pos = conf.getAtomPos(i); // get position
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // hydrogen neighbor necessary for H bond donation
dir = conf.getAtomPos(*nbrIdx) -
pos; // bond direction vector, IMPORTANT: no normalization,
// because operator() needs length of vector in case of
// donors
addVectElements(O, &cos_4, pos, dir);
interact++;
}
++nbrIdx;
}
break;
default:
break;
}
}
return interact;
}
unsigned int HBond::findDonorsUnfixed(
const RDKit::ROMol &mol, int confId,
const std::vector<unsigned int> &specials) {
using namespace HBondDetail;
const auto &conf = mol.getConformer(confId); // get conformer of molecule
RDKit::ROMol::ADJ_ITER nbrIdx, endNbrs;
RDGeom::Point3D pos, hbonddir, dir, plane;
unsigned int nbrs, nonhnbrs; // no of neighbors, no of non hydrogen neighbors
unsigned int interact = 0;
bool aromaticnbr;
for (unsigned int i = 0; i < d_nInteract; i++) { // loop over all atoms
if (std::find(specials.begin(), specials.end(), i) !=
specials.end()) { // check whether atom was already treated specially
interact++;
continue; // skip loop for this atom
}
const auto *atom = mol.getAtomWithIdx(i); // get ptr to atom
switch (
atom->getAtomicNum()) { // find atoms able to donate hydrogen bonds (O,
// N, of course only with attached hydrogen)
case 7: // Nitrogen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get position
nbrs = 0;
nonhnbrs = 0;
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() != 1) {
++nonhnbrs; // count non-hydrogen neighbors
}
++nbrs; // count neighbors
++nbrIdx;
}
nbrIdx -= nbrs;
if (nonhnbrs != nbrs) { // otherwise no hydrogens, no donation possible
switch (nbrs) { // number of neigbors
case 2: // eg. imines, heterocycles
if (nonhnbrs == 0) { //[N-]H2
addVectElements(N, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
} else { // primary imine, swapping of hydrogen is allowed
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() != 1) {
dir = pos - conf.getAtomPos(*nbrIdx); // the other bond
dir.normalize();
} else {
hbonddir = conf.getAtomPos(*nbrIdx) - pos; // hydrogen bond
}
++nbrIdx;
}
addVectElements(N, &cos_2, pos, hbonddir); // first hbond
interact++;
// let's swap hydrogen and lp:
plane = dir.crossProduct(hbonddir); // lp plane vector
plane =
plane.crossProduct(dir); // plane through other bond
// perpendicular to X-N-H plane
dir = plane *
hbonddir.dotProduct(plane); // projection of vector
// hbonddir on plane vector
hbonddir -= dir * 2; // mirroring of dir vector on plane
addVectElements(
N, &cos_2, pos,
hbonddir); // second hbond, hydrogen at other place
interact++;
}
break;
case 3: // amines, iminium ions
if (nonhnbrs ==
2) { // sec amines, no rotation, same as fixed bonds
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // hydrogen neighbor necessary for H bond donation
hbonddir =
conf.getAtomPos(*nbrIdx) - pos; // hbond direction
}
++nbrIdx;
}
addVectElements(N, &cos_2, pos, hbonddir);
interact++;
} else if (nonhnbrs == 1) { // primary amines, rotation
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // search for X-N bond
dir = pos - conf.getAtomPos(*nbrIdx); // X-N bond direction
dir.normalize();
}
++nbrIdx;
}
addVectElements(
N, &cos_2_rot, pos,
dir * (-bondlength[N])); // vector dir has typcial N-H
// bondlength, for approximate /
// average calculation of angle p
// (see operator())
interact++;
} else { // ammonia
addVectElements(N, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
}
break;
case 4:
if (nonhnbrs == 0) { // ammonium ion
addVectElements(N, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
} else if (nonhnbrs == 1) { // primary ammonium, rotation
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // search for X-N bond
dir = pos - conf.getAtomPos(*nbrIdx); // X-N bond direction
}
++nbrIdx;
}
addVectElements(
N, &cos_2_rot, pos,
dir * (-bondlength[N])); // vector dir has typcial N-H
// bondlength, for approximate /
// average calculation of angle p
// (see operator())
interact++;
} else { // secondary or tertiary ammonium, no rotation, same as
// fixed bonds
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // hydrogen neighbor necessary for H bond donation
dir = pos -
conf.getAtomPos(*nbrIdx); // hydrogen bond direction
addVectElements(N, &cos_2_rot, pos, dir);
interact++;
}
++nbrIdx;
}
}
break;
default: // more than four bonds: not possible with nitrogen
BOOST_LOG(rdErrorLog)
<< "HBond: Nitrogen atom bound to more than 4 neighbor atoms: Atom: "
<< i << std::endl;
}
}
break;
case 8: // Oxygen
boost::tie(nbrIdx, endNbrs) =
mol.getAtomNeighbors(atom); // get neighbors
pos = conf.getAtomPos(i); // get position
nbrs = 0;
nonhnbrs = 0;
aromaticnbr = false;
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() != 1) {
++nonhnbrs; // count non-hydrogen neighbors
}
if (mol.getAtomWithIdx(*nbrIdx)
->getIsAromatic()) { // check whether oxygen is bound to
// aromatic system
aromaticnbr = true;
}
++nbrs; // count neighbors
++nbrIdx;
}
nbrIdx -= nbrs;
if (nonhnbrs != nbrs) { // otherwise no hydrogen, no hydrogen bond
// donation possible
switch (nbrs) { // no of neighbors
case 1: // hydroxyl
addVectElements(O, &no_dep, pos, RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
break;
case 2:
if (nonhnbrs == 0) { // water
addVectElements(O, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
} else { // OH groups
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // search for X-O bond
dir = pos - conf.getAtomPos(*nbrIdx); // O-X bond direction
dir.normalize();
} else {
hbonddir =
conf.getAtomPos(*nbrIdx) - pos; // O-H bond direction
}
++nbrIdx;
}
if (aromaticnbr) { // phenolic oxygen, allows h to swap places
// in aromatic plane
addVectElements(O, &cos_4, pos, hbonddir); // first hbond
interact++;
// let's swap hydrogen and lp:
plane = dir.crossProduct(hbonddir); // lp plane vector
plane =
plane.crossProduct(dir); // plane through other bond
// perpendicular to X-N-H plane
dir = plane *
hbonddir.dotProduct(
plane); // projection of vector dir on plane vector
hbonddir -= dir * 2; // mirroring of dir vector on plane
addVectElements(
O, &cos_4, pos,
hbonddir); // second hbond, hydrogen at other place
interact++;
} else { // all other oxygens, flexible
addVectElements(
O, &cos_4_rot, pos,
dir * (-bondlength[O])); // vector dir has typcial O-H
// bondlength, for approximate /
// average calculation of angle
// p (see operator())
interact++;
}
}
break;
case 3: // positively charged oxygen
if (nonhnbrs == 0) { // oxonium
addVectElements(O, &no_dep, pos,
RDGeom::Point3D(0.0, 0.0, 0.0));
interact++;
} else if (nonhnbrs == 1) { // R[O+]H2
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() !=
1) { // search for X-O bond
dir = pos - conf.getAtomPos(*nbrIdx); // X-O bond direction
dir.normalize();
addVectElements(
O, &cos_4_rot, pos,
dir * (-bondlength[O])); // vector dir has typcial O-H
// bondlength, for approximate
// / average calculation of
// angle p (see operator())
interact++;
}
++nbrIdx;
}
} else { // R1R2[O+]H, non-flexible
while (nbrIdx != endNbrs) { // loop over neighbors
if (mol.getAtomWithIdx(*nbrIdx)->getAtomicNum() ==
1) { // search for O-H bond
dir = conf.getAtomPos(*nbrIdx) - pos; // O-H bond direction
addVectElements(O, &cos_4, pos, dir);
interact++;
}
++nbrIdx;
}
}
break;
case 0: // oxygen atom gas, no hydrogen bonding
break;
default:
BOOST_LOG(rdErrorLog)
<< "HBond: Oxygen atom bound to more than 3 neighbor atoms: Atom: "
<< i << std::endl;
}
}
break;
default:
break;
}
}
return interact;
}
double HBond::operator()(double x, double y, double z, double thres) const {
using namespace HBondDetail;
if (d_nInteract < 1) { // no interactions
return 0.0; // return 0.0
}
double res = 0.0;
if (d_DAprop == 'A') {
unsigned int minId = 0;
double minEne = std::numeric_limits<double>::max(); // minimal energy
double probeDirection[3]; // direction of probe
for (unsigned int i = 0, j = 0; i < d_nInteract;
i++, j += 3) { // calculation of energy contributions and searching
// for the favored probe direction ( direction of
// lowest energy contribution )
d_vectTargetProbe[j] = x - d_pos[j]; // vector of interaction
d_vectTargetProbe[j + 1] = y - d_pos[j + 1];
d_vectTargetProbe[j + 2] = z - d_pos[j + 2];
double dist2 =
d_vectTargetProbe[j] * d_vectTargetProbe[j] +
d_vectTargetProbe[j + 1] * d_vectTargetProbe[j + 1] +
d_vectTargetProbe[j + 2] *
d_vectTargetProbe[j +
2]; // calc of squared length of interaction
// std::cout << dist2 << std::endl;
if (dist2 < thres) {
double dis = sqrt(dist2);
double distN[3] = {d_vectTargetProbe[j] / dis,
d_vectTargetProbe[j + 1] / dis,
d_vectTargetProbe[j + 2] / dis};
dist2 = std::max(dist2, d_cutoff);
double t = angle(
distN[0], distN[1], distN[2], d_direction[j], d_direction[j + 1],
d_direction[j + 2]); // calc of angle between direction of hbond
// and target-probe direction
double eneTerm1 = rm[d_probetype][d_targettypes[i]];
eneTerm1 *= eneTerm1; // squared rm
eneTerm1 /= dist2; // division by squared distance
double eneTerm2 = eneTerm1 * eneTerm1; // fourth power of rm/distance
eneTerm1 = eneTerm2 * eneTerm1; // sixth power of rm/distance
eneTerm2 *= eneTerm2; // eigth power of rm/distance
eneTerm2 *= (4.0 * eneTerm1 - 3.0 * eneTerm2);
eneTerm2 *=
em[d_probetype][d_targettypes[i]]; // multiplication with em
d_eneContrib[i] = eneTerm2;
double t0, ti;
if (d_function[i] == &cos_acc ||
d_function[i] == &cos_2_0) { // only if dependent of ti and t0
double dotProd = d_vectTargetProbe[j] * d_plane[j] +
d_vectTargetProbe[j + 1] * d_plane[j + 1] +
d_vectTargetProbe[j + 2] * d_plane[j + 2];
double vectTargetProbeInPlane[3] = {
d_vectTargetProbe[j] -
d_plane[j] *
dotProd, // projection of targetProbe vector on lp plane
d_vectTargetProbe[j + 1] -
d_plane[j + 1] *
dotProd, // projection of targetProbe vector on lp plane
d_vectTargetProbe[j + 2] -
d_plane[j + 2] *
dotProd}; // projection of targetProbe vector on lp plane
normalize(vectTargetProbeInPlane[0], vectTargetProbeInPlane[1],
vectTargetProbeInPlane[2]);
// angles t_0 (out-of-lonepair-plane angle), t_i (in-lonepair-plane
// angle),
t0 = angle(distN[0], distN[1], distN[2], vectTargetProbeInPlane[0],
vectTargetProbeInPlane[1],
vectTargetProbeInPlane[2]); // out of plane angle
ti = angle(d_direction[j], d_direction[j + 1], d_direction[j + 2],
vectTargetProbeInPlane[0], vectTargetProbeInPlane[1],
vectTargetProbeInPlane[2]); // in plane angle
} else {
t0 = 0.0;
ti = 1.0;
}
d_eneContrib[i] *=
(*(d_function[i]))(t, t0, ti); // scaling of energy contribution
if (d_eneContrib[i] <
minEne) { // check whether most favored interaction
minEne = d_eneContrib[i];
minId = i;
}
} else {
d_eneContrib[i] = 0.0;
}
}
minId *= 3;
probeDirection[0] = -d_vectTargetProbe[minId]; // probe is directed to most
// favored interaction
probeDirection[1] = -d_vectTargetProbe[++minId];
probeDirection[2] = -d_vectTargetProbe[++minId];
normalize(probeDirection[0], probeDirection[1], probeDirection[2]);
for (unsigned int
i = 0,
j = 0;
i < d_nInteract;
i++,
j +=
3) { // scaling to take probe direction into account and adding up
double vectHydrogenTarget[3] = {
probeDirection[0] * bondlength[d_probetype] - d_vectTargetProbe[j],
probeDirection[1] * bondlength[d_probetype] -
d_vectTargetProbe[j + 1],
probeDirection[2] * bondlength[d_probetype] -
d_vectTargetProbe[j + 2]};
normalize(vectHydrogenTarget[0], vectHydrogenTarget[1],
vectHydrogenTarget[2]);
// p (angle between best hydrogen bond direction of probe and
// hydrogen-acceptor vector)
double p = angle(probeDirection[0], probeDirection[1], probeDirection[2],
vectHydrogenTarget[0], vectHydrogenTarget[1],
vectHydrogenTarget[2]);
res +=
d_eneContrib[i] *
cos_2(p, 0.0, 1.0); // scaling of energy contribution and adding up
}
} else { // d_DAprop='D'
unsigned int minId = 0;
double minEne = std::numeric_limits<double>::max(); // minimal energy
double probeDirection[3]; // direction of prope
for (unsigned int i = 0, j = 0; i < d_nInteract;
i++, j += 3) { // calculation of energy contributions and searching
// for the favored probe direction ( direction of
// lowest energy contribution )
d_vectTargetProbe[j] = x - d_pos[j]; // vector of interaction
d_vectTargetProbe[j + 1] = y - d_pos[j + 1];
d_vectTargetProbe[j + 2] = z - d_pos[j + 2];
double dist2 =
d_vectTargetProbe[j] * d_vectTargetProbe[j] +
d_vectTargetProbe[j + 1] * d_vectTargetProbe[j + 1] +
d_vectTargetProbe[j + 2] *
d_vectTargetProbe[j +
2]; // calc of squared length of interaction
// std::cout << dist2 << std::endl;
if (dist2 < thres) {
double dis = sqrt(dist2);
double distN[3] = {d_vectTargetProbe[j] / dis,
d_vectTargetProbe[j + 1] / dis,
d_vectTargetProbe[j + 2] / dis};
dist2 = std::max(dist2, d_cutoff);
double t = angle(
distN[0], distN[1], distN[2], d_direction[j], d_direction[j + 1],
d_direction[j + 2]); // calc of angle between direction of hbond
// and target-probe direction
double eneTerm1 = rm[d_probetype][d_targettypes[i]];
eneTerm1 *= eneTerm1; // squared rm
eneTerm1 /= dist2; // division by squared distance
double eneTerm2 = eneTerm1 * eneTerm1; // fourth power of rm/distance
eneTerm1 = eneTerm2 * eneTerm1; // sixth power of rm/distance
eneTerm2 *= eneTerm2; // eigth power of rm/distance
eneTerm2 *= (4.0 * eneTerm1 - 3.0 * eneTerm2);
eneTerm2 *=
em[d_probetype][d_targettypes[i]]; // multiplication with em
d_eneContrib[i] = eneTerm2;
d_eneContrib[i] *=
(*(d_function[i]))(t, 0.0, 1.0); // scaling of energy contribution
if (d_eneContrib[i] <
minEne) { // check whether most favored interaction
minEne = d_eneContrib[i];
minId = i;
}
} else {
d_eneContrib[i] = 0;
}
}
minId *= 3;
probeDirection[0] = -d_vectTargetProbe[minId]; // probe is directed to most
// favored interaction
probeDirection[1] = -d_vectTargetProbe[++minId];
probeDirection[2] = -d_vectTargetProbe[++minId];
normalize(probeDirection[0], probeDirection[1], probeDirection[2]);
for (unsigned int
i = 0,
j = 0;
i < d_nInteract;
i++,
j +=
3) { // scaling to take probe direction into account and adding up
double vectProbeHydrogen[3] = {
d_direction[j] * d_lengths[i] - d_vectTargetProbe[j],
d_direction[j + 1] * d_lengths[i] - d_vectTargetProbe[j + 1],
d_direction[j + 2] * d_lengths[i] - d_vectTargetProbe[j + 2]};
normalize(vectProbeHydrogen[0], vectProbeHydrogen[1],
vectProbeHydrogen[2]);
// p (angle between best hydrogen bond direction of probe and
// hydrogen-acceptor vector)
double p = angle(vectProbeHydrogen[0], vectProbeHydrogen[1],
vectProbeHydrogen[2], probeDirection[0],
probeDirection[1], probeDirection[2]);
res +=
d_eneContrib[i] *
cos_2(p, 0.0, 1.0); // scaling of energy contribution and adding up
}
}
return res;
}
void HBond::addVectElements(atomtype type,
double (*funct)(double, double, double),
const RDGeom::Point3D &pos,
const RDGeom::Point3D &dir,
const RDGeom::Point3D &plane) {
d_targettypes.push_back(type);
d_function.push_back(funct);
d_pos.push_back(pos.x);
d_pos.push_back(pos.y);
d_pos.push_back(pos.z);
double len = dir.length();
d_lengths.push_back(len);
d_direction.push_back(dir.x / len);
d_direction.push_back(dir.y / len);
d_direction.push_back(dir.z / len);
len = plane.length();
if (len > 1.e-16) {
d_plane.push_back(plane.x / len);
d_plane.push_back(plane.y / len);
d_plane.push_back(plane.z / len);
} else {
d_plane.push_back(0.0);
d_plane.push_back(0.0);
d_plane.push_back(0.0);
}
}
void HBond::normalize(double &x, double &y, double &z) const {
double temp = x * x + y * y + z * z;
temp = sqrt(temp);
x /= temp;
y /= temp;
z /= temp;
}
double HBond::angle(double x1, double y1, double z1, double x2, double y2,
double z2) const {
double dotProd = x1 * x2 + y1 * y2 + z1 * z2;
if (dotProd < -1.0) {
dotProd = -1.0;
} else if (dotProd > 1.0) {
dotProd = 1.0;
}
return acos(dotProd);
}
Hydrophilic::Hydrophilic(const RDKit::ROMol &mol, int confId, bool fixed,
double cutoff) {
d_hbondOH = HBond(mol, confId, "OH", fixed, cutoff);
d_hbondO = HBond(mol, confId, "O", fixed, cutoff);
}
double Hydrophilic::operator()(double x, double y, double z,
double thres) const {
double hbondO, hbondOH;
hbondO = d_hbondO(x, y, z, thres);
hbondOH = d_hbondOH(x, y, z, thres);
return std::min(hbondO, hbondOH);
}
void writeToCubeStream(const RDGeom::UniformRealValueGrid3D &grd,
std::ostream &outStrm, const RDKit::ROMol *mol,
int confid) {
PRECONDITION(outStrm, "bad stream");
const double bohr = 0.529177249;
int dimX = grd.getNumX(); //+2;
int dimY = grd.getNumY(); //+2;
int dimZ = grd.getNumZ(); //+2;
auto spacing = grd.getSpacing() / bohr;
auto offSet = grd.getOffset() / bohr;
auto nAtoms = mol ? mol->getNumAtoms() : 0u;
outStrm.setf(std::ios::left);
outStrm
<< "Gaussian cube format generated by RDKit\n*************************\n";
outStrm << std::setw(20) << std::setprecision(6) << nAtoms << std::setw(20)
<< std::setprecision(6) << offSet.x << std::setw(20)
<< std::setprecision(6) << offSet.y << std::setw(20)
<< std::setprecision(6) << offSet.z << std::endl;
outStrm << std::setw(20) << std::setprecision(6) << dimX << std::setw(20)
<< std::setprecision(6) << spacing << std::setw(20)
<< std::setprecision(6) << 0 << std::setw(20) << std::setprecision(6)
<< 0 << std::endl
<< std::setw(20) << std::setprecision(6) << dimY << std::setw(20)
<< std::setprecision(6) << 0 << std::setw(20) << std::setprecision(6)
<< spacing << std::setw(20) << std::setprecision(6) << 0 << std::endl
<< std::setw(20) << std::setprecision(6) << dimZ << std::setw(20)
<< std::setprecision(6) << std::setw(20) << std::setprecision(6) << 0
<< std::setw(20) << std::setprecision(6) << 0 << std::setw(20)
<< std::setprecision(6) << spacing << std::endl;
if (mol) {
// this will throw a ConformerException if confId does not exist
const auto &conf = mol->getConformer(confid);
for (const auto &atom : mol->atoms()) {
const auto &pt = conf.getAtomPos(atom->getIdx()) / bohr;
outStrm << std::setw(20) << std::setprecision(6) << std::left
<< atom->getAtomicNum() << std::setw(20) << std::setprecision(6)
<< atom->getAtomicNum() << std::setw(20) << std::setprecision(6)
<< pt.x << std::setw(20) << std::setprecision(6) << std::setw(20)
<< std::setprecision(6) << pt.y << std::setw(20)
<< std::setprecision(6) << pt.z << std::endl;
}
}
for (auto xi = 0u; xi < grd.getNumX(); ++xi) {
for (auto yi = 0u; yi < grd.getNumY(); ++yi) {
for (auto zi = 0u; zi < grd.getNumZ(); ++zi) {
outStrm << std::setw(20) << std::setprecision(6) << std::left
<< static_cast<double>(
grd.getVal(grd.getGridIndex(xi, yi, zi)));
// grd->d_numX-xi-1, grd->d_numY-yi-1, grd->d_numZ-zi-1
if ((zi + 1) % 8 == 0) {
outStrm << std::endl;
}
}
outStrm << std::endl;
}
outStrm << std::endl;
}
}
void writeToCubeFile(const RDGeom::UniformRealValueGrid3D &grd,
const std::string &filename, const RDKit::ROMol *mol,
int confid) {
std::ofstream ofStrm(filename.c_str());
writeToCubeStream(grd, ofStrm, mol, confid);
}
std::unique_ptr<RDKit::RWMol> readFromCubeStream(
RDGeom::UniformRealValueGrid3D &grd, std::istream &inStrm) {
PRECONDITION(inStrm, "bad stream");
constexpr double bohr = 0.529177249;
constexpr double spacingThreshold = 0.0001;
if (inStrm.eof()) {
return nullptr;
}
std::string string;
int nAtoms;
std::getline(inStrm, string);
std::getline(inStrm, string);
inStrm >> nAtoms;
double x, y, z;
inStrm >> x >> y >> z;
const RDGeom::Point3D offSet(x * bohr, y * bohr, z * bohr);
int dimX, dimY, dimZ;
double spacingX, spacingY, spacingZ, temp1, temp2;
inStrm >> dimX >> spacingX >> temp1 >> temp2;
inStrm >> dimY >> temp1 >> spacingY >> temp2;
inStrm >> dimZ >> temp1 >> temp2 >> spacingZ;
if ((fabs(spacingX - spacingY) > spacingThreshold) ||
(fabs(spacingX - spacingZ) > spacingThreshold)) {
std::ostringstream errout;
errout << "Same spacing in all directions needed";
throw RDKit::FileParseException(errout.str());
} else {
spacingX *= bohr;
grd = RDGeom::UniformRealValueGrid3D(spacingX * static_cast<double>(dimX),
spacingX * static_cast<double>(dimY),
spacingX * static_cast<double>(dimZ),
spacingX, &offSet);
}
std::unique_ptr<RDKit::RWMol> molecule;
if (nAtoms) {
molecule.reset(new RDKit::RWMol());
std::unique_ptr<RDKit::Conformer> conf(new RDKit::Conformer(nAtoms));
int atomNum;
for (auto i = 0; i < nAtoms; ++i) {
inStrm >> atomNum >> temp1 >> x >> y >> z;
RDKit::Atom atom(atomNum);
molecule->addAtom(&atom, true, false);
RDGeom::Point3D pos(x * bohr, y * bohr, z * bohr);
conf->setAtomPos(i, pos);
}
molecule->addConformer(conf.release(), false);
}
for (auto xi = 0; xi < dimX; ++xi) {
for (auto yi = 0; yi < dimY; ++yi) {
for (auto zi = 0; zi < dimZ; ++zi) {
double tempVal;
inStrm >> tempVal;
grd.setVal(grd.getGridIndex(xi, yi, zi), tempVal);
}
}
}
return molecule;
}
std::unique_ptr<RDKit::RWMol> readFromCubeFile(
RDGeom::UniformRealValueGrid3D &grd, const std::string &filename) {
std::ifstream ifStrm(filename.c_str());
if (ifStrm.bad()) {
std::ostringstream errout;
errout << "Bad input file " << filename;
throw RDKit::BadFileException(errout.str());
};
return readFromCubeStream(grd, ifStrm);
}
} // namespace RDMIF
Reference in New Issue View Git Blame Copy Permalink