1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.nonbonded;
39
40 import ffx.potential.parameters.ForceField;
41 import ffx.potential.parameters.VDWPairType;
42 import ffx.potential.parameters.VDWType;
43 import ffx.utilities.FFXProperty;
44
45 import java.util.Map;
46 import java.util.TreeMap;
47 import java.util.logging.Logger;
48
49 import static ffx.potential.parameters.ForceField.toEnumForm;
50 import static ffx.utilities.PropertyGroup.VanDerWaalsFunctionalForm;
51 import static java.lang.String.format;
52 import static org.apache.commons.math3.util.FastMath.pow;
53 import static org.apache.commons.math3.util.FastMath.sqrt;
54
55
56
57
58
59
60
61 public class VanDerWaalsForm {
62
63
64
65
66 private static final Logger logger = Logger.getLogger(VanDerWaalsForm.class.getName());
67
68
69
70
71 static final byte RADMIN = 0;
72
73
74
75 static final byte EPS = 1;
76
77
78
79
80 private static final double DEFAULT_GAMMA = 0.12;
81
82
83
84
85 private static final double DEFAULT_DELTA = 0.07;
86
87
88
89
90 private static final EPSILON_RULE DEFAULT_EPSILON_RULE = EPSILON_RULE.GEOMETRIC;
91
92
93
94
95 private static final RADIUS_RULE DEFAULT_RADIUS_RULE = RADIUS_RULE.ARITHMETIC;
96
97
98
99
100 private static final RADIUS_SIZE DEFAULT_RADIUS_SIZE = RADIUS_SIZE.RADIUS;
101
102
103
104
105 private static final RADIUS_TYPE DEFAULT_RADIUS_TYPE = RADIUS_TYPE.R_MIN;
106
107
108
109
110 private static final VDW_TYPE DEFAULT_VDW_TYPE = VDW_TYPE.LENNARD_JONES;
111
112
113
114
115 private static final double DEFAULT_VDW_12_SCALE = 0.0;
116
117
118
119
120 private static final double DEFAULT_VDW_13_SCALE = 0.0;
121
122
123
124
125 private static final double DEFAULT_VDW_14_SCALE = 1.0;
126
127
128
129
130 public static final double DEFAULT_VDW_TAPER = 0.9;
131
132
133
134
135 public static final double DEFAULT_VDW_CUTOFF = 12.0;
136
137
138
139
140 @FFXProperty(name = "halgren-gamma", propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "0.12",
141 description = """
142 Sets the value of the gamma parameter in Halgren’s buffered 14-7 vdw potential energy functional form.
143 In the absence of the gamma-halgren property, a default value of 0.12 is used."
144 """)
145
146 public final double gamma;
147
148
149
150
151 @FFXProperty(name = "halgren-delta", propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "0.07",
152 description = """
153 Sets the value of the delta parameter in Halgren’s buffered 14-7 vdw potential energy functional form.
154 In the absence of the delta-halgren property, a default value of 0.07 is used.
155 """)
156 public final double delta;
157
158
159
160
161 @FFXProperty(name = "vdwtype", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "LENNARD-JONES",
162 description = """
163 [LENNARD-JONES / BUFFERED-14-7]
164 Sets the functional form for the van der Waals potential energy term.
165 The text modifier gives the name of the functional form to be used.
166 The default in the absence of the vdwtype keyword is to use the standard two parameter Lennard-Jones function.
167 """)
168 public VDW_TYPE vdwType;
169
170
171
172
173 @FFXProperty(name = "epsilonrule", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "GEOMETRIC",
174 description = """
175 [GEOMETRIC / HHG / W-H]
176 Selects the combining rule used to derive the epsilon value for van der Waals interactions.
177 The default in the absence of the epsilonrule keyword is to use the geometric mean of the
178 individual epsilon values of the two atoms involved in the van der Waals interaction.
179 """)
180 public EPSILON_RULE epsilonRule;
181
182
183
184
185 @FFXProperty(name = "radiusrule", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "ARITHMETIC",
186 description = """
187 [ARITHMETIC / GEOMETRIC / CUBIC-MEAN]
188 Sets the functional form of the radius combining rule for heteroatomic van der Waals potential energy interactions.
189 The default in the absence of the radiusrule keyword is to use the arithmetic mean combining rule to get radii for heteroatomic interactions.
190 """)
191 public RADIUS_RULE radiusRule;
192
193
194
195
196 @FFXProperty(name = "radiussize", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "RADIUS",
197 description = """
198 [RADIUS / DIAMETER]
199 Determines whether the atom size values given in van der Waals parameters read from
200 VDW keyword statements are interpreted as atomic radius or diameter values.
201 The default in the absence of the radiussize keyword is to assume that vdw size parameters are given as radius values.
202 """)
203 public RADIUS_SIZE radiusSize;
204
205
206
207
208 @FFXProperty(name = "radiustype", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "R-MIN",
209 description = """
210 [R-MIN / SIGMA]
211 Determines whether atom size values given in van der Waals parameters read from VDW keyword
212 statements are interpreted as potential minimum (Rmin) or LJ-style sigma values.
213 The default in the absence of the radiustype keyword is to assume that vdw size parameters are given as Rmin values.
214 """)
215 public RADIUS_TYPE radiusType;
216
217
218
219
220 @FFXProperty(name = "vdw-12-scale", propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "0.0",
221 description = """
222 Provides a multiplicative scale factor that is applied to van der Waals potential
223 interactions between 1-2 connected atoms, i.e., atoms that are directly bonded.
224 The default value of 0.0 is used to omit 1-2 interactions,
225 if the vdw-12-scale property is not given in either the parameter file or the property file.
226 """)
227 protected double scale12;
228
229
230
231
232 @FFXProperty(name = "vdw-13-scale", propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "0.0",
233 description = """
234 Provides a multiplicative scale factor that is applied to van der Waals potential
235 interactions between 1-3 connected atoms, i.e., atoms separated by two covalent bonds.
236 The default value of 0.0 is used to omit 1-3 interactions, if the vdw-13-scale property
237 is not given in either the parameter file or the property file.
238 """)
239 protected double scale13;
240
241
242
243
244 @FFXProperty(name = "vdw-14-scale", propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "1.0",
245 description = """
246 Provides a multiplicative scale factor that is applied to van der Waals potential
247 interactions between 1-4 connected atoms, i.e., atoms separated by three covalent bonds.
248 The default value of 1.0 is used, if the vdw-14-scale keyword is not given in either
249 the parameter file or the property file.
250 """)
251 protected double scale14;
252
253
254
255
256 int maxClass;
257
258
259
260 double[] rMin;
261
262
263
264 double[] eps;
265
266
267
268 private final double[][] radEps;
269
270
271
272
273
274 private final double[][] radEps14;
275
276
277
278 private final VDWPowers vdwPowers;
279
280
281
282 final int dispersivePower;
283
284
285
286 private final int dispersivePower1;
287
288
289
290 final double t1n;
291
292
293
294 final int repDispPower;
295
296
297
298 private final int repDispPower1;
299
300
301
302 final double gamma1;
303
304
305
306
307
308
309 public VanDerWaalsForm(ForceField forceField) {
310
311
312 vdwType = DEFAULT_VDW_TYPE;
313 String value = forceField.getString("VDWTYPE", DEFAULT_VDW_TYPE.name());
314 try {
315 vdwType = VDW_TYPE.valueOf(toEnumForm(value));
316 } catch (Exception e) {
317 logger.info(format(" Unrecognized VDWTYPE %s; defaulting to %s.", value, vdwType));
318 }
319
320 switch (vdwType) {
321 case BUFFERED_14_7 -> vdwPowers = new Buffered_14_7();
322 case LENNARD_JONES -> vdwPowers = new LJ_6_12();
323 default -> vdwPowers = new VDWPowers();
324 }
325
326
327 epsilonRule = DEFAULT_EPSILON_RULE;
328 value = forceField.getString("EPSILONRULE", DEFAULT_EPSILON_RULE.name());
329 try {
330 epsilonRule = EPSILON_RULE.valueOf(toEnumForm(value));
331 } catch (Exception e) {
332 logger.info(format(" Unrecognized EPSILONRULE %s; defaulting to %s.", value, epsilonRule));
333 }
334
335
336 radiusRule = DEFAULT_RADIUS_RULE;
337 value = forceField.getString("RADIUSRULE", DEFAULT_RADIUS_RULE.name());
338 try {
339 radiusRule = RADIUS_RULE.valueOf(toEnumForm(value));
340 } catch (Exception e) {
341 logger.info(format(" Unrecognized RADIUSRULE %s; defaulting to %s.", value, radiusRule));
342 }
343
344
345 radiusSize = DEFAULT_RADIUS_SIZE;
346 value = forceField.getString("RADIUSSIZE", DEFAULT_RADIUS_SIZE.name());
347 try {
348 radiusSize = RADIUS_SIZE.valueOf(toEnumForm(value));
349 } catch (Exception e) {
350 logger.info(format(" Unrecognized RADIUSSIZE %s; defaulting to %s.", value, radiusSize));
351 }
352
353
354 radiusType = DEFAULT_RADIUS_TYPE;
355 value = forceField.getString("RADIUSTYPE", DEFAULT_RADIUS_TYPE.name());
356 try {
357 radiusType = RADIUS_TYPE.valueOf(toEnumForm(value));
358 } catch (Exception e) {
359 logger.info(format(" Unrecognized RADIUSTYPE %s; defaulting to %s", value, radiusType));
360 }
361
362
363 int repulsivePower;
364 switch (vdwType) {
365 case LENNARD_JONES -> {
366 repulsivePower = 12;
367 dispersivePower = 6;
368 delta = 0.0;
369 gamma = 0.0;
370 }
371 default -> {
372 repulsivePower = 14;
373 dispersivePower = 7;
374 delta = forceField.getDouble("DELTA-HALGREN", DEFAULT_DELTA);
375 gamma = forceField.getDouble("GAMMA-HALGREN", DEFAULT_GAMMA);
376 }
377 }
378
379 repDispPower = repulsivePower - dispersivePower;
380 dispersivePower1 = dispersivePower - 1;
381 repDispPower1 = repDispPower - 1;
382 double delta1 = 1.0 + delta;
383 t1n = pow(delta1, dispersivePower);
384 gamma1 = 1.0 + gamma;
385
386 scale12 = forceField.getDouble("VDW_12_SCALE", DEFAULT_VDW_12_SCALE);
387 scale13 = forceField.getDouble("VDW_13_SCALE", DEFAULT_VDW_13_SCALE);
388 scale14 = forceField.getDouble("VDW_14_SCALE", DEFAULT_VDW_14_SCALE);
389 double scale15 = forceField.getDouble("VDW_15_SCALE", 1.0);
390 if (scale15 != 1.0) {
391 logger.severe(" Van Der Waals 1-5 masking rules are not supported.");
392 }
393
394
395 if (scale12 > 1.0) {
396 scale12 = 1.0 / scale12;
397 }
398 if (scale13 > 1.0) {
399 scale13 = 1.0 / scale13;
400 }
401 if (scale14 > 1.0) {
402 scale14 = 1.0 / scale14;
403 }
404
405 Map<String, VDWType> map = forceField.getVDWTypes();
406 TreeMap<String, VDWType> vdwTypes = new TreeMap<>(map);
407 maxClass = 0;
408 for (VDWType currentType : vdwTypes.values()) {
409 if (currentType.atomClass > maxClass) {
410 maxClass = currentType.atomClass;
411 }
412 }
413 radEps = new double[maxClass + 1][2 * (maxClass + 1)];
414 radEps14 = new double[maxClass + 1][2 * (maxClass + 1)];
415
416
417 double radScale = switch (radiusSize) {
418 case DIAMETER -> 0.5;
419 default -> 1.0;
420 };
421 switch (radiusType) {
422 case SIGMA -> radScale *= 1.122462048309372981;
423 default -> {
424 }
425 }
426
427 rMin = new double[maxClass + 1];
428 eps = new double[maxClass + 1];
429
430
431 for (VDWType vdwi : vdwTypes.values()) {
432 int i = vdwi.atomClass;
433 double ri = radScale * vdwi.radius;
434 double e1 = vdwi.wellDepth;
435 rMin[i] = ri;
436 eps[i] = e1;
437 for (VDWType vdwj : vdwTypes.tailMap(vdwi.getKey()).values()) {
438 int j = vdwj.atomClass;
439 double rj = radScale * vdwj.radius;
440 double e2 = vdwj.wellDepth;
441 double radmin = getCombinedRadius(ri, rj, radiusRule);
442 double eps = getCombinedEps(e1, e2, ri, rj, epsilonRule);
443 if (radmin > 0) {
444 radEps[i][j * 2 + RADMIN] = 1.0 / radmin;
445 radEps[j][i * 2 + RADMIN] = 1.0 / radmin;
446 radEps14[i][j * 2 + RADMIN] = 1.0 / radmin;
447 radEps14[j][i * 2 + RADMIN] = 1.0 / radmin;
448 } else {
449 radEps[i][j * 2 + RADMIN] = 0.0;
450 radEps[j][i * 2 + RADMIN] = 0.0;
451 radEps14[i][j * 2 + RADMIN] = 0.0;
452 radEps14[j][i * 2 + RADMIN] = 0.0;
453 }
454 radEps[i][j * 2 + EPS] = eps;
455 radEps[j][i * 2 + EPS] = eps;
456 radEps14[i][j * 2 + EPS] = eps;
457 radEps14[j][i * 2 + EPS] = eps;
458 }
459 }
460
461
462 Map<String, VDWType> vdw14Types = forceField.getVDW14Types();
463 for (VDWType vdwi : vdwTypes.values()) {
464
465 VDWType vdw14 = forceField.getVDW14Type(vdwi.getKey());
466 if (vdw14 != null) {
467 vdwi = vdw14;
468 }
469 int i = vdwi.atomClass;
470 double ri = radScale * vdwi.radius;
471 double e1 = vdwi.wellDepth;
472
473 for (VDWType vdwj : vdw14Types.values()) {
474 int j = vdwj.atomClass;
475 double rj = radScale * vdwj.radius;
476 double e2 = vdwj.wellDepth;
477 double radmin = getCombinedRadius(ri, rj, radiusRule);
478 double eps = getCombinedEps(e1, e2, ri, rj, epsilonRule);
479 if (radmin > 0) {
480 radEps14[i][j * 2 + RADMIN] = 1.0 / radmin;
481 radEps14[j][i * 2 + RADMIN] = 1.0 / radmin;
482 } else {
483 radEps14[i][j * 2 + RADMIN] = 0.0;
484 radEps14[j][i * 2 + RADMIN] = 0.0;
485 }
486 radEps14[i][j * 2 + EPS] = eps;
487 radEps14[j][i * 2 + EPS] = eps;
488 }
489 }
490
491
492 Map<String, VDWPairType> vdwPairTypes = forceField.getVDWPairTypes();
493 for (VDWPairType vdwPairType : vdwPairTypes.values()) {
494 int i = vdwPairType.atomClasses[0];
495 int j = vdwPairType.atomClasses[1];
496 double radmin = vdwPairType.radius;
497 double eps = vdwPairType.wellDepth;
498 if (radmin > 0) {
499 radEps[i][j * 2 + RADMIN] = 1.0 / radmin;
500 radEps[j][i * 2 + RADMIN] = 1.0 / radmin;
501 radEps14[i][j * 2 + RADMIN] = 1.0 / radmin;
502 radEps14[j][i * 2 + RADMIN] = 1.0 / radmin;
503 } else {
504 radEps[i][j * 2 + RADMIN] = 0.0;
505 radEps[j][i * 2 + RADMIN] = 0.0;
506 radEps14[i][j * 2 + RADMIN] = 0.0;
507 radEps14[j][i * 2 + RADMIN] = 0.0;
508 }
509 radEps[i][j * 2 + EPS] = eps;
510 radEps[j][i * 2 + EPS] = eps;
511 radEps14[i][j * 2 + EPS] = eps;
512 radEps14[j][i * 2 + EPS] = eps;
513 }
514 }
515
516
517
518
519
520
521
522
523
524
525
526 public static double getCombinedEps(double ei, double ej, double ri, double rj,
527 EPSILON_RULE epsilonRule) {
528 double sei = sqrt(ei);
529 double sej = sqrt(ej);
530 switch (epsilonRule) {
531 case GEOMETRIC -> {
532 return sei * sej;
533 }
534 case W_H -> {
535
536 double rr = ri * rj;
537 double rr3 = rr * rr * rr;
538 double ri6 = pow(ri, 6);
539 double rj6 = pow(rj, 6);
540 return 2.0 * sei * sej * rr3 / (ri6 + rj6);
541 }
542 default -> {
543 return 4.0 * (ei * ej) / ((sei + sej) * (sei + sej));
544 }
545 }
546 }
547
548
549
550
551
552
553
554
555
556 public static double getCombinedRadius(double ri, double rj, RADIUS_RULE radiusRule) {
557 switch (radiusRule) {
558 case ARITHMETIC -> {
559 return ri + rj;
560 }
561 case GEOMETRIC -> {
562 return 2.0 * sqrt(ri) * sqrt(rj);
563 }
564 default -> {
565 double ri2 = ri * ri;
566 double ri3 = ri * ri2;
567 double rj2 = rj * rj;
568 double rj3 = rj * rj2;
569 return 2.0 * (ri3 + rj3) / (ri2 + rj2);
570 }
571 }
572 }
573
574
575
576
577
578
579
580
581 public double getCombinedEps(int class1, int class2) {
582 return radEps[class1][class2 * 2 + EPS];
583 }
584
585
586
587
588
589
590
591
592 public double getCombinedEps14(int class1, int class2) {
593 return radEps14[class1][class2 * 2 + EPS];
594 }
595
596
597
598
599
600
601
602
603 public double getCombinedInverseRmin(int class1, int class2) {
604 return radEps[class1][class2 * 2 + RADMIN];
605 }
606
607
608
609
610
611
612
613
614 public double getCombinedInverseRmin14(int class1, int class2) {
615 return radEps14[class1][class2 * 2 + RADMIN];
616 }
617
618
619
620
621
622
623 public double[] getEps() {
624 return eps;
625 }
626
627
628
629
630
631
632 public double[] getRmin() {
633 return rMin;
634 }
635
636
637
638
639
640
641 public double getScale14() {
642 return scale14;
643 }
644
645
646
647
648
649
650
651 double rhoDisp1(double rho) {
652 return vdwPowers.rhoDisp1(rho);
653 }
654
655
656
657
658
659
660
661 double rhoDelta1(double rhoDelta) {
662 return vdwPowers.rhoDisp1(rhoDelta);
663 }
664
665
666
667
668 public enum VDW_TYPE {
669 BUFFERED_14_7,
670 LENNARD_JONES
671 }
672
673
674
675
676 public enum RADIUS_RULE {
677 ARITHMETIC,
678 CUBIC_MEAN,
679 GEOMETRIC
680 }
681
682
683
684
685 public enum RADIUS_SIZE {
686 DIAMETER,
687 RADIUS
688 }
689
690
691
692
693 public enum RADIUS_TYPE {
694 R_MIN,
695 SIGMA
696 }
697
698
699
700
701 public enum EPSILON_RULE {
702 GEOMETRIC,
703 HHG,
704 W_H
705 }
706
707
708
709
710 private class VDWPowers {
711
712 public double rhoDelta1(double rhoDelta) {
713 return pow(rhoDelta, repDispPower1);
714 }
715
716 public double rhoDisp1(double rho) {
717 return pow(rho, dispersivePower1);
718 }
719 }
720
721
722
723
724 private class LJ_6_12 extends VDWPowers {
725
726 @Override
727 public double rhoDelta1(double rhoDelta) {
728 double rhoDelta2 = rhoDelta * rhoDelta;
729 return rhoDelta2 * rhoDelta2 * rhoDelta;
730 }
731
732 @Override
733 public double rhoDisp1(double rho) {
734 double rho2 = rho * rho;
735 return rho2 * rho2 * rho;
736 }
737 }
738
739
740
741
742 private class Buffered_14_7 extends VDWPowers {
743
744 @Override
745 public double rhoDelta1(double rhoDelta) {
746 double rhoDelta2 = rhoDelta * rhoDelta;
747 return rhoDelta2 * rhoDelta2 * rhoDelta2;
748 }
749
750 @Override
751 public double rhoDisp1(double rho) {
752 double rho2 = rho * rho;
753 return rho2 * rho2 * rho2;
754 }
755 }
756 }