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 edu.rit.pj.BarrierAction;
41 import edu.rit.pj.IntegerForLoop;
42 import edu.rit.pj.IntegerSchedule;
43 import edu.rit.pj.ParallelRegion;
44 import edu.rit.pj.ParallelTeam;
45 import edu.rit.pj.reduction.SharedDouble;
46 import edu.rit.pj.reduction.SharedInteger;
47 import ffx.crystal.Crystal;
48 import ffx.crystal.SymOp;
49 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
50 import ffx.numerics.atomic.AtomicDoubleArray3D;
51 import ffx.numerics.switching.MultiplicativeSwitch;
52 import ffx.potential.bonded.Atom;
53 import ffx.potential.bonded.Bond;
54 import ffx.potential.bonded.LambdaInterface;
55 import ffx.potential.extended.ExtendedSystem;
56 import ffx.potential.parameters.AtomType;
57 import ffx.potential.parameters.ForceField;
58 import ffx.potential.parameters.VDWType;
59 import ffx.utilities.FFXProperty;
60
61 import java.util.Arrays;
62 import java.util.List;
63 import java.util.logging.Level;
64 import java.util.logging.Logger;
65
66 import static ffx.potential.parameters.ForceField.toEnumForm;
67 import static ffx.utilities.PropertyGroup.NonBondedCutoff;
68 import static ffx.utilities.PropertyGroup.VanDerWaalsFunctionalForm;
69 import static java.lang.Double.isNaN;
70 import static java.lang.String.format;
71 import static java.util.Arrays.fill;
72 import static org.apache.commons.math3.util.FastMath.PI;
73 import static org.apache.commons.math3.util.FastMath.max;
74 import static org.apache.commons.math3.util.FastMath.min;
75 import static org.apache.commons.math3.util.FastMath.pow;
76 import static org.apache.commons.math3.util.FastMath.sqrt;
77
78
79
80
81
82
83
84
85
86
87 public class VanDerWaals implements MaskingInterface, LambdaInterface {
88
89 private static final Logger logger = Logger.getLogger(VanDerWaals.class.getName());
90 private static final byte HARD = 0;
91 private static final byte SOFT = 1;
92 private static final byte XX = 0;
93 private static final byte YY = 1;
94 private static final byte ZZ = 2;
95 private final boolean doLongRangeCorrection;
96
97
98
99 private final ParallelTeam parallelTeam;
100 private final int threadCount;
101 private final IntegerSchedule pairwiseSchedule;
102 private final SharedInteger sharedInteractions;
103 private final SharedDouble sharedEnergy;
104 private final SharedDouble shareddEdL;
105 private final SharedDouble sharedd2EdL2;
106 private final VanDerWaalsRegion vanDerWaalsRegion;
107 private final VanDerWaalsForm vdwForm;
108 @FFXProperty(name = "vdwindex", clazz = String.class, propertyGroup = VanDerWaalsFunctionalForm, defaultValue = "class",
109 description = """
110 [CLASS / TYPE]
111 Specifies whether van der Waals parameters are provided for atom classes or atom types.
112 While most force fields are indexed by atom class, in OPLS models the vdW values are indexed by atom type.
113 The default in the absence of the vdwindex property is to index vdW parameters by atom class.
114 """)
115 private final String vdwIndex;
116 @FFXProperty(name = "vdw-taper", propertyGroup = NonBondedCutoff, defaultValue = "0.9", description = """
117 Allows modification of the cutoff windows for van der Waals potential energy interactions. "
118 The default value in the absence of the vdw-taper keyword is to begin the cutoff window
119 at 0.9 of the vdw cutoff distance.
120 """)
121 private final double vdwTaper;
122
123 private final NonbondedCutoff nonbondedCutoff;
124 private final MultiplicativeSwitch multiplicativeSwitch;
125
126
127
128 private Crystal crystal;
129
130
131
132 private Atom[] atoms;
133
134
135
136 private int[] molecule;
137
138
139
140 private boolean[] neuralNetwork;
141
142
143
144 private ForceField forceField;
145
146
147
148 private boolean[] use = null;
149
150
151
152 private int nAtoms;
153
154
155
156 private int nSymm;
157
158
159
160 private boolean gradient;
161 private boolean lambdaTerm;
162 private boolean esvTerm;
163 private boolean[] isSoft;
164
165
166
167
168
169
170
171
172
173 private boolean[][] softCore;
174
175
176
177 private boolean intermolecularSoftcore = false;
178
179
180
181
182
183
184 private boolean intramolecularSoftcore = false;
185
186
187
188 private double lambda = 1.0;
189
190
191
192 private double vdwLambdaEnd = 1.0;
193
194
195
196 private double vdwLambdaExponent = 3.0;
197
198
199
200 private double vdwLambdaAlpha = 0.25;
201
202
203
204
205 private LambdaFactors[] lambdaFactors = null;
206
207
208
209 private double sc1 = 0.0;
210
211
212
213 private double sc2 = 1.0;
214 private double dsc1dL = 0.0;
215 private double dsc2dL = 0.0;
216 private double d2sc1dL2 = 0.0;
217 private double d2sc2dL2 = 0.0;
218
219
220
221 private ExtendedSystem esvSystem;
222
223
224
225
226 private double[] coordinates;
227
228
229
230 private double[][] reduced;
231
232 private double[] reducedXYZ;
233
234
235
236 private int[][][] neighborLists;
237
238
239
240 private int[] atomClass;
241
242
243
244
245 private int[] reductionIndex;
246
247 private int[][] mask12;
248 private int[][] mask13;
249 private int[][] mask14;
250
251
252
253
254 private double[] reductionValue;
255
256 private boolean reducedHydrogen;
257 private double longRangeCorrection;
258 private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
259
260
261
262 private AtomicDoubleArray3D grad;
263
264
265
266 private AtomicDoubleArray3D lambdaGrad;
267
268
269
270
271 private NeighborList neighborList;
272
273
274
275 private boolean forceNeighborListRebuild = true;
276
277 public VanDerWaals() {
278
279 doLongRangeCorrection = false;
280 parallelTeam = null;
281 threadCount = 0;
282 pairwiseSchedule = null;
283 sharedInteractions = null;
284 sharedEnergy = null;
285 shareddEdL = null;
286 sharedd2EdL2 = null;
287 vanDerWaalsRegion = null;
288 vdwForm = null;
289 nonbondedCutoff = null;
290 multiplicativeSwitch = null;
291 vdwIndex = null;
292 vdwTaper = VanDerWaalsForm.DEFAULT_VDW_TAPER;
293 }
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308 public VanDerWaals(
309 Atom[] atoms,
310 int[] molecule,
311 boolean[] neuralNetwork,
312 Crystal crystal,
313 ForceField forceField,
314 ParallelTeam parallelTeam,
315 double vdwCutoff,
316 double neighborListCutoff) {
317 this.atoms = atoms;
318 this.molecule = molecule;
319 this.neuralNetwork = neuralNetwork;
320 this.crystal = crystal;
321 this.parallelTeam = parallelTeam;
322 this.forceField = forceField;
323 nAtoms = atoms.length;
324 nSymm = crystal.spaceGroup.getNumberOfSymOps();
325 vdwForm = new VanDerWaalsForm(forceField);
326
327 vdwIndex = forceField.getString("VDWINDEX", "Class");
328 reducedHydrogen = forceField.getBoolean("REDUCE_HYDROGENS", true);
329
330
331 lambdaTerm = forceField.getBoolean("VDW_LAMBDATERM",
332 forceField.getBoolean("LAMBDATERM", false));
333 if (lambdaTerm) {
334 shareddEdL = new SharedDouble();
335 sharedd2EdL2 = new SharedDouble();
336 vdwLambdaAlpha = forceField.getDouble("VDW_LAMBDA_ALPHA", 0.25);
337 vdwLambdaExponent = forceField.getDouble("VDW_LAMBDA_EXPONENT", 3.0);
338 vdwLambdaEnd = forceField.getDouble("VDW_LAMBDA_END", 1.0);
339 if (vdwLambdaAlpha < 0.0) {
340 logger.warning(format(
341 " Invalid value %8.3g for vdw-lambda-alpha; must be greater than or equal to 0. Resetting to 0.25.",
342 vdwLambdaAlpha));
343 vdwLambdaAlpha = 0.25;
344 }
345 if (vdwLambdaExponent < 1.0) {
346 logger.warning(format(
347 " Invalid value %8.3g for vdw-lambda-exponent; must be greater than or equal to 1. Resetting to 3.",
348 vdwLambdaExponent));
349 vdwLambdaExponent = 3.0;
350 }
351
352 vdwLambdaEnd = forceField.getDouble("VDW_LAMBDA_END", 1.0);
353 if (vdwLambdaEnd < 0.0 || vdwLambdaEnd > 1.0) {
354 logger.warning(format(
355 " Invalid value %8.3g for vdw-lambda-end; must be between 0 and 1. Resetting to 1.", vdwLambdaEnd));
356 vdwLambdaEnd = 1.0;
357 }
358 intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
359 intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
360 } else {
361 shareddEdL = null;
362 sharedd2EdL2 = null;
363 }
364
365
366 threadCount = parallelTeam.getThreadCount();
367 sharedInteractions = new SharedInteger();
368 sharedEnergy = new SharedDouble();
369 doLongRangeCorrection = forceField.getBoolean("VDW_CORRECTION", false);
370 vanDerWaalsRegion = new VanDerWaalsRegion();
371
372
373 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
374 String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
375 try {
376 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
377 } catch (Exception e) {
378 logger.info(format(
379 " Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value, atomicDoubleArrayImpl));
380 }
381
382
383 initAtomArrays();
384
385
386
387
388
389
390 double taper = forceField.getDouble("VDW_TAPER", VanDerWaalsForm.DEFAULT_VDW_TAPER);
391 vdwTaper = taper * vdwCutoff;
392 double buff = 2.0;
393 nonbondedCutoff = new NonbondedCutoff(vdwCutoff, vdwTaper, buff);
394 multiplicativeSwitch = new MultiplicativeSwitch(vdwTaper, vdwCutoff);
395 neighborList = new NeighborList(null, this.crystal,
396 atoms, neighborListCutoff, buff, parallelTeam);
397 pairwiseSchedule = neighborList.getPairwiseSchedule();
398 neighborLists = new int[nSymm][][];
399
400
401 buildNeighborList(atoms);
402
403
404 neighborList.setDisableUpdates(forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
405
406 logger.info(toString());
407 }
408
409
410
411
412
413
414 @Override
415 public void applyMask(final int i, final boolean[] is14, final double[]... masks) {
416 double[] mask = masks[0];
417 var m12 = mask12[i];
418 for (int value : m12) {
419 mask[value] = vdwForm.scale12;
420 }
421 var m13 = mask13[i];
422 for (int value : m13) {
423 mask[value] = vdwForm.scale13;
424 }
425 var m14 = mask14[i];
426 for (int value : m14) {
427 mask[value] = vdwForm.scale14;
428 is14[value] = true;
429 }
430 }
431
432
433
434
435
436
437 public void attachExtendedSystem(ExtendedSystem system) {
438 if (system == null) {
439 logger.severe("Tried to attach null extended system.");
440 }
441 esvTerm = true;
442 esvSystem = system;
443
444
445 vdwLambdaAlpha = forceField.getDouble("VDW_LAMBDA_ALPHA", 0.05);
446 vdwLambdaExponent = forceField.getDouble("VDW_LAMBDA_EXPONENT", 1.0);
447 if (vdwLambdaExponent != 1.0) {
448 logger.warning(
449 format(
450 "ESVs are compatible only with a vdwLambdaExponent of unity!"
451 + " (found %g, resetting to 1.0)",
452 vdwLambdaExponent));
453 vdwLambdaExponent = 1.0;
454 }
455 if (vdwLambdaAlpha < 0.0) {
456 vdwLambdaAlpha = 0.05;
457 }
458
459 intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
460 intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
461
462 Atom[] atomsExt = esvSystem.getExtendedAtoms();
463 int[] moleculeExt = esvSystem.getExtendedMolecule();
464
465
466 boolean[] nn = new boolean[atomsExt.length];
467 Arrays.fill(nn, false);
468 for (Atom a : atomsExt) {
469 if (a.isNeuralNetwork()) {
470 logger.severe(
471 "The combination of an extended system and neural networks are not supported.");
472 }
473 }
474
475 setAtoms(atomsExt, moleculeExt, nn);
476 }
477
478
479
480
481
482
483 public ExtendedSystem getExtendedSystem() {
484 return esvSystem;
485 }
486
487
488
489
490
491
492 public void destroy() throws Exception {
493 if (neighborList != null) {
494 neighborList.destroy();
495 }
496 }
497
498
499
500
501
502
503
504
505
506 public double energy(boolean gradient, boolean print) {
507 this.gradient = gradient;
508
509 try {
510 parallelTeam.execute(vanDerWaalsRegion);
511 } catch (Exception e) {
512 String message = " Fatal exception expanding coordinates.\n";
513 logger.log(Level.SEVERE, message, e);
514 }
515 return sharedEnergy.get();
516 }
517
518
519
520
521
522 public boolean getLambdaTerm() {
523 return lambdaTerm;
524 }
525
526
527
528
529 public boolean getIntramolecularSoftcore() {
530 return intramolecularSoftcore;
531 }
532
533
534
535
536
537
538 public double getAlpha() {
539 return vdwLambdaAlpha;
540 }
541
542
543
544
545
546
547 public double getBeta() {
548 return vdwLambdaExponent;
549 }
550
551
552
553
554
555
556
557 public double getBuffer() {
558 return nonbondedCutoff.buff;
559 }
560
561
562
563
564
565
566 public boolean getDoLongRangeCorrection() {
567 return doLongRangeCorrection;
568 }
569
570
571
572
573
574
575
576 public double getEnergy() {
577 return sharedEnergy.get();
578 }
579
580
581
582
583
584
585
586 public int getInteractions() {
587 return sharedInteractions.get();
588 }
589
590
591
592
593 @Override
594 public double getLambda() {
595 return lambda;
596 }
597
598
599
600
601 @Override
602 public void setLambda(double lambda) {
603 assert (lambda >= 0.0 && lambda <= 1.0);
604 if (!lambdaTerm) {
605 return;
606 }
607
608 if (lambda > vdwLambdaEnd) {
609 lambda = 1.0;
610 } else {
611 lambda = lambda / vdwLambdaEnd;
612 }
613 this.lambda = lambda;
614
615
616 sc1 = vdwLambdaAlpha * (1.0 - lambda) * (1.0 - lambda);
617
618 dsc1dL = -2.0 * vdwLambdaAlpha * (1.0 - lambda);
619
620 d2sc1dL2 = 2.0 * vdwLambdaAlpha;
621
622 sc2 = pow(lambda, vdwLambdaExponent);
623
624 dsc2dL = vdwLambdaExponent * pow(lambda, vdwLambdaExponent - 1.0);
625 if (vdwLambdaExponent >= 2.0) {
626
627 d2sc2dL2 = vdwLambdaExponent * (vdwLambdaExponent - 1.0) * pow(lambda, vdwLambdaExponent - 2.0);
628 } else {
629 d2sc2dL2 = 0.0;
630 }
631
632
633 if (!esvTerm) {
634 for (LambdaFactors lf : lambdaFactors) {
635 lf.setFactors();
636 }
637 }
638
639 initSoftCore();
640
641
642 if (doLongRangeCorrection) {
643 longRangeCorrection = computeLongRangeCorrection();
644 } else {
645 longRangeCorrection = 0.0;
646 }
647 }
648
649
650
651
652
653
654 public int[][] getMask12() {
655 return mask12;
656 }
657
658
659
660
661
662
663 public int[][] getMask13() {
664 return mask13;
665 }
666
667
668
669
670
671
672 public int[][] getMask14() {
673 return mask14;
674 }
675
676
677
678
679
680
681 public NeighborList getNeighborList() {
682 return neighborList;
683 }
684
685
686
687
688
689
690 public int[][][] getNeighborLists() {
691 return neighborLists;
692 }
693
694
695
696
697
698
699 public NonbondedCutoff getNonbondedCutoff() {
700 return nonbondedCutoff;
701 }
702
703
704
705
706
707
708 public int[] getReductionIndex() {
709 return reductionIndex;
710 }
711
712
713
714
715
716
717 public VanDerWaalsForm getVDWForm() {
718 return vdwForm;
719 }
720
721
722
723
724 @Override
725 public double getd2EdL2() {
726 if (sharedd2EdL2 == null || !lambdaTerm) {
727 return 0.0;
728 }
729 return sharedd2EdL2.get();
730 }
731
732
733
734
735 @Override
736 public double getdEdL() {
737 if (shareddEdL == null || !lambdaTerm) {
738 return 0.0;
739 }
740 return shareddEdL.get();
741 }
742
743
744
745
746 @Override
747 public void getdEdXdL(double[] lambdaGradient) {
748 if (lambdaGrad == null || !lambdaTerm) {
749 return;
750 }
751 int index = 0;
752 for (int i = 0; i < nAtoms; i++) {
753 if (atoms[i].isActive()) {
754 lambdaGradient[index++] += lambdaGrad.getX(i);
755 lambdaGradient[index++] += lambdaGrad.getY(i);
756 lambdaGradient[index++] += lambdaGrad.getZ(i);
757 }
758 }
759 }
760
761
762
763
764
765
766 @Override
767 public void removeMask(final int i, final boolean[] is14, final double[]... masks) {
768 double[] mask = masks[0];
769 var m12 = mask12[i];
770 for (int value : m12) {
771 mask[value] = 1.0;
772 }
773 var m13 = mask13[i];
774 for (int value : m13) {
775 mask[value] = 1.0;
776 }
777 var m14 = mask14[i];
778 for (int value : m14) {
779 mask[value] = 1.0;
780 is14[value] = false;
781 }
782 }
783
784
785
786
787
788
789
790
791
792 public void setAtoms(Atom[] atoms, int[] molecule, boolean[] neuralNetwork) {
793 this.atoms = atoms;
794 this.nAtoms = atoms.length;
795 this.molecule = molecule;
796 this.neuralNetwork = neuralNetwork;
797
798 if (nAtoms != molecule.length) {
799 logger.warning("Atom and molecule arrays are of different lengths.");
800 throw new IllegalArgumentException();
801 }
802
803 initAtomArrays();
804 buildNeighborList(atoms);
805 }
806
807
808
809
810
811
812
813
814 public void setCrystal(Crystal crystal) {
815 this.crystal = crystal;
816 int newNSymm = crystal.spaceGroup.getNumberOfSymOps();
817 if (nSymm != newNSymm) {
818 nSymm = newNSymm;
819
820
821 if (reduced == null || reduced.length < nSymm) {
822 reduced = new double[nSymm][nAtoms * 3];
823 reducedXYZ = reduced[0];
824 neighborLists = new int[nSymm][][];
825 }
826 }
827 neighborList.setCrystal(crystal);
828 forceNeighborListRebuild = true;
829 try {
830 parallelTeam.execute(vanDerWaalsRegion);
831 } catch (Exception e) {
832 String message = " Fatal exception expanding coordinates.\n";
833 logger.log(Level.SEVERE, message, e);
834 }
835 }
836
837
838
839
840 @Override
841 public String toString() {
842 StringBuffer sb = new StringBuffer("\n Van der Waals\n");
843 if (multiplicativeSwitch.getSwitchStart() != Double.POSITIVE_INFINITY) {
844 sb.append(format(" Switch Start: %6.3f (A)\n",
845 multiplicativeSwitch.getSwitchStart()));
846 sb.append(format(" Cut-Off: %6.3f (A)\n",
847 multiplicativeSwitch.getSwitchEnd()));
848 } else {
849 sb.append(" Cut-Off: NONE\n");
850 }
851
852 sb.append(format(" Long-Range Correction: %6B\n", doLongRangeCorrection));
853 if (!reducedHydrogen) {
854 sb.append(format(" Reduce Hydrogen: %6B\n", reducedHydrogen));
855 }
856 if (lambdaTerm) {
857 sb.append(" Alchemical Parameters\n");
858 sb.append(format(" Softcore Alpha: %5.3f\n", vdwLambdaAlpha));
859 sb.append(format(" Lambda Exponent: %5.3f\n", vdwLambdaExponent));
860 sb.append(format(" Lambda End: %5.3f\n", vdwLambdaEnd));
861 }
862 return sb.toString();
863 }
864
865
866
867
868 private void initAtomArrays() {
869 if (esvTerm) {
870 atoms = esvSystem.getExtendedAtoms();
871 nAtoms = atoms.length;
872 }
873 if (atomClass == null || nAtoms > atomClass.length || lambdaTerm || esvTerm) {
874 atomClass = new int[nAtoms];
875 coordinates = new double[nAtoms * 3];
876 reduced = new double[nSymm][nAtoms * 3];
877 reducedXYZ = reduced[0];
878 reductionIndex = new int[nAtoms];
879 reductionValue = new double[nAtoms];
880 mask12 = new int[nAtoms][];
881 mask13 = new int[nAtoms][];
882 mask14 = new int[nAtoms][];
883 use = new boolean[nAtoms];
884 isSoft = new boolean[nAtoms];
885 softCore = new boolean[2][nAtoms];
886 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
887 if (lambdaTerm) {
888 lambdaGrad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
889 } else {
890 lambdaGrad = null;
891 }
892 }
893
894
895 fill(use, true);
896 fill(isSoft, false);
897 fill(softCore[HARD], false);
898 fill(softCore[SOFT], false);
899
900 lambdaFactors = new LambdaFactors[threadCount];
901 for (int i = 0; i < threadCount; i++) {
902 if (lambdaTerm) {
903 lambdaFactors[i] = new LambdaFactorsOST();
904 } else {
905 lambdaFactors[i] = new LambdaFactors();
906 }
907 }
908
909 for (int i = 0; i < nAtoms; i++) {
910 Atom ai = atoms[i];
911 assert (i == ai.getXyzIndex() - 1);
912 double[] xyz = ai.getXYZ(null);
913 int i3 = i * 3;
914 coordinates[i3 + XX] = xyz[XX];
915 coordinates[i3 + YY] = xyz[YY];
916 coordinates[i3 + ZZ] = xyz[ZZ];
917
918 VDWType vdwType = ai.getVDWType();
919 if (vdwType == null) {
920
921 AtomType atomType = ai.getAtomType();
922 if (atomType == null) {
923 logger.severe(ai.toString());
924 return;
925 }
926
927 if (vdwIndex.equalsIgnoreCase("Type")) {
928 atomClass[i] = atomType.type;
929 } else {
930 atomClass[i] = atomType.atomClass;
931 }
932 vdwType = forceField.getVDWType(Integer.toString(atomClass[i]));
933
934 if (vdwType == null) {
935 logger.info(" No VdW type for atom class " + atomClass[i]);
936 logger.severe(" No VdW type for atom " + ai);
937 return;
938 }
939 ai.setVDWType(vdwType);
940 }
941
942 atomClass[i] = vdwType.atomClass;
943
944 List<Bond> bonds = ai.getBonds();
945 int numBonds = bonds.size();
946 if (reducedHydrogen && vdwType.reductionFactor > 0.0 && numBonds == 1) {
947 Bond bond = bonds.get(0);
948 Atom heavyAtom = bond.get1_2(ai);
949
950 reductionIndex[i] = heavyAtom.getIndex() - 1;
951 reductionValue[i] = vdwType.reductionFactor;
952 } else {
953 reductionIndex[i] = i;
954 reductionValue[i] = 0.0;
955 }
956
957
958 List<Atom> n12 = ai.get12List();
959 mask12[i] = new int[n12.size()];
960 int j = 0;
961 for (Atom a12 : n12) {
962 mask12[i][j++] = a12.getIndex() - 1;
963 }
964
965
966 List<Atom> n13 = ai.get13List();
967 mask13[i] = new int[n13.size()];
968 j = 0;
969 for (Atom a13 : n13) {
970 mask13[i][j++] = a13.getIndex() - 1;
971 }
972
973
974 List<Atom> n14 = ai.get14List();
975 mask14[i] = new int[n14.size()];
976 j = 0;
977 for (Atom a14 : n14) {
978 mask14[i][j++] = a14.getIndex() - 1;
979 }
980 }
981 }
982
983
984
985
986
987
988 private void buildNeighborList(Atom[] atoms) {
989 neighborList.setAtoms(atoms);
990 forceNeighborListRebuild = true;
991 try {
992 parallelTeam.execute(vanDerWaalsRegion);
993 } catch (Exception e) {
994 String message = " Fatal exception expanding coordinates.\n";
995 logger.log(Level.SEVERE, message, e);
996 }
997 forceNeighborListRebuild = false;
998 }
999
1000
1001
1002
1003
1004
1005
1006
1007 private double computeLongRangeCorrection() {
1008
1009 if (esvTerm) {
1010 throw new UnsupportedOperationException();
1011 }
1012
1013
1014 if (crystal.aperiodic()) {
1015 return 0.0;
1016 }
1017
1018
1019
1020
1021
1022
1023 int stepsPerAngstrom = 2;
1024 double range = 100.0;
1025 int nDelta = (int) ((double) stepsPerAngstrom * (range - nonbondedCutoff.cut));
1026 double rDelta = (range - nonbondedCutoff.cut) / (double) nDelta;
1027 double offset = nonbondedCutoff.cut - 0.5 * rDelta;
1028 double oneMinLambda = 1.0 - lambda;
1029
1030
1031 int maxClass = vdwForm.maxClass;
1032 int[] radCount = new int[maxClass + 1];
1033 int[] softRadCount = new int[maxClass + 1];
1034 fill(radCount, 0);
1035 fill(softRadCount, 0);
1036 for (int i = 0; i < nAtoms; i++) {
1037 radCount[atomClass[i]]++;
1038 if (isSoft[i]) {
1039 softRadCount[atomClass[i]]++;
1040 }
1041 }
1042
1043 double total = 0.0;
1044
1045 for (int i = 1; i < maxClass + 1; i++) {
1046 for (int j = i; j < maxClass + 1; j++) {
1047 if (radCount[i] == 0 || radCount[j] == 0) {
1048 continue;
1049 }
1050 double irv = vdwForm.getCombinedInverseRmin(i, j);
1051 double ev = vdwForm.getCombinedEps(i, j);
1052 if (isNaN(irv) || irv == 0 || isNaN(ev)) {
1053 continue;
1054 }
1055 double sume = 0.0;
1056 for (int k = 1; k <= nDelta; k++) {
1057 double r = offset + (double) k * rDelta;
1058 double r2 = r * r;
1059 final double rho = r * irv;
1060 final double rho3 = rho * rho * rho;
1061 final double rhod = rho + vdwForm.delta;
1062 final double rhod3 = rhod * rhod * rhod;
1063 double t1, t2;
1064 switch (vdwForm.vdwType) {
1065 case LENNARD_JONES -> {
1066 final double rho6 = rho3 * rho3;
1067 final double rhod6 = rhod3 * rhod3;
1068 t1 = vdwForm.t1n / rhod6;
1069 t2 = vdwForm.gamma1 / (rho6 + vdwForm.gamma) - 2.0;
1070 }
1071 case BUFFERED_14_7 -> {
1072 final double rho7 = rho3 * rho3 * rho;
1073 final double rhod7 = rhod3 * rhod3 * rhod;
1074 t1 = vdwForm.t1n / rhod7;
1075 t2 = vdwForm.gamma1 / (rho7 + vdwForm.gamma) - 2.0;
1076 }
1077 default -> {
1078
1079 final double rhoN = pow(rho, vdwForm.dispersivePower);
1080 final double rhodN = pow(rhod, vdwForm.dispersivePower);
1081 t1 = vdwForm.t1n / rhodN;
1082 t2 = vdwForm.gamma1 / (rhoN + vdwForm.gamma) - 2.0;
1083 }
1084 }
1085 final double eij = ev * t1 * t2;
1086
1087
1088
1089
1090
1091 double taper = 1.0;
1092 if (r2 < nonbondedCutoff.off2) {
1093 double r3 = r * r2;
1094 double r4 = r2 * r2;
1095 double r5 = r2 * r3;
1096 taper = multiplicativeSwitch.taper(r, r2, r3, r4, r5);
1097 taper = 1.0 - taper;
1098 }
1099
1100 double jacobian = 4.0 * PI * r2;
1101 double e = jacobian * eij * taper;
1102 if (j != i) {
1103 sume += e;
1104 } else {
1105 sume += 0.5 * e;
1106 }
1107 }
1108 double trapezoid = rDelta * sume;
1109
1110
1111 total += radCount[i] * radCount[j] * trapezoid;
1112
1113 if (lambda < 1.0) {
1114 total -= (softRadCount[i] * radCount[j] + (radCount[i] - softRadCount[i]) * softRadCount[j])
1115 * oneMinLambda * trapezoid;
1116 }
1117 }
1118 }
1119
1120
1121 Crystal unitCell = crystal.getUnitCell();
1122 double asymmetricUnitVolume = unitCell.volume / unitCell.getNumSymOps();
1123 total /= asymmetricUnitVolume;
1124 logger.fine(format(" Long-range vdW correction: %16.8f", total));
1125 return total;
1126 }
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138 private void log(int i, int k, double minr, double r, double eij) {
1139 Atom ai = atoms[i];
1140 Atom ak = atoms[k];
1141 logger.info(format("VDW %6d-%s %d %6d-%s %d %10.4f %10.4f %10.4f",
1142 ai.getIndex(), ai.getAtomType().name, ai.getVDWType().atomClass,
1143 ak.getIndex(), ak.getAtomType().name, ak.getVDWType().atomClass,
1144 minr, r, eij));
1145 logger.info(format(" %d %s\n %d %s",
1146 i + 1, Arrays.toString(ai.getXYZ(null)),
1147 k + 1, Arrays.toString(ak.getXYZ(null))));
1148 }
1149
1150 private void initSoftCore() {
1151
1152 boolean rebuild = false;
1153
1154 for (int i = 0; i < nAtoms; i++) {
1155 boolean soft = atoms[i].applyLambda();
1156 if (soft != isSoft[i]) {
1157 isSoft[i] = soft;
1158 rebuild = true;
1159 }
1160 }
1161
1162 if (!rebuild) {
1163 return;
1164 }
1165
1166
1167 for (int i = 0; i < nAtoms; i++) {
1168 if (isSoft[i]) {
1169
1170 softCore[HARD][i] = true;
1171
1172 softCore[SOFT][i] = false;
1173 } else {
1174
1175 softCore[HARD][i] = false;
1176
1177 softCore[SOFT][i] = true;
1178 }
1179 }
1180 }
1181
1182
1183
1184
1185
1186
1187
1188 public static class LambdaFactors {
1189
1190 double sc1 = 0.0;
1191 double dsc1dL = 0.0;
1192 double d2sc1dL2 = 0.0;
1193 double sc2 = 1.0;
1194 double dsc2dL = 0.0;
1195 double d2sc2dL2 = 0.0;
1196
1197
1198
1199
1200 public void setFactors() {
1201 }
1202 }
1203
1204 public class LambdaFactorsOST extends LambdaFactors {
1205
1206 @Override
1207 public void setFactors() {
1208 sc1 = VanDerWaals.this.sc1;
1209 dsc1dL = VanDerWaals.this.dsc1dL;
1210 d2sc1dL2 = VanDerWaals.this.d2sc1dL2;
1211 sc2 = VanDerWaals.this.sc2;
1212 dsc2dL = VanDerWaals.this.dsc2dL;
1213 d2sc2dL2 = VanDerWaals.this.d2sc2dL2;
1214 }
1215 }
1216
1217 private class VanDerWaalsRegion extends ParallelRegion {
1218
1219 private final InitializationLoop[] initializationLoop;
1220 private final ExpandLoop[] expandLoop;
1221 private final VanDerWaalsLoop[] vanDerWaalsLoop;
1222 private final ReductionLoop[] reductionLoop;
1223 private final NeighborListBarrier neighborListAction;
1224
1225
1226
1227
1228 private long initLoopTotalTime, vdWLoopTotalTime, reductionLoopTimeTotal;
1229 private long neighborListTotalTime, vdwTimeTotal;
1230
1231 VanDerWaalsRegion() {
1232 initializationLoop = new InitializationLoop[threadCount];
1233 expandLoop = new ExpandLoop[threadCount];
1234 vanDerWaalsLoop = new VanDerWaalsLoop[threadCount];
1235 reductionLoop = new ReductionLoop[threadCount];
1236 neighborListAction = new NeighborListBarrier();
1237 }
1238
1239 @Override
1240 public void finish() {
1241 forceNeighborListRebuild = false;
1242 vdwTimeTotal += System.nanoTime();
1243
1244 if (logger.isLoggable(Level.FINE)) {
1245 logger.fine(format("\n Neighbor-List: %7.4f (sec)", neighborListTotalTime * 1e-9));
1246 logger.fine(format(" Van der Waals: %7.4f (sec)", (vdwTimeTotal - neighborListTotalTime) * 1e-9));
1247 logger.fine(" Thread Init Energy Reduce Total Counts");
1248 long initMax = 0;
1249 long vdwMax = 0;
1250 long reductionMax = 0;
1251 long initMin = Long.MAX_VALUE;
1252 long vdwMin = Long.MAX_VALUE;
1253 long reductionMin = Long.MAX_VALUE;
1254 int countMin = Integer.MAX_VALUE;
1255 int countMax = 0;
1256 for (int i = 0; i < threadCount; i++) {
1257 int count = vanDerWaalsLoop[i].getCount();
1258 long initTime = initializationLoop[i].time + expandLoop[i].time;
1259 long vdWTime = vanDerWaalsLoop[i].time;
1260 long reduceTime = reductionLoop[i].time;
1261 long totalTime = initTime + vdWTime + reduceTime;
1262 logger.fine(format(" %3d %7.4f %7.4f %7.4f %7.4f %10d", i,
1263 initTime * 1e-9, vdWTime * 1e-9, reduceTime * 1e-9, totalTime * 1e-9, count));
1264 initMax = max(initTime, initMax);
1265 vdwMax = max(vdWTime, vdwMax);
1266 reductionMax = max(reduceTime, reductionMax);
1267 countMax = max(countMax, count);
1268 initMin = min(initTime, initMin);
1269 vdwMin = min(vdWTime, vdwMin);
1270 reductionMin = min(reduceTime, reductionMin);
1271 countMin = min(countMin, count);
1272 }
1273 long totalMin = initMin + vdwMin + reductionMin;
1274 long totalMax = initMax + vdwMax + reductionMax;
1275 long totalActual = initLoopTotalTime + vdWLoopTotalTime + reductionLoopTimeTotal;
1276 logger.fine(format(" Min %7.4f %7.4f %7.4f %7.4f %10d",
1277 initMin * 1e-9, vdwMin * 1e-9, reductionMin * 1e-9, totalMin * 1e-9, countMin));
1278 logger.fine(format(" Max %7.4f %7.4f %7.4f %7.4f %10d",
1279 initMax * 1e-9, vdwMax * 1e-9, reductionMax * 1e-9, totalMax * 1e-9, countMax));
1280 logger.fine(format(" Delta %7.4f %7.4f %7.4f %7.4f %10d",
1281 (initMax - initMin) * 1e-9, (vdwMax - vdwMin) * 1e-9,
1282 (reductionMax - reductionMin) * 1e-9, (totalMax - totalMin) * 1e-9,
1283 (countMax - countMin)));
1284 logger.fine(format(" Actual %7.4f %7.4f %7.4f %7.4f %10d\n",
1285 initLoopTotalTime * 1e-9, vdWLoopTotalTime * 1e-9, reductionLoopTimeTotal * 1e-9,
1286 totalActual * 1e-9, sharedInteractions.get()));
1287 }
1288 }
1289
1290 @Override
1291 public void run() throws Exception {
1292 int threadIndex = getThreadIndex();
1293
1294
1295 if (initializationLoop[threadIndex] == null) {
1296 initializationLoop[threadIndex] = new InitializationLoop();
1297 expandLoop[threadIndex] = new ExpandLoop();
1298 vanDerWaalsLoop[threadIndex] = new VanDerWaalsLoop();
1299 reductionLoop[threadIndex] = new ReductionLoop();
1300 }
1301
1302
1303 try {
1304 if (threadIndex == 0) {
1305 initLoopTotalTime = -System.nanoTime();
1306 }
1307 execute(0, nAtoms - 1, initializationLoop[threadIndex]);
1308 execute(0, nAtoms - 1, expandLoop[threadIndex]);
1309 if (threadIndex == 0) {
1310 initLoopTotalTime += System.nanoTime();
1311 }
1312 } catch (RuntimeException ex) {
1313 logger.warning("Runtime exception expanding coordinates in thread: " + threadIndex);
1314 throw ex;
1315 } catch (Exception e) {
1316 String message = "Fatal exception expanding coordinates in thread: " + threadIndex + "\n";
1317 logger.log(Level.SEVERE, message, e);
1318 }
1319
1320
1321 barrier(neighborListAction);
1322 if (forceNeighborListRebuild) {
1323 return;
1324 }
1325
1326
1327 try {
1328 if (threadIndex == 0) {
1329 vdWLoopTotalTime = -System.nanoTime();
1330 }
1331 execute(0, nAtoms - 1, vanDerWaalsLoop[threadIndex]);
1332 if (threadIndex == 0) {
1333 vdWLoopTotalTime += System.nanoTime();
1334 }
1335 } catch (RuntimeException ex) {
1336 logger.warning("Runtime exception evaluating van der Waals energy in thread: " + threadIndex);
1337 throw ex;
1338 } catch (Exception e) {
1339 String message = "Fatal exception evaluating van der Waals energy in thread: " + threadIndex + "\n";
1340 logger.log(Level.SEVERE, message, e);
1341 }
1342
1343
1344 if (gradient || lambdaTerm) {
1345 try {
1346 if (threadIndex == 0) {
1347 reductionLoopTimeTotal = -System.nanoTime();
1348 }
1349 execute(0, nAtoms - 1, reductionLoop[threadIndex]);
1350 if (threadIndex == 0) {
1351 reductionLoopTimeTotal += System.nanoTime();
1352 }
1353 } catch (RuntimeException ex) {
1354 logger.warning(
1355 "Runtime exception reducing van der Waals gradient in thread: " + threadIndex);
1356 throw ex;
1357 } catch (Exception e) {
1358 String message =
1359 "Fatal exception reducing van der Waals gradient in thread: " + threadIndex + "\n";
1360 logger.log(Level.SEVERE, message, e);
1361 }
1362 }
1363 }
1364
1365
1366
1367
1368
1369
1370 @Override
1371 public void start() {
1372
1373 vdwTimeTotal = -System.nanoTime();
1374
1375
1376 if (doLongRangeCorrection) {
1377 longRangeCorrection = computeLongRangeCorrection();
1378 sharedEnergy.set(longRangeCorrection);
1379 } else {
1380 sharedEnergy.set(0.0);
1381 }
1382 sharedInteractions.set(0);
1383 if (lambdaTerm) {
1384 shareddEdL.set(0.0);
1385 sharedd2EdL2.set(0.0);
1386 }
1387 if (esvTerm) {
1388 esvSystem.initEsvVdw();
1389 lambdaFactors = new LambdaFactors[threadCount];
1390 for (int i = 0; i < threadCount; i++) {
1391 lambdaFactors[i] = new LambdaFactors();
1392 }
1393 }
1394
1395 grad.alloc(nAtoms);
1396 if (lambdaTerm) {
1397 lambdaGrad.alloc(nAtoms);
1398 }
1399 }
1400
1401
1402
1403
1404 private class InitializationLoop extends IntegerForLoop {
1405
1406 private int threadID;
1407 protected long time;
1408
1409 @Override
1410 public void run(int lb, int ub) {
1411 for (int i = lb, i3 = 3 * lb; i <= ub; i++, i3 += 3) {
1412 Atom atom = atoms[i];
1413 coordinates[i3 + XX] = atom.getX();
1414 coordinates[i3 + YY] = atom.getY();
1415 coordinates[i3 + ZZ] = atom.getZ();
1416 use[i] = atom.getUse();
1417
1418
1419 VDWType vdwType = atom.getVDWType();
1420 if (vdwType == null) {
1421 logger.severe(" No VdW type for atom " + atom);
1422 return;
1423 }
1424
1425
1426 atomClass[i] = vdwType.atomClass;
1427
1428
1429 List<Bond> bonds = atom.getBonds();
1430 int numBonds = bonds.size();
1431 if (reducedHydrogen && vdwType.reductionFactor > 0.0 && numBonds == 1) {
1432 Bond bond = bonds.get(0);
1433 Atom heavyAtom = bond.get1_2(atom);
1434
1435 reductionIndex[i] = heavyAtom.getIndex() - 1;
1436 reductionValue[i] = vdwType.reductionFactor;
1437 } else {
1438 reductionIndex[i] = i;
1439 reductionValue[i] = 0.0;
1440 }
1441
1442
1443 List<Atom> n12 = atom.get12List();
1444 if (mask12[i] == null || mask12[i].length != n12.size()) {
1445 mask12[i] = new int[n12.size()];
1446 }
1447 int j = 0;
1448 for (Atom a12 : n12) {
1449 mask12[i][j++] = a12.getIndex() - 1;
1450 }
1451
1452
1453 List<Atom> n13 = atom.get13List();
1454 if (mask13[i] == null || mask13[i].length != n13.size()) {
1455 mask13[i] = new int[n13.size()];
1456 }
1457 j = 0;
1458 for (Atom a13 : n13) {
1459 mask13[i][j++] = a13.getIndex() - 1;
1460 }
1461
1462
1463 List<Atom> n14 = atom.get14List();
1464 if (mask14[i] == null || mask14[i].length != n14.size()) {
1465 mask14[i] = new int[n14.size()];
1466 }
1467 j = 0;
1468 for (Atom a14 : n14) {
1469 mask14[i][j++] = a14.getIndex() - 1;
1470 }
1471 }
1472 if (gradient) {
1473 grad.reset(threadID, lb, ub);
1474 }
1475 if (lambdaTerm) {
1476 lambdaGrad.reset(threadID, lb, ub);
1477 }
1478 }
1479
1480 @Override
1481 public IntegerSchedule schedule() {
1482 return IntegerSchedule.fixed();
1483 }
1484
1485 @Override
1486 public void finish() {
1487 time += System.nanoTime();
1488 }
1489
1490 @Override
1491 public void start() {
1492 threadID = getThreadIndex();
1493 time = -System.nanoTime();
1494 }
1495 }
1496
1497 private class ExpandLoop extends IntegerForLoop {
1498
1499 private final double[] in = new double[3];
1500 private final double[] out = new double[3];
1501 private long time;
1502
1503 @Override
1504 public void finish() {
1505 time += System.nanoTime();
1506 }
1507
1508 @Override
1509 public void run(int lb, int ub) {
1510
1511
1512
1513
1514 for (int i = lb; i <= ub; i++) {
1515 int i3 = i * 3;
1516 int iX = i3 + XX;
1517 int iY = i3 + YY;
1518 int iZ = i3 + ZZ;
1519 double x = coordinates[iX];
1520 double y = coordinates[iY];
1521 double z = coordinates[iZ];
1522 int redIndex = reductionIndex[i];
1523 if (redIndex >= 0) {
1524 int r3 = redIndex * 3;
1525 double rx = coordinates[r3++];
1526 double ry = coordinates[r3++];
1527 double rz = coordinates[r3];
1528 double a = reductionValue[i];
1529 reducedXYZ[iX] = a * (x - rx) + rx;
1530 reducedXYZ[iY] = a * (y - ry) + ry;
1531 reducedXYZ[iZ] = a * (z - rz) + rz;
1532 in[0] = reducedXYZ[iX];
1533 in[1] = reducedXYZ[iY];
1534 in[2] = reducedXYZ[iZ];
1535 atoms[i].setRedXYZ(in);
1536 } else {
1537 reducedXYZ[iX] = x;
1538 reducedXYZ[iY] = y;
1539 reducedXYZ[iZ] = z;
1540 }
1541 }
1542
1543 List<SymOp> symOps = crystal.spaceGroup.symOps;
1544
1545 if (symOps.size() != nSymm) {
1546 String message = format(" Programming Error: nSymm %d != symOps.size %d", nSymm, symOps.size());
1547 logger.log(Level.WARNING, message);
1548 logger.log(Level.WARNING, " Replicates\n{0}", crystal.toString());
1549 logger.log(Level.WARNING, " Unit Cell\n{0}", crystal.getUnitCell().toString());
1550 }
1551
1552 double sp2 = crystal.getSpecialPositionCutoff();
1553 sp2 *= sp2;
1554 for (int iSymOp = 1; iSymOp < nSymm; iSymOp++) {
1555 SymOp symOp = symOps.get(iSymOp);
1556 double[] xyz = reduced[iSymOp];
1557 for (int i = lb; i <= ub; i++) {
1558 int i3 = i * 3;
1559 int iX = i3 + XX;
1560 int iY = i3 + YY;
1561 int iZ = i3 + ZZ;
1562 in[0] = reducedXYZ[iX];
1563 in[1] = reducedXYZ[iY];
1564 in[2] = reducedXYZ[iZ];
1565 crystal.applySymOp(in, out, symOp);
1566 xyz[iX] = out[0];
1567 xyz[iY] = out[1];
1568 xyz[iZ] = out[2];
1569
1570
1571 double dx = in[0] - out[0];
1572 double dy = in[1] - out[1];
1573 double dz = in[2] - out[2];
1574 double r2 = dx * dx + dy * dy + dz * dz;
1575 if (r2 < sp2 && logger.isLoggable(Level.FINEST)) {
1576 logger.log(Level.FINEST, " Atom may be at a special position: {0}", atoms[i].toString());
1577 }
1578 }
1579 }
1580 }
1581
1582 @Override
1583 public IntegerSchedule schedule() {
1584 return IntegerSchedule.fixed();
1585 }
1586
1587 @Override
1588 public void start() {
1589 time = -System.nanoTime();
1590 }
1591
1592 }
1593
1594
1595
1596
1597
1598
1599
1600
1601 private class VanDerWaalsLoop extends IntegerForLoop {
1602 protected long time;
1603 private final double[] dx_local;
1604 private final double[][] transOp;
1605 private int count;
1606 private double energy;
1607 private int threadID;
1608 private double dEdL;
1609 private double d2EdL2;
1610 private double[] mask;
1611 private boolean[] vdw14;
1612 private LambdaFactors lambdaFactorsLocal;
1613
1614 VanDerWaalsLoop() {
1615 super();
1616 dx_local = new double[3];
1617 transOp = new double[3][3];
1618 }
1619
1620 @Override
1621 public void finish() {
1622
1623 sharedEnergy.addAndGet(energy);
1624 sharedInteractions.addAndGet(count);
1625 if (lambdaTerm) {
1626 shareddEdL.addAndGet(dEdL);
1627 sharedd2EdL2.addAndGet(d2EdL2);
1628 }
1629 time += System.nanoTime();
1630 }
1631
1632 public int getCount() {
1633 return count;
1634 }
1635
1636 @Override
1637 public void run(int lb, int ub) {
1638 double e = 0.0;
1639 double[] xyzS = reduced[0];
1640
1641 int[][] list = neighborLists[0];
1642 double[] esvVdwPrefactori = new double[3];
1643 double[] esvVdwPrefactork = new double[3];
1644 for (int i = lb; i <= ub; i++) {
1645 if (!use[i]) {
1646 continue;
1647 }
1648
1649 final boolean esvi = esvTerm && esvSystem.isTitratingHydrogen(i);
1650 if (esvTerm) {
1651 esvSystem.getVdwPrefactor(i, esvVdwPrefactori);
1652 }
1653 int i3 = i * 3;
1654
1655 final double xi = reducedXYZ[i3++];
1656 final double yi = reducedXYZ[i3++];
1657 final double zi = reducedXYZ[i3];
1658
1659 final int redi = reductionIndex[i];
1660
1661 final double redv = reductionValue[i];
1662 final double rediv = 1.0 - redv;
1663 final int classI = atomClass[i];
1664
1665 final int moleculei = molecule[i];
1666
1667 boolean iNN = neuralNetwork[i];
1668
1669 double gxi = 0.0;
1670 double gyi = 0.0;
1671 double gzi = 0.0;
1672
1673 double gxredi = 0.0;
1674 double gyredi = 0.0;
1675 double gzredi = 0.0;
1676
1677 double lxi = 0.0;
1678 double lyi = 0.0;
1679 double lzi = 0.0;
1680
1681 double lxredi = 0.0;
1682 double lyredi = 0.0;
1683 double lzredi = 0.0;
1684
1685
1686 applyMask(i, vdw14, mask);
1687
1688 boolean[] softCorei = softCore[HARD];
1689 if (isSoft[i]) {
1690 softCorei = softCore[SOFT];
1691 }
1692
1693 final int[] neighbors = list[i];
1694 for (final int k : neighbors) {
1695
1696
1697 if (!use[k] || (iNN && neuralNetwork[k])) {
1698 continue;
1699 }
1700
1701 final boolean esvk = esvTerm && esvSystem.isTitratingHydrogen(k);
1702 if (esvTerm) {
1703 esvSystem.getVdwPrefactor(k, esvVdwPrefactork);
1704 }
1705 int k3 = k * 3;
1706
1707 final double xk = xyzS[k3++];
1708 final double yk = xyzS[k3++];
1709 final double zk = xyzS[k3];
1710
1711 dx_local[0] = xi - xk;
1712 dx_local[1] = yi - yk;
1713 dx_local[2] = zi - zk;
1714
1715 final double r2 = crystal.image(dx_local);
1716 int classK = atomClass[k];
1717 double irv = vdwForm.getCombinedInverseRmin(classI, classK);
1718 if (vdw14[k]) {
1719
1720 irv = vdwForm.getCombinedInverseRmin14(classI, classK);
1721 }
1722 if (r2 <= nonbondedCutoff.off2 && mask[k] > 0 && irv > 0) {
1723
1724 final double r = sqrt(r2);
1725
1726 boolean sameMolecule = (moleculei == molecule[k]);
1727 boolean soft = softCorei[k];
1728
1729
1730 if (isSoft[i] && isSoft[k]) {
1731 if (intermolecularSoftcore && !sameMolecule) {
1732 soft = true;
1733 } else if (intramolecularSoftcore && sameMolecule) {
1734 soft = true;
1735 }
1736 }
1737
1738 final double sc1, dsc1dL, d2sc1dL2;
1739 final double sc2, dsc2dL, d2sc2dL2;
1740 if (soft) {
1741
1742
1743
1744
1745 sc1 = lambdaFactorsLocal.sc1;
1746 dsc1dL = lambdaFactorsLocal.dsc1dL;
1747 d2sc1dL2 = lambdaFactorsLocal.d2sc1dL2;
1748 sc2 = lambdaFactorsLocal.sc2;
1749 dsc2dL = lambdaFactorsLocal.dsc2dL;
1750 d2sc2dL2 = lambdaFactorsLocal.d2sc2dL2;
1751 } else {
1752 sc1 = 0.0;
1753 dsc1dL = 0.0;
1754 d2sc1dL2 = 0.0;
1755 sc2 = 1.0;
1756 dsc2dL = 0.0;
1757 d2sc2dL2 = 0.0;
1758 }
1759 final double alpha = sc1;
1760 final double lambda5 = sc2;
1761
1762
1763
1764
1765
1766
1767 double ev = mask[k] * vdwForm.getCombinedEps(classI, classK);
1768 if (vdw14[k]) {
1769 ev = mask[k] * vdwForm.getCombinedEps14(classI, classK);
1770 }
1771 final double eps_lambda = ev * lambda5;
1772 final double rho = r * irv;
1773 final double rhoDisp1 = vdwForm.rhoDisp1(rho);
1774 final double rhoDisp = rhoDisp1 * rho;
1775 final double rhoDelta1 = vdwForm.rhoDelta1(rho + vdwForm.delta);
1776 final double rhoDelta = rhoDelta1 * (rho + vdwForm.delta);
1777 final double alphaRhoDelta = alpha + rhoDelta;
1778 final double alphaRhoDispGamma = alpha + rhoDisp + vdwForm.gamma;
1779 final double t1d = 1.0 / alphaRhoDelta;
1780 final double t2d = 1.0 / alphaRhoDispGamma;
1781 final double t1 = vdwForm.t1n * t1d;
1782 final double t2a = vdwForm.gamma1 * t2d;
1783 final double t2 = t2a - 2.0;
1784 double eik = eps_lambda * t1 * t2;
1785
1786
1787
1788
1789 double taper = 1.0;
1790 double dtaper = 0.0;
1791 if (r2 > nonbondedCutoff.cut2) {
1792 final double r3 = r2 * r;
1793 final double r4 = r2 * r2;
1794 final double r5 = r2 * r3;
1795 taper = multiplicativeSwitch.taper(r, r2, r3, r4, r5);
1796 dtaper = multiplicativeSwitch.dtaper(r, r2, r3, r4);
1797 }
1798 double esvik = 1.0;
1799 if (esvi || esvk) {
1800 esvik = esvVdwPrefactori[0] * esvVdwPrefactork[0];
1801 }
1802 e += eik * taper * esvik;
1803
1804 count++;
1805 if (!gradient && !soft) {
1806 continue;
1807 }
1808 final int redk = reductionIndex[k];
1809 final double red = reductionValue[k];
1810 final double redkv = 1.0 - red;
1811 final double dt1d_dr = vdwForm.repDispPower * rhoDelta1 * irv;
1812 final double dt2d_dr = vdwForm.dispersivePower * rhoDisp1 * irv;
1813 final double dt1_dr = t1 * dt1d_dr * t1d;
1814 final double dt2_dr = t2a * dt2d_dr * t2d;
1815 double dedr = -eps_lambda * (dt1_dr * t2 + t1 * dt2_dr);
1816 final double ir = 1.0 / r;
1817 final double drdx = dx_local[0] * ir;
1818 final double drdy = dx_local[1] * ir;
1819 final double drdz = dx_local[2] * ir;
1820 if (gradient) {
1821 final double dswitch = esvik * (eik * dtaper + dedr * taper);
1822 final double dedx = dswitch * drdx;
1823 final double dedy = dswitch * drdy;
1824 final double dedz = dswitch * drdz;
1825 gxi += dedx * redv;
1826 gyi += dedy * redv;
1827 gzi += dedz * redv;
1828 gxredi += dedx * rediv;
1829 gyredi += dedy * rediv;
1830 gzredi += dedz * rediv;
1831 grad.sub(threadID, k, red * dedx, red * dedy, red * dedz);
1832 grad.sub(threadID, redk, redkv * dedx, redkv * dedy, redkv * dedz);
1833
1834 if (esvi) {
1835 esvSystem.addVdwDeriv(i, eik * taper, esvVdwPrefactori, esvVdwPrefactork[0]);
1836 }
1837 if (esvk) {
1838 esvSystem.addVdwDeriv(k, eik * taper, esvVdwPrefactork, esvVdwPrefactori[0]);
1839 }
1840 }
1841 if (gradient && soft) {
1842 final double dt1 = -t1 * t1d * dsc1dL;
1843 final double dt2 = -t2a * t2d * dsc1dL;
1844 final double f1 = dsc2dL * t1 * t2;
1845 final double f2 = sc2 * dt1 * t2;
1846 final double f3 = sc2 * t1 * dt2;
1847 final double dedl = ev * (f1 + f2 + f3);
1848 dEdL += dedl * taper * esvik;
1849 final double t1d2 = -dsc1dL * t1d * t1d;
1850 final double t2d2 = -dsc1dL * t2d * t2d;
1851 final double d2t1 = -dt1 * t1d * dsc1dL - t1 * t1d * d2sc1dL2 - t1 * t1d2 * dsc1dL;
1852 final double d2t2 =
1853 -dt2 * t2d * dsc1dL - t2a * t2d * d2sc1dL2 - t2a * t2d2 * dsc1dL;
1854 final double df1 = d2sc2dL2 * t1 * t2 + dsc2dL * dt1 * t2 + dsc2dL * t1 * dt2;
1855 final double df2 = dsc2dL * dt1 * t2 + sc2 * d2t1 * t2 + sc2 * dt1 * dt2;
1856 final double df3 = dsc2dL * t1 * dt2 + sc2 * dt1 * dt2 + sc2 * t1 * d2t2;
1857 final double de2dl2 = ev * (df1 + df2 + df3);
1858 d2EdL2 += de2dl2 * taper * esvik;
1859 final double t11 = -dsc2dL * t2 * dt1_dr;
1860 final double t21 = -dsc2dL * t1 * dt2_dr;
1861 final double t13 = 2.0 * sc2 * t2 * dt1_dr * dsc1dL * t1d;
1862 final double t23 = 2.0 * sc2 * t1 * dt2_dr * dsc1dL * t2d;
1863 final double t12 = -sc2 * dt2 * dt1_dr;
1864 final double t22 = -sc2 * dt1 * dt2_dr;
1865 final double dedldr = ev * (t11 + t12 + t13 + t21 + t22 + t23);
1866 final double dswitch = esvik * (dedl * dtaper + dedldr * taper);
1867 final double dedldx = dswitch * drdx;
1868 final double dedldy = dswitch * drdy;
1869 final double dedldz = dswitch * drdz;
1870 lxi += dedldx * redv;
1871 lyi += dedldy * redv;
1872 lzi += dedldz * redv;
1873 lxredi += dedldx * rediv;
1874 lyredi += dedldy * rediv;
1875 lzredi += dedldz * rediv;
1876 if (lambdaTerm) {
1877 lambdaGrad.sub(threadID, k, red * dedldx, red * dedldy, red * dedldz);
1878 lambdaGrad.sub(threadID, redk, redkv * dedldx, redkv * dedldy, redkv * dedldz);
1879 }
1880 }
1881 }
1882 }
1883 if (gradient) {
1884 grad.add(threadID, i, gxi, gyi, gzi);
1885 grad.add(threadID, redi, gxredi, gyredi, gzredi);
1886 if (lambdaTerm) {
1887 lambdaGrad.add(threadID, i, lxi, lyi, lzi);
1888 lambdaGrad.add(threadID, redi, lxredi, lyredi, lzredi);
1889 }
1890 }
1891 removeMask(i, vdw14, mask);
1892 }
1893 energy += e;
1894 List<SymOp> symOps = crystal.spaceGroup.symOps;
1895 for (int iSymOp = 1; iSymOp < nSymm; iSymOp++) {
1896 e = 0.0;
1897 SymOp symOp = symOps.get(iSymOp);
1898
1899 crystal.getTransformationOperator(symOp, transOp);
1900 xyzS = reduced[iSymOp];
1901 list = neighborLists[iSymOp];
1902
1903 for (int i = lb; i <= ub; i++) {
1904 int i3 = i * 3;
1905 if (!use[i]) {
1906 continue;
1907 }
1908 final boolean esvi = esvTerm && esvSystem.isTitratingHydrogen(i);
1909 if (esvTerm) {
1910 esvSystem.getVdwPrefactor(i, esvVdwPrefactori);
1911 }
1912 final double xi = reducedXYZ[i3++];
1913 final double yi = reducedXYZ[i3++];
1914 final double zi = reducedXYZ[i3];
1915 final int redi = reductionIndex[i];
1916 final double redv = reductionValue[i];
1917 final double rediv = 1.0 - redv;
1918 final int classI = atomClass[i];
1919 double gxi = 0.0;
1920 double gyi = 0.0;
1921 double gzi = 0.0;
1922 double gxredi = 0.0;
1923 double gyredi = 0.0;
1924 double gzredi = 0.0;
1925 double lxi = 0.0;
1926 double lyi = 0.0;
1927 double lzi = 0.0;
1928 double lxredi = 0.0;
1929 double lyredi = 0.0;
1930 double lzredi = 0.0;
1931
1932 boolean[] softCorei = softCore[HARD];
1933 if (isSoft[i]) {
1934 softCorei = softCore[SOFT];
1935 }
1936
1937 final int[] neighbors = list[i];
1938 for (final int k : neighbors) {
1939 if (!use[k]) {
1940 continue;
1941 }
1942 final boolean esvk = esvTerm && esvSystem.isTitratingHydrogen(k);
1943 if (esvTerm) {
1944 esvSystem.getVdwPrefactor(k, esvVdwPrefactork);
1945 }
1946
1947 final double sc1, dsc1dL, d2sc1dL2;
1948 final double sc2, dsc2dL, d2sc2dL2;
1949 int k3 = k * 3;
1950 final double xk = xyzS[k3++];
1951 final double yk = xyzS[k3++];
1952 final double zk = xyzS[k3];
1953 dx_local[0] = xi - xk;
1954 dx_local[1] = yi - yk;
1955 dx_local[2] = zi - zk;
1956 final double r2 = crystal.image(dx_local);
1957 int classK = atomClass[k];
1958 final double irv = vdwForm.getCombinedInverseRmin(classI, classK);
1959 if (r2 <= nonbondedCutoff.off2 && irv > 0) {
1960 final double selfScale = (i == k) ? 0.5 : 1.0;
1961 final double r = sqrt(r2);
1962 boolean soft = isSoft[i] || softCorei[k];
1963 if (soft) {
1964 sc1 = lambdaFactorsLocal.sc1;
1965 dsc1dL = lambdaFactorsLocal.dsc1dL;
1966 d2sc1dL2 = lambdaFactorsLocal.d2sc1dL2;
1967 sc2 = lambdaFactorsLocal.sc2;
1968 dsc2dL = lambdaFactorsLocal.dsc2dL;
1969 d2sc2dL2 = lambdaFactorsLocal.d2sc2dL2;
1970 } else {
1971 sc1 = 0.0;
1972 dsc1dL = 0.0;
1973 d2sc1dL2 = 0.0;
1974 sc2 = 1.0;
1975 dsc2dL = 0.0;
1976 d2sc2dL2 = 0.0;
1977 }
1978 final double alpha = sc1;
1979 final double lambdaN = sc2;
1980 final double ev = vdwForm.getCombinedEps(classI, classK);
1981 final double eps_lambda = ev * lambdaN;
1982 final double rho = r * irv;
1983 final double rhoDisp1 = vdwForm.rhoDisp1(rho);
1984 final double rhoDisp = rhoDisp1 * rho;
1985 final double rhoDelta1 = vdwForm.rhoDelta1(rho + vdwForm.delta);
1986 final double rhoDelta = rhoDelta1 * (rho + vdwForm.delta);
1987 final double alphaRhoDelta = alpha + rhoDelta;
1988 final double alphaRhoDispGamma = alpha + rhoDisp + vdwForm.gamma;
1989 final double t1d = 1.0 / alphaRhoDelta;
1990 final double t2d = 1.0 / alphaRhoDispGamma;
1991 final double t1 = vdwForm.t1n * t1d;
1992 final double t2a = vdwForm.gamma1 * t2d;
1993 final double t2 = t2a - 2.0;
1994 double eik = eps_lambda * t1 * t2;
1995
1996
1997 double taper = 1.0;
1998 double dtaper = 0.0;
1999 if (r2 > nonbondedCutoff.cut2) {
2000 final double r3 = r2 * r;
2001 final double r4 = r2 * r2;
2002 final double r5 = r2 * r3;
2003 taper = multiplicativeSwitch.taper(r, r2, r3, r4, r5);
2004 dtaper = multiplicativeSwitch.dtaper(r, r2, r3, r4);
2005 }
2006
2007 double esvik = 1.0;
2008 if (esvi || esvk) {
2009 esvik = esvVdwPrefactori[0] * esvVdwPrefactork[0];
2010 }
2011 e += selfScale * eik * taper * esvik;
2012 count++;
2013 if (!gradient && !soft) {
2014 continue;
2015 }
2016 final int redk = reductionIndex[k];
2017 final double red = reductionValue[k];
2018 final double redkv = 1.0 - red;
2019 final double dt1d_dr = vdwForm.repDispPower * rhoDelta1 * irv;
2020 final double dt2d_dr = vdwForm.dispersivePower * rhoDisp1 * irv;
2021 final double dt1_dr = t1 * dt1d_dr * t1d;
2022 final double dt2_dr = t2a * dt2d_dr * t2d;
2023 double dedr = -eps_lambda * (dt1_dr * t2 + t1 * dt2_dr);
2024
2025 final double ir = 1.0 / r;
2026 double drdx = dx_local[0] * ir;
2027 double drdy = dx_local[1] * ir;
2028 double drdz = dx_local[2] * ir;
2029 dedr = selfScale * esvik * (eik * dtaper + dedr * taper);
2030 if (gradient) {
2031 double dedx = dedr * drdx;
2032 double dedy = dedr * drdy;
2033 double dedz = dedr * drdz;
2034 gxi += dedx * redv;
2035 gyi += dedy * redv;
2036 gzi += dedz * redv;
2037 gxredi += dedx * rediv;
2038 gyredi += dedy * rediv;
2039 gzredi += dedz * rediv;
2040
2041
2042 final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
2043 final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
2044 final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];
2045 grad.sub(threadID, k, red * dedxk, red * dedyk, red * dedzk);
2046 grad.sub(threadID, redk, redkv * dedxk, redkv * dedyk, redkv * dedzk);
2047
2048
2049 if (esvi) {
2050 esvSystem.addVdwDeriv(i, selfScale * eik * taper, esvVdwPrefactori,
2051 esvVdwPrefactork[0]);
2052 }
2053 if (esvk) {
2054 esvSystem.addVdwDeriv(k, selfScale * eik * taper, esvVdwPrefactork,
2055 esvVdwPrefactori[0]);
2056 }
2057 }
2058 if (gradient && lambdaTerm && soft) {
2059 double dt1 = -t1 * t1d * dsc1dL;
2060 double dt2 = -t2a * t2d * dsc1dL;
2061 double f1 = dsc2dL * t1 * t2;
2062 double f2 = sc2 * dt1 * t2;
2063 double f3 = sc2 * t1 * dt2;
2064 final double dedl = ev * (f1 + f2 + f3);
2065 dEdL += selfScale * dedl * taper * esvik;
2066 double t1d2 = -dsc1dL * t1d * t1d;
2067 double t2d2 = -dsc1dL * t2d * t2d;
2068 double d2t1 = -dt1 * t1d * dsc1dL - t1 * t1d * d2sc1dL2 - t1 * t1d2 * dsc1dL;
2069 double d2t2 = -dt2 * t2d * dsc1dL - t2a * t2d * d2sc1dL2 - t2a * t2d2 * dsc1dL;
2070 double df1 = d2sc2dL2 * t1 * t2 + dsc2dL * dt1 * t2 + dsc2dL * t1 * dt2;
2071 double df2 = dsc2dL * dt1 * t2 + sc2 * d2t1 * t2 + sc2 * dt1 * dt2;
2072 double df3 = dsc2dL * t1 * dt2 + sc2 * dt1 * dt2 + sc2 * t1 * d2t2;
2073 double de2dl2 = ev * (df1 + df2 + df3);
2074 d2EdL2 += selfScale * de2dl2 * taper * esvik;
2075 double t11 = -dsc2dL * t2 * dt1_dr;
2076 double t12 = -sc2 * dt2 * dt1_dr;
2077 double t13 = 2.0 * sc2 * t2 * dt1_dr * dsc1dL * t1d;
2078 double t21 = -dsc2dL * t1 * dt2_dr;
2079 double t22 = -sc2 * dt1 * dt2_dr;
2080 double t23 = 2.0 * sc2 * t1 * dt2_dr * dsc1dL * t2d;
2081 double dedldr = ev * (t11 + t12 + t13 + t21 + t22 + t23);
2082 dedldr = esvik * (dedl * dtaper + dedldr * taper);
2083 double dedldx = selfScale * dedldr * drdx;
2084 double dedldy = selfScale * dedldr * drdy;
2085 double dedldz = selfScale * dedldr * drdz;
2086 lxi += dedldx * redv;
2087 lyi += dedldy * redv;
2088 lzi += dedldz * redv;
2089 lxredi += dedldx * rediv;
2090 lyredi += dedldy * rediv;
2091 lzredi += dedldz * rediv;
2092
2093
2094 final double dedldxk =
2095 dedldx * transOp[0][0] + dedldy * transOp[1][0] + dedldz * transOp[2][0];
2096 final double dedldyk =
2097 dedldx * transOp[0][1] + dedldy * transOp[1][1] + dedldz * transOp[2][1];
2098 final double dedldzk =
2099 dedldx * transOp[0][2] + dedldy * transOp[1][2] + dedldz * transOp[2][2];
2100 lambdaGrad.sub(threadID, k, red * dedldxk, red * dedldyk, red * dedldzk);
2101 lambdaGrad.sub(threadID, redk, redkv * dedldxk, redkv * dedldyk, redkv * dedldzk);
2102 }
2103 }
2104 }
2105 if (gradient) {
2106 grad.add(threadID, i, gxi, gyi, gzi);
2107 grad.add(threadID, redi, gxredi, gyredi, gzredi);
2108 if (lambdaTerm) {
2109 lambdaGrad.add(threadID, i, lxi, lyi, lzi);
2110 lambdaGrad.add(threadID, redi, lxredi, lyredi, lzredi);
2111 }
2112 }
2113 }
2114 energy += e;
2115 }
2116 }
2117
2118 @Override
2119 public IntegerSchedule schedule() {
2120 return pairwiseSchedule;
2121 }
2122
2123 @Override
2124 public void start() {
2125 threadID = getThreadIndex();
2126 time = -System.nanoTime();
2127 energy = 0.0;
2128 count = 0;
2129 if (lambdaTerm) {
2130 dEdL = 0.0;
2131 d2EdL2 = 0.0;
2132 }
2133 lambdaFactorsLocal = lambdaFactors[threadID];
2134 if (lambdaFactorsLocal == null) {
2135 System.exit(1);
2136 }
2137
2138 if (mask == null || mask.length < nAtoms) {
2139 mask = new double[nAtoms];
2140 fill(mask, 1.0);
2141 vdw14 = new boolean[nAtoms];
2142 fill(vdw14, false);
2143 }
2144 }
2145
2146 }
2147
2148
2149
2150
2151 private class ReductionLoop extends IntegerForLoop {
2152
2153 protected long time;
2154 int threadID;
2155
2156 @Override
2157 public void finish() {
2158 time += System.nanoTime();
2159 }
2160
2161 @Override
2162 public void run(int lb, int ub) {
2163 if (gradient) {
2164 grad.reduce(lb, ub);
2165 for (int i = lb; i <= ub; i++) {
2166 Atom ai = atoms[i];
2167 ai.addToXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
2168 }
2169 }
2170 if (lambdaTerm) {
2171 lambdaGrad.reduce(lb, ub);
2172 }
2173 }
2174
2175 @Override
2176 public void start() {
2177 threadID = getThreadIndex();
2178 time = -System.nanoTime();
2179 }
2180
2181 }
2182
2183
2184
2185
2186 private class NeighborListBarrier extends BarrierAction {
2187
2188 @Override
2189 public void run() throws Exception {
2190 neighborListTotalTime = -System.nanoTime();
2191 neighborList.buildList(reduced, neighborLists, null, forceNeighborListRebuild, false);
2192 neighborListTotalTime += System.nanoTime();
2193 }
2194 }
2195 }
2196 }