1 //
2 // Copyright (C) 2013 Paolo Tosco
3 //
4 // Copyright (C) 2004-2006 Rational Discovery LLC
5 //
6 // @@ All Rights Reserved @@
7 // This file is part of the RDKit.
8 // The contents are covered by the terms of the BSD license
9 // which is included in the file license.txt, found at the root
10 // of the RDKit source tree.
11 //
12 #include <RDGeneral/export.h>
13 #ifndef __RD_MMFFPARAMS_H__
14 #define __RD_MMFFPARAMS_H__
15
16 #include <memory>
17 #include <RDGeneral/Invariant.h>
18 #include <cmath>
19 #include <string>
20 #include <vector>
21 #include <algorithm>
22 #include <map>
23 #include <iostream>
24 #include <cstdint>
25
26 #ifndef M_PI
27 #define M_PI 3.14159265358979323846
28 #endif
29
30 // binary searches are slightly faster than std::map;
31 // however when I moved to binary searches I had already
32 // written the code for std::map, so the two methods
33 // can be toggled defining RDKIT_MMFF_PARAMS_USE_STD_MAP
34
35 //#define RDKIT_MMFF_PARAMS_USE_STD_MAP 1
36
37 namespace ForceFields {
38 namespace MMFF {
39
40 const double DEG2RAD = M_PI / 180.0;
41 const double RAD2DEG = 180.0 / M_PI;
42 const double MDYNE_A_TO_KCAL_MOL = 143.9325;
isDoubleZero(const double x)43 inline bool isDoubleZero(const double x) {
44 return ((x < 1.0e-10) && (x > -1.0e-10));
45 }
clipToOne(double & x)46 inline void clipToOne(double &x) {
47 if (x > 1.0) {
48 x = 1.0;
49 } else if (x < -1.0) {
50 x = -1.0;
51 }
52 }
53
54 //! class to store MMFF atom type equivalence levels
55 class RDKIT_FORCEFIELD_EXPORT MMFFDef {
56 public:
57 std::uint8_t eqLevel[4];
58 };
59
60 //! class to store MMFF Properties
61 class RDKIT_FORCEFIELD_EXPORT MMFFProp {
62 public:
63 std::uint8_t atno;
64 std::uint8_t crd;
65 std::uint8_t val;
66 std::uint8_t pilp;
67 std::uint8_t mltb;
68 std::uint8_t arom;
69 std::uint8_t linh;
70 std::uint8_t sbmb;
71 };
72
73 //! class to store MMFF Partial Bond Charge Increments
74 class RDKIT_FORCEFIELD_EXPORT MMFFPBCI {
75 public:
76 double pbci;
77 double fcadj;
78 };
79
80 //! class to store MMFF bond-charge-increment parameters used to
81 //! construct MMFF partial atomic charges
82 class RDKIT_FORCEFIELD_EXPORT MMFFChg {
83 public:
84 double bci;
85 };
86
87 //! class to store MMFF parameters for bond stretching
88 class RDKIT_FORCEFIELD_EXPORT MMFFBond {
89 public:
90 double kb;
91 double r0;
92 };
93
94 //! class to store parameters for Herschbach-Laurie's version
95 //! of Badger's rule
96 class RDKIT_FORCEFIELD_EXPORT MMFFHerschbachLaurie {
97 public:
98 double a_ij;
99 double d_ij;
100 double dp_ij;
101 };
102
103 //! class to store covalent radius and Pauling electronegativity
104 //! values for MMFF bond stretching empirical rule
105 class RDKIT_FORCEFIELD_EXPORT MMFFCovRadPauEle {
106 public:
107 double r0;
108 double chi;
109 };
110
111 //! class to store MMFF parameters for angle bending
112 class RDKIT_FORCEFIELD_EXPORT MMFFAngle {
113 public:
114 double ka;
115 double theta0;
116 };
117
118 //! class to store MMFF parameters for stretch-bending
119 class RDKIT_FORCEFIELD_EXPORT MMFFStbn {
120 public:
121 double kbaIJK;
122 double kbaKJI;
123 };
124
125 //! class to store MMFF parameters for out-of-plane bending
126 class RDKIT_FORCEFIELD_EXPORT MMFFOop {
127 public:
128 double koop;
129 };
130
131 //! class to store MMFF parameters for torsions
132 class RDKIT_FORCEFIELD_EXPORT MMFFTor {
133 public:
134 double V1;
135 double V2;
136 double V3;
137 };
138
139 //! class to store MMFF parameters for non-bonded Van der Waals
140 class RDKIT_FORCEFIELD_EXPORT MMFFVdW {
141 public:
142 double alpha_i;
143 double N_i;
144 double A_i;
145 double G_i;
146 double R_star;
147 std::uint8_t DA;
148 };
149
150 class RDKIT_FORCEFIELD_EXPORT MMFFVdWRijstarEps {
151 public:
152 double R_ij_starUnscaled;
153 double epsilonUnscaled;
154 double R_ij_star;
155 double epsilon;
156 };
157
158 class RDKIT_FORCEFIELD_EXPORT MMFFAromCollection {
159 public:
160 //! Looks up the parameters for a particular key and returns them.
161 /*!
162 \return a pointer to the MMFFArom object, NULL on failure.
163 */
isMMFFAromatic(const unsigned int atomType)164 bool isMMFFAromatic(const unsigned int atomType) const {
165 return ((std::find(d_params.begin(), d_params.end(), atomType) !=
166 d_params.end())
167 ? true
168 : false);
169 }
170
171 MMFFAromCollection(const std::vector<std::uint8_t> *mmffArom = nullptr);
172 std::vector<std::uint8_t> d_params; //!< the aromatic type vector
173 };
174
175 class RDKIT_FORCEFIELD_EXPORT MMFFDefCollection {
176 public:
177 //! Looks up the parameters for a particular key and returns them.
178 /*!
179 \return a pointer to the MMFFDef object, NULL on failure.
180 */
operator()181 const MMFFDef *operator()(const unsigned int atomType) const {
182 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
183 std::map<const unsigned int, MMFFDef>::const_iterator res;
184 res = d_params.find(atomType);
185
186 return ((res != d_params.end()) ? &((*res).second) : NULL);
187 #else
188 return ((atomType && (atomType <= d_params.size()))
189 ? &d_params[atomType - 1]
190 : nullptr);
191 #endif
192 }
193
194 MMFFDefCollection(std::string mmffDef = "");
195
196 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
197 std::map<const unsigned int, MMFFDef> d_params; //!< the parameter map
198 #else
199 std::vector<MMFFDef> d_params; //!< the parameter vector
200 #endif
201 };
202
203 class RDKIT_FORCEFIELD_EXPORT MMFFPropCollection {
204 public:
205 //! Looks up the parameters for a particular key and returns them.
206 /*!
207 \return a pointer to the MMFFProp object, NULL on failure.
208 */
operator()209 const MMFFProp *operator()(const unsigned int atomType) const {
210 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
211 std::map<const unsigned int, MMFFProp>::const_iterator res;
212 res = d_params.find(atomType);
213
214 return ((res != d_params.end()) ? &((*res).second) : NULL);
215 #else
216 std::pair<std::vector<std::uint8_t>::const_iterator,
217 std::vector<std::uint8_t>::const_iterator>
218 bounds =
219 std::equal_range(d_iAtomType.begin(), d_iAtomType.end(), atomType);
220
221 return ((bounds.first != bounds.second)
222 ? &d_params[bounds.first - d_iAtomType.begin()]
223 : nullptr);
224 #endif
225 }
226
227 MMFFPropCollection(std::string mmffProp = "");
228 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
229 std::map<const unsigned int, MMFFProp> d_params; //!< the parameter map
230 #else
231 std::vector<MMFFProp> d_params;
232 std::vector<std::uint8_t> d_iAtomType; //!< the parameter vector
233 #endif
234 };
235
236 class RDKIT_FORCEFIELD_EXPORT MMFFPBCICollection {
237 public:
238 //! Looks up the parameters for a particular key and returns them.
239 /*!
240 \return a pointer to the MMFFPBCI object, NULL on failure.
241 */
operator()242 const MMFFPBCI *operator()(const unsigned int atomType) const {
243 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
244 std::map<const unsigned int, MMFFPBCI>::const_iterator res;
245 res = d_params.find(atomType);
246
247 return ((res != d_params.end()) ? &((*res).second) : NULL);
248 #else
249 return ((atomType && (atomType <= d_params.size()))
250 ? &d_params[atomType - 1]
251 : nullptr);
252 #endif
253 }
254
255 MMFFPBCICollection(std::string mmffPBCI = "");
256
257 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
258 std::map<const unsigned int, MMFFPBCI> d_params; //!< the parameter map
259 #else
260 std::vector<MMFFPBCI> d_params; //!< the parameter vector
261 #endif
262 };
263
264 class RDKIT_FORCEFIELD_EXPORT MMFFChgCollection {
265 public:
266 //! Looks up the parameters for a particular key and returns them.
267 /*!
268 \return a pointer to the MMFFChg object, NULL on failure.
269 */
getMMFFChgParams(const unsigned int bondType,const unsigned int iAtomType,const unsigned int jAtomType)270 const std::pair<int, const MMFFChg *> getMMFFChgParams(
271 const unsigned int bondType, const unsigned int iAtomType,
272 const unsigned int jAtomType) const {
273 int sign = -1;
274 const MMFFChg *mmffChgParams = nullptr;
275 unsigned int canIAtomType = iAtomType;
276 unsigned int canJAtomType = jAtomType;
277 if (iAtomType > jAtomType) {
278 canIAtomType = jAtomType;
279 canJAtomType = iAtomType;
280 sign = 1;
281 }
282 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
283 std::map<const unsigned int,
284 std::map<const unsigned int, MMFFChg>>::const_iterator res1;
285 std::map<const unsigned int, MMFFChg>::const_iterator res2;
286 res1 = d_params[bondType].find(canIAtomType);
287 if (res1 != d_params[bondType].end()) {
288 res2 = ((*res1).second).find(canJAtomType);
289 if (res2 != ((*res1).second).end()) {
290 mmffChgParams = &((*res2).second);
291 }
292 }
293 #else
294 std::pair<std::vector<std::uint8_t>::const_iterator,
295 std::vector<std::uint8_t>::const_iterator>
296 bounds;
297
298 bounds =
299 std::equal_range(d_iAtomType.begin(), d_iAtomType.end(), canIAtomType);
300 if (bounds.first != bounds.second) {
301 bounds = std::equal_range(
302 d_jAtomType.begin() + (bounds.first - d_iAtomType.begin()),
303 d_jAtomType.begin() + (bounds.second - d_iAtomType.begin()),
304 canJAtomType);
305 if (bounds.first != bounds.second) {
306 bounds = std::equal_range(
307 d_bondType.begin() + (bounds.first - d_jAtomType.begin()),
308 d_bondType.begin() + (bounds.second - d_jAtomType.begin()),
309 bondType);
310 if (bounds.first != bounds.second) {
311 mmffChgParams = &d_params[bounds.first - d_bondType.begin()];
312 }
313 }
314 }
315 #endif
316
317 return std::make_pair(sign, mmffChgParams);
318 }
319
320 MMFFChgCollection(std::string mmffChg = "");
321
322 //!< the parameter 3D-map
323 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
324 std::map<const unsigned int,
325 std::map<const unsigned int, std::map<const unsigned int, MMFFChg>>>
326 d_params; //!< the parameter 3D-map
327 #else
328 std::vector<MMFFChg> d_params; //! the parameter vector
329 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
330 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
331 std::vector<std::uint8_t> d_bondType; //! bond type vector for bond i-j
332 #endif
333 };
334
335 class RDKIT_FORCEFIELD_EXPORT MMFFBondCollection {
336 public:
337 //! Looks up the parameters for a particular key and returns them.
338 /*!
339 \return a pointer to the MMFFBond object, NULL on failure.
340 */
operator()341 const MMFFBond *operator()(const unsigned int bondType,
342 const unsigned int atomType,
343 const unsigned int nbrAtomType) const {
344 const MMFFBond *mmffBondParams = nullptr;
345 unsigned int canAtomType = atomType;
346 unsigned int canNbrAtomType = nbrAtomType;
347 if (atomType > nbrAtomType) {
348 canAtomType = nbrAtomType;
349 canNbrAtomType = atomType;
350 }
351 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
352 std::map<const unsigned int,
353 std::map<const unsigned int,
354 std::map<const unsigned int, MMFFBond>>>::const_iterator
355 res1;
356 std::map<const unsigned int,
357 std::map<const unsigned int, MMFFBond>>::const_iterator res2;
358 std::map<const unsigned int, MMFFBond>::const_iterator res3;
359 res1 = d_params.find(bondType);
360 if (res1 != d_params.end()) {
361 res2 = ((*res1).second).find(canAtomType);
362 if (res2 != ((*res1).second).end()) {
363 res3 = ((*res2).second).find(canNbrAtomType);
364 if (res3 != ((*res2).second).end()) {
365 mmffBondParams = &((*res3).second);
366 }
367 }
368 }
369 #else
370 std::pair<std::vector<std::uint8_t>::const_iterator,
371 std::vector<std::uint8_t>::const_iterator>
372 bounds;
373 bounds =
374 std::equal_range(d_iAtomType.begin(), d_iAtomType.end(), canAtomType);
375 if (bounds.first != bounds.second) {
376 bounds = std::equal_range(
377 d_jAtomType.begin() + (bounds.first - d_iAtomType.begin()),
378 d_jAtomType.begin() + (bounds.second - d_iAtomType.begin()),
379 canNbrAtomType);
380 if (bounds.first != bounds.second) {
381 bounds = std::equal_range(
382 d_bondType.begin() + (bounds.first - d_jAtomType.begin()),
383 d_bondType.begin() + (bounds.second - d_jAtomType.begin()),
384 bondType);
385 if (bounds.first != bounds.second) {
386 mmffBondParams = &d_params[bounds.first - d_bondType.begin()];
387 }
388 }
389 }
390 #endif
391
392 return mmffBondParams;
393 }
394
395 MMFFBondCollection(std::string mmffBond = "");
396 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
397 std::map<const unsigned int,
398 std::map<const unsigned int, std::map<const unsigned int, MMFFBond>>>
399 d_params; //!< the parameter 3D-map
400 #else
401 std::vector<MMFFBond> d_params; //!< the parameter vector
402 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
403 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
404 std::vector<std::uint8_t> d_bondType; //! bond type vector for bond i-j
405 #endif
406 };
407
408 class RDKIT_FORCEFIELD_EXPORT MMFFBndkCollection {
409 public:
410 //! Looks up the parameters for a particular key and returns them.
411 /*!
412 \return a pointer to the MMFFBndk object, NULL on failure.
413 */
operator()414 const MMFFBond *operator()(const int atomicNum,
415 const int nbrAtomicNum) const {
416 const MMFFBond *mmffBndkParams = nullptr;
417 unsigned int canAtomicNum = atomicNum;
418 unsigned int canNbrAtomicNum = nbrAtomicNum;
419 if (atomicNum > nbrAtomicNum) {
420 canAtomicNum = nbrAtomicNum;
421 canNbrAtomicNum = atomicNum;
422 }
423 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
424 std::map<const unsigned int,
425 std::map<const unsigned int, MMFFBond>>::const_iterator res1;
426 std::map<const unsigned int, MMFFBond>::const_iterator res2;
427 res1 = d_params.find(canAtomicNum);
428 if (res1 != d_params.end()) {
429 res2 = ((*res1).second).find(canNbrAtomicNum);
430 if (res2 != ((*res1).second).end()) {
431 mmffBndkParams = &((*res2).second);
432 }
433 }
434 #else
435 std::pair<std::vector<std::uint8_t>::const_iterator,
436 std::vector<std::uint8_t>::const_iterator>
437 bounds;
438 bounds = std::equal_range(d_iAtomicNum.begin(), d_iAtomicNum.end(),
439 canAtomicNum);
440 if (bounds.first != bounds.second) {
441 bounds = std::equal_range(
442 d_jAtomicNum.begin() + (bounds.first - d_iAtomicNum.begin()),
443 d_jAtomicNum.begin() + (bounds.second - d_iAtomicNum.begin()),
444 canNbrAtomicNum);
445 if (bounds.first != bounds.second) {
446 mmffBndkParams = &d_params[bounds.first - d_jAtomicNum.begin()];
447 }
448 }
449 #endif
450
451 return mmffBndkParams;
452 }
453
454 MMFFBndkCollection(std::string mmffBndk = "");
455 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
456 std::map<const unsigned int, std::map<const unsigned int, MMFFBond>>
457 d_params; //!< the parameter 2D-map
458 #else
459 std::vector<MMFFBond> d_params; //!< the parameter vector
460 std::vector<std::uint8_t> d_iAtomicNum; //! atomic number vector for atom i
461 std::vector<std::uint8_t> d_jAtomicNum; //! atomic number vector for atom j
462 #endif
463 };
464
465 class RDKIT_FORCEFIELD_EXPORT MMFFHerschbachLaurieCollection {
466 public:
467 //! Looks up the parameters for a particular key and returns them.
468 /*!
469 \return a pointer to the MMFFHerschbachLaurie object, NULL on failure.
470 */
operator()471 const MMFFHerschbachLaurie *operator()(const int iRow, const int jRow) const {
472 const MMFFHerschbachLaurie *mmffHerschbachLaurieParams = nullptr;
473 unsigned int canIRow = iRow;
474 unsigned int canJRow = jRow;
475 if (iRow > jRow) {
476 canIRow = jRow;
477 canJRow = iRow;
478 }
479 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
480 std::map<const unsigned int,
481 std::map<const unsigned int, MMFFHerschbachLaurie>>::const_iterator
482 res1;
483 std::map<const unsigned int, MMFFHerschbachLaurie>::const_iterator res2;
484 res1 = d_params.find(canIRow);
485 if (res1 != d_params.end()) {
486 res2 = ((*res1).second).find(canJRow);
487 if (res2 != ((*res1).second).end()) {
488 mmffHerschbachLaurieParams = &((*res2).second);
489 }
490 }
491 #else
492 std::pair<std::vector<std::uint8_t>::const_iterator,
493 std::vector<std::uint8_t>::const_iterator>
494 bounds;
495 bounds = std::equal_range(d_iRow.begin(), d_iRow.end(), canIRow);
496 if (bounds.first != bounds.second) {
497 bounds = std::equal_range(
498 d_jRow.begin() + (bounds.first - d_iRow.begin()),
499 d_jRow.begin() + (bounds.second - d_iRow.begin()), canJRow);
500 if (bounds.first != bounds.second) {
501 mmffHerschbachLaurieParams = &d_params[bounds.first - d_jRow.begin()];
502 }
503 }
504 #endif
505
506 return mmffHerschbachLaurieParams;
507 }
508
509 MMFFHerschbachLaurieCollection(std::string mmffHerschbachLaurie = "");
510
511 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
512 std::map<const unsigned int,
513 std::map<const unsigned int, MMFFHerschbachLaurie>>
514 d_params; //!< the parameter 2D-map
515 #else
516 std::vector<MMFFHerschbachLaurie> d_params; //!< the parameter vector
517 std::vector<std::uint8_t> d_iRow; //! periodic row number vector for atom i
518 std::vector<std::uint8_t> d_jRow; //! periodic row number vector for atom j
519 #endif
520 };
521
522 class RDKIT_FORCEFIELD_EXPORT MMFFCovRadPauEleCollection {
523 public:
524 //! Looks up the parameters for a particular key and returns them.
525 /*!
526 \return a pointer to the MMFFCovRadPauEle object, NULL on failure.
527 */
operator()528 const MMFFCovRadPauEle *operator()(const unsigned int atomicNum) const {
529 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
530 std::map<const unsigned int, MMFFCovRadPauEle>::const_iterator res;
531 res = d_params.find(atomicNum);
532
533 return ((res != d_params.end()) ? &((*res).second) : NULL);
534 #else
535 std::pair<std::vector<std::uint8_t>::const_iterator,
536 std::vector<std::uint8_t>::const_iterator>
537 bounds =
538 std::equal_range(d_atomicNum.begin(), d_atomicNum.end(), atomicNum);
539
540 return ((bounds.first != bounds.second)
541 ? &d_params[bounds.first - d_atomicNum.begin()]
542 : nullptr);
543 #endif
544 }
545
546 MMFFCovRadPauEleCollection(std::string mmffCovRadPauEle = "");
547 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
548 std::map<const unsigned int, MMFFCovRadPauEle>
549 d_params; //!< the parameter map
550 #else
551 std::vector<MMFFCovRadPauEle> d_params; //!< the parameter vector
552 std::vector<std::uint8_t> d_atomicNum; //!< the atomic number vector
553 #endif
554 };
555
556 class RDKIT_FORCEFIELD_EXPORT MMFFAngleCollection {
557 public:
558 //! Looks up the parameters for a particular key and returns them.
559 /*!
560 \return a pointer to the MMFFAngle object, NULL on failure.
561 */
operator()562 const MMFFAngle *operator()(const MMFFDefCollection *mmffDef,
563 const unsigned int angleType,
564 const unsigned int iAtomType,
565 const unsigned int jAtomType,
566 const unsigned int kAtomType) const {
567 const MMFFAngle *mmffAngleParams = nullptr;
568 unsigned int iter = 0;
569
570 // For bending of the i-j-k angle, a five-stage process based
571 // in the level combinations 1-1-1,2-2-2,3-2-3,4-2-4, and
572 // 5-2-5 is used. (MMFF.I, note 68, page 519)
573 // We skip 1-1-1 since Level 2 === Level 1
574 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
575 std::map<const unsigned int,
576 std::map<const unsigned int,
577 std::map<const unsigned int,
578 std::map<const unsigned int, MMFFAngle>>>>::
579 const_iterator res1;
580 std::map<const unsigned int,
581 std::map<const unsigned int,
582 std::map<const unsigned int, MMFFAngle>>>::const_iterator
583 res2;
584 std::map<const unsigned int,
585 std::map<const unsigned int, MMFFAngle>>::const_iterator res3;
586 std::map<const unsigned int, MMFFAngle>::const_iterator res4;
587 while ((iter < 4) && (!mmffAngleParams)) {
588 unsigned int canIAtomType = (*mmffDef)(iAtomType)->eqLevel[iter];
589 unsigned int canKAtomType = (*mmffDef)(kAtomType)->eqLevel[iter];
590 if (canIAtomType > canKAtomType) {
591 unsigned int temp = canKAtomType;
592 canKAtomType = canIAtomType;
593 canIAtomType = temp;
594 }
595 res1 = d_params.find(angleType);
596 if (res1 != d_params.end()) {
597 res2 = ((*res1).second).find(canIAtomType);
598 if (res2 != ((*res1).second).end()) {
599 res3 = ((*res2).second).find(jAtomType);
600 if (res3 != ((*res2).second).end()) {
601 res4 = ((*res3).second).find(canKAtomType);
602 if (res4 != ((*res3).second).end()) {
603 mmffAngleParams = &((*res4).second);
604 }
605 }
606 }
607 }
608 ++iter;
609 }
610 #else
611 std::pair<std::vector<std::uint8_t>::const_iterator,
612 std::vector<std::uint8_t>::const_iterator>
613 jBounds =
614 std::equal_range(d_jAtomType.begin(), d_jAtomType.end(), jAtomType);
615 std::pair<std::vector<std::uint8_t>::const_iterator,
616 std::vector<std::uint8_t>::const_iterator>
617 bounds;
618 if (jBounds.first != jBounds.second) {
619 while ((iter < 4) && (!mmffAngleParams)) {
620 unsigned int canIAtomType = (*mmffDef)(iAtomType)->eqLevel[iter];
621 unsigned int canKAtomType = (*mmffDef)(kAtomType)->eqLevel[iter];
622 if (canIAtomType > canKAtomType) {
623 unsigned int temp = canKAtomType;
624 canKAtomType = canIAtomType;
625 canIAtomType = temp;
626 }
627 bounds = std::equal_range(
628 d_iAtomType.begin() + (jBounds.first - d_jAtomType.begin()),
629 d_iAtomType.begin() + (jBounds.second - d_jAtomType.begin()),
630 canIAtomType);
631 if (bounds.first != bounds.second) {
632 bounds = std::equal_range(
633 d_kAtomType.begin() + (bounds.first - d_iAtomType.begin()),
634 d_kAtomType.begin() + (bounds.second - d_iAtomType.begin()),
635 canKAtomType);
636 if (bounds.first != bounds.second) {
637 bounds = std::equal_range(
638 d_angleType.begin() + (bounds.first - d_kAtomType.begin()),
639 d_angleType.begin() + (bounds.second - d_kAtomType.begin()),
640 angleType);
641 if (bounds.first != bounds.second) {
642 mmffAngleParams = &d_params[bounds.first - d_angleType.begin()];
643 }
644 }
645 }
646 ++iter;
647 }
648 }
649 #endif
650
651 return mmffAngleParams;
652 }
653
654 MMFFAngleCollection(std::string mmffAngle = "");
655 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
656 std::map<const unsigned int,
657 std::map<const unsigned int,
658 std::map<const unsigned int,
659 std::map<const unsigned int, MMFFAngle>>>>
660 d_params; //!< the parameter 4D-map
661 #else
662 std::vector<MMFFAngle> d_params; //!< the parameter vector
663 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
664 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
665 std::vector<std::uint8_t> d_kAtomType; //! atom type vector for atom k
666 std::vector<std::uint8_t> d_angleType; //! angle type vector for angle i-j-k
667 #endif
668 };
669
670 class RDKIT_FORCEFIELD_EXPORT MMFFStbnCollection {
671 public:
672 //! Looks up the parameters for a particular key and returns them.
673 /*!
674 \return a pointer to the MMFFStbn object, NULL on failure.
675 */
getMMFFStbnParams(const unsigned int stretchBendType,const unsigned int bondType1,const unsigned int bondType2,const unsigned int iAtomType,const unsigned int jAtomType,const unsigned int kAtomType)676 const std::pair<bool, const MMFFStbn *> getMMFFStbnParams(
677 const unsigned int stretchBendType, const unsigned int bondType1,
678 const unsigned int bondType2, const unsigned int iAtomType,
679 const unsigned int jAtomType, const unsigned int kAtomType) const {
680 const MMFFStbn *mmffStbnParams = nullptr;
681 bool swap = false;
682 unsigned int canIAtomType = iAtomType;
683 unsigned int canKAtomType = kAtomType;
684 unsigned int canStretchBendType = stretchBendType;
685 if (iAtomType > kAtomType) {
686 canIAtomType = kAtomType;
687 canKAtomType = iAtomType;
688 swap = true;
689 } else if (iAtomType == kAtomType) {
690 swap = (bondType1 < bondType2);
691 }
692 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
693 std::map<const unsigned int,
694 std::map<const unsigned int,
695 std::map<const unsigned int,
696 std::map<const unsigned int, MMFFStbn>>>>::
697 const_iterator res1;
698 std::map<const unsigned int,
699 std::map<const unsigned int,
700 std::map<const unsigned int, MMFFStbn>>>::const_iterator
701 res2;
702 std::map<const unsigned int,
703 std::map<const unsigned int, MMFFStbn>>::const_iterator res3;
704 std::map<const unsigned int, MMFFStbn>::const_iterator res4;
705 res1 = d_params.find(canStretchBendType);
706 if (res1 != d_params.end()) {
707 res2 = ((*res1).second).find(canIAtomType);
708 if (res2 != ((*res1).second).end()) {
709 res3 = ((*res2).second).find(jAtomType);
710 if (res3 != ((*res2).second).end()) {
711 res4 = ((*res3).second).find(canKAtomType);
712 if (res4 != ((*res3).second).end()) {
713 mmffStbnParams = &((*res4).second);
714 }
715 }
716 }
717 }
718 #else
719 std::pair<std::vector<std::uint8_t>::const_iterator,
720 std::vector<std::uint8_t>::const_iterator>
721 jBounds =
722 std::equal_range(d_jAtomType.begin(), d_jAtomType.end(), jAtomType);
723 std::pair<std::vector<std::uint8_t>::const_iterator,
724 std::vector<std::uint8_t>::const_iterator>
725 bounds;
726 if (jBounds.first != jBounds.second) {
727 bounds = std::equal_range(
728 d_iAtomType.begin() + (jBounds.first - d_jAtomType.begin()),
729 d_iAtomType.begin() + (jBounds.second - d_jAtomType.begin()),
730 canIAtomType);
731 if (bounds.first != bounds.second) {
732 bounds = std::equal_range(
733 d_kAtomType.begin() + (bounds.first - d_iAtomType.begin()),
734 d_kAtomType.begin() + (bounds.second - d_iAtomType.begin()),
735 canKAtomType);
736 if (bounds.first != bounds.second) {
737 bounds = std::equal_range(
738 d_stretchBendType.begin() + (bounds.first - d_kAtomType.begin()),
739 d_stretchBendType.begin() + (bounds.second - d_kAtomType.begin()),
740 canStretchBendType);
741 if (bounds.first != bounds.second) {
742 mmffStbnParams =
743 &d_params[bounds.first - d_stretchBendType.begin()];
744 }
745 }
746 }
747 }
748 #endif
749
750 return std::make_pair(swap, mmffStbnParams);
751 }
752
753 MMFFStbnCollection(std::string mmffStbn = "");
754 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
755 std::map<const unsigned int,
756 std::map<const unsigned int,
757 std::map<const unsigned int,
758 std::map<const unsigned int, MMFFStbn>>>>
759 d_params; //!< the parameter 4D-map
760 #else
761 std::vector<MMFFStbn> d_params; //!< the parameter vector
762 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
763 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
764 std::vector<std::uint8_t> d_kAtomType; //! atom type vector for atom k
765 std::vector<std::uint8_t>
766 d_stretchBendType; //! stretch-bend type vector for angle i-j-k
767 #endif
768 };
769
770 class RDKIT_FORCEFIELD_EXPORT MMFFDfsbCollection {
771 public:
772 //! Looks up the parameters for a particular key and returns them.
773 /*!
774 \return a pointer to the MMFFStbn object, NULL on failure.
775 */
getMMFFDfsbParams(const unsigned int periodicTableRow1,const unsigned int periodicTableRow2,const unsigned int periodicTableRow3)776 const std::pair<bool, const MMFFStbn *> getMMFFDfsbParams(
777 const unsigned int periodicTableRow1,
778 const unsigned int periodicTableRow2,
779 const unsigned int periodicTableRow3) const {
780 std::map<const unsigned int,
781 std::map<const unsigned int,
782 std::map<const unsigned int, MMFFStbn>>>::const_iterator
783 res1;
784 std::map<const unsigned int,
785 std::map<const unsigned int, MMFFStbn>>::const_iterator res2;
786 std::map<const unsigned int, MMFFStbn>::const_iterator res3;
787 const MMFFStbn *mmffDfsbParams = nullptr;
788 bool swap = false;
789 unsigned int canPeriodicTableRow1 = periodicTableRow1;
790 unsigned int canPeriodicTableRow3 = periodicTableRow3;
791 if (periodicTableRow1 > periodicTableRow3) {
792 canPeriodicTableRow1 = periodicTableRow3;
793 canPeriodicTableRow3 = periodicTableRow1;
794 swap = true;
795 }
796 res1 = d_params.find(canPeriodicTableRow1);
797 if (res1 != d_params.end()) {
798 res2 = ((*res1).second).find(periodicTableRow2);
799 if (res2 != ((*res1).second).end()) {
800 res3 = ((*res2).second).find(canPeriodicTableRow3);
801 if (res3 != ((*res2).second).end()) {
802 mmffDfsbParams = &((*res3).second);
803 }
804 }
805 }
806
807 return std::make_pair(swap, mmffDfsbParams);
808 }
809
810 MMFFDfsbCollection(std::string mmffDfsb = "");
811 std::map<const unsigned int,
812 std::map<const unsigned int, std::map<const unsigned int, MMFFStbn>>>
813 d_params; //!< the parameter 3D-map
814 };
815
816 class RDKIT_FORCEFIELD_EXPORT MMFFOopCollection {
817 public:
818 //! Looks up the parameters for a particular key and returns them.
819 /*!
820 \return a pointer to the MMFFOop object, NULL on failure.
821 */
operator()822 const MMFFOop *operator()(const MMFFDefCollection *mmffDef,
823 const unsigned int iAtomType,
824 const unsigned int jAtomType,
825 const unsigned int kAtomType,
826 const unsigned int lAtomType) const {
827 const MMFFOop *mmffOopParams = nullptr;
828 unsigned int iter = 0;
829 std::vector<unsigned int> canIKLAtomType(3);
830 // For out-of-plane bending ijk; I , where j is the central
831 // atom [cf. eq. (511, the five-stage protocol 1-1-1; 1, 2-2-2; 2,
832 // 3-2-3;3, 4-2-4;4, 5-2-5;5 is used. The final stage provides
833 // wild-card defaults for all except the central atom.
834 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
835 std::map<const unsigned int,
836 std::map<const unsigned int,
837 std::map<const unsigned int,
838 std::map<const unsigned int, MMFFOop>>>>::
839 const_iterator res1;
840 std::map<const unsigned int,
841 std::map<const unsigned int,
842 std::map<const unsigned int, MMFFOop>>>::const_iterator
843 res2;
844 std::map<const unsigned int,
845 std::map<const unsigned int, MMFFOop>>::const_iterator res3;
846 std::map<const unsigned int, MMFFOop>::const_iterator res4;
847 while ((iter < 4) && (!mmffOopParams)) {
848 canIKLAtomType[0] = (*mmffDef)(iAtomType)->eqLevel[iter];
849 unsigned int canJAtomType = jAtomType;
850 canIKLAtomType[1] = (*mmffDef)(kAtomType)->eqLevel[iter];
851 canIKLAtomType[2] = (*mmffDef)(lAtomType)->eqLevel[iter];
852 std::sort(canIKLAtomType.begin(), canIKLAtomType.end());
853 res1 = d_params.find(canIKLAtomType[0]);
854 if (res1 != d_params.end()) {
855 res2 = ((*res1).second).find(canJAtomType);
856 if (res2 != ((*res1).second).end()) {
857 res3 = ((*res2).second).find(canIKLAtomType[1]);
858 if (res3 != ((*res2).second).end()) {
859 res4 = ((*res3).second).find(canIKLAtomType[2]);
860 if (res4 != ((*res3).second).end()) {
861 mmffOopParams = &((*res4).second);
862 }
863 }
864 }
865 }
866 ++iter;
867 }
868 #else
869 std::pair<std::vector<std::uint8_t>::const_iterator,
870 std::vector<std::uint8_t>::const_iterator>
871 jBounds =
872 std::equal_range(d_jAtomType.begin(), d_jAtomType.end(), jAtomType);
873 std::pair<std::vector<std::uint8_t>::const_iterator,
874 std::vector<std::uint8_t>::const_iterator>
875 bounds;
876 if (jBounds.first != jBounds.second) {
877 while ((iter < 4) && (!mmffOopParams)) {
878 canIKLAtomType[0] = (*mmffDef)(iAtomType)->eqLevel[iter];
879 canIKLAtomType[1] = (*mmffDef)(kAtomType)->eqLevel[iter];
880 canIKLAtomType[2] = (*mmffDef)(lAtomType)->eqLevel[iter];
881 std::sort(canIKLAtomType.begin(), canIKLAtomType.end());
882 bounds = std::equal_range(
883 d_iAtomType.begin() + (jBounds.first - d_jAtomType.begin()),
884 d_iAtomType.begin() + (jBounds.second - d_jAtomType.begin()),
885 canIKLAtomType[0]);
886 if (bounds.first != bounds.second) {
887 bounds = std::equal_range(
888 d_kAtomType.begin() + (bounds.first - d_iAtomType.begin()),
889 d_kAtomType.begin() + (bounds.second - d_iAtomType.begin()),
890 canIKLAtomType[1]);
891 if (bounds.first != bounds.second) {
892 bounds = std::equal_range(
893 d_lAtomType.begin() + (bounds.first - d_kAtomType.begin()),
894 d_lAtomType.begin() + (bounds.second - d_kAtomType.begin()),
895 canIKLAtomType[2]);
896 if (bounds.first != bounds.second) {
897 mmffOopParams = &d_params[bounds.first - d_lAtomType.begin()];
898 }
899 }
900 }
901 ++iter;
902 }
903 }
904 #endif
905
906 return mmffOopParams;
907 }
908
909 MMFFOopCollection(const bool isMMFFs, std::string mmffOop = "");
910
911 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
912 std::map<const unsigned int,
913 std::map<const unsigned int,
914 std::map<const unsigned int,
915 std::map<const unsigned int, MMFFOop>>>>
916 d_params; //!< the parameter 4D-map
917 #else
918 std::vector<MMFFOop> d_params; //!< the parameter vector
919 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
920 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
921 std::vector<std::uint8_t> d_kAtomType; //! atom type vector for atom k
922 std::vector<std::uint8_t> d_lAtomType; //! atom type vector for atom l
923 #endif
924 };
925
926 class RDKIT_FORCEFIELD_EXPORT MMFFTorCollection {
927 public:
928 //! Looks up the parameters for a particular key and returns them.
929 /*!
930 \return a pointer to the MMFFTor object, NULL on failure.
931 */
getMMFFTorParams(const MMFFDefCollection * mmffDef,const std::pair<unsigned int,unsigned int> torType,const unsigned int iAtomType,const unsigned int jAtomType,const unsigned int kAtomType,const unsigned int lAtomType)932 const std::pair<const unsigned int, const MMFFTor *> getMMFFTorParams(
933 const MMFFDefCollection *mmffDef,
934 const std::pair<unsigned int, unsigned int> torType,
935 const unsigned int iAtomType, const unsigned int jAtomType,
936 const unsigned int kAtomType, const unsigned int lAtomType) const {
937 const MMFFTor *mmffTorParams = nullptr;
938 unsigned int iter = 0;
939 unsigned int iWildCard = 0;
940 unsigned int lWildCard = 0;
941 unsigned int canTorType = torType.first;
942 unsigned int maxIter = 5;
943 // For i-j-k-2 torsion interactions, a five-stage
944 // process based on level combinations 1-1-1-1, 2-2-2-2,
945 // 3-2-2-5, 5-2-2-3, and 5-2-2-5 is used, where stages 3
946 // and 4 correspond to "half-default" or "half-wild-card" entries.
947 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
948 std::map<
949 const unsigned int,
950 std::map<const unsigned int,
951 std::map<const unsigned int,
952 std::map<const unsigned int,
953 std::map<const unsigned int, MMFFTor>>>>>::
954 const_iterator res1;
955 std::map<const unsigned int,
956 std::map<const unsigned int,
957 std::map<const unsigned int,
958 std::map<const unsigned int, MMFFTor>>>>::
959 const_iterator res2;
960 std::map<const unsigned int,
961 std::map<const unsigned int,
962 std::map<const unsigned int, MMFFTor>>>::const_iterator
963 res3;
964 std::map<const unsigned int,
965 std::map<const unsigned int, MMFFTor>>::const_iterator res4;
966 std::map<const unsigned int, MMFFTor>::const_iterator res5;
967 #else
968 std::pair<std::vector<std::uint8_t>::const_iterator,
969 std::vector<std::uint8_t>::const_iterator>
970 jBounds;
971 std::pair<std::vector<std::uint8_t>::const_iterator,
972 std::vector<std::uint8_t>::const_iterator>
973 bounds;
974 #endif
975
976 while (((iter < maxIter) && ((!mmffTorParams) || (maxIter == 4))) ||
977 ((iter == 4) && (torType.first == 5) && torType.second)) {
978 // The rule of setting the torsion type to the value it had
979 // before being set to 5 as a last resort in case parameters
980 // could not be found is not mentioned in MMFF.IV; it was
981 // empirically discovered due to a number of tests in the
982 // MMFF validation suite otherwise failing
983 if ((maxIter == 5) && (iter == 4)) {
984 maxIter = 4;
985 iter = 0;
986 canTorType = torType.second;
987 }
988 iWildCard = iter;
989 lWildCard = iter;
990 if (iter == 1) {
991 iWildCard = 1;
992 lWildCard = 3;
993 } else if (iter == 2) {
994 iWildCard = 3;
995 lWildCard = 1;
996 }
997 unsigned int canIAtomType = (*mmffDef)(iAtomType)->eqLevel[iWildCard];
998 unsigned int canJAtomType = jAtomType;
999 unsigned int canKAtomType = kAtomType;
1000 unsigned int canLAtomType = (*mmffDef)(lAtomType)->eqLevel[lWildCard];
1001 if (canJAtomType > canKAtomType) {
1002 unsigned int temp = canKAtomType;
1003 canKAtomType = canJAtomType;
1004 canJAtomType = temp;
1005 temp = canLAtomType;
1006 canLAtomType = canIAtomType;
1007 canIAtomType = temp;
1008 } else if ((canJAtomType == canKAtomType) &&
1009 (canIAtomType > canLAtomType)) {
1010 unsigned int temp = canLAtomType;
1011 canLAtomType = canIAtomType;
1012 canIAtomType = temp;
1013 }
1014 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
1015 res1 = d_params.find(canTorType);
1016 if (res1 != d_params.end()) {
1017 res2 = ((*res1).second).find(canIAtomType);
1018 if (res2 != ((*res1).second).end()) {
1019 res3 = ((*res2).second).find(canJAtomType);
1020 if (res3 != ((*res2).second).end()) {
1021 res4 = ((*res3).second).find(canKAtomType);
1022 if (res4 != ((*res3).second).end()) {
1023 res5 = ((*res4).second).find(canLAtomType);
1024 if (res5 != ((*res4).second).end()) {
1025 mmffTorParams = &((*res5).second);
1026 if (maxIter == 4) {
1027 break;
1028 }
1029 }
1030 }
1031 }
1032 }
1033 }
1034 #else
1035 jBounds = std::equal_range(d_jAtomType.begin(), d_jAtomType.end(),
1036 canJAtomType);
1037 if (jBounds.first != jBounds.second) {
1038 bounds = std::equal_range(
1039 d_kAtomType.begin() + (jBounds.first - d_jAtomType.begin()),
1040 d_kAtomType.begin() + (jBounds.second - d_jAtomType.begin()),
1041 canKAtomType);
1042 if (bounds.first != bounds.second) {
1043 bounds = std::equal_range(
1044 d_iAtomType.begin() + (bounds.first - d_kAtomType.begin()),
1045 d_iAtomType.begin() + (bounds.second - d_kAtomType.begin()),
1046 canIAtomType);
1047 if (bounds.first != bounds.second) {
1048 bounds = std::equal_range(
1049 d_lAtomType.begin() + (bounds.first - d_iAtomType.begin()),
1050 d_lAtomType.begin() + (bounds.second - d_iAtomType.begin()),
1051 canLAtomType);
1052 if (bounds.first != bounds.second) {
1053 bounds = std::equal_range(
1054 d_torType.begin() + (bounds.first - d_lAtomType.begin()),
1055 d_torType.begin() + (bounds.second - d_lAtomType.begin()),
1056 canTorType);
1057 if (bounds.first != bounds.second) {
1058 mmffTorParams = &d_params[bounds.first - d_torType.begin()];
1059 if (maxIter == 4) {
1060 break;
1061 }
1062 }
1063 }
1064 }
1065 }
1066 }
1067 #endif
1068 ++iter;
1069 }
1070
1071 return std::make_pair(canTorType, mmffTorParams);
1072 }
1073
1074 MMFFTorCollection(const bool isMMFFs, std::string mmffTor = "");
1075 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
1076 std::map<
1077 const unsigned int,
1078 std::map<
1079 const unsigned int,
1080 std::map<const unsigned int,
1081 std::map<const unsigned int, std::map<const unsigned int,
1082 MMFFTor>>>>>
1083 d_params; //!< the parameter 5D-map
1084 #else
1085 std::vector<MMFFTor> d_params; //!< the parameter vector
1086 std::vector<std::uint8_t> d_iAtomType; //! atom type vector for atom i
1087 std::vector<std::uint8_t> d_jAtomType; //! atom type vector for atom j
1088 std::vector<std::uint8_t> d_kAtomType; //! atom type vector for atom k
1089 std::vector<std::uint8_t> d_lAtomType; //! atom type vector for atom l
1090 std::vector<std::uint8_t>
1091 d_torType; //! torsion type vector for angle i-j-k-l
1092 #endif
1093 };
1094
1095 class RDKIT_FORCEFIELD_EXPORT MMFFVdWCollection {
1096 public:
1097 //! gets a pointer to the singleton MMFFVdWCollection
1098 double power;
1099 double B;
1100 double Beta;
1101 double DARAD;
1102 double DAEPS;
1103 //! Looks up the parameters for a particular key and returns them.
1104 /*!
1105 \return a pointer to the MMFFVdW object, NULL on failure.
1106 */
operator()1107 const MMFFVdW *operator()(const unsigned int atomType) const {
1108 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
1109 std::map<const unsigned int, MMFFVdW>::const_iterator res;
1110 res = d_params.find(atomType);
1111
1112 return (res != d_params.end() ? &((*res).second) : NULL);
1113 #else
1114 std::pair<std::vector<std::uint8_t>::const_iterator,
1115 std::vector<std::uint8_t>::const_iterator>
1116 bounds =
1117 std::equal_range(d_atomType.begin(), d_atomType.end(), atomType);
1118
1119 return ((bounds.first != bounds.second)
1120 ? &d_params[bounds.first - d_atomType.begin()]
1121 : nullptr);
1122 #endif
1123 }
1124
1125 MMFFVdWCollection(std::string mmffVdW = "");
1126 #ifdef RDKIT_MMFF_PARAMS_USE_STD_MAP
1127 std::map<const unsigned int, MMFFVdW> d_params; //!< the parameter map
1128 #else
1129 std::vector<MMFFVdW> d_params; //!< the parameter vector
1130 std::vector<std::uint8_t> d_atomType; //! atom type vector
1131 #endif
1132 };
1133 } // namespace MMFF
1134 } // namespace ForceFields
1135
1136 #endif
1137