1 // $Id$
2 //
3 // Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC
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 #include <RDGeneral/test.h>
12 #include <iostream>
13 #include <iomanip>
14 #include <cmath>
15 #include <RDGeneral/Invariant.h>
16 #include <RDGeneral/utils.h>
17 #include <Geometry/point.h>
18
19 #include <ForceField/ForceField.h>
20 #include <ForceField/UFF/Params.h>
21 #include <ForceField/UFF/BondStretch.h>
22 #include <ForceField/UFF/AngleBend.h>
23 #include <ForceField/UFF/Nonbonded.h>
24 #include <ForceField/UFF/TorsionAngle.h>
25 #include <ForceField/UFF/DistanceConstraint.h>
26 #include <ForceField/UFF/AngleConstraint.h>
27 #include <ForceField/UFF/TorsionConstraint.h>
28 #include <ForceField/UFF/PositionConstraint.h>
29
30 #include <GraphMol/RDKitBase.h>
31 #include <GraphMol/FileParsers/FileParsers.h>
32 #include <GraphMol/ForceFieldHelpers/UFF/Builder.h>
33 #include <GraphMol/MolTransforms/MolTransforms.h>
34
35 using namespace RDKit;
36
test1()37 void test1() {
38 std::cerr << "-------------------------------------" << std::endl;
39 std::cerr << "Unit tests for force field basics." << std::endl;
40
41 ForceFields::ForceField ff;
42 TEST_ASSERT(ff.dimension() == 3);
43
44 RDGeom::Point3D p1(0, 0, 0), p2(1, 0, 0), p3(2, 0, 0), p4(0, 1, 0);
45 RDGeom::PointPtrVect &ps = ff.positions();
46 ps.push_back(&p1);
47 ps.push_back(&p2);
48 ps.push_back(&p3);
49 ps.push_back(&p4);
50
51 #if 0
52 RDGeom::Point3D f1,f2,f3,f4;
53 RDGeom::PointPtrVect &fs=ff.forces();
54 fs.push_back(&f1);
55 fs.push_back(&f2);
56 fs.push_back(&f3);
57 fs.push_back(&f4);
58 #endif
59 TEST_ASSERT(ff.positions().size() == 4);
60 // TEST_ASSERT(ff.forces().size()==4);
61
62 ff.initialize();
63
64 TEST_ASSERT(RDKit::feq(ff.distance(0, 1), 1.0));
65 TEST_ASSERT(RDKit::feq(ff.distance(1, 0), 1.0));
66 TEST_ASSERT(RDKit::feq(ff.distance(0, 0), 0.0));
67 TEST_ASSERT(RDKit::feq(ff.distance(0, 2), 2.0));
68 TEST_ASSERT(RDKit::feq(ff.distance(2, 0), 2.0));
69 TEST_ASSERT(RDKit::feq(ff.distance(0, 3), 1.0));
70 TEST_ASSERT(RDKit::feq(ff.distance(3, 0), 1.0));
71 TEST_ASSERT(RDKit::feq(ff.distance(3, 3), 0.0));
72 TEST_ASSERT(RDKit::feq(ff.distance(1, 2), 1.0));
73 TEST_ASSERT(RDKit::feq(ff.distance(2, 1), 1.0));
74
75 std::cerr << " done" << std::endl;
76 }
77
testUFF1()78 void testUFF1() {
79 std::cerr << "-------------------------------------" << std::endl;
80 std::cerr << "Unit tests for basics of UFF bond-stretch terms." << std::endl;
81
82 ForceFields::UFF::AtomicParams p1, p2;
83 double restLen, forceConstant;
84
85 // sp3 carbon:
86 p1.r1 = .757;
87 p1.Z1 = 1.912;
88 p1.GMP_Xi = 5.343;
89
90 // sp3 - sp3: checks basics
91 restLen = ForceFields::UFF::Utils::calcBondRestLength(1.0, &p1, &p1);
92 TEST_ASSERT(RDKit::feq(restLen, 1.514));
93
94 forceConstant =
95 ForceFields::UFF::Utils::calcBondForceConstant(restLen, &p1, &p1);
96 TEST_ASSERT(RDKit::feq(forceConstant, 699.5918));
97
98 // sp2 carbon:
99 p2.r1 = .732;
100 p2.Z1 = 1.912;
101 p2.GMP_Xi = 5.343;
102 // sp2 - sp2: checks rBO
103 restLen = ForceFields::UFF::Utils::calcBondRestLength(2.0, &p2, &p2);
104 TEST_ASSERT(RDKit::feq(restLen, 1.32883, 1e-5));
105
106 forceConstant =
107 ForceFields::UFF::Utils::calcBondForceConstant(restLen, &p2, &p2);
108 TEST_ASSERT(RDKit::feq(forceConstant, 1034.69, 1e-2));
109
110 // sp3 nitrogen:
111 p2.r1 = .700;
112 p2.Z1 = 2.544;
113 p2.GMP_Xi = 6.899;
114
115 // Csp3 - Nsp3: checks rEN
116 restLen = ForceFields::UFF::Utils::calcBondRestLength(1.0, &p1, &p2);
117 TEST_ASSERT(RDKit::feq(restLen, 1.451071, 1e-5));
118
119 forceConstant =
120 ForceFields::UFF::Utils::calcBondForceConstant(restLen, &p1, &p2);
121 TEST_ASSERT(RDKit::feq(forceConstant, 1057.27, 1e-2));
122
123 // amide bond: check we can reproduce values from the UFF paper:
124 // C_R:
125 p1.r1 = .729;
126 p1.Z1 = 1.912;
127 p1.GMP_Xi = 5.343;
128 // N_R:
129 p2.r1 = .699;
130 p2.Z1 = 2.544;
131 p2.GMP_Xi = 6.899;
132
133 restLen = ForceFields::UFF::Utils::calcBondRestLength(
134 ForceFields::UFF::Params::amideBondOrder, &p1, &p2);
135 TEST_ASSERT(RDKit::feq(restLen, 1.357, 1e-3)); // NOTE: the paper has 1.366
136
137 forceConstant =
138 ForceFields::UFF::Utils::calcBondForceConstant(restLen, &p1, &p2);
139 TEST_ASSERT(RDKit::feq(forceConstant, 1293., 1)); // NOTE: the paper has 1293
140
141 std::cerr << " done" << std::endl;
142 }
143
testUFF2()144 void testUFF2() {
145 std::cerr << "-------------------------------------" << std::endl;
146 std::cerr << "Unit tests for UFF bond-stretch terms." << std::endl;
147
148 ForceFields::ForceField ff;
149 RDGeom::Point3D p1(0, 0, 0), p2(1.514, 0, 0);
150 RDGeom::PointPtrVect &ps = ff.positions();
151 ps.push_back(&p1);
152 ps.push_back(&p2);
153
154 ForceFields::UFF::AtomicParams param1;
155 // sp3 carbon:
156 param1.r1 = .757;
157 param1.Z1 = 1.912;
158 param1.GMP_Xi = 5.343;
159
160 // C_3 - C_3, r0=1.514, k01=699.5918
161 ForceFields::ForceFieldContrib *bs;
162 bs = new ForceFields::UFF::BondStretchContrib(&ff, 0, 1, 1, ¶m1, ¶m1);
163 ff.contribs().push_back(ForceFields::ContribPtr(bs));
164 ff.initialize();
165
166 double *p, *g;
167 p = new double[6];
168 g = new double[6];
169 for (int i = 0; i < 6; i++) {
170 p[i] = 0.0;
171 g[i] = 0.0;
172 }
173
174 double E;
175 // edge case: zero bond length:
176 E = bs->getEnergy(p);
177 TEST_ASSERT(E > 0.0);
178 bs->getGrad(p, g);
179 for (int i = 0; i < 6; i++) {
180 TEST_ASSERT(fabs(g[i]) > 0.0);
181 }
182
183 p[0] = 0;
184 p[3] = 1.514;
185 for (int i = 0; i < 6; i++) {
186 g[i] = 0.0;
187 }
188 ff.initialize();
189 E = bs->getEnergy(p);
190 TEST_ASSERT(RDKit::feq(E, 0.0));
191 bs->getGrad(p, g);
192 for (int i = 0; i < 6; i++) {
193 TEST_ASSERT(RDKit::feq(g[i], 0.0));
194 }
195
196 (*ff.positions()[1])[0] = 1.814;
197 p[3] = 1.814;
198 ff.initialize();
199 E = bs->getEnergy(p);
200 TEST_ASSERT(RDKit::feq(E, 31.4816));
201 bs->getGrad(p, g);
202 TEST_ASSERT(RDKit::feq(g[0], -209.8775));
203 TEST_ASSERT(RDKit::feq(g[3], 209.8775));
204 TEST_ASSERT(RDKit::feq(g[1], 0.0));
205 TEST_ASSERT(RDKit::feq(g[2], 0.0));
206 TEST_ASSERT(RDKit::feq(g[4], 0.0));
207 TEST_ASSERT(RDKit::feq(g[5], 0.0));
208
209 // try a different axis:
210 for (int i = 0; i < 6; i++) {
211 g[i] = 0.0;
212 p[i] = 0.0;
213 }
214 ff.initialize();
215 (*ff.positions()[1])[0] = 0.0;
216 (*ff.positions()[1])[2] = 1.814;
217 p[5] = 1.814;
218 E = bs->getEnergy(p);
219 TEST_ASSERT(RDKit::feq(E, 31.4816));
220 bs->getGrad(p, g);
221 TEST_ASSERT(RDKit::feq(g[2], -209.8775));
222 TEST_ASSERT(RDKit::feq(g[5], 209.8775));
223 TEST_ASSERT(RDKit::feq(g[0], 0.0));
224 TEST_ASSERT(RDKit::feq(g[1], 0.0));
225 TEST_ASSERT(RDKit::feq(g[3], 0.0));
226 TEST_ASSERT(RDKit::feq(g[4], 0.0));
227
228 // try a bit of minimization
229 RDGeom::Point3D d;
230 ff.initialize();
231 (*ff.positions()[1])[2] = 0.0;
232 (*ff.positions()[1])[0] = 1.814;
233 ff.minimize(10, 1e-8);
234 d = *(RDGeom::Point3D *)ff.positions()[0] -
235 *(RDGeom::Point3D *)ff.positions()[1];
236 TEST_ASSERT(RDKit::feq(d.length(), 1.514, 1e-3));
237
238 // minimize in "3D"
239 ff.initialize();
240 (*ff.positions()[1])[2] = 1.1;
241 (*ff.positions()[1])[1] = 0.9;
242 (*ff.positions()[1])[0] = 1.00;
243 ff.minimize(10, 1e-8);
244 d = *(RDGeom::Point3D *)ff.positions()[0] -
245 *(RDGeom::Point3D *)ff.positions()[1];
246 TEST_ASSERT(RDKit::feq(d.length(), 1.514, 1e-3));
247
248 delete[] p;
249 delete[] g;
250 std::cerr << " done" << std::endl;
251 }
252
testUFF3()253 void testUFF3() {
254 std::cerr << "-------------------------------------" << std::endl;
255 std::cerr << "Unit tests for basics of UFF angle terms." << std::endl;
256
257 ForceFields::UFF::AtomicParams p1, p2, p3;
258 double restLen, forceConstant;
259
260 // sp3 carbon:
261 p3.r1 = .757;
262 p3.Z1 = 1.912;
263 p3.GMP_Xi = 5.343;
264 p3.theta0 = 109.47 * M_PI / 180.0;
265
266 // sp3 - sp3: checks basics
267 restLen = ForceFields::UFF::Utils::calcBondRestLength(1.0, &p3, &p3);
268 TEST_ASSERT(RDKit::feq(restLen, 1.514));
269
270 // C_3 - C_3 - C_3
271 forceConstant = ForceFields::UFF::Utils::calcAngleForceConstant(
272 p3.theta0, 1, 1, &p3, &p3, &p3);
273 // TEST_ASSERT(RDKit::feq(forceConstant,699.5918));
274
275 // amide bond bend:
276 // C_R - N_R - C_3
277 // C_R:
278 p1.r1 = .729;
279 p1.Z1 = 1.912;
280 p1.GMP_Xi = 5.343;
281 // N_R:
282 p2.r1 = .699;
283 p2.Z1 = 2.544;
284 p2.GMP_Xi = 6.899;
285 p2.theta0 = 120.0 * M_PI / 180.;
286 restLen = ForceFields::UFF::Utils::calcBondRestLength(
287 ForceFields::UFF::Params::amideBondOrder, &p1, &p2);
288 TEST_ASSERT(RDKit::feq(restLen, 1.357, 1e-3));
289 restLen = ForceFields::UFF::Utils::calcBondRestLength(1.0, &p2, &p3);
290 TEST_ASSERT(RDKit::feq(restLen, 1.450, 1e-3));
291
292 forceConstant = ForceFields::UFF::Utils::calcAngleForceConstant(
293 p2.theta0, ForceFields::UFF::Params::amideBondOrder, 1, &p1, &p2, &p3);
294 TEST_ASSERT(RDKit::feq(forceConstant, 211.0, 1e-1)); // paper has 105.5
295
296 std::cerr << " done" << std::endl;
297 }
298
testUFF4()299 void testUFF4() {
300 std::cerr << "-------------------------------------" << std::endl;
301 std::cerr << "Unit tests for UFF angle-bend terms." << std::endl;
302
303 ForceFields::ForceField ff;
304 RDGeom::Point3D p1(1.514, 0, 0), p2(0, 0, 0), p3(0.1, 1.5, 0);
305 RDGeom::PointPtrVect &ps = ff.positions();
306 ps.push_back(&p1);
307 ps.push_back(&p2);
308 ps.push_back(&p3);
309
310 ForceFields::UFF::AtomicParams param1;
311 // sp3 carbon:
312 param1.r1 = .757;
313 param1.Z1 = 1.912;
314 param1.GMP_Xi = 5.343;
315 // cheat to get the angle to 90 so that testing is easier:
316 param1.theta0 = 90.0 * M_PI / 180.;
317
318 // C_3 - C_3, r0=1.514, k01=699.5918
319 ForceFields::ForceFieldContrib *contrib;
320 contrib =
321 new ForceFields::UFF::BondStretchContrib(&ff, 0, 1, 1, ¶m1, ¶m1);
322 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
323 contrib =
324 new ForceFields::UFF::BondStretchContrib(&ff, 1, 2, 1, ¶m1, ¶m1);
325 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
326 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 1, 1, ¶m1,
327 ¶m1, ¶m1);
328 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
329
330 RDGeom::Point3D d, v1, v2;
331 double theta;
332 #if 1
333 // ------- ------- ------- ------- ------- ------- -------
334 // try a bit of minimization
335 ff.initialize();
336 ff.minimize(10, 1e-8, 1e-8);
337
338 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
339 *(RDGeom::Point3D *)ff.positions()[1];
340 v2 = *(RDGeom::Point3D *)ff.positions()[1] -
341 *(RDGeom::Point3D *)ff.positions()[2];
342 theta = v1.angleTo(v2);
343
344 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
345 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
346 TEST_ASSERT(RDKit::feq(theta, 90 * M_PI / 180., 1e-4));
347
348 // ------- ------- ------- ------- ------- ------- -------
349 // more complicated atomic coords:
350 p1.x = 1.3;
351 p1.y = 0.1;
352 p1.z = 0.1;
353 p2.x = -0.1;
354 p2.y = 0.05;
355 p2.z = -0.05;
356 p3.x = 0.1;
357 p3.y = 1.5;
358 p3.z = 0.05;
359 ff.initialize();
360 ff.minimize(10, 1e-8, 1e-8);
361
362 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
363 *(RDGeom::Point3D *)ff.positions()[1];
364 v2 = *(RDGeom::Point3D *)ff.positions()[1] -
365 *(RDGeom::Point3D *)ff.positions()[2];
366 theta = v1.angleTo(v2);
367
368 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
369 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
370 TEST_ASSERT(RDKit::feq(theta, 90 * M_PI / 180., 1e-4));
371
372 // ------- ------- ------- ------- ------- ------- -------
373 // try for the tetrahedral angle instead of 90:
374 param1.theta0 = 109.47 * M_PI / 180.;
375 ff.contribs().pop_back();
376 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 1, 1, ¶m1,
377 ¶m1, ¶m1);
378 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
379
380 p1.x = 1.3;
381 p1.y = 0.1;
382 p1.z = 0.1;
383 p2.x = -0.1;
384 p2.y = 0.05;
385 p2.z = -0.05;
386 p3.x = 0.1;
387 p3.y = 1.5;
388 p3.z = 0.05;
389
390 ff.initialize();
391 ff.minimize(100, 1e-8, 1e-8);
392 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
393 *(RDGeom::Point3D *)ff.positions()[1];
394 v2 = *(RDGeom::Point3D *)ff.positions()[2] -
395 *(RDGeom::Point3D *)ff.positions()[1];
396 theta = v1.angleTo(v2);
397 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
398 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
399 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
400
401 #endif
402
403 // ------- ------- ------- ------- ------- ------- -------
404 //
405 // Do a series of "special cases" (i.e. test the functional forms
406 // for linear, trigonal planar, square planar and octahedral)
407 //
408 // ------- ------- ------- ------- ------- ------- -------
409
410 // ------- ------- ------- ------- ------- ------- -------
411 // test a linear molecule:
412 param1.theta0 = M_PI;
413 // ff.contribs().pop_back();
414 // ff.contribs().pop_back();
415 ff.contribs().pop_back();
416 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 1, 1, ¶m1,
417 ¶m1, ¶m1, 2);
418 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
419
420 p1.x = 1.3;
421 p1.y = 0.1;
422 p1.z = 0.0;
423 p2.x = 0.0;
424 p2.y = 0.0;
425 p2.z = 0.0;
426 p3.x = -1.3;
427 p3.y = 0.1;
428 p3.z = 0.00;
429 ff.initialize();
430 ff.minimize(100, 1e-8, 1e-8);
431
432 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
433 *(RDGeom::Point3D *)ff.positions()[1];
434 v2 = *(RDGeom::Point3D *)ff.positions()[2] -
435 *(RDGeom::Point3D *)ff.positions()[1];
436 theta = v1.angleTo(v2);
437
438 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
439 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
440 std::cerr << "theta = " << theta << "; theta0 = " << param1.theta0
441 << std::endl;
442 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
443
444 // ------- ------- ------- ------- ------- ------- -------
445 // test n=3:
446 param1.theta0 = 120. * M_PI / 180.0;
447 // ff.contribs().pop_back();
448 // ff.contribs().pop_back();
449 ff.contribs().pop_back();
450 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 1, 1, ¶m1,
451 ¶m1, ¶m1, 3);
452 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
453
454 p1.x = 1.3;
455 p1.y = 0.1;
456 p1.z = 0.0;
457 p2.x = 0.0;
458 p2.y = 0.0;
459 p2.z = 0.0;
460 p3.x = -.3;
461 p3.y = -1.3;
462 p3.z = 0.00;
463
464 ff.initialize();
465 ff.minimize(100, 1e-8, 1e-8);
466
467 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
468 *(RDGeom::Point3D *)ff.positions()[1];
469 v2 = *(RDGeom::Point3D *)ff.positions()[2] -
470 *(RDGeom::Point3D *)ff.positions()[1];
471 theta = v1.angleTo(v2);
472
473 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
474 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
475 std::cerr << "theta = " << std::fixed << std::setprecision(6) << theta
476 << ", param1.theta0 = " << std::fixed << std::setprecision(6)
477 << param1.theta0 << std::endl;
478 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
479
480 // ------- ------- ------- ------- ------- ------- -------
481 // test n=4:
482 param1.theta0 = M_PI / 2.0;
483 // ff.contribs().pop_back();
484 // ff.contribs().pop_back();
485 ff.contribs().pop_back();
486 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 1, 1, ¶m1,
487 ¶m1, ¶m1, 4);
488 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
489
490 p1.x = 1.3;
491 p1.y = 0.1;
492 p1.z = 0.0;
493 p2.x = 0.0;
494 p2.y = 0.0;
495 p2.z = 0.0;
496 p3.x = -.3;
497 p3.y = -1.3;
498 p3.z = 0.00;
499
500 ff.initialize();
501 ff.minimize(100, 1e-8, 1e-8);
502
503 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
504 *(RDGeom::Point3D *)ff.positions()[1];
505 v2 = *(RDGeom::Point3D *)ff.positions()[2] -
506 *(RDGeom::Point3D *)ff.positions()[1];
507 theta = v1.angleTo(v2);
508
509 TEST_ASSERT(RDKit::feq(v1.length(), 1.514, 1e-3));
510 TEST_ASSERT(RDKit::feq(v2.length(), 1.514, 1e-3));
511 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
512
513 #if 0
514 std::cerr << " " << *ff.positions()[0] << std::endl;
515 std::cerr << " " << *ff.positions()[1] << std::endl;
516 std::cerr << " " << *ff.positions()[2] << std::endl;
517
518 std::cerr << "v1: " << v1 << std::endl;
519 std::cerr << "v2: " << v2 << std::endl;
520 std::cerr << "FINAL: " << v1.angleTo(v2) << " " << v1.signedAngleTo(v2) << std::endl;
521 #endif
522
523 std::cerr << " done" << std::endl;
524 }
525
testUFF5()526 void testUFF5() {
527 std::cerr << "-------------------------------------" << std::endl;
528 std::cerr << " Test Simple UFF molecule optimizations." << std::endl;
529
530 ForceFields::ForceField ff;
531 RDGeom::Point3D p1, p2, p3, p4, p5, p6;
532 RDGeom::PointPtrVect &ps = ff.positions();
533 ps.push_back(&p1);
534 ps.push_back(&p2);
535 ps.push_back(&p3);
536 ps.push_back(&p4);
537 ps.push_back(&p5);
538 ps.push_back(&p6);
539
540 ForceFields::UFF::AtomicParams param1, param2;
541 // sp2 carbon:
542 param1.r1 = .732;
543 param1.Z1 = 1.912;
544 param1.GMP_Xi = 5.343;
545 param1.theta0 = 120. * M_PI / 180.;
546
547 // H_1:
548 param2.r1 = 0.354;
549 param2.Z1 = 0.712;
550 param2.GMP_Xi = 4.528;
551
552 ForceFields::ForceFieldContrib *contrib;
553
554 // build ethylene:
555 contrib =
556 new ForceFields::UFF::BondStretchContrib(&ff, 0, 1, 2, ¶m1, ¶m1);
557 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
558 contrib =
559 new ForceFields::UFF::BondStretchContrib(&ff, 0, 2, 1, ¶m1, ¶m2);
560 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
561 contrib =
562 new ForceFields::UFF::BondStretchContrib(&ff, 0, 3, 1, ¶m1, ¶m2);
563 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
564 contrib =
565 new ForceFields::UFF::BondStretchContrib(&ff, 1, 4, 1, ¶m1, ¶m2);
566 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
567 contrib =
568 new ForceFields::UFF::BondStretchContrib(&ff, 1, 5, 1, ¶m1, ¶m2);
569 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
570
571 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 0, 2, 2, 1, ¶m1,
572 ¶m1, ¶m2, 3);
573 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
574 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 0, 3, 2, 1, ¶m1,
575 ¶m1, ¶m2, 3);
576 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
577 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 2, 0, 3, 1, 1, ¶m2,
578 ¶m1, ¶m2, 3);
579 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
580
581 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 4, 2, 1, ¶m1,
582 ¶m1, ¶m2, 3);
583 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
584 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 5, 2, 1, ¶m1,
585 ¶m1, ¶m2, 3);
586 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
587 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 4, 1, 5, 1, 1, ¶m2,
588 ¶m1, ¶m2, 3);
589 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
590
591 // dodge the fact that we're not using torsions yet by putting
592 // everything in the z=0 plane:
593 p1.x = -0.58;
594 p1.y = -0.33;
595 p1.z = 0.0;
596
597 p2.x = 0.58;
598 p2.y = 0.33;
599 p2.z = 0.0;
600
601 p3.x = -0.61;
602 p3.y = -1.43;
603 p3.z = 0.0;
604
605 p4.x = -1.54;
606 p4.y = 0.20;
607 p4.z = 0.0;
608
609 p5.x = 0.61;
610 p5.y = 1.43;
611 p5.z = 0.0;
612
613 p6.x = 1.54;
614 p6.y = -0.20;
615 p6.z = 0.0;
616
617 RDGeom::Point3D d, v1, v2;
618 double theta;
619 // ------- ------- ------- ------- ------- ------- -------
620 // try a bit of minimization
621 ff.initialize();
622 ff.minimize(10, 1e-8, 1e-8);
623
624 double CCDblBondLen =
625 ForceFields::UFF::Utils::calcBondRestLength(2, ¶m1, ¶m1);
626 double CHBondLen =
627 ForceFields::UFF::Utils::calcBondRestLength(1, ¶m1, ¶m2);
628
629 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
630 *(RDGeom::Point3D *)ff.positions()[1];
631 v2 = *(RDGeom::Point3D *)ff.positions()[0] -
632 *(RDGeom::Point3D *)ff.positions()[2];
633 theta = v1.angleTo(v2);
634 TEST_ASSERT(RDKit::feq(v1.length(), CCDblBondLen, 1e-3));
635 TEST_ASSERT(RDKit::feq(v2.length(), CHBondLen, 1e-3));
636 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
637 v2 = *(RDGeom::Point3D *)ff.positions()[0] -
638 *(RDGeom::Point3D *)ff.positions()[3];
639 theta = v1.angleTo(v2);
640 TEST_ASSERT(RDKit::feq(v2.length(), CHBondLen, 1e-3));
641 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
642
643 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
644 *(RDGeom::Point3D *)ff.positions()[2];
645 theta = v1.angleTo(v2);
646 TEST_ASSERT(RDKit::feq(theta, param1.theta0, 1e-4));
647
648 std::cerr << " done" << std::endl;
649 }
650
testUFF6()651 void testUFF6() {
652 std::cerr << "-------------------------------------" << std::endl;
653 std::cerr << "Unit tests for UFF nonbonded terms." << std::endl;
654
655 ForceFields::ForceField ff;
656 RDGeom::Point3D p1(0, 0, 0), p2(0.0, 0, 0);
657 RDGeom::PointPtrVect &ps = ff.positions();
658 ps.push_back(&p1);
659 ps.push_back(&p2);
660
661 ForceFields::UFF::AtomicParams param1;
662 // sp3 carbon:
663 param1.r1 = .757;
664 param1.Z1 = 1.912;
665 param1.GMP_Xi = 5.343;
666 param1.x1 = 3.851;
667 param1.D1 = 0.105;
668
669 ff.initialize();
670 ForceFields::ForceFieldContrib *contrib;
671 contrib = new ForceFields::UFF::vdWContrib(&ff, 0, 1, ¶m1, ¶m1);
672 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
673
674 // try a bit of minimization
675 RDGeom::Point3D d;
676 ff.initialize();
677
678 // edge case: our energy at zero length should be zero:
679 double E;
680 E = ff.calcEnergy();
681 TEST_ASSERT(RDKit::feq(E, 0.0));
682
683 (*ff.positions()[0])[0] = 0.0;
684 (*ff.positions()[1])[0] = 4.0;
685 ff.minimize(10, 1e-8, 1e-8);
686 d = *(RDGeom::Point3D *)ff.positions()[0] -
687 *(RDGeom::Point3D *)ff.positions()[1];
688 TEST_ASSERT(RDKit::feq(d.length(), 3.851, 1e-3));
689
690 // minimize in "3D"
691 ff.initialize();
692 (*ff.positions()[0])[0] = 0.0;
693 (*ff.positions()[0])[1] = 0.0;
694 (*ff.positions()[0])[2] = 0.0;
695 (*ff.positions()[1])[2] = 3.1;
696 (*ff.positions()[1])[1] = 0.9;
697 (*ff.positions()[1])[0] = 1.00;
698 ff.minimize(10, 1e-8, 1e-8);
699 d = *(RDGeom::Point3D *)ff.positions()[0] -
700 *(RDGeom::Point3D *)ff.positions()[1];
701 TEST_ASSERT(RDKit::feq(d.length(), 3.851, 1e-3));
702
703 std::cerr << " done" << std::endl;
704 }
705
testUFF7()706 void testUFF7() {
707 std::cerr << "-------------------------------------" << std::endl;
708 std::cerr << " Test UFF torsional terms." << std::endl;
709
710 ForceFields::ForceField ff;
711 RDGeom::Point3D p1, p2, p3, p4;
712 RDGeom::PointPtrVect &ps = ff.positions();
713 ps.push_back(&p1);
714 ps.push_back(&p2);
715 ps.push_back(&p3);
716 ps.push_back(&p4);
717
718 ForceFields::UFF::AtomicParams param1;
719 // sp3 carbon:
720 param1.r1 = .757;
721 param1.Z1 = 1.912;
722 param1.GMP_Xi = 5.343;
723 param1.x1 = 3.851;
724 param1.D1 = 0.105;
725 param1.V1 = 2.119;
726 param1.U1 = 2.0;
727
728 RDGeom::Point3D d, v1, v2;
729 double cosPhi;
730
731 ForceFields::ForceFieldContrib *contrib;
732 // ------- ------- ------- ------- ------- ------- -------
733 // Basic SP3 - SP3
734 // ------- ------- ------- ------- ------- ------- -------
735 contrib = new ForceFields::UFF::TorsionAngleContrib(
736 &ff, 0, 1, 2, 3, 1, 6, 6, RDKit::Atom::SP3, RDKit::Atom::SP3, ¶m1,
737 ¶m1);
738 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
739 #if 1
740 p1.x = 0;
741 p1.y = 1.5;
742 p1.z = 0;
743
744 p2.x = 0.0;
745 p2.y = 0.0;
746 p2.z = 0.0;
747
748 p3.x = 1.5;
749 p3.y = 0.0;
750 p3.z = 0.0;
751
752 p4.x = 1.5;
753 p4.y = 0.0;
754 p4.z = 1.5;
755
756 ff.initialize();
757 ff.minimize(10, 1e-8, 1e-8);
758 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
759 *(RDGeom::Point3D *)ff.positions()[0],
760 *(RDGeom::Point3D *)ff.positions()[1],
761 *(RDGeom::Point3D *)ff.positions()[2],
762 *(RDGeom::Point3D *)ff.positions()[3]);
763 TEST_ASSERT(RDKit::feq(cosPhi, 0.5, 1e-4));
764
765 // ------- ------- ------- ------- ------- ------- -------
766 // Basic SP2 - SP2
767 // ------- ------- ------- ------- ------- ------- -------
768 ff.contribs().pop_back();
769 contrib = new ForceFields::UFF::TorsionAngleContrib(
770 &ff, 0, 1, 2, 3, 1, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP2, ¶m1,
771 ¶m1);
772 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
773 p1.x = 0;
774 p1.y = 1.5;
775 p1.z = 0.1;
776
777 p2.x = 0.0;
778 p2.y = 0.0;
779 p2.z = 0.0;
780
781 p3.x = 1.5;
782 p3.y = 0.0;
783 p3.z = 0.0;
784
785 p4.x = 1.5;
786 p4.y = 0.2;
787 p4.z = 1.5;
788
789 ff.initialize();
790 ff.minimize(10, 1e-8, 1e-8);
791 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
792 *(RDGeom::Point3D *)ff.positions()[0],
793 *(RDGeom::Point3D *)ff.positions()[1],
794 *(RDGeom::Point3D *)ff.positions()[2],
795 *(RDGeom::Point3D *)ff.positions()[3]);
796 TEST_ASSERT(RDKit::feq(cosPhi, 1.0, 1e-4));
797
798 // ------- ------- ------- ------- ------- ------- -------
799 // Basic SP2 - SP3
800 // ------- ------- ------- ------- ------- ------- -------
801 ff.contribs().pop_back();
802 contrib = new ForceFields::UFF::TorsionAngleContrib(
803 &ff, 0, 1, 2, 3, 1, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, ¶m1,
804 ¶m1);
805 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
806 p1.x = 0;
807 p1.y = 1.5;
808 p1.z = 0.1;
809
810 p2.x = 0.0;
811 p2.y = 0.0;
812 p2.z = 0.0;
813
814 p3.x = 1.5;
815 p3.y = 0.0;
816 p3.z = 0.0;
817
818 p4.x = 1.5;
819 p4.y = 0.2;
820 p4.z = 1.5;
821
822 ff.initialize();
823 ff.minimize(100, 1e-8, 1e-8);
824 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
825 *(RDGeom::Point3D *)ff.positions()[0],
826 *(RDGeom::Point3D *)ff.positions()[1],
827 *(RDGeom::Point3D *)ff.positions()[2],
828 *(RDGeom::Point3D *)ff.positions()[3]);
829 TEST_ASSERT(RDKit::feq(cosPhi, 0.5, 1e-4));
830
831 // ------- ------- ------- ------- ------- ------- -------
832 // special case for group 6 - group 6 bonds:
833 // ------- ------- ------- ------- ------- ------- -------
834 ff.contribs().pop_back();
835 contrib = new ForceFields::UFF::TorsionAngleContrib(
836 &ff, 0, 1, 2, 3, 1, 8, 8, RDKit::Atom::SP3, RDKit::Atom::SP3, ¶m1,
837 ¶m1);
838 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
839 p1.x = 0;
840 p1.y = 1.5;
841 p1.z = 0.1;
842
843 p2.x = 0.0;
844 p2.y = 0.0;
845 p2.z = 0.0;
846
847 p3.x = 1.5;
848 p3.y = 0.0;
849 p3.z = 0.0;
850
851 p4.x = 1.5;
852 p4.y = 0.2;
853 p4.z = 1.5;
854
855 ff.initialize();
856 ff.minimize(100, 1e-8, 1e-8);
857 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
858 *(RDGeom::Point3D *)ff.positions()[0],
859 *(RDGeom::Point3D *)ff.positions()[1],
860 *(RDGeom::Point3D *)ff.positions()[2],
861 *(RDGeom::Point3D *)ff.positions()[3]);
862 TEST_ASSERT(RDKit::feq(cosPhi, 0.0, 1e-4));
863
864 // ------- ------- ------- ------- ------- ------- -------
865 // special case for SP3 group 6 - SP2 other group
866 // ------- ------- ------- ------- ------- ------- -------
867 ff.contribs().pop_back();
868 contrib = new ForceFields::UFF::TorsionAngleContrib(
869 &ff, 0, 1, 2, 3, 1, 8, 6, RDKit::Atom::SP3, RDKit::Atom::SP2, ¶m1,
870 ¶m1);
871 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
872 p1.x = 0;
873 p1.y = 1.5;
874 p1.z = 0.1;
875
876 p2.x = 0.0;
877 p2.y = 0.0;
878 p2.z = 0.0;
879
880 p3.x = 1.5;
881 p3.y = 0.0;
882 p3.z = 0.0;
883
884 p4.x = 1.5;
885 p4.y = 0.2;
886 p4.z = 1.5;
887
888 ff.initialize();
889 ff.minimize(100, 1e-8, 1e-8);
890 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
891 *(RDGeom::Point3D *)ff.positions()[0],
892 *(RDGeom::Point3D *)ff.positions()[1],
893 *(RDGeom::Point3D *)ff.positions()[2],
894 *(RDGeom::Point3D *)ff.positions()[3]);
895 TEST_ASSERT(RDKit::feq(cosPhi, 0.0, 1e-4));
896
897 #endif
898
899 // ------- ------- ------- ------- ------- ------- -------
900 // special case for (SP2 -) SP2 - SP3
901 // ------- ------- ------- ------- ------- ------- -------
902 ff.contribs().pop_back();
903 contrib = new ForceFields::UFF::TorsionAngleContrib(
904 &ff, 0, 1, 2, 3, 1, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, ¶m1,
905 ¶m1, true);
906 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
907 p1.x = 0;
908 p1.y = 1.5;
909 p1.z = 0.1;
910
911 p2.x = 0.0;
912 p2.y = 0.0;
913 p2.z = 0.0;
914
915 p3.x = 1.5;
916 p3.y = 0.0;
917 p3.z = 0.0;
918
919 p4.x = 1.5;
920 p4.y = 0.2;
921 p4.z = 1.5;
922
923 ff.initialize();
924 ff.minimize(100, 1e-8, 1e-8);
925 cosPhi = ForceFields::UFF::Utils::calculateCosTorsion(
926 *(RDGeom::Point3D *)ff.positions()[0],
927 *(RDGeom::Point3D *)ff.positions()[1],
928 *(RDGeom::Point3D *)ff.positions()[2],
929 *(RDGeom::Point3D *)ff.positions()[3]);
930 TEST_ASSERT(RDKit::feq(cosPhi, 0.5, 1e-4));
931
932 std::cerr << " done" << std::endl;
933 }
934
testUFFParams()935 void testUFFParams() {
936 std::cerr << "-------------------------------------" << std::endl;
937 std::cerr << " Test UFF Parameter objects" << std::endl;
938
939 ForceFields::UFF::ParamCollection *params =
940 ForceFields::UFF::ParamCollection::getParams();
941 TEST_ASSERT(params);
942
943 const ForceFields::UFF::AtomicParams *ptr;
944 ptr = (*params)("C_3");
945 TEST_ASSERT(ptr);
946 TEST_ASSERT(RDKit::feq(ptr->r1, 0.757));
947 TEST_ASSERT(RDKit::feq(ptr->theta0, 109.47 * M_PI / 180.));
948 TEST_ASSERT(RDKit::feq(ptr->x1, 3.851));
949 TEST_ASSERT(RDKit::feq(ptr->D1, 0.105));
950 TEST_ASSERT(RDKit::feq(ptr->zeta, 12.73));
951 TEST_ASSERT(RDKit::feq(ptr->Z1, 1.912));
952 TEST_ASSERT(RDKit::feq(ptr->V1, 2.119));
953 TEST_ASSERT(RDKit::feq(ptr->GMP_Xi, 5.343));
954 TEST_ASSERT(RDKit::feq(ptr->GMP_Hardness, 5.063));
955 TEST_ASSERT(RDKit::feq(ptr->GMP_Radius, 0.759));
956
957 ptr = (*params)("N_3");
958 TEST_ASSERT(ptr);
959
960 ptr = (*params)("C_5");
961 TEST_ASSERT(!ptr);
962 }
963
testUFF8()964 void testUFF8() {
965 std::cerr << "-------------------------------------" << std::endl;
966 std::cerr << " Test Simple UFF molecule optimization, part 2." << std::endl;
967
968 ForceFields::ForceField ff;
969 RDGeom::Point3D p1, p2, p3, p4, p5, p6;
970 RDGeom::PointPtrVect &ps = ff.positions();
971 ps.push_back(&p1);
972 ps.push_back(&p2);
973 ps.push_back(&p3);
974 ps.push_back(&p4);
975 ps.push_back(&p5);
976 ps.push_back(&p6);
977
978 ForceFields::UFF::ParamCollection *params =
979 ForceFields::UFF::ParamCollection::getParams();
980 const ForceFields::UFF::AtomicParams *param1, *param2;
981
982 // C_2 (sp2 carbon):
983 param1 = (*params)("C_2");
984 TEST_ASSERT(param1);
985 // H_:
986 param2 = (*params)("H_");
987 TEST_ASSERT(param2);
988 ForceFields::ForceFieldContrib *contrib;
989
990 // build ethylene:
991 // BONDS:
992 contrib =
993 new ForceFields::UFF::BondStretchContrib(&ff, 0, 1, 2, param1, param1);
994 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
995 contrib =
996 new ForceFields::UFF::BondStretchContrib(&ff, 0, 2, 1, param1, param2);
997 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
998 contrib =
999 new ForceFields::UFF::BondStretchContrib(&ff, 0, 3, 1, param1, param2);
1000 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1001 contrib =
1002 new ForceFields::UFF::BondStretchContrib(&ff, 1, 4, 1, param1, param2);
1003 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1004 contrib =
1005 new ForceFields::UFF::BondStretchContrib(&ff, 1, 5, 1, param1, param2);
1006 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1007
1008 // ANGLES:
1009 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 0, 2, 2, 1, param1,
1010 param1, param2, 3);
1011 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1012 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 0, 3, 2, 1, param1,
1013 param1, param2, 3);
1014 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1015 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 2, 0, 3, 1, 1, param2,
1016 param1, param2, 3);
1017 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1018 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 4, 2, 1, param1,
1019 param1, param2, 3);
1020 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1021 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 5, 2, 1, param1,
1022 param1, param2, 3);
1023 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1024 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 4, 1, 5, 1, 1, param2,
1025 param1, param2, 3);
1026 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1027
1028 // DIHEDRALS:
1029 contrib = new ForceFields::UFF::TorsionAngleContrib(
1030 &ff, 2, 0, 1, 4, 2, 6, 6, RDKit::Atom::SP3, RDKit::Atom::SP3, param1,
1031 param1);
1032 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1033 contrib = new ForceFields::UFF::TorsionAngleContrib(
1034 &ff, 2, 0, 1, 5, 2, 6, 6, RDKit::Atom::SP3, RDKit::Atom::SP3, param1,
1035 param1);
1036 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1037 contrib = new ForceFields::UFF::TorsionAngleContrib(
1038 &ff, 3, 0, 1, 4, 2, 6, 6, RDKit::Atom::SP3, RDKit::Atom::SP3, param1,
1039 param1);
1040 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1041 contrib = new ForceFields::UFF::TorsionAngleContrib(
1042 &ff, 3, 0, 1, 5, 2, 6, 6, RDKit::Atom::SP3, RDKit::Atom::SP3, param1,
1043 param1);
1044 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1045
1046 p1.x = -0.58;
1047 p1.y = -0.33;
1048 p1.z = 0.1;
1049
1050 p2.x = 0.58;
1051 p2.y = 0.33;
1052 p2.z = 0.1;
1053
1054 p3.x = -0.61;
1055 p3.y = -1.43;
1056 p3.z = 0.0;
1057
1058 p4.x = -1.54;
1059 p4.y = 0.20;
1060 p4.z = 0.0;
1061
1062 p5.x = 0.61;
1063 p5.y = 1.43;
1064 p5.z = 0.0;
1065
1066 p6.x = 1.54;
1067 p6.y = -0.20;
1068 p6.z = 0.0;
1069
1070 RDGeom::Point3D d, v1, v2;
1071 double theta;
1072 // ------- ------- ------- ------- ------- ------- -------
1073 // try a bit of minimization
1074 ff.initialize();
1075 ff.minimize(100, 1e-8, 1e-8);
1076
1077 double CCDblBondLen =
1078 ForceFields::UFF::Utils::calcBondRestLength(2, param1, param1);
1079 double CHBondLen =
1080 ForceFields::UFF::Utils::calcBondRestLength(1, param1, param2);
1081
1082 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
1083 *(RDGeom::Point3D *)ff.positions()[1];
1084 v2 = *(RDGeom::Point3D *)ff.positions()[0] -
1085 *(RDGeom::Point3D *)ff.positions()[2];
1086 theta = v1.angleTo(v2);
1087 TEST_ASSERT(RDKit::feq(v1.length(), CCDblBondLen, 1e-3));
1088 TEST_ASSERT(RDKit::feq(v2.length(), CHBondLen, 1e-3));
1089 TEST_ASSERT(RDKit::feq(theta, param1->theta0, 1e-4));
1090 v2 = *(RDGeom::Point3D *)ff.positions()[0] -
1091 *(RDGeom::Point3D *)ff.positions()[3];
1092 theta = v1.angleTo(v2);
1093 TEST_ASSERT(RDKit::feq(v2.length(), CHBondLen, 1e-3));
1094 TEST_ASSERT(RDKit::feq(theta, param1->theta0, 1e-4));
1095
1096 v1 = *(RDGeom::Point3D *)ff.positions()[0] -
1097 *(RDGeom::Point3D *)ff.positions()[2];
1098 theta = v1.angleTo(v2);
1099 TEST_ASSERT(RDKit::feq(theta, param1->theta0, 1e-4));
1100
1101 std::cerr << " done" << std::endl;
1102 }
1103
testUFFTorsionConflict()1104 void testUFFTorsionConflict() {
1105 std::cerr << "-------------------------------------" << std::endl;
1106 std::cerr << " Test UFF Torsion Conflicts." << std::endl;
1107
1108 ForceFields::ForceField ff;
1109 RDGeom::Point3D p1, p2, p3, p4, p5, p6, p7;
1110 RDGeom::PointPtrVect &ps = ff.positions();
1111 ps.push_back(&p1);
1112 ps.push_back(&p2);
1113 ps.push_back(&p3);
1114 ps.push_back(&p4);
1115 ps.push_back(&p5);
1116 ps.push_back(&p6);
1117 ps.push_back(&p7);
1118
1119 ForceFields::UFF::ParamCollection *params =
1120 ForceFields::UFF::ParamCollection::getParams();
1121 const ForceFields::UFF::AtomicParams *param1, *param2, *param3;
1122
1123 // C_2 (sp2 carbon):
1124 param1 = (*params)("C_2");
1125 TEST_ASSERT(param1);
1126 // H_:
1127 param2 = (*params)("H_");
1128 TEST_ASSERT(param2);
1129 // C_3 (sp3 carbon):
1130 param3 = (*params)("C_3");
1131 TEST_ASSERT(param3);
1132
1133 ForceFields::ForceFieldContrib *contrib;
1134
1135 // BONDS:
1136 contrib =
1137 new ForceFields::UFF::BondStretchContrib(&ff, 0, 1, 2, param1, param1);
1138 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1139 contrib =
1140 new ForceFields::UFF::BondStretchContrib(&ff, 1, 2, 1, param1, param3);
1141 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1142 contrib =
1143 new ForceFields::UFF::BondStretchContrib(&ff, 1, 3, 1, param1, param2);
1144 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1145 contrib =
1146 new ForceFields::UFF::BondStretchContrib(&ff, 2, 4, 1, param1, param2);
1147 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1148 contrib =
1149 new ForceFields::UFF::BondStretchContrib(&ff, 2, 5, 1, param1, param2);
1150 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1151 contrib =
1152 new ForceFields::UFF::BondStretchContrib(&ff, 2, 6, 1, param1, param2);
1153 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1154
1155 #if 1
1156 // ANGLES:
1157 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 2, 2.0, 1.0,
1158 param1, param1, param3);
1159 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1160 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 0, 1, 3, 2.0, 1.0,
1161 param1, param1, param2);
1162 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1163 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 2, 4, 1.0, 1.0,
1164 param1, param3, param2);
1165 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1166 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 2, 5, 1.0, 1.0,
1167 param1, param3, param2);
1168 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1169 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 1, 2, 6, 1.0, 1.0,
1170 param1, param3, param2);
1171 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1172 contrib = new ForceFields::UFF::AngleBendContrib(&ff, 2, 1, 3, 1.0, 1.0,
1173 param3, param1, param2);
1174 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1175 #endif
1176
1177 // DIHEDRALS:
1178 contrib = new ForceFields::UFF::TorsionAngleContrib(
1179 &ff, 0, 1, 2, 4, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1180 param3, true);
1181 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1182 contrib = new ForceFields::UFF::TorsionAngleContrib(
1183 &ff, 3, 1, 2, 4, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1184 param3, false);
1185 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1186
1187 contrib = new ForceFields::UFF::TorsionAngleContrib(
1188 &ff, 0, 1, 2, 5, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1189 param3, true);
1190 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1191 contrib = new ForceFields::UFF::TorsionAngleContrib(
1192 &ff, 3, 1, 2, 5, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1193 param3, false);
1194 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1195
1196 contrib = new ForceFields::UFF::TorsionAngleContrib(
1197 &ff, 0, 1, 2, 6, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1198 param3, true);
1199 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1200 contrib = new ForceFields::UFF::TorsionAngleContrib(
1201 &ff, 3, 1, 2, 6, 1.0, 6, 6, RDKit::Atom::SP2, RDKit::Atom::SP3, param1,
1202 param3, false);
1203 ff.contribs().push_back(ForceFields::ContribPtr(contrib));
1204
1205 p1.x = 0.5411;
1206 p1.y = -0.7741;
1207 p1.z = 0.0902;
1208
1209 p2.x = -0.5622;
1210 p2.y = -0.0368;
1211 p2.z = 0.1202;
1212
1213 p3.x = -0.5101;
1214 p3.y = 1.4485;
1215 p3.z = 0.0816;
1216
1217 p4.x = -1.5285;
1218 p4.y = -0.5341;
1219 p4.z = 0.1892;
1220
1221 p5.x = 0.5097;
1222 p5.y = 1.8065;
1223 p5.z = 0.1988;
1224
1225 p6.x = -1.1436;
1226 p6.y = 1.8781;
1227 p6.z = 0.8983;
1228
1229 p7.x = -0.9145;
1230 p7.y = 1.8185;
1231 p7.z = -0.8983;
1232
1233 RDGeom::Point3D d, v1, v2;
1234 // ------- ------- ------- ------- ------- ------- -------
1235 // try a bit of minimization
1236 ff.initialize();
1237 ff.minimize(100, 1e-8, 1e-8);
1238
1239 #if 1
1240 std::cerr.setf(std::ios_base::fixed, std::ios_base::floatfield);
1241 std::cerr.precision(4);
1242 std::cerr << "C " << *ff.positions()[0] << std::endl;
1243 std::cerr << "C " << *ff.positions()[1] << std::endl;
1244 std::cerr << "C " << *ff.positions()[2] << std::endl;
1245 std::cerr << "H " << *ff.positions()[3] << std::endl;
1246 std::cerr << "H " << *ff.positions()[4] << std::endl;
1247 std::cerr << "O " << *ff.positions()[5] << std::endl;
1248 std::cerr << "F " << *ff.positions()[6] << std::endl;
1249 #endif
1250 std::cerr << " done" << std::endl;
1251 }
1252
testUFFDistanceConstraints()1253 void testUFFDistanceConstraints() {
1254 std::cerr << "-------------------------------------" << std::endl;
1255 std::cerr << "Unit tests for UFF distance constraint terms." << std::endl;
1256
1257 ForceFields::ForceField ff;
1258 RDGeom::Point3D p1(0, 0, 0), p2(1.514, 0, 0);
1259 RDGeom::PointPtrVect &ps = ff.positions();
1260 ps.push_back(&p1);
1261 ps.push_back(&p2);
1262
1263 double *p, *g;
1264 p = new double[6];
1265 g = new double[6];
1266 for (int i = 0; i < 6; i++) {
1267 p[i] = 0.0;
1268 g[i] = 0.0;
1269 }
1270 p[0] = 0;
1271 p[3] = 1.40;
1272
1273 ff.initialize();
1274
1275 // C_3 - C_3, r0=1.514, k01=699.5918
1276 ForceFields::ForceFieldContrib *bs;
1277 bs = new ForceFields::UFF::DistanceConstraintContrib(&ff, 0, 1, 1.35, 1.55,
1278 1000.0);
1279 ff.contribs().push_back(ForceFields::ContribPtr(bs));
1280 double E;
1281 E = bs->getEnergy(p);
1282 TEST_ASSERT(RDKit::feq(E, 0.0));
1283 bs->getGrad(p, g);
1284 for (int i = 0; i < 6; i++) {
1285 TEST_ASSERT(RDKit::feq(g[i], 0.0));
1286 }
1287
1288 ff.initialize();
1289 (*ff.positions()[1])[0] = 1.20;
1290 p[3] = 1.20;
1291 E = bs->getEnergy(p);
1292 TEST_ASSERT(RDKit::feq(E, 11.25));
1293 bs->getGrad(p, g);
1294 TEST_ASSERT(RDKit::feq(g[0], 150.0));
1295 TEST_ASSERT(RDKit::feq(g[3], -150.0));
1296 TEST_ASSERT(RDKit::feq(g[1], 0.0));
1297 TEST_ASSERT(RDKit::feq(g[2], 0.0));
1298 TEST_ASSERT(RDKit::feq(g[4], 0.0));
1299 TEST_ASSERT(RDKit::feq(g[5], 0.0));
1300
1301 // try a bit of minimization
1302 RDGeom::Point3D d;
1303 ff.initialize();
1304 (*ff.positions()[1])[2] = 0.0;
1305 (*ff.positions()[1])[0] = 1.20;
1306 ff.minimize(10, 1e-8);
1307 d = *(RDGeom::Point3D *)ff.positions()[0] -
1308 *(RDGeom::Point3D *)ff.positions()[1];
1309 TEST_ASSERT(d.length() >= 1.35)
1310 TEST_ASSERT(d.length() <= 1.55)
1311
1312 ff.initialize();
1313 (*ff.positions()[1])[2] = 0.0;
1314 (*ff.positions()[1])[0] = 1.70;
1315 ff.minimize(10, 1e-8);
1316 d = *(RDGeom::Point3D *)ff.positions()[0] -
1317 *(RDGeom::Point3D *)ff.positions()[1];
1318 TEST_ASSERT(d.length() >= 1.35)
1319 TEST_ASSERT(d.length() <= 1.55)
1320
1321 delete[] p;
1322 delete[] g;
1323 std::cerr << " done" << std::endl;
1324 }
1325
testUFFAllConstraints()1326 void testUFFAllConstraints() {
1327 std::cerr << "-------------------------------------" << std::endl;
1328 std::cerr << "Unit tests for all UFF constraint terms." << std::endl;
1329
1330 std::string molBlock =
1331 "butane\n"
1332 " RDKit 3D\n"
1333 "butane\n"
1334 " 17 16 0 0 0 0 0 0 0 0999 V2000\n"
1335 " 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1336 " 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1337 " 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1338 " 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1339 " 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1340 " 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1341 " 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1342 " 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1343 " 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1344 " 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1345 " 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1346 " 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1347 " 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1348 " 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1349 " -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1350 " -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1351 " -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1352 " 1 2 1 0 0 0 0\n"
1353 " 1 15 1 0 0 0 0\n"
1354 " 1 16 1 0 0 0 0\n"
1355 " 1 17 1 0 0 0 0\n"
1356 " 2 3 1 0 0 0 0\n"
1357 " 2 4 1 0 0 0 0\n"
1358 " 2 14 1 0 0 0 0\n"
1359 " 4 5 1 0 0 0 0\n"
1360 " 4 6 1 0 0 0 0\n"
1361 " 4 7 1 0 0 0 0\n"
1362 " 7 8 1 0 0 0 0\n"
1363 " 7 9 1 0 0 0 0\n"
1364 " 7 13 1 0 0 0 0\n"
1365 " 9 10 1 0 0 0 0\n"
1366 " 9 11 1 0 0 0 0\n"
1367 " 9 12 1 0 0 0 0\n"
1368 "M END\n";
1369 RDKit::RWMol *mol;
1370 ForceFields::ForceField *field;
1371
1372 // distance constraints
1373 ForceFields::UFF::DistanceConstraintContrib *dc;
1374 mol = RDKit::MolBlockToMol(molBlock, true, false);
1375 TEST_ASSERT(mol);
1376 MolTransforms::setBondLength(mol->getConformer(), 1, 3, 2.0);
1377 field = RDKit::UFF::constructForceField(*mol);
1378 TEST_ASSERT(field);
1379 field->initialize();
1380 dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, 2.0, 2.0,
1381 1.0e5);
1382 field->contribs().push_back(ForceFields::ContribPtr(dc));
1383 field->minimize();
1384 TEST_ASSERT(RDKit::feq(
1385 MolTransforms::getBondLength(mol->getConformer(), 1, 3), 2.0, 0.1));
1386 delete field;
1387 field = RDKit::UFF::constructForceField(*mol);
1388 field->initialize();
1389 dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, true, -0.2,
1390 0.2, 1.0e5);
1391 field->contribs().push_back(ForceFields::ContribPtr(dc));
1392 field->minimize();
1393 TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.79);
1394 delete field;
1395 delete mol;
1396
1397 // angle constraints
1398 ForceFields::UFF::AngleConstraintContrib *ac;
1399 mol = RDKit::MolBlockToMol(molBlock, true, false);
1400 TEST_ASSERT(mol);
1401 MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 90.0);
1402 field = RDKit::UFF::constructForceField(*mol);
1403 TEST_ASSERT(field);
1404 field->initialize();
1405 ac = new ForceFields::UFF::AngleConstraintContrib(field, 1, 3, 6, 90.0, 90.0,
1406 100.0);
1407 field->contribs().push_back(ForceFields::ContribPtr(ac));
1408 field->minimize();
1409 TEST_ASSERT(RDKit::feq(
1410 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 90.0, 0.5));
1411 delete field;
1412 field = RDKit::UFF::constructForceField(*mol);
1413 field->initialize();
1414 ac = new ForceFields::UFF::AngleConstraintContrib(field, 1, 3, 6, true, -10.0,
1415 10.0, 100.0);
1416 field->contribs().push_back(ForceFields::ContribPtr(ac));
1417 field->minimize();
1418 TEST_ASSERT(RDKit::feq(
1419 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 100.0, 0.5));
1420 delete field;
1421 MolTransforms::setAngleDeg(mol->getConformer(), 1, 3, 6, 0.0);
1422 field = RDKit::UFF::constructForceField(*mol);
1423 field->initialize();
1424 ac = new ForceFields::UFF::AngleConstraintContrib(field, 1, 3, 6, false,
1425 -10.0, 10.0, 100.0);
1426 field->contribs().push_back(ForceFields::ContribPtr(ac));
1427 field->minimize();
1428 TEST_ASSERT(RDKit::feq(
1429 MolTransforms::getAngleDeg(mol->getConformer(), 1, 3, 6), 10.0, 0.5));
1430 delete field;
1431 delete mol;
1432
1433 // torsion constraints
1434 ForceFields::UFF::TorsionConstraintContrib *tc;
1435 mol = RDKit::MolBlockToMol(molBlock, true, false);
1436 TEST_ASSERT(mol);
1437 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, 15.0);
1438 field = RDKit::UFF::constructForceField(*mol);
1439 TEST_ASSERT(field);
1440 field->initialize();
1441 tc = new ForceFields::UFF::TorsionConstraintContrib(field, 1, 3, 6, 8, 10.0,
1442 10.0, 100.0);
1443 field->contribs().push_back(ForceFields::ContribPtr(tc));
1444 field->minimize();
1445 TEST_ASSERT(
1446 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1447 10.0, 0.5));
1448 delete field;
1449 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -30.0);
1450 field = RDKit::UFF::constructForceField(*mol);
1451 field->initialize();
1452 tc = new ForceFields::UFF::TorsionConstraintContrib(field, 1, 3, 6, 8, true,
1453 -1.0, 1.0, 100.0);
1454 field->contribs().push_back(ForceFields::ContribPtr(tc));
1455 field->minimize();
1456 TEST_ASSERT(
1457 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1458 -30.0, 1.5));
1459 delete field;
1460 MolTransforms::setDihedralDeg(mol->getConformer(), 1, 3, 6, 8, -10.0);
1461 field = RDKit::UFF::constructForceField(*mol);
1462 field->initialize();
1463 tc = new ForceFields::UFF::TorsionConstraintContrib(field, 1, 3, 6, 8, false,
1464 -10.0, 8.0, 100.0);
1465 field->contribs().push_back(ForceFields::ContribPtr(tc));
1466 field->minimize(500);
1467 TEST_ASSERT(
1468 RDKit::feq(MolTransforms::getDihedralDeg(mol->getConformer(), 1, 3, 6, 8),
1469 -9.0, 2.0));
1470 delete field;
1471 delete mol;
1472
1473 // position constraints
1474 ForceFields::UFF::PositionConstraintContrib *pc;
1475 mol = RDKit::MolBlockToMol(molBlock, true, false);
1476 TEST_ASSERT(mol);
1477 field = RDKit::UFF::constructForceField(*mol);
1478 TEST_ASSERT(field);
1479 field->initialize();
1480 RDGeom::Point3D p = mol->getConformer().getAtomPos(1);
1481 pc = new ForceFields::UFF::PositionConstraintContrib(field, 1, 0.3, 1.0e5);
1482 field->contribs().push_back(ForceFields::ContribPtr(pc));
1483 field->minimize();
1484 RDGeom::Point3D q = mol->getConformer().getAtomPos(1);
1485 TEST_ASSERT((p - q).length() < 0.3);
1486 delete field;
1487 delete mol;
1488
1489 // fixed atoms
1490 mol = RDKit::MolBlockToMol(molBlock, true, false);
1491 TEST_ASSERT(mol);
1492 field = RDKit::UFF::constructForceField(*mol);
1493 TEST_ASSERT(field);
1494 field->initialize();
1495 RDGeom::Point3D fp = mol->getConformer().getAtomPos(1);
1496 field->fixedPoints().push_back(1);
1497 field->minimize();
1498 RDGeom::Point3D fq = mol->getConformer().getAtomPos(1);
1499 TEST_ASSERT((fp - fq).length() < 0.01);
1500 delete field;
1501 delete mol;
1502
1503 std::cerr << " done" << std::endl;
1504 }
1505
testUFFCopy()1506 void testUFFCopy() {
1507 std::cerr << "-------------------------------------" << std::endl;
1508 std::cerr << "Unit tests for copying UFF ForceFields." << std::endl;
1509
1510 std::string molBlock =
1511 "butane\n"
1512 " RDKit 3D\n"
1513 "butane\n"
1514 " 17 16 0 0 0 0 0 0 0 0999 V2000\n"
1515 " 0.0000 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1516 " 1.4280 0.0000 0.0000 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1517 " 1.7913 -0.2660 0.9927 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1518 " 1.9040 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1519 " 1.5407 2.0271 0.3782 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1520 " 1.5407 1.5664 -1.3411 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1521 " 3.3320 1.3004 -0.3485 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1522 " 3.6953 1.5162 -1.3532 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1523 " 3.8080 0.0192 0.0649 C 0 0 0 0 0 0 0 0 0 0 0 0\n"
1524 " 3.4447 -0.7431 -0.6243 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1525 " 3.4447 -0.1966 1.0697 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1526 " 4.8980 0.0192 0.0649 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1527 " 3.6954 2.0627 0.3408 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1528 " 1.7913 -0.7267 -0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1529 " -0.3633 0.7267 0.7267 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1530 " -0.3633 -0.9926 0.2660 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1531 " -0.3633 0.2660 -0.9926 H 0 0 0 0 0 0 0 0 0 0 0 0\n"
1532 " 1 2 1 0 0 0 0\n"
1533 " 1 15 1 0 0 0 0\n"
1534 " 1 16 1 0 0 0 0\n"
1535 " 1 17 1 0 0 0 0\n"
1536 " 2 3 1 0 0 0 0\n"
1537 " 2 4 1 0 0 0 0\n"
1538 " 2 14 1 0 0 0 0\n"
1539 " 4 5 1 0 0 0 0\n"
1540 " 4 6 1 0 0 0 0\n"
1541 " 4 7 1 0 0 0 0\n"
1542 " 7 8 1 0 0 0 0\n"
1543 " 7 9 1 0 0 0 0\n"
1544 " 7 13 1 0 0 0 0\n"
1545 " 9 10 1 0 0 0 0\n"
1546 " 9 11 1 0 0 0 0\n"
1547 " 9 12 1 0 0 0 0\n"
1548 "M END\n";
1549 {
1550 RDKit::RWMol *mol = RDKit::MolBlockToMol(molBlock, true, false);
1551 TEST_ASSERT(mol);
1552 auto *cmol = new RDKit::RWMol(*mol);
1553 TEST_ASSERT(cmol);
1554
1555 ForceFields::ForceField *field = RDKit::UFF::constructForceField(*mol);
1556 TEST_ASSERT(field);
1557 field->initialize();
1558 auto *dc = new ForceFields::UFF::DistanceConstraintContrib(field, 1, 3, 2.0,
1559 2.0, 1.0e5);
1560 field->contribs().push_back(ForceFields::ContribPtr(dc));
1561 field->minimize();
1562 TEST_ASSERT(MolTransforms::getBondLength(mol->getConformer(), 1, 3) > 1.99);
1563
1564 auto *cfield = new ForceFields::ForceField(*field);
1565 cfield->positions().clear();
1566
1567 for (unsigned int i = 0; i < cmol->getNumAtoms(); i++) {
1568 cfield->positions().push_back(&cmol->getConformer().getAtomPos(i));
1569 }
1570 cfield->initialize();
1571 cfield->minimize();
1572 TEST_ASSERT(MolTransforms::getBondLength(cmol->getConformer(), 1, 3) >
1573 1.99);
1574 TEST_ASSERT(RDKit::feq(field->calcEnergy(), cfield->calcEnergy()));
1575
1576 const RDKit::Conformer &conf = mol->getConformer();
1577 const RDKit::Conformer &cconf = cmol->getConformer();
1578 for (unsigned int i = 0; i < mol->getNumAtoms(); i++) {
1579 RDGeom::Point3D p = conf.getAtomPos(i);
1580 RDGeom::Point3D cp = cconf.getAtomPos(i);
1581 TEST_ASSERT(RDKit::feq(p.x, cp.x));
1582 TEST_ASSERT(RDKit::feq(p.y, cp.y));
1583 TEST_ASSERT(RDKit::feq(p.z, cp.z));
1584 }
1585
1586 delete field;
1587 delete cfield;
1588 delete mol;
1589 delete cmol;
1590 }
1591 std::cerr << " done" << std::endl;
1592 }
1593
testUFFButaneScan()1594 void testUFFButaneScan() {
1595 std::cerr << "-------------------------------------" << std::endl;
1596 std::cerr << "Unit test for UFF butane scan." << std::endl;
1597
1598 auto molblock = R"(
1599 RDKit 3D
1600
1601 14 13 0 0 0 0 0 0 0 0999 V2000
1602 -1.5575 0.5684 -0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
1603 -0.6987 -0.6064 0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
1604 0.6987 -0.6064 -0.3102 C 0 0 0 0 0 0 0 0 0 0 0 0
1605 1.5575 0.5684 0.1320 C 0 0 0 0 0 0 0 0 0 0 0 0
1606 -1.6308 0.6129 -1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
1607 -2.5700 0.4662 0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
1608 -1.1524 1.5190 0.2272 H 0 0 0 0 0 0 0 0 0 0 0 0
1609 -1.2059 -1.5359 0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
1610 -0.6215 -0.6099 1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
1611 0.6215 -0.6099 -1.4037 H 0 0 0 0 0 0 0 0 0 0 0 0
1612 1.2059 -1.5359 -0.0253 H 0 0 0 0 0 0 0 0 0 0 0 0
1613 1.6308 0.6129 1.2232 H 0 0 0 0 0 0 0 0 0 0 0 0
1614 1.1524 1.5190 -0.2271 H 0 0 0 0 0 0 0 0 0 0 0 0
1615 2.5700 0.4662 -0.2715 H 0 0 0 0 0 0 0 0 0 0 0 0
1616 1 2 1 0
1617 2 3 1 0
1618 3 4 1 0
1619 1 5 1 0
1620 1 6 1 0
1621 1 7 1 0
1622 2 8 1 0
1623 2 9 1 0
1624 3 10 1 0
1625 3 11 1 0
1626 4 12 1 0
1627 4 13 1 0
1628 4 14 1 0
1629 M END
1630 )";
1631
1632 // torsion constraints
1633 std::unique_ptr<ROMol> mol(MolBlockToMol(molblock, true, false));
1634 TEST_ASSERT(mol);
1635 std::vector<std::pair<double, double>> diheEnergyVect;
1636 std::vector<std::pair<double, double>> expectedDiheEnergyVect{
1637 {179.9997, 1.7649}, {-175.1001, 1.8314}, {-170.1003, 2.0318},
1638 {-165.1004, 2.3527}, {-160.1005, 2.7726}, {-155.1005, 3.2633},
1639 {-150.1005, 3.7919}, {-145.1005, 4.3220}, {-140.1005, 4.8159},
1640 {-135.1004, 5.2372}, {-130.1003, 5.5532}, {-125.1001, 5.7385},
1641 {-119.9000, 5.7768}, {-114.8998, 5.6630}, {-109.8997, 5.4141},
1642 {-104.8996, 5.0555}, {-99.8995, 4.6243}, {-94.8995, 4.1667},
1643 {-89.8996, 3.7338}, {-84.8997, 3.3722}, {-79.8998, 3.1102},
1644 {-74.8999, 2.9561}, {-70.1000, 2.9112}, {-65.1001, 2.9777},
1645 {-60.1003, 3.1680}, {-55.1004, 3.4865}, {-50.1005, 3.9304},
1646 {-45.1006, 4.4873}, {-40.1007, 5.1363}, {-35.1007, 5.8490},
1647 {-30.1007, 6.5913}, {-25.1007, 7.3249}, {-20.1006, 8.0089},
1648 {-15.1005, 8.6016}, {-10.1004, 9.0626}, {-5.1002, 9.3573},
1649 {-0.1000, 9.4612}, {5.1002, 9.3573}, {10.1004, 9.0626},
1650 {15.1005, 8.6016}, {20.1006, 8.0089}, {25.1007, 7.3249},
1651 {30.1007, 6.5913}, {35.1007, 5.8490}, {40.1007, 5.1363},
1652 {45.1006, 4.4873}, {50.1005, 3.9304}, {55.1004, 3.4865},
1653 {60.1003, 3.1680}, {65.1001, 2.9777}, {70.1000, 2.9112},
1654 {74.8999, 2.9561}, {79.8998, 3.1102}, {84.8997, 3.3722},
1655 {89.8996, 3.7338}, {94.8995, 4.1667}, {99.8995, 4.6243},
1656 {104.8996, 5.0555}, {109.8997, 5.4141}, {114.8998, 5.6630},
1657 {119.9000, 5.7768}, {125.1001, 5.7385}, {130.1003, 5.5532},
1658 {135.1004, 5.2372}, {140.1005, 4.8159}, {145.1005, 4.3220},
1659 {150.1005, 3.7919}, {155.1005, 3.2633}, {160.1005, 2.7726},
1660 {165.1004, 2.3527}, {170.1003, 2.0318}, {175.1001, 1.8314}};
1661 std::unique_ptr<ForceFields::ForceField> globalFF(
1662 RDKit::UFF::constructForceField(*mol));
1663 TEST_ASSERT(globalFF);
1664 globalFF->initialize();
1665 std::unique_ptr<ROMol> mol2(new ROMol(*mol));
1666 const int start = -180;
1667 const int end = 180;
1668 const int incr = 5;
1669 MolTransforms::setDihedralDeg(mol2->getConformer(), 0, 1, 2, 3, start);
1670 for (int diheIn = start; diheIn < end; diheIn += incr) {
1671 std::unique_ptr<ForceFields::ForceField> localFF(
1672 RDKit::UFF::constructForceField(*mol2));
1673 TEST_ASSERT(localFF);
1674 localFF->initialize();
1675 auto *tc = new ForceFields::UFF::TorsionConstraintContrib(
1676 localFF.get(), 0, 1, 2, 3, static_cast<double>(diheIn) - 0.1,
1677 static_cast<double>(diheIn) + 0.1, 100.0);
1678 localFF->contribs().push_back(ForceFields::ContribPtr(tc));
1679 TEST_ASSERT(localFF->minimize(10000) == 0);
1680 std::vector<double> pos(3 * localFF->numPoints());
1681 TEST_ASSERT(localFF->numPoints() == globalFF->numPoints());
1682 auto posns = localFF->positions();
1683 for (unsigned int i = 0; i < localFF->numPoints(); ++i) {
1684 for (unsigned int j = 0; j < 3; ++j) {
1685 pos[i * 3 + j] = (*static_cast<RDGeom::Point3D *>(posns.at(i)))[j];
1686 }
1687 }
1688 auto diheOut =
1689 MolTransforms::getDihedralDeg(mol2->getConformer(), 0, 1, 2, 3);
1690 auto energy = globalFF->calcEnergy(pos.data());
1691 diheEnergyVect.emplace_back(diheOut, energy);
1692 }
1693 for (size_t i = 0; i < diheEnergyVect.size(); ++i) {
1694 TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).first,
1695 expectedDiheEnergyVect.at(i).first, 1.0));
1696 TEST_ASSERT(RDKit::feq(diheEnergyVect.at(i).second,
1697 expectedDiheEnergyVect.at(i).second, 0.2));
1698 }
1699
1700 std::cerr << " done" << std::endl;
1701 }
1702
main()1703 int main() {
1704 #if 1
1705 test1();
1706 testUFF1();
1707 testUFF2();
1708 testUFF3();
1709 testUFF4();
1710 testUFF5();
1711 testUFF6();
1712 testUFF7();
1713 testUFFParams();
1714 testUFF8();
1715 testUFFTorsionConflict();
1716
1717 #endif
1718 testUFFDistanceConstraints();
1719 testUFFAllConstraints();
1720 testUFFCopy();
1721 testUFFButaneScan();
1722 }
1723