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
master
rdkit/Code/ForceField/MMFF/testMMFFForceField.cpp
Ricardo Rodriguez 7b7a8a4e17 Refactor iostreams includes (#8846) 
* refactor iostreams includes

* restore ostream to MonomerInfo.cpp
2025-10-08 16:08:01 +02:00

1885 lines
74 KiB
C++
Raw Permalink Blame History

//
// Copyright (C) 2013-2024 Paolo Tosco 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 <RDGeneral/test.h>
#include <GraphMol/RDKitBase.h>
#include <GraphMol/FileParsers/MolSupplier.h>
#include <GraphMol/FileParsers/MolWriters.h>
#include <GraphMol/FileParsers/FileParsers.h>
#include <ForceField/MMFF/Params.h>
#include <ForceField/DistanceConstraints.h>
#include <ForceField/AngleConstraints.h>
#include <ForceField/MMFF/TorsionConstraint.h>
#include <ForceField/MMFF/PositionConstraint.h>
#include <GraphMol/ForceFieldHelpers/MMFF/AtomTyper.h>
#include <GraphMol/ForceFieldHelpers/MMFF/Builder.h>
#include <GraphMol/MolOps.h>
#include <GraphMol/DistGeomHelpers/Embedder.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include <ForceField/ForceField.h>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <string>
#include "testMMFFForceField.h"
#include <utility>
#define FCON_TOLERANCE 0.05
#define ENERGY_TOLERANCE 0.025
using namespace RDKit;
bool fexist(std::string filename) {
std::ifstream ifStream(filename.c_str());
return !ifStream.fail();
}
bool fgrep(std::fstream &fStream, std::string key, std::string &line) {
bool found = false;
std::streampos current = fStream.tellg();
while ((!found) &&
(!(std::getline(fStream, line).rdstate() & std::ifstream::failbit))) {
found = (line.find(key) != std::string::npos);
}
if (!found) {
fStream.seekg(current);
}
return found;
}
bool fgrep(std::fstream &rdkFStream, std::string key) {
std::string line;
return fgrep(rdkFStream, key, line);
}
bool getLineByNum(std::fstream &stream, std::streampos startPos, unsigned int n,
std::string &line) {
std::streampos current = stream.tellg();
stream.seekg(startPos);
bool fail;
unsigned int i = 0;
while ((i <= n) && (!(fail = std::getline(stream, line).rdstate() &
std::ifstream::failbit))) {
++i;
}
stream.seekg(current);
return (!fail);
}
bool sortBondStretchInstanceVec(BondStretchInstance *a,
BondStretchInstance *b) {
bool crit = false;
if (a->iAtomType == b->iAtomType) {
if (a->jAtomType == b->jAtomType) {
crit = (a->ffType < b->ffType);
} else {
crit = (a->jAtomType < b->jAtomType);
}
} else {
crit = (a->iAtomType < b->iAtomType);
}
return crit;
}
bool sortAngleBendInstanceVec(AngleBendInstance *a, AngleBendInstance *b) {
bool crit = false;
if (a->jAtomType == b->jAtomType) {
if (a->iAtomType == b->iAtomType) {
if (a->kAtomType == b->kAtomType) {
crit = (a->ffType < b->ffType);
} else {
crit = (a->kAtomType < b->kAtomType);
}
} else {
crit = (a->iAtomType < b->iAtomType);
}
} else {
crit = (a->jAtomType < b->jAtomType);
}
return crit;
}
bool sortStretchBendInstanceVec(StretchBendInstance *a,
StretchBendInstance *b) {
bool crit = false;
if (a->jAtomType == b->jAtomType) {
if (a->iAtomType == b->iAtomType) {
if (a->kAtomType == b->kAtomType) {
crit = (a->ffType < b->ffType);
} else {
crit = (a->kAtomType < b->kAtomType);
}
} else {
crit = (a->iAtomType < b->iAtomType);
}
} else {
crit = (a->jAtomType < b->jAtomType);
}
return crit;
}
bool sortOopBendInstanceVec(OopBendInstance *a, OopBendInstance *b) {
bool crit = false;
if (a->jAtomType == b->jAtomType) {
if (a->iAtomType == b->iAtomType) {
if (a->kAtomType == b->kAtomType) {
crit = (a->lAtomType < b->lAtomType);
} else {
crit = (a->kAtomType < b->kAtomType);
}
} else {
crit = (a->iAtomType < b->iAtomType);
}
} else {
crit = (a->jAtomType < b->jAtomType);
}
return crit;
}
bool sortTorsionInstanceVec(TorsionInstance *a, TorsionInstance *b) {
bool crit = false;
if (a->jAtomType == b->jAtomType) {
if (a->kAtomType == b->kAtomType) {
if (a->iAtomType == b->iAtomType) {
if (a->lAtomType == b->lAtomType) {
crit = (a->ffType < b->ffType);
} else {
crit = (a->lAtomType < b->lAtomType);
}
} else {
crit = (a->iAtomType < b->iAtomType);
}
} else {
crit = (a->kAtomType < b->kAtomType);
}
} else {
crit = (a->jAtomType < b->jAtomType);
}
return crit;
}
void fixBondStretchInstance(BondStretchInstance *bondStretchInstance) {
if (bondStretchInstance->iAtomType > bondStretchInstance->jAtomType) {
std::swap(bondStretchInstance->iAtomType, bondStretchInstance->jAtomType);
}
}
void fixAngleBendInstance(AngleBendInstance *angleBendInstance) {
if (angleBendInstance->iAtomType > angleBendInstance->kAtomType) {
std::swap(angleBendInstance->iAtomType, angleBendInstance->kAtomType);
}
}
void fixOopBendInstance(OopBendInstance *oopBendInstance) {
if (oopBendInstance->iAtomType > oopBendInstance->kAtomType) {
std::swap(oopBendInstance->iAtomType, oopBendInstance->kAtomType);
}
}
void fixTorsionInstance(TorsionInstance *torsionInstance) {
if (torsionInstance->jAtomType > torsionInstance->kAtomType) {
std::swap(torsionInstance->jAtomType, torsionInstance->kAtomType);
std::swap(torsionInstance->iAtomType, torsionInstance->lAtomType);
} else if (torsionInstance->jAtomType == torsionInstance->kAtomType) {
if (torsionInstance->iAtomType > torsionInstance->lAtomType) {
std::swap(torsionInstance->iAtomType, torsionInstance->lAtomType);
}
}
}
int mmffValidationSuite(int argc, char *argv[]) {
std::string arg;
std::string ffVariant = "";
std::vector<std::string> ffVec;
std::string verbosity;
std::string molFile = "";
std::string molType = "";
std::vector<std::string> molFileVec;
std::vector<std::string> molTypeVec;
std::string rdk = "";
std::string ref = "";
std::streampos rdkPos;
std::streampos rdkCurrent;
std::streampos refPos;
std::streampos refCurrent;
int iarg = 1;
unsigned int n = 0;
bool error = false;
bool fullTest = false;
int inc = 4;
bool testFailure = false;
while (iarg < argc) {
arg = argv[iarg];
error = (arg.at(0) != '-');
if (error) {
break;
}
error = ((arg == "-h") || (arg == "--help"));
if (error) {
break;
} else if (arg == "-f") {
error = ((iarg + 1) >= argc);
if (error) {
break;
}
ffVariant = argv[iarg + 1];
error = ((ffVariant != "MMFF94") && (ffVariant != "MMFF94s"));
if (error) {
break;
}
++iarg;
} else if ((arg == "-sdf") || (arg == "-smi")) {
molType = arg.substr(1);
if ((iarg + 1) < argc) {
molFile = argv[iarg + 1];
++iarg;
}
} else if (arg == "-L") {
fullTest = true;
inc = 1;
} else if (arg == "-l") {
error = ((iarg + 1) >= argc);
if (error) {
break;
}
rdk = argv[iarg + 1];
++iarg;
}
++iarg;
}
if (error) {
std::cerr << "Usage: testMMFFForceField [-L] [-f {MMFF94|MMFF94s}] "
"[{-sdf [<sdf_file>] | -smi [<smiles_file>]}] [-l <log_file>]"
<< std::endl;
return -1;
}
std::cerr << "-------------------------------------" << std::endl;
std::cerr << "Official MMFF validation suite." << std::endl;
std::string pathName = getenv("RDBASE");
pathName += "/Code/ForceField/MMFF/test_data/";
if (molFile != "") {
ffVariant = ((molFile.substr(0, 7) == "MMFF94s") ? "MMFF94s" : "MMFF94");
}
if (ffVariant == "") {
ffVec.push_back("MMFF94");
if (fullTest) {
ffVec.push_back("MMFF94s");
}
} else {
ffVec.push_back(ffVariant);
}
std::fstream rdkFStream;
std::fstream refFStream;
if (rdk == "") {
rdk = pathName + "testMMFFForceField.log";
}
bool firstTest = true;
if (molType == "") {
molTypeVec.push_back("sdf");
molTypeVec.push_back("smi");
} else {
molTypeVec.push_back(molType);
}
boost::logging::disable_logs("rdApp.error");
for (auto &molTypeIt : molTypeVec) {
bool checkEnergy = (molTypeIt == "sdf");
for (auto &ffIt : ffVec) {
ref = pathName + ffIt + "_reference.log";
molFileVec.clear();
if (molFile == "") {
molFileVec.push_back(pathName + ffIt + "_dative." + molTypeIt);
molFileVec.push_back(pathName + ffIt + "_hypervalent." + molTypeIt);
} else {
molFileVec.push_back(molFile);
}
for (auto &molFileIt : molFileVec) {
std::vector<ROMol *> molVec;
molVec.clear();
if (molTypeIt == "sdf") {
SDMolSupplier sdfSupplier(molFileIt, false, false);
for (size_t ii = 0; ii < sdfSupplier.length(); ++ii) {
molVec.push_back(sdfSupplier[ii]);
}
sdfSupplier.reset();
} else {
SmilesMolSupplier smiSupplier(molFileIt);
for (size_t ii = 0; ii < smiSupplier.length(); ++ii) {
molVec.push_back(smiSupplier[ii]);
}
smiSupplier.reset();
}
SDWriter *sdfWriter =
new SDWriter(molFileIt.substr(0, molFileIt.length() - 4) + "_min" +
((molTypeIt == "smi") ? "_from_SMILES" : "") + ".sdf");
ROMol *mol;
std::string molName;
std::vector<std::string> nameArray;
rdkFStream.open(rdk.c_str(),
(firstTest ? std::fstream::out | std::fstream::binary
: (std::fstream::out | std::fstream::binary |
std::fstream::app)));
std::string computingKey = ffIt + " energies for " + molFileIt;
std::cerr << std::endl
<< "Computing " << computingKey << "..." << std::endl;
for (size_t ii = 0; ii < computingKey.length(); ++ii) {
rdkFStream << "*";
}
rdkFStream << std::endl << computingKey << std::endl;
for (size_t ii = 0; ii < computingKey.length(); ++ii) {
rdkFStream << "*";
}
rdkFStream << std::endl << std::endl;
for (size_t ii = 0; ii < molVec.size(); ii += inc) {
mol = molVec[ii];
if (molTypeIt == "smi") {
DGeomHelpers::EmbedMolecule(*mol);
MolOps::addHs((RWMol &)*mol, false, true);
}
MMFF::sanitizeMMFFMol((RWMol &)(*mol));
if (mol->hasProp(common_properties::_Name)) {
mol->getProp(common_properties::_Name, molName);
rdkFStream << molName << std::endl;
nameArray.push_back(molName);
}
auto *mmffMolProperties = new MMFF::MMFFMolProperties(
*mol, ffIt, MMFF::MMFF_VERBOSITY_HIGH, rdkFStream);
if ((!mmffMolProperties) || (!(mmffMolProperties->isValid()))) {
std::cerr << molName + ": error setting up force-field"
<< std::endl;
continue;
} else {
ForceFields::ForceField *field = MMFF::constructForceField(
*mol, mmffMolProperties, 100.0, -1, false);
field->initialize();
field->minimize();
sdfWriter->write(*mol);
rdkFStream << "\nTOTAL MMFF ENERGY =" << std::right
<< std::setw(16) << std::fixed << std::setprecision(4)
<< field->calcEnergy() << "\n\n"
<< std::endl;
delete field;
}
delete mmffMolProperties;
}
for (auto m : molVec) {
delete m;
}
sdfWriter->close();
delete sdfWriter;
rdkFStream.close();
std::cerr << "Checking against " << ref << "..." << std::endl;
rdkFStream.open(rdk.c_str(), std::fstream::in | std::fstream::binary);
refFStream.open(ref.c_str(), std::fstream::in | std::fstream::binary);
bool found = false;
std::string skip;
std::string errorMsg;
std::string corruptedMsg;
std::string energyString;
std::vector<BondStretchInstance *> rdkBondStretchInstanceVec;
std::vector<BondStretchInstance *> refBondStretchInstanceVec;
std::vector<AngleBendInstance *> rdkAngleBendInstanceVec;
std::vector<AngleBendInstance *> refAngleBendInstanceVec;
std::vector<StretchBendInstance *> rdkStretchBendInstanceVec;
std::vector<StretchBendInstance *> refStretchBendInstanceVec;
std::vector<OopBendInstance *> rdkOopBendInstanceVec;
std::vector<OopBendInstance *> refOopBendInstanceVec;
std::vector<TorsionInstance *> rdkTorsionInstanceVec;
std::vector<TorsionInstance *> refTorsionInstanceVec;
double rdkEnergy;
double rdkEleEnergy;
double refEnergy;
double refVdWEnergy;
double refEleEnergy;
unsigned int n_failed = 0;
bool failed = false;
rdkFStream.seekg(0);
fgrep(rdkFStream, computingKey);
size_t ii;
for (ii = 0; ii < nameArray.size(); ii += inc) {
error = false;
failed = false;
errorMsg = rdk + ": Molecule not found";
found = fgrep(rdkFStream, nameArray[ii]);
if (!found) {
break;
}
errorMsg = ref + ": Molecule not found";
found = fgrep(refFStream, nameArray[ii]);
if (!found) {
break;
}
refPos = refFStream.tellg();
found = false;
bool refHaveBondStretch = false;
bool refHaveAngleBend = false;
bool refHaveStretchBend = false;
bool refHaveOopBend = false;
bool refHaveTorsion = false;
bool refHaveVdW = false;
bool refHaveEle = false;
std::string rdkLine;
std::string refLine;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
std::size_t occ = refLine.find("****");
found = (occ != std::string::npos);
if (!found) {
if (!refHaveVdW) {
occ = refLine.find("Net vdW");
refHaveVdW = (occ != std::string::npos);
if (refHaveVdW) {
std::stringstream refVdWStringStream(refLine);
refVdWStringStream >> skip >> skip >> refVdWEnergy;
}
}
if (!refHaveEle) {
occ = refLine.find("Electrostatic");
refHaveEle = (occ != std::string::npos);
if (refHaveEle) {
std::stringstream refEleStringStream(refLine);
refEleStringStream >> skip >> refEleEnergy;
}
}
if (!refHaveBondStretch) {
occ = refLine.find("B O N D S T R E T C H I N G");
refHaveBondStretch = (occ != std::string::npos);
}
if (!refHaveAngleBend) {
occ = refLine.find("A N G L E B E N D I N G");
refHaveAngleBend = (occ != std::string::npos);
}
if (!refHaveStretchBend) {
occ = refLine.find("S T R E T C H B E N D I N G");
refHaveStretchBend = (occ != std::string::npos);
}
if (!refHaveOopBend) {
occ = refLine.find("O U T - O F - P L A N E B E N D I N G");
refHaveOopBend = (occ != std::string::npos);
}
if (!refHaveTorsion) {
occ = refLine.find("T O R S I O N A L");
refHaveTorsion = (occ != std::string::npos);
}
}
}
refFStream.seekg(refPos);
errorMsg = "";
// --------------------------------------------------
// B O N D S T R E T C H I N G
// --------------------------------------------------
if (refHaveBondStretch) {
corruptedMsg = ": Bond stretching: corrupted input data\n";
found = fgrep(rdkFStream, "B O N D S T R E T C H I N G");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = false;
rdkBondStretchInstanceVec.clear();
rdkPos = rdkFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
std::ifstream::failbit))) {
found = (!(rdkLine.length()));
if (found) {
break;
}
std::stringstream rdkLineStream(rdkLine);
auto *bondStretchInstance = new BondStretchInstance();
bondStretchInstance->idx = n;
rdkLineStream >> skip >> skip >> skip >> skip >>
bondStretchInstance->iAtomType >>
bondStretchInstance->jAtomType >>
bondStretchInstance->ffType >> skip >> skip >> skip >> skip >>
bondStretchInstance->kb;
fixBondStretchInstance(bondStretchInstance);
rdkBondStretchInstanceVec.push_back(bondStretchInstance);
++n;
}
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::sort(rdkBondStretchInstanceVec.begin(),
rdkBondStretchInstanceVec.end(),
sortBondStretchInstanceVec);
found =
fgrep(rdkFStream, "TOTAL BOND STRETCH ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::stringstream rdkBondStretchStringStream(energyString);
rdkBondStretchStringStream >> skip >> skip >> skip >> skip >>
skip >> rdkEnergy;
found = fgrep(refFStream, "B O N D S T R E T C H I N G");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = fgrep(refFStream, "----------------");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = false;
refBondStretchInstanceVec.clear();
refPos = refFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
found = (!(refLine.length()));
if (found) {
break;
}
std::stringstream refLineStream(refLine);
auto *bondStretchInstance = new BondStretchInstance();
bondStretchInstance->idx = n;
refLineStream >> skip >> skip >> skip >> skip >>
bondStretchInstance->iAtomType >>
bondStretchInstance->jAtomType >>
bondStretchInstance->ffType >> skip >> skip >> skip >> skip >>
bondStretchInstance->kb;
fixBondStretchInstance(bondStretchInstance);
refBondStretchInstanceVec.push_back(bondStretchInstance);
++n;
}
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::sort(refBondStretchInstanceVec.begin(),
refBondStretchInstanceVec.end(),
sortBondStretchInstanceVec);
found = fgrep(refFStream, "TOTAL BOND STRAIN ENERGY", energyString);
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream refBondStretchStringStream(energyString);
refBondStretchStringStream >> skip >> skip >> skip >> skip >>
skip >> refEnergy;
error = (rdkBondStretchInstanceVec.size() !=
refBondStretchInstanceVec.size());
if (error) {
failed = true;
std::stringstream diff;
diff << "Bond stretching: expected "
<< refBondStretchInstanceVec.size()
<< " interactions, found "
<< rdkBondStretchInstanceVec.size() << "\n";
errorMsg += diff.str();
}
for (unsigned j = 0;
(!error) & (j < rdkBondStretchInstanceVec.size()); ++j) {
error =
((rdkBondStretchInstanceVec[j]->iAtomType !=
refBondStretchInstanceVec[j]->iAtomType) ||
(rdkBondStretchInstanceVec[j]->jAtomType !=
refBondStretchInstanceVec[j]->jAtomType) ||
(rdkBondStretchInstanceVec[j]->ffType !=
refBondStretchInstanceVec[j]->ffType) ||
(fabs(rdkBondStretchInstanceVec[j]->kb -
refBondStretchInstanceVec[j]->kb) > FCON_TOLERANCE));
if (error) {
failed = true;
bool haveLine;
haveLine =
getLineByNum(rdkFStream, rdkPos,
rdkBondStretchInstanceVec[j]->idx, rdkLine);
if (haveLine) {
haveLine =
getLineByNum(refFStream, refPos,
refBondStretchInstanceVec[j]->idx, refLine);
}
std::stringstream diff;
diff << "Bond stretching: found a difference\n";
if (haveLine) {
diff << "Expected:\n"
<< rdkLine << "\nFound:\n"
<< refLine << "\n";
}
errorMsg += diff.str();
}
delete rdkBondStretchInstanceVec[j];
delete refBondStretchInstanceVec[j];
}
error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Bond stretching: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
}
// --------------------------------------------------
// A N G L E B E N D I N G
// --------------------------------------------------
if (refHaveAngleBend) {
corruptedMsg = ": Angle bending: corrupted input data\n";
found = fgrep(rdkFStream, "A N G L E B E N D I N G");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = false;
rdkAngleBendInstanceVec.clear();
rdkPos = rdkFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
std::ifstream::failbit))) {
found = (!(rdkLine.length()));
if (found) {
break;
}
std::stringstream rdkLineStream(rdkLine);
auto *angleBendInstance = new AngleBendInstance();
angleBendInstance->idx = n;
rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
angleBendInstance->iAtomType >>
angleBendInstance->jAtomType >>
angleBendInstance->kAtomType >> angleBendInstance->ffType >>
skip >> skip >> skip >> skip >> angleBendInstance->ka;
fixAngleBendInstance(angleBendInstance);
rdkAngleBendInstanceVec.push_back(angleBendInstance);
++n;
}
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::sort(rdkAngleBendInstanceVec.begin(),
rdkAngleBendInstanceVec.end(), sortAngleBendInstanceVec);
found = fgrep(rdkFStream, "TOTAL ANGLE BEND ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::stringstream rdkAngleBendStringStream(energyString);
rdkAngleBendStringStream >> skip >> skip >> skip >> skip >> skip >>
rdkEnergy;
found = fgrep(refFStream, "A N G L E B E N D I N G");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = fgrep(refFStream, "----------------");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = false;
refAngleBendInstanceVec.clear();
refPos = refFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
found = (!(refLine.length()));
if (found) {
break;
}
std::stringstream refLineStream(refLine);
auto *angleBendInstance = new AngleBendInstance();
angleBendInstance->idx = n;
refLineStream >> skip >> skip >> skip >> skip >>
angleBendInstance->iAtomType >>
angleBendInstance->jAtomType >>
angleBendInstance->kAtomType >> angleBendInstance->ffType >>
skip >> skip >> skip >> skip >> angleBendInstance->ka;
fixAngleBendInstance(angleBendInstance);
refAngleBendInstanceVec.push_back(angleBendInstance);
++n;
}
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::sort(refAngleBendInstanceVec.begin(),
refAngleBendInstanceVec.end(), sortAngleBendInstanceVec);
found =
fgrep(refFStream, "TOTAL ANGLE STRAIN ENERGY", energyString);
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream refAngleBendStringStream(energyString);
refAngleBendStringStream >> skip >> skip >> skip >> skip >> skip >>
refEnergy;
error = (rdkAngleBendInstanceVec.size() !=
refAngleBendInstanceVec.size());
if (error) {
failed = true;
std::stringstream diff;
diff << "Angle bending: expected "
<< refAngleBendInstanceVec.size() << " interactions, found "
<< rdkAngleBendInstanceVec.size() << "\n";
errorMsg += diff.str();
}
for (unsigned j = 0;
(!error) & (j < rdkAngleBendInstanceVec.size()); ++j) {
error = ((rdkAngleBendInstanceVec[j]->iAtomType !=
refAngleBendInstanceVec[j]->iAtomType) ||
(rdkAngleBendInstanceVec[j]->jAtomType !=
refAngleBendInstanceVec[j]->jAtomType) ||
(rdkAngleBendInstanceVec[j]->kAtomType !=
refAngleBendInstanceVec[j]->kAtomType) ||
(rdkAngleBendInstanceVec[j]->ffType !=
refAngleBendInstanceVec[j]->ffType) ||
(fabs(rdkAngleBendInstanceVec[j]->ka -
refAngleBendInstanceVec[j]->ka) > FCON_TOLERANCE));
if (error) {
failed = true;
bool haveLine;
haveLine =
getLineByNum(rdkFStream, rdkPos,
rdkAngleBendInstanceVec[j]->idx, rdkLine);
if (haveLine) {
haveLine =
getLineByNum(refFStream, refPos,
refAngleBendInstanceVec[j]->idx, refLine);
}
std::stringstream diff;
diff << "Angle bending: found a difference\n";
if (haveLine) {
diff << "Expected:\n"
<< rdkLine << "\nFound:\n"
<< refLine << "\n";
}
errorMsg += diff.str();
}
delete rdkAngleBendInstanceVec[j];
delete refAngleBendInstanceVec[j];
}
error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Angle bending: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
}
// --------------------------------------------------
// S T R E T C H B E N D I N G
// --------------------------------------------------
if (refHaveStretchBend) {
corruptedMsg = ": Stretch-bending: corrupted input data\n";
found = fgrep(rdkFStream, "S T R E T C H B E N D I N G");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = false;
rdkStretchBendInstanceVec.clear();
rdkPos = rdkFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
std::ifstream::failbit))) {
found = (!(rdkLine.length()));
if (found) {
break;
}
std::stringstream rdkLineStream(rdkLine);
auto *stretchBendInstance = new StretchBendInstance();
stretchBendInstance->idx = n;
rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
stretchBendInstance->iAtomType >>
stretchBendInstance->jAtomType >>
stretchBendInstance->kAtomType >>
stretchBendInstance->ffType >> skip >> skip >> skip >> skip >>
skip >> stretchBendInstance->kba;
rdkStretchBendInstanceVec.push_back(stretchBendInstance);
++n;
}
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::sort(rdkStretchBendInstanceVec.begin(),
rdkStretchBendInstanceVec.end(),
sortStretchBendInstanceVec);
found =
fgrep(rdkFStream, "TOTAL STRETCH-BEND ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::stringstream rdkStretchBendStringStream(energyString);
rdkStretchBendStringStream >> skip >> skip >> skip >> skip >>
rdkEnergy;
found = fgrep(refFStream, "S T R E T C H B E N D I N G");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = fgrep(refFStream, "----------------");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = false;
refStretchBendInstanceVec.clear();
refPos = refFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
found = (!(refLine.length()));
if (found) {
break;
}
std::stringstream refLineStream(refLine);
auto *stretchBendInstance = new StretchBendInstance();
stretchBendInstance->idx = n;
refLineStream >> skip >> skip >> skip >> skip >>
stretchBendInstance->iAtomType >>
stretchBendInstance->jAtomType >>
stretchBendInstance->kAtomType >>
stretchBendInstance->ffType >> skip >> skip >> skip >> skip >>
stretchBendInstance->kba;
refStretchBendInstanceVec.push_back(stretchBendInstance);
++n;
}
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::sort(refStretchBendInstanceVec.begin(),
refStretchBendInstanceVec.end(),
sortStretchBendInstanceVec);
found = fgrep(refFStream, "TOTAL STRETCH-BEND STRAIN ENERGY",
energyString);
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream refStretchBendStringStream(energyString);
refStretchBendStringStream >> skip >> skip >> skip >> skip >>
skip >> refEnergy;
error = (rdkStretchBendInstanceVec.size() !=
refStretchBendInstanceVec.size());
if (error) {
failed = true;
std::stringstream diff;
diff << "Stretch-bending: expected "
<< refStretchBendInstanceVec.size()
<< " interactions, found "
<< rdkStretchBendInstanceVec.size() << "\n";
errorMsg += diff.str();
}
for (unsigned j = 0;
(!error) & (j < rdkStretchBendInstanceVec.size()); ++j) {
error =
((rdkStretchBendInstanceVec[j]->iAtomType !=
refStretchBendInstanceVec[j]->iAtomType) ||
(rdkStretchBendInstanceVec[j]->jAtomType !=
refStretchBendInstanceVec[j]->jAtomType) ||
(rdkStretchBendInstanceVec[j]->kAtomType !=
refStretchBendInstanceVec[j]->kAtomType) ||
(rdkStretchBendInstanceVec[j]->ffType !=
refStretchBendInstanceVec[j]->ffType) ||
(fabs(rdkStretchBendInstanceVec[j]->kba -
refStretchBendInstanceVec[j]->kba) > FCON_TOLERANCE));
if (error) {
failed = true;
bool haveLine;
haveLine =
getLineByNum(rdkFStream, rdkPos,
rdkStretchBendInstanceVec[j]->idx, rdkLine);
if (haveLine) {
haveLine =
getLineByNum(refFStream, refPos,
refStretchBendInstanceVec[j]->idx, refLine);
}
std::stringstream diff;
diff << "Stretch-bending: found a difference\n";
if (haveLine) {
diff << "Expected:\n"
<< rdkLine << "\nFound:\n"
<< refLine << "\n";
}
errorMsg += diff.str();
}
delete rdkStretchBendInstanceVec[j];
delete refStretchBendInstanceVec[j];
}
error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Stretch-bending: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
}
// --------------------------------------------------
// O U T - O F - P L A N E B E N D I N G
// --------------------------------------------------
if (refHaveOopBend) {
corruptedMsg = ": Out-of-plane bending: corrupted input data\n";
found =
fgrep(rdkFStream, "O U T - O F - P L A N E B E N D I N G");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = false;
rdkOopBendInstanceVec.clear();
rdkPos = rdkFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
std::ifstream::failbit))) {
found = (!(rdkLine.length()));
if (found) {
break;
}
std::stringstream rdkLineStream(rdkLine);
auto *oopBendInstance = new OopBendInstance();
oopBendInstance->idx = n;
rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
skip >> skip >> oopBendInstance->iAtomType >>
oopBendInstance->jAtomType >> oopBendInstance->kAtomType >>
oopBendInstance->lAtomType >> skip >> skip >>
oopBendInstance->koop;
fixOopBendInstance(oopBendInstance);
rdkOopBendInstanceVec.push_back(oopBendInstance);
++n;
}
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::sort(rdkOopBendInstanceVec.begin(),
rdkOopBendInstanceVec.end(), sortOopBendInstanceVec);
found = fgrep(rdkFStream, "TOTAL OUT-OF-PLANE BEND ENERGY",
energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::stringstream rdkOopBendStringStream(energyString);
rdkOopBendStringStream >> skip >> skip >> skip >> skip >> skip >>
rdkEnergy;
found =
fgrep(refFStream, "O U T - O F - P L A N E B E N D I N G");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = fgrep(refFStream, "----------------");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = false;
refOopBendInstanceVec.clear();
refPos = refFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
found = (!(refLine.length()));
if (found) {
break;
}
std::stringstream refLineStream(refLine);
auto *oopBendInstance = new OopBendInstance();
oopBendInstance->idx = n;
refLineStream >> skip >> skip >> skip >> skip >> skip >>
oopBendInstance->iAtomType >> oopBendInstance->jAtomType >>
oopBendInstance->kAtomType >> oopBendInstance->lAtomType >>
skip >> skip >> oopBendInstance->koop;
fixOopBendInstance(oopBendInstance);
refOopBendInstanceVec.push_back(oopBendInstance);
++n;
}
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::sort(refOopBendInstanceVec.begin(),
refOopBendInstanceVec.end(), sortOopBendInstanceVec);
found = fgrep(refFStream, "TOTAL OUT-OF-PLANE STRAIN ENERGY",
energyString);
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream refOopBendStringStream(energyString);
refOopBendStringStream >> skip >> skip >> skip >> skip >> skip >>
refEnergy;
error =
(rdkOopBendInstanceVec.size() != refOopBendInstanceVec.size());
if (error) {
failed = true;
std::stringstream diff;
diff << "Out-of-plane bending: expected "
<< refOopBendInstanceVec.size() << " interactions, found "
<< rdkOopBendInstanceVec.size() << "\n";
errorMsg += diff.str();
}
for (unsigned j = 0; (!error) & (j < rdkOopBendInstanceVec.size());
++j) {
error = ((rdkOopBendInstanceVec[j]->iAtomType !=
refOopBendInstanceVec[j]->iAtomType) ||
(rdkOopBendInstanceVec[j]->jAtomType !=
refOopBendInstanceVec[j]->jAtomType) ||
(rdkOopBendInstanceVec[j]->kAtomType !=
refOopBendInstanceVec[j]->kAtomType) ||
(rdkOopBendInstanceVec[j]->lAtomType !=
refOopBendInstanceVec[j]->lAtomType) ||
(fabs(rdkOopBendInstanceVec[j]->koop -
refOopBendInstanceVec[j]->koop) > FCON_TOLERANCE));
if (error) {
failed = true;
bool haveLine;
haveLine = getLineByNum(rdkFStream, rdkPos,
rdkOopBendInstanceVec[j]->idx, rdkLine);
if (haveLine) {
haveLine =
getLineByNum(refFStream, refPos,
refOopBendInstanceVec[j]->idx, refLine);
}
std::stringstream diff;
diff << "Out-of-plane bending: found a difference\n";
if (haveLine) {
diff << "Expected:\n"
<< rdkLine << "\nFound:\n"
<< refLine << "\n";
}
errorMsg += diff.str();
}
delete rdkOopBendInstanceVec[j];
delete refOopBendInstanceVec[j];
}
error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Out-of-plane bending: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
}
// --------------------------------------------------
// T O R S I O N A L
// --------------------------------------------------
if (refHaveTorsion) {
corruptedMsg = ": Torsional: corrupted input data\n";
found = fgrep(rdkFStream, "T O R S I O N A L");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = false;
rdkTorsionInstanceVec.clear();
rdkPos = rdkFStream.tellg();
n = 0;
while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
std::ifstream::failbit))) {
found = (!(rdkLine.length()));
if (found) {
break;
}
std::stringstream rdkLineStream(rdkLine);
auto *torsionInstance = new TorsionInstance();
torsionInstance->idx = n;
rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
skip >> skip >> torsionInstance->iAtomType >>
torsionInstance->jAtomType >> torsionInstance->kAtomType >>
torsionInstance->lAtomType >> torsionInstance->ffType >>
skip >> skip >> torsionInstance->V1 >> torsionInstance->V2 >>
torsionInstance->V3;
fixTorsionInstance(torsionInstance);
rdkTorsionInstanceVec.push_back(torsionInstance);
++n;
}
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::sort(rdkTorsionInstanceVec.begin(),
rdkTorsionInstanceVec.end(), sortTorsionInstanceVec);
found = fgrep(rdkFStream, "TOTAL TORSIONAL ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
std::stringstream rdkTorsionStringStream(energyString);
rdkTorsionStringStream >> skip >> skip >> skip >> skip >> rdkEnergy;
found = fgrep(refFStream, "T O R S I O N A L");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = fgrep(refFStream, "----------------");
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
found = false;
refTorsionInstanceVec.clear();
refPos = refFStream.tellg();
n = 0;
double refTorsionEnergy;
while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
std::ifstream::failbit))) {
found = (!(refLine.length()));
if (found) {
break;
}
std::stringstream refLineStream(refLine);
auto *torsionInstance = new TorsionInstance();
torsionInstance->idx = n;
refLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
torsionInstance->iAtomType >> torsionInstance->jAtomType >>
torsionInstance->kAtomType >> torsionInstance->lAtomType >>
torsionInstance->ffType >> skip >> refTorsionEnergy >>
torsionInstance->V1 >> torsionInstance->V2 >>
torsionInstance->V3;
if (!(ForceFields::MMFF::isDoubleZero(torsionInstance->V1) &&
ForceFields::MMFF::isDoubleZero(torsionInstance->V2) &&
ForceFields::MMFF::isDoubleZero(torsionInstance->V3) &&
ForceFields::MMFF::isDoubleZero(refTorsionEnergy))) {
fixTorsionInstance(torsionInstance);
refTorsionInstanceVec.push_back(torsionInstance);
} else {
delete torsionInstance;
}
++n;
}
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::sort(refTorsionInstanceVec.begin(),
refTorsionInstanceVec.end(), sortTorsionInstanceVec);
found =
fgrep(refFStream, "TOTAL TORSION STRAIN ENERGY", energyString);
if (!found) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream refTorsionStringStream(energyString);
refTorsionStringStream >> skip >> skip >> skip >> skip >> skip >>
refEnergy;
error =
(rdkTorsionInstanceVec.size() != refTorsionInstanceVec.size());
if (error) {
failed = true;
std::stringstream diff;
diff << "Torsional: expected " << refTorsionInstanceVec.size()
<< " interactions, found " << rdkTorsionInstanceVec.size()
<< "\n";
errorMsg += diff.str();
}
for (unsigned j = 0; (!error) & (j < rdkTorsionInstanceVec.size());
++j) {
error = ((rdkTorsionInstanceVec[j]->iAtomType !=
refTorsionInstanceVec[j]->iAtomType) ||
(rdkTorsionInstanceVec[j]->jAtomType !=
refTorsionInstanceVec[j]->jAtomType) ||
(rdkTorsionInstanceVec[j]->kAtomType !=
refTorsionInstanceVec[j]->kAtomType) ||
(rdkTorsionInstanceVec[j]->lAtomType !=
refTorsionInstanceVec[j]->lAtomType) ||
(rdkTorsionInstanceVec[j]->ffType !=
refTorsionInstanceVec[j]->ffType) ||
(fabs(rdkTorsionInstanceVec[j]->V1 -
refTorsionInstanceVec[j]->V1) > FCON_TOLERANCE) ||
(fabs(rdkTorsionInstanceVec[j]->V2 -
refTorsionInstanceVec[j]->V2) > FCON_TOLERANCE) ||
(fabs(rdkTorsionInstanceVec[j]->V3 -
refTorsionInstanceVec[j]->V3) > FCON_TOLERANCE));
if (error) {
failed = true;
bool haveLine;
haveLine = getLineByNum(rdkFStream, rdkPos,
rdkTorsionInstanceVec[j]->idx, rdkLine);
if (haveLine) {
haveLine =
getLineByNum(refFStream, refPos,
refTorsionInstanceVec[j]->idx, refLine);
}
std::stringstream diff;
diff << "Torsional: found a difference\n";
if (haveLine) {
diff << "Expected:\n"
<< rdkLine << "\nFound:\n"
<< refLine << "\n";
}
errorMsg += diff.str();
}
}
for (auto &ti : rdkTorsionInstanceVec) {
delete ti;
}
for (auto &ti : refTorsionInstanceVec) {
delete ti;
}
error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Torsional: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
}
// --------------------------------------------------
// V A N D E R W A A L S
// E L E C T R O S T A T I C
// --------------------------------------------------
corruptedMsg = ": Van der Waals: corrupted input data\n";
found = fgrep(rdkFStream, "V A N D E R W A A L S");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "TOTAL VAN DER WAALS ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
if (!refHaveVdW) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream rdkVdWStringStream(energyString);
rdkVdWStringStream >> skip >> skip >> skip >> skip >> skip >> skip >>
rdkEnergy;
corruptedMsg = ": Electrostatic: corrupted input data";
found = fgrep(rdkFStream, "E L E C T R O S T A T I C");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "----------------");
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
found = fgrep(rdkFStream, "TOTAL ELECTROSTATIC ENERGY", energyString);
if (!found) {
errorMsg += rdk + corruptedMsg;
break;
}
if (!refHaveEle) {
errorMsg += ref + corruptedMsg;
break;
}
std::stringstream rdkEleStringStream(energyString);
rdkEleStringStream >> skip >> skip >> skip >> skip >> rdkEleEnergy;
error = (fabs(rdkEnergy - refVdWEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Van der Waals: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refVdWEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
error = (fabs(rdkEleEnergy - refEleEnergy) > ENERGY_TOLERANCE);
if (error && checkEnergy) {
failed = true;
std::stringstream diff;
diff << "Electrostatic: energies differ\n"
"Expected "
<< std::fixed << std::setprecision(4) << refEleEnergy
<< ", found " << std::fixed << std::setprecision(4)
<< rdkEnergy << "\n";
errorMsg += diff.str();
}
if (failed) {
std::cerr << nameArray[ii] << "\n\n" << errorMsg << std::endl;
++n_failed;
}
}
if (!found) {
if (!failed) {
std::cerr << nameArray[ii] << "\n\n";
}
std::cerr << errorMsg << std::endl;
} else {
unsigned int n_passed = nameArray.size() - n_failed;
std::cerr << (n_failed ? "" : "All ") << n_passed << " test"
<< ((n_passed == 1) ? "" : "s") << " passed" << std::endl;
if (n_failed) {
std::cerr << n_failed << " test" << ((n_failed == 1) ? "" : "s")
<< " failed" << std::endl;
}
}
if (n_failed) {
testFailure = true;
}
rdkFStream.close();
refFStream.close();
firstTest = false;
}
}
}
TEST_ASSERT(!testFailure);
std::cerr << " done" << std::endl;
return 0;
}
void testMMFFAllConstraints() {
std::cerr << "-------------------------------------" << std::endl;
std::cerr << "Unit tests for all MMFF constraint terms." << std::endl;
std::string molBlock =
"butane\n"
" RDKit 3D\n"
"butane\n"
" 17 16 0 0 0 0 0 0 0 0999 V2000\n"
" 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1 2 1 0 0 0 0\n"
" 1 15 1 0 0 0 0\n"
" 1 16 1 0 0 0 0\n"
" 1 17 1 0 0 0 0\n"
" 2 3 1 0 0 0 0\n"
" 2 4 1 0 0 0 0\n"
" 2 14 1 0 0 0 0\n"
" 4 5 1 0 0 0 0\n"
" 4 6 1 0 0 0 0\n"
" 4 7 1 0 0 0 0\n"
" 7 8 1 0 0 0 0\n"
" 7 9 1 0 0 0 0\n"
" 7 13 1 0 0 0 0\n"
" 9 10 1 0 0 0 0\n"
" 9 11 1 0 0 0 0\n"
" 9 12 1 0 0 0 0\n"
"M END\n";
RDKit::RWMol *mol;
ForceFields::ForceField *field;
MMFF::MMFFMolProperties *mmffMolProperties;
// distance constraints
ForceFields::DistanceConstraintContribs *dc;
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
MolTransforms::setBondLength(mol->getConformer(), 1, 3, 2.0);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
dc = new ForceFields::DistanceConstraintContribs(field);
dc->addContrib(1, 3, 2.0, 2.0, 1.0e5);
field->contribs().emplace_back(dc);
field->minimize();
TEST_ASSERT(RDKit::feq(
MolTransforms::getBondLength(mol->getConformer(), 1, 3), 2.0, 0.1));
delete field;
delete mmffMolProperties;
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
field->initialize();
dc = new ForceFields::DistanceConstraintContribs(field);
dc->addContrib(1, 3, true, -0.2, 0.2, 1.0e5);
field->contribs().emplace_back(dc);
field->minimize();
TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.79);
delete field;
delete mmffMolProperties;
delete mol;
// angle constraints
ForceFields::AngleConstraintContribs *ac;
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 90.0);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
ac = new ForceFields::AngleConstraintContribs(field);
ac->addContrib(1, 3, 6, 90.0, 90.0, 100.0);
field->contribs().emplace_back(ac);
field->minimize();
TEST_ASSERT(RDKit::feq(
MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 90.0, 0.5));
delete field;
delete mmffMolProperties;
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
field->initialize();
ac = new ForceFields::AngleConstraintContribs(field);
ac->addContrib(1, 3, 6, true, -10.0, 10.0, 100.0);
field->contribs().emplace_back(ac);
field->minimize();
TEST_ASSERT(RDKit::feq(
MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 100.0, 0.5));
delete field;
delete mmffMolProperties;
delete mol;
// torsion constraints
ForceFields::MMFF::TorsionConstraintContrib *tc;
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, 15.0);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, 10.0,
10.0, 100.0);
field->contribs().push_back(ForceFields::ContribPtr(tc));
field->minimize();
std::cerr << MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8)
<< std::endl;
TEST_ASSERT(
RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
10.0, 0.5));
delete field;
delete mmffMolProperties;
MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -30.0);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
field->initialize();
tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, true,
-1.0, 1.0, 100.0);
field->contribs().push_back(ForceFields::ContribPtr(tc));
field->minimize();
TEST_ASSERT(
RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
-30.0, 1.5));
delete field;
delete mmffMolProperties;
MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -10.0);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
field->initialize();
tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, false,
-10.0, -8.0, 100.0);
field->contribs().push_back(ForceFields::ContribPtr(tc));
field->minimize(500);
TEST_ASSERT(
RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
-9.0, 2.0));
delete field;
delete mmffMolProperties;
delete mol;
// position constraints
ForceFields::MMFF::PositionConstraintContrib *pc;
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
RDGeom::Point3D p = mol->getConformer().getAtomPos(1);
pc = new ForceFields::MMFF::PositionConstraintContrib(field, 1, 0.3, 1.0e5);
field->contribs().push_back(ForceFields::ContribPtr(pc));
field->minimize();
RDGeom::Point3D q = mol->getConformer().getAtomPos(1);
TEST_ASSERT((p - q).length() < 0.3);
delete field;
delete mmffMolProperties;
delete mol;
// fixed atoms
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
field = RDKit::MMFF::constructForceField(*mol);
TEST_ASSERT(field);
field->initialize();
RDGeom::Point3D fp = mol->getConformer().getAtomPos(1);
field->fixedPoints().push_back(1);
field->minimize();
RDGeom::Point3D fq = mol->getConformer().getAtomPos(1);
TEST_ASSERT((fp - fq).length() < 0.01);
delete field;
delete mol;
std::cerr << " done" << std::endl;
}
void testMMFFTorsionConstraints() {
std::cerr << "-------------------------------------" << std::endl;
std::cerr << "Unit test for MMFF TorsionConstraint close to 0 degrees."
<< std::endl;
std::string molBlock = R"(
RDKit 3D
11 10 0 0 0 0 0 0 0 0999 V2000
-1.6091 -0.3348 -0.0457 N 0 0 0 0 0 0 0 0 0 0 0 0
-2.3935 0.2978 -0.1061 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.6708 -0.8686 0.8103 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.3560 0.4358 -0.0518 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.3583 1.1207 0.7968 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.2978 1.0122 -0.9759 H 0 0 0 0 0 0 0 0 0 0 0 0
0.8573 -0.4893 0.0451 C 0 0 0 0 0 0 0 0 0 0 0 0
0.8947 -1.1355 -0.8327 H 0 0 0 0 0 0 0 0 0 0 0 0
0.7774 -1.1052 0.9415 H 0 0 0 0 0 0 0 0 0 0 0 0
2.0396 0.2763 0.1127 O 0 0 0 0 0 0 0 0 0 0 0 0
2.1165 0.7905 -0.6941 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0 0 0 0
1 3 1 0 0 0 0
1 4 1 0 0 0 0
4 5 1 0 0 0 0
4 6 1 0 0 0 0
4 7 1 0 0 0 0
7 8 1 0 0 0 0
7 9 1 0 0 0 0
7 10 1 0 0 0 0
10 11 1 0 0 0 0
M END
)";
RDKit::RWMol *mol;
// torsion constraints
mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
std::vector<double> e;
for (double d : {0.0001, 0.0}) {
MolTransforms::setDihedralDeg(mol->getConformer(), 0, 3, 6, 9, d);
MMFF::MMFFMolProperties mmffMolProperties(*mol);
TEST_ASSERT(mmffMolProperties.isValid());
ForceFields::ForceField *field =
RDKit::MMFF::constructForceField(*mol, &mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
auto *tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 0, 3, 6,
9, d, d, 1.0e6);
field->contribs().push_back(ForceFields::ContribPtr(tc));
field->minimize();
e.push_back(field->calcEnergy());
TEST_ASSERT(RDKit::feq(
MolTransforms::getDihedralDeg(mol->getConformer(), 0, 3, 6, 9), d,
0.5));
delete field;
}
TEST_ASSERT(RDKit::feq(e[0], e[1], 0.5));
delete mol;
std::cerr << " done" << std::endl;
}
void testMMFFCopy() {
std::cerr << "-------------------------------------" << std::endl;
std::cerr << "Unit tests for copying MMFF ForceFields." << std::endl;
std::string molBlock =
"butane\n"
" RDKit 3D\n"
"butane\n"
" 17 16 0 0 0 0 0 0 0 0999 V2000\n"
" 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
" 1 2 1 0 0 0 0\n"
" 1 15 1 0 0 0 0\n"
" 1 16 1 0 0 0 0\n"
" 1 17 1 0 0 0 0\n"
" 2 3 1 0 0 0 0\n"
" 2 4 1 0 0 0 0\n"
" 2 14 1 0 0 0 0\n"
" 4 5 1 0 0 0 0\n"
" 4 6 1 0 0 0 0\n"
" 4 7 1 0 0 0 0\n"
" 7 8 1 0 0 0 0\n"
" 7 9 1 0 0 0 0\n"
" 7 13 1 0 0 0 0\n"
" 9 10 1 0 0 0 0\n"
" 9 11 1 0 0 0 0\n"
" 9 12 1 0 0 0 0\n"
"M END\n";
{
RDKit::RWMol *mol = RDKit::MolBlockToMol(molBlock, true, false);
TEST_ASSERT(mol);
auto *cmol = new RDKit::RWMol(*mol);
TEST_ASSERT(cmol);
MMFF::MMFFMolProperties *mmffMolProperties =
new MMFF::MMFFMolProperties(*mol);
TEST_ASSERT(mmffMolProperties);
TEST_ASSERT(mmffMolProperties->isValid());
ForceFields::ForceField *field =
RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
TEST_ASSERT(field);
field->initialize();
auto *dc = new ForceFields::DistanceConstraintContribs(field);
dc->addContrib(1, 3, 2.0, 2.0, 1.0e5);
field->contribs().emplace_back(dc);
field->minimize();
TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.99);
auto *cfield = new ForceFields::ForceField(*field);
cfield->positions().clear();
for (unsigned int i = 0; i < cmol->getNumAtoms(); i++) {
cfield->positions().push_back(&cmol->getConformer().getAtomPos(i));
}
cfield->initialize();
cfield->minimize();
TEST_ASSERT(MolTransforms::getBondLength(cmol->getConformer(), 1, 3) >
1.99);
TEST_ASSERT(RDKit::feq(field->calcEnergy(), cfield->calcEnergy()));
const RDKit::Conformer &conf = mol->getConformer();
const RDKit::Conformer &cconf = cmol->getConformer();
for (unsigned int i = 0; i < mol->getNumAtoms(); i++) {
RDGeom::Point3D p = conf.getAtomPos(i);
RDGeom::Point3D cp = cconf.getAtomPos(i);
TEST_ASSERT(RDKit::feq(p.x, cp.x));
TEST_ASSERT(RDKit::feq(p.y, cp.y));
TEST_ASSERT(RDKit::feq(p.z, cp.z));
}
delete field;
delete cfield;
delete mmffMolProperties;
delete mol;
delete cmol;
}
std::cerr << " done" << std::endl;
}
void testMMFFButaneScan() {
std::cerr << "-------------------------------------" << std::endl;
std::cerr << "Unit test for MMFF butane scan." << std::endl;
auto molblock = R"(
RDKit 3D
14 13 0 0 0 0 0 0 0 0999 V2000
-1.5575 0.5684 -0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.6987 -0.6064 0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
0.6987 -0.6064 -0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
1.5575 0.5684 0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.6308 0.6129 -1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.5700 0.4662 0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.1524 1.5190 0.2272 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.2059 -1.5359 0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
-0.6215 -0.6099 1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
0.6215 -0.6099 -1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
1.2059 -1.5359 -0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
1.6308 0.6129 1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
1.1524 1.5190 -0.2271 H 0 0 0 0 0 0 0 0 0 0 0 0
2.5700 0.4662 -0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0
2 3 1 0
3 4 1 0
1 5 1 0
1 6 1 0
1 7 1 0
2 8 1 0
2 9 1 0
3 10 1 0
3 11 1 0
4 12 1 0
4 13 1 0
4 14 1 0
M END
)";
// torsion constraints
std::unique_ptr<ROMol> mol(MolBlockToMol(molblock, true, false));
TEST_ASSERT(mol);
std::vector<std::pair<double, double>> diheEnergyVect;
std::vector<std::pair<double, double>> expectedDiheEnergyVect{
{-179.999, -5.07597}, {-175.1, -5.01299}, {-170.1, -4.82267},
{-165.1, -4.51646}, {-160.1, -4.11308}, {-155.101, -3.63773},
{-150.101, -3.12109}, {-145.101, -2.59801}, {-140.1, -2.10575},
{-135.1, -1.68158}, {-130.1, -1.35952}, {-125.1, -1.16635},
{-119.9, -1.11879}, {-114.9, -1.22194}, {-109.9, -1.45802},
{-104.9, -1.80199}, {-99.8996, -2.22002}, {-94.8995, -2.67302},
{-89.8996, -3.12055}, {-84.8996, -3.5255}, {-79.8997, -3.85916},
{-74.8998, -4.10358}, {-69.8999, -4.24974}, {-65.1, -4.29368},
{-60.1001, -4.23737}, {-55.1002, -4.0775}, {-50.1003, -3.81811},
{-45.1004, -3.46765}, {-40.1005, -3.0395}, {-35.1005, -2.55211},
{-30.1005, -2.02873}, {-25.1005, -1.49709}, {-20.1005, -0.988743},
{-15.1004, -0.538157}, {-10.1003, -0.180672}, {-5.10016, 0.051218},
{-0.100003, 0.13365}, {5.10016, 0.051218}, {10.1003, -0.180672},
{15.1004, -0.538157}, {20.1005, -0.988743}, {25.1005, -1.49709},
{30.1005, -2.02873}, {35.1005, -2.55211}, {40.1005, -3.0395},
{45.1004, -3.46765}, {50.1003, -3.81811}, {55.1002, -4.0775},
{60.1001, -4.23737}, {65.1, -4.29368}, {69.8999, -4.24974},
{74.8998, -4.10358}, {79.8997, -3.85916}, {84.8996, -3.5255},
{89.8996, -3.12055}, {94.8995, -2.67302}, {99.8996, -2.22002},
{104.9, -1.80199}, {109.9, -1.45802}, {114.9, -1.22194},
{119.9, -1.11879}, {125.1, -1.16635}, {130.1, -1.35952},
{135.1, -1.68158}, {140.1, -2.10575}, {145.101, -2.59801},
{150.101, -3.12109}, {155.101, -3.63773}, {160.1, -4.11308},
{165.1, -4.51646}, {170.1, -4.82267}, {175.1, -5.01299}};
MMFF::MMFFMolProperties mmffMolProperties(*mol);
TEST_ASSERT(mmffMolProperties.isValid());
std::unique_ptr<ForceFields::ForceField> globalFF(
RDKit::MMFF::constructForceField(*mol, &mmffMolProperties));
TEST_ASSERT(globalFF);
globalFF->initialize();
std::unique_ptr<ROMol> mol2(new ROMol(*mol));
const int start = -180;
const int end = 180;
const int incr = 5;
MolTransforms::setDihedralDeg(mol2->getConformer(), 0, 1, 2, 3, start);
for (int diheIn = start; diheIn < end; diheIn += incr) {
std::unique_ptr<ForceFields::ForceField> localFF(
RDKit::MMFF::constructForceField(*mol2, &mmffMolProperties));
TEST_ASSERT(localFF);
localFF->initialize();
auto *tc = new ForceFields::MMFF::TorsionConstraintContrib(
localFF.get(), 0, 1, 2, 3, static_cast<double>(diheIn) - 0.1,
static_cast<double>(diheIn) + 0.1, 100.0);
localFF->contribs().push_back(ForceFields::ContribPtr(tc));
TEST_ASSERT(localFF->minimize(10000) == 0);
std::vector<double> pos(3 * localFF->numPoints());
TEST_ASSERT(localFF->numPoints() == globalFF->numPoints());
auto posns = localFF->positions();
for (unsigned int i = 0; i < localFF->numPoints(); ++i) {
for (unsigned int j = 0; j < 3; ++j) {
pos[i * 3 + j] = (*static_cast<RDGeom::Point3D *>(posns.at(i)))[j];
}
}
auto diheOut =
MolTransforms::getDihedralDeg(mol2->getConformer(), 0, 1, 2, 3);
auto energy = globalFF->calcEnergy(pos.data());
diheEnergyVect.emplace_back(diheOut, energy);
}
for (size_t i = 0; i < diheEnergyVect.size(); ++i) {
TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).first,
expectedDiheEnergyVect.at(i).first, 1.0));
TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).second,
expectedDiheEnergyVect.at(i).second, 0.2));
}
std::cerr << " done" << std::endl;
}
int main(int argc, char *argv[]) {
if (!mmffValidationSuite(argc, argv)) {
testMMFFAllConstraints();
testMMFFTorsionConstraints();
testMMFFCopy();
testMMFFButaneScan();
}
return 0;
}
Reference in New Issue View Git Blame Copy Permalink