1 /*********************************************************************
2 forcefieldmmff94.cpp - MMFF94 force field
3 
4 Copyright (C) 2006-2008 by Tim Vandermeersch <tim.vandermeersch@gmail.com>
5 
6 This file is part of the Open Babel project.
7 For more information, see <http://openbabel.org/>
8 
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation version 2 of the License.
12 
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 GNU General Public License for more details.
17 ***********************************************************************/
18 
19 /*
20  * Source code layout:
21  * - Functions to calculate the actual interactions
22  * - Parse parameter files
23  * - Setup Functions
24  * - Validation functions
25  * - Calculate bond type, angle type, stretch-bend type, torsion type
26  * - Various tables & misc. functions
27  *
28  */
29 
30 #include <openbabel/babelconfig.h>
31 #include <openbabel/obconversion.h>
32 #include <openbabel/mol.h>
33 #include <openbabel/locale.h>
34 #include <openbabel/elements.h>
35 #include <openbabel/bond.h>
36 #include <openbabel/obiter.h>
37 #include <openbabel/ring.h>
38 
39 #ifndef M_PI
40 #define M_PI 3.14159265358979323846
41 #endif
42 
43 #include <iomanip>
44 #include "forcefieldmmff94.h"
45 
46 using namespace std;
47 
48 namespace OpenBabel
49 {
50   ////////////////////////////////////////////////////////////////////////////////
51   ////////////////////////////////////////////////////////////////////////////////
52   //
53   //  Functions to calculate the actual interactions
54   //
55   ////////////////////////////////////////////////////////////////////////////////
56   ////////////////////////////////////////////////////////////////////////////////
57 
Energy(bool gradients)58   double OBForceFieldMMFF94::Energy(bool gradients)
59   {
60     double energy;
61 
62     IF_OBFF_LOGLVL_MEDIUM
63       OBFFLog("\nE N E R G Y\n\n");
64 
65     if (gradients) {
66       ClearGradients();
67       energy  = E_Bond<true>();
68       energy += E_Angle<true>();
69       energy += E_StrBnd<true>();
70       energy += E_Torsion<true>();
71       energy += E_OOP<true>();
72       energy += E_VDW<true>();
73       energy += E_Electrostatic<true>();
74     } else {
75       energy  = E_Bond<false>();
76       energy += E_Angle<false>();
77       energy += E_StrBnd<false>();
78       energy += E_Torsion<false>();
79       energy += E_OOP<false>();
80       energy += E_VDW<false>();
81       energy += E_Electrostatic<false>();
82     }
83 
84     IF_OBFF_LOGLVL_MEDIUM {
85       snprintf(_logbuf, BUFF_SIZE, "\nTOTAL ENERGY = %8.5f %s\n", energy, GetUnit().c_str());
86       OBFFLog(_logbuf);
87     }
88 
89     return energy;
90   }
91 
92   //
93   // MMFF part I - page 494
94   //
95   //                   kb_ij                              7
96   // EB_ij = 143.9325 ------- /\r_ij^2 (1 + cs /\_rij + ---- cs^2 r_ij^2)
97   //                     2                               12
98   //
99   // kb_ij	force constant (md/A)
100   //
101   // /\r_ij 	r_ij - r0_ij (A)
102   //
103   // cs		cubic stretch constant = -2 A^(-1)
104   //
105   template<bool gradients>
Compute()106   inline void OBFFBondCalculationMMFF94::Compute()
107   {
108     if (OBForceField::IgnoreCalculation(idx_a, idx_b)) {
109       energy = 0.0;
110       return;
111     }
112 
113     double delta2;
114 
115     if (gradients) {
116       rab = OBForceField::VectorBondDerivative(pos_a, pos_b, force_a, force_b);
117       delta = rab - r0;
118       delta2 = delta * delta;
119 
120       const double dE = 143.9325 * kb * delta * (1.0 - 3.0 * delta + 14.0/3.0 * delta2);
121 
122       OBForceField::VectorSelfMultiply(force_a, dE);
123       OBForceField::VectorSelfMultiply(force_b, dE);
124     } else {
125       rab = OBForceField::VectorDistance(pos_a, pos_b);
126       delta = rab - r0;
127       delta2 = delta * delta;
128     }
129 
130     energy = kb * delta2 * (1.0 - 2.0 * delta + 7.0/3.0 * delta2);
131   }
132 
133   template<bool gradients>
E_Bond()134   double OBForceFieldMMFF94::E_Bond()
135   {
136     double energy = 0.0;
137 
138     IF_OBFF_LOGLVL_HIGH {
139       OBFFLog("\nB O N D   S T R E T C H I N G\n\n");
140       OBFFLog("ATOM TYPES   FF    BOND       IDEAL       FORCE\n");
141       OBFFLog(" I    J     CLASS  LENGTH     LENGTH     CONSTANT      DELTA      ENERGY\n");
142       OBFFLog("------------------------------------------------------------------------\n");
143     }
144 
145     #ifdef _OPENMP
146     #pragma omp parallel for reduction(+:energy)
147     #endif
148     for (int i = 0; i < _bondcalculations.size(); ++i) {
149       _bondcalculations[i].template Compute<gradients>();
150       energy += _bondcalculations[i].energy;
151 
152       #ifndef _OPENMP
153       if (gradients) {
154         AddGradient(_bondcalculations[i].force_a, _bondcalculations[i].idx_a);
155         AddGradient(_bondcalculations[i].force_b, _bondcalculations[i].idx_b);
156       }
157       #endif
158 
159       IF_OBFF_LOGLVL_HIGH {
160         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d      %d   %8.3f   %8.3f     %8.3f   %8.3f   %8.3f\n",
161                 atoi(_bondcalculations[i].a->GetType()), atoi(_bondcalculations[i].b->GetType()),
162                 _bondcalculations[i].bt, _bondcalculations[i].rab, _bondcalculations[i].r0,
163                 _bondcalculations[i].kb, _bondcalculations[i].delta,
164                 143.9325 * 0.5 * _bondcalculations[i].energy);
165         OBFFLog(_logbuf);
166       }
167     }
168 
169     #ifdef _OPENMP
170     for (int i = 0; i < _bondcalculations.size(); ++i) {
171       if (gradients) {
172         AddGradient(_bondcalculations[i].force_a, _bondcalculations[i].idx_a);
173         AddGradient(_bondcalculations[i].force_b, _bondcalculations[i].idx_b);
174       }
175     }
176     #endif
177 
178     IF_OBFF_LOGLVL_MEDIUM {
179       snprintf(_logbuf, BUFF_SIZE, "     TOTAL BOND STRETCHING ENERGY = %8.5f %s\n",  143.9325 * 0.5 * energy, GetUnit().c_str());
180       OBFFLog(_logbuf);
181     }
182 
183     return (143.9325 * 0.5 * energy);
184   }
185 
186   //
187   // MMFF part I - page 495
188   //
189   //                       ka_ijk
190   // EA_ijk = 0.438449325 -------- /\0_ijk^2 (1 + cs /\0_ijk)
191   //                         2
192   //
193   // ka_ijk	force constant (md A/rad^2)
194   //
195   // /\0_ijk 	0_ijk - 00_ijk (degrees)
196   //
197   // cs		cubic bend constant = -0.007 deg^-1 = -0.4 rad^-1
198   //
199   template<bool gradients>
Compute()200   inline void OBFFAngleCalculationMMFF94::Compute()
201   {
202     if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c)) {
203       energy = 0.0;
204       return;
205     }
206 
207     double delta2, dE;
208 
209     if (gradients) {
210       theta = OBForceField::VectorAngleDerivative(pos_a, pos_b, pos_c, force_a, force_b, force_c);
211 
212       if (!isfinite(theta))
213         theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us;
214 
215       delta = theta - theta0;
216 
217       if (linear) {
218         energy = 143.9325 * ka * (1.0 + cos((theta) * DEG_TO_RAD));
219         dE = -sin((theta) * DEG_TO_RAD) * 143.9325 * ka;
220       } else {
221         delta2 = delta * delta;
222         energy = 0.043844 * 0.5 * ka * delta2 * (1.0 - 0.007 * delta);
223         dE = RAD_TO_DEG * 0.043844 * ka * delta * (1.0 - 1.5 * 0.007 * delta);
224       }
225 
226       OBForceField::VectorSelfMultiply(force_a, dE);
227       OBForceField::VectorSelfMultiply(force_b, dE);
228       OBForceField::VectorSelfMultiply(force_c, dE);
229     } else {
230       theta = OBForceField::VectorAngle(pos_a, pos_b, pos_c);
231 
232       if (!isfinite(theta))
233         theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us;
234 
235       delta = theta - theta0;
236 
237       if (linear) {
238         energy = 143.9325 * ka * (1.0 + cos(theta * DEG_TO_RAD));
239       } else {
240         delta2 = delta * delta;
241         energy = 0.043844 * 0.5 * ka * delta2 * (1.0 - 0.007 * delta);
242       }
243 
244     }
245 
246   }
247 
248   template<bool gradients>
E_Angle()249   double OBForceFieldMMFF94::E_Angle()
250   {
251     double energy = 0.0;
252 
253     IF_OBFF_LOGLVL_HIGH {
254       OBFFLog("\nA N G L E   B E N D I N G\n\n");
255       OBFFLog("ATOM TYPES        FF    VALENCE     IDEAL      FORCE\n");
256       OBFFLog(" I    J    K     CLASS   ANGLE      ANGLE     CONSTANT      DELTA      ENERGY\n");
257       OBFFLog("-----------------------------------------------------------------------------\n");
258     }
259 
260     #ifdef _OPENMP
261     #pragma omp parallel for reduction(+:energy)
262     #endif
263     for (int i = 0; i < _anglecalculations.size(); ++i) {
264 
265       _anglecalculations[i].template Compute<gradients>();
266       energy += _anglecalculations[i].energy;
267 
268       #ifndef _OPENMP
269       if (gradients) {
270         AddGradient(_anglecalculations[i].force_a, _anglecalculations[i].idx_a);
271         AddGradient(_anglecalculations[i].force_b, _anglecalculations[i].idx_b);
272         AddGradient(_anglecalculations[i].force_c, _anglecalculations[i].idx_c);
273       }
274       #endif
275 
276       IF_OBFF_LOGLVL_HIGH {
277         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d   %2d      %d   %8.3f   %8.3f     %8.3f   %8.3f   %8.3f\n",
278                 atoi(_anglecalculations[i].a->GetType()), atoi(_anglecalculations[i].b->GetType()),
279                 atoi(_anglecalculations[i].c->GetType()), _anglecalculations[i].at,
280                 _anglecalculations[i].theta, _anglecalculations[i].theta0,
281                 _anglecalculations[i].ka, _anglecalculations[i].delta,
282                 _anglecalculations[i].energy);
283         OBFFLog(_logbuf);
284       }
285     }
286 
287     #ifdef _OPENMP
288     for (int i = 0; i < _anglecalculations.size(); ++i) {
289       if (gradients) {
290         AddGradient(_anglecalculations[i].force_a, _anglecalculations[i].idx_a);
291         AddGradient(_anglecalculations[i].force_b, _anglecalculations[i].idx_b);
292         AddGradient(_anglecalculations[i].force_c, _anglecalculations[i].idx_c);
293       }
294     }
295     #endif
296 
297     IF_OBFF_LOGLVL_MEDIUM {
298       snprintf(_logbuf, BUFF_SIZE, "     TOTAL ANGLE BENDING ENERGY = %8.5f %s\n", energy, GetUnit().c_str());
299       OBFFLog(_logbuf);
300     }
301 
302     return energy;
303   }
304 
305   //
306   // MMFF part I - page 495
307   //
308   // EBA_ijk = 2.51210 (kba_ijk /\r_ij + kba_kji /\r_kj) /\0_ijk
309   //
310   // kba_ijk	force constant (md/rad)
311   // kba_kji	force constant (md/rad)
312   //
313   // /\r_xx 	see above
314   // /\0_ijk 	see above
315   //
316   template<bool gradients>
Compute()317   inline void OBFFStrBndCalculationMMFF94::Compute()
318   {
319     if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c)) {
320       energy = 0.0;
321       return;
322     }
323 
324     if (gradients) {
325       theta = OBForceField::VectorAngleDerivative(pos_a, pos_b, pos_c,
326                                                   force_abc_a, force_abc_b, force_abc_c);
327       rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_ab_a, force_ab_b);
328       rbc = OBForceField::VectorDistanceDerivative(pos_b, pos_c, force_bc_b, force_bc_c);
329     } else {
330       theta = OBForceField::VectorAngle(pos_a, pos_b, pos_c);
331       rab = OBForceField::VectorDistance(pos_a, pos_b);
332       rbc = OBForceField::VectorDistance(pos_b, pos_c);
333     }
334 
335     if (!isfinite(theta))
336       theta = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us;
337 
338     delta_theta = theta - theta0;
339     delta_rab = rab - rab0;
340     delta_rbc = rbc - rbc0;
341     const double factor = RAD_TO_DEG * (kbaABC * delta_rab + kbaCBA * delta_rbc);
342 
343     energy = DEG_TO_RAD * factor * delta_theta;
344     if (gradients) {
345       //grada = 2.51210 * (kbaABC * rab_da * delta_theta + RAD_TO_DEG * theta_da * (kbaABC * delta_rab + kbaCBA * delta_rbc));
346       OBForceField::VectorSelfMultiply(force_ab_a, (kbaABC*delta_theta));
347       OBForceField::VectorSelfMultiply(force_abc_a, factor);
348       OBForceField::VectorAdd(force_ab_a, force_abc_a, force_ab_a);
349       OBForceField::VectorMultiply(force_ab_a, 2.51210, force_a);
350       //gradc = 2.51210 * (kbaCBA * rbc_dc * delta_theta + RAD_TO_DEG * theta_dc * (kbaABC * delta_rab + kbaCBA * delta_rbc));
351       OBForceField::VectorSelfMultiply(force_bc_c, (kbaCBA*delta_theta));
352       OBForceField::VectorSelfMultiply(force_abc_c, factor);
353       OBForceField::VectorAdd(force_bc_c, force_abc_c, force_bc_c);
354       OBForceField::VectorMultiply(force_bc_c, 2.51210, force_c);
355       //gradb = -grada - gradc;
356       OBForceField::VectorAdd(force_a, force_c, force_b);
357       OBForceField::VectorSelfMultiply(force_b, -1.0);
358     }
359   }
360 
361   template<bool gradients>
E_StrBnd()362   double OBForceFieldMMFF94::E_StrBnd()
363   {
364     double energy = 0.0;
365 
366     IF_OBFF_LOGLVL_HIGH {
367       OBFFLog("\nS T R E T C H   B E N D I N G\n\n");
368       OBFFLog("ATOM TYPES        FF    VALENCE     DELTA        FORCE CONSTANT\n");
369       OBFFLog(" I    J    K     CLASS   ANGLE      ANGLE        I J        J K      ENERGY\n");
370       OBFFLog("---------------------------------------------------------------------------\n");
371     }
372 
373     #ifdef _OPENMP
374     #pragma omp parallel for reduction(+:energy)
375     #endif
376     for (int i = 0; i < _strbndcalculations.size(); ++i) {
377 
378       _strbndcalculations[i].template Compute<gradients>();
379       energy += _strbndcalculations[i].energy;
380 
381       #ifndef _OPENMP
382       if (gradients) {
383         AddGradient(_strbndcalculations[i].force_a, _strbndcalculations[i].idx_a);
384         AddGradient(_strbndcalculations[i].force_b, _strbndcalculations[i].idx_b);
385         AddGradient(_strbndcalculations[i].force_c, _strbndcalculations[i].idx_c);
386       }
387       #endif
388 
389       IF_OBFF_LOGLVL_HIGH {
390         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d   %2d     %2d   %8.3f   %8.3f   %8.3f   %8.3f   %8.3f\n",
391                 atoi(_strbndcalculations[i].a->GetType()), atoi(_strbndcalculations[i].b->GetType()),
392                 atoi(_strbndcalculations[i].c->GetType()), _strbndcalculations[i].sbt,
393                 _strbndcalculations[i].theta, _strbndcalculations[i].delta_theta,
394                 _strbndcalculations[i].kbaABC, _strbndcalculations[i].kbaCBA,
395                 2.51210 * _strbndcalculations[i].energy);
396         OBFFLog(_logbuf);
397       }
398     }
399 
400     #ifdef _OPENMP
401     for (int i = 0; i < _strbndcalculations.size(); ++i) {
402       if (gradients) {
403         AddGradient(_strbndcalculations[i].force_a, _strbndcalculations[i].idx_a);
404         AddGradient(_strbndcalculations[i].force_b, _strbndcalculations[i].idx_b);
405         AddGradient(_strbndcalculations[i].force_c, _strbndcalculations[i].idx_c);
406       }
407     }
408     #endif
409 
410     IF_OBFF_LOGLVL_MEDIUM {
411       snprintf(_logbuf, BUFF_SIZE, "     TOTAL STRETCH BENDING ENERGY = %8.5f %s\n", 2.51210 * energy, GetUnit().c_str());
412       OBFFLog(_logbuf);
413     }
414 
415     return (2.51210 * energy);
416   }
417 
GetElementRow(OBAtom * atom)418   int OBForceFieldMMFF94::GetElementRow(OBAtom *atom)
419   {
420     int row;
421 
422     row = 0;
423 
424     if (atom->GetAtomicNum() > 2)
425       row++;
426     if (atom->GetAtomicNum() > 10)
427       row++;
428     if (atom->GetAtomicNum() > 18)
429       row++;
430     if (atom->GetAtomicNum() > 36)
431       row++;
432     if (atom->GetAtomicNum() > 54)
433       row++;
434     if (atom->GetAtomicNum() > 86)
435       row++;
436 
437     return row;
438   }
439 
440   //
441   // MMFF part I - page 495
442   //
443   // ET_ijkl = 0.5 ( V1 (1 + cos(0_ijkl)) + V2 (1 - cos(2 0_ijkl)) + V3 (1 + cos(3 0_ijkl)) )
444   //
445   // V1		force constant (md/rad)
446   // V2		force constant (md/rad)
447   // V3		force constant (md/rad)
448   //
449   // 0_ijkl 	torsion angle (degrees)
450   //
451   template<bool gradients>
Compute()452   inline void OBFFTorsionCalculationMMFF94::Compute()
453   {
454     if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c, idx_d)) {
455       energy = 0.0;
456       return;
457     }
458 
459     double cosine, cosine2, cosine3;
460     double phi1, phi2, phi3;
461     double dE, sine, sine2, sine3;
462 
463     if (gradients) {
464       tor = OBForceField::VectorTorsionDerivative(pos_a, pos_b, pos_c, pos_d,
465                                                   force_a, force_b, force_c, force_d);
466       if (!isfinite(tor))
467         tor = 1.0e-3;
468 
469       sine = sin(DEG_TO_RAD * tor);
470       sine2 = sin(2.0 * DEG_TO_RAD * tor);
471       sine3 = sin(3.0 * DEG_TO_RAD * tor);
472 
473       dE = 0.5 * (v1 * sine - 2.0 * v2 * sine2 + 3.0 * v3 * sine3); // MMFF
474 
475       OBForceField::VectorSelfMultiply(force_a, dE);
476       OBForceField::VectorSelfMultiply(force_b, dE);
477       OBForceField::VectorSelfMultiply(force_c, dE);
478       OBForceField::VectorSelfMultiply(force_d, dE);
479     } else {
480       tor = OBForceField::VectorTorsion(pos_a, pos_b, pos_c, pos_d);
481       if (!isfinite(tor))
482         tor = 1.0e-3;
483     }
484 
485     cosine = cos(DEG_TO_RAD * tor);
486     cosine2 = cos(DEG_TO_RAD * 2 * tor);
487     cosine3 = cos(DEG_TO_RAD * 3 * tor);
488 
489     phi1 = 1.0 + cosine;
490     phi2 = 1.0 - cosine2;
491     phi3 = 1.0 + cosine3;
492 
493     energy = (v1 * phi1 + v2 * phi2 + v3 * phi3);
494 
495   }
496 
497   template<bool gradients>
E_Torsion()498   double OBForceFieldMMFF94::E_Torsion()
499   {
500     double energy = 0.0;
501 
502     IF_OBFF_LOGLVL_HIGH {
503       OBFFLog("\nT O R S I O N A L\n\n");
504       OBFFLog("ATOM TYPES             FF     TORSION       FORCE CONSTANT\n");
505       OBFFLog(" I    J    K    L     CLASS    ANGLE         V1   V2   V3     ENERGY\n");
506       OBFFLog("--------------------------------------------------------------------\n");
507     }
508 
509     #ifdef _OPENMP
510     #pragma omp parallel for reduction(+:energy)
511     #endif
512     for (int i = 0; i < _torsioncalculations.size(); ++i) {
513 
514       _torsioncalculations[i].template Compute<gradients>();
515       energy += _torsioncalculations[i].energy;
516 
517       #ifndef _OPENMP
518       if (gradients) {
519         AddGradient(_torsioncalculations[i].force_a, _torsioncalculations[i].idx_a);
520         AddGradient(_torsioncalculations[i].force_b, _torsioncalculations[i].idx_b);
521         AddGradient(_torsioncalculations[i].force_c, _torsioncalculations[i].idx_c);
522         AddGradient(_torsioncalculations[i].force_d, _torsioncalculations[i].idx_d);
523       }
524       #endif
525 
526       IF_OBFF_LOGLVL_HIGH {
527         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d   %2d   %2d      %d   %8.3f   %6.3f   %6.3f   %6.3f   %8.3f\n",
528                 atoi(_torsioncalculations[i].a->GetType()), atoi(_torsioncalculations[i].b->GetType()),
529                 atoi(_torsioncalculations[i].c->GetType()), atoi(_torsioncalculations[i].d->GetType()),
530                 _torsioncalculations[i].tt, _torsioncalculations[i].tor, _torsioncalculations[i].v1,
531                 _torsioncalculations[i].v2, _torsioncalculations[i].v3, 0.5 * _torsioncalculations[i].energy);
532         OBFFLog(_logbuf);
533       }
534 
535     }
536 
537     #ifdef _OPENMP
538     for (int i = 0; i < _torsioncalculations.size(); ++i) {
539       if (gradients) {
540         AddGradient(_torsioncalculations[i].force_a, _torsioncalculations[i].idx_a);
541         AddGradient(_torsioncalculations[i].force_b, _torsioncalculations[i].idx_b);
542         AddGradient(_torsioncalculations[i].force_c, _torsioncalculations[i].idx_c);
543         AddGradient(_torsioncalculations[i].force_d, _torsioncalculations[i].idx_d);
544       }
545     }
546     #endif
547 
548     IF_OBFF_LOGLVL_MEDIUM {
549       snprintf(_logbuf, BUFF_SIZE, "     TOTAL TORSIONAL ENERGY = %8.5f %s\n", 0.5 * energy, GetUnit().c_str());
550       OBFFLog(_logbuf);
551     }
552 
553     return (0.5 * energy);
554   }
555 
556   //						//
557   //  a						//
558   //   \  					//
559   //    b---d      plane = a-b-c		//
560   //   / 					//
561   //  c						//
562   //						//
563   template<bool gradients>
Compute()564   void OBFFOOPCalculationMMFF94::Compute()
565   {
566     if (OBForceField::IgnoreCalculation(idx_a, idx_b, idx_c, idx_d)) {
567       energy = 0.0;
568       return;
569     }
570 
571     double angle2, dE;
572 
573     if (gradients) {
574       angle = OBForceField::VectorOOPDerivative(pos_a, pos_b, pos_c, pos_d,
575                                                 force_a, force_b, force_c, force_d);
576 
577       dE =  (-1.0 * RAD_TO_DEG * 0.043844 * angle * koop) / cos(angle * DEG_TO_RAD);
578 
579       OBForceField::VectorSelfMultiply(force_a, dE);
580       OBForceField::VectorSelfMultiply(force_b, dE);
581       OBForceField::VectorSelfMultiply(force_c, dE);
582       OBForceField::VectorSelfMultiply(force_d, dE);
583     } else {
584       angle = OBForceField::VectorOOP(pos_a, pos_b, pos_c, pos_d);
585     }
586 
587     if (!isfinite(angle))
588       angle = 0.0; // doesn't explain why GetAngle is returning NaN but solves it for us;
589 
590     angle2 = angle * angle;
591     energy = koop * angle2;
592 
593   }
594 
595   template<bool gradients>
E_OOP()596   double OBForceFieldMMFF94::E_OOP()
597   {
598     double energy = 0.0;
599 
600     IF_OBFF_LOGLVL_HIGH {
601       OBFFLog("\nO U T - O F - P L A N E   B E N D I N G\n\n");
602       OBFFLog("ATOM TYPES             FF       OOP     FORCE\n");
603       OBFFLog(" I    J    K    L     CLASS    ANGLE   CONSTANT     ENERGY\n");
604       OBFFLog("----------------------------------------------------------\n");
605     }
606 
607     #ifdef _OPENMP
608     #pragma omp parallel for reduction(+:energy)
609     #endif
610     for (int i = 0; i < _oopcalculations.size(); ++i) {
611 
612       _oopcalculations[i].template Compute<gradients>();
613       energy += _oopcalculations[i].energy;
614 
615       #ifndef _OPENMP
616       if (gradients) {
617         AddGradient(_oopcalculations[i].force_a, _oopcalculations[i].idx_a);
618         AddGradient(_oopcalculations[i].force_b, _oopcalculations[i].idx_b);
619         AddGradient(_oopcalculations[i].force_c, _oopcalculations[i].idx_c);
620         AddGradient(_oopcalculations[i].force_d, _oopcalculations[i].idx_d);
621       }
622       #endif
623 
624       IF_OBFF_LOGLVL_HIGH {
625         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d   %2d   %2d      0   %8.3f   %8.3f     %8.3f\n",
626                 atoi(_oopcalculations[i].a->GetType()), atoi(_oopcalculations[i].b->GetType()),
627                 atoi(_oopcalculations[i].c->GetType()), atoi(_oopcalculations[i].d->GetType()),
628                 _oopcalculations[i].angle, _oopcalculations[i].koop,
629                 0.043844 * 0.5 * _oopcalculations[i].energy);
630         OBFFLog(_logbuf);
631       }
632     }
633 
634     #ifdef _OPENMP
635     for (int i = 0; i < _oopcalculations.size(); ++i) {
636       if (gradients) {
637         AddGradient(_oopcalculations[i].force_a, _oopcalculations[i].idx_a);
638         AddGradient(_oopcalculations[i].force_b, _oopcalculations[i].idx_b);
639         AddGradient(_oopcalculations[i].force_c, _oopcalculations[i].idx_c);
640         AddGradient(_oopcalculations[i].force_d, _oopcalculations[i].idx_d);
641       }
642     }
643     #endif
644 
645     IF_OBFF_LOGLVL_MEDIUM {
646       snprintf(_logbuf, BUFF_SIZE, "     TOTAL OUT-OF-PLANE BENDING ENERGY = %8.5f %s\n", 0.043844 * 0.5 * energy, GetUnit().c_str());
647       OBFFLog(_logbuf);
648     }
649 
650     return (0.043844 * 0.5 * energy);
651   }
652 
653   template<bool gradients>
Compute()654   inline void OBFFVDWCalculationMMFF94::Compute()
655   {
656     if (OBForceField::IgnoreCalculation(idx_a, idx_b)) {
657       energy = 0.0;
658       return;
659     }
660 
661     if (gradients) {
662       rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_a, force_b);
663     } else {
664       rab = OBForceField::VectorDistance(pos_a, pos_b);
665     }
666 
667     const double rab7 = rab*rab*rab*rab*rab*rab*rab;
668 
669     double erep = (1.07 * R_AB) / (rab + 0.07 * R_AB); //***
670     double erep7 = erep*erep*erep*erep*erep*erep*erep;
671 
672     double eattr = (((1.12 * R_AB7) / (rab7 + 0.12 * R_AB7)) - 2.0);
673 
674     energy = epsilon * erep7 * eattr;
675 
676     if (gradients) {
677       const double q = rab / R_AB;
678       const double q6 = q*q*q*q*q*q;
679       const double q7 = q6 * q;
680       erep = 1.07 / (q + 0.07);
681       erep7 = erep*erep*erep*erep*erep*erep*erep;
682       const double term = q7 + 0.12;
683       const double term2 = term * term;
684       eattr = (-7.84 * q6) / term2 + ((-7.84 / term) + 14) / (q + 0.07);
685       const double dE = (epsilon / R_AB) * erep7 * eattr;
686       OBForceField::VectorSelfMultiply(force_a, dE);
687       OBForceField::VectorSelfMultiply(force_b, dE);
688     }
689   }
690 
691   template<bool gradients>
E_VDW()692   double OBForceFieldMMFF94::E_VDW()
693   {
694     double energy = 0.0;
695 
696     IF_OBFF_LOGLVL_HIGH {
697       OBFFLog("\nV A N   D E R   W A A L S\n\n");
698       OBFFLog("ATOM TYPES\n");
699       OBFFLog(" I    J        Rij       R*IJ    EPSILON    ENERGY\n");
700       OBFFLog("--------------------------------------------------\n");
701       //       XX   XX     -000.000  -000.000  -000.000  -000.000
702     }
703 
704     #ifdef _OPENMP
705     #pragma omp parallel for reduction(+:energy)
706     #endif
707     for (int i = 0; i < _vdwcalculations.size(); ++i) {
708       // Cut-off check
709       if (_cutoff)
710         if (!_vdwpairs.BitIsSet(_vdwcalculations[i].pairIndex))
711           continue;
712 
713       _vdwcalculations[i].template Compute<gradients>();
714       energy += _vdwcalculations[i].energy;
715 
716       #ifndef _OPENMP
717       if (gradients) {
718         AddGradient(_vdwcalculations[i].force_a, _vdwcalculations[i].idx_a);
719         AddGradient(_vdwcalculations[i].force_b, _vdwcalculations[i].idx_b);
720       }
721       #endif
722 
723       IF_OBFF_LOGLVL_HIGH {
724         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d     %8.3f  %8.3f  %8.3f  %8.3f\n",
725                 atoi(_vdwcalculations[i].a->GetType()), atoi(_vdwcalculations[i].b->GetType()),
726                 _vdwcalculations[i].rab, _vdwcalculations[i].R_AB, _vdwcalculations[i].epsilon, _vdwcalculations[i].energy);
727         OBFFLog(_logbuf);
728       }
729 
730     }
731 
732     #ifdef _OPENMP
733     for (int i = 0; i < _vdwcalculations.size(); ++i) {
734       // Cut-off check
735       if (_cutoff)
736         if (!_vdwpairs.BitIsSet(i))
737           continue;
738 
739       if (gradients) {
740         AddGradient(_vdwcalculations[i].force_a, _vdwcalculations[i].idx_a);
741         AddGradient(_vdwcalculations[i].force_b, _vdwcalculations[i].idx_b);
742       }
743     }
744     #endif
745 
746     IF_OBFF_LOGLVL_MEDIUM {
747       snprintf(_logbuf, BUFF_SIZE, "     TOTAL VAN DER WAALS ENERGY = %8.5f %s\n", energy, GetUnit().c_str());
748       OBFFLog(_logbuf);
749     }
750 
751     return energy;
752   }
753 
754   template<bool gradients>
Compute()755   inline void OBFFElectrostaticCalculationMMFF94::Compute()
756   {
757     if (OBForceField::IgnoreCalculation(idx_a, idx_b)) {
758       energy = 0.0;
759       return;
760     }
761 
762     if (gradients) {
763       rab = OBForceField::VectorDistanceDerivative(pos_a, pos_b, force_a, force_b);
764       rab += 0.05; // ??
765       const double rab2 = rab * rab;
766       const double dE = -qq / rab2;
767       OBForceField::VectorSelfMultiply(force_a, dE);
768       OBForceField::VectorSelfMultiply(force_b, dE);
769     } else {
770       rab = OBForceField::VectorDistance(pos_a, pos_b);
771       rab += 0.05; // ??
772     }
773 
774     energy = qq / rab;
775   }
776 
777   template<bool gradients>
E_Electrostatic()778   double OBForceFieldMMFF94::E_Electrostatic()
779   {
780     double energy = 0.0;
781 
782     IF_OBFF_LOGLVL_HIGH {
783       OBFFLog("\nE L E C T R O S T A T I C   I N T E R A C T I O N S\n\n");
784       OBFFLog("ATOM TYPES\n");
785       OBFFLog(" I    J        Rij        Qi         Qj        ENERGY\n");
786       OBFFLog("-----------------------------------------------------\n");
787       //       XX   XX     XXXXXXXX   XXXXXXXX   XXXXXXXX   XXXXXXXX
788     }
789 
790     #ifdef _OPENMP
791     #pragma omp parallel for reduction(+:energy)
792     #endif
793     for (int i = 0; i < _electrostaticcalculations.size(); ++i) {
794       // Cut-off check
795       if (_cutoff)
796         if (!_elepairs.BitIsSet(_electrostaticcalculations[i].pairIndex))
797           continue;
798 
799       _electrostaticcalculations[i].template Compute<gradients>();
800       energy += _electrostaticcalculations[i].energy;
801 
802       #ifndef _OPENMP
803       if (gradients) {
804         AddGradient(_electrostaticcalculations[i].force_a, _electrostaticcalculations[i].idx_a);
805         AddGradient(_electrostaticcalculations[i].force_b, _electrostaticcalculations[i].idx_b);
806       }
807       #endif
808 
809       IF_OBFF_LOGLVL_HIGH {
810         snprintf(_logbuf, BUFF_SIZE, "%2d   %2d   %8.3f  %8.3f  %8.3f  %8.3f\n",
811                 atoi(_electrostaticcalculations[i].a->GetType()), atoi(_electrostaticcalculations[i].b->GetType()),
812                 _electrostaticcalculations[i].rab, _electrostaticcalculations[i].a->GetPartialCharge(),
813                 _electrostaticcalculations[i].b->GetPartialCharge(), _electrostaticcalculations[i].energy);
814         OBFFLog(_logbuf);
815       }
816     }
817 
818     #ifdef _OPENMP
819     for (int i = 0; i < _electrostaticcalculations.size(); ++i) {
820       // Cut-off check
821       if (_cutoff)
822         if (!_elepairs.BitIsSet(i))
823           continue;
824 
825       if (gradients) {
826         AddGradient(_electrostaticcalculations[i].force_a, _electrostaticcalculations[i].idx_a);
827         AddGradient(_electrostaticcalculations[i].force_b, _electrostaticcalculations[i].idx_b);
828       }
829     }
830     #endif
831 
832     IF_OBFF_LOGLVL_MEDIUM {
833       snprintf(_logbuf, BUFF_SIZE, "     TOTAL ELECTROSTATIC ENERGY = %8.5f %s\n", energy, GetUnit().c_str());
834       OBFFLog(_logbuf);
835     }
836 
837     return energy;
838   }
839 
840   //
841   // OBForceFieldMMFF member functions
842   //
843   //***********************************************
844   //Make a global instance
845   OBForceFieldMMFF94 theForceFieldMMFF94("MMFF94", false);
846   OBForceFieldMMFF94 theForceFieldMMFF94s("MMFF94s", false);
847   //***********************************************
848 
~OBForceFieldMMFF94()849   OBForceFieldMMFF94::~OBForceFieldMMFF94()
850   {
851   }
852 
operator =(OBForceFieldMMFF94 & src)853   OBForceFieldMMFF94 &OBForceFieldMMFF94::operator=(OBForceFieldMMFF94 &src)
854   {
855     _mol = src._mol;
856     _init = src._init;
857     return *this;
858   }
859 
860   ////////////////////////////////////////////////////////////////////////////////
861   ////////////////////////////////////////////////////////////////////////////////
862   //
863   //  Parse parameter files
864   //
865   ////////////////////////////////////////////////////////////////////////////////
866   ////////////////////////////////////////////////////////////////////////////////
867 
ParseParamFile()868   bool OBForceFieldMMFF94::ParseParamFile()
869   {
870     // Set the locale for number parsing to avoid locale issues: PR#1785463
871     obLocale.SetLocale();
872 
873     vector<string> vs;
874     char buffer[80];
875 
876     // open data/_parFile
877     ifstream ifs;
878     if (OpenDatafile(ifs, _parFile).length() == 0) {
879       obErrorLog.ThrowError(__FUNCTION__, "Cannot open parameter file", obError);
880       return false;
881     }
882 
883     while (ifs.getline(buffer, 80)) {
884       if (EQn(buffer, "#", 1)) continue;
885 
886       tokenize(vs, buffer);
887       if (vs.size() < 2)
888         continue;
889 
890       if (vs[0] == "prop")
891         ParseParamProp(vs[1]);
892       if (vs[0] == "def")
893         ParseParamDef(vs[1]);
894       if (vs[0] == "bond")
895         ParseParamBond(vs[1]);
896       if (vs[0] == "ang")
897         ParseParamAngle(vs[1]);
898       if (vs[0] == "bndk")
899         ParseParamBndk(vs[1]);
900       if (vs[0] == "chg")
901         ParseParamCharge(vs[1]);
902       if (vs[0] == "dfsb")
903         ParseParamDfsb(vs[1]);
904       if (vs[0] == "oop")
905         ParseParamOOP(vs[1]);
906       if (vs[0] == "pbci")
907         ParseParamPbci(vs[1]);
908       if (vs[0] == "stbn")
909         ParseParamStrBnd(vs[1]);
910       if (vs[0] == "tor")
911         ParseParamTorsion(vs[1]);
912       if (vs[0] == "vdw")
913         ParseParamVDW(vs[1]);
914     }
915 
916     if (ifs)
917       ifs.close();
918 
919     // return the locale to the original one
920     obLocale.RestoreLocale();
921     return true;
922   }
923 
ParseParamBond(std::string & filename)924   bool OBForceFieldMMFF94::ParseParamBond(std::string &filename)
925   {
926     vector<string> vs;
927     char buffer[80];
928 
929     OBFFParameter parameter;
930 
931     // open data/mmffbond.par
932     ifstream ifs;
933     if (OpenDatafile(ifs, filename).length() == 0) {
934       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffbond.par", obError);
935       return false;
936     }
937 
938     while (ifs.getline(buffer, 80)) {
939       if (EQn(buffer, "*", 1)) continue;
940       if (EQn(buffer, "$", 1)) continue;
941 
942       tokenize(vs, buffer);
943 
944       parameter.clear();
945       parameter._ipar.push_back(atoi(vs[0].c_str()));  // FF class
946       parameter.a = atoi(vs[1].c_str());
947       parameter.b = atoi(vs[2].c_str());
948       parameter._dpar.push_back(atof(vs[3].c_str()));  // kb
949       parameter._dpar.push_back(atof(vs[4].c_str()));  // r0
950       _ffbondparams.push_back(parameter);
951     }
952 
953     if (ifs)
954       ifs.close();
955 
956     return 0;
957   }
958 
ParseParamBndk(std::string & filename)959   bool OBForceFieldMMFF94::ParseParamBndk(std::string &filename)
960   {
961     vector<string> vs;
962     char buffer[80];
963 
964     OBFFParameter parameter;
965 
966     // open data/mmffbndk.par
967     ifstream ifs;
968     if (OpenDatafile(ifs, filename).length() == 0) {
969       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffbndk.par", obError);
970       return false;
971     }
972 
973     while (ifs.getline(buffer, 80)) {
974       if (EQn(buffer, "*", 1)) continue;
975       if (EQn(buffer, "$", 1)) continue;
976 
977       tokenize(vs, buffer);
978 
979       parameter.clear();
980       parameter.a = atoi(vs[0].c_str());
981       parameter.b = atoi(vs[1].c_str());
982       parameter._dpar.push_back(atof(vs[2].c_str()));  // r0-ref
983       parameter._dpar.push_back(atof(vs[3].c_str()));  // kb-ref
984       _ffbndkparams.push_back(parameter);
985     }
986 
987     if (ifs)
988       ifs.close();
989 
990     return 0;
991   }
992 
ParseParamAngle(std::string & filename)993   bool OBForceFieldMMFF94::ParseParamAngle(std::string &filename)
994   {
995     vector<string> vs;
996     char buffer[80];
997 
998     OBFFParameter parameter;
999 
1000     // open data/mmffang.par
1001     ifstream ifs;
1002     if (OpenDatafile(ifs, filename).length() == 0) {
1003       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffang.par", obError);
1004       return false;
1005     }
1006 
1007     while (ifs.getline(buffer, 80)) {
1008       if (EQn(buffer, "*", 1)) continue;
1009       if (EQn(buffer, "$", 1)) continue;
1010 
1011       tokenize(vs, buffer);
1012 
1013       parameter.clear();
1014       parameter._ipar.push_back(atoi(vs[0].c_str()));  // FF class
1015       parameter.a = atoi(vs[1].c_str());
1016       parameter.b = atoi(vs[2].c_str());
1017       parameter.c = atoi(vs[3].c_str());
1018       parameter._dpar.push_back(atof(vs[4].c_str()));  // ka
1019       parameter._dpar.push_back(atof(vs[5].c_str()));  // theta0
1020       _ffangleparams.push_back(parameter);
1021     }
1022 
1023     if (ifs)
1024       ifs.close();
1025 
1026     return 0;
1027   }
1028 
ParseParamStrBnd(std::string & filename)1029   bool OBForceFieldMMFF94::ParseParamStrBnd(std::string &filename)
1030   {
1031     vector<string> vs;
1032     char buffer[80];
1033 
1034     OBFFParameter parameter;
1035 
1036     // open data/mmffstbn.par
1037     ifstream ifs;
1038     if (OpenDatafile(ifs, filename).length() == 0) {
1039       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffstbn.par", obError);
1040       return false;
1041     }
1042 
1043     while (ifs.getline(buffer, 80)) {
1044       if (EQn(buffer, "*", 1)) continue;
1045       if (EQn(buffer, "$", 1)) continue;
1046 
1047       tokenize(vs, buffer);
1048 
1049       parameter.clear();
1050       parameter._ipar.push_back(atoi(vs[0].c_str()));  // FF class
1051       parameter.a = atoi(vs[1].c_str());
1052       parameter.b = atoi(vs[2].c_str());
1053       parameter.c = atoi(vs[3].c_str());
1054       parameter._dpar.push_back(atof(vs[4].c_str()));  // kbaIJK
1055       parameter._dpar.push_back(atof(vs[5].c_str()));  // kbaKJI
1056       _ffstrbndparams.push_back(parameter);
1057     }
1058 
1059     if (ifs)
1060       ifs.close();
1061 
1062     return 0;
1063   }
1064 
ParseParamDfsb(std::string & filename)1065   bool OBForceFieldMMFF94::ParseParamDfsb(std::string &filename)
1066   {
1067     vector<string> vs;
1068     char buffer[80];
1069 
1070     OBFFParameter parameter;
1071 
1072     // open data/mmffdfsb.par
1073     ifstream ifs;
1074     if (OpenDatafile(ifs, filename).length() == 0) {
1075       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffdfsb.par", obError);
1076       return false;
1077     }
1078 
1079     while (ifs.getline(buffer, 80)) {
1080       if (EQn(buffer, "*", 1)) continue;
1081       if (EQn(buffer, "$", 1)) continue;
1082 
1083       tokenize(vs, buffer);
1084 
1085       parameter.clear();
1086       parameter.a = atoi(vs[0].c_str());
1087       parameter.b = atoi(vs[1].c_str());
1088       parameter.c = atoi(vs[2].c_str());
1089       parameter._dpar.push_back(atof(vs[3].c_str()));  // kbaIJK
1090       parameter._dpar.push_back(atof(vs[4].c_str()));  // kbaKJI
1091       _ffdfsbparams.push_back(parameter);
1092     }
1093 
1094     if (ifs)
1095       ifs.close();
1096 
1097     return 0;
1098   }
1099 
ParseParamOOP(std::string & filename)1100   bool OBForceFieldMMFF94::ParseParamOOP(std::string &filename)
1101   {
1102     vector<string> vs;
1103     char buffer[80];
1104 
1105     OBFFParameter parameter;
1106 
1107     // open data/mmffoop.par
1108     ifstream ifs;
1109     if (OpenDatafile(ifs, filename).length() == 0) {
1110       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffoop.par", obError);
1111       return false;
1112     }
1113 
1114     while (ifs.getline(buffer, 80)) {
1115       if (EQn(buffer, "*", 1)) continue;
1116       if (EQn(buffer, "$", 1)) continue;
1117 
1118       tokenize(vs, buffer);
1119 
1120       parameter.clear();
1121       parameter.a = atoi(vs[0].c_str());
1122       parameter.b = atoi(vs[1].c_str());
1123       parameter.c = atoi(vs[2].c_str());
1124       parameter.d = atoi(vs[3].c_str());
1125       parameter._dpar.push_back(atof(vs[4].c_str()));  // koop
1126       _ffoopparams.push_back(parameter);
1127     }
1128 
1129     if (ifs)
1130       ifs.close();
1131 
1132     return 0;
1133   }
1134 
ParseParamTorsion(std::string & filename)1135   bool OBForceFieldMMFF94::ParseParamTorsion(std::string &filename)
1136   {
1137     vector<string> vs;
1138     char buffer[80];
1139 
1140     OBFFParameter parameter;
1141 
1142     // open data/mmfftor.par
1143     ifstream ifs;
1144     if (OpenDatafile(ifs, filename).length() == 0) {
1145       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmfftor.par", obError);
1146       return false;
1147     }
1148 
1149     while (ifs.getline(buffer, 80)) {
1150       if (EQn(buffer, "*", 1)) continue;
1151       if (EQn(buffer, "$", 1)) continue;
1152 
1153       tokenize(vs, buffer);
1154 
1155       parameter.clear();
1156       parameter._ipar.push_back(atoi(vs[0].c_str()));  // FF class
1157       parameter.a = atoi(vs[1].c_str());
1158       parameter.b = atoi(vs[2].c_str());
1159       parameter.c = atoi(vs[3].c_str());
1160       parameter.d = atoi(vs[4].c_str());
1161       parameter._dpar.push_back(atof(vs[5].c_str()));  // v1
1162       parameter._dpar.push_back(atof(vs[6].c_str()));  // v2
1163       parameter._dpar.push_back(atof(vs[7].c_str()));  // v3
1164       _fftorsionparams.push_back(parameter);
1165     }
1166 
1167     if (ifs)
1168       ifs.close();
1169 
1170     return 0;
1171   }
1172 
ParseParamVDW(std::string & filename)1173   bool OBForceFieldMMFF94::ParseParamVDW(std::string &filename)
1174   {
1175     vector<string> vs;
1176     char buffer[80];
1177 
1178     OBFFParameter parameter;
1179 
1180     // open data/mmffvdw.par
1181     ifstream ifs;
1182     if (OpenDatafile(ifs, filename).length() == 0) {
1183       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffvdw.par", obError);
1184       return false;
1185     }
1186 
1187     while (ifs.getline(buffer, 80)) {
1188       if (EQn(buffer, "*", 1)) continue;
1189       if (EQn(buffer, "$", 1)) continue;
1190 
1191       tokenize(vs, buffer);
1192 
1193       parameter.clear();
1194       parameter.a = atoi(vs[0].c_str());
1195       parameter._dpar.push_back(atof(vs[1].c_str()));  // alpha-i
1196       parameter._dpar.push_back(atof(vs[2].c_str()));  // N-i
1197       parameter._dpar.push_back(atof(vs[3].c_str()));  // A-i
1198       parameter._dpar.push_back(atof(vs[4].c_str()));  // G-i
1199       if (EQn(vs[5].c_str(), "-", 1))
1200         parameter._ipar.push_back(0);
1201       else if (EQn(vs[5].c_str(), "D", 1))
1202         parameter._ipar.push_back(1);  // hydrogen bond donor
1203       else if (EQn(vs[5].c_str(), "A", 1))
1204         parameter._ipar.push_back(2);  // hydrogen bond acceptor
1205       _ffvdwparams.push_back(parameter);
1206     }
1207 
1208     if (ifs)
1209       ifs.close();
1210 
1211     return 0;
1212   }
1213 
ParseParamCharge(std::string & filename)1214   bool OBForceFieldMMFF94::ParseParamCharge(std::string &filename)
1215   {
1216     vector<string> vs;
1217     char buffer[80];
1218 
1219     OBFFParameter parameter;
1220 
1221     // open data/mmffchg.par
1222     ifstream ifs;
1223     if (OpenDatafile(ifs, filename).length() == 0) {
1224       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffchg.par", obError);
1225       return false;
1226     }
1227 
1228     while (ifs.getline(buffer, 80)) {
1229       if (EQn(buffer, "*", 1)) continue;
1230       if (EQn(buffer, "$", 1)) continue;
1231 
1232       tokenize(vs, buffer);
1233 
1234       parameter.clear();
1235       parameter._ipar.push_back(atoi(vs[0].c_str()));  // FF class
1236       parameter.a = atoi(vs[1].c_str());
1237       parameter.b = atoi(vs[2].c_str());
1238       parameter._dpar.push_back(atof(vs[3].c_str()));  // bci
1239       _ffchgparams.push_back(parameter);
1240     }
1241 
1242     if (ifs)
1243       ifs.close();
1244 
1245     return 0;
1246   }
1247 
ParseParamPbci(std::string & filename)1248   bool OBForceFieldMMFF94::ParseParamPbci(std::string &filename)
1249   {
1250     vector<string> vs;
1251     char buffer[80];
1252 
1253     OBFFParameter parameter;
1254 
1255     // open data/mmffpbci.par
1256     ifstream ifs;
1257     if (OpenDatafile(ifs, filename).length() == 0) {
1258       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffpbci", obError);
1259       return false;
1260     }
1261 
1262     while (ifs.getline(buffer, 80)) {
1263       if (EQn(buffer, "*", 1)) continue;
1264       if (EQn(buffer, "$", 1)) continue;
1265 
1266       tokenize(vs, buffer);
1267 
1268       parameter.clear();
1269       parameter.a = atoi(vs[1].c_str());
1270       parameter._dpar.push_back(atof(vs[2].c_str()));  // pbci
1271       parameter._dpar.push_back(atof(vs[3].c_str()));  // fcadj
1272       _ffpbciparams.push_back(parameter);
1273     }
1274 
1275     if (ifs)
1276       ifs.close();
1277 
1278     return 0;
1279   }
1280 
ParseParamProp(std::string & filename)1281   bool OBForceFieldMMFF94::ParseParamProp(std::string &filename)
1282   {
1283     vector<string> vs;
1284     char buffer[80];
1285 
1286     OBFFParameter parameter;
1287 
1288     // open data/mmffprop.par
1289     ifstream ifs;
1290     if (OpenDatafile(ifs, filename).length() == 0) {
1291       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffprop.par", obError);
1292       return false;
1293     }
1294 
1295     while (ifs.getline(buffer, 80)) {
1296       if (EQn(buffer, "*", 1)) continue;
1297       if (EQn(buffer, "$", 1)) continue;
1298 
1299       tokenize(vs, buffer);
1300 
1301       parameter.clear();
1302       parameter.a = atoi(vs[0].c_str()); // atom type
1303       parameter._ipar.push_back(atoi(vs[1].c_str()));  // at.no
1304       parameter._ipar.push_back(atoi(vs[2].c_str()));  // crd
1305       parameter._ipar.push_back(atoi(vs[3].c_str()));  // val
1306       parameter._ipar.push_back(atoi(vs[4].c_str()));  // pilp
1307       parameter._ipar.push_back(atoi(vs[5].c_str()));  // mltb
1308       parameter._ipar.push_back(atoi(vs[6].c_str()));  // arom
1309       parameter._ipar.push_back(atoi(vs[7].c_str()));  // linh
1310       parameter._ipar.push_back(atoi(vs[8].c_str()));  // sbmb
1311 
1312       if (parameter._ipar[3])
1313         _ffpropPilp.SetBitOn(parameter.a);
1314       if (parameter._ipar[5])
1315         _ffpropArom.SetBitOn(parameter.a);
1316       if (parameter._ipar[6])
1317         _ffpropLin.SetBitOn(parameter.a);
1318       if (parameter._ipar[7])
1319         _ffpropSbmb.SetBitOn(parameter.a);
1320 
1321       _ffpropparams.push_back(parameter);
1322     }
1323 
1324     if (ifs)
1325       ifs.close();
1326 
1327     return 0;
1328   }
1329 
ParseParamDef(std::string & filename)1330   bool OBForceFieldMMFF94::ParseParamDef(std::string &filename)
1331   {
1332     vector<string> vs;
1333     char buffer[80];
1334 
1335     OBFFParameter parameter;
1336 
1337     // open data/mmffdef.par
1338     ifstream ifs;
1339     if (OpenDatafile(ifs, filename).length() == 0) {
1340       obErrorLog.ThrowError(__FUNCTION__, "Cannot open mmffdef.par", obError);
1341       return false;
1342     }
1343 
1344     while (ifs.getline(buffer, 80)) {
1345       if (EQn(buffer, "*", 1)) continue;
1346       if (EQn(buffer, "$", 1)) continue;
1347 
1348       tokenize(vs, buffer);
1349 
1350       parameter.clear();
1351       parameter._ipar.push_back(atoi(vs[1].c_str()));  // level 1
1352       parameter._ipar.push_back(atoi(vs[2].c_str()));  // level 2
1353       parameter._ipar.push_back(atoi(vs[3].c_str()));  // level 3
1354       parameter._ipar.push_back(atoi(vs[4].c_str()));  // level 4
1355       parameter._ipar.push_back(atoi(vs[5].c_str()));  // level 5
1356       _ffdefparams.push_back(parameter);
1357     }
1358 
1359     if (ifs)
1360       ifs.close();
1361 
1362     return 0;
1363   }
1364 
1365   ////////////////////////////////////////////////////////////////////////////////
1366   ////////////////////////////////////////////////////////////////////////////////
1367   //
1368   //  Setup Functions
1369   //
1370   ////////////////////////////////////////////////////////////////////////////////
1371   ////////////////////////////////////////////////////////////////////////////////
1372 
1373   // The MMFF94 article doesn't seem to include information about how
1374   // aromaticity is perceived. This function was written by studying the
1375   // MMFF_opti.log file, trail-and-error and using the MMFF94 validation
1376   // set to check the results (If all atom types are assigned correctly,
1377   // aromatic rings are probably detected correctly)
PerceiveAromatic()1378   bool OBForceFieldMMFF94::PerceiveAromatic()
1379   {
1380     bool done = false; // not done actually....
1381     OBAtom *ringatom;
1382     OBBond *ringbond;
1383     vector<OBRing*> vr;
1384     vr = _mol.GetSSSR();
1385 
1386     vector<OBRing*>::iterator ri;
1387     vector<int>::iterator rj;
1388     int n, index, ringsize, first_rj, prev_rj, pi_electrons, c60;
1389     for (ri = vr.begin();ri != vr.end();++ri) { // for each ring
1390       ringsize = (*ri)->Size();
1391 
1392       n = 1;
1393       pi_electrons = 0;
1394       c60 = 0; // we have a special case to get c60 right (all atom type 37)
1395       for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) { // for each ring atom
1396         index = *rj;
1397         ringatom = _mol.GetAtom(index);
1398 
1399         // is the bond to the previous ring atom double?
1400         if (n > 1) {
1401           ringbond = _mol.GetBond(prev_rj, index);
1402           if (!ringbond) {
1403             prev_rj = index;
1404             continue;
1405           }
1406           if (ringbond->GetBondOrder() == 2) {
1407             pi_electrons += 2;
1408             prev_rj = index;
1409             n++;
1410             continue;
1411           }
1412           prev_rj = index;
1413         } else {
1414           prev_rj = index;
1415           first_rj = index;
1416         }
1417 
1418         // does the current ring atom have a exocyclic double bond?
1419         FOR_NBORS_OF_ATOM (nbr, ringatom) {
1420           if ((*ri)->IsInRing(nbr->GetIdx()))
1421             continue;
1422 
1423           if (!nbr->IsAromatic()) {
1424             if (ringatom->GetAtomicNum() == OBElements::Carbon && ringatom->IsInRingSize(5)
1425                 && ringatom->IsInRingSize(6) && nbr->GetAtomicNum() == OBElements::Carbon && nbr->IsInRingSize(5)
1426                 && nbr->IsInRingSize(6)) {
1427               c60++;
1428             } else {
1429               continue;
1430             }
1431           }
1432 
1433           ringbond = _mol.GetBond(nbr->GetIdx(), index);
1434           if (!ringbond) {
1435             continue;
1436           }
1437           if (ringbond->GetBondOrder() == 2)
1438             pi_electrons++;
1439         }
1440 
1441         // is the atom N, O or S in 5 rings
1442         if (ringsize == 5 &&
1443             ringatom->GetIdx() == (*ri)->GetRootAtom())
1444           pi_electrons += 2;
1445 
1446         n++;
1447 
1448       } // for each ring atom
1449 
1450       // is the bond from the first to the last atom double?
1451       ringbond = _mol.GetBond(first_rj, index);
1452       if (ringbond) {
1453         if (ringbond->GetBondOrder() == 2)
1454           pi_electrons += 2;
1455       }
1456 
1457       if (((pi_electrons == 6) && ((ringsize == 5) || (ringsize == 6)))
1458         || ((pi_electrons == 5) && (c60 == 5))) {
1459         // mark ring atoms as aromatic
1460         for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) {
1461           if (!_mol.GetAtom(*rj)->IsAromatic())
1462             done = true;
1463           _mol.GetAtom(*rj)->SetAromatic();
1464         }
1465         // mark all ring bonds as aromatic
1466         FOR_BONDS_OF_MOL (bond, _mol)
1467           if((*ri)->IsMember(&*bond))
1468             bond->SetAromatic();
1469       }
1470     }
1471 
1472     return done;
1473   }
1474 
1475   // Symbolic atom typing is skipped
1476   //
1477   // atom typing is based on:
1478   //   MMFF94 I - Table III
1479   //   MMFF94 V - Table I
1480   //
GetType(OBAtom * atom)1481   int OBForceFieldMMFF94::GetType(OBAtom *atom)
1482   {
1483     OBBond *bond;
1484     int oxygenCount, nitrogenCount, sulphurCount, doubleBondTo;
1485     ////////////////////////////////
1486     // Aromatic Atoms
1487     ////////////////////////////////
1488     if (atom->IsAromatic()) {
1489       if (atom->IsInRingSize(5)) {
1490         bool IsAromatic = false;
1491         vector<OBAtom*> alphaPos, betaPos;
1492         vector<OBAtom*> alphaAtoms, betaAtoms;
1493 
1494         if (atom->GetAtomicNum() == OBElements::Sulfur) {
1495           return 44; // Aromatic 5-ring sulfur with pi lone pair (STHI)
1496         }
1497         if (atom->GetAtomicNum() == OBElements::Oxygen) {
1498           return 59; // Aromatic 5-ring oxygen with pi lone pair (OFUR)
1499         }
1500         if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1501           FOR_NBORS_OF_ATOM (nbr, atom) {
1502             if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1503               return 82; // N-oxide nitrogen in 5-ring alpha position,
1504               // N-oxide nitrogen in 5-ring beta position,
1505               // N-oxide nitrogen in other 5-ring  position,
1506               // (N5AX, N5BX, N5OX)
1507             }
1508           }
1509         }
1510         FOR_NBORS_OF_ATOM (nbr, atom) {
1511           if (!((_mol.GetBond(atom, &*nbr))->IsAromatic()) || !nbr->IsInRingSize(5))
1512             continue;
1513 
1514           if (IsInSameRing(atom, &*nbr)) {
1515             alphaPos.push_back(&*nbr);
1516           }
1517 
1518           FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
1519             if (nbrNbr->GetIdx() == atom->GetIdx())
1520               continue;
1521             if (!((_mol.GetBond(&*nbr, &*nbrNbr))->IsAromatic()) || !nbrNbr->IsInRingSize(5))
1522               continue;
1523 
1524             IsAromatic = true;
1525 
1526             if (IsInSameRing(atom, &*nbrNbr)) {
1527               betaPos.push_back(&*nbrNbr);
1528             }
1529           }
1530         }
1531 
1532         if (IsAromatic) {
1533 
1534 
1535           for (unsigned int i = 0; i < alphaPos.size(); i++) {
1536             if (alphaPos[i]->GetAtomicNum() == OBElements::Sulfur) {
1537               alphaAtoms.push_back(alphaPos[i]);
1538             } else if (alphaPos[i]->GetAtomicNum() == OBElements::Oxygen) {
1539               alphaAtoms.push_back(alphaPos[i]);
1540             } else if (alphaPos[i]->GetAtomicNum() == OBElements::Nitrogen && (alphaPos[i]->GetExplicitDegree() == 3)) {
1541               bool IsNOxide = false;
1542               FOR_NBORS_OF_ATOM (nbr, alphaPos[i]) {
1543                 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1544                   IsNOxide = true;
1545                 }
1546               }
1547 
1548               if (!IsNOxide) {
1549                 alphaAtoms.push_back(alphaPos[i]);
1550               }
1551             }
1552           }
1553           for (unsigned int i = 0; i < betaPos.size(); i++) {
1554             if (betaPos[i]->GetAtomicNum() == OBElements::Sulfur) {
1555               betaAtoms.push_back(betaPos[i]);
1556             } else if (betaPos[i]->GetAtomicNum() == OBElements::Oxygen) {
1557               betaAtoms.push_back(betaPos[i]);
1558             } else if (betaPos[i]->GetAtomicNum() == OBElements::Nitrogen && (betaPos[i]->GetExplicitDegree() == 3)) {
1559               bool IsNOxide = false;
1560               FOR_NBORS_OF_ATOM (nbr, betaPos[i]) {
1561                 if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1562                   IsNOxide = true;
1563                 }
1564               }
1565 
1566               if (!IsNOxide) {
1567                 betaAtoms.push_back(betaPos[i]);
1568               }
1569             }
1570           }
1571           if (!betaAtoms.size()) {
1572             nitrogenCount = 0;
1573             FOR_NBORS_OF_ATOM (nbr, atom) {
1574               if (nbr->GetAtomicNum() == OBElements::Nitrogen && (nbr->GetExplicitDegree() == 3)) {
1575                 if ((nbr->GetExplicitValence() == 4) && nbr->IsAromatic()) {
1576                   nitrogenCount++;
1577                 } else if ((nbr->GetExplicitValence() == 3) && !nbr->IsAromatic()) {
1578                   nitrogenCount++;
1579                 }
1580               }
1581             }
1582             if (nitrogenCount >= 2) {
1583               return 80; // Aromatic carbon between N's in imidazolium (CIM+)
1584             }
1585           }
1586           if (!alphaAtoms.size() && !betaAtoms.size()) {
1587             if (atom->GetAtomicNum() == OBElements::Carbon) {
1588               bool c60 = true; // special case to ensure c60 is typed correctly -- Paolo Tosco
1589               FOR_NBORS_OF_ATOM (nbr, atom) {
1590                 if (!(nbr->GetAtomicNum() == OBElements::Carbon && nbr->IsAromatic() && nbr->IsInRingSize(6)))
1591                   c60 = false;
1592               }
1593               if (c60)
1594                 return 37; // correct atom type for c in c60 (all atoms symmetric)
1595               // there is no S:, O:, or N:
1596               // this is the case for anions with only carbon and nitrogen in the ring
1597               return 78; // General carbon in 5-membered aromatic ring (C5)
1598             } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1599               if (atom->GetExplicitDegree() == 3) {
1600                 // this is the N: atom
1601                 return 39; // Aromatic 5 ring nitrogen with pi lone pair (NPYL)
1602               } else {
1603                 // again, no S:, O:, or N:
1604                 return 76; // Nitrogen in 5-ring aromatic anion (N5M)
1605               }
1606             }
1607           }
1608           if (alphaAtoms.size() == 2) {
1609             if (atom->GetAtomicNum() == OBElements::Carbon && IsInSameRing(alphaAtoms[0], alphaAtoms[1])) {
1610               if (alphaAtoms[0]->GetAtomicNum() == OBElements::Nitrogen && alphaAtoms[1]->GetAtomicNum() == OBElements::Nitrogen) {
1611                 if ((alphaAtoms[0]->GetExplicitDegree() == 3) && (alphaAtoms[1]->GetExplicitDegree() == 3)) {
1612                   return 80; // Aromatic carbon between N's in imidazolium (CIM+)
1613                 }
1614               }
1615             }
1616           }
1617           if (alphaAtoms.size() && !betaAtoms.size()) {
1618             if (atom->GetAtomicNum() == OBElements::Carbon) {
1619               return 63; // Aromatic 5-ring C, alpha to N:, O:, or S: (C5A)
1620             } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1621               if (atom->GetExplicitDegree() == 3) {
1622                 return 81; // Posivite nitrogen in 5-ring alpha position (N5A+)
1623               } else {
1624                 return 65; // Aromatic 5-ring N, alpha to N:, O:, or S: (N5A)
1625               }
1626             }
1627           }
1628           if (!alphaAtoms.size() && betaAtoms.size()) {
1629             if (atom->GetAtomicNum() == OBElements::Carbon) {
1630               return 64; // Aromatic 5-ring C, beta to N:, O:, or S: (C5B)
1631             } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1632               if (atom->GetExplicitDegree() == 3) {
1633                 return 81; // Posivite nitrogen in 5-ring beta position (N5B+)
1634               } else {
1635                 return 66; // Aromatic 5-ring N, beta to N:, O:, or S: (N5B)
1636               }
1637             }
1638           }
1639           if (alphaAtoms.size() && betaAtoms.size()) {
1640             for (unsigned int i = 0; i < alphaAtoms.size(); i++) {
1641               for (unsigned int j = 0; j < betaAtoms.size(); j++) {
1642                 if (!IsInSameRing(alphaAtoms[i], betaAtoms[j])) {
1643                   if (atom->GetAtomicNum() == OBElements::Carbon) {
1644                     return 78; // General carbon in 5-membered aromatic ring (C5)
1645                   } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1646                     return 79; // General nitrogen in 5-membered aromatic ring (N5)
1647                   }
1648                 }
1649               }
1650             }
1651             for (unsigned int i = 0; i < alphaAtoms.size(); i++) {
1652               if (alphaAtoms[i]->GetAtomicNum() == OBElements::Sulfur || alphaAtoms[i]->GetAtomicNum() == OBElements::Oxygen) {
1653                 if (atom->GetAtomicNum() == OBElements::Carbon) {
1654                   return 63; // Aromatic 5-ring C, alpha to N:, O:, or S: (C5A)
1655                 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1656                   return 65; // Aromatic 5-ring N, alpha to N:, O:, or S: (N5A)
1657                 }
1658               }
1659             }
1660             for (unsigned int i = 0; i < betaAtoms.size(); i++) {
1661               if (betaAtoms[i]->GetAtomicNum() == OBElements::Sulfur || betaAtoms[i]->GetAtomicNum() == OBElements::Oxygen) {
1662                 if (atom->GetAtomicNum() == OBElements::Carbon) {
1663                   return 64; // Aromatic 5-ring C, beta to N:, O:, or S: (C5B)
1664                 } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1665                   return 66; // Aromatic 5-ring N, beta to N:, O:, or S: (N5B)
1666                 }
1667               }
1668             }
1669 
1670             if (atom->GetAtomicNum() == OBElements::Carbon) {
1671               return 78; // General carbon in 5-membered aromatic ring (C5)
1672             } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1673               return 79; // General nitrogen in 5-membered aromatic ring (N5)
1674             }
1675           }
1676         }
1677       }
1678 
1679       if (atom->IsInRingSize(6)) {
1680 
1681         if (atom->GetAtomicNum() == OBElements::Carbon) {
1682           return 37; // Aromatic carbon, e.g., in benzene (CB)
1683         } else if (atom->GetAtomicNum() == OBElements::Nitrogen) {
1684           FOR_NBORS_OF_ATOM (nbr, atom) {
1685             if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1686               return 69; // Pyridinium N-oxide nitrogen (NPOX)
1687             }
1688           }
1689 
1690           if (atom->GetExplicitDegree() == 3) {
1691             return 58; // Aromatic nitrogen in pyridinium (NPD+)
1692           } else {
1693             return 38; // Aromatic nitrogen with sigma lone pair (NPYD)
1694           }
1695         }
1696       }
1697     }
1698 
1699     ////////////////////////////////
1700     // Hydrogen
1701     ////////////////////////////////
1702     if (atom->GetAtomicNum() == 1) {
1703       FOR_NBORS_OF_ATOM (nbr, atom) {
1704         if (nbr->GetAtomicNum() == OBElements::Carbon) {
1705           return 5; // Hydrogen attatched to carbon (HC)
1706         }
1707         if (nbr->GetAtomicNum() == 14) {
1708           return 5; // Hydrogen attatched to silicon (HSI)
1709         }
1710         if (nbr->GetAtomicNum() == OBElements::Oxygen) {
1711           if (nbr->GetExplicitValence() == 3) {
1712             if (nbr->GetExplicitDegree() == 3) {
1713               return 50; // Hydrogen on oxonium oxygen (HO+)
1714             } else {
1715               return 52; // Hydrogen on oxenium oxygen (HO=+)
1716             }
1717           }
1718 
1719           int hydrogenCount = 0;
1720           FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
1721             if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen) {
1722               hydrogenCount++;
1723               continue;
1724             }
1725             if (nbrNbr->GetAtomicNum() == OBElements::Carbon) {
1726               if (nbrNbr->IsAromatic()) {
1727                 return 29; // phenol
1728               }
1729 
1730               FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
1731                 if (nbrNbrNbr->GetIdx() == nbr->GetIdx())
1732                   continue;
1733 
1734                 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr);
1735                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1736                   if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) {
1737                     return 24; // Hydroxyl hydrogen in carboxylic acids (HOCO)
1738                   }
1739                   if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen) {
1740                     return 29; // Enolic or phenolic hydroxyl hydrogen,
1741                     // Hydroxyl hydrogen in HO-C=N moiety (HOCC, HOCN)
1742                   }
1743                 }
1744               }
1745             }
1746             if (nbrNbr->GetAtomicNum() == OBElements::Phosphorus) {
1747               return 24; // Hydroxyl hydrogen in H-O-P moiety (HOP)
1748             }
1749             if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) {
1750               return 33; // Hydrogen on oxygen attached to sulfur (HOS)
1751             }
1752 
1753           }
1754           if (hydrogenCount == 2) {
1755             return 31; // Hydroxyl hydrogen in water (HOH)
1756           }
1757 
1758           return 21; // Hydroxyl hydrogen in alcohols, Generic hydroxyl hydrogen (HOR, HO)
1759         }
1760         if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
1761           switch (GetType(&*nbr)) {
1762           case 81:
1763             return 36; // Hydrogen on imidazolium nitrogen (HIM+)
1764           case 68:
1765             return 23; // Hydrogen on N in N-oxide (HNOX)
1766           case 67:
1767             return 23; // Hydrogen on N in N-oxide (HNOX)
1768           case 62:
1769             return 23; // Generic hydrogen on sp3 nitrogen, e.g., in amines (HNR)
1770           case 56:
1771             return 36; // Hydrogen on guanimdinium nitrogen (HGD+)
1772           case 55:
1773             return 36; // Hydrogen on amidinium nitrogen (HNN+)
1774           case 43:
1775             return 28; // Hydrogen on NSO, NSO2, or NSO3 nitrogen, Hydrogen on N triply bonded to C (HNSO, HNC%)
1776           case 39:
1777             return 23; // Hydrogen on nitrogen in pyrrole (HPYL)
1778           case 8:
1779             return 23; // Generic hydrogen on sp3 nitrogen, e.g., in amines, Hydrogen on nitrogen in ammonia (HNR, H3N)
1780           }
1781 
1782           if (nbr->GetExplicitValence() == 4) {
1783             if (nbr->GetExplicitDegree() == 2) {
1784               return 28; // Hydrogen on N triply bonded to C (HNC%)
1785             } else {
1786               return 36; // Hydrogen on pyridinium nitrogen, Hydrogen on protonated imine nitrogen (HPD+, HNC+)
1787             }
1788           }
1789 
1790           if (nbr->GetExplicitDegree() == 2) {
1791             FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
1792               if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen)
1793                 continue;
1794 
1795               bond = _mol.GetBond(&*nbr, &*nbrNbr);
1796               if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1797                 if (nbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbr->GetAtomicNum() == OBElements::Nitrogen) {
1798                   return 27; // Hydrogen on imine nitrogen, Hydrogen on azo nitrogen (HN=C, HN=N)
1799                 }
1800 
1801                 return 28; // Generic hydrogen on sp2 nitrogen (HSP2)
1802               }
1803             }
1804           }
1805 
1806           FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
1807             if (nbrNbr->GetAtomicNum() == OBElements::Hydrogen)
1808               continue;
1809 
1810             if (nbrNbr->GetAtomicNum() == OBElements::Carbon) {
1811               if (nbrNbr->IsAromatic()) {
1812                 return 28; // deloc. lp pair
1813               }
1814 
1815               FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
1816                 if (nbrNbrNbr->GetIdx() == nbr->GetIdx())
1817                   continue;
1818 
1819                 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr);
1820                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1821                   if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen || nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen || nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) {
1822                     return 28; // Hydrogen on amide nitrogen, Hydrogen on thioamide nitrogen,
1823                     // Hydrogen on enamine nitrogen, Hydrogen in H-N-C=N moiety (HNCO, HNCS, HNCC, HNCN)
1824                   }
1825                 }
1826               }
1827             }
1828             if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen) {
1829               FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
1830                 if (nbrNbrNbr->GetIdx() == nbr->GetIdx())
1831                   continue;
1832 
1833                 bond = _mol.GetBond(&*nbrNbr, &*nbrNbrNbr);
1834                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1835                   if (nbrNbrNbr->GetAtomicNum() == OBElements::Carbon || nbrNbrNbr->GetAtomicNum() == OBElements::Nitrogen) {
1836                     return 28; // Hydrogen in H-N-N=C moiety, Hydrogen in H-N-N=N moiety (HNNC, HNNN)
1837                   }
1838                 }
1839               }
1840             }
1841             if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) {
1842               FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
1843                 if (nbrNbrNbr->GetIdx() == nbr->GetIdx())
1844                   continue;
1845 
1846                 if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen || (nbrNbrNbr->GetExplicitDegree() == 1)) {
1847                   return 28; // Hydrogen on NSO, NSO2 or NSO3 nitrogen (HNSO)
1848                 }
1849               }
1850             }
1851           }
1852 
1853           return 23; // Generic hydrogen on sp3 nitrogen e.g., in amines,
1854           // Hydrogen on nitrogen in pyrrole, Hydrogen in ammonia,
1855           // Hydrogen on N in N-oxide (HNR, HPYL, H3N, HNOX)
1856         }
1857         if (nbr->GetAtomicNum() == OBElements::Sulfur || nbr->GetAtomicNum() == OBElements::Phosphorus) {
1858           return 71; // Hydrogen attached to sulfur, Hydrogen attached to >S= sulfur doubly bonded to N,
1859           // Hydrogen attached to phosphorus (HS, HS=N, HP)
1860         }
1861       }
1862     }
1863 
1864     ////////////////////////////////
1865     // Lithium
1866     ////////////////////////////////
1867     if (atom->GetAtomicNum() == 3) {
1868       // 0 neighbours
1869       if (atom->GetExplicitDegree() == 0) {
1870         return 92; // Lithium cation (LI+)
1871       }
1872     }
1873 
1874     ////////////////////////////////
1875     // Carbon
1876     ////////////////////////////////
1877     if (atom->GetAtomicNum() == 6) {
1878       // 4 neighbours
1879       if (atom->GetExplicitDegree() == 4) {
1880         if (atom->IsInRingSize(3)) {
1881           return 22; // Aliphatic carbon in 3-membered ring (CR3R)
1882         }
1883 
1884         if (atom->IsInRingSize(4)) {
1885           return 20; // Aliphatic carbon in 4-membered ring (CR4R)
1886         }
1887 
1888         return 1; // Alkyl carbon (CR)
1889       }
1890       // 3 neighbours
1891       if (atom->GetExplicitDegree() == 3) {
1892         int N2count = 0;
1893         int N3count = 0;
1894         int N3fcharge = 0;
1895         oxygenCount = sulphurCount = doubleBondTo = 0;
1896 
1897         FOR_NBORS_OF_ATOM (nbr, atom) {
1898           bond = _mol.GetBond(&*nbr, atom);
1899           if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1900             doubleBondTo = nbr->GetAtomicNum();
1901           }
1902 
1903           if (nbr->GetExplicitDegree() == 1) {
1904             if (nbr->GetAtomicNum() == OBElements::Oxygen) {
1905               oxygenCount++;
1906             } else if (nbr->GetAtomicNum() == OBElements::Sulfur) {
1907               sulphurCount++;
1908             }
1909           } else if (nbr->GetExplicitDegree() == 3) {
1910             if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
1911               N3fcharge += nbr->GetFormalCharge();
1912               N3count++;
1913             }
1914           } else if ((nbr->GetExplicitDegree() == 2) && !bond->IsAromatic() && bond->GetBondOrder() == 2) {
1915             if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
1916               N2count++;
1917             }
1918           }
1919         }
1920         if ((N3count >= 2) && (doubleBondTo == 7 || (!(doubleBondTo == 6) && atom->GetExplicitDegree() == 3
1921           && N3fcharge == 1)) && !N2count && !oxygenCount && !sulphurCount) {
1922           // N3==C--N3
1923           return 57; // Guanidinium carbon, Carbon in +N=C-N: resonance structures (CGD+, CNN+)
1924         }
1925         if ((oxygenCount == 2) || (sulphurCount == 2)) {
1926           // O1-?-C-?-O1 or S1-?-C-?-S1
1927           return 41; // Carbon in carboxylate anion, Carbon in thiocarboxylate anion (CO2M, CS2M)
1928         }
1929         if (atom->IsInRingSize(4) && (doubleBondTo == 6)) {
1930 	        return 30; // Olefinic carbon in 4-membered ring (CR4E)
1931         }
1932         if ((doubleBondTo ==  7) || (doubleBondTo ==  8) ||
1933             (doubleBondTo == 15) || (doubleBondTo == 16)) {
1934           // C==N, C==O, C==P, C==S
1935           return 3; // Generic carbonyl carbon, Imine-type carbon, Guanidine carbon,
1936           // Ketone or aldehyde carbonyl carbon, Amide carbonyl carbon,
1937           // Carboxylic acid or ester carbonyl carbon, Carbamate carbonyl carbon,
1938           // Carbonic acid or ester carbonyl carbon, Thioester carbonyl (double
1939           // bonded to O or S), Thioamide carbon (double bonded to S), Carbon
1940           // in >C=SO2, Sulfinyl carbon in >C=S=O, Thiocarboxylic acid or ester
1941           // carbon, Carbon doubly bonded to P (C=O, C=N, CGD, C=OR, C=ON, COO,
1942           // COON, COOO, C=OS, C=S, C=SN, CSO2, CS=O, CSS, C=P)
1943         }
1944 
1945         return 2; // Vinylic Carbon, Generic sp2 carbon (C=C, CSP2)
1946 
1947       }
1948       // 2 neighbours
1949       if (atom->GetExplicitDegree() == 2) {
1950         return 4; // Acetylenic carbon, Allenic caron (CSP, =C=)
1951       }
1952       // 1 neighbours
1953       if (atom->GetExplicitDegree() == 1) {
1954         return 60; // Isonitrile carbon (C%-)
1955       }
1956     }
1957 
1958     ////////////////////////////////
1959     // Nitrogen
1960     ////////////////////////////////
1961     if (atom->GetAtomicNum() == 7) {
1962       // 4 neighbours
1963       if (atom->GetExplicitDegree() == 4) {
1964         FOR_NBORS_OF_ATOM (nbr, atom) {
1965           if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1966             return 68; // sp3-hybridized N-oxide nitrogen (N3OX)
1967           }
1968         }
1969 
1970         return 34; // Quaternary nitrogen (NR+)
1971       }
1972       // 3 neighbours
1973       if (atom->GetExplicitDegree() == 3) {
1974         if (atom->GetExplicitValence() >= 4) { // > because we accept -N(=O)=O as a valid nitro group
1975           oxygenCount = nitrogenCount = doubleBondTo = 0;
1976 
1977           FOR_NBORS_OF_ATOM (nbr, atom) {
1978             if (nbr->GetAtomicNum() == OBElements::Oxygen && (nbr->GetExplicitDegree() == 1)) {
1979               oxygenCount++;
1980             }
1981             if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
1982               bond = _mol.GetBond(&*nbr, atom);
1983               if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1984                 doubleBondTo = 7;
1985               }
1986             }
1987             if (nbr->GetAtomicNum() == OBElements::Carbon) {
1988               bond = _mol.GetBond(&*nbr, atom);
1989               if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
1990                 FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
1991                   if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 3)) {
1992                     nitrogenCount++;
1993                   }
1994                 }
1995               }
1996             }
1997           }
1998 
1999           if (oxygenCount == 1) {
2000             return 67; // sp2-hybridized N-oxide nitrogen (N2OX)
2001           }
2002           if (oxygenCount >= 2) {
2003             return 45; // Nitrogen in nitro group, Nitrogen in nitrate group (NO2, NO3)
2004           }
2005 
2006           if (nitrogenCount == 1) {
2007             return 54; // Iminium nitrogen (N+=C)
2008           }
2009           if (nitrogenCount == 2) {
2010             return 55; // Either nitrogen in N+=C-N: (NCN+)
2011           }
2012           if (nitrogenCount == 3) {
2013             return 56; // Guanidinium nitrogen (NGD+)
2014           }
2015 
2016           if (doubleBondTo == 7) {
2017             return 54; // Positivly charged nitrogen doubly bonded to nitrogen (N+=N)
2018           }
2019         }
2020 
2021         if (atom->GetExplicitValence() == 3) {
2022           bool IsAmide = false;
2023           bool IsSulfonAmide = false;
2024           bool IsNNNorNNC = false;
2025           int tripleBondTo = 0;
2026           doubleBondTo = 0;
2027 
2028           FOR_NBORS_OF_ATOM (nbr, atom) {
2029             if (nbr->GetAtomicNum() == OBElements::Sulfur || nbr->GetAtomicNum() == OBElements::Phosphorus) {
2030               oxygenCount = 0;
2031 
2032               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2033                 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) {
2034                   oxygenCount++;
2035                 }
2036               }
2037               if (oxygenCount >= 2) {
2038                 IsSulfonAmide = true;
2039                 //return 43; // Sulfonamide nitrogen (NSO2, NSO3)
2040               }
2041             }
2042           }
2043 
2044           FOR_NBORS_OF_ATOM (nbr, atom) {
2045             if (nbr->GetAtomicNum() == OBElements::Carbon) {
2046               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2047                 bond = _mol.GetBond(&*nbr, &*nbrNbr);
2048                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2 && (nbrNbr->GetAtomicNum() == OBElements::Oxygen || nbrNbr->GetAtomicNum() == OBElements::Sulfur)) {
2049                   IsAmide = true;
2050                   //return 10; // Amide nitrogen, Thioamide nitrogen (NC=O, NC=S)
2051                 }
2052               }
2053             }
2054           }
2055 
2056           FOR_NBORS_OF_ATOM (nbr, atom) {
2057             if (nbr->GetAtomicNum() == OBElements::Carbon) {
2058               int N2count = 0;
2059               int N3count = 0;
2060               oxygenCount = sulphurCount = 0;
2061 
2062               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2063                 bond = _mol.GetBond(&*nbr, &*nbrNbr);
2064                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
2065                   doubleBondTo = nbrNbr->GetAtomicNum();
2066                 }
2067                 if (bond->IsAromatic()) {
2068                   if ((nbrNbr->GetAtomicNum() == 7) || (nbrNbr->GetAtomicNum() == 6)) {
2069                     doubleBondTo = nbrNbr->GetAtomicNum();
2070                   }
2071                 }
2072                 if (!bond->IsAromatic() && bond->GetBondOrder() == 3) {
2073                   tripleBondTo = nbrNbr->GetAtomicNum();
2074                 }
2075                 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 3)) {
2076                   int nbrOxygen = 0;
2077                   FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
2078                     if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) {
2079                       nbrOxygen++;
2080                     }
2081                   }
2082                   if (nbrOxygen < 2) {
2083                     N3count++;
2084                   }
2085                 }
2086                 if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen && (nbrNbr->GetExplicitDegree() == 2) && (bond->GetBondOrder() == 2 || bond->IsAromatic())) {
2087                   N2count++;
2088                 }
2089                 if (nbrNbr->IsAromatic()) {
2090                   if (nbrNbr->GetAtomicNum() == OBElements::Oxygen) {
2091                     oxygenCount++;
2092                   }
2093                   if (nbrNbr->GetAtomicNum() == OBElements::Sulfur) {
2094                     sulphurCount++;
2095                   }
2096                 }
2097               }
2098               if (N3count == 3) {
2099                 return 56; // Guanidinium nitrogen (NGD+)
2100               }
2101 
2102               if (!IsAmide && !IsSulfonAmide && !oxygenCount && !sulphurCount && nbr->IsAromatic()) {
2103                 return 40;
2104               }
2105 
2106               if ((N3count == 2) && (doubleBondTo == 7) && !N2count) {
2107                 return 55; // Either nitrogen in N+=C-N: (NCN+)
2108               }
2109             }
2110 
2111             if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
2112               nitrogenCount = 0;
2113               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2114                 bond = _mol.GetBond(&*nbr, &*nbrNbr);
2115                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
2116                   if (nbrNbr->GetAtomicNum() == OBElements::Carbon) {
2117                     oxygenCount = sulphurCount = 0;
2118                     FOR_NBORS_OF_ATOM (nbrNbrNbr, &*nbrNbr) {
2119                       if (nbrNbrNbr->GetAtomicNum() == OBElements::Oxygen) {
2120                         oxygenCount++;
2121                       }
2122                       if (nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) {
2123                         sulphurCount++;
2124                       }
2125                       if (nbrNbrNbr->GetAtomicNum() == OBElements::Sulfur) {
2126                         nitrogenCount++;
2127                       }
2128                     }
2129                     if (!oxygenCount && !sulphurCount && (nitrogenCount == 1)) {
2130                       bool bondToAromC = false;
2131                       FOR_NBORS_OF_ATOM (nbr2, atom) {
2132                         if (nbr2->IsAromatic() && nbr2->GetAtomicNum() == OBElements::Carbon && nbr2->IsInRingSize(6)) {
2133                           bondToAromC = true;
2134                         }
2135                       }
2136                       if (!bondToAromC) {
2137                         IsNNNorNNC = true;
2138                       }
2139                     }
2140                   }
2141                   if (nbrNbr->GetAtomicNum() == OBElements::Nitrogen) {
2142                     bool bondToAromC = false;
2143                     FOR_NBORS_OF_ATOM (nbr2, atom) {
2144                       if (nbr2->IsAromatic() && nbr2->GetAtomicNum() == OBElements::Carbon && nbr2->IsInRingSize(6)) {
2145                         bondToAromC = true;
2146                       }
2147                     }
2148                     if (!bondToAromC) {
2149                       IsNNNorNNC = true;
2150                     }
2151                   }
2152                 }
2153               }
2154             }
2155           }
2156 
2157           if (IsSulfonAmide) {
2158             return 43; // Sulfonamide nitrogen (NSO2, NSO3)
2159           }
2160           if (IsAmide) {
2161             return 10; // Amide nitrogen, Thioamide nitrogen (NC=O, NC=S)
2162           }
2163 
2164           if ((doubleBondTo ==  6) || (doubleBondTo == 7) ||(doubleBondTo == 15) || (tripleBondTo == 6)) {
2165             return 40; // Enamine or aniline nitrogen (deloc. lp), Nitrogen in N-C=N with deloc. lp,
2166             // Nitrogen in N-C=N with deloc. lp, Nitrogen attached to C-C triple bond
2167             // (NC=C, NC=N, NC=P, NC%C)
2168           }
2169           if (tripleBondTo == 7) {
2170             return 43; // Nitrogen attached to cyano group (NC%N)
2171           }
2172           if (IsNNNorNNC) {
2173             return 10; // Nitrogen in N-N=C moiety with deloc. lp
2174             // Nitrogen in N-N=N moiety with deloc. lp (NN=C, NN=N)
2175           }
2176 
2177           return 8; // Amine nitrogen (NR)
2178         }
2179       }
2180       // 2 neighbours
2181       if (atom->GetExplicitDegree() == 2) {
2182         if (atom->GetExplicitValence() >= 4) {
2183           bool isNbrCarbon = false;
2184           bool isBondTriple = false;
2185           FOR_NBORS_OF_ATOM (nbr, atom) {
2186             if (!isNbrCarbon)
2187               isNbrCarbon = nbr->GetAtomicNum() == OBElements::Carbon;
2188             bond = _mol.GetBond(&*nbr, atom);
2189             if (!isBondTriple)
2190               isBondTriple = !bond->IsAromatic() && bond->GetBondOrder() == 3;
2191           }
2192           if (isBondTriple && isNbrCarbon)
2193             return 61; // Isonitrile nitrogen (NR%)
2194 
2195           return 53; // Central nitrogen in C=N=N or N=N=N (=N=)
2196         }
2197 
2198         if (atom->GetExplicitValence() == 3) {
2199           doubleBondTo = 0;
2200 
2201           FOR_NBORS_OF_ATOM (nbr, atom) {
2202             bond = _mol.GetBond(&*nbr, atom);
2203             if (nbr->GetAtomicNum() == OBElements::Oxygen && !bond->IsAromatic() && bond->GetBondOrder() == 2 && (nbr->GetExplicitDegree() == 1)) {
2204               return 46; // Nitrogen in nitroso group (N=O)
2205             }
2206             if ((nbr->GetAtomicNum() == OBElements::Carbon || nbr->GetAtomicNum() == OBElements::Nitrogen) && !bond->IsAromatic() && bond->GetBondOrder() == 2) {
2207               return 9; // Iminie nitrogen, Azo-group nitrogen (N=C, N=N)
2208             }
2209           }
2210           FOR_NBORS_OF_ATOM (nbr, atom) {
2211             if (nbr->GetAtomicNum() == OBElements::Sulfur) {
2212               oxygenCount = 0;
2213 
2214               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2215                 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) {
2216                   oxygenCount++;
2217                 }
2218               }
2219               if (oxygenCount >= 2) {
2220                 return 43; // Sulfonamide nitrogen (NSO2, NSO3)
2221               }
2222             }
2223           }
2224         }
2225 
2226         if (atom->GetExplicitValence() >= 2) { // Bug reported by Paolo Tosco
2227           oxygenCount = sulphurCount = 0;
2228 
2229           FOR_NBORS_OF_ATOM (nbr, atom) {
2230             if (nbr->GetAtomicNum() == OBElements::Sulfur) {
2231               FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2232                 if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) {
2233                   oxygenCount++;
2234                 }
2235               }
2236               if (oxygenCount == 1) {
2237                 return 48; // Divalent nitrogen replacing monovalent O in SO2 group (NSO)
2238               }
2239             }
2240           }
2241 
2242           return 62; // Anionic divalent nitrogen (NM)
2243         }
2244       }
2245       // 1 neighbours
2246       if (atom->GetExplicitDegree() == 1) {
2247        	FOR_NBORS_OF_ATOM (nbr, atom) {
2248           bond = _mol.GetBond(&*nbr, atom);
2249           if (!bond->IsAromatic() && bond->GetBondOrder() == 3) {
2250             FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2251               if (atom != &*nbrNbr && !(nbr->GetAtomicNum() == OBElements::Nitrogen
2252                 && nbrNbr->GetAtomicNum() == OBElements::Nitrogen && nbrNbr->GetExplicitDegree() == 2)) {
2253                 return 42; // Triply bonded nitrogen (NSP)
2254               }
2255             }
2256           }
2257           if (nbr->GetAtomicNum() == OBElements::Nitrogen && (nbr->GetExplicitDegree() == 2)) {
2258             return 47; // Terminal nitrogen in azido group (NAZT)
2259           }
2260         }
2261       }
2262       return 8; // generic amine nitrogen
2263     }
2264 
2265     ////////////////////////////////
2266     // Oxygen
2267     ////////////////////////////////
2268     if (atom->GetAtomicNum() == 8) {
2269       // 3 neighbours
2270       if (atom->GetExplicitDegree() == 3) {
2271         return 49; // Oxonium oxygen (O+)
2272       }
2273       // 2 neighbours
2274       if (atom->GetExplicitDegree() == 2) {
2275         int hydrogenCount = 0;
2276         FOR_NBORS_OF_ATOM (nbr, atom) {
2277           if (nbr->GetAtomicNum() == OBElements::Hydrogen) {
2278             hydrogenCount++;
2279           }
2280         }
2281 
2282         if (hydrogenCount == 2) {
2283           // H--O--H
2284           return 70; // Oxygen in water (OH2)
2285         }
2286         if (atom->GetExplicitValence() == 3) {
2287           return 51; // Oxenium oxygen (O=+)
2288         }
2289 
2290         return 6; // Generic divalent oxygen, Ether oxygen, Carboxylic acid or ester oxygen,
2291         // Enolic or phenolic oxygen, Oxygen in -O-C=N- moiety, Divalent oxygen in
2292         // thioacid or ester, Divalent nitrate "ether" oxygen, Divalent oxygen in
2293         // sulfate group, Divalent oxygen in sulfite group, One of two divalent
2294         // oxygens attached to sulfur, Divalent oxygen in R(RO)S=O, Other divalent
2295         // oxygen attached to sulfur, Divalent oxygen in phosphate group, Divalent
2296         // oxygen in phosphite group, Divalent oxygen (one of two oxygens attached
2297         // to P), Other divalent oxygen (-O-, OR, OC=O, OC=C, OC=N, OC=S, ONO2,
2298         // ON=O, OSO3, OSO2, OSO, OS=O, -OS, OPO3, OPO2, OPO, -OP)
2299 
2300         // 59 ar
2301       }
2302       // 1 neighbour
2303       if (atom->GetExplicitDegree() == 1) {
2304         oxygenCount = sulphurCount = 0;
2305 
2306         FOR_NBORS_OF_ATOM (nbr, atom) {
2307           bond = _mol.GetBond(&*nbr, atom);
2308 
2309           if (nbr->GetAtomicNum() == OBElements::Carbon || nbr->GetAtomicNum() == OBElements::Nitrogen) {
2310             FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2311               if (nbrNbr->GetAtomicNum() == OBElements::Oxygen && (nbrNbr->GetExplicitDegree() == 1)) {
2312                 oxygenCount++;
2313               }
2314               if (nbrNbr->GetAtomicNum() == OBElements::Sulfur && (nbrNbr->GetExplicitDegree() == 1)) {
2315                 sulphurCount++;
2316               }
2317             }
2318           }
2319           // O---H
2320           if (nbr->GetAtomicNum() == OBElements::Hydrogen) {
2321             return 35;
2322           }
2323           // O-?-C
2324           if (nbr->GetAtomicNum() == OBElements::Carbon) {
2325             if (oxygenCount == 2) {
2326               // O-?-C-?-O
2327               return 32; // Oxygen in carboxylate group (O2CM)
2328             }
2329             if (!bond->IsAromatic() && bond->GetBondOrder() == 1) {
2330               // O--C
2331               return 35; // Oxide oxygen on sp3 carbon, Oxide oxygen on sp2 carbon (OM, OM2)
2332             } else {
2333               // O==C
2334               return 7; // Generic carbonyl oxygen, Carbonyl oxygen in amides,
2335               // Carbonyl oxygen in aldehydes and ketones, Carbonyl
2336               // oxygen in acids or esters (O=C, O=CN, O=CR, O=CO)
2337             }
2338           }
2339           // O-?-N
2340           if (nbr->GetAtomicNum() == OBElements::Nitrogen) {
2341             if (oxygenCount >= 2) {
2342               // O-?-N-?-O
2343               return 32; // Oxygen in nitro group, Nitro-group oxygen in nitrate,
2344               // Nitrate anion oxygen (O2N, O2NO, O3N)
2345             }
2346             if (!bond->IsAromatic() && bond->GetBondOrder() == 1) {
2347 	      if ((nbr->GetExplicitDegree() == 2) || (nbr->GetExplicitValence() == 3))
2348 	      // O(-)--N
2349 	        return 35;
2350 	      else
2351               // O--N
2352                 return 32; // Oxygen in N-oxides (ONX)
2353             } else {
2354               // sometimes ONX bonds are labelled as double bonds
2355               // (e.g. by MOE, noticed by Paolo Tosco)
2356               int ndab = 0;
2357               FOR_BONDS_OF_ATOM(bond2, &*nbr) {
2358                 if (bond2->GetBondOrder() == 2 || bond2->GetBondOrder() == 5)
2359                   ndab++;
2360               }
2361               if (ndab + nbr->GetExplicitDegree() == 5)
2362                 // O--N
2363                 return 32; // Oxygen in N-oxides (ONX)
2364               else
2365                 // O==N
2366                 return 7; // Nitroso oxygen (O=N)
2367             }
2368           }
2369           // O-?-S
2370           if (nbr->GetAtomicNum() == OBElements::Sulfur) {
2371             if (sulphurCount == 1) {
2372               // O1-?-S-?-S1
2373               return 32; // Terminal oxygen in thiosulfinate anion (OSMS)
2374             }
2375             if (!bond->IsAromatic() && bond->GetBondOrder() == 1) {
2376               // O--S
2377               return 32; // Single terminal oxygen on sulfur, One of 2 terminal O's on sulfur,
2378               // One of 3 terminal O's on sulfur, Terminal O in sulfate anion,
2379               // (O-S, O2S, O3S, O4S)
2380             } else {
2381               // O==S
2382 
2383               // are all sulfur nbr atoms carbon?
2384               bool isSulfoxide = true;
2385               int oxygenBoundToSulfur = 0;
2386               FOR_NBORS_OF_ATOM (nbr2, &*nbr) {
2387                 if (atom == &*nbr2)
2388                   continue;
2389 
2390                 if (nbr2->GetAtomicNum() == OBElements::Oxygen)
2391                   ++oxygenBoundToSulfur;
2392               }
2393               FOR_NBORS_OF_ATOM (nbr2, &*nbr) {
2394                 if (atom == &*nbr2)
2395                   continue;
2396                 OBBond* bond = nbr->GetBond(&*nbr2);
2397                 if (!bond->IsAromatic() && bond->GetBondOrder() == 2
2398                   && nbr2->GetAtomicNum() == OBElements::Carbon && oxygenBoundToSulfur == 1)
2399                   isSulfoxide = false; // O=S on sulfur doubly bonded to, e.g., C (O=S=)
2400 
2401                 if ((nbr2->GetAtomicNum() == OBElements::Oxygen && nbr2->GetExplicitDegree() == 1)
2402 		  || (nbr2->GetAtomicNum() == OBElements::Nitrogen && nbr2->GetExplicitDegree() == 2))
2403                   isSulfoxide = false;
2404               }
2405 
2406               if (isSulfoxide)
2407                 return 7; // Doubly bonded sulfoxide oxygen (O=S)
2408               else
2409                 return 32; // (O2S, O2S=C, O3S, O4S)
2410             }
2411           }
2412 
2413           return 32; // Oxygen in phosphine oxide, One of 2 terminal O's on sulfur,
2414           // One of 3 terminal O's on sulfur, One of 4 terminal O's on sulfur,
2415           // Oxygen in perchlorate anion (OP, O2P, O3P, O4P, O4Cl)
2416         }
2417       }
2418     }
2419 
2420     ////////////////////////////////
2421     // Flourine
2422     ////////////////////////////////
2423     if (atom->GetAtomicNum() == 9) {
2424       // 1 neighbour
2425       if (atom->GetExplicitDegree() == 1) {
2426         return 11; // Fluorine (F)
2427       }
2428       // 0 neighbours
2429       if (atom->GetExplicitDegree() == 0) {
2430         return 89; // Fluoride anion (F-)
2431       }
2432     }
2433 
2434     ////////////////////////////////
2435     // Sodium
2436     ////////////////////////////////
2437     if (atom->GetAtomicNum() == 11) {
2438       return 93; // Sodium cation (NA+)
2439     }
2440 
2441     ////////////////////////////////
2442     // Magnesium
2443     ////////////////////////////////
2444     if (atom->GetAtomicNum() == 12) {
2445       return 99; // Dipositive magnesium cation (MG+2)
2446     }
2447 
2448     ////////////////////////////////
2449     // Silicon
2450     ////////////////////////////////
2451     if (atom->GetAtomicNum() == 14) {
2452       return 19; // Silicon (SI)
2453     }
2454 
2455     ////////////////////////////////
2456     // Phosphorus
2457     ////////////////////////////////
2458     if (atom->GetAtomicNum() == 15) {
2459       if (atom->GetExplicitDegree() == 4) {
2460         return 25; // Phosphate group phosphorus, Phosphorus with 3 attached oxygens,
2461         // Phosphorus with 2 attached oxygens, Phosphine oxide phosphorus,
2462         // General tetracoordinate phosphorus (PO4, PO3, PO2, PO, PTET)
2463       }
2464       if (atom->GetExplicitDegree() == 3) {
2465         return 26; // Phosphorus in phosphines (P)
2466       }
2467       if (atom->GetExplicitDegree() == 2) {
2468         return 75; // Phosphorus doubly bonded to C (-P=C)
2469       }
2470     }
2471 
2472     ////////////////////////////////
2473     // Sulfur
2474     ////////////////////////////////
2475     if (atom->GetAtomicNum() == 16) {
2476       // 4 neighbours
2477       if (atom->GetExplicitDegree() == 4) {
2478         return 18; // Sulfone sulfur, Sulfonamide sulfur, Sulfonate group sulfur,
2479         // Sulfate group sulfur, Sulfur in nitrogen analog of sulfone
2480         // (SO2, SO2N, SO3, SO4, SNO)
2481       }
2482       // 3 neighbours
2483       if (atom->GetExplicitDegree() == 3) {
2484         oxygenCount = sulphurCount = doubleBondTo = 0;
2485 
2486         FOR_NBORS_OF_ATOM (nbr, atom) {
2487           bond = _mol.GetBond(&*nbr, atom);
2488           if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
2489             if (nbr->GetAtomicNum() == 6)
2490 	            doubleBondTo = 6;
2491           }
2492 
2493           if (nbr->GetExplicitDegree() == 1) {
2494             if (nbr->GetAtomicNum() == OBElements::Oxygen) {
2495               oxygenCount++;
2496             } else if (nbr->GetAtomicNum() == OBElements::Sulfur) {
2497               sulphurCount++;
2498             }
2499           }
2500         }
2501 
2502         if (oxygenCount == 2) {
2503           if (doubleBondTo == 6) {
2504             return 18; // Sulfone sulfur, doubly bonded to carbon (=SO2)
2505           }
2506           return 73; // Sulfur in anionic sulfinate group (SO2M)
2507         }
2508         if (oxygenCount && sulphurCount)
2509           return 73; // Tricoordinate sulfur in anionic thiosulfinate group (SSOM)
2510 
2511         //if ((doubleBondTo == 6) || (doubleBondTo == 8))
2512         return 17; // Sulfur doubly bonded to carbon, Sulfoxide sulfur (S=C, S=O)
2513       }
2514       // 2 neighbours
2515       if (atom->GetExplicitDegree() == 2) {
2516         doubleBondTo = 0;
2517 
2518         FOR_NBORS_OF_ATOM (nbr, atom) {
2519           if (nbr->GetAtomicNum() == OBElements::Oxygen) {
2520             bond = _mol.GetBond(&*nbr, atom);
2521             if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
2522               doubleBondTo = 8;
2523             }
2524           }
2525         }
2526 
2527         if (doubleBondTo == 8)
2528           return 74; // Sulfinyl sulfur, e.g., in C=S=O (=S=O)
2529 
2530         return 15; // Thiol, sulfide, or disulfide sulfor (S)
2531       }
2532       // 1 neighbour
2533       if (atom->GetExplicitDegree() == 1) {
2534         sulphurCount = doubleBondTo = 0;
2535 
2536         FOR_NBORS_OF_ATOM (nbr, atom) {
2537           FOR_NBORS_OF_ATOM (nbrNbr, &*nbr) {
2538             if (nbrNbr->GetAtomicNum() == OBElements::Sulfur && (nbrNbr->GetExplicitDegree() == 1)) {
2539               sulphurCount++;
2540             }
2541           }
2542           bond = _mol.GetBond(&*nbr, atom);
2543           if (!bond->IsAromatic() && bond->GetBondOrder() == 2) {
2544             doubleBondTo = nbr->GetAtomicNum();
2545           }
2546         }
2547 
2548         if ((doubleBondTo == 6) && (sulphurCount != 2)) {
2549           return 16; // Sulfur doubly bonded to carbon (S=C)
2550         }
2551 
2552         return 72; // Terminal sulfur bonded to P, Anionic terminal sulfur,
2553         // Terminal sulfur in thiosulfinate group (S-P, SM, SSMO)
2554       }
2555 
2556       // 44 ar
2557     }
2558 
2559     ////////////////////////////////
2560     // Clorine
2561     ////////////////////////////////
2562     if (atom->GetAtomicNum() == 17) {
2563       // 4 neighbour
2564       if (atom->GetExplicitDegree() == 4) {
2565         oxygenCount = 0;
2566 
2567         FOR_NBORS_OF_ATOM (nbr, atom) {
2568           if (nbr->GetAtomicNum() == OBElements::Oxygen) {
2569             oxygenCount++;
2570           }
2571         }
2572         if (oxygenCount == 4)
2573           return 77; // Perchlorate anion chlorine (CLO4)
2574       }
2575       // 1 neighbour
2576       if (atom->GetExplicitDegree() == 1) {
2577         return 12; // Chlorine (CL)
2578       }
2579       // 0 neighbours
2580       if (atom->GetExplicitDegree() == 0) {
2581         return 90; // Chloride anion (CL-)
2582       }
2583     }
2584 
2585     ////////////////////////////////
2586     // Potasium
2587     ////////////////////////////////
2588     if (atom->GetAtomicNum() == 19) {
2589       return 94; // Potasium cation (K+)
2590     }
2591 
2592     ////////////////////////////////
2593     // Calcium
2594     ////////////////////////////////
2595     if (atom->GetAtomicNum() == 20) {
2596       // 0 neighbours
2597       if (atom->GetExplicitDegree() == 0) {
2598         return 96; // Dipositive calcium cation (CA+2)
2599       }
2600     }
2601 
2602     ////////////////////////////////
2603     // Iron
2604     ////////////////////////////////
2605     if (atom->GetAtomicNum() == 26) {
2606       if (atom->GetFormalCharge() == 2)
2607         return 87; // Dipositive iron (FE+2)
2608       else
2609         return 88; // Tripositive iron (FE+3)
2610     }
2611 
2612     ////////////////////////////////
2613     // Copper
2614     ////////////////////////////////
2615     if (atom->GetAtomicNum() == 29) {
2616       if (atom->GetFormalCharge() == 1)
2617         return 97; // Monopositive copper cation (CU+1)
2618       else
2619         return 98; // Dipositive copper cation (CU+2)
2620     }
2621 
2622     ////////////////////////////////
2623     // Zinc
2624     ////////////////////////////////
2625     if (atom->GetAtomicNum() == 30) {
2626       return 95; // Dipositive zinc cation (ZN+2)
2627     }
2628 
2629     ////////////////////////////////
2630     // Bromine
2631     ////////////////////////////////
2632     if (atom->GetAtomicNum() == 35) {
2633       // 1 neighbour
2634       if (atom->GetExplicitDegree() == 1) {
2635         return 13; // Bromine (BR)
2636       }
2637       // 0 neighbours
2638       if (atom->GetExplicitDegree() == 0) {
2639         return 91; // Bromide anion (BR-)
2640       }
2641     }
2642 
2643     ////////////////////////////////
2644     // Iodine
2645     ////////////////////////////////
2646     if (atom->GetAtomicNum() == 53) {
2647       // 1 neighbour
2648       if (atom->GetExplicitDegree() == 1) {
2649         return 14; // Iodine (I)
2650       }
2651     }
2652 
2653 
2654 
2655     return 0;
2656   }
2657 
SetTypes()2658   bool OBForceFieldMMFF94::SetTypes()
2659   {
2660     char type[4];
2661 
2662     _mol.SetAtomTypesPerceived();
2663 
2664     // mark all atoms and bonds as non-aromatic
2665     _mol.SetAromaticPerceived();
2666     FOR_BONDS_OF_MOL (bond, _mol)
2667       bond->SetAromatic(false);
2668     FOR_ATOMS_OF_MOL (atom, _mol)
2669       atom->SetAromatic(false);
2670 
2671     // It might be needed to run this function more than once...
2672     bool done = true;
2673     while (done) {
2674       done = PerceiveAromatic();
2675     }
2676 
2677     FOR_ATOMS_OF_MOL (atom, _mol) {
2678       snprintf(type, 3, "%d", GetType(&*atom));
2679       atom->SetType(type);
2680     }
2681 
2682     PrintTypes();
2683 
2684     return true;
2685   }
2686 
SetupCalculations()2687   bool OBForceFieldMMFF94::SetupCalculations()
2688   {
2689     OBFFParameter *parameter;
2690     OBAtom *a, *b, *c, *d;
2691     int type_a, type_b, type_c, type_d;
2692     bool found;
2693     int order;
2694 
2695     IF_OBFF_LOGLVL_LOW
2696       OBFFLog("\nS E T T I N G   U P   C A L C U L A T I O N S\n\n");
2697 
2698     //
2699     // Bond Calculations
2700     //
2701     // no "step-down" procedure
2702     // MMFF part V - page 625 (empirical rule)
2703     //
2704     IF_OBFF_LOGLVL_LOW
2705       OBFFLog("SETTING UP BOND CALCULATIONS...\n");
2706 
2707     OBFFBondCalculationMMFF94 bondcalc;
2708     int bondtype;
2709 
2710     _bondcalculations.clear();
2711 
2712     FOR_BONDS_OF_MOL(bond, _mol) {
2713       a = bond->GetBeginAtom();
2714       b = bond->GetEndAtom();
2715 
2716       // skip this bond if the atoms are ignored
2717       if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) )
2718         continue;
2719 
2720       // if there are any groups specified, check if the two bond atoms are in a single intraGroup
2721       if (HasGroups()) {
2722         bool validBond = false;
2723         for (unsigned int i=0; i < _intraGroup.size(); ++i) {
2724           if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx())) {
2725             validBond = true;
2726             break;
2727           }
2728         }
2729         if (!validBond)
2730           continue;
2731       }
2732 
2733       bondtype = GetBondType(a, b);
2734 
2735       parameter = GetTypedParameter2Atom(bondtype, atoi(a->GetType()), atoi(b->GetType()), _ffbondparams); // from mmffbond.par
2736       if (parameter == nullptr) {
2737         parameter = GetParameter2Atom(a->GetAtomicNum(), b->GetAtomicNum(), _ffbndkparams); // from mmffbndk.par - emperical rules
2738         if (parameter == nullptr) {
2739           IF_OBFF_LOGLVL_LOW {
2740             // This should never happen
2741             snprintf(_logbuf, BUFF_SIZE, "    COULD NOT FIND PARAMETERS FOR BOND %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx());
2742             OBFFLog(_logbuf);
2743           }
2744           return false;
2745         } else {
2746           IF_OBFF_LOGLVL_LOW {
2747             snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR BOND STRETCHING %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx());
2748             OBFFLog(_logbuf);
2749           }
2750 
2751           double rr, rr2, rr4, rr6;
2752           bondcalc.a = a;
2753           bondcalc.b = b;
2754           bondcalc.r0 = GetRuleBondLength(a, b);
2755 
2756           rr = parameter->_dpar[0] / bondcalc.r0; // parameter->_dpar[0]  = r0-ref
2757           rr2 = rr * rr;
2758           rr4 = rr2 * rr2;
2759           rr6 = rr4 * rr2;
2760 
2761           bondcalc.kb = parameter->_dpar[1] * rr6; // parameter->_dpar[1]  = kb-ref
2762           bondcalc.bt = bondtype;
2763           bondcalc.SetupPointers();
2764 
2765           _bondcalculations.push_back(bondcalc);
2766         }
2767       } else {
2768         bondcalc.a = a;
2769         bondcalc.b = b;
2770         bondcalc.kb = parameter->_dpar[0];
2771         bondcalc.r0 = parameter->_dpar[1];
2772         bondcalc.bt = bondtype;
2773         bondcalc.SetupPointers();
2774 
2775         _bondcalculations.push_back(bondcalc);
2776       }
2777     }
2778 
2779     //
2780     // Angle Calculations
2781     //
2782     // MMFF part I - page 513 ("step-down" prodedure)
2783     // MMFF part I - page 519 (reference 68 is actually a footnote)
2784     // MMFF part V - page 627 (empirical rule)
2785     //
2786     // First try and find an exact match, if this fails, step down using the equivalences from mmffdef.par
2787     // five-stage protocol: 1-1-1, 2-2-2, 3-2-3, 4-2-4, 5-2-5
2788     // If this fails, use empirical rules
2789     // Since 1-1-1 = 2-2-2, we will only try 1-1-1 before going to 3-2-3
2790     //
2791     // Stretch-Bend Calculations
2792     //
2793     IF_OBFF_LOGLVL_LOW
2794       OBFFLog("SETTING UP ANGLE & STRETCH-BEND CALCULATIONS...\n");
2795 
2796     OBFFAngleCalculationMMFF94 anglecalc;
2797     OBFFStrBndCalculationMMFF94 strbndcalc;
2798     int angletype, strbndtype, bondtype1, bondtype2;
2799 
2800     _anglecalculations.clear();
2801     _strbndcalculations.clear();
2802 
2803     FOR_ANGLES_OF_MOL(angle, _mol) {
2804       b = _mol.GetAtom((*angle)[0] + 1);
2805       a = _mol.GetAtom((*angle)[1] + 1);
2806       c = _mol.GetAtom((*angle)[2] + 1);
2807 
2808       type_a = atoi(a->GetType());
2809       type_b = atoi(b->GetType());
2810       type_c = atoi(c->GetType());
2811 
2812       // skip this angle if the atoms are ignored
2813       if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) || _constraints.IsIgnored(c->GetIdx()) )
2814         continue;
2815 
2816       // if there are any groups specified, check if the three angle atoms are in a single intraGroup
2817       if (HasGroups()) {
2818         bool validAngle = false;
2819         for (unsigned int i=0; i < _intraGroup.size(); ++i) {
2820           if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) &&
2821               _intraGroup[i].BitIsSet(c->GetIdx())) {
2822             validAngle = true;
2823             break;
2824           }
2825         }
2826         if (!validAngle)
2827           continue;
2828       }
2829 
2830       angletype = GetAngleType(a, b, c);
2831       strbndtype = GetStrBndType(a, b, c);
2832       bondtype1 = GetBondType(a, b);
2833       bondtype2 = GetBondType(b, c);
2834 
2835       if (HasLinSet(type_b)) {
2836         anglecalc.linear = true;
2837       } else {
2838         anglecalc.linear = false;
2839       }
2840 
2841       // try exact match
2842       parameter = GetTypedParameter3Atom(angletype, type_a, type_b, type_c, _ffangleparams);
2843       if (parameter == nullptr) // try 3-2-3
2844         parameter = GetTypedParameter3Atom(angletype, EqLvl3(type_a), type_b, EqLvl3(type_c), _ffangleparams);
2845       if (parameter == nullptr) // try 4-2-4
2846         parameter = GetTypedParameter3Atom(angletype, EqLvl4(type_a), type_b, EqLvl4(type_c), _ffangleparams);
2847       if (parameter == nullptr) // try 5-2-5
2848         parameter = GetTypedParameter3Atom(angletype, EqLvl5(type_a), type_b, EqLvl5(type_c), _ffangleparams);
2849 
2850       if (parameter) {
2851         anglecalc.ka = parameter->_dpar[0];
2852         anglecalc.theta0 = parameter->_dpar[1];
2853         strbndcalc.theta0 = parameter->_dpar[1]; // **
2854       } else {
2855         IF_OBFF_LOGLVL_LOW {
2856           snprintf(_logbuf, BUFF_SIZE, "   USING DEFAULT ANGLE FOR %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx());
2857           snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR ANGLE BENDING %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx());
2858           OBFFLog(_logbuf);
2859         }
2860 
2861         anglecalc.ka = 0.0;
2862         anglecalc.theta0 = 120.0;
2863 
2864         if (GetCrd(type_b) == 4)
2865           anglecalc.theta0 = 109.45;
2866 
2867         if ((GetCrd(type_b) == 2) && b->GetAtomicNum() == OBElements::Oxygen)
2868           anglecalc.theta0 = 105.0;
2869 
2870         if (b->GetAtomicNum() > 10)
2871           anglecalc.theta0 = 95.0;
2872 
2873         if (HasLinSet(type_b))
2874           anglecalc.theta0 = 180.0;
2875 
2876         if ((GetCrd(type_b) == 3) && (GetVal(type_b) == 3) && !GetMltb(type_b)) {
2877           if (b->GetAtomicNum() == OBElements::Nitrogen) {
2878             anglecalc.theta0 = 107.0;
2879           } else {
2880             anglecalc.theta0 = 92.0;
2881           }
2882         }
2883 
2884         if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c))
2885           anglecalc.theta0 = 60.0;
2886 
2887         if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c))
2888           anglecalc.theta0 = 90.0;
2889 
2890         strbndcalc.theta0 = anglecalc.theta0; // **
2891       }
2892 
2893       // empirical rule for 0-b-0 and standard angles
2894       if (anglecalc.ka == 0.0) {
2895         IF_OBFF_LOGLVL_LOW {
2896           snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR ANGLE BENDING FORCE CONSTANT %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx());
2897           OBFFLog(_logbuf);
2898         }
2899 
2900         double beta, Za, Zc, Cb, r0ab, r0bc, theta, theta2, D, rr, rr2;
2901         Za = GetZParam(a);
2902         Cb = GetCParam(b); // Fixed typo -- PR#2741658
2903         Zc = GetZParam(c);
2904 
2905         r0ab = GetBondLength(a, b);
2906         r0bc = GetBondLength(b, c);
2907         rr = r0ab + r0bc;
2908         rr2 = rr * rr;
2909         D = (r0ab - r0bc) / rr2;
2910 
2911         theta = anglecalc.theta0;
2912         theta2 = theta * theta;
2913 
2914         beta = 1.75;
2915         if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c))
2916           beta = 0.85 * beta;
2917         if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c))
2918           beta = 0.05 * beta;
2919 
2920         // Theta2 is in Degrees^2, but parameters are expecting radians
2921         // PR#2741669
2922         anglecalc.ka = (beta * Za * Cb * Zc * exp(-2 * D)) / (rr * theta2 * DEG_TO_RAD * DEG_TO_RAD);
2923       }
2924 
2925       anglecalc.a = a;
2926       anglecalc.b = b;
2927       anglecalc.c = c;
2928       anglecalc.at = angletype;
2929 
2930       anglecalc.SetupPointers();
2931       _anglecalculations.push_back(anglecalc);
2932 
2933       if (anglecalc.linear)
2934         continue;
2935 
2936       parameter = GetTypedParameter3Atom(strbndtype, type_a, type_b, type_c, _ffstrbndparams);
2937       if (parameter == nullptr) {
2938         int rowa, rowb, rowc;
2939 
2940         // This is not a real empirical rule...
2941         //IF_OBFF_LOGLVL_LOW {
2942         //  snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR STRETCH-BENDING FORCE CONSTANT %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx());
2943         //  OBFFLog(_logbuf);
2944         //}
2945 
2946         rowa = GetElementRow(a);
2947         rowb = GetElementRow(b);
2948         rowc = GetElementRow(c);
2949 
2950         parameter = GetParameter3Atom(rowa, rowb, rowc, _ffdfsbparams);
2951 
2952         if (parameter == nullptr) {
2953           // This should never happen
2954           IF_OBFF_LOGLVL_LOW {
2955             snprintf(_logbuf, BUFF_SIZE, "    COULD NOT FIND PARAMETERS FOR STRETCH-BEND %d-%d-%d (IDX)...\n", a->GetIdx(), b->GetIdx(), c->GetIdx());
2956             OBFFLog(_logbuf);
2957           }
2958           return false;
2959         }
2960 
2961         if (rowa == parameter->a) {
2962           strbndcalc.kbaABC = parameter->_dpar[0];
2963           strbndcalc.kbaCBA = parameter->_dpar[1];
2964         } else {
2965           strbndcalc.kbaABC = parameter->_dpar[1];
2966           strbndcalc.kbaCBA = parameter->_dpar[0];
2967         }
2968       } else {
2969         if (type_a == parameter->a) {
2970           strbndcalc.kbaABC = parameter->_dpar[0];
2971           strbndcalc.kbaCBA = parameter->_dpar[1];
2972         } else {
2973           strbndcalc.kbaABC = parameter->_dpar[1];
2974           strbndcalc.kbaCBA = parameter->_dpar[0];
2975         }
2976       }
2977 
2978       strbndcalc.rab0 = GetBondLength(a, b);
2979       strbndcalc.rbc0 = GetBondLength(b ,c);
2980       strbndcalc.a = a;
2981       strbndcalc.b = b;
2982       strbndcalc.c = c;
2983       strbndcalc.sbt = strbndtype;
2984       strbndcalc.SetupPointers();
2985       // Set the pointers to addresses in the anglecalc, find the matching bondcalcs and do the same.
2986       // This should improve performance by not calculating all this twice. We could do the same
2987       // for torsion and angles since the bond lengths are calculated for bond stretching first.
2988       //bool found_angle = false;
2989       /*
2990         for (unsigned int ai = 0; ai < _anglecalculations.size(); ++ai) {
2991         if ( (_anglecalculations[ai].a->GetIdx() == a->GetIdx()) &&
2992         (_anglecalculations[ai].b->GetIdx() == b->GetIdx()) &&
2993         (_anglecalculations[ai].c->GetIdx() == c->GetIdx()) ) {
2994         strbndcalc.theta = &(_anglecalculations[ai].theta);
2995         strbndcalc.force_abc_a = _anglecalculations[ai].force_a;
2996         strbndcalc.force_abc_b = _anglecalculations[ai].force_b;
2997         strbndcalc.force_abc_c = _anglecalculations[ai].force_c;
2998         found_angle = true;
2999         break;
3000         } else if ( (_anglecalculations[ai].a->GetIdx() == c->GetIdx()) &&
3001         (_anglecalculations[ai].b->GetIdx() == b->GetIdx()) &&
3002         (_anglecalculations[ai].c->GetIdx() == a->GetIdx()) ) {
3003         strbndcalc.theta = &(_anglecalculations[ai].theta);
3004         strbndcalc.force_abc_a = _anglecalculations[ai].force_c;
3005         strbndcalc.force_abc_b = _anglecalculations[ai].force_b;
3006         strbndcalc.force_abc_c = _anglecalculations[ai].force_a;
3007         found_angle = true;
3008         break;
3009         }
3010         }
3011       */
3012 
3013       /*
3014         vector<OBFFAngleCalculationMMFF94>::iterator ai;
3015         for (ai = _anglecalculations.begin(); ai != _anglecalculations.end(); ++ai) {
3016         if ( (((*ai).a)->GetIdx() == a->GetIdx()) && (((*ai).b)->GetIdx() == b->GetIdx()) && (((*ai).c)->GetIdx() == c->GetIdx()) ) {
3017         strbndcalc.theta = (*ai).theta;
3018         cout << "theta prt       = " << (*ai).theta << endl;
3019         cout << "delta prt       = " << &((*ai).delta) << endl;
3020         cout << "GetThetaPointer = " << ai->GetThetaPointer() << endl;
3021         strbndcalc.force_abc_a = (*ai).force_a;
3022         strbndcalc.force_abc_b = (*ai).force_b;
3023         strbndcalc.force_abc_c = (*ai).force_c;
3024         found_angle = true;
3025         break;
3026         } else if ( (((*ai).a)->GetIdx() == c->GetIdx()) && (((*ai).b)->GetIdx() == b->GetIdx()) && (((*ai).c)->GetIdx() == a->GetIdx()) ) {
3027         strbndcalc.theta = (*ai).theta;
3028         strbndcalc.force_abc_a = (*ai).force_c;
3029         strbndcalc.force_abc_b = (*ai).force_b;
3030         strbndcalc.force_abc_c = (*ai).force_a;
3031         found_angle = true;
3032         break;
3033         }
3034         }
3035         if (!found_angle) // didn't find matching angle, shouldn't happen, but continue to be safe
3036         continue;
3037 
3038 
3039         bool found_rab = false;
3040         bool found_rbc = false;
3041         vector<OBFFBondCalculationMMFF94>::iterator bi;
3042         for (bi = _bondcalculations.begin(); bi != _bondcalculations.end(); ++bi) {
3043         // find rab
3044         if ( (((*bi).a)->GetIdx() == a->GetIdx()) && (((*bi).b)->GetIdx() == b->GetIdx()) ) {
3045         strbndcalc.rab = &((*bi).rab);
3046         strbndcalc.force_ab_a = (*bi).force_a;
3047         strbndcalc.force_ab_b = (*bi).force_b;
3048         found_rab = true;
3049         } else if ( (((*bi).a)->GetIdx() == b->GetIdx()) && (((*bi).b)->GetIdx() == a->GetIdx()) ) {
3050         strbndcalc.rab = &((*bi).rab);
3051         strbndcalc.force_ab_a = (*bi).force_b;
3052         strbndcalc.force_ab_b = (*bi).force_a;
3053         found_rab = true;
3054         }
3055         // find rbc
3056         if ( (((*bi).a)->GetIdx() == b->GetIdx()) && (((*bi).b)->GetIdx() == c->GetIdx()) ) {
3057         strbndcalc.rbc = &(bondcalc.rab);
3058         strbndcalc.force_ab_a = (*bi).force_a;
3059         strbndcalc.force_ab_b = (*bi).force_b;
3060         found_rbc = true;
3061         } else if ( (((*bi).a)->GetIdx() == c->GetIdx()) && (((*bi).b)->GetIdx() == b->GetIdx()) ) {
3062         strbndcalc.rbc = &(bondcalc.rab);
3063         strbndcalc.force_ab_a = (*bi).force_b;
3064         strbndcalc.force_ab_b = (*bi).force_a;
3065         found_rbc = true;
3066         }
3067 
3068         if (found_rab && found_rbc)
3069         break;
3070         }
3071         if (!found_rab || !found_rbc) // didn't find matching bond, or atoms overlap
3072         continue;
3073       */
3074       _strbndcalculations.push_back(strbndcalc);
3075 
3076     }
3077 
3078     //
3079     // Torsion Calculations
3080     //
3081     // MMFF part I - page 513 ("step-down" prodedure)
3082     // MMFF part I - page 519 (reference 68 is actually a footnote)
3083     // MMFF part IV - page 631 (empirical rule)
3084     //
3085     // First try and find an exact match, if this fails, step down using the equivalences from mmffdef.par
3086     // five-stage protocol: 1-1-1-1, 2-2-2-2, 3-2-2-5, 5-2-2-3, 5-2-2-5
3087     // If this fails, use empirical rules
3088     // Since 1-1-1-1 = 2-2-2-2, we will only try 1-1-1-1 before going to 3-2-2-5
3089     //
3090     IF_OBFF_LOGLVL_LOW
3091       OBFFLog("SETTING UP TORSION CALCULATIONS...\n");
3092 
3093     OBFFTorsionCalculationMMFF94 torsioncalc;
3094     int torsiontype;
3095 
3096     _torsioncalculations.clear();
3097 
3098     FOR_TORSIONS_OF_MOL(t, _mol) {
3099       a = _mol.GetAtom((*t)[0] + 1);
3100       b = _mol.GetAtom((*t)[1] + 1);
3101       c = _mol.GetAtom((*t)[2] + 1);
3102       d = _mol.GetAtom((*t)[3] + 1);
3103 
3104       type_a = atoi(a->GetType());
3105       type_b = atoi(b->GetType());
3106       type_c = atoi(c->GetType());
3107       type_d = atoi(d->GetType());
3108 
3109       // skip this torsion if the atoms are ignored
3110       if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) ||
3111            _constraints.IsIgnored(c->GetIdx()) || _constraints.IsIgnored(d->GetIdx()) )
3112         continue;
3113 
3114       // if there are any groups specified, check if the four torsion atoms are in a single intraGroup
3115       if (HasGroups()) {
3116         bool validTorsion = false;
3117         for (unsigned int i=0; i < _intraGroup.size(); ++i) {
3118           if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) &&
3119               _intraGroup[i].BitIsSet(c->GetIdx()) && _intraGroup[i].BitIsSet(d->GetIdx())) {
3120             validTorsion = true;
3121             break;
3122           }
3123         }
3124         if (!validTorsion)
3125           continue;
3126       }
3127 
3128       torsiontype = GetTorsionType(a, b, c, d);
3129       // CXT = MC*(J*MA**3 + K*MA**2 + I*MA + L) + TTijkl  MC = 6, MA = 136
3130       order = (type_c*2515456 + type_b*18496 + type_d*136 + type_a)
3131         - (type_b*2515456 + type_c*18496 + type_a*136 + type_d);
3132 
3133       if (order >= 0) {
3134         // try exact match
3135         parameter = GetTypedParameter4Atom(torsiontype, type_a, type_b, type_c, type_d, _fftorsionparams);
3136         if (parameter == nullptr) // try 3-2-2-5
3137           parameter = GetTypedParameter4Atom(torsiontype, EqLvl3(type_a), type_b, type_c, EqLvl5(type_d), _fftorsionparams);
3138         if (parameter == nullptr) // try 5-2-2-3
3139           parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_a), type_b, type_c, EqLvl3(type_d), _fftorsionparams);
3140         if (parameter == nullptr) // try 5-2-2-5
3141           parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_a), type_b, type_c, EqLvl5(type_d), _fftorsionparams);
3142       } else {
3143         // try exact match
3144         parameter = GetTypedParameter4Atom(torsiontype, type_d, type_c, type_b, type_a, _fftorsionparams);
3145         if (parameter == nullptr) // try 3-2-2-5
3146           parameter = GetTypedParameter4Atom(torsiontype, EqLvl3(type_d), type_c, type_b, EqLvl5(type_a), _fftorsionparams);
3147         if (parameter == nullptr) // try 5-2-2-3
3148           parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_d), type_c, type_b, EqLvl3(type_a), _fftorsionparams);
3149         if (parameter == nullptr) // try 5-2-2-5
3150           parameter = GetTypedParameter4Atom(torsiontype, EqLvl5(type_d), type_c, type_b, EqLvl5(type_a), _fftorsionparams);
3151       }
3152 
3153       if (parameter) {
3154         torsioncalc.v1 = parameter->_dpar[0];
3155         torsioncalc.v2 = parameter->_dpar[1];
3156         torsioncalc.v3 = parameter->_dpar[2];
3157       } else {
3158         bool found_rule = false;
3159 
3160         //IF_OBFF_LOGLVL_LOW {
3161         //  snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR TORSION FORCE CONSTANT %d-%d-%d-%d (IDX)...\n",
3162         //    a->GetIdx(), b->GetIdx(), c->GetIdx(), d->GetIdx());
3163         //  OBFFLog(_logbuf);
3164         //}
3165 
3166         // rule (a) page 631
3167         if (HasLinSet(type_b) || HasLinSet(type_c))
3168           continue;
3169 
3170         // rule (b) page 631
3171         if (b->GetBond(c)->IsAromatic()) {
3172           double Ub, Uc, pi_bc, beta;
3173           Ub = GetUParam(b);
3174           Uc = GetUParam(c);
3175 
3176           if (!HasPilpSet(type_b) && !HasPilpSet(type_c))
3177             pi_bc = 0.5;
3178           else
3179             pi_bc = 0.3;
3180 
3181           if (((GetVal(type_b) == 3) && (GetVal(type_c) == 4)) ||
3182               ((GetVal(type_b) == 4) && (GetVal(type_c) == 3)))
3183             beta = 3.0;
3184           else
3185             beta = 6.0;
3186 
3187           torsioncalc.v1 = 0.0;
3188           torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc);
3189           torsioncalc.v3 = 0.0;
3190           found_rule = true;
3191         } else {
3192           // rule (c) page 631
3193        	  double Ub, Uc, pi_bc, beta;
3194           Ub = GetUParam(b);
3195           Uc = GetUParam(c);
3196           OBBond *bond = a->GetBond(b);
3197           if (((GetMltb(type_b) == 2) && (GetMltb(type_c) == 2)) && !bond->IsAromatic() && bond->GetBondOrder() == 2)
3198             pi_bc = 1.0;
3199           else
3200             pi_bc = 0.4;
3201 
3202           beta = 6.0;
3203           torsioncalc.v1 = 0.0;
3204           torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc);
3205           torsioncalc.v3 = 0.0;
3206           found_rule = true;
3207         }
3208 
3209         // rule (d) page 632
3210         if (!found_rule)
3211           if (((GetCrd(type_b) == 4) && (GetCrd(type_c) == 4))) {
3212             double Vb, Vc;
3213             Vb = GetVParam(b);
3214             Vc = GetVParam(c);
3215 
3216             torsioncalc.v1 = 0.0;
3217             torsioncalc.v2 = 0.0;
3218             torsioncalc.v3 = sqrt(Vb * Vc) / 9.0;
3219             found_rule = true;
3220           }
3221 
3222         // rule (e) page 632
3223         if (!found_rule)
3224           if (((GetCrd(type_b) == 4) && (GetCrd(type_c) != 4))) {
3225             if (GetCrd(type_c) == 3) // case (1)
3226               if ((GetVal(type_c) == 4) || (GetVal(type_c) == 34) || (GetMltb(type_c) != 0))
3227                 continue;
3228 
3229             if (GetCrd(type_c) == 2) // case (2)
3230               if ((GetVal(type_c) == 3) || (GetMltb(type_c) != 0))
3231                 continue;
3232 
3233             // case (3) saturated bonds -- see rule (h)
3234           }
3235 
3236         // rule (f) page 632
3237         if (!found_rule)
3238           if (((GetCrd(type_b) != 4) && (GetCrd(type_c) == 4))) {
3239             if (GetCrd(type_b) == 3) // case (1)
3240               if ((GetVal(type_b) == 4) || (GetVal(type_b) == 34) || (GetMltb(type_b) != 0))
3241                 continue;
3242 
3243             if (GetCrd(type_b) == 2) // case (2)
3244               if ((GetVal(type_b) == 3) || (GetMltb(type_b) != 0))
3245                 continue;
3246 
3247             // case (3) saturated bonds
3248           }
3249 
3250         // rule (g) page 632
3251         if (!found_rule) {
3252           OBBond *bond = b->GetBond(c);
3253           if (!bond->IsAromatic() && bond->GetBondOrder() == 1 && (
3254             (GetMltb(type_b) && GetMltb(type_c)) ||
3255             (GetMltb(type_b) && HasPilpSet(type_c)) ||
3256             (GetMltb(type_c) && HasPilpSet(type_b)))) {
3257             if (HasPilpSet(type_b) && HasPilpSet(type_c)) // case (1)
3258               continue;
3259 
3260             double Ub, Uc, pi_bc, beta;
3261             Ub = GetUParam(b);
3262             Uc = GetUParam(c);
3263             beta = 6.0;
3264 
3265             if (HasPilpSet(type_b) && GetMltb(type_c)) { // case (2)
3266               if (GetMltb(type_c) == 1)
3267                 pi_bc = 0.5;
3268               else if ((GetElementRow(b) == 1) && (GetElementRow(c) == 1))
3269                 pi_bc = 0.3;
3270               else
3271                 pi_bc = 0.15;
3272               found_rule = true;
3273             }
3274 
3275             if (HasPilpSet(type_c) && GetMltb(type_b)) { // case (3)
3276               if (GetMltb(type_b) == 1)
3277                 pi_bc = 0.5;
3278               else if ((GetElementRow(b) == 1) && (GetElementRow(c) == 1))
3279                 pi_bc = 0.3;
3280               else
3281                 pi_bc = 0.15;
3282               found_rule = true;
3283             }
3284 
3285             if (!found_rule)
3286               if (((GetMltb(type_b) == 1) || (GetMltb(type_c) == 1)) && (b->GetAtomicNum() != OBElements::Carbon || c->GetAtomicNum() != OBElements::Carbon)) {
3287                 pi_bc = 0.4;
3288                 found_rule = true;
3289               }
3290 
3291             if (!found_rule)
3292               pi_bc = 0.15;
3293 
3294             torsioncalc.v1 = 0.0;
3295             torsioncalc.v2 = beta * pi_bc * sqrt(Ub * Uc);
3296             torsioncalc.v3 = 0.0;
3297             found_rule = true;
3298           }
3299         }
3300 
3301         // rule (h) page 632
3302         if (!found_rule) {
3303           if ((b->GetAtomicNum() == OBElements::Oxygen || b->GetAtomicNum() == OBElements::Sulfur) && (c->GetAtomicNum() == OBElements::Oxygen || c->GetAtomicNum() == OBElements::Sulfur)) {
3304             double Wb, Wc;
3305 
3306             if (b->GetAtomicNum() == OBElements::Oxygen) {
3307               Wb = 2.0;
3308             }
3309             else {
3310               Wb = 8.0;
3311             }
3312 
3313             if (c->GetAtomicNum() == OBElements::Oxygen) {
3314               Wc = 2.0;
3315             }
3316             else {
3317               Wc = 8.0;
3318             }
3319 
3320             torsioncalc.v1 = 0.0;
3321             torsioncalc.v2 = -sqrt(Wb * Wc);
3322             torsioncalc.v3 = 0.0;
3323           } else {
3324             double Vb, Vc, Nbc;
3325             Vb = GetVParam(b);
3326             Vc = GetVParam(c);
3327 
3328             IF_OBFF_LOGLVL_LOW {
3329               snprintf(_logbuf, BUFF_SIZE, "   USING EMPIRICAL RULE FOR TORSION FORCE CONSTANT %d-%d-%d-%d (IDX)...\n",
3330                       a->GetIdx(), b->GetIdx(), c->GetIdx(), d->GetIdx());
3331               OBFFLog(_logbuf);
3332             }
3333 
3334             Nbc = GetCrd(type_b) * GetCrd(type_c);
3335 
3336             torsioncalc.v1 = 0.0;
3337             torsioncalc.v2 = 0.0;
3338             torsioncalc.v3 = sqrt(Vb * Vc) / Nbc;
3339           }
3340         }
3341       }
3342 
3343       torsioncalc.a = a;
3344       torsioncalc.b = b;
3345       torsioncalc.c = c;
3346       torsioncalc.d = d;
3347       torsioncalc.SetupPointers();
3348       torsioncalc.tt = torsiontype;
3349 
3350       _torsioncalculations.push_back(torsioncalc);
3351     }
3352 
3353     //
3354     // Out-Of-Plane Calculations
3355     //
3356     IF_OBFF_LOGLVL_LOW
3357       OBFFLog("SETTING UP OOP CALCULATIONS...\n");
3358 
3359     OBFFOOPCalculationMMFF94 oopcalc;
3360 
3361     _oopcalculations.clear();
3362 
3363     FOR_ATOMS_OF_MOL(atom, _mol) {
3364       b = (OBAtom*) &*atom;
3365 
3366       found = false;
3367 
3368       type_b = atoi(b->GetType());
3369 
3370       for (unsigned int idx=0; idx < _ffoopparams.size(); idx++) {
3371         if (type_b == _ffoopparams[idx].b) {
3372           a = nullptr;
3373           c = nullptr;
3374           d = nullptr;
3375 
3376           FOR_NBORS_OF_ATOM(nbr, b) {
3377             if (a ==nullptr)
3378               a = (OBAtom*) &*nbr;
3379             else if (c == nullptr)
3380               c = (OBAtom*) &*nbr;
3381             else
3382               d = (OBAtom*) &*nbr;
3383           }
3384 
3385           if (a == nullptr || c == nullptr || d == nullptr)
3386             break;
3387 
3388           type_a = atoi(a->GetType());
3389           type_c = atoi(c->GetType());
3390           type_d = atoi(d->GetType());
3391 
3392           // skip this oop if the atoms are ignored
3393           if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) ||
3394                _constraints.IsIgnored(c->GetIdx()) || _constraints.IsIgnored(d->GetIdx()) )
3395             continue;
3396 
3397           // if there are any groups specified, check if the four oop atoms are in a single intraGroup
3398           if (HasGroups()) {
3399             bool validOOP = false;
3400             for (unsigned int i=0; i < _intraGroup.size(); ++i) {
3401               if (_intraGroup[i].BitIsSet(a->GetIdx()) && _intraGroup[i].BitIsSet(b->GetIdx()) &&
3402                   _intraGroup[i].BitIsSet(c->GetIdx()) && _intraGroup[i].BitIsSet(d->GetIdx())) {
3403                 validOOP = true;
3404                 break;
3405               }
3406             }
3407             if (!validOOP)
3408               continue;
3409           }
3410 
3411           if (((type_a == _ffoopparams[idx].a) && (type_c == _ffoopparams[idx].c) && (type_d == _ffoopparams[idx].d)) ||
3412               ((type_c == _ffoopparams[idx].a) && (type_a == _ffoopparams[idx].c) && (type_d == _ffoopparams[idx].d)) ||
3413               ((type_c == _ffoopparams[idx].a) && (type_d == _ffoopparams[idx].c) && (type_a == _ffoopparams[idx].d)) ||
3414               ((type_d == _ffoopparams[idx].a) && (type_c == _ffoopparams[idx].c) && (type_a == _ffoopparams[idx].d)) ||
3415               ((type_a == _ffoopparams[idx].a) && (type_d == _ffoopparams[idx].c) && (type_c == _ffoopparams[idx].d)) ||
3416               ((type_d == _ffoopparams[idx].a) && (type_a == _ffoopparams[idx].c) && (type_c == _ffoopparams[idx].d)))
3417             {
3418               found = true;
3419 
3420               oopcalc.koop = _ffoopparams[idx]._dpar[0];
3421 
3422               // A-B-CD || C-B-AD  PLANE = ABC
3423               oopcalc.a = a;
3424               oopcalc.b = b;
3425               oopcalc.c = c;
3426               oopcalc.d = d;
3427 
3428               oopcalc.SetupPointers();
3429               _oopcalculations.push_back(oopcalc);
3430 
3431               // C-B-DA || D-B-CA  PLANE BCD
3432               oopcalc.a = d;
3433               oopcalc.d = a;
3434 
3435               oopcalc.SetupPointers();
3436               _oopcalculations.push_back(oopcalc);
3437 
3438               // A-B-DC || D-B-AC  PLANE ABD
3439               oopcalc.a = a;
3440               oopcalc.c = d;
3441               oopcalc.d = c;
3442 
3443               oopcalc.SetupPointers();
3444               _oopcalculations.push_back(oopcalc);
3445             }
3446 
3447           if ((_ffoopparams[idx].a == 0) && (_ffoopparams[idx].c == 0) && (_ffoopparams[idx].d == 0) && !found) // *-XX-*-*
3448             {
3449               oopcalc.koop = _ffoopparams[idx]._dpar[0];
3450 
3451               // A-B-CD || C-B-AD  PLANE = ABC
3452               oopcalc.a = a;
3453               oopcalc.b = b;
3454               oopcalc.c = c;
3455               oopcalc.d = d;
3456 
3457               oopcalc.SetupPointers();
3458               _oopcalculations.push_back(oopcalc);
3459 
3460               // C-B-DA || D-B-CA  PLANE BCD
3461               oopcalc.a = d;
3462               oopcalc.d = a;
3463 
3464               oopcalc.SetupPointers();
3465               _oopcalculations.push_back(oopcalc);
3466 
3467               // A-B-DC || D-B-AC  PLANE ABD
3468               oopcalc.a = a;
3469               oopcalc.c = d;
3470               oopcalc.d = c;
3471 
3472               oopcalc.SetupPointers();
3473               _oopcalculations.push_back(oopcalc);
3474             }
3475         }
3476       }
3477     }
3478 
3479     //
3480     // VDW Calculations
3481     //
3482     IF_OBFF_LOGLVL_LOW
3483       OBFFLog("SETTING UP VAN DER WAALS CALCULATIONS...\n");
3484 
3485     OBFFVDWCalculationMMFF94 vdwcalc;
3486 
3487     _vdwcalculations.clear();
3488 
3489     int pairIndex = -1;
3490     FOR_PAIRS_OF_MOL(p, _mol) {
3491       ++pairIndex;
3492       a = _mol.GetAtom((*p)[0]);
3493       b = _mol.GetAtom((*p)[1]);
3494 
3495       // skip this vdw if the atoms are ignored
3496       if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) )
3497         continue;
3498 
3499       // if there are any groups specified, check if the two atoms are in a single _interGroup or if
3500       // two two atoms are in one of the _interGroups pairs.
3501       if (HasGroups()) {
3502         bool validVDW = false;
3503         for (unsigned int i=0; i < _interGroup.size(); ++i) {
3504           if (_interGroup[i].BitIsSet(a->GetIdx()) && _interGroup[i].BitIsSet(b->GetIdx())) {
3505             validVDW = true;
3506             break;
3507           }
3508         }
3509         if (!validVDW) {
3510           for (unsigned int i=0; i < _interGroups.size(); ++i) {
3511             if (_interGroups[i].first.BitIsSet(a->GetIdx()) && _interGroups[i].second.BitIsSet(b->GetIdx())) {
3512               validVDW = true;
3513               break;
3514             }
3515             if (_interGroups[i].first.BitIsSet(b->GetIdx()) && _interGroups[i].second.BitIsSet(a->GetIdx())) {
3516               validVDW = true;
3517               break;
3518             }
3519           }
3520         }
3521 
3522         if (!validVDW)
3523           continue;
3524       }
3525 
3526       OBFFParameter *parameter_a, *parameter_b;
3527       parameter_a = GetParameter1Atom(atoi(a->GetType()), _ffvdwparams);
3528       parameter_b = GetParameter1Atom(atoi(b->GetType()), _ffvdwparams);
3529       if (parameter_a == nullptr || parameter_b == nullptr) {
3530         IF_OBFF_LOGLVL_LOW {
3531           snprintf(_logbuf, BUFF_SIZE, "   COULD NOT FIND VAN DER WAALS PARAMETERS FOR %d-%d (IDX)...\n", a->GetIdx(), b->GetIdx());
3532           OBFFLog(_logbuf);
3533         }
3534 
3535         return false;
3536       }
3537 
3538       vdwcalc.a = a;
3539       vdwcalc.alpha_a = parameter_a->_dpar[0];
3540       vdwcalc.Na = parameter_a->_dpar[1];
3541       vdwcalc.Aa = parameter_a->_dpar[2];
3542       vdwcalc.Ga = parameter_a->_dpar[3];
3543       vdwcalc.aDA = parameter_a->_ipar[0];
3544 
3545       vdwcalc.b = b;
3546       vdwcalc.alpha_b = parameter_b->_dpar[0];
3547       vdwcalc.Nb = parameter_b->_dpar[1];
3548       vdwcalc.Ab = parameter_b->_dpar[2];
3549       vdwcalc.Gb = parameter_b->_dpar[3];
3550       vdwcalc.bDA = parameter_b->_ipar[0];
3551 
3552       //these calculations only need to be done once for each pair,
3553       //we do them now and save them for later use
3554       double R_AA, R_BB, R_AB6, g_AB, g_AB2;
3555       double R_AB2, R_AB4, /*R_AB7,*/ sqrt_a, sqrt_b;
3556 
3557       R_AA = vdwcalc.Aa * pow(vdwcalc.alpha_a, 0.25);
3558       R_BB = vdwcalc.Ab * pow(vdwcalc.alpha_b, 0.25);
3559       sqrt_a = sqrt(vdwcalc.alpha_a / vdwcalc.Na);
3560       sqrt_b = sqrt(vdwcalc.alpha_b / vdwcalc.Nb);
3561 
3562       if (vdwcalc.aDA == 1) { // hydrogen bond donor
3563         vdwcalc.R_AB = 0.5 * (R_AA + R_BB);
3564         R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB;
3565         R_AB4 = R_AB2 * R_AB2;
3566         R_AB6 = R_AB4 * R_AB2;
3567 
3568         if (vdwcalc.bDA == 2) { // hydrogen bond acceptor
3569           vdwcalc.epsilon = 0.5 * (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6);
3570           // R_AB is scaled to 0.8 for D-A interactions. The value used in the calculation of epsilon is not scaled.
3571           vdwcalc.R_AB = 0.8 * vdwcalc.R_AB;
3572         } else
3573           vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6);
3574 
3575         R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB;
3576         R_AB4 = R_AB2 * R_AB2;
3577         R_AB6 = R_AB4 * R_AB2;
3578         vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB;
3579       } else if (vdwcalc.bDA == 1) { // hydrogen bond donor
3580         vdwcalc.R_AB = 0.5 * (R_AA + R_BB);
3581        	R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB;
3582         R_AB4 = R_AB2 * R_AB2;
3583         R_AB6 = R_AB4 * R_AB2;
3584 
3585         if (vdwcalc.aDA == 2) { // hydrogen bond acceptor
3586           vdwcalc.epsilon = 0.5 * (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6);
3587           // R_AB is scaled to 0.8 for D-A interactions. The value used in the calculation of epsilon is not scaled.
3588           vdwcalc.R_AB = 0.8 * vdwcalc.R_AB;
3589         } else
3590           vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6);
3591 
3592         R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB;
3593         R_AB4 = R_AB2 * R_AB2;
3594         R_AB6 = R_AB4 * R_AB2;
3595         vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB;
3596       } else {
3597         g_AB = (R_AA - R_BB) / ( R_AA + R_BB);
3598         g_AB2 = g_AB * g_AB;
3599         vdwcalc.R_AB =  0.5 * (R_AA + R_BB) * (1.0 + 0.2 * (1.0 - exp(-12.0 * g_AB2)));
3600         R_AB2 = vdwcalc.R_AB * vdwcalc.R_AB;
3601         R_AB4 = R_AB2 * R_AB2;
3602         R_AB6 = R_AB4 * R_AB2;
3603         vdwcalc.R_AB7 = R_AB6 * vdwcalc.R_AB;
3604         vdwcalc.epsilon = (181.16 * vdwcalc.Ga * vdwcalc.Gb * vdwcalc.alpha_a * vdwcalc.alpha_b) / (sqrt_a + sqrt_b) * (1.0 / R_AB6);
3605       }
3606 
3607       vdwcalc.pairIndex = pairIndex;
3608       vdwcalc.SetupPointers();
3609       _vdwcalculations.push_back(vdwcalc);
3610     }
3611 
3612     //
3613     // Electrostatic Calculations
3614     //
3615     IF_OBFF_LOGLVL_LOW
3616       OBFFLog("SETTING UP ELECTROSTATIC CALCULATIONS...\n");
3617 
3618     OBFFElectrostaticCalculationMMFF94 elecalc;
3619 
3620     _electrostaticcalculations.clear();
3621 
3622     pairIndex = -1;
3623     FOR_PAIRS_OF_MOL(p, _mol) {
3624       ++pairIndex;
3625       a = _mol.GetAtom((*p)[0]);
3626       b = _mol.GetAtom((*p)[1]);
3627 
3628       // skip this ele if the atoms are ignored
3629       if ( _constraints.IsIgnored(a->GetIdx()) || _constraints.IsIgnored(b->GetIdx()) )
3630         continue;
3631 
3632       // if there are any groups specified, check if the two atoms are in a single _interGroup or if
3633       // two two atoms are in one of the _interGroups pairs.
3634       if (HasGroups()) {
3635         bool validEle = false;
3636         for (unsigned int i=0; i < _interGroup.size(); ++i) {
3637           if (_interGroup[i].BitIsSet(a->GetIdx()) && _interGroup[i].BitIsSet(b->GetIdx())) {
3638             validEle = true;
3639             break;
3640           }
3641         }
3642         if (!validEle) {
3643           for (unsigned int i=0; i < _interGroups.size(); ++i) {
3644             if (_interGroups[i].first.BitIsSet(a->GetIdx()) && _interGroups[i].second.BitIsSet(b->GetIdx())) {
3645               validEle = true;
3646               break;
3647             }
3648             if (_interGroups[i].first.BitIsSet(b->GetIdx()) && _interGroups[i].second.BitIsSet(a->GetIdx())) {
3649               validEle = true;
3650               break;
3651             }
3652           }
3653         }
3654 
3655         if (!validEle)
3656           continue;
3657       }
3658 
3659       elecalc.qq = 332.0716 * a->GetPartialCharge() * b->GetPartialCharge() / _epsilon;
3660 
3661       if (elecalc.qq) {
3662         elecalc.a = &*a;
3663         elecalc.b = &*b;
3664 
3665         // 1-4 scaling
3666         if (a->IsOneFour(b))
3667           elecalc.qq *= 0.75;
3668 
3669         elecalc.pairIndex = pairIndex;
3670         elecalc.SetupPointers();
3671         _electrostaticcalculations.push_back(elecalc);
3672       }
3673     }
3674 
3675     return true;
3676   }
3677 
SetupPointers()3678   bool OBForceFieldMMFF94::SetupPointers()
3679   {
3680     for (unsigned int i = 0; i < _bondcalculations.size(); ++i)
3681       _bondcalculations[i].SetupPointers();
3682     for (unsigned int i = 0; i < _anglecalculations.size(); ++i)
3683       _anglecalculations[i].SetupPointers();
3684     for (unsigned int i = 0; i < _strbndcalculations.size(); ++i)
3685       _strbndcalculations[i].SetupPointers();
3686     for (unsigned int i = 0; i < _torsioncalculations.size(); ++i)
3687       _torsioncalculations[i].SetupPointers();
3688     for (unsigned int i = 0; i < _oopcalculations.size(); ++i)
3689       _oopcalculations[i].SetupPointers();
3690     for (unsigned int i = 0; i < _vdwcalculations.size(); ++i)
3691       _vdwcalculations[i].SetupPointers();
3692     for (unsigned int i = 0; i < _electrostaticcalculations.size(); ++i)
3693       _electrostaticcalculations[i].SetupPointers();
3694 
3695     return true;
3696   }
3697 
3698 
3699   // we set the the formal charge with SetPartialCharge because formal charges
3700   // in MMFF94 are not always and integer
SetFormalCharges()3701   bool OBForceFieldMMFF94::SetFormalCharges()
3702   {
3703     _mol.SetAutomaticPartialCharge(false);
3704 
3705     FOR_ATOMS_OF_MOL (atom, _mol) {
3706       int type = atoi(atom->GetType());
3707       atom->SetPartialCharge(0.0);
3708 
3709       bool done = false;
3710       switch (type) {
3711       case 34:
3712       case 49:
3713       case 51:
3714       case 54:
3715       case 58:
3716       case 92:
3717       case 93:
3718       case 94:
3719       case 97:
3720         atom->SetPartialCharge(1.0);
3721         done = true;
3722         break;
3723       case 35:
3724       case 62:
3725       case 89:
3726       case 90:
3727       case 91:
3728         atom->SetPartialCharge(-1.0);
3729         done = true;
3730         break;
3731       case 55:
3732         atom->SetPartialCharge(0.5);
3733         done = true;
3734         break;
3735       case 87:
3736       case 95:
3737       case 96:
3738       case 98:
3739       case 99:
3740         atom->SetPartialCharge(2.0);
3741         done = true;
3742         break;
3743       case 88:
3744         atom->SetPartialCharge(3.0);
3745         done = true;
3746         break;
3747         //case 98:
3748         //  atom->SetPartialCharge(3.0);
3749       default:
3750         break;
3751       }
3752 
3753       if (done)
3754         continue;
3755 
3756       if (type == 56) {
3757         int n_count = 0;
3758         int temp_type;
3759         FOR_ATOMS_OF_MOL (atom2, _mol) {
3760           temp_type = atoi(atom2->GetType());
3761           if (temp_type == 56 || temp_type == 81)
3762             ++n_count;
3763         }
3764         atom->SetPartialCharge((double)((n_count + 1) / 3) / (double)n_count);
3765       } else if (type == 32) {
3766         int o_count = 0;
3767         bool sulfonamide = false;
3768         bool sulfone_s_c = false;
3769         int s_count = 0;
3770 
3771         FOR_NBORS_OF_ATOM(nbr, &*atom) {
3772           FOR_NBORS_OF_ATOM(nbr2, &*nbr) {
3773             if (nbr2->GetAtomicNum() == OBElements::Oxygen && (nbr2->GetExplicitDegree() == 1))
3774               o_count++;
3775             if (nbr2->GetAtomicNum() == OBElements::Sulfur && (nbr2->GetExplicitDegree() == 1))
3776               s_count++;
3777             if (nbr2->GetAtomicNum() == OBElements::Nitrogen && !nbr2->IsAromatic())
3778               sulfonamide = true;
3779             OBBond *bond = nbr->GetBond(&*nbr2);
3780             if (nbr2->GetAtomicNum() == OBElements::Carbon && !bond->IsAromatic() && bond->GetBondOrder() == 2)
3781               sulfone_s_c = true;
3782           }
3783 
3784           if (nbr->GetAtomicNum() == OBElements::Carbon)
3785             atom->SetPartialCharge(-0.5); // O2CM
3786 
3787           if (nbr->GetAtomicNum() == OBElements::Nitrogen && (o_count == 3))
3788             atom->SetPartialCharge(-1.0 / o_count);  // O3N
3789 
3790           if (nbr->GetAtomicNum() == OBElements::Sulfur && !sulfonamide) {
3791             if (((o_count + s_count) == 2) && (nbr->GetExplicitDegree() == 3)
3792                 && (nbr->GetExplicitValence() >= 3) && !sulfone_s_c) {
3793               atom->SetPartialCharge(-0.5); // O2S (sulfinate)
3794             }
3795             else if ((o_count + s_count) == 3) {
3796               atom->SetPartialCharge(-1.0 / 3.0); // O3S
3797             }
3798             else if ((o_count + s_count) == 4) {
3799               atom->SetPartialCharge(-0.5); // O4S
3800             }
3801           }
3802 
3803           if (nbr->GetAtomicNum() == OBElements::Phosphorus) {
3804             if ((o_count + s_count) == 2) {
3805               atom->SetPartialCharge(-0.5); // O2P
3806             }
3807             else if ((o_count + s_count) == 3) {
3808               atom->SetPartialCharge(-2.0 / 3.0); // O3P
3809             }
3810             else if ((o_count + s_count) == 4) {
3811               atom->SetPartialCharge(-0.25); // O4P
3812             }
3813           }
3814 
3815           if (atoi(nbr->GetType()) == 77)
3816             atom->SetPartialCharge(-0.25); // O4CL
3817         }
3818       } else if (type == 61) {
3819         FOR_BONDS_OF_ATOM(bond, &*atom) {
3820           OBAtom *nbr = bond->GetNbrAtom(&*atom);
3821           if (!bond->IsAromatic() && bond->GetBondOrder() == 3 && nbr->GetAtomicNum() == OBElements::Nitrogen)
3822             atom->SetPartialCharge(1.0);
3823         }
3824       } else if (type == 72) {
3825         int s_count = 0;
3826 
3827         FOR_NBORS_OF_ATOM(nbr, &*atom) {
3828           if (nbr->GetAtomicNum() == OBElements::Sulfur)
3829             s_count++;
3830 
3831           if (nbr->GetAtomicNum() == OBElements::Phosphorus || nbr->GetAtomicNum() == OBElements::Sulfur) {
3832             FOR_NBORS_OF_ATOM(nbr2, &*nbr)
3833               if ((nbr2->GetAtomicNum() == OBElements::Sulfur || nbr2->GetAtomicNum() == OBElements::Oxygen) && (nbr2->GetExplicitDegree() == 1) && (atom->GetIdx() != nbr2->GetIdx()))
3834                 atom->SetPartialCharge(-0.5);
3835           } else
3836             atom->SetPartialCharge(-1.0);
3837 
3838           if (nbr->GetAtomicNum() == OBElements::Carbon)
3839             FOR_NBORS_OF_ATOM(nbr2, &*nbr)
3840               if (nbr2->GetAtomicNum() == OBElements::Sulfur && (nbr2->GetExplicitDegree() == 1) && (atom->GetIdx() != nbr2->GetIdx()))
3841                 atom->SetPartialCharge(-0.5); // SSMO
3842 
3843           if (s_count >= 2)
3844             atom->SetPartialCharge(-0.5); // SSMO
3845         }
3846       } else if (type == 76) {
3847        	vector<OBRing*> vr;
3848         vr = _mol.GetSSSR();
3849         vector<OBRing*>::iterator ri;
3850         vector<int>::iterator rj;
3851         int n_count;
3852 
3853         for (ri = vr.begin();ri != vr.end();++ri) { // for each ring
3854           n_count = 0;
3855 
3856           if ((*ri)->IsAromatic() && (*ri)->IsMember(&*atom) && ((*ri)->Size() == 5)) {
3857             for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) // for each ring atom
3858               if (_mol.GetAtom(*rj)->GetAtomicNum() == OBElements::Nitrogen)
3859                 n_count++;
3860 
3861             if (n_count > 1)
3862               atom->SetPartialCharge(-1.0 / n_count);
3863           }
3864         }
3865       } else if (type == 81) {
3866         atom->SetPartialCharge(1.0);
3867 
3868         vector<OBRing*> vr;
3869         vr = _mol.GetSSSR();
3870         vector<OBRing*>::iterator ri;
3871         vector<int>::iterator rj;
3872         for (ri = vr.begin();ri != vr.end();++ri) // for each ring
3873           if ((*ri)->IsAromatic() && (*ri)->IsMember(&*atom) && ((*ri)->Size() == 5)) {
3874             int n_count = 0;
3875             for(rj = (*ri)->_path.begin();rj != (*ri)->_path.end();++rj) // for each ring atom
3876               if (_mol.GetAtom(*rj)->GetAtomicNum() == OBElements::Nitrogen && (_mol.GetAtom(*rj)->GetExplicitDegree() == 3))
3877                 n_count++;
3878 
3879             if (n_count) // coverity defensive testing
3880               atom->SetPartialCharge(1.0 / n_count); // NIM+
3881 
3882             FOR_NBORS_OF_ATOM(nbr, &*atom)
3883               FOR_NBORS_OF_ATOM(nbr2, &*nbr)
3884               if (atoi(nbr2->GetType()) == 56)
3885                 atom->SetPartialCharge(1.0 / 3.0);
3886 
3887             FOR_NBORS_OF_ATOM(nbr, &*atom)
3888               FOR_NBORS_OF_ATOM(nbr2, &*nbr)
3889               if (atoi(nbr2->GetType()) == 55)
3890                 atom->SetPartialCharge(1.0 / (1.0 + n_count));
3891           }
3892       }
3893 
3894     }
3895 
3896     PrintFormalCharges();
3897 
3898     return true;
3899   }
3900 
SetPartialCharges()3901   bool OBForceFieldMMFF94::SetPartialCharges()
3902   {
3903     vector<double> charges(_mol.NumAtoms()+1, 0);
3904     double M, Wab, factor, q0a, q0b, Pa, Pb;
3905 
3906     FOR_ATOMS_OF_MOL (atom, _mol) {
3907       int type = atoi(atom->GetType());
3908 
3909       switch (type) {
3910       case 32:
3911       case 35:
3912       case 72:
3913         factor = 0.5;
3914         break;
3915       case 62:
3916       case 76:
3917         factor = 0.25;
3918         break;
3919       default:
3920         factor = 0.0;
3921         break;
3922       }
3923 
3924       M = GetCrd(type);
3925       q0a = atom->GetPartialCharge();
3926 
3927       // charge sharing
3928       if (!factor)
3929         FOR_NBORS_OF_ATOM (nbr, &*atom)
3930           if (nbr->GetPartialCharge() < 0.0)
3931             q0a += nbr->GetPartialCharge() / (2.0 * (double)(nbr->GetExplicitDegree()));
3932 
3933       // needed for SEYWUO, positive charge sharing?
3934       if (type == 62)
3935         FOR_NBORS_OF_ATOM (nbr, &*atom)
3936           if (nbr->GetPartialCharge() > 0.0)
3937             q0a -= nbr->GetPartialCharge() / 2.0;
3938 
3939       q0b = 0.0;
3940       Wab = 0.0;
3941       Pa = Pb = 0.0;
3942       FOR_NBORS_OF_ATOM (nbr, &*atom) {
3943         int nbr_type = atoi(nbr->GetType());
3944 
3945         q0b += nbr->GetPartialCharge();
3946 
3947         bool bci_found = false;
3948         for (unsigned int idx=0; idx < _ffchgparams.size(); idx++)
3949           if (GetBondType(&*atom, &*nbr) == _ffchgparams[idx]._ipar[0]) {
3950             if ((type == _ffchgparams[idx].a) && (nbr_type == _ffchgparams[idx].b)) {
3951               Wab += -_ffchgparams[idx]._dpar[0];
3952               bci_found = true;
3953             } else if  ((type == _ffchgparams[idx].b) && (nbr_type == _ffchgparams[idx].a)) {
3954               Wab += _ffchgparams[idx]._dpar[0];
3955               bci_found = true;
3956             }
3957 	  }
3958 
3959         if (!bci_found) {
3960           for (unsigned int idx=0; idx < _ffpbciparams.size(); idx++) {
3961             if (type == _ffpbciparams[idx].a)
3962               Pa = _ffpbciparams[idx]._dpar[0];
3963             if (nbr_type == _ffpbciparams[idx].a)
3964               Pb = _ffpbciparams[idx]._dpar[0];
3965           }
3966           Wab += Pa - Pb;
3967         }
3968       }
3969       if (factor)
3970         charges[atom->GetIdx()] = (1.0 - M * factor) * q0a + factor * q0b + Wab;
3971       else
3972         charges[atom->GetIdx()] = q0a + Wab;
3973     }
3974 
3975     FOR_ATOMS_OF_MOL (atom, _mol)
3976       atom->SetPartialCharge(charges[atom->GetIdx()]);
3977 
3978     PrintPartialCharges();
3979 
3980     return true;
3981   }
3982 
3983   ////////////////////////////////////////////////////////////////////////////////
3984   ////////////////////////////////////////////////////////////////////////////////
3985   //
3986   //  Validation functions
3987   //
3988   ////////////////////////////////////////////////////////////////////////////////
3989   ////////////////////////////////////////////////////////////////////////////////
3990 
3991   // used to validate the implementation
Validate()3992   bool OBForceFieldMMFF94::Validate ()
3993   {
3994     OBConversion conv;
3995     OBFormat *format_in = conv.FindFormat("mol2");
3996     vector<string> vs;
3997     vector<int> types;
3998     vector<double> fcharges, pcharges;
3999     vector<double> bond_lengths;
4000     char buffer[150];
4001     bool molfound, atomfound, bondfound, fchgfound, pchgfound;
4002     double etot, ebond, eangle, eoop, estbn, etor, evdw, eeq;
4003     double termcount; //1=bond, 2=angle, 3=strbnd, 4=torsion, 5=oop
4004     int n = 0;
4005 
4006     if (!format_in || !conv.SetInFormat(format_in)) {
4007       obErrorLog.ThrowError(__FUNCTION__, "Could not set mol2 input format", obError);
4008       return false;
4009     }
4010 
4011     ifstream ifs, ifs2;
4012     ofstream ofs;
4013 
4014     ifs.open("MMFF94_dative.mol2");
4015     if (!ifs) {
4016       obErrorLog.ThrowError(__FUNCTION__, "Could not open ./MMFF94_dative.mol2", obError);
4017       return false;
4018     }
4019 
4020     ifs2.open("MMFF94_opti.log");
4021     if (!ifs2) {
4022       obErrorLog.ThrowError(__FUNCTION__, "Coulg not open ./MMFF_opti.log", obError);
4023       return false;
4024     }
4025 
4026     ofs.open("MMFF94_openbabel.log");
4027     if (!ofs) {
4028       obErrorLog.ThrowError(__FUNCTION__, "Coulg not open ./MMFF_openbabel.log", obError);
4029       return false;
4030     }
4031 
4032     if (!_init) {
4033       ParseParamFile();
4034       _init = true;
4035     }
4036 
4037 
4038     SetLogFile(&ofs);
4039     SetLogLevel(OBFF_LOGLVL_HIGH);
4040 
4041     for (unsigned int c=1;; c++) {
4042       _mol.Clear();
4043       types.clear();
4044       fcharges.clear();
4045       pcharges.clear();
4046       bond_lengths.clear();
4047 
4048       if (!conv.Read(&_mol, &ifs))
4049         break;
4050       if (_mol.Empty())
4051         break;
4052 
4053       _ncoords = _mol.NumAtoms() * 3;
4054       _gradientPtr = new double[_ncoords];
4055 
4056       SetTypes();
4057 
4058       if ((c == 98) || (c == 692)) // CUDPAS & VUWXUG
4059         continue;
4060 
4061       termcount = 0;
4062       molfound = false;
4063       atomfound = false;
4064       bondfound = false;
4065       fchgfound = false;
4066       pchgfound = false;
4067 
4068       while (ifs2.getline(buffer, 150)) {
4069         tokenize(vs, buffer);
4070         if (vs.size() == 0) {
4071           bondfound = false;
4072           continue;
4073         }
4074 
4075         string str(buffer);
4076         if (string::npos != str.find(_mol.GetTitle(),0))
4077           molfound = true;
4078 
4079         if (atomfound) {
4080           if (n) {
4081             types.push_back(atoi(vs[2].c_str()));
4082             types.push_back(atoi(vs[5].c_str()));
4083             types.push_back(atoi(vs[8].c_str()));
4084             types.push_back(atoi(vs[11].c_str()));
4085           } else {
4086             if (vs.size() > 2)
4087               types.push_back(atoi(vs[2].c_str()));
4088             if (vs.size() > 5)
4089               types.push_back(atoi(vs[5].c_str()));
4090             if (vs.size() > 8)
4091               types.push_back(atoi(vs[8].c_str()));
4092 
4093             atomfound = false;
4094           }
4095           n--;
4096         }
4097 
4098         if (fchgfound) {
4099           if (n) {
4100             fcharges.push_back(atof(vs[2].c_str()));
4101             fcharges.push_back(atof(vs[5].c_str()));
4102             fcharges.push_back(atof(vs[8].c_str()));
4103             fcharges.push_back(atof(vs[11].c_str()));
4104           } else {
4105             if (vs.size() > 2)
4106               fcharges.push_back(atof(vs[2].c_str()));
4107             if (vs.size() > 5)
4108               fcharges.push_back(atof(vs[5].c_str()));
4109             if (vs.size() > 8)
4110               fcharges.push_back(atof(vs[8].c_str()));
4111 
4112             fchgfound = false;
4113           }
4114           n--;
4115         }
4116 
4117         if (pchgfound) {
4118           if (n) {
4119             pcharges.push_back(atof(vs[2].c_str()));
4120             pcharges.push_back(atof(vs[5].c_str()));
4121             pcharges.push_back(atof(vs[8].c_str()));
4122             pcharges.push_back(atof(vs[11].c_str()));
4123           } else {
4124             if (vs.size() > 2)
4125               pcharges.push_back(atof(vs[2].c_str()));
4126             if (vs.size() > 5)
4127               pcharges.push_back(atof(vs[5].c_str()));
4128             if (vs.size() > 8)
4129               pcharges.push_back(atof(vs[8].c_str()));
4130 
4131             pchgfound = false;
4132           }
4133           n--;
4134         }
4135 
4136         if (molfound && EQn(buffer, " ATOM NAME  TYPE", 16)) {
4137           atomfound = true;
4138           n = _mol.NumAtoms() / 4;
4139         }
4140         if (molfound && EQn(buffer, "   ATOM   FCHARGE", 17)) {
4141           fchgfound = true;
4142           n = _mol.NumAtoms() / 4;
4143         }
4144         if (molfound && EQn(buffer, "   ATOM    CHARGE", 17)) {
4145           pchgfound = true;
4146           n = _mol.NumAtoms() / 4;
4147         }
4148 
4149         if (bondfound)
4150           bond_lengths.push_back(atof(vs[7].c_str()));
4151 
4152         if (molfound) {
4153           if (EQn(buffer, " Total ENERGY", 13))
4154             etot = atof(vs[3].c_str());
4155           if (EQn(buffer, " Bond Stretching", 16))
4156             ebond = atof(vs[2].c_str());
4157           if (EQn(buffer, " Angle Bending", 14))
4158             eangle = atof(vs[2].c_str());
4159           if (EQn(buffer, " Out-of-Plane Bending", 21))
4160             eoop = atof(vs[2].c_str());
4161           if (EQn(buffer, " Stretch-Bend", 13))
4162             estbn = atof(vs[1].c_str());
4163           if (EQn(buffer, "     Total Torsion", 18))
4164             etor = atof(vs[2].c_str());
4165           if (EQn(buffer, "     Net vdW", 12))
4166             evdw = atof(vs[2].c_str());
4167           if (EQn(buffer, " Electrostatic", 14))
4168             eeq = atof(vs[1].c_str());
4169           if (EQn(buffer, " ---------------------", 22) && (termcount == 0)) {
4170             termcount++;
4171             bondfound = true;
4172           }
4173           if (EQn(buffer, " OPTIMOL>  # read next", 22))
4174             break;
4175         }
4176 
4177 
4178 
4179       } // while (getline)
4180 
4181       vector<int>::iterator i;
4182       vector<double>::iterator di;
4183       unsigned int ni;
4184       bool failed;
4185 
4186       cout << "--------------------------------------------------------------------------------" << endl;
4187       cout << "                                                                                " << endl;
4188       cout << "  VALIDATE MOLECULE " << c << ": " << _mol.GetTitle() << endl;
4189       cout << "                                                                                " << endl;
4190       cout << "IDX  HYB  AROM  OB_TYPE  LOG_TYPE       RESULT                                  " << endl;
4191       cout << "----------------------------------------------                                  " << endl;
4192 
4193       ni = 1;
4194       failed = false;
4195       for (i = types.begin(); i != types.end();++i) {
4196         if (ni > _mol.NumAtoms())
4197           continue;
4198 
4199         if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) ||
4200              (atoi(_mol.GetAtom(ni)->GetType()) == 97)
4201              ) continue;
4202 
4203         if (atoi(_mol.GetAtom(ni)->GetType()) == (*i))
4204           snprintf(_logbuf, BUFF_SIZE, "%2d   %3d  %4d    %3d      %3d          PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetHyb(),
4205                   _mol.GetAtom(ni)->IsAromatic(), atoi(_mol.GetAtom(ni)->GetType()), *i);
4206         else {
4207           snprintf(_logbuf, BUFF_SIZE, "%2d   %3d  %4d    %3d      %3d      XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetHyb(),
4208                   _mol.GetAtom(ni)->IsAromatic(), atoi(_mol.GetAtom(ni)->GetType()), *i);
4209           failed = true;
4210         }
4211 
4212         cout << _logbuf << endl;
4213 
4214         ni++;
4215       }
4216 
4217       if (failed) {
4218         cout << "Could not successfully assign atom types" << endl;
4219         return false;
4220         //continue;
4221       }
4222 
4223       SetFormalCharges();
4224       cout << endl;
4225       cout << "IDX  OB_FCARGE  LOG_FCHARGE       RESULT" << endl;
4226       cout << "----------------------------------------" << endl;
4227 
4228       ni = 1;
4229       for (di = fcharges.begin(); di != fcharges.end(); ++di) {
4230         if (ni > _mol.NumAtoms())
4231           continue;
4232 
4233         if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) ||
4234              (atoi(_mol.GetAtom(ni)->GetType()) == 97)
4235              ) continue;
4236 
4237         if (fabs((*di) - _mol.GetAtom(ni)->GetPartialCharge()) <= 0.001)
4238           snprintf(_logbuf, BUFF_SIZE, "%2d   %7.4f     %7.4f          PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di);
4239         else {
4240           snprintf(_logbuf, BUFF_SIZE, "%2d   %7.4f     %7.4f      XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di);
4241           failed = true;
4242         }
4243 
4244         cout << _logbuf << endl;
4245 
4246         ni++;
4247       }
4248 
4249       if (failed) {
4250         cout << "Could not successfully assign formal charges" << endl;
4251         //return false;
4252         continue;
4253       }
4254 
4255       SetPartialCharges();
4256       cout << endl;
4257       cout << "IDX  OB_PCARGE  LOG_PCHARGE       RESULT" << endl;
4258       cout << "----------------------------------------" << endl;
4259 
4260       ni = 1;
4261       for (di = pcharges.begin(); di != pcharges.end(); ++di) {
4262         if (ni > _mol.NumAtoms())
4263           continue;
4264 
4265         if ( (atoi(_mol.GetAtom(ni)->GetType()) == 87) ||
4266              (atoi(_mol.GetAtom(ni)->GetType()) == 97)
4267              ) continue;
4268 
4269         if (fabs((*di) - _mol.GetAtom(ni)->GetPartialCharge()) <= 0.001)
4270           snprintf(_logbuf, BUFF_SIZE, "%2d   %7.4f     %7.4f          PASSED", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di);
4271         else {
4272           snprintf(_logbuf, BUFF_SIZE, "%2d   %7.4f     %7.4f      XXX FAILED XXX", _mol.GetAtom(ni)->GetIdx(), _mol.GetAtom(ni)->GetPartialCharge(), *di);
4273           failed = true;
4274         }
4275 
4276         cout << _logbuf << endl;
4277 
4278         ni++;
4279       }
4280 
4281       if (failed) {
4282         cout << "Could not successfully assign partial charges" << endl;
4283         //return false;
4284         continue;
4285       }
4286 
4287 
4288 
4289       if (!SetupCalculations()) {
4290         cout << "Could not setup calculations (missing parameters...)" << endl;
4291         return false;
4292         //continue;
4293       }
4294 
4295       double delta;
4296       cout << endl;
4297       cout << "TERM                     OB ENERGY     LOG ENERGY         DELTA" << endl;
4298       cout << "---------------------------------------------------------------" << endl;
4299 
4300       delta = (E_Bond() - ebond);
4301       snprintf(_logbuf, BUFF_SIZE, "Bond Stretching        %11.5f    %11.5f   %11.5f", E_Bond(), ebond, delta);
4302       cout << _logbuf << endl;
4303 
4304       delta = (E_Angle() - eangle);
4305       snprintf(_logbuf, BUFF_SIZE, "Angle Bending          %11.5f    %11.5f   %11.5f", E_Angle(), eangle, delta);
4306       cout << _logbuf << endl;
4307 
4308       delta = (E_StrBnd() - estbn);
4309       snprintf(_logbuf, BUFF_SIZE, "Stretch-Bending        %11.5f    %11.5f   %11.5f", E_StrBnd(), estbn, delta);
4310       cout << _logbuf << endl;
4311 
4312       delta = (E_OOP() - eoop);
4313       snprintf(_logbuf, BUFF_SIZE, "Out-Of-Plane Bending   %11.5f    %11.5f   %11.5f", E_OOP(), eoop, delta);
4314       cout << _logbuf << endl;
4315 
4316       delta = (E_Torsion() - etor);
4317       snprintf(_logbuf, BUFF_SIZE, "Torsional              %11.5f    %11.5f   %11.5f", E_Torsion(), etor, delta);
4318       cout << _logbuf << endl;
4319 
4320       delta = (E_VDW() - evdw);
4321       snprintf(_logbuf, BUFF_SIZE, "Van der Waals          %11.5f    %11.5f   %11.5f", E_VDW(), evdw, delta);
4322       cout << _logbuf << endl;
4323 
4324       delta = (E_Electrostatic() - eeq);
4325       snprintf(_logbuf, BUFF_SIZE, "Electrostatic          %11.5f    %11.5f   %11.5f", E_Electrostatic(), eeq, delta);
4326       cout << _logbuf << endl;
4327 
4328       cout << endl;
4329       delta = (Energy() - etot);
4330       snprintf(_logbuf, BUFF_SIZE, "Total ENERGY           %11.5f    %11.5f   %11.5f", Energy(), etot, delta);
4331       cout << _logbuf << endl;
4332 
4333     } // for (unsigned int c;; c++ )
4334 
4335     if (ifs)
4336       ifs.close();
4337     if (ifs2)
4338       ifs2.close();
4339 
4340     return true;
4341   }
4342 
ValidateGradients()4343   bool OBForceFieldMMFF94::ValidateGradients ()
4344   {
4345     vector3 numgrad, anagrad, err;
4346     int coordIdx;
4347 
4348     bool passed = true; // set to false if any component fails
4349 
4350     cout << "----------------------------------------------------------------------------------------" << endl;
4351     cout << "                                                                                        " << endl;
4352     cout << "  VALIDATE GRADIENTS : " << _mol.GetTitle() << endl;
4353     cout << "                                                                                        " << endl;
4354     cout << "                                                                                        " << endl;
4355     cout << "ATOM IDX      NUMERICAL GRADIENT           ANALYTICAL GRADIENT        REL. ERROR (%)   " << endl;
4356     cout << "----------------------------------------------------------------------------------------" << endl;
4357     //     "XX       (000.000, 000.000, 000.000)  (000.000, 000.000, 000.000)  (00.00, 00.00, 00.00)"
4358 
4359     FOR_ATOMS_OF_MOL (a, _mol) {
4360       coordIdx = (a->GetIdx() - 1) * 3;
4361 
4362       // OBFF_ENERGY
4363       numgrad = NumericalDerivative(&*a, OBFF_ENERGY);
4364       Energy(); // compute
4365       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4366       err = ValidateGradientError(numgrad, anagrad);
4367 
4368       snprintf(_logbuf, BUFF_SIZE, "%2d       (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", a->GetIdx(), numgrad.x(), numgrad.y(), numgrad.z(),
4369               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4370       OBFFLog(_logbuf);
4371 
4372       // OBFF_EBOND
4373       numgrad = NumericalDerivative(&*a, OBFF_EBOND);
4374       ClearGradients();
4375       E_Bond(); // compute
4376       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4377       err = ValidateGradientError(numgrad, anagrad);
4378 
4379       snprintf(_logbuf, BUFF_SIZE, "    bond    (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4380               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4381       OBFFLog(_logbuf);
4382       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4383         passed = false;
4384 
4385       // OBFF_EANGLE
4386       numgrad = NumericalDerivative(&*a, OBFF_EANGLE);
4387       ClearGradients();
4388       E_Angle(); // compute
4389       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4390       err = ValidateGradientError(numgrad, anagrad);
4391 
4392       snprintf(_logbuf, BUFF_SIZE, "    angle   (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4393               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4394       OBFFLog(_logbuf);
4395       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4396         passed = false;
4397 
4398       // OBFF_ESTRBND
4399       numgrad = NumericalDerivative(&*a, OBFF_ESTRBND);
4400       ClearGradients();
4401       E_StrBnd(); // compute
4402       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4403       err = ValidateGradientError(numgrad, anagrad);
4404 
4405       snprintf(_logbuf, BUFF_SIZE, "    strbnd  (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4406               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4407       OBFFLog(_logbuf);
4408       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4409         passed = false;
4410 
4411       // OBFF_ETORSION
4412       numgrad = NumericalDerivative(&*a, OBFF_ETORSION);
4413       ClearGradients();
4414       E_Torsion(); // compute
4415       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4416       err = ValidateGradientError(numgrad, anagrad);
4417 
4418       snprintf(_logbuf, BUFF_SIZE, "    torsion (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4419               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4420       OBFFLog(_logbuf);
4421       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4422         passed = false;
4423 
4424       // OBFF_EOOP
4425       numgrad = NumericalDerivative(&*a, OBFF_EOOP);
4426       ClearGradients();
4427       E_OOP(); // compute
4428       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4429       err = ValidateGradientError(numgrad, anagrad);
4430 
4431       snprintf(_logbuf, BUFF_SIZE, "    oop     (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4432               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4433       OBFFLog(_logbuf);
4434       // disable OOP gradient validation for now -- some small errors, but nothing major
4435       //      if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4436       //        passed = false;
4437 
4438       // OBFF_EVDW
4439       numgrad = NumericalDerivative(&*a, OBFF_EVDW);
4440       ClearGradients();
4441       E_VDW(); // compute
4442       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4443       err = ValidateGradientError(numgrad, anagrad);
4444 
4445       snprintf(_logbuf, BUFF_SIZE, "    vdw     (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4446               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4447       OBFFLog(_logbuf);
4448       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4449         passed = false;
4450 
4451       // OBFF_EELECTROSTATIC
4452       numgrad = NumericalDerivative(&*a, OBFF_EELECTROSTATIC);
4453       ClearGradients();
4454       E_Electrostatic(); // compute
4455       anagrad.Set(_gradientPtr[coordIdx], _gradientPtr[coordIdx+1], _gradientPtr[coordIdx+2]);
4456       err = ValidateGradientError(numgrad, anagrad);
4457 
4458       snprintf(_logbuf, BUFF_SIZE, "    electro (%7.3f, %7.3f, %7.3f)  (%7.3f, %7.3f, %7.3f)  (%5.2f, %5.2f, %5.2f)\n", numgrad.x(), numgrad.y(), numgrad.z(),
4459               anagrad.x(), anagrad.y(), anagrad.z(), err.x(), err.y(), err.z());
4460       OBFFLog(_logbuf);
4461       if (err.x() > 5.0 || err.y() > 5.0 || err.z() > 5.0)
4462         passed = false;
4463     }
4464 
4465     return passed; // did we pass every single component?
4466   }
4467 
4468   ////////////////////////////////////////////////////////////////////////////////
4469   ////////////////////////////////////////////////////////////////////////////////
4470   //
4471   //  Calculate bond type, angle type, stretch-bend type, torsion type
4472   //
4473   ////////////////////////////////////////////////////////////////////////////////
4474   ////////////////////////////////////////////////////////////////////////////////
4475 
4476   //
4477   // MMFF part V - page 620
4478   //
4479   // BTij is 1 when:
4480   // a) single bond between atoms i and j, both i and j are not aromatic and both types have sbmb set in mmffprop.par, or
4481   // b) bewtween two aromatic atoms, but the bond is not aromatic (e.g. connecting bond in biphenyl)
4482   //
GetBondType(OBAtom * a,OBAtom * b)4483   int OBForceFieldMMFF94::GetBondType(OBAtom* a, OBAtom* b)
4484   {
4485     OBBond *bond = _mol.GetBond(a, b);
4486     if (bond->GetBondOrder() != 1 || bond->IsAromatic())
4487       return 0;
4488 
4489     if (HasAromSet(atoi(a->GetType())) && HasAromSet(atoi(b->GetType())))
4490       return 1;
4491 
4492     if (HasSbmbSet(atoi(a->GetType())) && HasSbmbSet(atoi(b->GetType())))
4493       return 1;
4494 
4495     return 0;
4496   }
4497 
GetAngleType(OBAtom * a,OBAtom * b,OBAtom * c)4498   int OBForceFieldMMFF94::GetAngleType(OBAtom* a, OBAtom* b, OBAtom *c)
4499   {
4500     int sumbondtypes;
4501 
4502     sumbondtypes = GetBondType(a,b) + GetBondType(b, c);
4503 
4504     if (a->IsInRingSize(3) && b->IsInRingSize(3) && c->IsInRingSize(3) && IsInSameRing(a, c))
4505       switch (sumbondtypes) {
4506       case 0:
4507         return 3;
4508       case 1:
4509         return 5;
4510       case 2:
4511         return 6;
4512       }
4513 
4514     if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && IsInSameRing(a, c))
4515       switch (sumbondtypes) {
4516       case 0:
4517         return 4;
4518       case 1:
4519         return 7;
4520       case 2:
4521         return 8;
4522       }
4523 
4524     return sumbondtypes;
4525   }
4526 
GetStrBndType(OBAtom * a,OBAtom * b,OBAtom * c)4527   int OBForceFieldMMFF94::GetStrBndType(OBAtom* a, OBAtom* b, OBAtom *c)
4528   {
4529     int btab, btbc, atabc;
4530     bool inverse;
4531 
4532     btab = GetBondType(a, b);
4533     btbc = GetBondType(b, c);
4534     atabc = GetAngleType(a, b, c);
4535 
4536     if (atoi(a->GetType()) <= atoi(c->GetType()))
4537       inverse = false;
4538     else
4539       inverse = true;
4540 
4541     switch (atabc) {
4542     case 0:
4543       return 0;
4544 
4545     case 1:
4546       if (btab) {
4547         if (!inverse) {
4548           return 1;
4549         } else {
4550           return 2;
4551         }
4552       }
4553       if (btbc) {
4554         if (!inverse) {
4555           return 2;
4556         } else {
4557           return 1;
4558         }
4559       }
4560 
4561     case 2:
4562       return 3;
4563 
4564     case 3:
4565       return 5;
4566 
4567     case 4:
4568       return 4;
4569 
4570     case 5:
4571       if (btab) {
4572         if (!inverse) {
4573           return 6;
4574         } else {
4575           return 7;
4576         }
4577       }
4578       if (btbc) {
4579         if (!inverse) {
4580           return 7;
4581         } else {
4582           return 6;
4583         }
4584       }
4585 
4586     case 6:
4587       return 8;
4588 
4589     case 7:
4590       if (btab) {
4591         if (!inverse) {
4592           return 9;
4593         } else {
4594           return 10;
4595         }
4596       }
4597       if (btbc) {
4598         if (!inverse) {
4599           return 10;
4600         } else {
4601           return 9;
4602         }
4603       }
4604 
4605     case 8:
4606       return 11;
4607     }
4608 
4609     return 0;
4610   }
4611 
4612   //
4613   // MMFF part IV - page 609
4614   //
4615   // TTijkl = 1 when BTjk = 1
4616   // TTijkl = 2 when BTjk = 0 but BTij and/or BTkl = 1
4617   // TTijkl = 4 when i, j, k and l are all members of the same four-membered ring
4618   // TTijkl = 5 when i, j, k and l are members of a five-membered ring and at least one is a sp3-hybridized carbon (MMFF atom type 1)
4619   //
GetTorsionType(OBAtom * a,OBAtom * b,OBAtom * c,OBAtom * d)4620   int OBForceFieldMMFF94::GetTorsionType(OBAtom* a, OBAtom* b, OBAtom *c, OBAtom *d)
4621   {
4622     int btab, btbc, btcd;
4623 
4624     btab = GetBondType(a, b);
4625     btbc = GetBondType(b, c);
4626     btcd = GetBondType(c, d);
4627 
4628     if (btbc == 1)
4629       return 1;
4630 
4631     if (a->IsInRingSize(4) && b->IsInRingSize(4) && c->IsInRingSize(4) && d->IsInRingSize(4))
4632       if (IsInSameRing(a,b) && IsInSameRing(b,c) && IsInSameRing(c,d))
4633         return 4;
4634 
4635     OBBond *bond = _mol.GetBond(b, c);
4636     if (bond->GetBondOrder() == 1 && !bond->IsAromatic()) {
4637       if (btab || btcd)
4638         return 2;
4639       /*
4640         unsigned int order1 = GetCXT(0, atoi(d->GetType()), atoi(c->GetType()), atoi(b->GetType()), atoi(a->GetType()));
4641         unsigned int order2 = GetCXT(0, atoi(a->GetType()), atoi(b->GetType()), atoi(c->GetType()), atoi(d->GetType()));
4642 
4643         cout << "GetTorsionType(" << a->GetType() << ", " << b->GetType() << ", " << c->GetType() << ", " << d->GetType() << ")" << endl;
4644         cout << "    order1 = " << order1 << endl;
4645         cout << "    order2 = " << order2 << endl;
4646         cout << "    btab = " << btab << endl;
4647         cout << "    btbc = " << btbc << endl;
4648         cout << "    btcd = " << btcd << endl;
4649       */
4650     }
4651 
4652     if (a->IsInRingSize(5) && b->IsInRingSize(5) && c->IsInRingSize(5) && d->IsInRingSize(5)) {
4653       vector<OBRing*> vr;
4654       vr = _mol.GetSSSR();
4655 
4656       if( !((atoi(a->GetType()) == 1) || (atoi(b->GetType()) == 1) || (atoi(c->GetType()) == 1) || (atoi(d->GetType()) == 1)) )
4657         return 0;
4658 
4659       vector<OBRing*>::iterator ri;
4660       vector<int>::iterator rj;
4661       for (ri = vr.begin();ri != vr.end();++ri) { // for each ring
4662         if ((*ri)->IsAromatic())
4663           continue;
4664 
4665         if ((*ri)->Size() != 5)
4666           continue;
4667 
4668         if (!(*ri)->IsMember(a) || !(*ri)->IsMember(b) || !(*ri)->IsMember(c) || !(*ri)->IsMember(d))
4669           continue;
4670 
4671         return 5;
4672       }
4673     }
4674 
4675 
4676     return 0;
4677   }
4678 
4679   // CXB = MC * (I * MA + J) + BTij
GetCXB(int type,int a,int b)4680   unsigned int OBForceFieldMMFF94::GetCXB(int type, int a, int b)
4681   {
4682     unsigned int cxb;
4683     cxb = 2 * (a * 136 + b) + type;
4684     return cxb;
4685   }
4686 
4687   // CXA = MC * (J * MA^2 + I * MA + K) + ATijk
GetCXA(int type,int a,int b,int c)4688   unsigned int OBForceFieldMMFF94::GetCXA(int type, int a, int b, int c)
4689   {
4690     unsigned int cxa;
4691     cxa = 9 * (b * 18496 + a * 136 + c) + type;
4692     return cxa;
4693   }
4694 
4695   // CXS = MC * (J * MA^2 + I * MA + K) + STijk
GetCXS(int type,int a,int b,int c)4696   unsigned int OBForceFieldMMFF94::GetCXS(int type, int a, int b, int c)
4697   {
4698     unsigned int cxs;
4699     cxs = 12 * (b * 18496 + a * 136 + c) + type;
4700     return cxs;
4701   }
4702 
4703   // CXO = J * MA^3 + I * MA^2 + K * MA + L
GetCXO(int a,int b,int c,int d)4704   unsigned int OBForceFieldMMFF94::GetCXO(int a, int b, int c, int d)
4705   {
4706     unsigned int cxo;
4707     cxo = b * 2515456 + a * 18496 + c * 136 + d;
4708     return cxo;
4709   }
4710 
4711   // CXT = MC * (J * MA^3 + K * MA^2 + I * MA + L) + TTijkl
GetCXT(int type,int a,int b,int c,int d)4712   unsigned int OBForceFieldMMFF94::GetCXT(int type, int a, int b, int c, int d)
4713   {
4714     unsigned int cxt;
4715     cxt = 6 * (b * 2515456 + c * 18496 + a * 136 + d) + type;
4716     return cxt;
4717   }
4718 
4719   // CXQ = MC * (I * MA + J) + BTij
GetCXQ(int type,int a,int b)4720   unsigned int OBForceFieldMMFF94::GetCXQ(int type, int a, int b)
4721   {
4722     unsigned int cxq;
4723     cxq = 2 * (a * 136 + b) + type;
4724     return cxq;
4725   }
4726 
4727   ////////////////////////////////////////////////////////////////////////////////
4728   ////////////////////////////////////////////////////////////////////////////////
4729   //
4730   //  Various tables & misc. functions
4731   //
4732   ////////////////////////////////////////////////////////////////////////////////
4733   ////////////////////////////////////////////////////////////////////////////////
4734 
4735   // MMFF part V - TABLE I
HasLinSet(int atomtype)4736   bool OBForceFieldMMFF94::HasLinSet(int atomtype)
4737   {
4738     return _ffpropLin.BitIsSet(atomtype);
4739   }
4740 
4741   // MMFF part V - TABLE I
HasPilpSet(int atomtype)4742   bool OBForceFieldMMFF94::HasPilpSet(int atomtype)
4743   {
4744     return _ffpropPilp.BitIsSet(atomtype);
4745   }
4746 
4747   // MMFF part V - TABLE I
HasAromSet(int atomtype)4748   bool OBForceFieldMMFF94::HasAromSet(int atomtype)
4749   {
4750     return _ffpropArom.BitIsSet(atomtype);
4751   }
4752 
4753   // MMFF part V - TABLE I
HasSbmbSet(int atomtype)4754   bool OBForceFieldMMFF94::HasSbmbSet(int atomtype)
4755   {
4756     return _ffpropSbmb.BitIsSet(atomtype);
4757   }
4758 
4759   // MMFF part V - TABLE I
GetCrd(int atomtype)4760   int OBForceFieldMMFF94::GetCrd(int atomtype)
4761   {
4762     OBFFParameter *par;
4763 
4764     par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par
4765     if (par)
4766       return par->_ipar[1];
4767 
4768     return 0;
4769   }
4770 
4771   // MMFF part V - TABLE I
GetVal(int atomtype)4772   int OBForceFieldMMFF94::GetVal(int atomtype)
4773   {
4774     OBFFParameter *par;
4775 
4776     par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par
4777     if (par)
4778       return par->_ipar[2];
4779 
4780     return 0;
4781   }
4782 
4783   // MMFF part V - TABLE I
GetMltb(int atomtype)4784   int OBForceFieldMMFF94::GetMltb(int atomtype)
4785   {
4786     OBFFParameter *par;
4787 
4788     par = GetParameter1Atom(atomtype, _ffpropparams); // from mmffprop.par
4789     if (par)
4790       return par->_ipar[4];
4791 
4792     return 0;
4793   }
4794 
4795   // MMFF part I - TABLE IV
EqLvl2(int type)4796   int OBForceFieldMMFF94::EqLvl2(int type)
4797   {
4798     for (unsigned int idx=0; idx < _ffdefparams.size(); idx++)
4799       if (_ffdefparams[idx]._ipar[0] == type)
4800         return _ffdefparams[idx]._ipar[1];
4801 
4802     return type;
4803   }
4804 
4805   // MMFF part I - TABLE IV
EqLvl3(int type)4806   int OBForceFieldMMFF94::EqLvl3(int type)
4807   {
4808     for (unsigned int idx=0; idx < _ffdefparams.size(); idx++)
4809       if (_ffdefparams[idx]._ipar[0] == type)
4810         return _ffdefparams[idx]._ipar[2];
4811 
4812     return type;
4813   }
4814 
4815   // MMFF part I - TABLE IV
EqLvl4(int type)4816   int OBForceFieldMMFF94::EqLvl4(int type)
4817   {
4818     for (unsigned int idx=0; idx < _ffdefparams.size(); idx++)
4819       if (_ffdefparams[idx]._ipar[0] == type)
4820         return _ffdefparams[idx]._ipar[3];
4821 
4822     return type;
4823   }
4824 
4825   // MMFF part I - TABLE IV
EqLvl5(int type)4826   int OBForceFieldMMFF94::EqLvl5(int type)
4827   {
4828     for (unsigned int idx=0; idx < _ffdefparams.size(); idx++)
4829       if (_ffdefparams[idx]._ipar[0] == type)
4830         return _ffdefparams[idx]._ipar[4];
4831 
4832     return type;
4833   }
4834 
4835   // MMFF part V - TABLE VI
GetZParam(OBAtom * atom)4836   double OBForceFieldMMFF94::GetZParam(OBAtom* atom)
4837   {
4838     switch (atom->GetAtomicNum()) {
4839     case OBElements::Hydrogen:
4840       return 1.395;
4841     case OBElements::Carbon:
4842       return 2.494;
4843     case OBElements::Nitrogen:
4844       return 2.711;
4845     case OBElements::Oxygen:
4846       return 3.045;
4847     case OBElements::Fluorine:
4848       return 2.847;
4849     case OBElements::Silicon:
4850       return 2.350;
4851     case OBElements::Phosphorus:
4852       return 2.350;
4853     case OBElements::Sulfur:
4854       return 2.980;
4855     case OBElements::Chlorine:
4856       return 2.909;
4857     case OBElements::Bromine:
4858       return 3.017;
4859     case OBElements::Iodine:
4860       return 3.086;
4861     }
4862 
4863     return 0.0;
4864   }
4865 
4866   // MMFF part V - TABLE VI
GetCParam(OBAtom * atom)4867   double OBForceFieldMMFF94::GetCParam(OBAtom* atom)
4868   {
4869     switch (atom->GetAtomicNum()) {
4870     case OBElements::Boron:
4871       return 0.704;
4872     case OBElements::Carbon:
4873       return 1.016;
4874     case OBElements::Nitrogen:
4875       return 1.113;
4876     case OBElements::Oxygen:
4877       return 1.337;
4878     case OBElements::Silicon:
4879       return 0.811;
4880     case OBElements::Phosphorus:
4881       return 1.068;
4882     case OBElements::Sulfur:
4883       return 1.249;
4884     case OBElements::Chlorine:
4885       return 1.078;
4886     case OBElements::Arsenic:
4887       return 0.825;
4888     }
4889 
4890     return 0.0;
4891   }
4892 
4893   // MMFF part V - TABLE X
GetUParam(OBAtom * atom)4894   double OBForceFieldMMFF94::GetUParam(OBAtom* atom)
4895   {
4896     switch (atom->GetAtomicNum()) {
4897     case OBElements::Carbon:
4898       return 2.0;
4899     case OBElements::Nitrogen:
4900       return 2.0;
4901     case OBElements::Oxygen:
4902       return 2.0;
4903     case OBElements::Silicon:
4904       return 1.25;
4905     case OBElements::Phosphorus:
4906       return 1.25;
4907     case OBElements::Sulfur:
4908       return 1.25;
4909     }
4910 
4911     return 0.0;
4912   }
4913 
4914   // MMFF part V - TABLE X
GetVParam(OBAtom * atom)4915   double OBForceFieldMMFF94::GetVParam(OBAtom* atom)
4916   {
4917     switch (atom->GetAtomicNum()) {
4918     case OBElements::Carbon:
4919       return 2.12;
4920     case OBElements::Nitrogen:
4921       return 1.5;
4922     case OBElements::Oxygen:
4923       return 0.2;
4924     case OBElements::Silicon:
4925       return 1.22;
4926     case OBElements::Phosphorus:
4927       return 2.4;
4928     case OBElements::Sulfur:
4929       return 0.49;
4930     }
4931 
4932     return 0.0;
4933   }
4934 
4935   // R Blom and A Haaland, J. Mol. Struct., 128, 21-27 (1985)
GetCovalentRadius(OBAtom * a)4936   double OBForceFieldMMFF94::GetCovalentRadius(OBAtom* a) {
4937 
4938     switch (a->GetAtomicNum()) {
4939     case 1:
4940       return 0.33; // corrected value from MMFF part V
4941     case 5:
4942       return 0.81;
4943     case 6:
4944       return 0.77; // corrected value from MMFF part V
4945     case 7:
4946       return 0.73;
4947     case 8:
4948       return 0.72;
4949     case 9:
4950       return 0.74;
4951     case 13:
4952       return 1.22;
4953     case 14:
4954       return 1.15;
4955     case 15:
4956       return 1.09;
4957     case 16:
4958       return 1.03;
4959     case 17:
4960       return 1.01;
4961     case 31:
4962       return 1.19;
4963     case 32:
4964       return 1.20;
4965     case 33:
4966       return 1.20;
4967     case 34:
4968       return 1.16;
4969     case 35:
4970       return 1.15;
4971     case 44:
4972       return 1.46;
4973     case 50:
4974       return 1.40;
4975     case 51:
4976       return 1.41;
4977     case 52:
4978       return 1.35;
4979     case 53:
4980       return 1.33;
4981     case 81:
4982       return 1.51;
4983     case 82:
4984       return 1.53;
4985     case 83:
4986       return 1.55;
4987     default:
4988       return OBElements::GetCovalentRad(a->GetAtomicNum());
4989     }
4990   }
4991 
GetBondLength(OBAtom * a,OBAtom * b)4992   double OBForceFieldMMFF94::GetBondLength(OBAtom* a, OBAtom* b)
4993   {
4994     OBFFParameter *parameter;
4995     double rab;
4996 
4997     parameter = GetTypedParameter2Atom(GetBondType(a, b), atoi(a->GetType()), atoi(b->GetType()), _ffbondparams);
4998     if (parameter == nullptr)
4999       rab = GetRuleBondLength(a, b);
5000     else
5001       rab = parameter->_dpar[1];
5002 
5003     return rab;
5004   }
5005 
5006   // MMFF part V - page 625
GetRuleBondLength(OBAtom * a,OBAtom * b)5007   double OBForceFieldMMFF94::GetRuleBondLength(OBAtom* a, OBAtom* b)
5008   {
5009     double r0ab, r0a, r0b, c, Xa, Xb;
5010     int Ha, Hb, BOab;
5011     r0a = GetCovalentRadius(a);
5012     r0b = GetCovalentRadius(b);
5013     Xa = OBElements::GetAllredRochowElectroNeg(a->GetAtomicNum());
5014     Xb = OBElements::GetAllredRochowElectroNeg(b->GetAtomicNum());
5015 
5016 
5017     if (a->GetAtomicNum() == OBElements::Hydrogen)
5018       r0a = 0.33;
5019     if (b->GetAtomicNum() == OBElements::Hydrogen)
5020       r0b = 0.33;
5021 
5022     if (a->GetAtomicNum() == OBElements::Hydrogen || b->GetAtomicNum() == OBElements::Hydrogen)
5023       c = 0.050;
5024     else
5025       c = 0.085;
5026 
5027     if (GetMltb(atoi(a->GetType()) == 3))
5028       Ha = 1;
5029     else if ((GetMltb(atoi(a->GetType())) == 1) || (GetMltb(atoi(a->GetType())) == 2))
5030       Ha = 2;
5031     else
5032       Ha = 3;
5033 
5034     if (GetMltb(atoi(b->GetType()) == 3))
5035       Hb = 1;
5036     else if ((GetMltb(atoi(b->GetType())) == 1) || (GetMltb(atoi(b->GetType())) == 2))
5037       Hb = 2;
5038     else
5039       Hb = 3;
5040 
5041     BOab = a->GetBond(b)->GetBondOrder();
5042     if ((GetMltb(atoi(a->GetType())) == 1) && (GetMltb(atoi(b->GetType())) == 1))
5043       BOab = 4;
5044     if ((GetMltb(atoi(a->GetType())) == 1) && (GetMltb(atoi(b->GetType())) == 2))
5045       BOab = 5;
5046     if ((GetMltb(atoi(a->GetType())) == 2) && (GetMltb(atoi(b->GetType())) == 1))
5047       BOab = 5;
5048     if (a->GetBond(b)->IsAromatic()) {
5049       if (!HasPilpSet(atoi(a->GetType())) && !HasPilpSet(atoi(b->GetType()))) {
5050         BOab = 4;
5051       } else {
5052         BOab = 5;
5053       }
5054     }
5055 
5056     switch (BOab) {
5057     case 5:
5058       r0a -= 0.04;
5059       r0b -= 0.04;
5060       break;
5061     case 4:
5062       r0a -= 0.075;
5063       r0b -= 0.075;
5064       break;
5065     case 3:
5066       r0a -= 0.17;
5067       r0b -= 0.17;
5068       break;
5069     case 2:
5070       r0a -= 0.10;
5071       r0b -= 0.10;
5072       break;
5073     case 1:
5074       if (Ha == 1)
5075         r0a -= 0.08;
5076       if (Ha == 2)
5077         r0a -= 0.03;
5078       if (Hb == 1)
5079         r0b -= 0.08;
5080       if (Hb == 2)
5081         r0b -= 0.03;
5082     }
5083 
5084     /*
5085       cout << "Ha=" << Ha << "  Hb=" << Hb << "  BOab=" << BOab << endl;
5086       cout << "r0a=" << r0a << "  Xa=" << Xa << endl;
5087       cout << "r0b=" << r0b << "  Xb=" << Xb << endl;
5088       cout << "r0a + r0b=" << r0a +r0b << endl;
5089       cout << "c=" << c << "  |Xa-Xb|=" << fabs(Xa-Xb) << "  |Xa-Xb|^1.4=" << pow(fabs(Xa-Xb), 1.4) << endl;
5090     */
5091     r0ab = r0a + r0b - c * pow(fabs(Xa - Xb), 1.4) - 0.008;
5092 
5093     return r0ab;
5094   }
5095 
GetParameter1Atom(int a,std::vector<OBFFParameter> & parameter)5096   OBFFParameter* OBForceFieldMMFF94::GetParameter1Atom(int a, std::vector<OBFFParameter> &parameter)
5097   {
5098     OBFFParameter *par;
5099 
5100     for (unsigned int idx=0; idx < parameter.size(); idx++)
5101       if (a == parameter[idx].a) {
5102         par = &parameter[idx];
5103         return par;
5104       }
5105 
5106     return nullptr;
5107   }
5108 
GetParameter2Atom(int a,int b,std::vector<OBFFParameter> & parameter)5109   OBFFParameter* OBForceFieldMMFF94::GetParameter2Atom(int a, int b, std::vector<OBFFParameter> &parameter)
5110   {
5111     OBFFParameter *par;
5112 
5113     for (unsigned int idx=0; idx < parameter.size(); idx++)
5114       if (((a == parameter[idx].a) && (b == parameter[idx].b)) ||
5115           ((a == parameter[idx].b) && (b == parameter[idx].a)))
5116         {
5117           par = &parameter[idx];
5118           return par;
5119         }
5120 
5121     return nullptr;
5122   }
5123 
GetParameter3Atom(int a,int b,int c,std::vector<OBFFParameter> & parameter)5124   OBFFParameter* OBForceFieldMMFF94::GetParameter3Atom(int a, int b, int c, std::vector<OBFFParameter> &parameter)
5125   {
5126     OBFFParameter *par;
5127 
5128     for (unsigned int idx=0; idx < parameter.size(); idx++)
5129       if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c)) ||
5130           ((a == parameter[idx].c) && (b == parameter[idx].b) && (c == parameter[idx].a)))
5131         {
5132           par = &parameter[idx];
5133           return par;
5134         }
5135 
5136     return nullptr;
5137   }
5138 
GetTypedParameter2Atom(int ffclass,int a,int b,std::vector<OBFFParameter> & parameter)5139   OBFFParameter* OBForceFieldMMFF94::GetTypedParameter2Atom(int ffclass, int a, int b, std::vector<OBFFParameter> &parameter)
5140   {
5141     OBFFParameter *par;
5142 
5143     for (unsigned int idx=0; idx < parameter.size(); idx++)
5144       if (((a == parameter[idx].a) && (b == parameter[idx].b) && (ffclass == parameter[idx]._ipar[0])) ||
5145           ((a == parameter[idx].b) && (b == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0])))
5146         {
5147           par = &parameter[idx];
5148           return par;
5149         }
5150 
5151     return nullptr;
5152   }
5153 
GetTypedParameter3Atom(int ffclass,int a,int b,int c,std::vector<OBFFParameter> & parameter)5154   OBFFParameter* OBForceFieldMMFF94::GetTypedParameter3Atom(int ffclass, int a, int b, int c, std::vector<OBFFParameter> &parameter)
5155   {
5156     OBFFParameter *par;
5157 
5158     for (unsigned int idx=0; idx < parameter.size(); idx++)
5159       if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c) && (ffclass == parameter[idx]._ipar[0])) ||
5160           ((a == parameter[idx].c) && (b == parameter[idx].b) && (c == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0])) )
5161         {
5162           par = &parameter[idx];
5163           return par;
5164         }
5165 
5166     return nullptr;
5167   }
5168 
GetTypedParameter4Atom(int ffclass,int a,int b,int c,int d,std::vector<OBFFParameter> & parameter)5169   OBFFParameter* OBForceFieldMMFF94::GetTypedParameter4Atom(int ffclass, int a, int b, int c, int d, std::vector<OBFFParameter> &parameter)
5170   {
5171     OBFFParameter *par;
5172 
5173     for (unsigned int idx=0; idx < parameter.size(); idx++)
5174       if (((a == parameter[idx].a) && (b == parameter[idx].b) && (c == parameter[idx].c) &&
5175            (d == parameter[idx].d) && (ffclass == parameter[idx]._ipar[0]))
5176           /* || ((a == parameter[idx].d) && (b == parameter[idx].c) && (c == parameter[idx].b) &&
5177              (d == parameter[idx].a) && (ffclass == parameter[idx]._ipar[0])) */ )
5178         {
5179           par = &parameter[idx];
5180           return par;
5181         }
5182 
5183     return nullptr;
5184   }
5185 
5186 } // end namespace OpenBabel
5187 
5188 //! \file forcefieldmmff94.cpp
5189 //! \brief MMFF94 force field
5190