1 //
2 // Copyright (C) 2004-2018 Greg Landrum and Rational Discovery LLC
3 //
4 // @@ All Rights Reserved @@
5 // This file is part of the RDKit.
6 // The contents are covered by the terms of the BSD license
7 // which is included in the file license.txt, found at the root
8 // of the RDKit source tree.
9 //
10
11 #include <RDGeneral/test.h>
12 #include <iostream>
13 #include <RDGeneral/Invariant.h>
14 #include <RDGeneral/RDLog.h>
15 #include <RDGeneral/utils.h>
16
17 #include <GraphMol/RDKitBase.h>
18 #include <GraphMol/SmilesParse/SmilesParse.h>
19 #include <GraphMol/SmilesParse/SmilesWrite.h>
20 #include <GraphMol/FileParsers/FileParsers.h>
21 #include <GraphMol/FileParsers/MolSupplier.h>
22 #include <GraphMol/FileParsers/MolWriters.h>
23
24 #include <GraphMol/ForceFieldHelpers/FFConvenience.h>
25 #include <GraphMol/ForceFieldHelpers/UFF/AtomTyper.h>
26 #include <GraphMol/ForceFieldHelpers/UFF/Builder.h>
27 #include <GraphMol/ForceFieldHelpers/UFF/UFF.h>
28 #include <ForceField/ForceField.h>
29 #include <GraphMol/DistGeomHelpers/Embedder.h>
30
31 using namespace RDKit;
32 #if 1
testUFFTyper1()33 void testUFFTyper1() {
34 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
35 BOOST_LOG(rdErrorLog) << " Test UFF atom labels." << std::endl;
36
37 ROMol *mol;
38 std::string key;
39
40 mol = SmilesToMol("[SiH3]CC(=O)NC");
41 TEST_ASSERT(mol);
42
43 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
44 TEST_ASSERT(key == "Si3");
45 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
46 TEST_ASSERT(key == "C_3");
47 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
48 TEST_ASSERT(key == "C_R");
49 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(3));
50 TEST_ASSERT(key == "O_R");
51 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(4));
52 TEST_ASSERT(key == "N_R");
53
54 delete mol;
55 mol = SmilesToMol("CC(=O)C");
56 TEST_ASSERT(mol);
57
58 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
59 TEST_ASSERT(key == "C_3");
60 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
61 TEST_ASSERT(key == "C_2");
62 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
63 TEST_ASSERT(key == "O_2");
64 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(3));
65 TEST_ASSERT(key == "C_3");
66
67 delete mol;
68 mol = SmilesToMol("C(=O)S");
69 TEST_ASSERT(mol);
70
71 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
72 TEST_ASSERT(key == "C_2");
73 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
74 TEST_ASSERT(key == "O_2");
75 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
76 TEST_ASSERT(key == "S_3+2");
77
78 delete mol;
79 mol = SmilesToMol("SCS(=O)S(=O)(=O)O");
80 TEST_ASSERT(mol);
81
82 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
83 TEST_ASSERT(key == "S_3+2");
84 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
85 TEST_ASSERT(key == "C_3");
86 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
87 TEST_ASSERT(key == "S_3+4");
88 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(4));
89 TEST_ASSERT(key == "S_3+6");
90
91 delete mol;
92 mol = SmilesToMol("PCP(O)CP(=O)(=O)");
93 TEST_ASSERT(mol);
94
95 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
96 TEST_ASSERT(key == "P_3+3");
97 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
98 TEST_ASSERT(key == "C_3");
99 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
100 TEST_ASSERT(key == "P_3+3");
101 // FIX: I am not 100% convinced that this should be SP2
102 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(5));
103 TEST_ASSERT(key == "P_3+5");
104
105 delete mol;
106 mol = SmilesToMol("C(F)(Cl)(Br)I");
107 TEST_ASSERT(mol);
108
109 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
110 TEST_ASSERT(key == "C_3");
111 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
112 TEST_ASSERT(key == "F_");
113 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
114 TEST_ASSERT(key == "Cl");
115 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(3));
116 TEST_ASSERT(key == "Br");
117 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(4));
118 TEST_ASSERT(key == "I_");
119
120 delete mol;
121 mol = SmilesToMol("[Li].[Na].[K].[Rb].[Cs]");
122 TEST_ASSERT(mol);
123
124 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
125 TEST_ASSERT(key == "Li");
126 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
127 TEST_ASSERT(key == "Na");
128 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
129 TEST_ASSERT(key == "K_");
130 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(3));
131 TEST_ASSERT(key == "Rb");
132 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(4));
133 TEST_ASSERT(key == "Cs");
134 delete mol;
135
136 BOOST_LOG(rdErrorLog) << " done" << std::endl;
137 }
138
testUFFTyper2()139 void testUFFTyper2() {
140 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
141 BOOST_LOG(rdErrorLog) << " Test UFF atom typer." << std::endl;
142
143 ROMol *mol, *mol2;
144 std::string key;
145
146 mol = SmilesToMol("[SiH3]CC(=O)NC");
147 TEST_ASSERT(mol);
148
149 UFF::AtomicParamVect types;
150 bool foundAll;
151 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
152 TEST_ASSERT(foundAll);
153 TEST_ASSERT(types.size() == mol->getNumAtoms());
154 for (UFF::AtomicParamVect::const_iterator it = types.begin();
155 it != types.end(); it++) {
156 TEST_ASSERT((*it));
157 }
158 mol2 = MolOps::addHs(*mol);
159 delete mol;
160 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol2);
161 TEST_ASSERT(foundAll);
162 TEST_ASSERT(types.size() == mol2->getNumAtoms());
163 for (UFF::AtomicParamVect::const_iterator it = types.begin();
164 it != types.end(); it++) {
165 TEST_ASSERT((*it));
166 }
167 delete mol2;
168
169 // connected with sf.net bug 2094445
170 mol = SmilesToMol("[SiH2]=C");
171 TEST_ASSERT(mol);
172 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
173 TEST_ASSERT(key == "Si3");
174 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
175 TEST_ASSERT(key == "C_2");
176 delete mol;
177
178 mol = SmilesToMol("[AlH]=C");
179 TEST_ASSERT(mol);
180 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
181 TEST_ASSERT(key == "Al3");
182 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
183 TEST_ASSERT(key == "C_2");
184 delete mol;
185
186 mol = SmilesToMol("[Mg]=C");
187 TEST_ASSERT(mol);
188 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
189 TEST_ASSERT(key == "Mg3+2");
190 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
191 TEST_ASSERT(key == "C_2");
192 delete mol;
193
194 mol = SmilesToMol("[SiH3][Si]([SiH3])=C");
195 TEST_ASSERT(mol);
196 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(0));
197 TEST_ASSERT(key == "Si3");
198 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(1));
199 TEST_ASSERT(key == "Si3");
200 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(2));
201 TEST_ASSERT(key == "Si3");
202 key = UFF::Tools::getAtomLabel(mol->getAtomWithIdx(3));
203 TEST_ASSERT(key == "C_2");
204 delete mol;
205
206 BOOST_LOG(rdErrorLog) << " done" << std::endl;
207 }
208
testUFFBuilder1()209 void testUFFBuilder1() {
210 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
211 BOOST_LOG(rdErrorLog) << " Testing UFF builder tools." << std::endl;
212
213 ROMol *mol, *mol2;
214 std::string key;
215
216 UFF::AtomicParamVect types;
217 bool foundAll;
218 ForceFields::ForceField *field;
219 boost::shared_array<std::uint8_t> nbrMat;
220
221 mol = SmilesToMol("CC(O)C");
222 auto *conf = new Conformer(mol->getNumAtoms());
223 int cid = static_cast<int>(mol->addConformer(conf, true));
224 TEST_ASSERT(mol);
225 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
226 TEST_ASSERT(foundAll);
227 TEST_ASSERT(types.size() == mol->getNumAtoms());
228 field = new ForceFields::ForceField();
229 // add the atomic positions:
230 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
231 field->positions().push_back(&((conf->getAtomPos(i))));
232 }
233
234 UFF::Tools::addBonds(*mol, types, field);
235
236 TEST_ASSERT(field->contribs().size() == 3);
237
238 nbrMat = UFF::Tools::buildNeighborMatrix(*mol);
239 // the neighbor matrix is an upper triangular matrix
240 // position indices are as follows:
241 // 0 1 2 3
242 // 4 5 6
243 // 7 8
244 // 9
245 TEST_ASSERT(UFF::Tools::twoBitCellPos(mol->getNumAtoms(), 1, 1) == 4);
246 TEST_ASSERT(UFF::Tools::twoBitCellPos(mol->getNumAtoms(), 2, 1) == 5);
247 TEST_ASSERT(UFF::Tools::twoBitCellPos(mol->getNumAtoms(), 1, 2) == 5);
248 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 0) == UFF::Tools::RELATION_1_X);
249 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 1) == UFF::Tools::RELATION_1_2);
250 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 2) == UFF::Tools::RELATION_1_3);
251 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 3) == UFF::Tools::RELATION_1_3);
252
253 UFF::Tools::addAngles(*mol, types, field);
254 TEST_ASSERT(field->contribs().size() == 6);
255
256 // there are no non-bonded terms here:
257 UFF::Tools::addNonbonded(*mol, cid, types, field, nbrMat);
258 TEST_ASSERT(field->contribs().size() == 6);
259
260 // and no torsions either, until we add Hs:
261 UFF::Tools::addTorsions(*mol, types, field);
262 TEST_ASSERT(field->contribs().size() == 6);
263
264 delete mol;
265 delete field;
266 mol = SmilesToMol("CCOC");
267 auto *conf2 = new Conformer(mol->getNumAtoms());
268 cid = static_cast<int>(mol->addConformer(conf2, true));
269 TEST_ASSERT(mol);
270 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
271 TEST_ASSERT(foundAll);
272 TEST_ASSERT(types.size() == mol->getNumAtoms());
273 field = new ForceFields::ForceField();
274 // add the atomic positions:
275 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
276 field->positions().push_back(&(conf2->getAtomPos(i)));
277 }
278
279 UFF::Tools::addBonds(*mol, types, field);
280
281 TEST_ASSERT(field->contribs().size() == 3);
282
283 nbrMat = UFF::Tools::buildNeighborMatrix(*mol);
284 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 0) == UFF::Tools::RELATION_1_X);
285 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 1) == UFF::Tools::RELATION_1_2);
286 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 2) == UFF::Tools::RELATION_1_3);
287 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 3) == UFF::Tools::RELATION_1_X);
288
289 UFF::Tools::addAngles(*mol, types, field);
290 TEST_ASSERT(field->contribs().size() == 5);
291 UFF::Tools::addNonbonded(*mol, cid, types, field, nbrMat);
292 TEST_ASSERT(field->contribs().size() == 6);
293 UFF::Tools::addTorsions(*mol, types, field);
294 TEST_ASSERT(field->contribs().size() == 7);
295
296 delete mol;
297 delete field;
298 mol = SmilesToMol("CO");
299 auto *conf3 = new Conformer(mol->getNumAtoms());
300 cid = static_cast<int>(mol->addConformer(conf3, true));
301 TEST_ASSERT(mol);
302 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
303 TEST_ASSERT(foundAll);
304 TEST_ASSERT(types.size() == mol->getNumAtoms());
305
306 field = new ForceFields::ForceField();
307 // add the atomic positions:
308 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
309 field->positions().push_back(&(conf3->getAtomPos(i)));
310 }
311
312 UFF::Tools::addBonds(*mol, types, field);
313
314 TEST_ASSERT(field->contribs().size() == 1);
315
316 nbrMat = UFF::Tools::buildNeighborMatrix(*mol);
317 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 0) == UFF::Tools::RELATION_1_X);
318 TEST_ASSERT(UFF::Tools::getTwoBitCell(nbrMat, 1) == UFF::Tools::RELATION_1_2);
319
320 UFF::Tools::addAngles(*mol, types, field);
321 TEST_ASSERT(field->contribs().size() == 1);
322 UFF::Tools::addNonbonded(*mol, cid, types, field, nbrMat);
323 TEST_ASSERT(field->contribs().size() == 1);
324 UFF::Tools::addTorsions(*mol, types, field);
325 TEST_ASSERT(field->contribs().size() == 1);
326
327 mol2 = MolOps::addHs(*mol);
328 TEST_ASSERT(mol2->getNumAtoms() == 6);
329 delete field;
330
331 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol2);
332 TEST_ASSERT(foundAll);
333 TEST_ASSERT(types.size() == mol2->getNumAtoms());
334 auto *conf4 = new Conformer(mol2->getNumAtoms());
335 cid = static_cast<int>(mol2->addConformer(conf4, true));
336
337 field = new ForceFields::ForceField();
338 // add the atomic positions:
339 for (unsigned int i = 0; i < mol2->getNumAtoms(); ++i) {
340 field->positions().push_back(&(conf4->getAtomPos(i)));
341 }
342
343 UFF::Tools::addBonds(*mol2, types, field);
344 TEST_ASSERT(field->contribs().size() == 5);
345
346 nbrMat = UFF::Tools::buildNeighborMatrix(*mol2);
347 UFF::Tools::addAngles(*mol2, types, field);
348 TEST_ASSERT(field->contribs().size() == 12);
349 UFF::Tools::addNonbonded(*mol2, cid, types, field, nbrMat);
350 TEST_ASSERT(field->contribs().size() == 15);
351 UFF::Tools::addTorsions(*mol2, types, field);
352 TEST_ASSERT(field->contribs().size() == 18);
353 delete mol2;
354
355 delete mol;
356 delete field;
357
358 BOOST_LOG(rdErrorLog) << " done" << std::endl;
359 }
360
testUFFBuilder2()361 void testUFFBuilder2() {
362 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
363 BOOST_LOG(rdErrorLog) << " Testing UFF builder+minimization." << std::endl;
364
365 std::string pathName = getenv("RDBASE");
366 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
367 {
368 RWMol *mol = MolFileToMol(pathName + "/small1.mol", false);
369 TEST_ASSERT(mol);
370 MolOps::sanitizeMol(*mol);
371
372 ForceFields::ForceField *field;
373 field = UFF::constructForceField(*mol);
374 TEST_ASSERT(field);
375 field->initialize();
376 int needMore = field->minimize();
377 TEST_ASSERT(!needMore);
378 // std::cout << MolToMolBlock(mol);
379 delete mol;
380 delete field;
381 }
382
383 { // make sure the confId argument works
384 RWMol *mol = MolFileToMol(pathName + "/small1.mol", false);
385 TEST_ASSERT(mol);
386 MolOps::sanitizeMol(*mol);
387 auto *newConf = new Conformer(mol->getConformer());
388 newConf->setId(111);
389 mol->addConformer(newConf, false);
390 RDGeom::Point3D p0 = mol->getConformer().getAtomPos(0);
391 RDGeom::Point3D p1 = mol->getConformer().getAtomPos(1);
392
393 ForceFields::ForceField *field;
394 field = UFF::constructForceField(*mol, 100, 111);
395 TEST_ASSERT(field);
396 field->initialize();
397 int needMore = field->minimize();
398 TEST_ASSERT(!needMore);
399 // std::cout << MolToMolBlock(mol);
400
401 RDGeom::Point3D np0 = mol->getConformer().getAtomPos(0);
402 RDGeom::Point3D np1 = mol->getConformer().getAtomPos(1);
403 TEST_ASSERT(feq(p0.x, np0.x));
404 TEST_ASSERT(feq(p0.y, np0.y));
405 TEST_ASSERT(feq(p0.z, np0.z));
406 TEST_ASSERT(feq(p1.x, np1.x));
407 TEST_ASSERT(feq(p1.y, np1.y));
408 TEST_ASSERT(feq(p1.z, np1.z));
409
410 delete mol;
411 delete field;
412 }
413
414 {
415 RWMol *mol = MolFileToMol(pathName + "/small2.mol", false);
416 TEST_ASSERT(mol);
417 MolOps::sanitizeMol(*mol);
418
419 ForceFields::ForceField *field;
420 field = UFF::constructForceField(*mol);
421 TEST_ASSERT(field);
422 field->initialize();
423 int needMore = field->minimize(150);
424 TEST_ASSERT(!needMore);
425 // std::cout << MolToMolBlock(mol);
426 delete mol;
427 delete field;
428 }
429 {
430 RWMol *mol = MolFileToMol(pathName + "/benzene.mol", false);
431 TEST_ASSERT(mol);
432 MolOps::sanitizeMol(*mol);
433
434 ForceFields::ForceField *field;
435 field = UFF::constructForceField(*mol);
436 TEST_ASSERT(field);
437 field->initialize();
438 int needMore = field->minimize();
439 TEST_ASSERT(!needMore);
440 // std::cout << MolToMolBlock(mol);
441 delete mol;
442 delete field;
443 }
444 {
445 RWMol *mol = MolFileToMol(pathName + "/toluene.mol", false);
446 TEST_ASSERT(mol);
447 MolOps::sanitizeMol(*mol);
448
449 ForceFields::ForceField *field;
450 field = UFF::constructForceField(*mol);
451 TEST_ASSERT(field);
452 field->initialize();
453 int needMore = field->minimize();
454 TEST_ASSERT(!needMore);
455 // std::cout << MolToMolBlock(mol);
456 delete mol;
457 delete field;
458 }
459 {
460 RWMol *mol = MolFileToMol(pathName + "/complex1.mol", false);
461 TEST_ASSERT(mol);
462 MolOps::sanitizeMol(*mol);
463
464 ForceFields::ForceField *field;
465 field = UFF::constructForceField(*mol);
466 TEST_ASSERT(field);
467 field->initialize();
468 int needMore = field->minimize();
469 TEST_ASSERT(!needMore);
470 // std::cout << MolToMolBlock(mol);
471 delete mol;
472 delete field;
473 }
474
475 { // test the convenience function
476 RWMol *mol = MolFileToMol(pathName + "/small1.mol", false);
477 TEST_ASSERT(mol);
478 MolOps::sanitizeMol(*mol);
479
480 UFF::AtomicParamVect types;
481 bool foundAll;
482 boost::shared_array<std::uint8_t> nbrMat;
483 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
484
485 ForceFields::ForceField *field;
486 field = new ForceFields::ForceField();
487
488 // add the atomic positions:
489 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
490 field->positions().push_back(&((mol->getConformer().getAtomPos(i))));
491 }
492
493 UFF::Tools::addBonds(*mol, types, field);
494 nbrMat = UFF::Tools::buildNeighborMatrix(*mol);
495 UFF::Tools::addAngles(*mol, types, field);
496 UFF::Tools::addTorsions(*mol, types, field);
497 // std::cout << field->contribs().size() << std::endl;
498 UFF::Tools::addNonbonded(*mol, 0, types, field, nbrMat);
499 delete field;
500
501 field = UFF::constructForceField(*mol);
502 field->initialize();
503 field->minimize();
504 delete field;
505
506 auto *mol2 = new RWMol(*mol);
507
508 field = UFF::constructForceField(*mol2);
509 TEST_ASSERT(field);
510 field->initialize();
511 int needMore = field->minimize();
512 TEST_ASSERT(!needMore);
513 delete field;
514
515 needMore = UFF::UFFOptimizeMolecule(*mol2).first;
516 TEST_ASSERT(!needMore);
517 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
518 const RDGeom::Point3D p1 = mol->getConformer().getAtomPos(i);
519 const RDGeom::Point3D p2 = mol2->getConformer().getAtomPos(i);
520 TEST_ASSERT(feq(p1.x, p2.x));
521 TEST_ASSERT(feq(p1.y, p2.y));
522 TEST_ASSERT(feq(p1.z, p2.z));
523 }
524
525 delete mol;
526 delete mol2;
527 }
528
529 { // test the convenience function for all confs
530 RWMol *mol = MolFileToMol(pathName + "/complex1.mol", false);
531 TEST_ASSERT(mol);
532 MolOps::sanitizeMol(*mol);
533 auto *mol2 = new RWMol(*mol);
534 auto *newConf = new Conformer(mol->getConformer());
535 newConf->setId(111);
536 mol->addConformer(newConf, false);
537
538 ForceFields::ForceField *field;
539 field = UFF::constructForceField(*mol, 100, 111);
540 TEST_ASSERT(field);
541 field->initialize();
542 int needMore = field->minimize();
543 TEST_ASSERT(!needMore);
544 delete field;
545
546 // the first conf is the same as above,
547 // but we add a second that's already minimized
548 newConf = new Conformer(mol->getConformer(111));
549 newConf->setId(112);
550 mol2->addConformer(newConf, false);
551
552 std::vector<std::pair<int, double>> res;
553 UFF::UFFOptimizeMoleculeConfs(*mol2, res);
554 TEST_ASSERT(res.size() == 2);
555 TEST_ASSERT(!res[0].first);
556 TEST_ASSERT(!res[1].first);
557 // we expect the energy to go down at least a little bit.
558 TEST_ASSERT(res[1].second < res[0].second);
559
560 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
561 const RDGeom::Point3D p1 = mol->getConformer(111).getAtomPos(i);
562 const RDGeom::Point3D p2 = mol2->getConformer(0).getAtomPos(i);
563 TEST_ASSERT(feq(p1.x, p2.x));
564 TEST_ASSERT(feq(p1.y, p2.y));
565 TEST_ASSERT(feq(p1.z, p2.z));
566 }
567
568 delete mol;
569 delete mol2;
570 }
571
572 BOOST_LOG(rdErrorLog) << " done" << std::endl;
573 }
574
575 #ifdef RDK_TEST_MULTITHREADED
testUFFBatch()576 void testUFFBatch() {}
577 #else
testUFFBatch()578 void testUFFBatch() {
579 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
580 BOOST_LOG(rdErrorLog)
581 << " Testing bulk UFF (regression to check that things run)."
582 << std::endl;
583
584 ROMol *mol;
585 std::string key;
586
587 ForceFields::ForceField *field;
588
589 std::string pathName = getenv("RDBASE");
590 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
591 SDMolSupplier suppl(pathName + "/bulk.sdf", false);
592
593 int count = 0;
594 mol = suppl.next();
595 while (mol && !suppl.atEnd()) {
596 count++;
597 MolOps::sanitizeMol(*(RWMol *)mol);
598 std::string origMolBlock = MolToMolBlock(*mol);
599
600 BOOST_LOG(rdErrorLog) << "Mol:" << count << std::endl;
601 try {
602 field = UFF::constructForceField(*mol);
603 } catch (...) {
604 field = nullptr;
605 }
606 if (field) {
607 field->initialize();
608 int failed = field->minimize(500);
609 if (failed) {
610 BOOST_LOG(rdErrorLog)
611 << " not converged (code=" << failed << ")" << std::endl;
612 std::cout << origMolBlock << "$$$$" << std::endl;
613 std::cout << MolToMolBlock(*mol) << "$$$$" << std::endl;
614 }
615 delete field;
616 }
617 delete mol;
618 mol = suppl.next();
619 }
620 delete mol;
621
622 BOOST_LOG(rdErrorLog) << " done" << std::endl;
623 }
624 #endif
625
testUFFBuilderSpecialCases()626 void testUFFBuilderSpecialCases() {
627 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
628 BOOST_LOG(rdErrorLog) << " Testing UFF special cases." << std::endl;
629
630 RWMol *mol;
631 std::string key;
632 int needMore;
633 RDGeom::Point3D v1, v2;
634 ForceFields::ForceField *field;
635
636 std::string pathName = getenv("RDBASE");
637 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
638 // ----------
639 // Trigonal bipyramid
640 // ----------
641 mol = MolFileToMol(pathName + "/tbp.mol", false);
642 TEST_ASSERT(mol);
643 MolOps::sanitizeMol(*mol);
644
645 const Conformer &conf = mol->getConformer();
646 field = UFF::constructForceField(*mol);
647 TEST_ASSERT(field);
648 field->initialize();
649 needMore = field->minimize(200, 1e-8, 1e-4);
650 TEST_ASSERT(!needMore);
651 v1 = conf.getAtomPos(0).directionVector(conf.getAtomPos(1));
652 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(2));
653 TEST_ASSERT(feq(v1.dotProduct(v2), -1.0, 1e-3));
654 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(3));
655 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
656 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(4));
657 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
658 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(5));
659 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
660
661 v1 = conf.getAtomPos(0).directionVector(conf.getAtomPos(2));
662 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(3));
663 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
664 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(4));
665 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
666 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(5));
667 TEST_ASSERT(feq(v1.dotProduct(v2), 0.0, 1e-3));
668
669 v1 = conf.getAtomPos(0).directionVector(conf.getAtomPos(3));
670 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(4));
671 TEST_ASSERT(feq(v1.dotProduct(v2), -0.5, 1e-3));
672 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(5));
673 TEST_ASSERT(feq(v1.dotProduct(v2), -0.5, 1e-3));
674
675 v1 = conf.getAtomPos(0).directionVector(conf.getAtomPos(4));
676 v2 = conf.getAtomPos(0).directionVector(conf.getAtomPos(5));
677 TEST_ASSERT(feq(v1.dotProduct(v2), -0.5, 1e-3));
678
679 delete mol;
680 delete field;
681
682 BOOST_LOG(rdErrorLog) << " done" << std::endl;
683 }
684
testIssue239()685 void testIssue239() {
686 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
687 BOOST_LOG(rdErrorLog) << " Testing Issue239." << std::endl;
688
689 RWMol *mol;
690 int needMore;
691 (void)needMore; // add test later
692 ForceFields::ForceField *field;
693 double e1, e2;
694
695 std::string pathName = getenv("RDBASE");
696 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
697 mol = MolFileToMol(pathName + "/Issue239.mol", false);
698 TEST_ASSERT(mol);
699 MolOps::sanitizeMol(*mol);
700
701 field = UFF::constructForceField(*mol);
702 TEST_ASSERT(field);
703 field->initialize();
704 needMore = field->minimize(200, 1e-6, 1e-3);
705 e1 = field->calcEnergy();
706 needMore = field->minimize(200, 1e-6, 1e-3);
707 e2 = field->calcEnergy();
708 TEST_ASSERT(feq(e2, e1, 0.1));
709
710 delete mol;
711 delete field;
712
713 BOOST_LOG(rdErrorLog) << " done" << std::endl;
714 }
715 #endif
716
testCalcEnergyPassedCoords()717 void testCalcEnergyPassedCoords() {
718 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
719 BOOST_LOG(rdErrorLog) << " Testing calcEnergy with passed coords."
720 << std::endl;
721
722 RWMol *mol;
723 ForceFields::ForceField *field;
724 double e1, e2, e3;
725
726 std::string pathName = getenv("RDBASE");
727 pathName += "/Code/GraphMol/ForceFieldHelpers/MMFF/test_data";
728 mol = MolFileToMol(pathName + "/Issue239.mol", false);
729 TEST_ASSERT(mol);
730 MolOps::sanitizeMol(*mol);
731
732 field = UFF::constructForceField(*mol);
733 TEST_ASSERT(field);
734 field->initialize();
735 auto *savedPos = new double[3 * field->numPoints()];
736 size_t i = 0;
737 for (const auto pptr : field->positions()) {
738 for (size_t j = 0; j < 3; ++j) {
739 savedPos[i++] = (*pptr)[j];
740 }
741 }
742 e1 = field->calcEnergy();
743 field->minimize(10000, 1.0e-6, 1.0e-3);
744 e2 = field->calcEnergy();
745 TEST_ASSERT(e2 < e1);
746 e3 = field->calcEnergy(savedPos);
747 TEST_ASSERT(feq(e3, e1, 0.01));
748
749 delete[] savedPos;
750 delete mol;
751 delete field;
752
753 BOOST_LOG(rdErrorLog) << " done" << std::endl;
754 }
755
testCalcGrad()756 void testCalcGrad() {
757 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
758 BOOST_LOG(rdErrorLog) << " Testing calcGrad." << std::endl;
759
760 RWMol *mol;
761 ForceFields::ForceField *field;
762
763 std::string pathName = getenv("RDBASE");
764 pathName += "/Code/GraphMol/ForceFieldHelpers/MMFF/test_data";
765 mol = MolFileToMol(pathName + "/Issue239.mol", false);
766 TEST_ASSERT(mol);
767 MolOps::sanitizeMol(*mol);
768
769 field = UFF::constructForceField(*mol);
770 TEST_ASSERT(field);
771 field->initialize();
772 size_t l = 3 * field->numPoints();
773 auto *savedPos = new double[l];
774 auto *grad1 = new double[l];
775 auto *grad2 = new double[l];
776 size_t i = 0;
777 for (const auto pptr : field->positions()) {
778 for (size_t j = 0; j < 3; ++j) {
779 savedPos[i++] = (*pptr)[j];
780 }
781 }
782 TEST_ASSERT(i == l);
783
784 std::memset(grad1, 0, l * sizeof(double));
785 field->calcGrad(grad1);
786 for (i = 0; i < l; ++i) {
787 TEST_ASSERT(!feq(grad1[i], 0.0, 0.001));
788 }
789
790 field->minimize(10000, 1.0e-6, 1.0e-3);
791 std::memset(grad2, 0, l * sizeof(double));
792 field->calcGrad(grad2);
793 for (i = 0; i < l; ++i) {
794 TEST_ASSERT(feq(grad2[i], 0.0, 0.001));
795 }
796
797 field->initialize();
798 std::memset(grad2, 0, l * sizeof(double));
799 field->calcGrad(savedPos, grad2);
800 for (i = 0; i < l; ++i) {
801 TEST_ASSERT(feq(grad1[i], grad2[i], 0.001));
802 }
803
804 delete[] savedPos;
805 delete[] grad1;
806 delete[] grad2;
807 delete mol;
808 delete field;
809
810 BOOST_LOG(rdErrorLog) << " done" << std::endl;
811 }
testIssue242()812 void testIssue242() {
813 #if 0
814 // FIX: Changes to the forcefield (connected to Issue 408) have
815 // made it so that this particular problem no longer manifests
816 // in this molecule/embedding. A new test case is needed.
817 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
818 BOOST_LOG(rdErrorLog) << " Testing Issue242." << std::endl;
819
820 RWMol *mol,*mol2;
821 int needMore;
822 ForceFields::ForceField *field=0,*field2=0;
823 std::string mb1,mb2;
824 double e1,e2;
825
826 std::string pathName=getenv("RDBASE");
827 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
828
829 mol = MolFileToMol(pathName+"/Issue242.mol");
830 TEST_ASSERT(mol);
831
832 mol2 = MolFileToMol(pathName+"/Issue242.mol");
833 TEST_ASSERT(mol2);
834
835 TEST_ASSERT(DGeomHelpers::EmbedMolecule(*mol2,30,2300)>=0);
836 mb1 = MolToMolBlock(*mol);
837 mb2 = MolToMolBlock(*mol2);
838
839 //BOOST_LOG(rdInfoLog) << "\nMol1\n" << mb1 << std::endl;
840 //BOOST_LOG(rdInfoLog) << "\nMol2\n" << mb2 << std::endl;
841
842
843 field=UFF::constructForceField(*mol);
844 TEST_ASSERT(field);
845 field->initialize();
846 field2=UFF::constructForceField(*mol2);
847 TEST_ASSERT(field2);
848 field2->initialize();
849 e1 = field->calcEnergy();
850 e2 = field2->calcEnergy();
851 BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
852 BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
853 //TEST_ASSERT(feq(e2,e1,0.1));
854
855 needMore = field->minimize(600,1e-4);
856 TEST_ASSERT(!needMore)
857 needMore = field2->minimize(600,1e-4);
858 TEST_ASSERT(!needMore)
859 e1 = field->calcEnergy();
860 e2 = field2->calcEnergy();
861 BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
862 BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
863 TEST_ASSERT(feq(e2,e1,1.0));
864
865 needMore = field->minimize(600,1e-4);
866 TEST_ASSERT(!needMore)
867 needMore = field2->minimize(600,1e-4);
868 TEST_ASSERT(!needMore)
869 e1 = field->calcEnergy();
870 e2 = field2->calcEnergy();
871 BOOST_LOG(rdInfoLog) << "rE1: " << e1 << std::endl;
872 BOOST_LOG(rdInfoLog) << "rE2: " << e2 << std::endl;
873 TEST_ASSERT(feq(e2,e1,1.0));
874
875 delete mol;
876 delete mol2;
877 delete field;
878 delete field2;
879
880 mol = MolFileToMol(pathName+"/Issue242-2.mol");
881 TEST_ASSERT(mol);
882
883 mol2 = MolFileToMol(pathName+"/Issue242-2.mol");
884 TEST_ASSERT(mol2);
885
886 TEST_ASSERT(DGeomHelpers::EmbedMolecule(*mol2,30,2370)>=0);
887 mb1 = MolToMolBlock(*mol);
888 mb2 = MolToMolBlock(*mol2);
889
890 //std::cout << mb1 << "------\n";
891 //std::cout << mb2 << "------\n";
892
893 field=UFF::constructForceField(*mol);
894 TEST_ASSERT(field);
895 field->initialize();
896 field2=UFF::constructForceField(*mol2);
897 TEST_ASSERT(field2);
898 field2->initialize();
899 e1 = field->calcEnergy();
900 e2 = field2->calcEnergy();
901 BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
902 BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
903 //TEST_ASSERT(feq(e2,e1,0.1));
904
905 needMore = field->minimize(200,1e-4);
906 needMore = field2->minimize(200,1e-4);
907 e1 = field->calcEnergy();
908 e2 = field2->calcEnergy();
909 BOOST_LOG(rdInfoLog) << "E1: " << e1 << std::endl;
910 BOOST_LOG(rdInfoLog) << "E2: " << e2 << std::endl;
911 TEST_ASSERT(feq(e2,e1,0.1));
912
913 delete mol;
914 delete mol2;
915 delete field;
916 delete field2;
917
918
919 BOOST_LOG(rdErrorLog) << " done" << std::endl;
920 #endif
921 }
922
testSFIssue1653802()923 void testSFIssue1653802() {
924 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
925 BOOST_LOG(rdErrorLog) << " Testing SFIssue1653802." << std::endl;
926
927 RWMol *mol;
928 int needMore;
929 ForceFields::ForceField *field;
930
931 std::string pathName = getenv("RDBASE");
932 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
933
934 mol = MolFileToMol(pathName + "/cyclobutadiene.mol", false);
935 TEST_ASSERT(mol);
936 MolOps::sanitizeMol(*mol);
937
938 UFF::AtomicParamVect types;
939 bool foundAll;
940 boost::shared_array<std::uint8_t> nbrMat;
941 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
942 TEST_ASSERT(foundAll);
943 TEST_ASSERT(types.size() == mol->getNumAtoms());
944 field = new ForceFields::ForceField();
945 // add the atomic positions:
946 for (unsigned int i = 0; i < mol->getNumAtoms(); ++i) {
947 field->positions().push_back(&((mol->getConformer().getAtomPos(i))));
948 }
949
950 UFF::Tools::addBonds(*mol, types, field);
951 TEST_ASSERT(field->contribs().size() == 8);
952
953 nbrMat = UFF::Tools::buildNeighborMatrix(*mol);
954 UFF::Tools::addAngles(*mol, types, field);
955 TEST_ASSERT(field->contribs().size() == 20);
956 UFF::Tools::addTorsions(*mol, types, field);
957 // std::cout << field->contribs().size() << std::endl;
958 TEST_ASSERT(field->contribs().size() == 36);
959 UFF::Tools::addNonbonded(*mol, 0, types, field, nbrMat);
960 delete field;
961
962 field = UFF::constructForceField(*mol);
963 field->initialize();
964 needMore = field->minimize(200, 1e-6, 1e-3);
965 TEST_ASSERT(!needMore);
966
967 delete mol;
968 delete field;
969 BOOST_LOG(rdErrorLog) << " done" << std::endl;
970 }
971
testSFIssue2378119()972 void testSFIssue2378119() {
973 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
974 BOOST_LOG(rdErrorLog) << " Testing SFIssue2378119." << std::endl;
975
976 std::string pathName = getenv("RDBASE");
977 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
978 {
979 RWMol *mol = MolFileToMol(pathName + "/Issue2378119.mol");
980 TEST_ASSERT(mol);
981 ForceFields::ForceField *field = UFF::constructForceField(*mol);
982 TEST_ASSERT(field);
983 field->initialize();
984 double e1 = field->calcEnergy();
985 TEST_ASSERT(e1 > 0.0 && e1 < 1e12);
986
987 int needMore = field->minimize(200, 1e-6, 1e-3);
988 TEST_ASSERT(!needMore);
989 double e2 = field->calcEnergy();
990 TEST_ASSERT(e2 < e1);
991
992 delete mol;
993 delete field;
994 }
995 {
996 RWMol *mol = MolFileToMol(pathName + "/Issue2378119.2.mol");
997 TEST_ASSERT(mol);
998 ForceFields::ForceField *field = UFF::constructForceField(*mol);
999 TEST_ASSERT(field);
1000 field->initialize();
1001 double e1 = field->calcEnergy();
1002 TEST_ASSERT(e1 == 0.0);
1003
1004 int needMore = field->minimize(200, 1e-6, 1e-3);
1005 TEST_ASSERT(!needMore);
1006 double e2 = field->calcEnergy();
1007 TEST_ASSERT(e2 == e1);
1008
1009 delete mol;
1010 delete field;
1011 }
1012 {
1013 RWMol *mol = MolFileToMol(pathName + "/Issue2378119.2.mol");
1014 TEST_ASSERT(mol);
1015 ForceFields::ForceField *field =
1016 UFF::constructForceField(*mol, 100.0, -1, false);
1017 TEST_ASSERT(field);
1018 field->initialize();
1019 double e1 = field->calcEnergy();
1020 TEST_ASSERT(e1 > 0.0 && e1 < 1e12);
1021
1022 int needMore = field->minimize(200, 1e-6, 1e-3);
1023 TEST_ASSERT(!needMore);
1024 double e2 = field->calcEnergy();
1025 TEST_ASSERT(e2 < e1);
1026
1027 delete mol;
1028 delete field;
1029 }
1030 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1031 }
1032
testMissingParams()1033 void testMissingParams() {
1034 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1035 BOOST_LOG(rdErrorLog) << " Testing handling missing parameters."
1036 << std::endl;
1037
1038 {
1039 UFF::AtomicParamVect types;
1040 bool foundAll;
1041
1042 ROMol *mol = SmilesToMol("[Cu](C)(C)(C)(C)C");
1043 TEST_ASSERT(mol);
1044
1045 ROMol *mol2 = MolOps::addHs(*mol);
1046 delete mol;
1047
1048 TEST_ASSERT(DGeomHelpers::EmbedMolecule(*mol2) >= 0);
1049
1050 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol2);
1051 TEST_ASSERT(!foundAll);
1052 TEST_ASSERT(types.size() == mol2->getNumAtoms());
1053 TEST_ASSERT(!types[0]);
1054
1055 // make sure we can optimize anyway:
1056 ForceFields::ForceField *field = UFF::constructForceField(*mol2, types);
1057 TEST_ASSERT(field);
1058 field->initialize();
1059 double e1 = field->calcEnergy();
1060 int needMore = field->minimize();
1061 TEST_ASSERT(needMore);
1062 double e2 = field->calcEnergy();
1063 TEST_ASSERT(e2 < e1);
1064 // std::cerr<<" DE: "<<e1<<" -> "<<e2<<std::endl;
1065 delete mol2;
1066 delete field;
1067 }
1068
1069 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1070 }
testSFIssue3009337()1071 void testSFIssue3009337() {
1072 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1073 BOOST_LOG(rdErrorLog) << " Testing SFIssue3009337." << std::endl;
1074
1075 std::string pathName = getenv("RDBASE");
1076 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
1077 {
1078 RWMol *mol = MolFileToMol(pathName + "/Issue3009337.mol", true, false);
1079 TEST_ASSERT(mol);
1080 ForceFields::ForceField *field = UFF::constructForceField(*mol);
1081 TEST_ASSERT(field);
1082 field->initialize();
1083 double e1 = field->calcEnergy();
1084 TEST_ASSERT(e1 > 0.0 && e1 < 1e12);
1085
1086 int needMore = field->minimize(200, 1e-6, 1e-3);
1087 TEST_ASSERT(!needMore);
1088 double e2 = field->calcEnergy();
1089 TEST_ASSERT(e2 < e1);
1090
1091 delete mol;
1092 delete field;
1093 }
1094 {
1095 RWMol *mol = MolFileToMol(pathName + "/Issue3009337.2.mol", true, false);
1096 TEST_ASSERT(mol);
1097 ForceFields::ForceField *field = UFF::constructForceField(*mol);
1098 TEST_ASSERT(field);
1099 field->initialize();
1100 double e1 = field->calcEnergy();
1101 TEST_ASSERT(e1 > 0.0 && e1 < 1e12);
1102
1103 int needMore = field->minimize(200, 1e-6, 1e-3);
1104 TEST_ASSERT(!needMore);
1105 double e2 = field->calcEnergy();
1106 TEST_ASSERT(e2 < e1);
1107
1108 delete mol;
1109 delete field;
1110 }
1111 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1112 }
testGitHubIssue62()1113 void testGitHubIssue62() {
1114 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1115 BOOST_LOG(rdErrorLog) << " Testing GitHubIssue62." << std::endl;
1116
1117 std::string pathName = getenv("RDBASE");
1118 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
1119 {
1120 double energyValues[] = {
1121 38.687, 174.698, 337.986, 115.248, 2.482, 1.918, 10.165, 99.492,
1122 41.016, 267.236, 15.747, 203.398, 206.852, 20.044, 218.879, 79.614};
1123 SmilesMolSupplier smiSupplier(pathName + "/Issue62.smi");
1124 SDWriter *sdfWriter = new SDWriter(pathName + "/Issue62.sdf");
1125 for (unsigned int i = 0; i < smiSupplier.length(); ++i) {
1126 auto *tmp = smiSupplier[i];
1127 ROMol *mol = MolOps::addHs(*tmp);
1128 delete tmp;
1129 TEST_ASSERT(mol);
1130 std::string molName = "";
1131 if (mol->hasProp(common_properties::_Name)) {
1132 mol->getProp(common_properties::_Name, molName);
1133 }
1134 DGeomHelpers::EmbedMolecule(*mol);
1135 ForceFields::ForceField *field = UFF::constructForceField(*mol);
1136 TEST_ASSERT(field);
1137 field->initialize();
1138 int needMore = field->minimize(200, 1.e-6, 1.e-3);
1139 TEST_ASSERT(!needMore);
1140 sdfWriter->write(*mol);
1141 double e = field->calcEnergy();
1142 BOOST_LOG(rdErrorLog) << molName << " " << e << std::endl;
1143 TEST_ASSERT(fabs(e - energyValues[i]) < 1.);
1144 delete field;
1145 delete mol;
1146 }
1147 sdfWriter->close();
1148 delete sdfWriter;
1149 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1150 }
1151 }
testUFFParamGetters()1152 void testUFFParamGetters() {
1153 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1154 BOOST_LOG(rdErrorLog) << " Test UFF force-field parameter getters."
1155 << std::endl;
1156 {
1157 ROMol *mol = SmilesToMol("c1ccccc1CCNN");
1158 TEST_ASSERT(mol);
1159 ROMol *molH = MolOps::addHs(*mol);
1160 delete mol;
1161 TEST_ASSERT(molH);
1162 ForceFields::UFF::UFFBond uffBondStretchParams;
1163 TEST_ASSERT(
1164 UFF::getUFFBondStretchParams(*molH, 6, 7, uffBondStretchParams));
1165 TEST_ASSERT(((int)std::round(uffBondStretchParams.kb * 1000) == 699592) &&
1166 ((int)std::round(uffBondStretchParams.r0 * 1000) == 1514));
1167 TEST_ASSERT(
1168 !UFF::getUFFBondStretchParams(*molH, 0, 7, uffBondStretchParams));
1169 ForceFields::UFF::UFFAngle uffAngleBendParams;
1170 TEST_ASSERT(UFF::getUFFAngleBendParams(*molH, 6, 7, 8, uffAngleBendParams));
1171 TEST_ASSERT(((int)std::round(uffAngleBendParams.ka * 1000) == 303297) &&
1172 ((int)std::round(uffAngleBendParams.theta0 * 1000) == 109470));
1173 TEST_ASSERT(
1174 !UFF::getUFFAngleBendParams(*molH, 0, 7, 8, uffAngleBendParams));
1175 ForceFields::UFF::UFFTor uffTorsionParams;
1176 TEST_ASSERT(UFF::getUFFTorsionParams(*molH, 6, 7, 8, 9, uffTorsionParams));
1177 TEST_ASSERT(((int)std::round(uffTorsionParams.V * 1000) == 976));
1178 TEST_ASSERT(!UFF::getUFFTorsionParams(*molH, 0, 7, 8, 9, uffTorsionParams));
1179 ForceFields::UFF::UFFInv uffInversionParams;
1180 TEST_ASSERT(
1181 UFF::getUFFInversionParams(*molH, 6, 5, 4, 0, uffInversionParams));
1182 TEST_ASSERT(((int)std::round(uffInversionParams.K * 1000) == 2000));
1183 TEST_ASSERT(
1184 !UFF::getUFFInversionParams(*molH, 6, 5, 4, 1, uffInversionParams));
1185 ForceFields::UFF::UFFVdW uffVdWParams;
1186 TEST_ASSERT(UFF::getUFFVdWParams(*molH, 0, 9, uffVdWParams));
1187 TEST_ASSERT(((int)std::round(uffVdWParams.x_ij * 1000) == 3754) &&
1188 ((int)std::round(uffVdWParams.D_ij * 1000) == 85));
1189 delete molH;
1190 }
1191 }
1192
1193 #ifdef RDK_TEST_MULTITHREADED
1194 namespace {
runblock_uff(const std::vector<ROMol * > & mols,const std::vector<double> & energies,unsigned int count,unsigned int idx)1195 void runblock_uff(const std::vector<ROMol *> &mols,
1196 const std::vector<double> &energies, unsigned int count,
1197 unsigned int idx) {
1198 for (unsigned int rep = 0; rep < 200; ++rep) {
1199 for (unsigned int i = 0; i < mols.size(); ++i) {
1200 if (i % count != idx) {
1201 continue;
1202 }
1203 ROMol *mol = mols[i];
1204 ForceFields::ForceField *field = nullptr;
1205 if (!(rep % 100)) {
1206 BOOST_LOG(rdErrorLog) << "Rep: " << rep << " Mol:" << i << std::endl;
1207 }
1208 try {
1209 field = UFF::constructForceField(*mol);
1210 } catch (...) {
1211 field = nullptr;
1212 }
1213 TEST_ASSERT(field);
1214 field->initialize();
1215 int failed = field->minimize(500);
1216 TEST_ASSERT(!failed);
1217 double eng = field->calcEnergy();
1218 TEST_ASSERT(feq(eng, energies[i]));
1219 delete field;
1220 }
1221 }
1222 }
1223 } // namespace
1224 #include <thread>
1225 #include <future>
testUFFMultiThread()1226 void testUFFMultiThread() {
1227 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1228 BOOST_LOG(rdErrorLog) << " Test UFF multithreading" << std::endl;
1229
1230 // ForceFields::ForceField *field;
1231
1232 std::string pathName = getenv("RDBASE");
1233 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
1234 SDMolSupplier suppl(pathName + "/bulk.sdf");
1235 std::vector<ROMol *> mols;
1236 while (!suppl.atEnd() && mols.size() < 100) {
1237 ROMol *mol = nullptr;
1238 try {
1239 mol = suppl.next();
1240 } catch (...) {
1241 continue;
1242 }
1243 if (!mol) {
1244 continue;
1245 }
1246 mols.push_back(mol);
1247 }
1248
1249 std::cerr << "generating reference data" << std::endl;
1250 std::vector<double> energies(mols.size(), 0.0);
1251 for (unsigned int i = 0; i < mols.size(); ++i) {
1252 ROMol mol(*mols[i]);
1253 ForceFields::ForceField *field = nullptr;
1254 try {
1255 field = UFF::constructForceField(mol);
1256 } catch (...) {
1257 field = nullptr;
1258 }
1259 TEST_ASSERT(field);
1260 field->initialize();
1261 int failed = field->minimize(500);
1262 TEST_ASSERT(!failed);
1263 energies[i] = field->calcEnergy();
1264 delete field;
1265 }
1266
1267 std::vector<std::future<void>> tg;
1268
1269 std::cerr << "processing" << std::endl;
1270 unsigned int count = 4;
1271 for (unsigned int i = 0; i < count; ++i) {
1272 std::cerr << " launch :" << i << std::endl;
1273 std::cerr.flush();
1274 tg.emplace_back(
1275 std::async(std::launch::async, runblock_uff, mols, energies, count, i));
1276 }
1277 for (auto &fut : tg) {
1278 fut.get();
1279 }
1280
1281 for (auto *mol : mols) {
1282 delete mol;
1283 }
1284 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1285 }
1286
testUFFMultiThread2()1287 void testUFFMultiThread2() {
1288 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1289 BOOST_LOG(rdErrorLog) << " Test UFF multithreading2" << std::endl;
1290
1291 std::string pathName = getenv("RDBASE");
1292 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
1293 SDMolSupplier suppl(pathName + "/bulk.sdf");
1294 ROMol *m = suppl[4];
1295 TEST_ASSERT(m);
1296 auto *om = new ROMol(*m);
1297 for (unsigned int i = 0; i < 200; ++i) {
1298 m->addConformer(new Conformer(m->getConformer()), true);
1299 }
1300 std::vector<std::pair<int, double>> res;
1301
1302 UFF::UFFOptimizeMolecule(*om);
1303 UFF::UFFOptimizeMoleculeConfs(*m, res, 0);
1304 for (unsigned int i = 1; i < res.size(); ++i) {
1305 TEST_ASSERT(!res[i].first);
1306 TEST_ASSERT(feq(res[i].second, res[0].second, .00001));
1307 }
1308 for (unsigned int i = 0; i < m->getNumAtoms(); ++i) {
1309 RDGeom::Point3D p0 = om->getConformer().getAtomPos(i);
1310 RDGeom::Point3D np0 = m->getConformer().getAtomPos(i);
1311 TEST_ASSERT(feq(p0.x, np0.x));
1312 TEST_ASSERT(feq(p0.y, np0.y));
1313 TEST_ASSERT(feq(p0.z, np0.z));
1314 np0 =
1315 m->getConformer(11).getAtomPos(i); // pick some random other conformer
1316 TEST_ASSERT(feq(p0.x, np0.x));
1317 TEST_ASSERT(feq(p0.y, np0.y));
1318 TEST_ASSERT(feq(p0.z, np0.z));
1319 }
1320 delete m;
1321 delete om;
1322 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1323 }
1324
testUFFMultiThread3()1325 void testUFFMultiThread3() {
1326 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1327 BOOST_LOG(rdErrorLog) << " Test UFF multithreading3" << std::endl;
1328
1329 std::string pathName = getenv("RDBASE");
1330 pathName += "/Code/GraphMol/ForceFieldHelpers/UFF/test_data";
1331 SDMolSupplier suppl(pathName + "/bulk.sdf");
1332 ROMol *m = suppl[4];
1333 TEST_ASSERT(m);
1334 auto *om = new ROMol(*m);
1335 for (unsigned int i = 0; i < 200; ++i) {
1336 m->addConformer(new Conformer(m->getConformer()), true);
1337 }
1338 std::vector<std::pair<int, double>> res;
1339
1340 ForceFields::ForceField *omField = UFF::constructForceField(*om);
1341 TEST_ASSERT(omField);
1342 omField->initialize();
1343 ForceFields::ForceField *mField = UFF::constructForceField(*m);
1344 TEST_ASSERT(mField);
1345 mField->initialize();
1346
1347 ForceFieldsHelper::OptimizeMolecule(*omField);
1348 ForceFieldsHelper::OptimizeMoleculeConfs(*m, *mField, res, 0);
1349 for (unsigned int i = 1; i < res.size(); ++i) {
1350 TEST_ASSERT(!res[i].first);
1351 TEST_ASSERT(feq(res[i].second, res[0].second, .00001));
1352 }
1353 for (unsigned int i = 0; i < m->getNumAtoms(); ++i) {
1354 RDGeom::Point3D p0 = om->getConformer().getAtomPos(i);
1355 RDGeom::Point3D np0 = m->getConformer().getAtomPos(i);
1356 TEST_ASSERT(feq(p0.x, np0.x));
1357 TEST_ASSERT(feq(p0.y, np0.y));
1358 TEST_ASSERT(feq(p0.z, np0.z));
1359 np0 =
1360 m->getConformer(11).getAtomPos(i); // pick some random other conformer
1361 TEST_ASSERT(feq(p0.x, np0.x));
1362 TEST_ASSERT(feq(p0.y, np0.y));
1363 TEST_ASSERT(feq(p0.z, np0.z));
1364 }
1365 delete m;
1366 delete om;
1367 delete mField;
1368 delete omField;
1369 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1370 }
1371 #endif
1372
testGitHubIssue613()1373 void testGitHubIssue613() {
1374 BOOST_LOG(rdErrorLog) << "-------------------------------------" << std::endl;
1375 BOOST_LOG(rdErrorLog) << " Test Github Issue 613: UFF Atom type not "
1376 "properly assigned to lanthanides."
1377 << std::endl;
1378 {
1379 ROMol *mol =
1380 SmilesToMol("[Eu+3]123456.[Cl]1.[Cl]2.[Cl]3.[Cl]4.[Cl]5.[Cl]6");
1381 TEST_ASSERT(mol);
1382 mol->getAtomWithIdx(0)->setHybridization(Atom::SP3D2);
1383
1384 UFF::AtomicParamVect types;
1385 bool foundAll;
1386 boost::tie(types, foundAll) = UFF::getAtomTypes(*mol);
1387 TEST_ASSERT(foundAll);
1388 TEST_ASSERT(types.size() == mol->getNumAtoms());
1389
1390 ForceFields::UFF::ParamCollection *params =
1391 ForceFields::UFF::ParamCollection::getParams();
1392 const ForceFields::UFF::AtomicParams *ap = (*params)("Eu6+3");
1393 TEST_ASSERT(ap);
1394 TEST_ASSERT(ap->r1 == types[0]->r1);
1395 TEST_ASSERT(ap->theta0 == types[0]->theta0);
1396 delete mol;
1397 }
1398 BOOST_LOG(rdErrorLog) << " done" << std::endl;
1399 }
1400
1401 //-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
1402 //
1403 //-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
main()1404 int main() {
1405 RDLog::InitLogs();
1406 // we get a ton of warnings here about missing Hs... disable them
1407 boost::logging::disable_logs("rdApp.warning");
1408
1409 #if 1
1410 testUFFTyper1();
1411 testUFFTyper2();
1412 testUFFBuilder1();
1413 testUFFBuilder2();
1414 testUFFBatch();
1415 testUFFBuilderSpecialCases();
1416 testIssue239();
1417 testCalcEnergyPassedCoords();
1418 testCalcGrad();
1419 testIssue242();
1420 testSFIssue1653802();
1421 testSFIssue2378119();
1422 testUFFParamGetters();
1423 testMissingParams();
1424 testSFIssue3009337();
1425 #ifdef RDK_TEST_MULTITHREADED
1426 testUFFMultiThread();
1427 testUFFMultiThread2();
1428 testUFFMultiThread3();
1429 #endif
1430 testGitHubIssue62();
1431 #endif
1432 testGitHubIssue613();
1433 }
1434