1 // $Id$
2 //
3 // Copyright (C) 2013 Paolo Tosco
4 //
5 // @@ All Rights Reserved @@
6 // This file is part of the RDKit.
7 // The contents are covered by the terms of the BSD license
8 // which is included in the file license.txt, found at the root
9 // of the RDKit source tree.
10 //
11
12 #include <RDGeneral/test.h>
13 #include <GraphMol/RDKitBase.h>
14 #include <GraphMol/FileParsers/MolSupplier.h>
15 #include <GraphMol/FileParsers/MolWriters.h>
16 #include <GraphMol/FileParsers/FileParsers.h>
17 #include <ForceField/MMFF/Params.h>
18 #include <ForceField/MMFF/DistanceConstraint.h>
19 #include <ForceField/MMFF/AngleConstraint.h>
20 #include <ForceField/MMFF/TorsionConstraint.h>
21 #include <ForceField/MMFF/PositionConstraint.h>
22 #include <GraphMol/ForceFieldHelpers/MMFF/AtomTyper.h>
23 #include <GraphMol/ForceFieldHelpers/MMFF/Builder.h>
24 #include <GraphMol/MolOps.h>
25 #include <GraphMol/DistGeomHelpers/Embedder.h>
26 #include <GraphMol/MolTransforms/MolTransforms.h>
27 #include <ForceField/ForceField.h>
28 #include <iostream>
29 #include <iomanip>
30 #include <fstream>
31 #include <sstream>
32 #include <string>
33 #include "testMMFFForceField.h"
34
35 #define FCON_TOLERANCE 0.05
36 #define ENERGY_TOLERANCE 0.025
37
38 using namespace RDKit;
39
fexist(std::string filename)40 bool fexist(std::string filename) {
41 std::ifstream ifStream(filename.c_str());
42
43 return (ifStream ? true : false);
44 }
45
fgrep(std::fstream & fStream,std::string key,std::string & line)46 bool fgrep(std::fstream &fStream, std::string key, std::string &line) {
47 bool found = false;
48 std::streampos current = fStream.tellg();
49 while ((!found) &&
50 (!(std::getline(fStream, line).rdstate() & std::ifstream::failbit))) {
51 found = (line.find(key) != std::string::npos);
52 }
53 if (!found) {
54 fStream.seekg(current);
55 }
56
57 return found;
58 }
59
fgrep(std::fstream & rdkFStream,std::string key)60 bool fgrep(std::fstream &rdkFStream, std::string key) {
61 std::string line;
62
63 return fgrep(rdkFStream, key, line);
64 }
65
getLineByNum(std::fstream & stream,std::streampos startPos,unsigned int n,std::string & line)66 bool getLineByNum(std::fstream &stream, std::streampos startPos, unsigned int n,
67 std::string &line) {
68 std::streampos current = stream.tellg();
69 stream.seekg(startPos);
70 bool fail;
71 unsigned int i = 0;
72 while ((i <= n) && (!(fail = std::getline(stream, line).rdstate() &
73 std::ifstream::failbit))) {
74 ++i;
75 }
76 stream.seekg(current);
77
78 return (!fail);
79 }
80
sortBondStretchInstanceVec(BondStretchInstance * a,BondStretchInstance * b)81 bool sortBondStretchInstanceVec(BondStretchInstance *a,
82 BondStretchInstance *b) {
83 bool crit = false;
84
85 if (a->iAtomType == b->iAtomType) {
86 if (a->jAtomType == b->jAtomType) {
87 crit = (a->ffType < b->ffType);
88 } else {
89 crit = (a->jAtomType < b->jAtomType);
90 }
91 } else {
92 crit = (a->iAtomType < b->iAtomType);
93 }
94
95 return crit;
96 }
97
sortAngleBendInstanceVec(AngleBendInstance * a,AngleBendInstance * b)98 bool sortAngleBendInstanceVec(AngleBendInstance *a, AngleBendInstance *b) {
99 bool crit = false;
100
101 if (a->jAtomType == b->jAtomType) {
102 if (a->iAtomType == b->iAtomType) {
103 if (a->kAtomType == b->kAtomType) {
104 crit = (a->ffType < b->ffType);
105 } else {
106 crit = (a->kAtomType < b->kAtomType);
107 }
108 } else {
109 crit = (a->iAtomType < b->iAtomType);
110 }
111 } else {
112 crit = (a->jAtomType < b->jAtomType);
113 }
114
115 return crit;
116 }
117
sortStretchBendInstanceVec(StretchBendInstance * a,StretchBendInstance * b)118 bool sortStretchBendInstanceVec(StretchBendInstance *a,
119 StretchBendInstance *b) {
120 bool crit = false;
121
122 if (a->jAtomType == b->jAtomType) {
123 if (a->iAtomType == b->iAtomType) {
124 if (a->kAtomType == b->kAtomType) {
125 crit = (a->ffType < b->ffType);
126 } else {
127 crit = (a->kAtomType < b->kAtomType);
128 }
129 } else {
130 crit = (a->iAtomType < b->iAtomType);
131 }
132 } else {
133 crit = (a->jAtomType < b->jAtomType);
134 }
135
136 return crit;
137 }
138
sortOopBendInstanceVec(OopBendInstance * a,OopBendInstance * b)139 bool sortOopBendInstanceVec(OopBendInstance *a, OopBendInstance *b) {
140 bool crit = false;
141
142 if (a->jAtomType == b->jAtomType) {
143 if (a->iAtomType == b->iAtomType) {
144 if (a->kAtomType == b->kAtomType) {
145 crit = (a->lAtomType < b->lAtomType);
146 } else {
147 crit = (a->kAtomType < b->kAtomType);
148 }
149 } else {
150 crit = (a->iAtomType < b->iAtomType);
151 }
152 } else {
153 crit = (a->jAtomType < b->jAtomType);
154 }
155
156 return crit;
157 }
158
sortTorsionInstanceVec(TorsionInstance * a,TorsionInstance * b)159 bool sortTorsionInstanceVec(TorsionInstance *a, TorsionInstance *b) {
160 bool crit = false;
161
162 if (a->jAtomType == b->jAtomType) {
163 if (a->kAtomType == b->kAtomType) {
164 if (a->iAtomType == b->iAtomType) {
165 if (a->lAtomType == b->lAtomType) {
166 crit = (a->ffType < b->ffType);
167 } else {
168 crit = (a->lAtomType < b->lAtomType);
169 }
170 } else {
171 crit = (a->iAtomType < b->iAtomType);
172 }
173 } else {
174 crit = (a->kAtomType < b->kAtomType);
175 }
176 } else {
177 crit = (a->jAtomType < b->jAtomType);
178 }
179
180 return crit;
181 }
182
fixBondStretchInstance(BondStretchInstance * bondStretchInstance)183 void fixBondStretchInstance(BondStretchInstance *bondStretchInstance) {
184 if (bondStretchInstance->iAtomType > bondStretchInstance->jAtomType) {
185 unsigned int temp;
186 temp = bondStretchInstance->iAtomType;
187 bondStretchInstance->iAtomType = bondStretchInstance->jAtomType;
188 bondStretchInstance->jAtomType = temp;
189 }
190 }
191
fixAngleBendInstance(AngleBendInstance * angleBendInstance)192 void fixAngleBendInstance(AngleBendInstance *angleBendInstance) {
193 if (angleBendInstance->iAtomType > angleBendInstance->kAtomType) {
194 unsigned int temp;
195 temp = angleBendInstance->iAtomType;
196 angleBendInstance->iAtomType = angleBendInstance->kAtomType;
197 angleBendInstance->kAtomType = temp;
198 }
199 }
200
fixOopBendInstance(OopBendInstance * oopBendInstance)201 void fixOopBendInstance(OopBendInstance *oopBendInstance) {
202 if (oopBendInstance->iAtomType > oopBendInstance->kAtomType) {
203 unsigned int temp;
204 temp = oopBendInstance->iAtomType;
205 oopBendInstance->iAtomType = oopBendInstance->kAtomType;
206 oopBendInstance->kAtomType = temp;
207 }
208 }
209
fixTorsionInstance(TorsionInstance * torsionInstance)210 void fixTorsionInstance(TorsionInstance *torsionInstance) {
211 if (torsionInstance->jAtomType > torsionInstance->kAtomType) {
212 unsigned int temp;
213 temp = torsionInstance->jAtomType;
214 torsionInstance->jAtomType = torsionInstance->kAtomType;
215 torsionInstance->kAtomType = temp;
216 temp = torsionInstance->iAtomType;
217 torsionInstance->iAtomType = torsionInstance->lAtomType;
218 torsionInstance->lAtomType = temp;
219 } else if (torsionInstance->jAtomType == torsionInstance->kAtomType) {
220 unsigned int temp;
221 if (torsionInstance->iAtomType > torsionInstance->lAtomType) {
222 temp = torsionInstance->iAtomType;
223 torsionInstance->iAtomType = torsionInstance->lAtomType;
224 torsionInstance->lAtomType = temp;
225 }
226 }
227 }
228
mmffValidationSuite(int argc,char * argv[])229 int mmffValidationSuite(int argc, char *argv[]) {
230 std::string arg;
231 std::string ffVariant = "";
232 std::vector<std::string> ffVec;
233 std::string verbosity;
234 std::string molFile = "";
235 std::string molType = "";
236 std::vector<std::string> molFileVec;
237 std::vector<std::string> molTypeVec;
238 std::string rdk = "";
239 std::string ref = "";
240 std::streampos rdkPos;
241 std::streampos rdkCurrent;
242 std::streampos refPos;
243 std::streampos refCurrent;
244 int iarg = 1;
245 unsigned int n = 0;
246 bool error = false;
247 bool fullTest = false;
248 int inc = 4;
249 bool testFailure = false;
250 while (iarg < argc) {
251 arg = argv[iarg];
252 error = (arg.at(0) != '-');
253 if (error) {
254 break;
255 }
256 error = ((arg == "-h") || (arg == "--help"));
257 if (error) {
258 break;
259 } else if (arg == "-f") {
260 error = ((iarg + 1) >= argc);
261 if (error) {
262 break;
263 }
264 ffVariant = argv[iarg + 1];
265 error = ((ffVariant != "MMFF94") && (ffVariant != "MMFF94s"));
266 if (error) {
267 break;
268 }
269 ++iarg;
270 } else if ((arg == "-sdf") || (arg == "-smi")) {
271 molType = arg.substr(1);
272 if ((iarg + 1) < argc) {
273 molFile = argv[iarg + 1];
274 ++iarg;
275 }
276 } else if (arg == "-L") {
277 fullTest = true;
278 inc = 1;
279 } else if (arg == "-l") {
280 error = ((iarg + 1) >= argc);
281 if (error) {
282 break;
283 }
284 rdk = argv[iarg + 1];
285 ++iarg;
286 }
287 ++iarg;
288 }
289 if (error) {
290 std::cerr << "Usage: testMMFFForceField [-L] [-f {MMFF94|MMFF94s}] "
291 "[{-sdf [<sdf_file>] | -smi [<smiles_file>]}] [-l <log_file>]"
292 << std::endl;
293
294 return -1;
295 }
296 std::cerr << "-------------------------------------" << std::endl;
297 std::cerr << "Official MMFF validation suite." << std::endl;
298 std::string pathName = getenv("RDBASE");
299 pathName += "/Code/ForceField/MMFF/test_data/";
300 if (molFile != "") {
301 ffVariant = ((molFile.substr(0, 7) == "MMFF94s") ? "MMFF94s" : "MMFF94");
302 }
303 if (ffVariant == "") {
304 ffVec.push_back("MMFF94");
305 if (fullTest) {
306 ffVec.push_back("MMFF94s");
307 }
308 } else {
309 ffVec.push_back(ffVariant);
310 }
311 std::fstream rdkFStream;
312 std::fstream refFStream;
313 if (rdk == "") {
314 rdk = pathName + "testMMFFForceField.log";
315 }
316 bool firstTest = true;
317 if (molType == "") {
318 molTypeVec.push_back("sdf");
319 molTypeVec.push_back("smi");
320 } else {
321 molTypeVec.push_back(molType);
322 }
323 boost::logging::disable_logs("rdApp.error");
324 for (auto &molTypeIt : molTypeVec) {
325 bool checkEnergy = (molTypeIt == "sdf");
326 for (auto &ffIt : ffVec) {
327 ref = pathName + ffIt + "_reference.log";
328 molFileVec.clear();
329 if (molFile == "") {
330 molFileVec.push_back(pathName + ffIt + "_dative." + molTypeIt);
331 molFileVec.push_back(pathName + ffIt + "_hypervalent." + molTypeIt);
332 } else {
333 molFileVec.push_back(molFile);
334 }
335 for (auto &molFileIt : molFileVec) {
336 std::vector<ROMol *> molVec;
337 molVec.clear();
338 if (molTypeIt == "sdf") {
339 SDMolSupplier sdfSupplier(molFileIt, false, false);
340 for (size_t ii = 0; ii < sdfSupplier.length(); ++ii) {
341 molVec.push_back(sdfSupplier[ii]);
342 }
343 sdfSupplier.reset();
344 } else {
345 SmilesMolSupplier smiSupplier(molFileIt);
346 for (size_t ii = 0; ii < smiSupplier.length(); ++ii) {
347 molVec.push_back(smiSupplier[ii]);
348 }
349 smiSupplier.reset();
350 }
351 SDWriter *sdfWriter =
352 new SDWriter(molFileIt.substr(0, molFileIt.length() - 4) + "_min" +
353 ((molTypeIt == "smi") ? "_from_SMILES" : "") + ".sdf");
354 ROMol *mol;
355 std::string molName;
356 std::vector<std::string> nameArray;
357
358 rdkFStream.open(rdk.c_str(),
359 (firstTest ? std::fstream::out | std::fstream::binary
360 : (std::fstream::out | std::fstream::binary |
361 std::fstream::app)));
362 std::string computingKey = ffIt + " energies for " + molFileIt;
363 std::cerr << std::endl
364 << "Computing " << computingKey << "..." << std::endl;
365 for (size_t ii = 0; ii < computingKey.length(); ++ii) {
366 rdkFStream << "*";
367 }
368 rdkFStream << std::endl << computingKey << std::endl;
369 for (size_t ii = 0; ii < computingKey.length(); ++ii) {
370 rdkFStream << "*";
371 }
372 rdkFStream << std::endl << std::endl;
373 for (size_t ii = 0; ii < molVec.size(); ii += inc) {
374 mol = molVec[ii];
375 if (molTypeIt == "smi") {
376 DGeomHelpers::EmbedMolecule(*mol);
377 MolOps::addHs((RWMol &)*mol, false, true);
378 }
379 MMFF::sanitizeMMFFMol((RWMol &)(*mol));
380 if (mol->hasProp(common_properties::_Name)) {
381 mol->getProp(common_properties::_Name, molName);
382 rdkFStream << molName << std::endl;
383 nameArray.push_back(molName);
384 }
385 auto *mmffMolProperties = new MMFF::MMFFMolProperties(
386 *mol, ffIt, MMFF::MMFF_VERBOSITY_HIGH, rdkFStream);
387 if ((!mmffMolProperties) || (!(mmffMolProperties->isValid()))) {
388 std::cerr << molName + ": error setting up force-field"
389 << std::endl;
390 continue;
391 } else {
392 ForceFields::ForceField *field = MMFF::constructForceField(
393 *mol, mmffMolProperties, 100.0, -1, false);
394 field->initialize();
395 field->minimize();
396 sdfWriter->write(*mol);
397 rdkFStream << "\nTOTAL MMFF ENERGY =" << std::right
398 << std::setw(16) << std::fixed << std::setprecision(4)
399 << field->calcEnergy() << "\n\n"
400 << std::endl;
401 delete field;
402 }
403 delete mmffMolProperties;
404 }
405 for (auto m : molVec) {
406 delete m;
407 }
408 sdfWriter->close();
409 delete sdfWriter;
410 rdkFStream.close();
411 std::cerr << "Checking against " << ref << "..." << std::endl;
412 rdkFStream.open(rdk.c_str(), std::fstream::in | std::fstream::binary);
413 refFStream.open(ref.c_str(), std::fstream::in | std::fstream::binary);
414
415 bool found = false;
416 std::string skip;
417 std::string errorMsg;
418 std::string corruptedMsg;
419 std::string energyString;
420 std::vector<BondStretchInstance *> rdkBondStretchInstanceVec;
421 std::vector<BondStretchInstance *> refBondStretchInstanceVec;
422 std::vector<AngleBendInstance *> rdkAngleBendInstanceVec;
423 std::vector<AngleBendInstance *> refAngleBendInstanceVec;
424 std::vector<StretchBendInstance *> rdkStretchBendInstanceVec;
425 std::vector<StretchBendInstance *> refStretchBendInstanceVec;
426 std::vector<OopBendInstance *> rdkOopBendInstanceVec;
427 std::vector<OopBendInstance *> refOopBendInstanceVec;
428 std::vector<TorsionInstance *> rdkTorsionInstanceVec;
429 std::vector<TorsionInstance *> refTorsionInstanceVec;
430 double rdkEnergy;
431 double rdkEleEnergy;
432 double refEnergy;
433 double refVdWEnergy;
434 double refEleEnergy;
435 unsigned int n_failed = 0;
436 bool failed = false;
437 rdkFStream.seekg(0);
438 fgrep(rdkFStream, computingKey);
439 size_t ii;
440 for (ii = 0; ii < nameArray.size(); ii += inc) {
441 error = false;
442 failed = false;
443 errorMsg = rdk + ": Molecule not found";
444 found = fgrep(rdkFStream, nameArray[ii]);
445 if (!found) {
446 break;
447 }
448 errorMsg = ref + ": Molecule not found";
449 found = fgrep(refFStream, nameArray[ii]);
450 if (!found) {
451 break;
452 }
453 refPos = refFStream.tellg();
454 found = false;
455 bool refHaveBondStretch = false;
456 bool refHaveAngleBend = false;
457 bool refHaveStretchBend = false;
458 bool refHaveOopBend = false;
459 bool refHaveTorsion = false;
460 bool refHaveVdW = false;
461 bool refHaveEle = false;
462 std::string rdkLine;
463 std::string refLine;
464 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
465 std::ifstream::failbit))) {
466 std::size_t occ = refLine.find("****");
467 found = (occ != std::string::npos);
468 if (!found) {
469 if (!refHaveVdW) {
470 occ = refLine.find("Net vdW");
471 refHaveVdW = (occ != std::string::npos);
472 if (refHaveVdW) {
473 std::stringstream refVdWStringStream(refLine);
474 refVdWStringStream >> skip >> skip >> refVdWEnergy;
475 }
476 }
477 if (!refHaveEle) {
478 occ = refLine.find("Electrostatic");
479 refHaveEle = (occ != std::string::npos);
480 if (refHaveEle) {
481 std::stringstream refEleStringStream(refLine);
482 refEleStringStream >> skip >> refEleEnergy;
483 }
484 }
485 if (!refHaveBondStretch) {
486 occ = refLine.find("B O N D S T R E T C H I N G");
487 refHaveBondStretch = (occ != std::string::npos);
488 }
489 if (!refHaveAngleBend) {
490 occ = refLine.find("A N G L E B E N D I N G");
491 refHaveAngleBend = (occ != std::string::npos);
492 }
493 if (!refHaveStretchBend) {
494 occ = refLine.find("S T R E T C H B E N D I N G");
495 refHaveStretchBend = (occ != std::string::npos);
496 }
497 if (!refHaveOopBend) {
498 occ = refLine.find("O U T - O F - P L A N E B E N D I N G");
499 refHaveOopBend = (occ != std::string::npos);
500 }
501 if (!refHaveTorsion) {
502 occ = refLine.find("T O R S I O N A L");
503 refHaveTorsion = (occ != std::string::npos);
504 }
505 }
506 }
507 refFStream.seekg(refPos);
508 errorMsg = "";
509
510 // --------------------------------------------------
511 // B O N D S T R E T C H I N G
512 // --------------------------------------------------
513 if (refHaveBondStretch) {
514 corruptedMsg = ": Bond stretching: corrupted input data\n";
515 found = fgrep(rdkFStream, "B O N D S T R E T C H I N G");
516 if (!found) {
517 errorMsg += rdk + corruptedMsg;
518 break;
519 }
520 found = fgrep(rdkFStream, "----------------");
521 if (!found) {
522 errorMsg += rdk + corruptedMsg;
523 break;
524 }
525 found = false;
526 rdkBondStretchInstanceVec.clear();
527 rdkPos = rdkFStream.tellg();
528 n = 0;
529 while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
530 std::ifstream::failbit))) {
531 found = (!(rdkLine.length()));
532 if (found) {
533 break;
534 }
535 std::stringstream rdkLineStream(rdkLine);
536 auto *bondStretchInstance = new BondStretchInstance();
537 bondStretchInstance->idx = n;
538 rdkLineStream >> skip >> skip >> skip >> skip >>
539 bondStretchInstance->iAtomType >>
540 bondStretchInstance->jAtomType >>
541 bondStretchInstance->ffType >> skip >> skip >> skip >> skip >>
542 bondStretchInstance->kb;
543 fixBondStretchInstance(bondStretchInstance);
544 rdkBondStretchInstanceVec.push_back(bondStretchInstance);
545 ++n;
546 }
547 if (!found) {
548 errorMsg += rdk + corruptedMsg;
549 break;
550 }
551 std::sort(rdkBondStretchInstanceVec.begin(),
552 rdkBondStretchInstanceVec.end(),
553 sortBondStretchInstanceVec);
554 found =
555 fgrep(rdkFStream, "TOTAL BOND STRETCH ENERGY", energyString);
556 if (!found) {
557 errorMsg += rdk + corruptedMsg;
558 break;
559 }
560 std::stringstream rdkBondStretchStringStream(energyString);
561 rdkBondStretchStringStream >> skip >> skip >> skip >> skip >>
562 skip >> rdkEnergy;
563 found = fgrep(refFStream, "B O N D S T R E T C H I N G");
564 if (!found) {
565 errorMsg += ref + corruptedMsg;
566 break;
567 }
568 found = fgrep(refFStream, "----------------");
569 if (!found) {
570 errorMsg += ref + corruptedMsg;
571 break;
572 }
573 found = false;
574 refBondStretchInstanceVec.clear();
575 refPos = refFStream.tellg();
576 n = 0;
577 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
578 std::ifstream::failbit))) {
579 found = (!(refLine.length()));
580 if (found) {
581 break;
582 }
583 std::stringstream refLineStream(refLine);
584 auto *bondStretchInstance = new BondStretchInstance();
585 bondStretchInstance->idx = n;
586 refLineStream >> skip >> skip >> skip >> skip >>
587 bondStretchInstance->iAtomType >>
588 bondStretchInstance->jAtomType >>
589 bondStretchInstance->ffType >> skip >> skip >> skip >> skip >>
590 bondStretchInstance->kb;
591 fixBondStretchInstance(bondStretchInstance);
592 refBondStretchInstanceVec.push_back(bondStretchInstance);
593 ++n;
594 }
595 if (!found) {
596 errorMsg += ref + corruptedMsg;
597 break;
598 }
599 std::sort(refBondStretchInstanceVec.begin(),
600 refBondStretchInstanceVec.end(),
601 sortBondStretchInstanceVec);
602 found = fgrep(refFStream, "TOTAL BOND STRAIN ENERGY", energyString);
603 if (!found) {
604 errorMsg += ref + corruptedMsg;
605 break;
606 }
607 std::stringstream refBondStretchStringStream(energyString);
608 refBondStretchStringStream >> skip >> skip >> skip >> skip >>
609 skip >> refEnergy;
610 error = (rdkBondStretchInstanceVec.size() !=
611 refBondStretchInstanceVec.size());
612 if (error) {
613 failed = true;
614 std::stringstream diff;
615 diff << "Bond stretching: expected "
616 << refBondStretchInstanceVec.size()
617 << " interactions, found "
618 << rdkBondStretchInstanceVec.size() << "\n";
619 errorMsg += diff.str();
620 }
621 for (unsigned j = 0;
622 (!error) & (j < rdkBondStretchInstanceVec.size()); ++j) {
623 error =
624 ((rdkBondStretchInstanceVec[j]->iAtomType !=
625 refBondStretchInstanceVec[j]->iAtomType) ||
626 (rdkBondStretchInstanceVec[j]->jAtomType !=
627 refBondStretchInstanceVec[j]->jAtomType) ||
628 (rdkBondStretchInstanceVec[j]->ffType !=
629 refBondStretchInstanceVec[j]->ffType) ||
630 (fabs(rdkBondStretchInstanceVec[j]->kb -
631 refBondStretchInstanceVec[j]->kb) > FCON_TOLERANCE));
632 if (error) {
633 failed = true;
634 bool haveLine;
635 haveLine =
636 getLineByNum(rdkFStream, rdkPos,
637 rdkBondStretchInstanceVec[j]->idx, rdkLine);
638 if (haveLine) {
639 haveLine =
640 getLineByNum(refFStream, refPos,
641 refBondStretchInstanceVec[j]->idx, refLine);
642 }
643 std::stringstream diff;
644 diff << "Bond stretching: found a difference\n";
645 if (haveLine) {
646 diff << "Expected:\n"
647 << rdkLine << "\nFound:\n"
648 << refLine << "\n";
649 }
650 errorMsg += diff.str();
651 }
652 delete rdkBondStretchInstanceVec[j];
653 delete refBondStretchInstanceVec[j];
654 }
655 error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
656 if (error && checkEnergy) {
657 failed = true;
658 std::stringstream diff;
659 diff << "Bond stretching: energies differ\n"
660 "Expected "
661 << std::fixed << std::setprecision(4) << refEnergy
662 << ", found " << std::fixed << std::setprecision(4)
663 << rdkEnergy << "\n";
664 errorMsg += diff.str();
665 }
666 }
667
668 // --------------------------------------------------
669 // A N G L E B E N D I N G
670 // --------------------------------------------------
671 if (refHaveAngleBend) {
672 corruptedMsg = ": Angle bending: corrupted input data\n";
673 found = fgrep(rdkFStream, "A N G L E B E N D I N G");
674 if (!found) {
675 errorMsg += rdk + corruptedMsg;
676 break;
677 }
678 found = fgrep(rdkFStream, "----------------");
679 if (!found) {
680 errorMsg += rdk + corruptedMsg;
681 break;
682 }
683 found = false;
684 rdkAngleBendInstanceVec.clear();
685 rdkPos = rdkFStream.tellg();
686 n = 0;
687 while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
688 std::ifstream::failbit))) {
689 found = (!(rdkLine.length()));
690 if (found) {
691 break;
692 }
693 std::stringstream rdkLineStream(rdkLine);
694 auto *angleBendInstance = new AngleBendInstance();
695 angleBendInstance->idx = n;
696 rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
697 angleBendInstance->iAtomType >>
698 angleBendInstance->jAtomType >>
699 angleBendInstance->kAtomType >> angleBendInstance->ffType >>
700 skip >> skip >> skip >> skip >> angleBendInstance->ka;
701 fixAngleBendInstance(angleBendInstance);
702 rdkAngleBendInstanceVec.push_back(angleBendInstance);
703 ++n;
704 }
705 if (!found) {
706 errorMsg += rdk + corruptedMsg;
707 break;
708 }
709 std::sort(rdkAngleBendInstanceVec.begin(),
710 rdkAngleBendInstanceVec.end(), sortAngleBendInstanceVec);
711 found = fgrep(rdkFStream, "TOTAL ANGLE BEND ENERGY", energyString);
712 if (!found) {
713 errorMsg += rdk + corruptedMsg;
714 break;
715 }
716 std::stringstream rdkAngleBendStringStream(energyString);
717 rdkAngleBendStringStream >> skip >> skip >> skip >> skip >> skip >>
718 rdkEnergy;
719 found = fgrep(refFStream, "A N G L E B E N D I N G");
720 if (!found) {
721 errorMsg += ref + corruptedMsg;
722 break;
723 }
724 found = fgrep(refFStream, "----------------");
725 if (!found) {
726 errorMsg += ref + corruptedMsg;
727 break;
728 }
729 found = false;
730 refAngleBendInstanceVec.clear();
731 refPos = refFStream.tellg();
732 n = 0;
733 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
734 std::ifstream::failbit))) {
735 found = (!(refLine.length()));
736 if (found) {
737 break;
738 }
739 std::stringstream refLineStream(refLine);
740 auto *angleBendInstance = new AngleBendInstance();
741 angleBendInstance->idx = n;
742 refLineStream >> skip >> skip >> skip >> skip >>
743 angleBendInstance->iAtomType >>
744 angleBendInstance->jAtomType >>
745 angleBendInstance->kAtomType >> angleBendInstance->ffType >>
746 skip >> skip >> skip >> skip >> angleBendInstance->ka;
747 fixAngleBendInstance(angleBendInstance);
748 refAngleBendInstanceVec.push_back(angleBendInstance);
749 ++n;
750 }
751 if (!found) {
752 errorMsg += ref + corruptedMsg;
753 break;
754 }
755 std::sort(refAngleBendInstanceVec.begin(),
756 refAngleBendInstanceVec.end(), sortAngleBendInstanceVec);
757 found =
758 fgrep(refFStream, "TOTAL ANGLE STRAIN ENERGY", energyString);
759 if (!found) {
760 errorMsg += ref + corruptedMsg;
761 break;
762 }
763 std::stringstream refAngleBendStringStream(energyString);
764 refAngleBendStringStream >> skip >> skip >> skip >> skip >> skip >>
765 refEnergy;
766 error = (rdkAngleBendInstanceVec.size() !=
767 refAngleBendInstanceVec.size());
768 if (error) {
769 failed = true;
770 std::stringstream diff;
771 diff << "Angle bending: expected "
772 << refAngleBendInstanceVec.size() << " interactions, found "
773 << rdkAngleBendInstanceVec.size() << "\n";
774 errorMsg += diff.str();
775 }
776 for (unsigned j = 0;
777 (!error) & (j < rdkAngleBendInstanceVec.size()); ++j) {
778 error = ((rdkAngleBendInstanceVec[j]->iAtomType !=
779 refAngleBendInstanceVec[j]->iAtomType) ||
780 (rdkAngleBendInstanceVec[j]->jAtomType !=
781 refAngleBendInstanceVec[j]->jAtomType) ||
782 (rdkAngleBendInstanceVec[j]->kAtomType !=
783 refAngleBendInstanceVec[j]->kAtomType) ||
784 (rdkAngleBendInstanceVec[j]->ffType !=
785 refAngleBendInstanceVec[j]->ffType) ||
786 (fabs(rdkAngleBendInstanceVec[j]->ka -
787 refAngleBendInstanceVec[j]->ka) > FCON_TOLERANCE));
788 if (error) {
789 failed = true;
790 bool haveLine;
791 haveLine =
792 getLineByNum(rdkFStream, rdkPos,
793 rdkAngleBendInstanceVec[j]->idx, rdkLine);
794 if (haveLine) {
795 haveLine =
796 getLineByNum(refFStream, refPos,
797 refAngleBendInstanceVec[j]->idx, refLine);
798 }
799 std::stringstream diff;
800 diff << "Angle bending: found a difference\n";
801 if (haveLine) {
802 diff << "Expected:\n"
803 << rdkLine << "\nFound:\n"
804 << refLine << "\n";
805 }
806 errorMsg += diff.str();
807 }
808 delete rdkAngleBendInstanceVec[j];
809 delete refAngleBendInstanceVec[j];
810 }
811 error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
812 if (error && checkEnergy) {
813 failed = true;
814 std::stringstream diff;
815 diff << "Angle bending: energies differ\n"
816 "Expected "
817 << std::fixed << std::setprecision(4) << refEnergy
818 << ", found " << std::fixed << std::setprecision(4)
819 << rdkEnergy << "\n";
820 errorMsg += diff.str();
821 }
822 }
823
824 // --------------------------------------------------
825 // S T R E T C H B E N D I N G
826 // --------------------------------------------------
827 if (refHaveStretchBend) {
828 corruptedMsg = ": Stretch-bending: corrupted input data\n";
829 found = fgrep(rdkFStream, "S T R E T C H B E N D I N G");
830 if (!found) {
831 errorMsg += rdk + corruptedMsg;
832 break;
833 }
834 found = fgrep(rdkFStream, "----------------");
835 if (!found) {
836 errorMsg += rdk + corruptedMsg;
837 break;
838 }
839 found = false;
840 rdkStretchBendInstanceVec.clear();
841 rdkPos = rdkFStream.tellg();
842 n = 0;
843 while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
844 std::ifstream::failbit))) {
845 found = (!(rdkLine.length()));
846 if (found) {
847 break;
848 }
849 std::stringstream rdkLineStream(rdkLine);
850 auto *stretchBendInstance = new StretchBendInstance();
851 stretchBendInstance->idx = n;
852 rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
853 stretchBendInstance->iAtomType >>
854 stretchBendInstance->jAtomType >>
855 stretchBendInstance->kAtomType >>
856 stretchBendInstance->ffType >> skip >> skip >> skip >> skip >>
857 skip >> stretchBendInstance->kba;
858 rdkStretchBendInstanceVec.push_back(stretchBendInstance);
859 ++n;
860 }
861 if (!found) {
862 errorMsg += rdk + corruptedMsg;
863 break;
864 }
865 std::sort(rdkStretchBendInstanceVec.begin(),
866 rdkStretchBendInstanceVec.end(),
867 sortStretchBendInstanceVec);
868 found =
869 fgrep(rdkFStream, "TOTAL STRETCH-BEND ENERGY", energyString);
870 if (!found) {
871 errorMsg += rdk + corruptedMsg;
872 break;
873 }
874 std::stringstream rdkStretchBendStringStream(energyString);
875 rdkStretchBendStringStream >> skip >> skip >> skip >> skip >>
876 rdkEnergy;
877 found = fgrep(refFStream, "S T R E T C H B E N D I N G");
878 if (!found) {
879 errorMsg += ref + corruptedMsg;
880 break;
881 }
882 found = fgrep(refFStream, "----------------");
883 if (!found) {
884 errorMsg += ref + corruptedMsg;
885 break;
886 }
887 found = false;
888 refStretchBendInstanceVec.clear();
889 refPos = refFStream.tellg();
890 n = 0;
891 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
892 std::ifstream::failbit))) {
893 found = (!(refLine.length()));
894 if (found) {
895 break;
896 }
897 std::stringstream refLineStream(refLine);
898 auto *stretchBendInstance = new StretchBendInstance();
899 stretchBendInstance->idx = n;
900 refLineStream >> skip >> skip >> skip >> skip >>
901 stretchBendInstance->iAtomType >>
902 stretchBendInstance->jAtomType >>
903 stretchBendInstance->kAtomType >>
904 stretchBendInstance->ffType >> skip >> skip >> skip >> skip >>
905 stretchBendInstance->kba;
906 refStretchBendInstanceVec.push_back(stretchBendInstance);
907 ++n;
908 }
909 if (!found) {
910 errorMsg += ref + corruptedMsg;
911 break;
912 }
913 std::sort(refStretchBendInstanceVec.begin(),
914 refStretchBendInstanceVec.end(),
915 sortStretchBendInstanceVec);
916 found = fgrep(refFStream, "TOTAL STRETCH-BEND STRAIN ENERGY",
917 energyString);
918 if (!found) {
919 errorMsg += ref + corruptedMsg;
920 break;
921 }
922 std::stringstream refStretchBendStringStream(energyString);
923 refStretchBendStringStream >> skip >> skip >> skip >> skip >>
924 skip >> refEnergy;
925 error = (rdkStretchBendInstanceVec.size() !=
926 refStretchBendInstanceVec.size());
927 if (error) {
928 failed = true;
929 std::stringstream diff;
930 diff << "Stretch-bending: expected "
931 << refStretchBendInstanceVec.size()
932 << " interactions, found "
933 << rdkStretchBendInstanceVec.size() << "\n";
934 errorMsg += diff.str();
935 }
936 for (unsigned j = 0;
937 (!error) & (j < rdkStretchBendInstanceVec.size()); ++j) {
938 error =
939 ((rdkStretchBendInstanceVec[j]->iAtomType !=
940 refStretchBendInstanceVec[j]->iAtomType) ||
941 (rdkStretchBendInstanceVec[j]->jAtomType !=
942 refStretchBendInstanceVec[j]->jAtomType) ||
943 (rdkStretchBendInstanceVec[j]->kAtomType !=
944 refStretchBendInstanceVec[j]->kAtomType) ||
945 (rdkStretchBendInstanceVec[j]->ffType !=
946 refStretchBendInstanceVec[j]->ffType) ||
947 (fabs(rdkStretchBendInstanceVec[j]->kba -
948 refStretchBendInstanceVec[j]->kba) > FCON_TOLERANCE));
949 if (error) {
950 failed = true;
951 bool haveLine;
952 haveLine =
953 getLineByNum(rdkFStream, rdkPos,
954 rdkStretchBendInstanceVec[j]->idx, rdkLine);
955 if (haveLine) {
956 haveLine =
957 getLineByNum(refFStream, refPos,
958 refStretchBendInstanceVec[j]->idx, refLine);
959 }
960 std::stringstream diff;
961 diff << "Stretch-bending: found a difference\n";
962 if (haveLine) {
963 diff << "Expected:\n"
964 << rdkLine << "\nFound:\n"
965 << refLine << "\n";
966 }
967 errorMsg += diff.str();
968 }
969 delete rdkStretchBendInstanceVec[j];
970 delete refStretchBendInstanceVec[j];
971 }
972 error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
973 if (error && checkEnergy) {
974 failed = true;
975 std::stringstream diff;
976 diff << "Stretch-bending: energies differ\n"
977 "Expected "
978 << std::fixed << std::setprecision(4) << refEnergy
979 << ", found " << std::fixed << std::setprecision(4)
980 << rdkEnergy << "\n";
981 errorMsg += diff.str();
982 }
983 }
984
985 // --------------------------------------------------
986 // O U T - O F - P L A N E B E N D I N G
987 // --------------------------------------------------
988 if (refHaveOopBend) {
989 corruptedMsg = ": Out-of-plane bending: corrupted input data\n";
990 found =
991 fgrep(rdkFStream, "O U T - O F - P L A N E B E N D I N G");
992 if (!found) {
993 errorMsg += rdk + corruptedMsg;
994 break;
995 }
996 found = fgrep(rdkFStream, "----------------");
997 if (!found) {
998 errorMsg += rdk + corruptedMsg;
999 break;
1000 }
1001 found = false;
1002 rdkOopBendInstanceVec.clear();
1003 rdkPos = rdkFStream.tellg();
1004 n = 0;
1005 while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
1006 std::ifstream::failbit))) {
1007 found = (!(rdkLine.length()));
1008 if (found) {
1009 break;
1010 }
1011 std::stringstream rdkLineStream(rdkLine);
1012 auto *oopBendInstance = new OopBendInstance();
1013 oopBendInstance->idx = n;
1014 rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
1015 skip >> skip >> oopBendInstance->iAtomType >>
1016 oopBendInstance->jAtomType >> oopBendInstance->kAtomType >>
1017 oopBendInstance->lAtomType >> skip >> skip >>
1018 oopBendInstance->koop;
1019 fixOopBendInstance(oopBendInstance);
1020 rdkOopBendInstanceVec.push_back(oopBendInstance);
1021 ++n;
1022 }
1023 if (!found) {
1024 errorMsg += rdk + corruptedMsg;
1025 break;
1026 }
1027 std::sort(rdkOopBendInstanceVec.begin(),
1028 rdkOopBendInstanceVec.end(), sortOopBendInstanceVec);
1029 found = fgrep(rdkFStream, "TOTAL OUT-OF-PLANE BEND ENERGY",
1030 energyString);
1031 if (!found) {
1032 errorMsg += rdk + corruptedMsg;
1033 break;
1034 }
1035 std::stringstream rdkOopBendStringStream(energyString);
1036 rdkOopBendStringStream >> skip >> skip >> skip >> skip >> skip >>
1037 rdkEnergy;
1038 found =
1039 fgrep(refFStream, "O U T - O F - P L A N E B E N D I N G");
1040 if (!found) {
1041 errorMsg += ref + corruptedMsg;
1042 break;
1043 }
1044 found = fgrep(refFStream, "----------------");
1045 if (!found) {
1046 errorMsg += ref + corruptedMsg;
1047 break;
1048 }
1049 found = false;
1050 refOopBendInstanceVec.clear();
1051 refPos = refFStream.tellg();
1052 n = 0;
1053 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
1054 std::ifstream::failbit))) {
1055 found = (!(refLine.length()));
1056 if (found) {
1057 break;
1058 }
1059 std::stringstream refLineStream(refLine);
1060 auto *oopBendInstance = new OopBendInstance();
1061 oopBendInstance->idx = n;
1062 refLineStream >> skip >> skip >> skip >> skip >> skip >>
1063 oopBendInstance->iAtomType >> oopBendInstance->jAtomType >>
1064 oopBendInstance->kAtomType >> oopBendInstance->lAtomType >>
1065 skip >> skip >> oopBendInstance->koop;
1066 fixOopBendInstance(oopBendInstance);
1067 refOopBendInstanceVec.push_back(oopBendInstance);
1068 ++n;
1069 }
1070 if (!found) {
1071 errorMsg += ref + corruptedMsg;
1072 break;
1073 }
1074 std::sort(refOopBendInstanceVec.begin(),
1075 refOopBendInstanceVec.end(), sortOopBendInstanceVec);
1076 found = fgrep(refFStream, "TOTAL OUT-OF-PLANE STRAIN ENERGY",
1077 energyString);
1078 if (!found) {
1079 errorMsg += ref + corruptedMsg;
1080 break;
1081 }
1082 std::stringstream refOopBendStringStream(energyString);
1083 refOopBendStringStream >> skip >> skip >> skip >> skip >> skip >>
1084 refEnergy;
1085 error =
1086 (rdkOopBendInstanceVec.size() != refOopBendInstanceVec.size());
1087 if (error) {
1088 failed = true;
1089 std::stringstream diff;
1090 diff << "Out-of-plane bending: expected "
1091 << refOopBendInstanceVec.size() << " interactions, found "
1092 << rdkOopBendInstanceVec.size() << "\n";
1093 errorMsg += diff.str();
1094 }
1095 for (unsigned j = 0; (!error) & (j < rdkOopBendInstanceVec.size());
1096 ++j) {
1097 error = ((rdkOopBendInstanceVec[j]->iAtomType !=
1098 refOopBendInstanceVec[j]->iAtomType) ||
1099 (rdkOopBendInstanceVec[j]->jAtomType !=
1100 refOopBendInstanceVec[j]->jAtomType) ||
1101 (rdkOopBendInstanceVec[j]->kAtomType !=
1102 refOopBendInstanceVec[j]->kAtomType) ||
1103 (rdkOopBendInstanceVec[j]->lAtomType !=
1104 refOopBendInstanceVec[j]->lAtomType) ||
1105 (fabs(rdkOopBendInstanceVec[j]->koop -
1106 refOopBendInstanceVec[j]->koop) > FCON_TOLERANCE));
1107 if (error) {
1108 failed = true;
1109 bool haveLine;
1110 haveLine = getLineByNum(rdkFStream, rdkPos,
1111 rdkOopBendInstanceVec[j]->idx, rdkLine);
1112 if (haveLine) {
1113 haveLine =
1114 getLineByNum(refFStream, refPos,
1115 refOopBendInstanceVec[j]->idx, refLine);
1116 }
1117 std::stringstream diff;
1118 diff << "Out-of-plane bending: found a difference\n";
1119 if (haveLine) {
1120 diff << "Expected:\n"
1121 << rdkLine << "\nFound:\n"
1122 << refLine << "\n";
1123 }
1124 errorMsg += diff.str();
1125 }
1126 delete rdkOopBendInstanceVec[j];
1127 delete refOopBendInstanceVec[j];
1128 }
1129 error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
1130 if (error && checkEnergy) {
1131 failed = true;
1132 std::stringstream diff;
1133 diff << "Out-of-plane bending: energies differ\n"
1134 "Expected "
1135 << std::fixed << std::setprecision(4) << refEnergy
1136 << ", found " << std::fixed << std::setprecision(4)
1137 << rdkEnergy << "\n";
1138 errorMsg += diff.str();
1139 }
1140 }
1141
1142 // --------------------------------------------------
1143 // T O R S I O N A L
1144 // --------------------------------------------------
1145 if (refHaveTorsion) {
1146 corruptedMsg = ": Torsional: corrupted input data\n";
1147 found = fgrep(rdkFStream, "T O R S I O N A L");
1148 if (!found) {
1149 errorMsg += rdk + corruptedMsg;
1150 break;
1151 }
1152 found = fgrep(rdkFStream, "----------------");
1153 if (!found) {
1154 errorMsg += rdk + corruptedMsg;
1155 break;
1156 }
1157 found = false;
1158 rdkTorsionInstanceVec.clear();
1159 rdkPos = rdkFStream.tellg();
1160 n = 0;
1161 while ((!found) && (!(std::getline(rdkFStream, rdkLine).rdstate() &
1162 std::ifstream::failbit))) {
1163 found = (!(rdkLine.length()));
1164 if (found) {
1165 break;
1166 }
1167 std::stringstream rdkLineStream(rdkLine);
1168 auto *torsionInstance = new TorsionInstance();
1169 torsionInstance->idx = n;
1170 rdkLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
1171 skip >> skip >> torsionInstance->iAtomType >>
1172 torsionInstance->jAtomType >> torsionInstance->kAtomType >>
1173 torsionInstance->lAtomType >> torsionInstance->ffType >>
1174 skip >> skip >> torsionInstance->V1 >> torsionInstance->V2 >>
1175 torsionInstance->V3;
1176 fixTorsionInstance(torsionInstance);
1177 rdkTorsionInstanceVec.push_back(torsionInstance);
1178 ++n;
1179 }
1180 if (!found) {
1181 errorMsg += rdk + corruptedMsg;
1182 break;
1183 }
1184 std::sort(rdkTorsionInstanceVec.begin(),
1185 rdkTorsionInstanceVec.end(), sortTorsionInstanceVec);
1186 found = fgrep(rdkFStream, "TOTAL TORSIONAL ENERGY", energyString);
1187 if (!found) {
1188 errorMsg += rdk + corruptedMsg;
1189 break;
1190 }
1191 std::stringstream rdkTorsionStringStream(energyString);
1192 rdkTorsionStringStream >> skip >> skip >> skip >> skip >> rdkEnergy;
1193 found = fgrep(refFStream, "T O R S I O N A L");
1194 if (!found) {
1195 errorMsg += ref + corruptedMsg;
1196 break;
1197 }
1198 found = fgrep(refFStream, "----------------");
1199 if (!found) {
1200 errorMsg += ref + corruptedMsg;
1201 break;
1202 }
1203 found = false;
1204 refTorsionInstanceVec.clear();
1205 refPos = refFStream.tellg();
1206 n = 0;
1207 double refTorsionEnergy;
1208 while ((!found) && (!(std::getline(refFStream, refLine).rdstate() &
1209 std::ifstream::failbit))) {
1210 found = (!(refLine.length()));
1211 if (found) {
1212 break;
1213 }
1214 std::stringstream refLineStream(refLine);
1215 auto *torsionInstance = new TorsionInstance();
1216 torsionInstance->idx = n;
1217 refLineStream >> skip >> skip >> skip >> skip >> skip >> skip >>
1218 torsionInstance->iAtomType >> torsionInstance->jAtomType >>
1219 torsionInstance->kAtomType >> torsionInstance->lAtomType >>
1220 torsionInstance->ffType >> skip >> refTorsionEnergy >>
1221 torsionInstance->V1 >> torsionInstance->V2 >>
1222 torsionInstance->V3;
1223 if (!(ForceFields::MMFF::isDoubleZero(torsionInstance->V1) &&
1224 ForceFields::MMFF::isDoubleZero(torsionInstance->V2) &&
1225 ForceFields::MMFF::isDoubleZero(torsionInstance->V3) &&
1226 ForceFields::MMFF::isDoubleZero(refTorsionEnergy))) {
1227 fixTorsionInstance(torsionInstance);
1228 refTorsionInstanceVec.push_back(torsionInstance);
1229 } else {
1230 delete torsionInstance;
1231 }
1232 ++n;
1233 }
1234 if (!found) {
1235 errorMsg += ref + corruptedMsg;
1236 break;
1237 }
1238 std::sort(refTorsionInstanceVec.begin(),
1239 refTorsionInstanceVec.end(), sortTorsionInstanceVec);
1240 found =
1241 fgrep(refFStream, "TOTAL TORSION STRAIN ENERGY", energyString);
1242 if (!found) {
1243 errorMsg += ref + corruptedMsg;
1244 break;
1245 }
1246 std::stringstream refTorsionStringStream(energyString);
1247 refTorsionStringStream >> skip >> skip >> skip >> skip >> skip >>
1248 refEnergy;
1249 error =
1250 (rdkTorsionInstanceVec.size() != refTorsionInstanceVec.size());
1251 if (error) {
1252 failed = true;
1253 std::stringstream diff;
1254 diff << "Torsional: expected " << refTorsionInstanceVec.size()
1255 << " interactions, found " << rdkTorsionInstanceVec.size()
1256 << "\n";
1257 errorMsg += diff.str();
1258 }
1259 for (unsigned j = 0; (!error) & (j < rdkTorsionInstanceVec.size());
1260 ++j) {
1261 error = ((rdkTorsionInstanceVec[j]->iAtomType !=
1262 refTorsionInstanceVec[j]->iAtomType) ||
1263 (rdkTorsionInstanceVec[j]->jAtomType !=
1264 refTorsionInstanceVec[j]->jAtomType) ||
1265 (rdkTorsionInstanceVec[j]->kAtomType !=
1266 refTorsionInstanceVec[j]->kAtomType) ||
1267 (rdkTorsionInstanceVec[j]->lAtomType !=
1268 refTorsionInstanceVec[j]->lAtomType) ||
1269 (rdkTorsionInstanceVec[j]->ffType !=
1270 refTorsionInstanceVec[j]->ffType) ||
1271 (fabs(rdkTorsionInstanceVec[j]->V1 -
1272 refTorsionInstanceVec[j]->V1) > FCON_TOLERANCE) ||
1273 (fabs(rdkTorsionInstanceVec[j]->V2 -
1274 refTorsionInstanceVec[j]->V2) > FCON_TOLERANCE) ||
1275 (fabs(rdkTorsionInstanceVec[j]->V3 -
1276 refTorsionInstanceVec[j]->V3) > FCON_TOLERANCE));
1277 if (error) {
1278 failed = true;
1279 bool haveLine;
1280 haveLine = getLineByNum(rdkFStream, rdkPos,
1281 rdkTorsionInstanceVec[j]->idx, rdkLine);
1282 if (haveLine) {
1283 haveLine =
1284 getLineByNum(refFStream, refPos,
1285 refTorsionInstanceVec[j]->idx, refLine);
1286 }
1287 std::stringstream diff;
1288 diff << "Torsional: found a difference\n";
1289 if (haveLine) {
1290 diff << "Expected:\n"
1291 << rdkLine << "\nFound:\n"
1292 << refLine << "\n";
1293 }
1294 errorMsg += diff.str();
1295 }
1296 }
1297 for (auto &ti : rdkTorsionInstanceVec) {
1298 delete ti;
1299 }
1300 for (auto &ti : refTorsionInstanceVec) {
1301 delete ti;
1302 }
1303 error = (fabs(rdkEnergy - refEnergy) > ENERGY_TOLERANCE);
1304 if (error && checkEnergy) {
1305 failed = true;
1306 std::stringstream diff;
1307 diff << "Torsional: energies differ\n"
1308 "Expected "
1309 << std::fixed << std::setprecision(4) << refEnergy
1310 << ", found " << std::fixed << std::setprecision(4)
1311 << rdkEnergy << "\n";
1312 errorMsg += diff.str();
1313 }
1314 }
1315
1316 // --------------------------------------------------
1317 // V A N D E R W A A L S
1318 // E L E C T R O S T A T I C
1319 // --------------------------------------------------
1320 corruptedMsg = ": Van der Waals: corrupted input data\n";
1321 found = fgrep(rdkFStream, "V A N D E R W A A L S");
1322 if (!found) {
1323 errorMsg += rdk + corruptedMsg;
1324 break;
1325 }
1326 found = fgrep(rdkFStream, "----------------");
1327 if (!found) {
1328 errorMsg += rdk + corruptedMsg;
1329 break;
1330 }
1331 found = fgrep(rdkFStream, "TOTAL VAN DER WAALS ENERGY", energyString);
1332 if (!found) {
1333 errorMsg += rdk + corruptedMsg;
1334 break;
1335 }
1336 if (!refHaveVdW) {
1337 errorMsg += ref + corruptedMsg;
1338 break;
1339 }
1340 std::stringstream rdkVdWStringStream(energyString);
1341 rdkVdWStringStream >> skip >> skip >> skip >> skip >> skip >> skip >>
1342 rdkEnergy;
1343 corruptedMsg = ": Electrostatic: corrupted input data";
1344 found = fgrep(rdkFStream, "E L E C T R O S T A T I C");
1345 if (!found) {
1346 errorMsg += rdk + corruptedMsg;
1347 break;
1348 }
1349 found = fgrep(rdkFStream, "----------------");
1350 if (!found) {
1351 errorMsg += rdk + corruptedMsg;
1352 break;
1353 }
1354 found = fgrep(rdkFStream, "TOTAL ELECTROSTATIC ENERGY", energyString);
1355 if (!found) {
1356 errorMsg += rdk + corruptedMsg;
1357 break;
1358 }
1359 if (!refHaveEle) {
1360 errorMsg += ref + corruptedMsg;
1361 break;
1362 }
1363 std::stringstream rdkEleStringStream(energyString);
1364 rdkEleStringStream >> skip >> skip >> skip >> skip >> rdkEleEnergy;
1365 error = (fabs(rdkEnergy - refVdWEnergy) > ENERGY_TOLERANCE);
1366 if (error && checkEnergy) {
1367 failed = true;
1368 std::stringstream diff;
1369 diff << "Van der Waals: energies differ\n"
1370 "Expected "
1371 << std::fixed << std::setprecision(4) << refVdWEnergy
1372 << ", found " << std::fixed << std::setprecision(4)
1373 << rdkEnergy << "\n";
1374 errorMsg += diff.str();
1375 }
1376 error = (fabs(rdkEleEnergy - refEleEnergy) > ENERGY_TOLERANCE);
1377 if (error && checkEnergy) {
1378 failed = true;
1379 std::stringstream diff;
1380 diff << "Electrostatic: energies differ\n"
1381 "Expected "
1382 << std::fixed << std::setprecision(4) << refEleEnergy
1383 << ", found " << std::fixed << std::setprecision(4)
1384 << rdkEnergy << "\n";
1385 errorMsg += diff.str();
1386 }
1387 if (failed) {
1388 std::cerr << nameArray[ii] << "\n\n" << errorMsg << std::endl;
1389 ++n_failed;
1390 }
1391 }
1392 if (!found) {
1393 if (!failed) {
1394 std::cerr << nameArray[ii] << "\n\n";
1395 }
1396 std::cerr << errorMsg << std::endl;
1397 } else {
1398 unsigned int n_passed = nameArray.size() - n_failed;
1399 std::cerr << (n_failed ? "" : "All ") << n_passed << " test"
1400 << ((n_passed == 1) ? "" : "s") << " passed" << std::endl;
1401 if (n_failed) {
1402 std::cerr << n_failed << " test" << ((n_failed == 1) ? "" : "s")
1403 << " failed" << std::endl;
1404 }
1405 }
1406 if (n_failed) {
1407 testFailure = true;
1408 }
1409 rdkFStream.close();
1410 refFStream.close();
1411 firstTest = false;
1412 }
1413 }
1414 }
1415
1416 TEST_ASSERT(!testFailure);
1417 std::cerr << " done" << std::endl;
1418 return 0;
1419 }
1420
testMMFFAllConstraints()1421 void testMMFFAllConstraints() {
1422 std::cerr << "-------------------------------------" << std::endl;
1423 std::cerr << "Unit tests for all MMFF constraint terms." << std::endl;
1424
1425 std::string molBlock =
1426 "butane\n"
1427 " RDKit 3D\n"
1428 "butane\n"
1429 " 17 16 0 0 0 0 0 0 0 0999 V2000\n"
1430 " 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1431 " 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1432 " 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1433 " 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1434 " 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1435 " 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1436 " 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1437 " 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1438 " 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1439 " 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1440 " 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1441 " 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1442 " 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1443 " 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1444 " -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1445 " -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1446 " -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1447 " 1 2 1 0 0 0 0\n"
1448 " 1 15 1 0 0 0 0\n"
1449 " 1 16 1 0 0 0 0\n"
1450 " 1 17 1 0 0 0 0\n"
1451 " 2 3 1 0 0 0 0\n"
1452 " 2 4 1 0 0 0 0\n"
1453 " 2 14 1 0 0 0 0\n"
1454 " 4 5 1 0 0 0 0\n"
1455 " 4 6 1 0 0 0 0\n"
1456 " 4 7 1 0 0 0 0\n"
1457 " 7 8 1 0 0 0 0\n"
1458 " 7 9 1 0 0 0 0\n"
1459 " 7 13 1 0 0 0 0\n"
1460 " 9 10 1 0 0 0 0\n"
1461 " 9 11 1 0 0 0 0\n"
1462 " 9 12 1 0 0 0 0\n"
1463 "M END\n";
1464 RDKit::RWMol *mol;
1465 ForceFields::ForceField *field;
1466 MMFF::MMFFMolProperties *mmffMolProperties;
1467
1468 // distance constraints
1469 ForceFields::MMFF::DistanceConstraintContrib *dc;
1470 mol = RDKit::MolBlockToMol(molBlock, true, false);
1471 TEST_ASSERT(mol);
1472 MolTransforms::setBondLength(mol->getConformer(), 1, 3, 2.0);
1473 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1474 TEST_ASSERT(mmffMolProperties);
1475 TEST_ASSERT(mmffMolProperties->isValid());
1476 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1477 TEST_ASSERT(field);
1478 field->initialize();
1479 dc = new ForceFields::MMFF::DistanceConstraintContrib(field, 1, 3, 2.0, 2.0,
1480 1.0e5);
1481 field->contribs().push_back(ForceFields::ContribPtr(dc));
1482 field->minimize();
1483 TEST_ASSERT(RDKit::feq(
1484 MolTransforms::getBondLength(mol->getConformer(), 1, 3), 2.0, 0.1));
1485 delete field;
1486 delete mmffMolProperties;
1487 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1488 TEST_ASSERT(mmffMolProperties);
1489 TEST_ASSERT(mmffMolProperties->isValid());
1490 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1491 field->initialize();
1492 dc = new ForceFields::MMFF::DistanceConstraintContrib(field, 1, 3, true, -0.2,
1493 0.2, 1.0e5);
1494 field->contribs().push_back(ForceFields::ContribPtr(dc));
1495 field->minimize();
1496 TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.79);
1497 delete field;
1498 delete mmffMolProperties;
1499 delete mol;
1500
1501 // angle constraints
1502 ForceFields::MMFF::AngleConstraintContrib *ac;
1503 mol = RDKit::MolBlockToMol(molBlock, true, false);
1504 TEST_ASSERT(mol);
1505 MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 90.0);
1506 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1507 TEST_ASSERT(mmffMolProperties);
1508 TEST_ASSERT(mmffMolProperties->isValid());
1509 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1510 TEST_ASSERT(field);
1511 field->initialize();
1512 ac = new ForceFields::MMFF::AngleConstraintContrib(field, 1, 3, 6, 90.0, 90.0,
1513 100.0);
1514 field->contribs().push_back(ForceFields::ContribPtr(ac));
1515 field->minimize();
1516 TEST_ASSERT(RDKit::feq(
1517 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 90.0, 0.5));
1518 delete field;
1519 delete mmffMolProperties;
1520 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1521 TEST_ASSERT(mmffMolProperties);
1522 TEST_ASSERT(mmffMolProperties->isValid());
1523 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1524 field->initialize();
1525 ac = new ForceFields::MMFF::AngleConstraintContrib(field, 1, 3, 6, true,
1526 -10.0, 10.0, 100.0);
1527 field->contribs().push_back(ForceFields::ContribPtr(ac));
1528 field->minimize();
1529 TEST_ASSERT(RDKit::feq(
1530 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 100.0, 0.5));
1531 delete field;
1532 delete mmffMolProperties;
1533 MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 0.0);
1534 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1535 TEST_ASSERT(mmffMolProperties);
1536 TEST_ASSERT(mmffMolProperties->isValid());
1537 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1538 field->initialize();
1539 ac = new ForceFields::MMFF::AngleConstraintContrib(field, 1, 3, 6, false,
1540 -10.0, 10.0, 100.0);
1541 field->contribs().push_back(ForceFields::ContribPtr(ac));
1542 field->minimize();
1543 TEST_ASSERT(RDKit::feq(
1544 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 10.0, 0.5));
1545 delete field;
1546 delete mmffMolProperties;
1547 delete mol;
1548
1549 // torsion constraints
1550 ForceFields::MMFF::TorsionConstraintContrib *tc;
1551 mol = RDKit::MolBlockToMol(molBlock, true, false);
1552 TEST_ASSERT(mol);
1553 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, 15.0);
1554 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1555 TEST_ASSERT(mmffMolProperties);
1556 TEST_ASSERT(mmffMolProperties->isValid());
1557 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1558 TEST_ASSERT(field);
1559 field->initialize();
1560 tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, 10.0,
1561 10.0, 100.0);
1562 field->contribs().push_back(ForceFields::ContribPtr(tc));
1563 field->minimize();
1564 std::cerr << MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8)
1565 << std::endl;
1566 TEST_ASSERT(
1567 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1568 10.0, 0.5));
1569 delete field;
1570 delete mmffMolProperties;
1571 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -30.0);
1572 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1573 TEST_ASSERT(mmffMolProperties);
1574 TEST_ASSERT(mmffMolProperties->isValid());
1575 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1576 field->initialize();
1577 tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, true,
1578 -1.0, 1.0, 100.0);
1579 field->contribs().push_back(ForceFields::ContribPtr(tc));
1580 field->minimize();
1581 TEST_ASSERT(
1582 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1583 -30.0, 1.5));
1584 delete field;
1585 delete mmffMolProperties;
1586 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -10.0);
1587 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1588 TEST_ASSERT(mmffMolProperties);
1589 TEST_ASSERT(mmffMolProperties->isValid());
1590 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1591 field->initialize();
1592 tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 1, 3, 6, 8, false,
1593 -10.0, -8.0, 100.0);
1594 field->contribs().push_back(ForceFields::ContribPtr(tc));
1595 field->minimize(500);
1596 TEST_ASSERT(
1597 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1598 -9.0, 2.0));
1599 delete field;
1600 delete mmffMolProperties;
1601 delete mol;
1602
1603 // position constraints
1604 ForceFields::MMFF::PositionConstraintContrib *pc;
1605 mol = RDKit::MolBlockToMol(molBlock, true, false);
1606 TEST_ASSERT(mol);
1607 mmffMolProperties = new MMFF::MMFFMolProperties(*mol);
1608 TEST_ASSERT(mmffMolProperties);
1609 TEST_ASSERT(mmffMolProperties->isValid());
1610 field = RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1611 TEST_ASSERT(field);
1612 field->initialize();
1613 RDGeom::Point3D p = mol->getConformer().getAtomPos(1);
1614 pc = new ForceFields::MMFF::PositionConstraintContrib(field, 1, 0.3, 1.0e5);
1615 field->contribs().push_back(ForceFields::ContribPtr(pc));
1616 field->minimize();
1617 RDGeom::Point3D q = mol->getConformer().getAtomPos(1);
1618 TEST_ASSERT((p - q).length() < 0.3);
1619 delete field;
1620 delete mmffMolProperties;
1621 delete mol;
1622
1623 // fixed atoms
1624 mol = RDKit::MolBlockToMol(molBlock, true, false);
1625 TEST_ASSERT(mol);
1626 field = RDKit::MMFF::constructForceField(*mol);
1627 TEST_ASSERT(field);
1628 field->initialize();
1629 RDGeom::Point3D fp = mol->getConformer().getAtomPos(1);
1630 field->fixedPoints().push_back(1);
1631 field->minimize();
1632 RDGeom::Point3D fq = mol->getConformer().getAtomPos(1);
1633 TEST_ASSERT((fp - fq).length() < 0.01);
1634 delete field;
1635 delete mol;
1636
1637 std::cerr << " done" << std::endl;
1638 }
1639
testMMFFTorsionConstraints()1640 void testMMFFTorsionConstraints() {
1641 std::cerr << "-------------------------------------" << std::endl;
1642 std::cerr << "Unit test for MMFF TorsionConstraint close to 0 degrees."
1643 << std::endl;
1644
1645 std::string molBlock = R"(
1646 RDKit 3D
1647
1648 11 10 0 0 0 0 0 0 0 0999 V2000
1649 -1.6091 -0.3348 -0.0457 N 0 0 0 0 0 0 0 0 0 0 0 0
1650 -2.3935 0.2978 -0.1061 H 0 0 0 0 0 0 0 0 0 0 0 0
1651 -1.6708 -0.8686 0.8103 H 0 0 0 0 0 0 0 0 0 0 0 0
1652 -0.3560 0.4358 -0.0518 C 0 0 0 0 0 0 0 0 0 0 0 0
1653 -0.3583 1.1207 0.7968 H 0 0 0 0 0 0 0 0 0 0 0 0
1654 -0.2978 1.0122 -0.9759 H 0 0 0 0 0 0 0 0 0 0 0 0
1655 0.8573 -0.4893 0.0451 C 0 0 0 0 0 0 0 0 0 0 0 0
1656 0.8947 -1.1355 -0.8327 H 0 0 0 0 0 0 0 0 0 0 0 0
1657 0.7774 -1.1052 0.9415 H 0 0 0 0 0 0 0 0 0 0 0 0
1658 2.0396 0.2763 0.1127 O 0 0 0 0 0 0 0 0 0 0 0 0
1659 2.1165 0.7905 -0.6941 H 0 0 0 0 0 0 0 0 0 0 0 0
1660 1 2 1 0 0 0 0
1661 1 3 1 0 0 0 0
1662 1 4 1 0 0 0 0
1663 4 5 1 0 0 0 0
1664 4 6 1 0 0 0 0
1665 4 7 1 0 0 0 0
1666 7 8 1 0 0 0 0
1667 7 9 1 0 0 0 0
1668 7 10 1 0 0 0 0
1669 10 11 1 0 0 0 0
1670 M END
1671 )";
1672 RDKit::RWMol *mol;
1673
1674 // torsion constraints
1675 mol = RDKit::MolBlockToMol(molBlock, true, false);
1676 TEST_ASSERT(mol);
1677 std::vector<double> e;
1678 for (double d : {0.0001, 0.0}) {
1679 MolTransforms::setDihedralDeg(mol->getConformer(), 0, 3, 6, 9, d);
1680 MMFF::MMFFMolProperties mmffMolProperties(*mol);
1681 TEST_ASSERT(mmffMolProperties.isValid());
1682 ForceFields::ForceField *field =
1683 RDKit::MMFF::constructForceField(*mol, &mmffMolProperties);
1684 TEST_ASSERT(field);
1685 field->initialize();
1686 auto *tc = new ForceFields::MMFF::TorsionConstraintContrib(field, 0, 3, 6,
1687 9, d, d, 1.0e6);
1688 field->contribs().push_back(ForceFields::ContribPtr(tc));
1689 field->minimize();
1690 e.push_back(field->calcEnergy());
1691 TEST_ASSERT(RDKit::feq(
1692 MolTransforms::getDihedralDeg(mol->getConformer(), 0, 3, 6, 9), d,
1693 0.5));
1694 delete field;
1695 }
1696 TEST_ASSERT(RDKit::feq(e[0], e[1], 0.5));
1697 delete mol;
1698
1699 std::cerr << " done" << std::endl;
1700 }
1701
testMMFFCopy()1702 void testMMFFCopy() {
1703 std::cerr << "-------------------------------------" << std::endl;
1704 std::cerr << "Unit tests for copying MMFF ForceFields." << std::endl;
1705
1706 std::string molBlock =
1707 "butane\n"
1708 " RDKit 3D\n"
1709 "butane\n"
1710 " 17 16 0 0 0 0 0 0 0 0999 V2000\n"
1711 " 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1712 " 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1713 " 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1714 " 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1715 " 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1716 " 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1717 " 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1718 " 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1719 " 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1720 " 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1721 " 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1722 " 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1723 " 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1724 " 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1725 " -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1726 " -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1727 " -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1728 " 1 2 1 0 0 0 0\n"
1729 " 1 15 1 0 0 0 0\n"
1730 " 1 16 1 0 0 0 0\n"
1731 " 1 17 1 0 0 0 0\n"
1732 " 2 3 1 0 0 0 0\n"
1733 " 2 4 1 0 0 0 0\n"
1734 " 2 14 1 0 0 0 0\n"
1735 " 4 5 1 0 0 0 0\n"
1736 " 4 6 1 0 0 0 0\n"
1737 " 4 7 1 0 0 0 0\n"
1738 " 7 8 1 0 0 0 0\n"
1739 " 7 9 1 0 0 0 0\n"
1740 " 7 13 1 0 0 0 0\n"
1741 " 9 10 1 0 0 0 0\n"
1742 " 9 11 1 0 0 0 0\n"
1743 " 9 12 1 0 0 0 0\n"
1744 "M END\n";
1745
1746 {
1747 RDKit::RWMol *mol = RDKit::MolBlockToMol(molBlock, true, false);
1748 TEST_ASSERT(mol);
1749 auto *cmol = new RDKit::RWMol(*mol);
1750 TEST_ASSERT(cmol);
1751
1752 MMFF::MMFFMolProperties *mmffMolProperties =
1753 new MMFF::MMFFMolProperties(*mol);
1754 TEST_ASSERT(mmffMolProperties);
1755 TEST_ASSERT(mmffMolProperties->isValid());
1756 ForceFields::ForceField *field =
1757 RDKit::MMFF::constructForceField(*mol, mmffMolProperties);
1758 TEST_ASSERT(field);
1759 field->initialize();
1760 auto *dc = new ForceFields::MMFF::DistanceConstraintContrib(
1761 field, 1, 3, 2.0, 2.0, 1.0e5);
1762 field->contribs().push_back(ForceFields::ContribPtr(dc));
1763 field->minimize();
1764 TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.99);
1765
1766 auto *cfield = new ForceFields::ForceField(*field);
1767 cfield->positions().clear();
1768
1769 for (unsigned int i = 0; i < cmol->getNumAtoms(); i++) {
1770 cfield->positions().push_back(&cmol->getConformer().getAtomPos(i));
1771 }
1772 cfield->initialize();
1773 cfield->minimize();
1774 TEST_ASSERT(MolTransforms::getBondLength(cmol->getConformer(), 1, 3) >
1775 1.99);
1776 TEST_ASSERT(RDKit::feq(field->calcEnergy(), cfield->calcEnergy()));
1777
1778 const RDKit::Conformer &conf = mol->getConformer();
1779 const RDKit::Conformer &cconf = cmol->getConformer();
1780 for (unsigned int i = 0; i < mol->getNumAtoms(); i++) {
1781 RDGeom::Point3D p = conf.getAtomPos(i);
1782 RDGeom::Point3D cp = cconf.getAtomPos(i);
1783 TEST_ASSERT(RDKit::feq(p.x, cp.x));
1784 TEST_ASSERT(RDKit::feq(p.y, cp.y));
1785 TEST_ASSERT(RDKit::feq(p.z, cp.z));
1786 }
1787
1788 delete field;
1789 delete cfield;
1790 delete mmffMolProperties;
1791 delete mol;
1792 delete cmol;
1793 }
1794 std::cerr << " done" << std::endl;
1795 }
1796
testMMFFButaneScan()1797 void testMMFFButaneScan() {
1798 std::cerr << "-------------------------------------" << std::endl;
1799 std::cerr << "Unit test for MMFF butane scan." << std::endl;
1800
1801 auto molblock = R"(
1802 RDKit 3D
1803
1804 14 13 0 0 0 0 0 0 0 0999 V2000
1805 -1.5575 0.5684 -0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
1806 -0.6987 -0.6064 0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
1807 0.6987 -0.6064 -0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
1808 1.5575 0.5684 0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
1809 -1.6308 0.6129 -1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
1810 -2.5700 0.4662 0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
1811 -1.1524 1.5190 0.2272 H 0 0 0 0 0 0 0 0 0 0 0 0
1812 -1.2059 -1.5359 0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
1813 -0.6215 -0.6099 1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
1814 0.6215 -0.6099 -1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
1815 1.2059 -1.5359 -0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
1816 1.6308 0.6129 1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
1817 1.1524 1.5190 -0.2271 H 0 0 0 0 0 0 0 0 0 0 0 0
1818 2.5700 0.4662 -0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
1819 1 2 1 0
1820 2 3 1 0
1821 3 4 1 0
1822 1 5 1 0
1823 1 6 1 0
1824 1 7 1 0
1825 2 8 1 0
1826 2 9 1 0
1827 3 10 1 0
1828 3 11 1 0
1829 4 12 1 0
1830 4 13 1 0
1831 4 14 1 0
1832 M END
1833 )";
1834
1835 // torsion constraints
1836 std::unique_ptr<ROMol> mol(MolBlockToMol(molblock, true, false));
1837 TEST_ASSERT(mol);
1838 std::vector<std::pair<double, double>> diheEnergyVect;
1839 std::vector<std::pair<double, double>> expectedDiheEnergyVect{
1840 {-179.999, -5.07597}, {-175.1, -5.01299}, {-170.1, -4.82267},
1841 {-165.1, -4.51646}, {-160.1, -4.11308}, {-155.101, -3.63773},
1842 {-150.101, -3.12109}, {-145.101, -2.59801}, {-140.1, -2.10575},
1843 {-135.1, -1.68158}, {-130.1, -1.35952}, {-125.1, -1.16635},
1844 {-119.9, -1.11879}, {-114.9, -1.22194}, {-109.9, -1.45802},
1845 {-104.9, -1.80199}, {-99.8996, -2.22002}, {-94.8995, -2.67302},
1846 {-89.8996, -3.12055}, {-84.8996, -3.5255}, {-79.8997, -3.85916},
1847 {-74.8998, -4.10358}, {-69.8999, -4.24974}, {-65.1, -4.29368},
1848 {-60.1001, -4.23737}, {-55.1002, -4.0775}, {-50.1003, -3.81811},
1849 {-45.1004, -3.46765}, {-40.1005, -3.0395}, {-35.1005, -2.55211},
1850 {-30.1005, -2.02873}, {-25.1005, -1.49709}, {-20.1005, -0.988743},
1851 {-15.1004, -0.538157}, {-10.1003, -0.180672}, {-5.10016, 0.051218},
1852 {-0.100003, 0.13365}, {5.10016, 0.051218}, {10.1003, -0.180672},
1853 {15.1004, -0.538157}, {20.1005, -0.988743}, {25.1005, -1.49709},
1854 {30.1005, -2.02873}, {35.1005, -2.55211}, {40.1005, -3.0395},
1855 {45.1004, -3.46765}, {50.1003, -3.81811}, {55.1002, -4.0775},
1856 {60.1001, -4.23737}, {65.1, -4.29368}, {69.8999, -4.24974},
1857 {74.8998, -4.10358}, {79.8997, -3.85916}, {84.8996, -3.5255},
1858 {89.8996, -3.12055}, {94.8995, -2.67302}, {99.8996, -2.22002},
1859 {104.9, -1.80199}, {109.9, -1.45802}, {114.9, -1.22194},
1860 {119.9, -1.11879}, {125.1, -1.16635}, {130.1, -1.35952},
1861 {135.1, -1.68158}, {140.1, -2.10575}, {145.101, -2.59801},
1862 {150.101, -3.12109}, {155.101, -3.63773}, {160.1, -4.11308},
1863 {165.1, -4.51646}, {170.1, -4.82267}, {175.1, -5.01299}};
1864 MMFF::MMFFMolProperties mmffMolProperties(*mol);
1865 TEST_ASSERT(mmffMolProperties.isValid());
1866 std::unique_ptr<ForceFields::ForceField> globalFF(
1867 RDKit::MMFF::constructForceField(*mol, &mmffMolProperties));
1868 TEST_ASSERT(globalFF);
1869 globalFF->initialize();
1870 std::unique_ptr<ROMol> mol2(new ROMol(*mol));
1871 const int start = -180;
1872 const int end = 180;
1873 const int incr = 5;
1874 MolTransforms::setDihedralDeg(mol2->getConformer(), 0, 1, 2, 3, start);
1875 for (int diheIn = start; diheIn < end; diheIn += incr) {
1876 std::unique_ptr<ForceFields::ForceField> localFF(
1877 RDKit::MMFF::constructForceField(*mol2, &mmffMolProperties));
1878 TEST_ASSERT(localFF);
1879 localFF->initialize();
1880 auto *tc = new ForceFields::MMFF::TorsionConstraintContrib(
1881 localFF.get(), 0, 1, 2, 3, static_cast<double>(diheIn) - 0.1,
1882 static_cast<double>(diheIn) + 0.1, 100.0);
1883 localFF->contribs().push_back(ForceFields::ContribPtr(tc));
1884 TEST_ASSERT(localFF->minimize(10000) == 0);
1885 std::vector<double> pos(3 * localFF->numPoints());
1886 TEST_ASSERT(localFF->numPoints() == globalFF->numPoints());
1887 auto posns = localFF->positions();
1888 for (unsigned int i = 0; i < localFF->numPoints(); ++i) {
1889 for (unsigned int j = 0; j < 3; ++j) {
1890 pos[i * 3 + j] = (*static_cast<RDGeom::Point3D *>(posns.at(i)))[j];
1891 }
1892 }
1893 auto diheOut =
1894 MolTransforms::getDihedralDeg(mol2->getConformer(), 0, 1, 2, 3);
1895 auto energy = globalFF->calcEnergy(pos.data());
1896 diheEnergyVect.emplace_back(diheOut, energy);
1897 }
1898 for (size_t i = 0; i < diheEnergyVect.size(); ++i) {
1899 TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).first,
1900 expectedDiheEnergyVect.at(i).first, 1.0));
1901 TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).second,
1902 expectedDiheEnergyVect.at(i).second, 0.2));
1903 }
1904
1905 std::cerr << " done" << std::endl;
1906 }
1907
main(int argc,char * argv[])1908 int main(int argc, char *argv[]) {
1909 if (!mmffValidationSuite(argc, argv)) {
1910 testMMFFAllConstraints();
1911 testMMFFTorsionConstraints();
1912 testMMFFCopy();
1913 testMMFFButaneScan();
1914 }
1915
1916 return 0;
1917 }
1918