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