View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * The Van der Waals class computes Van der Waals interaction in parallel using a
80   * {@link ffx.potential.nonbonded.NeighborList} for any {@link ffx.crystal.Crystal}. The repulsive
81   * power (e.g. 12), attractive power (e.g. 6) and buffering (e.g. for the AMOEBA buffered-14-7) can
82   * all be specified such that both Lennard-Jones and AMOEBA are supported.
83   *
84   * @author Michael J. Schnieders
85   * @since 1.0
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     * The ParallelTeam used to evaluate van der Waals interactions.
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    * Boundary conditions and crystal symmetry.
127    */
128   private Crystal crystal;
129   /**
130    * An array of all atoms in the system.
131    */
132   private Atom[] atoms;
133   /**
134    * Specification of the molecular index for each atom.
135    */
136   private int[] molecule;
137   /**
138    * Flag to indicate the atom is treated by a neural network.
139    */
140   private boolean[] neuralNetwork;
141   /**
142    * The Force Field that defines the Van der Waals interactions.
143    */
144   private ForceField forceField;
145   /**
146    * An array of whether each atom in the system should be used in the calculations.
147    */
148   private boolean[] use = null;
149   /**
150    * A local convenience variable equal to atoms.length.
151    */
152   private int nAtoms;
153   /**
154    * A local convenience variable equal to the number of crystal symmetry operators.
155    */
156   private int nSymm;
157   /**
158    * **************************************************************** Lambda variables.
159    */
160   private boolean gradient;
161   private boolean lambdaTerm;
162   private boolean esvTerm;
163   private boolean[] isSoft;
164   /**
165    * There are 2 softCore arrays of length nAtoms.
166    *
167    * <p>The first is used for atoms in the outer loop that are hard. This mask equals: false for
168    * inner loop hard atoms true for inner loop soft atoms
169    *
170    * <p>The second is used for atoms in the outer loop that are soft. This mask equals: true for
171    * inner loop hard atoms false for inner loop soft atoms
172    */
173   private boolean[][] softCore;
174   /**
175    * Turn on inter-molecular softcore interactions using molecular index.
176    */
177   private boolean intermolecularSoftcore = false;
178 
179   // *************************************************************************
180   // Coordinate arrays.
181   /**
182    * Turn on intra-molecular softcore interactions using molecular index.
183    */
184   private boolean intramolecularSoftcore = false;
185   /**
186    * Current value of the lambda state variable.
187    */
188   private double lambda = 1.0;
189   /**
190    * The vdW energy is fully turned on when this value is reached.
191    */
192   private double vdwLambdaEnd = 1.0;
193   /**
194    * Exponent on lambda (beta).
195    */
196   private double vdwLambdaExponent = 3.0;
197   /**
198    * Offset in Angstroms (alpha).
199    */
200   private double vdwLambdaAlpha = 0.25;
201 
202   /**
203    * Polymorphic inner class to set sc1,sc2,dsc1,etc only when necessary. [nThreads]
204    */
205   private LambdaFactors[] lambdaFactors = null;
206   /**
207    * alpha * (1 - lambda)^2
208    */
209   private double sc1 = 0.0;
210   /**
211    * lambda^lambdaExponent
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    * Generalized extended system variables.
220    */
221   private ExtendedSystem esvSystem;
222 
223   /**
224    * A local copy of atomic coordinates, including reductions on the hydrogen atoms.
225    */
226   private double[] coordinates;
227   /**
228    * Reduced coordinates of size: [nSymm][nAtoms * 3]
229    */
230   private double[][] reduced;
231 
232   private double[] reducedXYZ;
233   /**
234    * Neighbor lists for each atom. Size: [nSymm][nAtoms][nNeighbors]
235    */
236   private int[][][] neighborLists;
237   /**
238    * A local reference to the atom class of each atom in the system.
239    */
240   private int[] atomClass;
241   /**
242    * Hydrogen atom vdW sites are located toward their heavy atom relative to their nucleus. This is a
243    * look-up that gives the heavy atom index for each hydrogen.
244    */
245   private int[] reductionIndex;
246 
247   private int[][] mask12;
248   private int[][] mask13;
249   private int[][] mask14;
250   /**
251    * Each hydrogen vdW site is located a fraction of the way from the heavy atom nucleus to the
252    * hydrogen nucleus (~0.9).
253    */
254   private double[] reductionValue;
255 
256   private boolean reducedHydrogen;
257   private double longRangeCorrection;
258   private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
259   /**
260    * Cartesian coordinate gradient.
261    */
262   private AtomicDoubleArray3D grad;
263   /**
264    * Lambda derivative of the Cartesian coordinate gradient.
265    */
266   private AtomicDoubleArray3D lambdaGrad;
267   /**
268    * The neighbor-list includes 1-2 and 1-3 interactions, which are masked out in the Van der Waals
269    * energy code. The AMOEBA force field includes 1-4 interactions fully.
270    */
271   private NeighborList neighborList;
272   /**
273    * Force building of the neighbor list.
274    */
275   private boolean forceNeighborListRebuild = true;
276 
277   public VanDerWaals() {
278     // Empty constructor for use with VanDerWaalsTornado
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    * The VanDerWaals class constructor.
297    *
298    * @param atoms              the Atom array to do Van Der Waals calculations on.
299    * @param molecule           the molecule number for each atom.
300    * @param neuralNetwork      an array of flags to indicate the atom is treated by a neural network.
301    * @param crystal            The boundary conditions.
302    * @param forceField         the ForceField parameters to apply.
303    * @param parallelTeam       The parallel environment.
304    * @param vdwCutoff          a double.
305    * @param neighborListCutoff a double.
306    * @since 1.0
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     // Lambda parameters.
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     // Parallel constructs.
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     // Define how force arrays will be accumulated.
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     // Allocate coordinate arrays and set up reduction indices and values.
383     initAtomArrays();
384 
385     /*
386      Define the multiplicative switch, which sets vdW energy to zero
387      at the cutoff distance using a window that begin at 90% of the
388      cutoff distance.
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     // Reduce and expand the coordinates of the asymmetric unit. Then build the first neighbor-list.
401     buildNeighborList(atoms);
402 
403     // Then, optionally, prevent that neighbor list from ever updating.
404     neighborList.setDisableUpdates(forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
405 
406     logger.info(toString());
407   }
408 
409   /**
410    * {@inheritDoc}
411    *
412    * <p>Apply masking rules for 1-2, 1-3 and 1-4 interactions.
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    * attachExtendedSystem.
434    *
435    * @param system a {@link ffx.potential.extended.ExtendedSystem} object.
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     // Launch shared lambda/esvLambda initializers if missed (ie. !lambdaTerm) in constructor.
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     // The combination of an extended system and neural networks are not supported.
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    * Get the ExtendedSystem instance.
480    *
481    * @return The ExtendedSystem is returned.
482    */
483   public ExtendedSystem getExtendedSystem() {
484     return esvSystem;
485   }
486 
487   /**
488    * destroy.
489    *
490    * @throws java.lang.Exception if any.
491    */
492   public void destroy() throws Exception {
493     if (neighborList != null) {
494       neighborList.destroy();
495     }
496   }
497 
498   /**
499    * The energy routine may be called repeatedly.
500    *
501    * @param gradient If true, gradients with respect to atomic coordinates are computed.
502    * @param print    If true, there is verbose printing.
503    * @return The energy.
504    * @since 1.0
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    * If true, there are alchemical atoms impacted by the lambda state variable.
520    * @return True if there are alchemical atoms.
521    */
522   public boolean getLambdaTerm() {
523     return lambdaTerm;
524   }
525 
526   /**
527    * If true, intra-molecular interactions are annihilated by the lambda state variable.
528    */
529   public boolean getIntramolecularSoftcore() {
530     return intramolecularSoftcore;
531   }
532 
533   /**
534    * getAlpha.
535    *
536    * @return a double.
537    */
538   public double getAlpha() {
539     return vdwLambdaAlpha;
540   }
541 
542   /**
543    * getBeta.
544    *
545    * @return a double.
546    */
547   public double getBeta() {
548     return vdwLambdaExponent;
549   }
550 
551   /**
552    * Get the buffer size.
553    *
554    * @return The buffer.
555    * @since 1.0
556    */
557   public double getBuffer() {
558     return nonbondedCutoff.buff;
559   }
560 
561   /**
562    * Return use of the long-range vdW correction.
563    *
564    * @return True if it is on.
565    */
566   public boolean getDoLongRangeCorrection() {
567     return doLongRangeCorrection;
568   }
569 
570   /**
571    * Get the total Van der Waals potential energy.
572    *
573    * @return The energy.
574    * @since 1.0
575    */
576   public double getEnergy() {
577     return sharedEnergy.get();
578   }
579 
580   /**
581    * Get the number of interacting pairs.
582    *
583    * @return The interaction count.
584    * @since 1.0
585    */
586   public int getInteractions() {
587     return sharedInteractions.get();
588   }
589 
590   /**
591    * {@inheritDoc}
592    */
593   @Override
594   public double getLambda() {
595     return lambda;
596   }
597 
598   /**
599    * {@inheritDoc}
600    */
601   @Override
602   public void setLambda(double lambda) {
603     assert (lambda >= 0.0 && lambda <= 1.0);
604     if (!lambdaTerm) {
605       return;
606     }
607     // Scale lambda if vdwLambdaEnd is less than 1.0.
608     if (lambda > vdwLambdaEnd) {
609       lambda = 1.0;
610     } else {
611       lambda = lambda / vdwLambdaEnd;
612     }
613     this.lambda = lambda;
614 
615     // Softcore constant in vdW denominator sc1 = alpha * (1-L)^2
616     sc1 = vdwLambdaAlpha * (1.0 - lambda) * (1.0 - lambda);
617     // d(sc1)/dL = -2 * alpha * (1-L)
618     dsc1dL = -2.0 * vdwLambdaAlpha * (1.0 - lambda);
619     // d^2(sc1)/dL^2 = 2 * alpha
620     d2sc1dL2 = 2.0 * vdwLambdaAlpha;
621     // Multiplicative vdW function: sc2 = L^beta
622     sc2 = pow(lambda, vdwLambdaExponent);
623     // d(sc2)/dL = beta * L^(beta-1)
624     dsc2dL = vdwLambdaExponent * pow(lambda, vdwLambdaExponent - 1.0);
625     if (vdwLambdaExponent >= 2.0) {
626       // d^2(sc2)/dL^2 = beta * (beta-1) * L^(beta-2)
627       d2sc2dL2 = vdwLambdaExponent * (vdwLambdaExponent - 1.0) * pow(lambda, vdwLambdaExponent - 2.0);
628     } else {
629       d2sc2dL2 = 0.0;
630     }
631 
632     // If LambdaFactors are in OST mode, update them now.
633     if (!esvTerm) {
634       for (LambdaFactors lf : lambdaFactors) {
635         lf.setFactors();
636       }
637     }
638 
639     initSoftCore();
640 
641     // Redo the long range correction.
642     if (doLongRangeCorrection) {
643       longRangeCorrection = computeLongRangeCorrection();
644     } else {
645       longRangeCorrection = 0.0;
646     }
647   }
648 
649   /**
650    * Getter for the field <code>bondMask</code>.
651    *
652    * @return an array of {@link int} objects.
653    */
654   public int[][] getMask12() {
655     return mask12;
656   }
657 
658   /**
659    * Getter for the field <code>angleMask</code>.
660    *
661    * @return an array of {@link int} objects.
662    */
663   public int[][] getMask13() {
664     return mask13;
665   }
666 
667   /**
668    * Getter for the field <code>torsionMask</code>.
669    *
670    * @return an array of {@link int} objects.
671    */
672   public int[][] getMask14() {
673     return mask14;
674   }
675 
676   /**
677    * Allow sharing the of the VanDerWaals NeighborList with ParticleMeshEwald.
678    *
679    * @return The NeighborList.
680    */
681   public NeighborList getNeighborList() {
682     return neighborList;
683   }
684 
685   /**
686    * Getter for the field <code>neighborLists</code>.
687    *
688    * @return an array of int.
689    */
690   public int[][][] getNeighborLists() {
691     return neighborLists;
692   }
693 
694   /**
695    * Get details of the non-bonded cutoff.
696    *
697    * @return a {@link ffx.potential.nonbonded.NonbondedCutoff} object.
698    */
699   public NonbondedCutoff getNonbondedCutoff() {
700     return nonbondedCutoff;
701   }
702 
703   /**
704    * Get the reduction index.
705    *
706    * @return an array of {@link int} objects.
707    */
708   public int[] getReductionIndex() {
709     return reductionIndex;
710   }
711 
712   /**
713    * getVDWForm.
714    *
715    * @return a {@link ffx.potential.nonbonded.VanDerWaalsForm} object.
716    */
717   public VanDerWaalsForm getVDWForm() {
718     return vdwForm;
719   }
720 
721   /**
722    * {@inheritDoc}
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    * {@inheritDoc}
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    * {@inheritDoc}
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    * {@inheritDoc}
763    *
764    * <p>Remove the masking rules for 1-2, 1-3 and 1-4 interactions.
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    * Setter for the field <code>atoms</code>.
786    *
787    * @param atoms         an array of {@link ffx.potential.bonded.Atom} objects.
788    * @param molecule      an array of {@link int} objects.
789    * @param neuralNetwork an array of flags to indicate if the atom is treated by a neural
790    *                      network.
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    * If the crystal being passed in is not equal to the current crystal, then some Van der Waals data
809    * structures may need to updated. If <code>nSymm</code> has changed, update arrays dimensioned by
810    * nSymm. Finally, rebuild the neighbor-lists.
811    *
812    * @param crystal The new crystal instance defining the symmetry and boundary conditions.
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       // Allocate memory if necessary.
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    * {@inheritDoc}
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    * Allocate coordinate arrays and set up reduction indices and values.
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     // Initialize all atoms to be used in the energy.
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         // Find the vdW parameters from the AtomType.
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         // Atom indexes start at 1
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       // Collect 1-2 interactions.
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       // Collect 1-3 interactions.
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       // Collect 1-4 interactions.
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    * buildNeighborList.
985    *
986    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
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    * Computes the long range van der Waals correction to the energy via numerical integration
1002    *
1003    * @return Long range correction (kcal/mol)
1004    * @see "M. P. Allen and D. J. Tildesley, "Computer Simulation of Liquids, 2nd Ed.", Oxford
1005    * University Press, 2017, Section 2.8"
1006    */
1007   private double computeLongRangeCorrection() {
1008     // Need to treat esvLambda chain terms below before you can do this.
1009     if (esvTerm) {
1010       throw new UnsupportedOperationException();
1011     }
1012 
1013     // Only applicable under periodic boundary conditions.
1014     if (crystal.aperiodic()) {
1015       return 0.0;
1016     }
1017 
1018     /*
1019      Integrate to maxR = 100 Angstroms or ~33 sigma.
1020      nDelta integration steps: (100.0 - cutoff) * 2.0
1021      rDelta integration step size: (100.0 - cutoff) / nDelta
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     // Count the number of classes and their frequencies
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     // Loop over vdW classes.
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               // Arbitrary power of "N" for the attractive dispersion.
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            Apply one minus the multiplicative switch if the
1088            interaction distance is less than the end of the
1089            switching window.
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         // Normal correction
1111         total += radCount[i] * radCount[j] * trapezoid;
1112         // Correct for softCore vdW that are being turned off.
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     // Divide by the volume of the asymmetric unit.
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    * Log the Van der Waals interaction.
1130    *
1131    * @param i    Atom i.
1132    * @param k    Atom j.
1133    * @param minr The minimum energy vdW separation distance.
1134    * @param r    The distance rij.
1135    * @param eij  The interaction energy.
1136    * @since 1.0
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     // Initialize the softcore atom masks.
1167     for (int i = 0; i < nAtoms; i++) {
1168       if (isSoft[i]) {
1169         // Outer loop atom hard, inner loop atom soft.
1170         softCore[HARD][i] = true;
1171         // Both soft: full intramolecular ligand interactions.
1172         softCore[SOFT][i] = false;
1173       } else {
1174         // Both hard: full interaction between atoms.
1175         softCore[HARD][i] = false;
1176         // Outer loop atom soft, inner loop atom hard.
1177         softCore[SOFT][i] = true;
1178       }
1179     }
1180   }
1181 
1182   /**
1183    * The trick: The setFactors(i,k) method is called every time through the inner VdW loop, avoiding
1184    * an "if (esv)" branch statement. A plain OST run will have an object of type LambdaFactorsOST
1185    * instead, which contains an empty version of setFactors(i,k). The OST version sets new factors
1186    * only on lambda updates, in setLambda().
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      * Overriden by the OST version which updates only during setLambda().
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      * Timing variables.
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       // Log timings.
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       // Locally initialize the Loops to help with NUMA?
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       // Initialize and expand coordinates.
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       // Build the neighbor-list (if necessary) using reduced coordinates.
1321       barrier(neighborListAction);
1322       if (forceNeighborListRebuild) {
1323         return;
1324       }
1325 
1326       // Compute Van der Waals energy and gradient.
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       // Reduce derivatives.
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      * {@inheritDoc}
1367      *
1368      * @since 1.0
1369      */
1370     @Override
1371     public void start() {
1372       // Start timing.
1373       vdwTimeTotal = -System.nanoTime();
1374 
1375       // Initialize the shared variables.
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      * Update the local coordinate array and initialize reduction variables.
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           // Load the vdw type.
1419           VDWType vdwType = atom.getVDWType();
1420           if (vdwType == null) {
1421             logger.severe(" No VdW type for atom " + atom);
1422             return;
1423           }
1424 
1425           // Update the atom class.
1426           atomClass[i] = vdwType.atomClass;
1427 
1428           // Set reduction values.
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             // Atom indexes start at 1
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           // Collect 1-2 interactions.
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           // Collect 1-3 interactions.
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           // Collect 1-4 interactions.
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          Set up the local coordinate array for the asymmetric unit,
1512          applying reduction factors to the hydrogen atoms.
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             // Check if the atom is at a special position.
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      * The Van der Waals loop class contains methods and thread local variables used to evaluate the
1596      * Van der Waals energy and gradients with respect to atomic coordinates.
1597      *
1598      * @author Michael J. Schnieders
1599      * @since 1.0
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         // Reduce the energy, interaction count and gradients into the shared variables.
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         // neighborLists array: [nSymm][nAtoms][nNeighbors]
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           // Flag to indicate if atom i is effected by an extended system variable.
1649           final boolean esvi = esvTerm && esvSystem.isTitratingHydrogen(i);
1650           if (esvTerm) {
1651             esvSystem.getVdwPrefactor(i, esvVdwPrefactori);
1652           }
1653           int i3 = i * 3;
1654           // Atomic coordinates of atom i, including application of reduction factors.
1655           final double xi = reducedXYZ[i3++];
1656           final double yi = reducedXYZ[i3++];
1657           final double zi = reducedXYZ[i3];
1658           // If this hydrogen atomic position is reduced, this is the index of its heavy atom.
1659           final int redi = reductionIndex[i];
1660           // If this hydrogen atomic position is reduced, this is the reduction factor.
1661           final double redv = reductionValue[i];
1662           final double rediv = 1.0 - redv;
1663           final int classI = atomClass[i];
1664           // Molecule index for atom i.
1665           final int moleculei = molecule[i];
1666           // Neural network identity.
1667           boolean iNN = neuralNetwork[i];
1668           // Zero out local gradient accumulation variables.
1669           double gxi = 0.0;
1670           double gyi = 0.0;
1671           double gzi = 0.0;
1672           // Zero out local gradient accumulation variables for reduced atoms.
1673           double gxredi = 0.0;
1674           double gyredi = 0.0;
1675           double gzredi = 0.0;
1676           // Zero out local variables for accumulation of dU/dX/dL.
1677           double lxi = 0.0;
1678           double lyi = 0.0;
1679           double lzi = 0.0;
1680           // Zero out local variables for accumulation of dU/dX/dL for reduced atoms.
1681           double lxredi = 0.0;
1682           double lyredi = 0.0;
1683           double lzredi = 0.0;
1684           // Collect information about special vdW 1-4 interactions,
1685           // and fill the mask array due to application of 1-2, 1-3 and 1-4 vdW scale factors.
1686           applyMask(i, vdw14, mask);
1687           // Default is that the outer loop atom is hard.
1688           boolean[] softCorei = softCore[HARD];
1689           if (isSoft[i]) {
1690             softCorei = softCore[SOFT];
1691           }
1692           // Loop over the neighbor list.
1693           final int[] neighbors = list[i];
1694           for (final int k : neighbors) {
1695             // Check that atom k is in use.
1696             // Do not compute vdW interactions between asymmetric unit neural network atoms.
1697             if (!use[k] || (iNN && neuralNetwork[k])) {
1698               continue;
1699             }
1700             // Flag to indicate if atom k is effected by an extended system variable.
1701             final boolean esvk = esvTerm && esvSystem.isTitratingHydrogen(k);
1702             if (esvTerm) {
1703               esvSystem.getVdwPrefactor(k, esvVdwPrefactork);
1704             }
1705             int k3 = k * 3;
1706             // Atomic coordinates of atom k, including application of reduction factors.
1707             final double xk = xyzS[k3++];
1708             final double yk = xyzS[k3++];
1709             final double zk = xyzS[k3];
1710             // Compute the separation vector between atom i and k.
1711             dx_local[0] = xi - xk;
1712             dx_local[1] = yi - yk;
1713             dx_local[2] = zi - zk;
1714             // Apply the minimum image convention (if periodic).
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               // Some force fields such as CHARMM have special 1-4 vdW parameters.
1720               irv = vdwForm.getCombinedInverseRmin14(classI, classK);
1721             }
1722             if (r2 <= nonbondedCutoff.off2 && mask[k] > 0 && irv > 0) {
1723               // Compute the separation distance.
1724               final double r = sqrt(r2);
1725               // Check if i and k are part of the same molecule.
1726               boolean sameMolecule = (moleculei == molecule[k]);
1727               boolean soft = softCorei[k];
1728               // If both atoms are softcore, respect the intermolecularSoftcore and
1729               // intramolecularSoftcore flags.
1730               if (isSoft[i] && isSoft[k]) {
1731                 if (intermolecularSoftcore && !sameMolecule) {
1732                   soft = true;
1733                 } else if (intramolecularSoftcore && sameMolecule) {
1734                   soft = true;
1735                 }
1736               }
1737               // Hide these global variable names for thread safety.
1738               final double sc1, dsc1dL, d2sc1dL2;
1739               final double sc2, dsc2dL, d2sc2dL2;
1740               if (soft) {
1741                 /*
1742                 The setFactors(i,k) method is empty unless ESVs are present.
1743                 If OST lambda present, lambdaFactors will already have been updated during setLambda().
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                Calculate Van der Waals interaction energy. Notation from
1763                "The structure, thermodynamics, and solubility of organic
1764                crystals from simulation with a polarizable force field."
1765                J. Chem. Theory Comput. 8, 1721–1736 (2012).
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                Apply a multiplicative switch if the interaction
1787                distance is greater than the beginning of the taper.
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               // log(i, k, ev, r, eik * taper * esvik);
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                 // Assign this gradient to attached ESVs.
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           // Compute the total transformation operator: R = ToCart * Rot * ToFrac.
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             // Default is that the outer loop atom is hard.
1932             boolean[] softCorei = softCore[HARD];
1933             if (isSoft[i]) {
1934               softCorei = softCore[SOFT];
1935             }
1936             // Loop over the neighbor list.
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               // Hide these global variable names for thread safety.
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                 // Apply a multiplicative switch if the interaction distance is greater than
1996                 // the beginning of the taper.
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                   // Apply the transpose of the transformation operator.
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                   // Assign this gradient to attached ESVs.
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                   // Apply the transpose of the transformation operator.
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      * Reduce Van der Waals gradient.
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      * Build the NeighborList.
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 }