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-2021.
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;
39  
40  import static ffx.potential.nonbonded.ParticleMeshEwald.DEFAULT_EWALD_CUTOFF;
41  import static ffx.potential.nonbonded.VanDerWaals.DEFAULT_VDW_CUTOFF;
42  import static ffx.potential.parameters.ForceField.toEnumForm;
43  import static ffx.potential.parsers.XYZFileFilter.isXYZ;
44  import static ffx.utilities.KeywordGroup.NonBondedCutoff;
45  import static ffx.utilities.KeywordGroup.PotentialFunctionSelection;
46  import static java.lang.Double.isInfinite;
47  import static java.lang.Double.isNaN;
48  import static java.lang.String.format;
49  import static java.util.Arrays.fill;
50  import static java.util.Arrays.sort;
51  import static org.apache.commons.io.FilenameUtils.removeExtension;
52  import static org.apache.commons.math3.util.FastMath.max;
53  import static org.apache.commons.math3.util.FastMath.min;
54  import static org.apache.commons.math3.util.FastMath.sqrt;
55  
56  import edu.rit.pj.IntegerForLoop;
57  import edu.rit.pj.IntegerSchedule;
58  import edu.rit.pj.ParallelRegion;
59  import edu.rit.pj.ParallelTeam;
60  import edu.rit.pj.reduction.SharedDouble;
61  import ffx.crystal.Crystal;
62  import ffx.crystal.CrystalPotential;
63  import ffx.crystal.LatticeSystem;
64  import ffx.crystal.NCSCrystal;
65  import ffx.crystal.ReplicatesCrystal;
66  import ffx.crystal.SpaceGroup;
67  import ffx.crystal.SpaceGroupDefinitions;
68  import ffx.crystal.SymOp;
69  import ffx.numerics.Constraint;
70  import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
71  import ffx.numerics.atomic.AtomicDoubleArray3D;
72  import ffx.numerics.switching.ConstantSwitch;
73  import ffx.numerics.switching.UnivariateFunctionFactory;
74  import ffx.numerics.switching.UnivariateSwitchingFunction;
75  import ffx.potential.bonded.Angle;
76  import ffx.potential.bonded.AngleTorsion;
77  import ffx.potential.bonded.Atom;
78  import ffx.potential.bonded.Atom.Resolution;
79  import ffx.potential.bonded.Bond;
80  import ffx.potential.bonded.BondedTerm;
81  import ffx.potential.bonded.ImproperTorsion;
82  import ffx.potential.bonded.LambdaInterface;
83  import ffx.potential.bonded.MSNode;
84  import ffx.potential.bonded.MultiResidue;
85  import ffx.potential.bonded.OutOfPlaneBend;
86  import ffx.potential.bonded.PiOrbitalTorsion;
87  import ffx.potential.bonded.RelativeSolvation;
88  import ffx.potential.bonded.RelativeSolvation.SolvationLibrary;
89  import ffx.potential.bonded.Residue;
90  import ffx.potential.bonded.RestraintBond;
91  import ffx.potential.bonded.RestraintTorsion;
92  import ffx.potential.bonded.StretchBend;
93  import ffx.potential.bonded.StretchTorsion;
94  import ffx.potential.bonded.Torsion;
95  import ffx.potential.bonded.TorsionTorsion;
96  import ffx.potential.bonded.UreyBradley;
97  import ffx.potential.constraint.CcmaConstraint;
98  import ffx.potential.constraint.SettleConstraint;
99  import ffx.potential.extended.ExtendedSystem;
100 import ffx.potential.nonbonded.COMRestraint;
101 import ffx.potential.nonbonded.CoordRestraint;
102 import ffx.potential.nonbonded.GeneralizedKirkwood;
103 import ffx.potential.nonbonded.NCSRestraint;
104 import ffx.potential.nonbonded.ParticleMeshEwald;
105 import ffx.potential.nonbonded.RestrainGroups;
106 import ffx.potential.nonbonded.VanDerWaals;
107 import ffx.potential.nonbonded.VanDerWaalsTornado;
108 import ffx.potential.parameters.AngleType.AngleMode;
109 import ffx.potential.parameters.BondType;
110 import ffx.potential.parameters.ForceField;
111 import ffx.potential.parameters.ForceField.ELEC_FORM;
112 import ffx.potential.parameters.TorsionType;
113 import ffx.potential.utils.ConvexHullOps;
114 import ffx.potential.utils.EnergyException;
115 import ffx.potential.utils.PotentialsFunctions;
116 import ffx.potential.utils.PotentialsUtils;
117 import ffx.utilities.Constants;
118 import ffx.utilities.FFXKeyword;
119 import java.io.File;
120 import java.time.LocalDateTime;
121 import java.time.format.DateTimeFormatter;
122 import java.util.ArrayList;
123 import java.util.Arrays;
124 import java.util.Collections;
125 import java.util.HashSet;
126 import java.util.List;
127 import java.util.Set;
128 import java.util.logging.Level;
129 import java.util.logging.Logger;
130 import java.util.stream.Collectors;
131 import java.util.stream.Stream;
132 import org.apache.commons.configuration2.CompositeConfiguration;
133 
134 /**
135  * Compute the potential energy and derivatives of a molecular system described by a force field.
136  *
137  * @author Michael J. Schnieders
138  * @since 1.0
139  */
140 public class ForceFieldEnergy implements CrystalPotential, LambdaInterface {
141 
142   /** Default tolerance for numerical methods of solving constraints. */
143   public static final double DEFAULT_CONSTRAINT_TOLERANCE = 1E-4;
144   /** A Logger for the ForceFieldEnergy class. */
145   private static final Logger logger = Logger.getLogger(ForceFieldEnergy.class.getName());
146   /** Convert from nanoseconds to seconds. */
147   private static final double toSeconds = 1.0e-9;
148   /** The MolecularAssembly associated with this force field energy. */
149   protected final MolecularAssembly molecularAssembly;
150   /**
151    * If the absolute value of a gradient component is greater than "maxDebugGradient", verbose
152    * logging results.
153    */
154   final double maxDebugGradient;
155   /** The Parallel Java ParallelTeam instance. */
156   private final ParallelTeam parallelTeam;
157   /** An array of Coordinate Restraint terms. */
158   private final List<CoordRestraint> coordRestraints;
159   /** An NCS restraint term. */
160   private final NCSRestraint ncsRestraint;
161   /** A coordinate restraint term. */
162   private final CoordRestraint autoCoordRestraint;
163   /** A Center-of-Mass restraint term. */
164   private final COMRestraint comRestraint;
165   /** Restrain groups */
166   private final RestrainGroups restrainGroups;
167   /** Non-Bonded van der Waals energy. */
168   private final VanDerWaals vanderWaals;
169 
170   private final List<Constraint> constraints;
171 
172   @FFXKeyword(name = "vdw-cutoff", keywordGroup = NonBondedCutoff, defaultValue = "12.0",
173       description =
174           "Sets the cutoff distance value in Angstroms for van der Waals potential energy interactions. "
175               + "The energy for any pair of van der Waals sites beyond the cutoff distance will be set to zero. "
176               + "Other properties can be used to define the smoothing scheme near the cutoff distance. "
177               + "The default cutoff distance in the absence of the vdw-cutoff keyword is infinite for nonperiodic "
178               + "systems and 12.0 for periodic systems.")
179   private double vdwCutoff;
180 
181   @FFXKeyword(name = "ewald-cutoff", keywordGroup = NonBondedCutoff, defaultValue = "7.0",
182       description =
183           "Sets the value in Angstroms of the real-space distance cutoff for use during Ewald summation "
184               + "By default, in the absence of the ewald-cutoff property, a value of 7.0 is used.")
185   private double ewaldCutoff;
186 
187   @FFXKeyword(name = "gk-cutoff", keywordGroup = NonBondedCutoff, defaultValue = "12.0",
188       description =
189           "Sets the value in Angstroms of the generalized Kirkwood distance cutoff for use during implicit solvent simulations. "
190               + "By default, in the absence of the gk-cutoff property, no cutoff is used under aperiodic boundary conditions and the vdw-cutoff is used under PBC.")
191   private double gkCutoff;
192 
193   private static final double DEFAULT_LIST_BUFFER = 2.0;
194 
195   @FFXKeyword(name = "list-buffer", keywordGroup = NonBondedCutoff, defaultValue = "2.0",
196       description = "Sets the size of the neighbor list buffer in Angstroms for potential energy functions. "
197           + "This value is added to the actual cutoff distance to determine which pairs will be kept on the neighbor list. This buffer value is used for all potential function neighbor lists. "
198           + "The default value in the absence of the list-buffer keyword is 2.0 Angstroms.")
199   private double listBuffer;
200 
201   /** 2.0 times the neighbor list cutoff. */
202   private final double cutOff2;
203   /** The non-bonded cut-off plus buffer distance (Angstroms). */
204   private final double cutoffPlusBuffer;
205   /** An array of Stretch-Torsion terms. */
206   private final StretchTorsion[] stretchTorsions;
207   /** An array of Angle-Torsion terms. */
208   private final AngleTorsion[] angleTorsions;
209   /** Particle-Mesh Ewald electrostatic energy. */
210   private final ParticleMeshEwald particleMeshEwald;
211   /** Number of Angle-Torsion terms in the system. */
212   private final int nAngleTorsions;
213   /** Number of Stretch-Torsion terms in the system. */
214   private final int nStretchTorsions;
215   /** Original state of the Bond energy term flag. */
216   private final boolean bondTermOrig;
217   /** Original state of the Angle energy term flag. */
218   private final boolean angleTermOrig;
219   /** Original state of the Stretch-Bend energy term flag. */
220   private final boolean stretchBendTermOrig;
221   /** Original state of the Urey-Bradley energy term flag. */
222   private final boolean ureyBradleyTermOrig;
223   /** Original state of the Out-of-Plane Bend energy term flag. */
224   private final boolean outOfPlaneBendTermOrig;
225   /** Original state of the Torsion energy term flag. */
226   private final boolean torsionTermOrig;
227   /** Original state of the Stretch-Torsion energy term flag. */
228   private final boolean stretchTorsionTermOrig;
229   /** Original state of the Angle-Torsion energy term flag. */
230   private final boolean angleTorsionTermOrig;
231   /** Original state of the Improper Torsion energy term flag. */
232   private final boolean improperTorsionTermOrig;
233   /** Original state of the Pi-Orbital Torsion energy term flag. */
234   private final boolean piOrbitalTorsionTermOrig;
235   /** Original state of the Torsion-Torsion energy term flag. */
236   private final boolean torsionTorsionTermOrig;
237   /** Original state of the Restraint Bond energy term flag. */
238   private final boolean restraintBondTermOrig;
239   /** Original state of the van der Waals energy term flag. */
240   private final boolean vanderWaalsTermOrig;
241   /** Original state of the multipole energy term flag. */
242   private final boolean multipoleTermOrig;
243   /** Original state of the polarization energy term flag. */
244   private final boolean polarizationTermOrig;
245   /** Original state of the GK energy term flag. */
246   private final boolean generalizedKirkwoodTermOrig;
247   private final boolean rTorsTermOrig;
248   /** Flag to indicate hydrogen bonded terms should be scaled up. */
249   private final boolean rigidHydrogens;
250   /** Indicates application of lambda scaling to all Torsion based energy terms. */
251   private final boolean lambdaTorsions;
252   /** Relative solvation term (TODO: needs further testing). */
253   private final RelativeSolvation relativeSolvation;
254   private final boolean relativeSolvationTerm;
255   private final Platform platform = Platform.FFX;
256 
257   /** Indicates use of the Lambda state variable. */
258   @FFXKeyword(name = "lambdaterm", clazz = Boolean.class,
259       keywordGroup = PotentialFunctionSelection, defaultValue = "false",
260       description = "Specifies use of the Lambda state variable."
261   )
262   protected boolean lambdaTerm;
263   /** Current value of the Lambda state variable. */
264   private double lambda = 1.0;
265 
266 
267   /** Optimization scaling value to use for each degree of freedom. */
268   protected double[] optimizationScaling = null;
269   /** Indicates only bonded energy terms effected by Lambda should be evaluated. */
270   boolean lambdaBondedTerms = false;
271   /** Indicates all bonded energy terms should be evaluated if lambdaBondedTerms is true. */
272   boolean lambdaAllBondedTerms = false;
273 
274   /** Flag to indicate proper shutdown of the ForceFieldEnergy. */
275   boolean destroyed = false;
276   /**
277    * An instance of the STATE enumeration to specify calculation of slowly varying energy terms, fast
278    * varying or both.
279    */
280   private STATE state = STATE.BOTH;
281   /** The array of Atoms being evaluated. */
282   private Atom[] atoms;
283   /** An array of Bond terms. */
284   private Bond[] bonds;
285   /** An array of Angle terms. */
286   private Angle[] angles;
287   /** An array of Stretch-Bend terms. */
288   private StretchBend[] stretchBends;
289   /** An array of Urey-Bradley terms. */
290   private UreyBradley[] ureyBradleys;
291   /** An array of Out of Plane Bend terms. */
292   private OutOfPlaneBend[] outOfPlaneBends;
293   /** An array of Torsion-Torsion terms. */
294   private TorsionTorsion[] torsionTorsions;
295   /** An array of Improper Torsion terms. */
296   private ImproperTorsion[] improperTorsions;
297   /** An array of Torsion terms. */
298   private Torsion[] torsions;
299   /** An array of Pi-Orbital Torsion terms. */
300   private PiOrbitalTorsion[] piOrbitalTorsions;
301   /** An array of Bond Restraint terms. */
302   private RestraintBond[] restraintBonds;
303   protected RestraintTorsion[] rTors;
304   /** Number of atoms in the system. */
305   private int nAtoms;
306   /** Number of bond terms in the system. */
307   private int nBonds;
308   /** Number of angle terms in the system. */
309   private int nAngles;
310   /** Number of stretch-bend terms in the system. */
311   private int nStretchBends;
312   /** Number of Urey-Bradley terms in the system. */
313   private int nUreyBradleys;
314   /** Number of Out of Plane Bend terms in the system. */
315   private int nOutOfPlaneBends;
316   /** Number of Improper Torsion terms in the system. */
317   private int nImproperTorsions;
318   /** Number of Pi-Orbital Torsion terms in the system. */
319   private int nPiOrbitalTorsions;
320   /** Number of Torsion-Torsion terms in the system. */
321   private int nTorsionTorsions;
322   /** Number of Torsion terms in the system. */
323   private int nTorsions;
324   /** Number of Restraint Bond terms in the system. */
325   private int nRestraintBonds = 0;
326   /** Number of Restraint Bond terms in the system. */
327   private int nRestrainGroups = 0;
328   /** Number of van der Waals interactions evaluated. */
329   private int nVanDerWaalInteractions;
330   /** Number of electrostatic interactions evaluated. */
331   private int nPermanentInteractions;
332   /** Number of implicit solvent interactions evaluated. */
333   private int nGKInteractions;
334   private int nRestTors = 0;
335   /** The boundary conditions used when evaluating the force field energy. */
336   private Crystal crystal;
337   /** A Parallel Java Region used to evaluate Bonded energy values. */
338   private BondedRegion bondedRegion;
339 
340   /** Specifies use of the bond stretch potential. */
341   @FFXKeyword(name = "bondterm", clazz = Boolean.class,
342       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
343       description = "Specifies use of the bond stretch potential."
344   )
345   private boolean bondTerm;
346 
347   /** Specifies use of the angle bend potential. */
348   @FFXKeyword(name = "angleterm", clazz = Boolean.class,
349       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
350       description = "Specifies use of the angle bend potential.")
351   private boolean angleTerm;
352 
353   /** Evaluate Stretch-Bend energy terms. */
354   @FFXKeyword(name = "strbndterm", clazz = Boolean.class,
355       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
356       description = "Specifies use of the stretch-bend potential.")
357   private boolean stretchBendTerm;
358 
359   /** Evaluate Urey-Bradley energy terms. */
360   @FFXKeyword(name = "ureyterm", clazz = Boolean.class,
361       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
362       description = "Specifies use of the Urey-Bradley potential.")
363   private boolean ureyBradleyTerm;
364 
365   /** Evaluate Out of Plane Bend energy terms. */
366   @FFXKeyword(name = "opbendterm", clazz = Boolean.class,
367       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
368       description = "Specifies use of the out-of-plane potential.")
369   private boolean outOfPlaneBendTerm;
370 
371   /** Evaluate Torsion energy terms. */
372   @FFXKeyword(name = "torsionterm", clazz = Boolean.class,
373       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
374       description = "Specifies use of the torsional potential.")
375   private boolean torsionTerm;
376 
377   /** Evaluate Stretch-Torsion energy terms. */
378   @FFXKeyword(name = "strtorterm", clazz = Boolean.class,
379       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
380       description = "Specifies use of the stretch-torsion potential.")
381   private boolean stretchTorsionTerm;
382 
383   /** Evaluate Angle-Torsion energy terms. */
384   @FFXKeyword(name = "angtorterm", clazz = Boolean.class,
385       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
386       description = "Specifies use of the angle-torsion potential.")
387   private boolean angleTorsionTerm;
388 
389   /** Evaluate Improper Torsion energy terms. */
390   @FFXKeyword(name = "imptorterm", clazz = Boolean.class,
391       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
392       description = "Specifies use of the improper torsion potential.")
393   private boolean improperTorsionTerm;
394 
395   /** Evaluate Pi-Orbital Torsion energy terms. */
396   @FFXKeyword(name = "pitorsterm", clazz = Boolean.class,
397       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
398       description = "Specifies use of the pi-system torsion potential.")
399   private boolean piOrbitalTorsionTerm;
400 
401   /** Evaluate Torsion-Torsion energy terms. */
402   @FFXKeyword(name = "tortorterm", clazz = Boolean.class,
403       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
404       description = "Specifies use of the pi-system torsion potential.")
405   private boolean torsionTorsionTerm;
406 
407   /** Evaluate van der Waals energy term. */
408   @FFXKeyword(name = "vdwterm", clazz = Boolean.class,
409       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
410       description = "Specifies use of the vdw der Waals potential. "
411           + "If set to false, all non-bonded terms are turned off.")
412   private boolean vanderWaalsTerm;
413 
414   /** Evaluate permanent multipole electrostatics energy term. */
415   @FFXKeyword(name = "mpoleterm", clazz = Boolean.class,
416       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
417       description = "Specifies use of the fixed charge electrostatic potential. "
418           + "Setting mpoleterm to false also turns off polarization and generalized Kirkwood, "
419           + "overriding the polarizeterm and gkterm properties.")
420   private boolean multipoleTerm;
421 
422   /** Evaluate polarization energy term. */
423   @FFXKeyword(name = "polarizeterm", clazz = Boolean.class,
424       keywordGroup = PotentialFunctionSelection, defaultValue = "true",
425       description = "Specifies use of the polarizable electrostatic potential. "
426           + "Setting polarizeterm to false overrides the polarization property.")
427   private boolean polarizationTerm;
428 
429   /** Evaluate generalized Kirkwood energy term. */
430   @FFXKeyword(name = "gkterm", clazz = Boolean.class,
431       keywordGroup = PotentialFunctionSelection, defaultValue = "false",
432       description = "Specifies use of generalized Kirkwood electrostatics.")
433   private boolean generalizedKirkwoodTerm;
434 
435   /** Evaluate COM energy term. */
436   private boolean comTerm;
437   /** Original state of the COM energy term flag. */
438   private boolean comTermOrig;
439   /** Evaluate NCS energy term. */
440   private boolean ncsTerm;
441   /** Original state of the NCS energy term flag. */
442   private boolean ncsTermOrig;
443   /** Evaluate Restrain energy term. */
444   private boolean restrainTerm;
445   /** Evaluate Restraint Bond energy terms. */
446   private boolean restraintBondTerm;
447   private boolean rTorsTerm;
448   /** Scale factor for increasing the strength of bonded terms involving hydrogen atoms. */
449   private double rigidScale;
450   /** The total Bond term energy. */
451   private double bondEnergy;
452   /** The RMSD of all Bond distances relative to their equilibrium values. */
453   private double bondRMSD;
454   /** The total Angle term energy. */
455   private double angleEnergy;
456   /** The RMSD of all Angle bends relative to their equilibrium values. */
457   private double angleRMSD;
458   /** The total Stretch-Bend term energy. */
459   private double stretchBendEnergy;
460   /** The total Urey-Bradley term energy. */
461   private double ureyBradleyEnergy;
462   /** The total Out-of-Plane term energy. */
463   private double outOfPlaneBendEnergy;
464   /** The total Torsion term energy. */
465   private double torsionEnergy;
466   /** The total Stretch-Torsion term energy. */
467   private double stretchTorsionEnergy;
468   /** The total Angle-Torsion term energy. */
469   private double angleTorsionEnergy;
470   /** The total Improper Torsion term energy. */
471   private double improperTorsionEnergy;
472   /** The total Pi-Orbital Torsion term energy. */
473   private double piOrbitalTorsionEnergy;
474   /** The total Torsion-Torsion term energy. */
475   private double torsionTorsionEnergy;
476   /** The total Restraint Bond term energy. */
477   private double restraintBondEnergy;
478   /** The total energy for all Bonded Energy terms. */
479   private double totalBondedEnergy;
480   /** The total van der Waals energy. */
481   private double vanDerWaalsEnergy;
482   /** The permanent Multipole energy. */
483   private double permanentMultipoleEnergy;
484   /** The permanent multipole real space energy. */
485   private double permanentRealSpaceEnergy;
486   /** The total polarization energy. */
487   private double polarizationEnergy;
488   /** The total electrostatic energy. */
489   private double totalMultipoleEnergy;
490   /** The total non-bonded energy (van der Waals plus electrostatic). */
491   private double totalNonBondedEnergy;
492   /** The solvation energy (GK). */
493   private double solvationEnergy;
494   /** The total NCS Energy. */
495   private double ncsEnergy;
496   /** The total Restrain Energy. */
497   private double restrainEnergy;
498   /** The total COM Restraint Energy. */
499   private double comRestraintEnergy;
500   /** The total COM Restraint Energy. */
501   private double restrainGroupEnergy;
502   private double rTorsEnergy;
503   /** The total system energy. */
504   private double totalEnergy;
505   /** Time to evaluate Bond terms. */
506   private long bondTime;
507   /** Time to evaluate Angle terms. */
508   private long angleTime;
509   /** Time to evaluate Stretch-Bend terms. */
510   private long stretchBendTime;
511   /** Time to evaluate Urey-Bradley terms. */
512   private long ureyBradleyTime;
513   /** Time to evaluate Out-Of-Plane Bend terms. */
514   private long outOfPlaneBendTime;
515   /** Time to evaluate Torsion terms. */
516   private long torsionTime;
517   /** Time to evaluate Angle-Torsion terms. */
518   private long angleTorsionTime;
519   /** Time to evaluate Stretch-Torsion terms. */
520   private long stretchTorsionTime;
521   /** Time to evaluate Pi-Orbital Torsion terms. */
522   private long piOrbitalTorsionTime;
523   /** Time to evaluate Improper Torsion terms. */
524   private long improperTorsionTime;
525   /** Time to evaluate Torsion-Torsion terms. */
526   private long torsionTorsionTime;
527   /** Time to evaluate van der Waals term. */
528   private long vanDerWaalsTime;
529   /** Time to evaluate electrostatics term. */
530   private long electrostaticTime;
531   /** Time to evaluate Restraint Bond term. */
532   private long restraintBondTime;
533   private long rTorsTime;
534   /** Original state of the Restrain energy term flag. */
535   private boolean restrainTermOrig;
536   /** Time to evaluate Center of Mass restraint term. */
537   private long comRestraintTime;
538   /** Time to evaluate restrain group term. */
539   private long restrainGroupTime;
540   /** Time to evaluate all energy terms. */
541   private long totalTime;
542   /** Evaluate Restrain Group energy term. */
543   private boolean restrainGroupTerm;
544   /** Value of each degree of freedom. */
545   private double[] xyz;
546   /** Original state of the Restrain Group energy term flag. */
547   private boolean restrainGroupTermOrig;
548   /** Indicate resolution of this ForceFieldEnergy (TODO: needs further testing). */
549   private Resolution resolution = Resolution.AMOEBA;
550   /** Constant pH extended system (TODO: needs further testing). */
551   private ExtendedSystem esvSystem = null;
552 
553   private boolean esvTerm;
554   private double esvBias;
555   /** Time to evaluate NCS term. */
556   private long ncsTime;
557   private int nRelativeSolvations;
558   /** Time to evaluate coordinate restraint term. */
559   private long coordRestraintTime;
560 
561   private double relativeSolvationEnergy;
562   /** Enable verbose printing if large energy gradient components are observed. */
563   private boolean printOnFailure;
564 
565   /**
566    * Constructor for ForceFieldEnergy.
567    *
568    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
569    */
570   protected ForceFieldEnergy(MolecularAssembly molecularAssembly) {
571     this(molecularAssembly, null);
572   }
573 
574   /**
575    * Constructor for ForceFieldEnergy.
576    *
577    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
578    * @param restraints list of {@link ffx.potential.nonbonded.CoordRestraint} objects.
579    */
580   protected ForceFieldEnergy(MolecularAssembly molecularAssembly, List<CoordRestraint> restraints) {
581     this(molecularAssembly, restraints, ParallelTeam.getDefaultThreadCount());
582   }
583 
584   /**
585    * Constructor for ForceFieldEnergy.
586    *
587    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
588    * @param restraints a {@link java.util.List} object.
589    * @param numThreads a int.
590    */
591   protected ForceFieldEnergy(
592       MolecularAssembly molecularAssembly, List<CoordRestraint> restraints, int numThreads) {
593     // Get a reference to the sorted atom array.
594     this.molecularAssembly = molecularAssembly;
595     molecularAssembly.renameWaterProtons();
596     atoms = molecularAssembly.getAtomArray();
597     nAtoms = atoms.length;
598     xyz = new double[nAtoms * 3];
599 
600     // Check that atom ordering is correct and count the number of active atoms.
601     for (int i = 0; i < nAtoms; i++) {
602       int index = atoms[i].getXyzIndex() - 1;
603       assert (i == index);
604     }
605 
606     // Enforce that the number of threads be less than or equal to the number of atoms.
607     int nThreads = min(nAtoms, numThreads);
608     parallelTeam = new ParallelTeam(nThreads);
609 
610     ForceField forceField = molecularAssembly.getForceField();
611     CompositeConfiguration properties = forceField.getProperties();
612     String name = forceField.toString().toUpperCase();
613 
614     logger.info(format(" Constructing Force Field %s", name));
615     logger.info(format("\n SMP threads:                        %10d", nThreads));
616 
617     bondTerm = forceField.getBoolean("BONDTERM", true);
618     angleTerm = forceField.getBoolean("ANGLETERM", true);
619     stretchBendTerm = forceField.getBoolean("STRBNDTERM", true);
620     ureyBradleyTerm = forceField.getBoolean("UREYTERM", true);
621     outOfPlaneBendTerm = forceField.getBoolean("OPBENDTERM", true);
622     torsionTerm = forceField.getBoolean("TORSIONTERM", true);
623     stretchTorsionTerm = forceField.getBoolean("STRTORTERM", true);
624     angleTorsionTerm = forceField.getBoolean("ANGTORTERM", true);
625     piOrbitalTorsionTerm = forceField.getBoolean("PITORSTERM", true);
626     torsionTorsionTerm = forceField.getBoolean("TORTORTERM", true);
627     improperTorsionTerm = forceField.getBoolean("IMPROPERTERM", true);
628     vanderWaalsTerm = forceField.getBoolean("VDWTERM", true);
629 
630     boolean tornadoVM = forceField.getBoolean("tornado", false);
631     if (vanderWaalsTerm && !tornadoVM) {
632       multipoleTerm = forceField.getBoolean("MPOLETERM", true);
633       if (multipoleTerm) {
634         String polarizeString = forceField.getString("POLARIZATION", "NONE");
635         boolean defaultPolarizeTerm = !polarizeString.equalsIgnoreCase("NONE");
636         polarizationTerm = forceField.getBoolean("POLARIZETERM", defaultPolarizeTerm);
637         generalizedKirkwoodTerm = forceField.getBoolean("GKTERM", false);
638       } else {
639         // If multipole electrostatics is turned off, turn off all electrostatics.
640         polarizationTerm = false;
641         generalizedKirkwoodTerm = false;
642       }
643     } else {
644       // If van der Waals is turned off, turn off all non-bonded terms.
645       multipoleTerm = false;
646       polarizationTerm = false;
647       generalizedKirkwoodTerm = false;
648     }
649     restraintBondTerm = false;
650     lambdaTerm = forceField.getBoolean("LAMBDATERM", false);
651     restrainTerm = forceField.getBoolean("RESTRAINTERM", false);
652     comTerm = forceField.getBoolean("COMRESTRAINTERM", false);
653     lambdaTorsions = forceField.getBoolean("TORSION_LAMBDATERM", false);
654     printOnFailure = forceField.getBoolean("PRINT_ON_FAILURE", false);
655 
656     if (properties.containsKey("restrain-groups")) {
657       restrainGroupTerm = true;
658     }
659 
660     // For RESPA
661     bondTermOrig = bondTerm;
662     angleTermOrig = angleTerm;
663     stretchBendTermOrig = stretchBendTerm;
664     ureyBradleyTermOrig = ureyBradleyTerm;
665     outOfPlaneBendTermOrig = outOfPlaneBendTerm;
666     torsionTermOrig = torsionTerm;
667     angleTorsionTermOrig = angleTorsionTerm;
668     stretchTorsionTermOrig = stretchTorsionTerm;
669     improperTorsionTermOrig = improperTorsionTerm;
670     piOrbitalTorsionTermOrig = piOrbitalTorsionTerm;
671     torsionTorsionTermOrig = torsionTorsionTerm;
672     restraintBondTermOrig = restraintBondTerm;
673     vanderWaalsTermOrig = vanderWaalsTerm;
674     multipoleTermOrig = multipoleTerm;
675     polarizationTermOrig = polarizationTerm;
676     generalizedKirkwoodTermOrig = generalizedKirkwoodTerm;
677     ncsTermOrig = ncsTerm;
678     restrainTermOrig = restrainTerm;
679     comTermOrig = comTerm;
680     restrainGroupTermOrig = restrainGroupTerm;
681 
682     // Determine the unit cell dimensions and Spacegroup
683     String spacegroup;
684     double a, b, c, alpha, beta, gamma;
685     boolean aperiodic;
686     try {
687       spacegroup = forceField.getString("SPACEGROUP", "P 1");
688       SpaceGroup sg = SpaceGroupDefinitions.spaceGroupFactory(spacegroup);
689       LatticeSystem latticeSystem = sg.latticeSystem;
690       a = forceField.getDouble("A_AXIS");
691       aperiodic = false;
692       b = forceField.getDouble("B_AXIS", latticeSystem.getDefaultBAxis(a));
693       c = forceField.getDouble("C_AXIS", latticeSystem.getDefaultCAxis(a, b));
694       alpha = forceField.getDouble("ALPHA", latticeSystem.getDefaultAlpha());
695       beta = forceField.getDouble("BETA", latticeSystem.getDefaultBeta());
696       gamma = forceField.getDouble("GAMMA", latticeSystem.getDefaultGamma());
697       if (!sg.latticeSystem.validParameters(a, b, c, alpha, beta, gamma)) {
698         logger.severe(" Check lattice parameters.");
699       }
700       if (a == 1.0 && b == 1.0 && c == 1.0) {
701         String message = " A-, B-, and C-axis values equal to 1.0.";
702         logger.info(message);
703         throw new Exception(message);
704       }
705     } catch (Exception e) {
706       logger.info(" The system will be treated as aperiodic.");
707       aperiodic = true;
708 
709       long convTime = -System.nanoTime();
710       double maxr = 10.0;
711       for (int i = 0; i < nAtoms - 1; i++) {
712         Atom ai = atoms[i];
713         for (int j = 1; j < nAtoms; j++) {
714           Atom aj = atoms[j];
715           double dx = ai.getX() - aj.getX();
716           double dy = ai.getY() - aj.getY();
717           double dz = ai.getZ() - aj.getZ();
718           double r = sqrt(dx * dx + dy * dy + dz * dz);
719           maxr = max(r, maxr);
720         }
721       }
722       convTime += System.nanoTime();
723       if (nAtoms > 10) {
724         long hullTime = -System.nanoTime();
725         double maxC = ConvexHullOps.maxDist(ConvexHullOps.constructHull(atoms));
726         maxC = Math.max(10.0, maxC);
727         hullTime += System.nanoTime();
728         double diff = maxr - maxC;
729         if (Math.abs(diff) > 1.0E-5) {
730           logger.warning(
731               String.format(
732                   " Max particle-particle distance by convex hull %12.6g "
733                       + "disagrees with max distance by brute-force loop %12.6g: difference %12.6g!",
734                   maxC, maxr, diff));
735         }
736         logger.fine(
737             String.format(
738                 " Time for convex hull calculation: %12.6g sec. By O(n^2) " + "loop: %12.6g sec.",
739                 hullTime * Constants.NS2SEC, convTime * Constants.NS2SEC));
740       } // At N < 4, the 3D convex hull may not even be defined properly, nevermind be more
741       // efficient.
742 
743       // Turn off reciprocal space calculations.
744       forceField.addProperty("EWALD_ALPHA", "0.0");
745 
746       // Specify some dummy values for the crystal.
747       spacegroup = "P1";
748       a = 2.0 * maxr;
749       b = 2.0 * maxr;
750       c = 2.0 * maxr;
751       alpha = 90.0;
752       beta = 90.0;
753       gamma = 90.0;
754     }
755     Crystal unitCell = new Crystal(a, b, c, alpha, beta, gamma, spacegroup);
756     unitCell.setAperiodic(aperiodic);
757 
758     // First set cutoffs assuming PBC.
759     vdwCutoff = DEFAULT_VDW_CUTOFF;
760     ewaldCutoff = DEFAULT_EWALD_CUTOFF;
761     gkCutoff = vdwCutoff;
762     if (unitCell.aperiodic()) {
763       // Default to no cutoffs for aperiodic systems.
764       vdwCutoff = Double.POSITIVE_INFINITY;
765       ewaldCutoff = Double.POSITIVE_INFINITY;
766       gkCutoff = Double.POSITIVE_INFINITY;
767     }
768     // Load user specified cutoffs.
769     vdwCutoff = forceField.getDouble("VDW_CUTOFF", vdwCutoff);
770     ewaldCutoff = forceField.getDouble("EWALD_CUTOFF", ewaldCutoff);
771     gkCutoff = forceField.getDouble("GK_CUTOFF", gkCutoff);
772 
773     // Neighbor list cutoff is the maximum of the vdW, Ewald and GK cutoffs (i.e. to use a single list).
774     double nlistCutoff = vdwCutoff;
775     if (multipoleTerm) {
776       nlistCutoff = max(vdwCutoff, ewaldCutoff);
777     }
778     if (generalizedKirkwoodTerm) {
779       nlistCutoff = max(nlistCutoff, gkCutoff);
780     }
781 
782     // Check for a frozen neighbor list.
783     boolean disabledNeighborUpdates = forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false);
784     if (disabledNeighborUpdates) {
785       logger.info(format(" Neighbor list updates disabled; interactions will "
786           + "only be calculated between atoms that started the simulation "
787           + "within a radius of %9.3g Angstroms of each other.", nlistCutoff));
788       // The individual cutoffs are infinity (i.e. they are limited by interactions in the list).
789       vdwCutoff = Double.POSITIVE_INFINITY;
790       ewaldCutoff = Double.POSITIVE_INFINITY;
791       gkCutoff = Double.POSITIVE_INFINITY;
792     }
793 
794     listBuffer = forceField.getDouble("LIST_BUFFER", DEFAULT_LIST_BUFFER);
795     cutoffPlusBuffer = nlistCutoff + listBuffer;
796     cutOff2 = 2.0 * cutoffPlusBuffer;
797     unitCell = configureNCS(forceField, unitCell);
798 
799     // If necessary, create a ReplicatesCrystal.
800     if (!aperiodic) {
801       this.crystal = ReplicatesCrystal.replicatesCrystalFactory(unitCell, cutOff2);
802       logger.info(format("\n Density:                                %6.3f (g/cc)",
803           crystal.getDensity(molecularAssembly.getMass())));
804       logger.info(crystal.toString());
805     } else {
806       this.crystal = unitCell;
807     }
808 
809     if (!unitCell.aperiodic() && unitCell.spaceGroup.number == 1) {
810       ncsTerm = forceField.getBoolean("NCSTERM", false);
811       ncsTermOrig = ncsTerm;
812     } else {
813       ncsTerm = false;
814       ncsTermOrig = false;
815     }
816 
817     rigidHydrogens = forceField.getBoolean("RIGID_HYDROGENS", false);
818     rigidScale = forceField.getDouble("RIGID_SCALE", 10.0);
819 
820     nRelativeSolvations = 0;
821     String relSolvLibrary = forceField.getString("RELATIVE_SOLVATION", "NONE").toUpperCase();
822     SolvationLibrary library = SolvationLibrary.valueOf(relSolvLibrary);
823     if (library == SolvationLibrary.AUTO) {
824       if (generalizedKirkwoodTerm && name.toUpperCase().contains("OPLS")) {
825         library = SolvationLibrary.MACCALLUM_TIP4P;
826         // Change when we have good Generalized Born numbers of our own for OPLS.
827       } else if (generalizedKirkwoodTerm) {
828         library = SolvationLibrary.GK;
829       } else {
830         library = SolvationLibrary.NONE;
831       }
832     }
833     if (library != SolvationLibrary.NONE) {
834       relativeSolvationTerm = true;
835       relativeSolvation = new RelativeSolvation(library, forceField);
836     } else {
837       relativeSolvationTerm = false;
838       relativeSolvation = null;
839     }
840 
841     boolean checkAllNodeCharges = forceField.getBoolean("CHECK_ALL_NODE_CHARGES", false);
842 
843     if (rigidScale <= 1.0) {
844       rigidScale = 1.0;
845     }
846 
847     logger.info("\n Bonded Terms");
848     if (rigidHydrogens && rigidScale > 1.0) {
849       logger.info(format("  Rigid hydrogens:                   %10.2f", rigidScale));
850     }
851 
852     // Collect, count, pack and sort bonds.
853     if (bondTerm) {
854       List<Bond> bond = molecularAssembly.getBondList();
855       nBonds = bond.size();
856       bonds = bond.toArray(new Bond[0]);
857       sort(bonds);
858       if (nBonds > 0) {
859         logger.info(format("  Bonds:                             %10d", nBonds));
860       }
861     } else {
862       nBonds = 0;
863       bonds = null;
864     }
865 
866     // Collect, count, pack and sort angles.
867     if (angleTerm) {
868       List<Angle> angle = molecularAssembly.getAngleList();
869       nAngles = angle.size();
870       angles = angle.toArray(new Angle[0]);
871       sort(angles);
872       if (nAngles > 0) {
873         logger.info(format("  Angles:                            %10d", nAngles));
874       }
875     } else {
876       nAngles = 0;
877       angles = null;
878     }
879 
880     // Collect, count, pack and sort stretch-bends.
881     if (stretchBendTerm) {
882       List<StretchBend> stretchBend = molecularAssembly.getStretchBendList();
883       nStretchBends = stretchBend.size();
884       stretchBends = stretchBend.toArray(new StretchBend[0]);
885       sort(stretchBends);
886       if (nStretchBends > 0) {
887         logger.info(format("  Stretch-Bends:                     %10d", nStretchBends));
888       }
889     } else {
890       nStretchBends = 0;
891       stretchBends = null;
892     }
893 
894     // Collect, count, pack and sort Urey-Bradleys.
895     if (ureyBradleyTerm) {
896       List<UreyBradley> ureyBradley = molecularAssembly.getUreyBradleyList();
897       nUreyBradleys = ureyBradley.size();
898       ureyBradleys = ureyBradley.toArray(new UreyBradley[0]);
899       sort(ureyBradleys);
900       if (nUreyBradleys > 0) {
901         logger.info(format("  Urey-Bradleys:                     %10d", nUreyBradleys));
902       }
903     } else {
904       nUreyBradleys = 0;
905       ureyBradleys = null;
906     }
907 
908     // Set a multiplier on the force constants of bonded terms containing hydrogens.
909     if (rigidHydrogens) {
910       if (bonds != null) {
911         for (Bond bond : bonds) {
912           if (bond.containsHydrogen()) {
913             bond.setRigidScale(rigidScale);
914           }
915         }
916       }
917       if (angles != null) {
918         for (Angle angle : angles) {
919           if (angle.containsHydrogen()) {
920             angle.setRigidScale(rigidScale);
921           }
922         }
923       }
924       if (stretchBends != null) {
925         for (StretchBend stretchBend : stretchBends) {
926           if (stretchBend.containsHydrogen()) {
927             stretchBend.setRigidScale(rigidScale);
928           }
929         }
930       }
931       if (ureyBradleys != null) {
932         for (UreyBradley ureyBradley : ureyBradleys) {
933           if (ureyBradley.containsHydrogen()) {
934             ureyBradley.setRigidScale(rigidScale);
935           }
936         }
937       }
938     }
939 
940     // Collect, count, pack and sort out-of-plane bends.
941     if (outOfPlaneBendTerm) {
942       List<OutOfPlaneBend> outOfPlaneBend = molecularAssembly.getOutOfPlaneBendList();
943       nOutOfPlaneBends = outOfPlaneBend.size();
944       outOfPlaneBends = outOfPlaneBend.toArray(new OutOfPlaneBend[0]);
945       sort(outOfPlaneBends);
946       if (nOutOfPlaneBends > 0) {
947         logger.info(format("  Out-of-Plane Bends:                %10d", nOutOfPlaneBends));
948       }
949     } else {
950       nOutOfPlaneBends = 0;
951       outOfPlaneBends = null;
952     }
953 
954     double torsionScale = forceField.getDouble("TORSION_SCALE", 1.0);
955     if (torsionScale != 1.0) {
956       forceField.setTorsionScale(torsionScale);
957     }
958 
959     // Collect, count, pack and sort torsions.
960     if (torsionTerm) {
961       List<Torsion> torsion = molecularAssembly.getTorsionList();
962       nTorsions = torsion.size();
963       torsions = torsion.toArray(new Torsion[0]);
964       if (nTorsions > 0) {
965         if (torsionScale == 1.0) {
966           logger.info(format("  Torsions:                          %10d", nTorsions));
967         } else {
968           logger.info(format("  Torsions (%5.2f):                  %10d", torsionScale, nTorsions));
969         }
970       }
971     } else {
972       nTorsions = 0;
973       torsions = null;
974     }
975 
976     // Collect, count, pack and sort stretch torsions.
977     if (stretchTorsionTerm) {
978       List<StretchTorsion> stretchTorsion = molecularAssembly.getStretchTorsionList();
979       nStretchTorsions = stretchTorsion.size();
980       stretchTorsions = stretchTorsion.toArray(new StretchTorsion[0]);
981       if (nStretchTorsions > 0) {
982         if (torsionScale == 1.0) {
983           logger.info(format("  Stretch-Torsions:                  %10d", nStretchTorsions));
984         } else {
985           logger.info(
986               format("  Stretch-Torsions (%5.2f):          %10d", torsionScale, nStretchTorsions));
987         }
988       }
989     } else {
990       nStretchTorsions = 0;
991       stretchTorsions = null;
992     }
993 
994     // Collect, count, pack and sort angle torsions.
995     if (angleTorsionTerm) {
996       List<AngleTorsion> angleTorsion = molecularAssembly.getAngleTorsionList();
997       nAngleTorsions = angleTorsion.size();
998       angleTorsions = angleTorsion.toArray(new AngleTorsion[0]);
999       if (nAngleTorsions > 0) {
1000         if (torsionScale == 1.0) {
1001           logger.info(format("  Angle-Torsions:                    %10d", nAngleTorsions));
1002         } else {
1003           logger.info(
1004               format("  Angle-Torsions (%5.2f):            %10d", torsionScale, nAngleTorsions));
1005         }
1006       }
1007     } else {
1008       nAngleTorsions = 0;
1009       angleTorsions = null;
1010     }
1011 
1012     // Collect, count, pack and sort pi-orbital torsions.
1013     if (piOrbitalTorsionTerm) {
1014       List<PiOrbitalTorsion> piOrbitalTorsion = molecularAssembly.getPiOrbitalTorsionList();
1015       nPiOrbitalTorsions = piOrbitalTorsion.size();
1016       piOrbitalTorsions = piOrbitalTorsion.toArray(new PiOrbitalTorsion[0]);
1017       if (nPiOrbitalTorsions > 0) {
1018         if (torsionScale == 1.0) {
1019           logger.info(format("  Pi-Orbital Torsions:               %10d", nPiOrbitalTorsions));
1020         } else {
1021           logger.info(
1022               format(
1023                   "  Pi-Orbital Torsions (%5.2f):       %10d", torsionScale, nPiOrbitalTorsions));
1024         }
1025       }
1026     } else {
1027       nPiOrbitalTorsions = 0;
1028       piOrbitalTorsions = null;
1029     }
1030 
1031     // Collect, count, pack and sort torsion-torsions.
1032     if (torsionTorsionTerm) {
1033       List<TorsionTorsion> torsionTorsion = molecularAssembly.getTorsionTorsionList();
1034       nTorsionTorsions = torsionTorsion.size();
1035       torsionTorsions = torsionTorsion.toArray(new TorsionTorsion[0]);
1036       if (nTorsionTorsions > 0) {
1037         logger.info(format("  Torsion-Torsions:                  %10d", nTorsionTorsions));
1038       }
1039     } else {
1040       nTorsionTorsions = 0;
1041       torsionTorsions = null;
1042     }
1043 
1044     // Collect, count, pack and sort improper torsions.
1045     if (improperTorsionTerm) {
1046       List<ImproperTorsion> improperTorsion = molecularAssembly.getImproperTorsionList();
1047       nImproperTorsions = improperTorsion.size();
1048       improperTorsions = improperTorsion.toArray(new ImproperTorsion[0]);
1049       if (nImproperTorsions > 0) {
1050         logger.info(format("  Improper Torsions:                 %10d", nImproperTorsions));
1051       }
1052     } else {
1053       nImproperTorsions = 0;
1054       improperTorsions = null;
1055     }
1056 
1057     int[] molecule = molecularAssembly.getMoleculeNumbers();
1058     if (vanderWaalsTerm) {
1059       logger.info("\n Non-Bonded Terms");
1060       if (!tornadoVM) {
1061         vanderWaals = new VanDerWaals(
1062             atoms, molecule, crystal, forceField, parallelTeam, vdwCutoff, nlistCutoff);
1063       } else {
1064         vanderWaals = new VanDerWaalsTornado(atoms, crystal, forceField, vdwCutoff);
1065       }
1066     } else {
1067       vanderWaals = null;
1068     }
1069 
1070     if (multipoleTerm) {
1071       ELEC_FORM form;
1072       if (name.contains("OPLS") || name.contains("AMBER") || name.contains("CHARMM")) {
1073         form = ELEC_FORM.FIXED_CHARGE;
1074       } else {
1075         form = ELEC_FORM.PAM;
1076       }
1077 
1078       particleMeshEwald =
1079           new ParticleMeshEwald(
1080               atoms,
1081               molecule,
1082               forceField,
1083               crystal,
1084               vanderWaals.getNeighborList(),
1085               form,
1086               ewaldCutoff,
1087               gkCutoff,
1088               parallelTeam);
1089       double charge = molecularAssembly.getCharge(checkAllNodeCharges);
1090       logger.info(format("\n  Overall system charge:             %10.3f", charge));
1091     } else {
1092       particleMeshEwald = null;
1093     }
1094 
1095     if (ncsTerm) {
1096       String sg = forceField.getString("NCSGROUP", "P 1");
1097       Crystal ncsCrystal = new Crystal(a, b, c, alpha, beta, gamma, sg);
1098       ncsRestraint = new NCSRestraint(atoms, forceField, ncsCrystal);
1099     } else {
1100       ncsRestraint = null;
1101     }
1102 
1103     coordRestraints = new ArrayList<>();
1104     if (restraints != null) {
1105       coordRestraints.addAll(restraints);
1106     }
1107 
1108     if (restrainTerm) {
1109       this.autoCoordRestraint = new CoordRestraint(atoms, forceField);
1110       coordRestraints.add(autoCoordRestraint);
1111     } else {
1112       autoCoordRestraint = null;
1113     }
1114 
1115     if (!coordRestraints.isEmpty()) {
1116       restrainTerm = true;
1117       restrainTermOrig = restrainTerm;
1118       logger.log(Level.FINE, " restrainTerm set true");
1119     }
1120 
1121     if (comTerm) {
1122       comRestraint = new COMRestraint(molecularAssembly);
1123     } else {
1124       comRestraint = null;
1125     }
1126 
1127     if (restrainGroupTerm) {
1128       restrainGroups = new RestrainGroups(molecularAssembly);
1129       nRestrainGroups = restrainGroups.getNumberOfRestraints();
1130     } else {
1131       restrainGroups = null;
1132     }
1133 
1134     String[] restrainTorsionCos = properties.getStringArray("restrain-torsion-cos");
1135     double torsUnits = forceField.getDouble("TORSIONUNIT", TorsionType.DEFAULT_TORSION_UNIT);
1136     List<RestraintTorsion> rTorsList = new ArrayList<>(restrainTorsionCos.length);
1137     for (String rtc : restrainTorsionCos) {
1138       String[] toks = rtc.split("\\s+");
1139       int nTok = toks.length;
1140       if (nTok < 7) {
1141         throw new IllegalArgumentException(format(
1142             "restrain-torsion-cos record %s must have 4 atom indices, amplitude, phase shift and periodicity!",
1143             rtc));
1144       }
1145       int ai1 = Integer.parseInt(toks[0]) - 1;
1146       int ai2 = Integer.parseInt(toks[1]) - 1;
1147       int ai3 = Integer.parseInt(toks[2]) - 1;
1148       int ai4 = Integer.parseInt(toks[3]) - 1;
1149       Atom a1 = atoms[ai1];
1150       Atom a2 = atoms[ai2];
1151       Atom a3 = atoms[ai3];
1152       Atom a4 = atoms[ai4];
1153       int[] atomClasses = new int[] {-1, -1, -1, -1};
1154 
1155       int startTerms = 4;
1156       boolean lamEnabled = false;
1157       boolean revLam = false;
1158       if (toks[4].equalsIgnoreCase("lambda")) {
1159         ++startTerms;
1160         lamEnabled = true;
1161         if (toks[5].toLowerCase().startsWith("reverse")) {
1162           ++startTerms;
1163           revLam = true;
1164         }
1165       }
1166       if (!lambdaTerm) {
1167         lamEnabled = false;
1168       }
1169 
1170       int nTerms = (nTok - startTerms) / 3;
1171       assert (nTok - startTerms) % 3 == 0;
1172       double[] amp = new double[nTerms];
1173       double[] phase = new double[nTerms];
1174       int[] period = new int[nTerms];
1175 
1176       for (int i = 0; i < nTerms; i++) {
1177         int i0 = startTerms + (3 * i);
1178         amp[i] = Double.parseDouble(toks[i0]);
1179         phase[i] = Double.parseDouble(toks[i0 + 1]);
1180         period[i] = Integer.parseInt(toks[i0 + 2]);
1181       }
1182 
1183       // Lambda-enabled torsion restraints require lambda term to be true.
1184       assert !lamEnabled || lambdaTerm;
1185       TorsionType tType = new TorsionType(atomClasses, amp, phase, period);
1186       rTorsList.add(new RestraintTorsion(a1, a2, a3, a4, tType, lamEnabled, revLam, torsUnits));
1187     }
1188 
1189     nRestTors = rTorsList.size();
1190     if (nRestTors > 0) {
1191       logger.info(format(" Adding %4d cosine-based torsion restraints.", nRestTors));
1192       rTorsTerm = true;
1193       rTorsTermOrig = true;
1194       rTors = rTorsList.toArray(new RestraintTorsion[0]);
1195       rTorsEnergy = 0;
1196       rTorsTime = 0;
1197     } else {
1198       rTorsTermOrig = false;
1199     }
1200 
1201     bondedRegion = new BondedRegion();
1202 
1203     maxDebugGradient = forceField.getDouble("MAX_DEBUG_GRADIENT", Double.POSITIVE_INFINITY);
1204 
1205     molecularAssembly.setPotential(this);
1206 
1207     // Add restrain-bond records. If no restrain-distance records exist, the empty array will be
1208     // returned.
1209     String[] bondRestraints = molecularAssembly.getProperties().getStringArray("restrain-distance");
1210     for (String bondRest : bondRestraints) {
1211       try {
1212         String[] toks = bondRest.split("\\s+");
1213         if (toks.length < 2) {
1214           throw new IllegalArgumentException(
1215               format(" restrain-distance value %s could not be parsed!", bondRest));
1216         }
1217         // Internally, everything starts with 0, but restrain distance starts at 1, so that 1 has to
1218         // be subtracted
1219         int at1 = Integer.parseInt(toks[0]) - 1;
1220         int at2 = Integer.parseInt(toks[1]) - 1;
1221 
1222         double forceConst = 100.0;
1223         double flatBottomRadius = 0;
1224         Atom a1 = atoms[at1];
1225         Atom a2 = atoms[at2];
1226 
1227         if (toks.length > 2) {
1228           forceConst = Double.parseDouble(toks[2]);
1229         }
1230         double dist;
1231         switch (toks.length) {
1232           case 2:
1233           case 3:
1234             double[] xyz1 = new double[3];
1235             xyz1 = a1.getXYZ(xyz1);
1236             double[] xyz2 = new double[3];
1237             xyz2 = a2.getXYZ(xyz2);
1238             // Current distance between restrained atoms
1239             dist = crystal.minDistOverSymOps(xyz1, xyz2);
1240             break;
1241           case 4:
1242             dist = Double.parseDouble(toks[3]);
1243             break;
1244           case 5:
1245           default:
1246             double minDist = Double.parseDouble(toks[3]);
1247             double maxDist = Double.parseDouble(toks[4]);
1248             dist = 0.5 * (minDist + maxDist);
1249             flatBottomRadius = 0.5 * Math.abs(maxDist - minDist);
1250             break;
1251         }
1252 
1253         UnivariateSwitchingFunction switchF;
1254         double lamStart = RestraintBond.DEFAULT_RB_LAM_START;
1255         double lamEnd = RestraintBond.DEFAULT_RB_LAM_END;
1256         if (toks.length > 5) {
1257           int offset = 5;
1258           if (toks[5].matches("^[01](?:\\.[0-9]*)?")) {
1259             offset = 6;
1260             lamStart = Double.parseDouble(toks[5]);
1261             if (toks[6].matches("^[01](?:\\.[0-9]*)?")) {
1262               offset = 7;
1263               lamEnd = Double.parseDouble(toks[6]);
1264             }
1265           }
1266           switchF = UnivariateFunctionFactory.parseUSF(toks, offset);
1267         } else {
1268           switchF = new ConstantSwitch();
1269         }
1270 
1271         setRestraintBond(a1, a2, dist, forceConst, flatBottomRadius, lamStart, lamEnd, switchF);
1272       } catch (Exception ex) {
1273         logger.info(format(" Exception in parsing restrain-distance: %s", ex.toString()));
1274       }
1275     }
1276 
1277     String constraintStrings =
1278         forceField.getString("CONSTRAIN", forceField.getString("RATTLE", null));
1279     if (constraintStrings != null) {
1280       constraints = new ArrayList<>();
1281       logger.info(format(" Experimental: parsing constraints option %s", constraintStrings));
1282 
1283       Set<Bond> numericBonds = new HashSet<>(1);
1284       Set<Angle> numericAngles = new HashSet<>(1);
1285 
1286       // Totally empty constrain option: constrain only X-H bonds. No other options applied.
1287       if (constraintStrings.isEmpty() || constraintStrings.matches("^\\s*$")) {
1288         // Assume constraining only X-H bonds (i.e. RIGID-HYDROGEN).
1289         logger.info(" Constraining X-H bonds.");
1290         numericBonds =
1291             Arrays.stream(bonds)
1292                 .filter(
1293                     (Bond bond) ->
1294                         bond.getAtom(0).getAtomicNumber() == 1
1295                             || bond.getAtom(1).getAtomicNumber() == 1)
1296                 .collect(Collectors.toSet());
1297       } else {
1298         String[] constraintToks = constraintStrings.split("\\s+");
1299 
1300         // First, accumulate SETTLE constraints.
1301         for (String tok : constraintToks) {
1302           if (tok.equalsIgnoreCase("WATER")) {
1303             logger.info(" Constraining waters to be rigid based on angle & bonds.");
1304             // XYZ files, in particular, have waters mislabeled as generic Molecules.
1305             // First, find any such mislabeled water.
1306             Stream<MSNode> settleStream =
1307                 molecularAssembly.getMolecules().stream()
1308                     .filter((MSNode m) -> m.getAtomList().size() == 3)
1309                     .filter(
1310                         (MSNode m) -> {
1311                           List<Atom> atoms = m.getAtomList();
1312                           Atom O = null;
1313                           List<Atom> H = new ArrayList<>(2);
1314                           for (Atom at : atoms) {
1315                             int atN = at.getAtomicNumber();
1316                             if (atN == 8) {
1317                               O = at;
1318                             } else if (atN == 1) {
1319                               H.add(at);
1320                             }
1321                           }
1322                           return O != null && H.size() == 2;
1323                         });
1324             // Now concatenate the stream with the properly labeled waters.
1325             settleStream = Stream.concat(settleStream, molecularAssembly.getWater().stream());
1326             // Map them into new Settle constraints and collect.
1327             List<SettleConstraint> settleConstraints =
1328                 settleStream
1329                     .map((MSNode m) -> m.getAngleList().get(0))
1330                     .map(SettleConstraint::settleFactory)
1331                     .collect(Collectors.toList());
1332             constraints.addAll(settleConstraints);
1333 
1334           } else if (tok.equalsIgnoreCase("DIATOMIC")) {
1335             logger.severe(" Diatomic distance constraints not yet implemented properly.");
1336           } else if (tok.equalsIgnoreCase("TRIATOMIC")) {
1337             logger.severe(
1338                 " Triatomic SETTLE constraints for non-water molecules not yet implemented properly.");
1339           }
1340         }
1341 
1342         // Second, accumulate bond/angle constraints.
1343         for (String tok : constraintToks) {
1344           if (tok.equalsIgnoreCase("BONDS")) {
1345             numericBonds = new HashSet<>(Arrays.asList(bonds));
1346           } else if (tok.equalsIgnoreCase("ANGLES")) {
1347             numericAngles = new HashSet<>(Arrays.asList(angles));
1348           }
1349         }
1350       }
1351 
1352       // Remove bonds that are already dealt with via angles.
1353       for (Angle angle : numericAngles) {
1354         numericBonds.removeAll(angle.getBondList());
1355       }
1356 
1357       // Remove already-constrained angles and bonds (e.g. SETTLE-constrained ones).
1358       List<Angle> ccmaAngles =
1359           numericAngles.stream()
1360               .filter((Angle ang) -> !ang.isConstrained())
1361               .collect(Collectors.toList());
1362       List<Bond> ccmaBonds =
1363           numericBonds.stream()
1364               .filter((Bond bond) -> !bond.isConstrained())
1365               .collect(Collectors.toList());
1366 
1367       CcmaConstraint ccmaConstraint =
1368           CcmaConstraint.ccmaFactory(
1369               ccmaBonds, ccmaAngles, atoms, getMass(), CcmaConstraint.DEFAULT_CCMA_NONZERO_CUTOFF);
1370       constraints.add(ccmaConstraint);
1371 
1372       logger.info(format(" Added %d constraints.", constraints.size()));
1373     } else {
1374       constraints = Collections.emptyList();
1375     }
1376 
1377     if (lambdaTerm) {
1378       this.setLambda(1.0);
1379     }
1380   }
1381 
1382   /**
1383    * Static factory method to create a ForceFieldEnergy, possibly via FFX or OpenMM implementations.
1384    *
1385    * @param assembly To create FFE over
1386    * @return a {@link ffx.potential.ForceFieldEnergy} object.
1387    */
1388   public static ForceFieldEnergy energyFactory(MolecularAssembly assembly) {
1389     return energyFactory(assembly, null);
1390   }
1391 
1392   /**
1393    * Static factory method to create a ForceFieldEnergy, possibly via FFX or OpenMM implementations.
1394    *
1395    * @param assembly To create FFE over
1396    * @param restraints Harmonic restraints
1397    * @return a {@link ffx.potential.ForceFieldEnergy} object.
1398    */
1399   public static ForceFieldEnergy energyFactory(
1400       MolecularAssembly assembly, List<CoordRestraint> restraints) {
1401     return energyFactory(assembly, restraints, ParallelTeam.getDefaultThreadCount());
1402   }
1403 
1404   /**
1405    * Static factory method to create a ForceFieldEnergy, possibly via FFX or OpenMM implementations.
1406    *
1407    * @param assembly To create FFE over
1408    * @param restraints Harmonic restraints
1409    * @param numThreads Number of threads to use for FFX energy
1410    * @return A ForceFieldEnergy on some Platform
1411    */
1412   public static ForceFieldEnergy energyFactory(
1413       MolecularAssembly assembly, List<CoordRestraint> restraints, int numThreads) {
1414     ForceField ffield = assembly.getForceField();
1415     String platformString = toEnumForm(ffield.getString("PLATFORM", "FFX"));
1416     try {
1417       Platform platform = Platform.valueOf(platformString);
1418       switch (platform) {
1419         case OMM:
1420         case OMM_REF: // Should be split from the code once we figure out how to specify a kernel.
1421         case OMM_CUDA:
1422           try {
1423             return new ForceFieldEnergyOpenMM(assembly, platform, restraints, numThreads);
1424           } catch (Exception ex) {
1425             logger.warning(format(" Exception creating ForceFieldEnergyOpenMM: %s", ex));
1426             ex.printStackTrace();
1427 
1428             ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1429             if (ffxEnergy == null) {
1430               ffxEnergy = new ForceFieldEnergy(assembly, restraints, numThreads);
1431               assembly.setPotential(ffxEnergy);
1432             }
1433             return ffxEnergy;
1434           }
1435         case OMM_OPENCL:
1436         case OMM_OPTCPU:
1437           logger.warning(format(" Platform %s not supported; defaulting to FFX", platform));
1438         case FFX:
1439         default:
1440           ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1441           if (ffxEnergy == null) {
1442             ffxEnergy = new ForceFieldEnergy(assembly, restraints, numThreads);
1443             assembly.setPotential(ffxEnergy);
1444           }
1445           return ffxEnergy;
1446       }
1447     } catch (IllegalArgumentException | NullPointerException ex) {
1448       logger.warning(
1449           format(" String %s did not match a known energy implementation", platformString));
1450       ForceFieldEnergy ffxEnergy = assembly.getPotentialEnergy();
1451       if (ffxEnergy == null) {
1452         ffxEnergy = new ForceFieldEnergy(assembly, restraints, numThreads);
1453         assembly.setPotential(ffxEnergy);
1454       }
1455       return ffxEnergy;
1456     }
1457   }
1458 
1459   /**
1460    * Applies constraints to positions
1461    *
1462    * @param xPrior Prior coodinates.
1463    * @param xNew New coordinates.
1464    */
1465   public void applyAllConstraintPositions(double[] xPrior, double[] xNew) {
1466     applyAllConstraintPositions(xPrior, xNew, DEFAULT_CONSTRAINT_TOLERANCE);
1467   }
1468 
1469   /**
1470    * Applies constraints to positions
1471    *
1472    * @param xPrior Prior coodinates.
1473    * @param xNew New coordinates.
1474    * @param tol Constraint Tolerance.
1475    */
1476   public void applyAllConstraintPositions(double[] xPrior, double[] xNew, double tol) {
1477     if (xPrior == null) {
1478       xPrior = Arrays.copyOf(xNew, xNew.length);
1479     }
1480     for (Constraint constraint : constraints) {
1481       constraint.applyConstraintToStep(xPrior, xNew, getMass(), tol);
1482     }
1483   }
1484 
1485   /**
1486    * Overwrites current esvSystem if present. Multiple ExtendedSystems is possible but unnecessary;
1487    * add all ESVs to one system (per FFE, at least).
1488    *
1489    * @param system a {@link ffx.potential.extended.ExtendedSystem} object.
1490    */
1491   public void attachExtendedSystem(ExtendedSystem system) {
1492     if (system == null) {
1493       throw new IllegalArgumentException();
1494     }
1495     esvTerm = true;
1496     esvSystem = system;
1497     if (vanderWaalsTerm) {
1498       if (vanderWaals == null) {
1499         logger.warning("Null VdW during ESV setup.");
1500       } else {
1501         vanderWaals.attachExtendedSystem(system);
1502       }
1503     }
1504     if (multipoleTerm) {
1505       if (particleMeshEwald == null) {
1506         logger.warning("Null PME during ESV setup.");
1507       }
1508       particleMeshEwald.attachExtendedSystem(system);
1509     }
1510     if (crystal != null) {
1511       crystal.setSpecialPositionCutoff(0.0);
1512     }
1513     // Check if reInit is needed, now that we init with fully protonated residues ASD/GLD/HIS
1514     reInit();
1515   }
1516 
1517   /**
1518    * {@inheritDoc}
1519    *
1520    * <p>Returns true if lambda term is not enabled for this ForceFieldEnergy.
1521    */
1522   @Override
1523   public boolean dEdLZeroAtEnds() {
1524     // This may actually be true even with softcored atoms.
1525     // For now, serves the purpose of reporting true when nothing is softcored.
1526     return !lambdaTerm;
1527   }
1528 
1529   /**
1530    * Frees up assets associated with this ForceFieldEnergy, such as worker Threads.
1531    *
1532    * @return If successful in freeing up assets.
1533    */
1534   public boolean destroy() {
1535     if (destroyed) {
1536       // This regularly occurs with Repex OST, as multiple OrthogonalSpaceTempering objects wrap a
1537       // single FFE.
1538       logger.fine(format(" This ForceFieldEnergy is already destroyed: %s", this.toString()));
1539       return true;
1540     } else {
1541       try {
1542         if (parallelTeam != null) {
1543           parallelTeam.shutdown();
1544         }
1545         if (vanderWaals != null) {
1546           vanderWaals.destroy();
1547         }
1548         if (particleMeshEwald != null) {
1549           particleMeshEwald.destroy();
1550         }
1551         molecularAssembly.finishDestruction();
1552         destroyed = true;
1553         return true;
1554       } catch (Exception ex) {
1555         logger.warning(format(" Exception in shutting down a ForceFieldEnergy: %s", ex));
1556         logger.info(Utilities.stackTraceToString(ex));
1557         return false;
1558       }
1559     }
1560   }
1561 
1562   /**
1563    * energy.
1564    *
1565    * @return a double.
1566    */
1567   public double energy() {
1568     return energy(false, false);
1569   }
1570 
1571   /**
1572    * {@inheritDoc}
1573    *
1574    * <p>energy
1575    */
1576   public double energy(boolean gradient, boolean print) {
1577 
1578     try {
1579       bondTime = 0;
1580       angleTime = 0;
1581       stretchBendTime = 0;
1582       ureyBradleyTime = 0;
1583       outOfPlaneBendTime = 0;
1584       torsionTime = 0;
1585       stretchTorsionTime = 0;
1586       angleTorsionTime = 0;
1587       piOrbitalTorsionTime = 0;
1588       torsionTorsionTime = 0;
1589       improperTorsionTime = 0;
1590       vanDerWaalsTime = 0;
1591       electrostaticTime = 0;
1592       restraintBondTime = 0;
1593       ncsTime = 0;
1594       coordRestraintTime = 0;
1595       rTorsTime = 0;
1596       totalTime = System.nanoTime();
1597 
1598       // Zero out the potential energy of each bonded term.
1599       bondEnergy = 0.0;
1600       angleEnergy = 0.0;
1601       stretchBendEnergy = 0.0;
1602       ureyBradleyEnergy = 0.0;
1603       outOfPlaneBendEnergy = 0.0;
1604       torsionEnergy = 0.0;
1605       angleTorsionEnergy = 0.0;
1606       stretchTorsionEnergy = 0.0;
1607       piOrbitalTorsionEnergy = 0.0;
1608       torsionTorsionEnergy = 0.0;
1609       improperTorsionEnergy = 0.0;
1610       totalBondedEnergy = 0.0;
1611 
1612       // Zero out potential energy of restraint terms
1613       restraintBondEnergy = 0.0;
1614       ncsEnergy = 0.0;
1615       rTorsEnergy = 0.0;
1616       restrainEnergy = 0.0;
1617 
1618       // Zero out bond and angle RMSDs.
1619       bondRMSD = 0.0;
1620       angleRMSD = 0.0;
1621 
1622       // Zero out the potential energy of each non-bonded term.
1623       vanDerWaalsEnergy = 0.0;
1624       permanentMultipoleEnergy = 0.0;
1625       permanentRealSpaceEnergy = 0.0;
1626       polarizationEnergy = 0.0;
1627       totalMultipoleEnergy = 0.0;
1628       totalNonBondedEnergy = 0.0;
1629 
1630       // Zero out the solvation energy.
1631       solvationEnergy = 0.0;
1632 
1633       // Zero out the relative solvation energy (sequence optimization)
1634       relativeSolvationEnergy = 0.0;
1635       nRelativeSolvations = 0;
1636 
1637       esvBias = 0.0;
1638 
1639       // Zero out the total potential energy.
1640       totalEnergy = 0.0;
1641 
1642       // Zero out the Cartesian coordinate gradient for each atom.
1643       //            if (gradient) {
1644       //                for (int i = 0; i < nAtoms; i++) {
1645       //                    atoms[i].setXYZGradient(0.0, 0.0, 0.0);
1646       //                    atoms[i].setLambdaXYZGradient(0.0, 0.0, 0.0);
1647       //                }
1648       //            }
1649 
1650       // Computed the bonded energy terms in parallel.
1651       try {
1652         bondedRegion.setGradient(gradient);
1653         parallelTeam.execute(bondedRegion);
1654       } catch (RuntimeException ex) {
1655         logger.warning("Runtime exception during bonded term calculation.");
1656         throw ex;
1657       } catch (Exception ex) {
1658         logger.info(Utilities.stackTraceToString(ex));
1659         logger.severe(ex.toString());
1660       }
1661 
1662       if (!lambdaBondedTerms) {
1663         // Compute restraint terms.
1664         if (ncsTerm) {
1665           ncsTime = -System.nanoTime();
1666           ncsEnergy = ncsRestraint.residual(gradient, print);
1667           ncsTime += System.nanoTime();
1668         }
1669         if (restrainTerm && !coordRestraints.isEmpty()) {
1670           coordRestraintTime = -System.nanoTime();
1671           for (CoordRestraint restraint : coordRestraints) {
1672             restrainEnergy += restraint.residual(gradient, print);
1673           }
1674           coordRestraintTime += System.nanoTime();
1675         }
1676         if (comTerm) {
1677           comRestraintTime = -System.nanoTime();
1678           comRestraintEnergy = comRestraint.residual(gradient, print);
1679           comRestraintTime += System.nanoTime();
1680         }
1681         if (restrainGroupTerm) {
1682           restrainGroupTime = -System.nanoTime();
1683           restrainGroupEnergy = restrainGroups.energy(gradient);
1684           restrainGroupTime += System.nanoTime();
1685         }
1686         // Compute non-bonded terms.
1687         if (vanderWaalsTerm) {
1688           vanDerWaalsTime = -System.nanoTime();
1689           vanDerWaalsEnergy = vanderWaals.energy(gradient, print);
1690           nVanDerWaalInteractions = this.vanderWaals.getInteractions();
1691           vanDerWaalsTime += System.nanoTime();
1692         }
1693         if (multipoleTerm) {
1694           electrostaticTime = -System.nanoTime();
1695           totalMultipoleEnergy = particleMeshEwald.energy(gradient, print);
1696           permanentMultipoleEnergy = particleMeshEwald.getPermanentEnergy();
1697           permanentRealSpaceEnergy = particleMeshEwald.getPermRealEnergy();
1698           polarizationEnergy = particleMeshEwald.getPolarizationEnergy();
1699           nPermanentInteractions = particleMeshEwald.getInteractions();
1700           solvationEnergy = particleMeshEwald.getSolvationEnergy();
1701           nGKInteractions = particleMeshEwald.getGKInteractions();
1702           electrostaticTime += System.nanoTime();
1703         }
1704       }
1705 
1706       if (relativeSolvationTerm) {
1707         List<Residue> residuesList = molecularAssembly.getResidueList();
1708         for (Residue residue : residuesList) {
1709           if (residue instanceof MultiResidue) {
1710             Atom refAtom = residue.getSideChainAtoms().get(0);
1711             if (refAtom != null && refAtom.getUse()) {
1712               // Reasonably confident that it should be -=,
1713               // as we are trying to penalize residues with strong solvation energy.
1714               double thisSolvation = relativeSolvation.getSolvationEnergy(residue, false);
1715               relativeSolvationEnergy -= thisSolvation;
1716               if (thisSolvation != 0) {
1717                 nRelativeSolvations++;
1718               }
1719             }
1720           }
1721         }
1722       }
1723 
1724       totalTime = System.nanoTime() - totalTime;
1725 
1726       totalBondedEnergy =
1727           bondEnergy
1728               + restraintBondEnergy
1729               + angleEnergy
1730               + stretchBendEnergy
1731               + ureyBradleyEnergy
1732               + outOfPlaneBendEnergy
1733               + torsionEnergy
1734               + angleTorsionEnergy
1735               + stretchTorsionEnergy
1736               + piOrbitalTorsionEnergy
1737               + improperTorsionEnergy
1738               + torsionTorsionEnergy
1739               + ncsEnergy
1740               + restrainEnergy
1741               + restrainGroupEnergy
1742               + rTorsEnergy;
1743       totalNonBondedEnergy = vanDerWaalsEnergy + totalMultipoleEnergy + relativeSolvationEnergy;
1744       totalEnergy = totalBondedEnergy + totalNonBondedEnergy + solvationEnergy;
1745       if (esvTerm) {
1746         esvBias = esvSystem.getBiasEnergy();
1747         totalEnergy += esvBias;
1748       }
1749 
1750       if (print) {
1751         StringBuilder sb = new StringBuilder();
1752         if (gradient) {
1753           sb.append("\n Computed Potential Energy and Atomic Coordinate Gradients\n");
1754         } else {
1755           sb.append("\n Computed Potential Energy\n");
1756         }
1757         sb.append(this);
1758         logger.info(sb.toString());
1759       }
1760       return totalEnergy;
1761     } catch (EnergyException ex) {
1762       if (printOnFailure) {
1763         printFailure();
1764       }
1765       if (ex.doCauseSevere()) {
1766         logger.info(Utilities.stackTraceToString(ex));
1767         logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
1768       } else {
1769         logger.log(Level.INFO, format(" Exception in energy calculation:\n %s", ex));
1770       }
1771       throw ex;
1772     }
1773   }
1774 
1775   /** {@inheritDoc} */
1776   @Override
1777   public double energy(double[] x) {
1778     return energy(x, false);
1779   }
1780 
1781   /** {@inheritDoc} */
1782   @Override
1783   public double energy(double[] x, boolean verbose) {
1784     assert Arrays.stream(x).allMatch(Double::isFinite);
1785 
1786     // Unscale the coordinates.
1787     unscaleCoordinates(x);
1788 
1789     // Set coordinates.
1790     setCoordinates(x);
1791 
1792     double e = this.energy(false, verbose);
1793 
1794     // Rescale the coordinates.
1795     scaleCoordinates(x);
1796 
1797     return e;
1798   }
1799 
1800   /** {@inheritDoc} */
1801   @Override
1802   public double energyAndGradient(double[] x, double[] g) {
1803     return energyAndGradient(x, g, false);
1804   }
1805 
1806   /** {@inheritDoc} */
1807   @Override
1808   public double energyAndGradient(double[] x, double[] g, boolean verbose) {
1809     assert Arrays.stream(x).allMatch(Double::isFinite);
1810 
1811     // Un-scale the coordinates.
1812     unscaleCoordinates(x);
1813 
1814     // Set coordinates.
1815     setCoordinates(x);
1816     double e = energy(true, verbose);
1817 
1818     // Try block already exists inside energy(boolean, boolean), so only
1819     // need to try-catch fillGradient.
1820     try {
1821       fillGradient(g);
1822 
1823       // Scale the coordinates and gradients.
1824       scaleCoordinatesAndGradient(x, g);
1825 
1826       if (maxDebugGradient < Double.MAX_VALUE) {
1827         boolean extremeGrad =
1828             Arrays.stream(g)
1829                 .anyMatch((double gi) -> (gi > maxDebugGradient || gi < -maxDebugGradient));
1830         if (extremeGrad) {
1831           File origFile = molecularAssembly.getFile();
1832           String timeString =
1833               LocalDateTime.now().format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
1834 
1835           String filename =
1836               format(
1837                   "%s-LARGEGRAD-%s.pdb",
1838                   removeExtension(molecularAssembly.getFile().getName()), timeString);
1839           PotentialsFunctions ef = new PotentialsUtils();
1840           filename = ef.versionFile(filename);
1841 
1842           logger.warning(
1843               format(
1844                   " Excessively large gradient detected; printing snapshot to file %s", filename));
1845           ef.saveAsPDB(molecularAssembly, new File(filename));
1846           molecularAssembly.setFile(origFile);
1847         }
1848       }
1849       return e;
1850     } catch (EnergyException ex) {
1851       if (printOnFailure) {
1852         printFailure();
1853       }
1854       if (ex.doCauseSevere()) {
1855         logger.info(Utilities.stackTraceToString(ex));
1856         logger.log(Level.SEVERE, " Error in calculating energies or gradients", ex);
1857       } else {
1858         logger.log(Level.INFO, format(" Exception in energy calculation:\n %s", ex));
1859       }
1860       throw ex;
1861     }
1862   }
1863 
1864   /**
1865    * Save coordinates when an EnergyException is caught.
1866    */
1867   private void printFailure() {
1868     String timeString = LocalDateTime.now()
1869         .format(DateTimeFormatter.ofPattern("yyyy_MM_dd-HH_mm_ss"));
1870     File file = molecularAssembly.getFile();
1871     String ext = "pdb";
1872     if (isXYZ(file)) {
1873       ext = "xyz";
1874     }
1875     String filename = format("%s-ERROR-%s.%s", removeExtension(file.getName()), timeString, ext);
1876     PotentialsFunctions ef = new PotentialsUtils();
1877     filename = ef.versionFile(filename);
1878     logger.info(format(" Writing on-error snapshot to file %s", filename));
1879     ef.save(molecularAssembly, new File(filename));
1880   }
1881 
1882   /**
1883    * {@inheritDoc}
1884    *
1885    * <p>Returns an array of acceleration values for active atoms.
1886    */
1887   @Override
1888   public double[] getAcceleration(double[] acceleration) {
1889     int n = getNumberOfVariables();
1890     if (acceleration == null || acceleration.length < n) {
1891       acceleration = new double[n];
1892     }
1893     int index = 0;
1894     double[] a = new double[3];
1895     for (int i = 0; i < nAtoms; i++) {
1896       if (atoms[i].isActive()) {
1897         atoms[i].getAcceleration(a);
1898         acceleration[index++] = a[0];
1899         acceleration[index++] = a[1];
1900         acceleration[index++] = a[2];
1901       }
1902     }
1903     return acceleration;
1904   }
1905 
1906   /**
1907    * Getter for the field <code>angleEnergy</code>.
1908    *
1909    * @return a double.
1910    */
1911   public double getAngleEnergy() {
1912     return angleEnergy;
1913   }
1914 
1915   /**
1916    * Getter for the field <code>angleTorsionEnergy</code>.
1917    *
1918    * @return a double.
1919    */
1920   public double getAngleTorsionEnergy() {
1921     return angleTorsionEnergy;
1922   }
1923 
1924   /**
1925    * Getter for the field <code>angleTorsions</code>.
1926    *
1927    * @return an array of {@link ffx.potential.bonded.AngleTorsion} objects.
1928    */
1929   public AngleTorsion[] getAngleTorsions() {
1930     return angleTorsions;
1931   }
1932 
1933   /**
1934    * Getter for the field <code>angles</code>. Both normal and in-plane angles are returned.
1935    *
1936    * @return an array of {@link ffx.potential.bonded.Angle} objects.
1937    */
1938   public Angle[] getAngles() {
1939     return angles;
1940   }
1941 
1942   /**
1943    * Getter for the field <code>angles</code> with only <code>AngleMode</code> angles.
1944    *
1945    * @param angleMode Only angles of this mode will be returned.
1946    * @return an array of {@link ffx.potential.bonded.Angle} objects.
1947    */
1948   public Angle[] getAngles(AngleMode angleMode) {
1949     if (angles == null || angles.length < 1) {
1950       return null;
1951     }
1952     int nAngles = angles.length;
1953     List<Angle> angleList = new ArrayList<>();
1954     // Sort all normal angles from in-plane angles
1955     for (int i = 0; i < nAngles; i++) {
1956       if (angles[i].getAngleMode() == angleMode) {
1957         angleList.add(angles[i]);
1958       }
1959     }
1960     nAngles = angleList.size();
1961     if (nAngles < 1) {
1962       return null;
1963     }
1964     return angleList.toArray(new Angle[0]);
1965   }
1966 
1967   /**
1968    * Getter for the field <code>bondEnergy</code>.
1969    *
1970    * @return a double.
1971    */
1972   public double getBondEnergy() {
1973     return bondEnergy;
1974   }
1975 
1976   /**
1977    * Getter for the field <code>bonds</code>.
1978    *
1979    * @return an array of {@link ffx.potential.bonded.Bond} objects.
1980    */
1981   public Bond[] getBonds() {
1982     return bonds;
1983   }
1984 
1985   /**
1986    * getCavitationEnergy.
1987    *
1988    * @return a double.
1989    */
1990   public double getCavitationEnergy() {
1991     return particleMeshEwald.getCavitationEnergy();
1992   }
1993 
1994   /**
1995    * Returns a copy of the list of constraints this ForceFieldEnergy has.
1996    *
1997    * @return Copied list of constraints.
1998    */
1999   @Override
2000   public List<Constraint> getConstraints() {
2001     return constraints.isEmpty() ? Collections.emptyList() : new ArrayList<>(constraints);
2002   }
2003 
2004   /**
2005    * Getter for the field <code>coordRestraints</code>.
2006    *
2007    * @return a {@link java.util.List} object.
2008    */
2009   public List<CoordRestraint> getCoordRestraints() {
2010     return new ArrayList<>(coordRestraints);
2011   }
2012 
2013   /** {@inheritDoc} */
2014   @Override
2015   public double[] getCoordinates(double[] x) {
2016     int n = getNumberOfVariables();
2017     if (x == null || x.length < n) {
2018       x = new double[n];
2019     }
2020     int index = 0;
2021     for (int i = 0; i < nAtoms; i++) {
2022       Atom a = atoms[i];
2023       if (a.isActive() && !a.isBackground()) {
2024         x[index++] = a.getX();
2025         x[index++] = a.getY();
2026         x[index++] = a.getZ();
2027       }
2028     }
2029     return x;
2030   }
2031 
2032   /**
2033    * {@inheritDoc}
2034    *
2035    * <p>Getter for the field <code>crystal</code>.
2036    */
2037   @Override
2038   public Crystal getCrystal() {
2039     return crystal;
2040   }
2041 
2042   /**
2043    * {@inheritDoc}
2044    *
2045    * <p>Set the boundary conditions for this calculation.
2046    */
2047   @Override
2048   public void setCrystal(Crystal crystal) {
2049     setCrystal(crystal, false);
2050   }
2051 
2052   /**
2053    * Getter for the field <code>cutoffPlusBuffer</code>.
2054    *
2055    * @return a double.
2056    */
2057   public double getCutoffPlusBuffer() {
2058     return cutoffPlusBuffer;
2059   }
2060 
2061   /**
2062    * getDispersionEnergy.
2063    *
2064    * @return a double.
2065    */
2066   public double getDispersionEnergy() {
2067     return particleMeshEwald.getDispersionEnergy();
2068   }
2069 
2070   /**
2071    * getElectrostaticEnergy.
2072    *
2073    * @return a double.
2074    */
2075   public double getElectrostaticEnergy() {
2076     return totalMultipoleEnergy;
2077   }
2078 
2079   /**
2080    * getEnergyComponent.
2081    *
2082    * @param component a {@link ffx.potential.PotentialComponent} object.
2083    * @return a double.
2084    */
2085   public double getEnergyComponent(PotentialComponent component) {
2086     switch (component) {
2087       case Topology:
2088       case ForceFieldEnergy:
2089         return totalEnergy;
2090       case VanDerWaals:
2091         return vanDerWaalsEnergy;
2092       case Multipoles:
2093         return particleMeshEwald.getTotalMultipoleEnergy();
2094       case Permanent:
2095         return particleMeshEwald.getPermanentEnergy();
2096       case Induced:
2097         return particleMeshEwald.getPolarizationEnergy();
2098       case PermanentRealSpace:
2099         return particleMeshEwald.getPermRealEnergy();
2100       case PermanentSelf:
2101         return particleMeshEwald.getPermSelfEnergy();
2102       case PermanentReciprocal:
2103         return particleMeshEwald.getPermRecipEnergy();
2104       case InducedRealSpace:
2105         return particleMeshEwald.getIndRealEnergy();
2106       case InducedSelf:
2107         return particleMeshEwald.getIndSelfEnergy();
2108       case InducedReciprocal:
2109         return particleMeshEwald.getIndRecipEnergy();
2110       case GeneralizedKirkwood:
2111         return particleMeshEwald.getSolvationEnergy();
2112       case Bonded:
2113         return totalBondedEnergy;
2114       case Bond:
2115         return bondEnergy;
2116       case Angle:
2117         return angleEnergy;
2118       case Torsion:
2119         return torsionEnergy;
2120       case StretchBend:
2121         return stretchBendEnergy;
2122       case OutOfPlaneBend:
2123         return outOfPlaneBendEnergy;
2124       case PiOrbitalTorsion:
2125         return piOrbitalTorsionEnergy;
2126       case TorsionTorsion:
2127         return torsionTorsionEnergy;
2128       case UreyBradley:
2129         return ureyBradleyEnergy;
2130       case RestraintBond:
2131         return restraintBondEnergy;
2132       case ImproperTorsion:
2133         return improperTorsionEnergy;
2134       case NCS:
2135         return ncsEnergy;
2136       case Restrain:
2137         return restrainEnergy;
2138       case pHMD:
2139       case Bias:
2140       case DiscretizeBias:
2141       case ModelBias:
2142       case pHBias:
2143         return (esvTerm) ? esvSystem.getEnergyComponent(component) : 0.0;
2144       case XRay:
2145       default:
2146         throw new AssertionError(component.name());
2147     }
2148   }
2149 
2150   /** {@inheritDoc} */
2151   @Override
2152   public STATE getEnergyTermState() {
2153     return state;
2154   }
2155 
2156   /**
2157    * {@inheritDoc}
2158    *
2159    * <p>This method is for the RESPA integrator only.
2160    */
2161   @Override
2162   public void setEnergyTermState(STATE state) {
2163     this.state = state;
2164     switch (state) {
2165       case FAST:
2166         bondTerm = bondTermOrig;
2167         angleTerm = angleTermOrig;
2168         stretchBendTerm = stretchBendTermOrig;
2169         ureyBradleyTerm = ureyBradleyTermOrig;
2170         outOfPlaneBendTerm = outOfPlaneBendTermOrig;
2171         torsionTerm = torsionTermOrig;
2172         stretchTorsionTerm = stretchTorsionTermOrig;
2173         angleTorsionTerm = angleTorsionTermOrig;
2174         piOrbitalTorsionTerm = piOrbitalTorsionTermOrig;
2175         torsionTorsionTerm = torsionTorsionTermOrig;
2176         improperTorsionTerm = improperTorsionTermOrig;
2177         restraintBondTerm = restraintBondTermOrig;
2178         ncsTerm = ncsTermOrig;
2179         restrainTerm = restrainTermOrig;
2180         restrainGroupTerm = restrainGroupTermOrig;
2181         comTerm = comTermOrig;
2182         vanderWaalsTerm = false;
2183         multipoleTerm = false;
2184         polarizationTerm = false;
2185         generalizedKirkwoodTerm = false;
2186         break;
2187       case SLOW:
2188         vanderWaalsTerm = vanderWaalsTermOrig;
2189         multipoleTerm = multipoleTermOrig;
2190         polarizationTerm = polarizationTermOrig;
2191         generalizedKirkwoodTerm = generalizedKirkwoodTermOrig;
2192         bondTerm = false;
2193         angleTerm = false;
2194         stretchBendTerm = false;
2195         ureyBradleyTerm = false;
2196         outOfPlaneBendTerm = false;
2197         torsionTerm = false;
2198         stretchTorsionTerm = false;
2199         angleTorsionTerm = false;
2200         piOrbitalTorsionTerm = false;
2201         torsionTorsionTerm = false;
2202         improperTorsionTerm = false;
2203         restraintBondTerm = false;
2204         ncsTerm = false;
2205         restrainTerm = false;
2206         comTerm = false;
2207         restrainGroupTerm = false;
2208         break;
2209       default:
2210         bondTerm = bondTermOrig;
2211         angleTerm = angleTermOrig;
2212         stretchBendTerm = stretchBendTermOrig;
2213         ureyBradleyTerm = ureyBradleyTermOrig;
2214         outOfPlaneBendTerm = outOfPlaneBendTermOrig;
2215         torsionTerm = torsionTermOrig;
2216         stretchTorsionTerm = stretchTorsionTermOrig;
2217         angleTorsionTerm = angleTorsionTermOrig;
2218         piOrbitalTorsionTerm = piOrbitalTorsionTermOrig;
2219         torsionTorsionTerm = torsionTorsionTermOrig;
2220         improperTorsionTerm = improperTorsionTermOrig;
2221         restraintBondTerm = restraintBondTermOrig;
2222         ncsTerm = ncsTermOrig;
2223         restrainTermOrig = restrainTerm;
2224         comTermOrig = comTerm;
2225         restrainGroupTermOrig = restrainGroupTerm;
2226         vanderWaalsTerm = vanderWaalsTermOrig;
2227         multipoleTerm = multipoleTermOrig;
2228         polarizationTerm = polarizationTermOrig;
2229         generalizedKirkwoodTerm = generalizedKirkwoodTermOrig;
2230     }
2231   }
2232 
2233   /**
2234    * getEsvBiasEnergy.
2235    *
2236    * @return a double.
2237    */
2238   public double getEsvBiasEnergy() {
2239     return esvBias;
2240   }
2241 
2242   /**
2243    * getExtendedSystem.
2244    *
2245    * @return a {@link ffx.potential.extended.ExtendedSystem} object.
2246    */
2247   public ExtendedSystem getExtendedSystem() {
2248     return esvSystem;
2249   }
2250 
2251   /**
2252    * getGK.
2253    *
2254    * @return a {@link ffx.potential.nonbonded.GeneralizedKirkwood} object.
2255    */
2256   public GeneralizedKirkwood getGK() {
2257     if (particleMeshEwald != null) {
2258       return particleMeshEwald.getGK();
2259     } else {
2260       return null;
2261     }
2262   }
2263 
2264   /**
2265    * getGKEnergy.
2266    *
2267    * @return a double.
2268    */
2269   public double getGKEnergy() {
2270     return particleMeshEwald.getGKEnergy();
2271   }
2272 
2273   /**
2274    * getGradient
2275    *
2276    * @param g an array of double.
2277    * @return an array of {@link double} objects.
2278    */
2279   public double[] getGradient(double[] g) {
2280     return fillGradient(g);
2281   }
2282 
2283   /**
2284    * Getter for the field <code>improperTorsionEnergy</code>.
2285    *
2286    * @return a double.
2287    */
2288   public double getImproperTorsionEnergy() {
2289     return improperTorsionEnergy;
2290   }
2291 
2292   /**
2293    * Getter for the field <code>improperTorsions</code>.
2294    *
2295    * @return an array of {@link ffx.potential.bonded.ImproperTorsion} objects.
2296    */
2297   public ImproperTorsion[] getImproperTorsions() {
2298     return improperTorsions;
2299   }
2300 
2301   /** {@inheritDoc} */
2302   @Override
2303   public double getLambda() {
2304     return lambda;
2305   }
2306 
2307   /** {@inheritDoc} */
2308   @Override
2309   public void setLambda(double lambda) {
2310     if (lambdaTerm) {
2311       if (lambda <= 1.0 && lambda >= 0.0) {
2312         this.lambda = lambda;
2313         if (vanderWaalsTerm) {
2314           vanderWaals.setLambda(lambda);
2315         }
2316         if (multipoleTerm) {
2317           particleMeshEwald.setLambda(lambda);
2318         }
2319         if (restraintBondTerm && restraintBonds != null) {
2320           for (RestraintBond restraintBond : restraintBonds) {
2321             restraintBond.setLambda(lambda);
2322           }
2323         }
2324         if (ncsTerm && ncsRestraint != null) {
2325           ncsRestraint.setLambda(lambda);
2326         }
2327         if (restrainTerm && !coordRestraints.isEmpty()) {
2328           for (CoordRestraint restraint : coordRestraints) {
2329             restraint.setLambda(lambda);
2330           }
2331         }
2332         if (comTerm && comRestraint != null) {
2333           comRestraint.setLambda(lambda);
2334         }
2335         if (rTorsTerm) {
2336           for (RestraintTorsion rt : rTors) {
2337             rt.setLambda(lambda);
2338           }
2339         }
2340         if (lambdaTorsions) {
2341           for (int i = 0; i < nTorsions; i++) {
2342             torsions[i].setLambda(lambda);
2343           }
2344           for (int i = 0; i < nPiOrbitalTorsions; i++) {
2345             piOrbitalTorsions[i].setLambda(lambda);
2346           }
2347           for (int i = 0; i < nTorsionTorsions; i++) {
2348             torsionTorsions[i].setLambda(lambda);
2349           }
2350         }
2351       } else {
2352         String message = format("Lambda value %8.3f is not in the range [0..1].", lambda);
2353         logger.warning(message);
2354       }
2355     } else {
2356       logger.fine(" Attempting to set a lambda value on a ForceFieldEnergy with lambdaterm false.");
2357     }
2358   }
2359 
2360   /** {@inheritDoc} */
2361   @Override
2362   public double[] getMass() {
2363     int n = getNumberOfVariables();
2364     double[] mass = new double[n];
2365     int index = 0;
2366     for (int i = 0; i < nAtoms; i++) {
2367       Atom a = atoms[i];
2368       if (a.isActive()) {
2369         double m = a.getMass();
2370         mass[index++] = m;
2371         mass[index++] = m;
2372         mass[index++] = m;
2373       }
2374     }
2375     return mass;
2376   }
2377 
2378   /** {@inheritDoc} */
2379   @Override
2380   public int getNumberOfVariables() {
2381     int nActive = 0;
2382     for (int i = 0; i < nAtoms; i++) {
2383       if (atoms[i].isActive() && !atoms[i].isBackground()) {
2384         nActive++;
2385       }
2386     }
2387     return nActive * 3;
2388   }
2389 
2390   /**
2391    * getNumberofAngleTorsions.
2392    *
2393    * @return a int.
2394    */
2395   public int getNumberofAngleTorsions() {
2396     return nAngleTorsions;
2397   }
2398 
2399   /**
2400    * getNumberofAngles.
2401    *
2402    * @return a int.
2403    */
2404   public int getNumberofAngles() {
2405     return nAngles;
2406   }
2407 
2408   /**
2409    * getNumberofBonds.
2410    *
2411    * @return a int.
2412    */
2413   public int getNumberofBonds() {
2414     return nBonds;
2415   }
2416 
2417   /**
2418    * getNumberofImproperTorsions.
2419    *
2420    * @return a int.
2421    */
2422   public int getNumberofImproperTorsions() {
2423     return nImproperTorsions;
2424   }
2425 
2426   /**
2427    * getNumberofOutOfPlaneBends.
2428    *
2429    * @return a int.
2430    */
2431   public int getNumberofOutOfPlaneBends() {
2432     return nOutOfPlaneBends;
2433   }
2434 
2435   /**
2436    * getNumberofPiOrbitalTorsions.
2437    *
2438    * @return a int.
2439    */
2440   public int getNumberofPiOrbitalTorsions() {
2441     return nPiOrbitalTorsions;
2442   }
2443 
2444   /**
2445    * getNumberofStretchBends.
2446    *
2447    * @return a int.
2448    */
2449   public int getNumberofStretchBends() {
2450     return nStretchBends;
2451   }
2452 
2453   /**
2454    * getNumberofStretchTorsions.
2455    *
2456    * @return a int.
2457    */
2458   public int getNumberofStretchTorsions() {
2459     return nStretchTorsions;
2460   }
2461 
2462   /**
2463    * getNumberofTorsionTorsions.
2464    *
2465    * @return a int.
2466    */
2467   public int getNumberofTorsionTorsions() {
2468     return nTorsionTorsions;
2469   }
2470 
2471   /**
2472    * getNumberofTorsions.
2473    *
2474    * @return a int.
2475    */
2476   public int getNumberofTorsions() {
2477     return nTorsions;
2478   }
2479 
2480   /**
2481    * getNumberofUreyBradleys.
2482    *
2483    * @return a int.
2484    */
2485   public int getNumberofUreyBradleys() {
2486     return nUreyBradleys;
2487   }
2488 
2489   /**
2490    * Getter for the field <code>outOfPlaneBendEnergy</code>.
2491    *
2492    * @return a double.
2493    */
2494   public double getOutOfPlaneBendEnergy() {
2495     return outOfPlaneBendEnergy;
2496   }
2497 
2498   /**
2499    * Getter for the field <code>outOfPlaneBends</code>.
2500    *
2501    * @return an array of {@link ffx.potential.bonded.OutOfPlaneBend} objects.
2502    */
2503   public OutOfPlaneBend[] getOutOfPlaneBends() {
2504     return outOfPlaneBends;
2505   }
2506 
2507   /**
2508    * getPDBHeaderString
2509    *
2510    * @return a {@link java.lang.String} object.
2511    */
2512   public String getPDBHeaderString() {
2513     energy(false, false);
2514     StringBuilder sb = new StringBuilder();
2515     sb.append("REMARK   3  CALCULATED POTENTIAL ENERGY\n");
2516     if (bondTerm) {
2517       sb.append(
2518           format("REMARK   3   %s %g (%d)\n", "BOND STRETCHING            : ", bondEnergy, nBonds));
2519       sb.append(format("REMARK   3   %s %g\n", "BOND RMSD                  : ", bondRMSD));
2520     }
2521     if (angleTerm) {
2522       sb.append(
2523           format(
2524               "REMARK   3   %s %g (%d)\n", "ANGLE BENDING              : ", angleEnergy, nAngles));
2525       sb.append(format("REMARK   3   %s %g\n", "ANGLE RMSD                 : ", angleRMSD));
2526     }
2527     if (stretchBendTerm) {
2528       sb.append(
2529           format(
2530               "REMARK   3   %s %g (%d)\n",
2531               "STRETCH-BEND               : ", stretchBendEnergy, nStretchBends));
2532     }
2533     if (ureyBradleyTerm) {
2534       sb.append(
2535           format(
2536               "REMARK   3   %s %g (%d)\n",
2537               "UREY-BRADLEY               : ", ureyBradleyEnergy, nUreyBradleys));
2538     }
2539     if (outOfPlaneBendTerm) {
2540       sb.append(
2541           format(
2542               "REMARK   3   %s %g (%d)\n",
2543               "OUT-OF-PLANE BEND          : ", outOfPlaneBendEnergy, nOutOfPlaneBends));
2544     }
2545     if (torsionTerm) {
2546       sb.append(
2547           format(
2548               "REMARK   3   %s %g (%d)\n",
2549               "TORSIONAL ANGLE            : ", torsionEnergy, nTorsions));
2550     }
2551     if (piOrbitalTorsionTerm) {
2552       sb.append(
2553           format(
2554               "REMARK   3   %s %g (%d)\n",
2555               "PI-ORBITAL TORSION         : ", piOrbitalTorsionEnergy, nPiOrbitalTorsions));
2556     }
2557     if (torsionTorsionTerm) {
2558       sb.append(
2559           format(
2560               "REMARK   3   %s %g (%d)\n",
2561               "TORSION-TORSION            : ", torsionTorsionEnergy, nTorsionTorsions));
2562     }
2563     if (improperTorsionTerm) {
2564       sb.append(
2565           format(
2566               "REMARK   3   %s %g (%d)\n",
2567               "IMPROPER TORSION           : ", improperTorsionEnergy, nImproperTorsions));
2568     }
2569     if (restraintBondTerm) {
2570       sb.append(
2571           format(
2572               "REMARK   3   %s %g (%d)\n",
2573               "RESTRAINT BOND STRETCHING            : ", restraintBondEnergy, nRestraintBonds));
2574     }
2575 
2576     if (ncsTerm) {
2577       sb.append(
2578           format("REMARK   3   %s %g (%d)\n", "NCS RESTRAINT              : ", ncsEnergy, nAtoms));
2579     }
2580 
2581     if (restrainTerm && !coordRestraints.isEmpty()) {
2582       int nRests = 0;
2583       for (CoordRestraint restraint : coordRestraints) {
2584         nRests += restraint.getNumAtoms();
2585       }
2586       sb.append(
2587           format(
2588               "REMARK   3   %s %g (%d)\n",
2589               "COORDINATE RESTRAINTS      : ", restrainEnergy, nRests));
2590     }
2591 
2592     if (comTerm) {
2593       sb.append(
2594           format(
2595               "REMARK   3   %s %g (%d)\n",
2596               "COM RESTRAINT              : ", comRestraintEnergy, nAtoms));
2597     }
2598 
2599     if (vanderWaalsTerm) {
2600       sb.append(
2601           format(
2602               "REMARK   3   %s %g (%d)\n",
2603               "VAN DER WAALS              : ", vanDerWaalsEnergy, nVanDerWaalInteractions));
2604     }
2605     if (multipoleTerm) {
2606       sb.append(
2607           format(
2608               "REMARK   3   %s %g (%d)\n",
2609               "ATOMIC MULTIPOLES          : ", permanentMultipoleEnergy, nPermanentInteractions));
2610     }
2611     if (polarizationTerm) {
2612       sb.append(
2613           format(
2614               "REMARK   3   %s %g (%d)\n",
2615               "POLARIZATION               : ", polarizationEnergy, nPermanentInteractions));
2616     }
2617     sb.append(format("REMARK   3   %s %g\n", "TOTAL POTENTIAL (KCAL/MOL) : ", totalEnergy));
2618     int nsymm = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
2619     if (nsymm > 1) {
2620       sb.append(
2621           format("REMARK   3   %s %g\n", "UNIT CELL POTENTIAL        : ", totalEnergy * nsymm));
2622     }
2623     if (crystal.getUnitCell() != crystal) {
2624       nsymm = crystal.spaceGroup.getNumberOfSymOps();
2625       sb.append(
2626           format("REMARK   3   %s %g\n", "REPLICATES CELL POTENTIAL  : ", totalEnergy * nsymm));
2627     }
2628     sb.append("REMARK   3\n");
2629 
2630     return sb.toString();
2631   }
2632 
2633   /**
2634    * Getter for the field <code>parallelTeam</code>.
2635    *
2636    * @return a {@link edu.rit.pj.ParallelTeam} object.
2637    */
2638   public ParallelTeam getParallelTeam() {
2639     return parallelTeam;
2640   }
2641 
2642   /**
2643    * getPermanentInteractions.
2644    *
2645    * @return a int.
2646    */
2647   public int getPermanentInteractions() {
2648     return nPermanentInteractions;
2649   }
2650 
2651   /**
2652    * Getter for the field <code>permanentMultipoleEnergy</code>.
2653    *
2654    * @return a double.
2655    */
2656   public double getPermanentMultipoleEnergy() {
2657     return permanentMultipoleEnergy;
2658   }
2659 
2660   /**
2661    * Getter for the field <code>permanentRealSpaceEnergy</code>.
2662    *
2663    * @return a double.
2664    */
2665   public double getPermanentRealSpaceEnergy() {
2666     return permanentRealSpaceEnergy;
2667   }
2668 
2669   /**
2670    * getPermanentReciprocalMpoleEnergy.
2671    *
2672    * @return a double.
2673    */
2674   public double getPermanentReciprocalMpoleEnergy() {
2675     return particleMeshEwald.getPermRecipEnergy();
2676   }
2677 
2678   /**
2679    * getPermanentReciprocalSelfEnergy.
2680    *
2681    * @return a double.
2682    */
2683   public double getPermanentReciprocalSelfEnergy() {
2684     return particleMeshEwald.getPermSelfEnergy();
2685   }
2686 
2687   /**
2688    * Getter for the field <code>piOrbitalTorsionEnergy</code>.
2689    *
2690    * @return a double.
2691    */
2692   public double getPiOrbitalTorsionEnergy() {
2693     return piOrbitalTorsionEnergy;
2694   }
2695 
2696   /**
2697    * Getter for the field <code>piOrbitalTorsions</code>.
2698    *
2699    * @return an array of {@link ffx.potential.bonded.PiOrbitalTorsion} objects.
2700    */
2701   public PiOrbitalTorsion[] getPiOrbitalTorsions() {
2702     return piOrbitalTorsions;
2703   }
2704 
2705   /**
2706    * Gets the Platform associated with this force field energy. For the reference platform, always
2707    * returns FFX.
2708    *
2709    * @return A Platform.
2710    */
2711   public Platform getPlatform() {
2712     return platform;
2713   }
2714 
2715   /**
2716    * getPmeCartNode.
2717    *
2718    * @return a {@link ParticleMeshEwald} object.
2719    */
2720   public ParticleMeshEwald getPmeCartNode() {
2721     return (particleMeshEwald instanceof ParticleMeshEwald)
2722         ? (ParticleMeshEwald) particleMeshEwald
2723         : null;
2724   }
2725 
2726   /**
2727    * getPmeNode.
2728    *
2729    * @return a {@link ParticleMeshEwald} object.
2730    */
2731   public ParticleMeshEwald getPmeNode() {
2732     return particleMeshEwald;
2733   }
2734 
2735   /**
2736    * Getter for the field <code>polarizationEnergy</code>.
2737    *
2738    * @return a double.
2739    */
2740   public double getPolarizationEnergy() {
2741     return polarizationEnergy;
2742   }
2743 
2744 
2745   /**
2746    * {@inheritDoc}
2747    *
2748    * <p>Returns an array of previous acceleration values for active atoms.
2749    */
2750   @Override
2751   public double[] getPreviousAcceleration(double[] previousAcceleration) {
2752     int n = getNumberOfVariables();
2753     if (previousAcceleration == null || previousAcceleration.length < n) {
2754       previousAcceleration = new double[n];
2755     }
2756     int index = 0;
2757     double[] a = new double[3];
2758     for (int i = 0; i < nAtoms; i++) {
2759       if (atoms[i].isActive()) {
2760         atoms[i].getPreviousAcceleration(a);
2761         previousAcceleration[index++] = a[0];
2762         previousAcceleration[index++] = a[1];
2763         previousAcceleration[index++] = a[2];
2764       }
2765     }
2766     return previousAcceleration;
2767   }
2768 
2769   /**
2770    * Getter for the field <code>relativeSolvationEnergy</code>.
2771    *
2772    * @return a double.
2773    */
2774   public double getRelativeSolvationEnergy() {
2775     return relativeSolvationEnergy;
2776   }
2777 
2778   /** {@inheritDoc} */
2779   @Override
2780   public double[] getScaling() {
2781     return optimizationScaling;
2782   }
2783 
2784   /** {@inheritDoc} */
2785   @Override
2786   public void setScaling(double[] scaling) {
2787     optimizationScaling = scaling;
2788   }
2789 
2790   /**
2791    * Getter for the field <code>solvationEnergy</code>.
2792    *
2793    * @return a double.
2794    */
2795   public double getSolvationEnergy() {
2796     return solvationEnergy;
2797   }
2798 
2799   /**
2800    * getSolvationInteractions.
2801    *
2802    * @return a int.
2803    */
2804   public int getSolvationInteractions() {
2805     return nGKInteractions;
2806   }
2807 
2808   /**
2809    * getStrenchBendEnergy.
2810    *
2811    * @return a double.
2812    */
2813   public double getStrenchBendEnergy() {
2814     return stretchBendEnergy;
2815   }
2816 
2817   /**
2818    * Getter for the field <code>stretchBends</code>.
2819    *
2820    * @return an array of {@link ffx.potential.bonded.StretchBend} objects.
2821    */
2822   public StretchBend[] getStretchBends() {
2823     return stretchBends;
2824   }
2825 
2826   /**
2827    * Getter for the field <code>stretchTorsionEnergy</code>.
2828    *
2829    * @return a double.
2830    */
2831   public double getStretchTorsionEnergy() {
2832     return stretchTorsionEnergy;
2833   }
2834 
2835   /**
2836    * Getter for the field <code>stretchTorsions</code>.
2837    *
2838    * @return an array of {@link ffx.potential.bonded.StretchTorsion} objects.
2839    */
2840   public StretchTorsion[] getStretchTorsions() {
2841     return stretchTorsions;
2842   }
2843 
2844   /**
2845    * Getter for the field <code>torsionEnergy</code>.
2846    *
2847    * @return a double.
2848    */
2849   public double getTorsionEnergy() {
2850     return torsionEnergy;
2851   }
2852 
2853   /**
2854    * Getter for the field <code>torsionTorsionEnergy</code>.
2855    *
2856    * @return a double.
2857    */
2858   public double getTorsionTorsionEnergy() {
2859     return torsionTorsionEnergy;
2860   }
2861 
2862   /**
2863    * Getter for the field <code>torsionTorsions</code>.
2864    *
2865    * @return an array of {@link ffx.potential.bonded.TorsionTorsion} objects.
2866    */
2867   public TorsionTorsion[] getTorsionTorsions() {
2868     return torsionTorsions;
2869   }
2870 
2871   /**
2872    * Getter for the field <code>torsions</code>.
2873    *
2874    * @return an array of {@link ffx.potential.bonded.Torsion} objects.
2875    */
2876   public Torsion[] getTorsions() {
2877     return torsions;
2878   }
2879 
2880   /**
2881    * getTotalElectrostaticEnergy.
2882    *
2883    * @return a double.
2884    */
2885   public double getTotalElectrostaticEnergy() {
2886     return totalMultipoleEnergy + solvationEnergy;
2887   }
2888 
2889   /** {@inheritDoc} */
2890   @Override
2891   public double getTotalEnergy() {
2892     return totalEnergy;
2893   }
2894 
2895   /**
2896    * Getter for the field <code>ureyBradleyEnergy</code>.
2897    *
2898    * @return a double.
2899    */
2900   public double getUreyBradleyEnergy() {
2901     return ureyBradleyEnergy;
2902   }
2903 
2904   /**
2905    * Getter for the field <code>ureyBradleys</code>.
2906    *
2907    * @return an array of {@link ffx.potential.bonded.UreyBradley} objects.
2908    */
2909   public UreyBradley[] getUreyBradleys() {
2910     return ureyBradleys;
2911   }
2912 
2913   /**
2914    * Getter for the field <code>vanDerWaalsEnergy</code>.
2915    *
2916    * @return a double.
2917    */
2918   public double getVanDerWaalsEnergy() {
2919     return vanDerWaalsEnergy;
2920   }
2921 
2922   /**
2923    * getVanDerWaalsInteractions.
2924    *
2925    * @return a int.
2926    */
2927   public int getVanDerWaalsInteractions() {
2928     return nVanDerWaalInteractions;
2929   }
2930 
2931   /**
2932    * {@inheritDoc}
2933    *
2934    * <p>Return a reference to each variables type.
2935    */
2936   @Override
2937   public VARIABLE_TYPE[] getVariableTypes() {
2938     int n = getNumberOfVariables();
2939     VARIABLE_TYPE[] type = new VARIABLE_TYPE[n];
2940     int i = 0;
2941     for (int j = 0; j < nAtoms; j++) {
2942       if (atoms[j].isActive()) {
2943         type[i++] = VARIABLE_TYPE.X;
2944         type[i++] = VARIABLE_TYPE.Y;
2945         type[i++] = VARIABLE_TYPE.Z;
2946       }
2947     }
2948     return type;
2949   }
2950 
2951   /**
2952    * getVdwNode.
2953    *
2954    * @return a {@link ffx.potential.nonbonded.VanDerWaals} object.
2955    */
2956   public VanDerWaals getVdwNode() {
2957     return vanderWaals;
2958   }
2959 
2960   /**
2961    * {@inheritDoc}
2962    *
2963    * <p>Returns an array of velocity values for active atoms.
2964    */
2965   @Override
2966   public double[] getVelocity(double[] velocity) {
2967     int n = getNumberOfVariables();
2968     if (velocity == null || velocity.length < n) {
2969       velocity = new double[n];
2970     }
2971     int index = 0;
2972     double[] v = new double[3];
2973     for (int i = 0; i < nAtoms; i++) {
2974       Atom a = atoms[i];
2975       if (a.isActive()) {
2976         a.getVelocity(v);
2977         velocity[index++] = v[0];
2978         velocity[index++] = v[1];
2979         velocity[index++] = v[2];
2980       }
2981     }
2982     return velocity;
2983   }
2984 
2985   /** {@inheritDoc} */
2986   @Override
2987   public double getd2EdL2() {
2988     double d2EdLambda2 = 0.0;
2989     if (!lambdaBondedTerms) {
2990       if (vanderWaalsTerm) {
2991         d2EdLambda2 = vanderWaals.getd2EdL2();
2992       }
2993       if (multipoleTerm) {
2994         d2EdLambda2 += particleMeshEwald.getd2EdL2();
2995       }
2996       if (restraintBondTerm) {
2997         for (int i = 0; i < nRestraintBonds; i++) {
2998           d2EdLambda2 += restraintBonds[i].getd2EdL2();
2999         }
3000       }
3001       if (ncsTerm && ncsRestraint != null) {
3002         d2EdLambda2 += ncsRestraint.getd2EdL2();
3003       }
3004       if (restrainTerm && !coordRestraints.isEmpty()) {
3005         for (CoordRestraint restraint : coordRestraints) {
3006           d2EdLambda2 += restraint.getd2EdL2();
3007         }
3008       }
3009       if (comTerm && comRestraint != null) {
3010         d2EdLambda2 += comRestraint.getd2EdL2();
3011       }
3012       if (lambdaTorsions) {
3013         for (int i = 0; i < nTorsions; i++) {
3014           d2EdLambda2 += torsions[i].getd2EdL2();
3015         }
3016         for (int i = 0; i < nPiOrbitalTorsions; i++) {
3017           d2EdLambda2 += piOrbitalTorsions[i].getd2EdL2();
3018         }
3019         for (int i = 0; i < nTorsionTorsions; i++) {
3020           d2EdLambda2 += torsionTorsions[i].getd2EdL2();
3021         }
3022       }
3023     }
3024     return d2EdLambda2;
3025   }
3026 
3027   /** {@inheritDoc} */
3028   @Override
3029   public double getdEdL() {
3030     double dEdLambda = 0.0;
3031     if (!lambdaBondedTerms) {
3032       if (vanderWaalsTerm) {
3033         dEdLambda = vanderWaals.getdEdL();
3034       }
3035       if (multipoleTerm) {
3036         dEdLambda += particleMeshEwald.getdEdL();
3037       }
3038       if (restraintBondTerm) {
3039         for (int i = 0; i < nRestraintBonds; i++) {
3040           dEdLambda += restraintBonds[i].getdEdL();
3041         }
3042       }
3043       if (ncsTerm && ncsRestraint != null) {
3044         dEdLambda += ncsRestraint.getdEdL();
3045       }
3046       if (restrainTerm && !coordRestraints.isEmpty()) {
3047         for (CoordRestraint restraint : coordRestraints) {
3048           dEdLambda += restraint.getdEdL();
3049         }
3050       }
3051       if (comTerm && comRestraint != null) {
3052         dEdLambda += comRestraint.getdEdL();
3053       }
3054       if (lambdaTorsions) {
3055         for (int i = 0; i < nTorsions; i++) {
3056           dEdLambda += torsions[i].getdEdL();
3057         }
3058         for (int i = 0; i < nPiOrbitalTorsions; i++) {
3059           dEdLambda += piOrbitalTorsions[i].getdEdL();
3060         }
3061         for (int i = 0; i < nTorsionTorsions; i++) {
3062           dEdLambda += torsionTorsions[i].getdEdL();
3063         }
3064       }
3065     }
3066     return dEdLambda;
3067   }
3068 
3069   /** {@inheritDoc} */
3070   @Override
3071   public void getdEdXdL(double[] gradients) {
3072     if (!lambdaBondedTerms) {
3073       if (vanderWaalsTerm) {
3074         vanderWaals.getdEdXdL(gradients);
3075       }
3076       if (multipoleTerm) {
3077         particleMeshEwald.getdEdXdL(gradients);
3078       }
3079       if (restraintBondTerm) {
3080         for (int i = 0; i < nRestraintBonds; i++) {
3081           restraintBonds[i].getdEdXdL(gradients);
3082         }
3083       }
3084       if (ncsTerm && ncsRestraint != null) {
3085         ncsRestraint.getdEdXdL(gradients);
3086       }
3087       if (restrainTerm && !coordRestraints.isEmpty()) {
3088         for (CoordRestraint restraint : coordRestraints) {
3089           restraint.getdEdXdL(gradients);
3090         }
3091       }
3092       if (comTerm && comRestraint != null) {
3093         comRestraint.getdEdXdL(gradients);
3094       }
3095       if (lambdaTorsions) {
3096         double[] grad = new double[3];
3097         int index = 0;
3098         for (int i = 0; i < nAtoms; i++) {
3099           Atom a = atoms[i];
3100           if (a.isActive()) {
3101             a.getLambdaXYZGradient(grad);
3102             gradients[index++] += grad[0];
3103             gradients[index++] += grad[1];
3104             gradients[index++] += grad[2];
3105           }
3106         }
3107       }
3108     }
3109   }
3110 
3111   /** Need to remove degrees of freedom that are lost to prevent heating. */
3112   @Deprecated
3113   public void reInit() {
3114     int[] molecule;
3115     atoms = molecularAssembly.getAtomArray();
3116     molecule = molecularAssembly.getMoleculeNumbers();
3117     nAtoms = atoms.length;
3118 
3119     xyz = new double[nAtoms * 3];
3120     getCoordinates(xyz);
3121 
3122     // Check that atom ordering is correct and count number of Active atoms.
3123     for (int i = 0; i < nAtoms; i++) {
3124       Atom atom = atoms[i];
3125       int index = atom.getXyzIndex() - 1;
3126       if (index != i) {
3127         atom.setXyzIndex(i + 1);
3128       }
3129     }
3130 
3131     // Collect, count, pack and sort bonds.
3132     if (bondTerm) {
3133       List<Bond> bond = molecularAssembly.getBondList();
3134       nBonds = 0;
3135       for (Bond r : bond) {
3136         if (keep(r)) {
3137           nBonds++;
3138         }
3139       }
3140       if (nBonds > bonds.length) {
3141         bonds = new Bond[nBonds];
3142       }
3143       fill(bonds, null);
3144       nBonds = 0;
3145       for (Bond r : bond) {
3146         if (keep(r)) {
3147           bonds[nBonds++] = r;
3148         }
3149       }
3150       sort(bonds, 0, nBonds);
3151       if (nBonds > 0 && logger.isLoggable(Level.FINEST)) {
3152         logger.finest(format("  Bonds:                             %10d", nBonds));
3153       }
3154     } else {
3155       nBonds = 0;
3156       bonds = null;
3157     }
3158 
3159     // Collect, count, pack and sort angles.
3160     if (angleTerm) {
3161       List<Angle> angle = molecularAssembly.getAngleList();
3162       nAngles = 0;
3163       for (Angle r : angle) {
3164         if (keep(r)) {
3165           nAngles++;
3166         }
3167       }
3168       if (nAngles > angles.length) {
3169         angles = new Angle[nAngles];
3170       }
3171       fill(angles, null);
3172       nAngles = 0;
3173       for (Angle r : angle) {
3174         if (keep(r)) {
3175           angles[nAngles++] = r;
3176         }
3177       }
3178 
3179       sort(angles, 0, nAngles);
3180       if (nAngles > 0 && logger.isLoggable(Level.FINEST)) {
3181         logger.finest(format("  Angles:                            %10d", nAngles));
3182       }
3183     } else {
3184       nAngles = 0;
3185       angles = null;
3186     }
3187 
3188     // Collect, count, pack and sort stretch-bends.
3189     if (stretchBendTerm) {
3190       List<StretchBend> stretchBend = molecularAssembly.getStretchBendList();
3191       nStretchBends = 0;
3192       for (StretchBend r : stretchBend) {
3193         if (keep(r)) {
3194           nStretchBends++;
3195         }
3196       }
3197       if (nStretchBends > stretchBends.length) {
3198         stretchBends = new StretchBend[nStretchBends];
3199       }
3200       fill(stretchBends, null);
3201       nStretchBends = 0;
3202       for (StretchBend r : stretchBend) {
3203         if (keep(r)) {
3204           stretchBends[nStretchBends++] = r;
3205         }
3206       }
3207       sort(stretchBends, 0, nStretchBends);
3208       if (nStretchBends > 0 && logger.isLoggable(Level.FINEST)) {
3209         logger.finest(format("  Stretch-Bends:                     %10d", nStretchBends));
3210       }
3211     } else {
3212       nStretchBends = 0;
3213       stretchBends = null;
3214     }
3215 
3216     // Collect, count, pack and sort Urey-Bradleys.
3217     if (ureyBradleyTerm) {
3218       List<UreyBradley> ureyBradley = molecularAssembly.getUreyBradleyList();
3219       nUreyBradleys = 0;
3220       for (UreyBradley r : ureyBradley) {
3221         if (keep(r)) {
3222           nUreyBradleys++;
3223         }
3224       }
3225       if (nUreyBradleys > ureyBradleys.length) {
3226         ureyBradleys = new UreyBradley[nUreyBradleys];
3227       }
3228       fill(ureyBradleys, null);
3229       nUreyBradleys = 0;
3230       for (UreyBradley r : ureyBradley) {
3231         if (keep(r)) {
3232           ureyBradleys[nUreyBradleys++] = r;
3233         }
3234       }
3235       sort(ureyBradleys, 0, nUreyBradleys);
3236       if (nUreyBradleys > 0 && logger.isLoggable(Level.FINEST)) {
3237         logger.finest(format("  Urey-Bradleys:                     %10d", nUreyBradleys));
3238       }
3239     } else {
3240       nUreyBradleys = 0;
3241       ureyBradleys = null;
3242     }
3243 
3244     // Set a multiplier on the force constants of bonded terms containing hydrogens.
3245     if (rigidHydrogens) {
3246       if (bonds != null) {
3247         for (Bond bond : bonds) {
3248           if (bond.containsHydrogen()) {
3249             bond.setRigidScale(rigidScale);
3250           }
3251         }
3252       }
3253       if (angles != null) {
3254         for (Angle angle : angles) {
3255           if (angle.containsHydrogen()) {
3256             angle.setRigidScale(rigidScale);
3257           }
3258         }
3259       }
3260       if (stretchBends != null) {
3261         for (StretchBend stretchBend : stretchBends) {
3262           if (stretchBend.containsHydrogen()) {
3263             stretchBend.setRigidScale(rigidScale);
3264           }
3265         }
3266       }
3267       if (ureyBradleys != null) {
3268         for (UreyBradley ureyBradley : ureyBradleys) {
3269           if (ureyBradley.containsHydrogen()) {
3270             ureyBradley.setRigidScale(rigidScale);
3271           }
3272         }
3273       }
3274     }
3275 
3276     // Collect, count, pack and sort out-of-plane bends.
3277     if (outOfPlaneBendTerm) {
3278       List<OutOfPlaneBend> outOfPlaneBend = molecularAssembly.getOutOfPlaneBendList();
3279       nOutOfPlaneBends = 0;
3280       for (OutOfPlaneBend r : outOfPlaneBend) {
3281         if (keep(r)) {
3282           nOutOfPlaneBends++;
3283         }
3284       }
3285       if (nOutOfPlaneBends > outOfPlaneBends.length) {
3286         outOfPlaneBends = new OutOfPlaneBend[nOutOfPlaneBends];
3287       }
3288       fill(outOfPlaneBends, null);
3289       nOutOfPlaneBends = 0;
3290       for (OutOfPlaneBend r : outOfPlaneBend) {
3291         if (keep(r)) {
3292           outOfPlaneBends[nOutOfPlaneBends++] = r;
3293         }
3294       }
3295       sort(outOfPlaneBends, 0, nOutOfPlaneBends);
3296       if (nOutOfPlaneBends > 0 && logger.isLoggable(Level.FINEST)) {
3297         logger.finest(format("  Out-of-Plane Bends:                %10d", nOutOfPlaneBends));
3298       }
3299     } else {
3300       nOutOfPlaneBends = 0;
3301       outOfPlaneBends = null;
3302     }
3303 
3304     // Collect, count, pack and sort torsions.
3305     if (torsionTerm) {
3306       List<Torsion> torsion = molecularAssembly.getTorsionList();
3307       nTorsions = 0;
3308       for (Torsion r : torsion) {
3309         if (keep(r)) {
3310           nTorsions++;
3311         }
3312       }
3313       if (nTorsions >= torsions.length) {
3314         torsions = new Torsion[nTorsions];
3315       }
3316       fill(torsions, null);
3317       nTorsions = 0;
3318       for (Torsion r : torsion) {
3319         if (keep(r)) {
3320           torsions[nTorsions++] = r;
3321         }
3322       }
3323       // Arrays.sort(torsions);
3324       if (nTorsions > 0 && logger.isLoggable(Level.FINEST)) {
3325         logger.finest(format("  Torsions:                          %10d", nTorsions));
3326       }
3327     } else {
3328       nTorsions = 0;
3329       torsions = null;
3330     }
3331 
3332     // Collect, count, pack and sort pi-orbital torsions.
3333     if (piOrbitalTorsionTerm) {
3334       List<PiOrbitalTorsion> piOrbitalTorsion = molecularAssembly.getPiOrbitalTorsionList();
3335       nPiOrbitalTorsions = 0;
3336       for (PiOrbitalTorsion r : piOrbitalTorsion) {
3337         if (keep(r)) {
3338           nPiOrbitalTorsions++;
3339         }
3340       }
3341       if (nPiOrbitalTorsions >= piOrbitalTorsions.length) {
3342         piOrbitalTorsions = new PiOrbitalTorsion[nPiOrbitalTorsions];
3343       }
3344       fill(piOrbitalTorsions, null);
3345       nPiOrbitalTorsions = 0;
3346       for (PiOrbitalTorsion r : piOrbitalTorsion) {
3347         if (keep(r)) {
3348           piOrbitalTorsions[nPiOrbitalTorsions++] = r;
3349         }
3350       }
3351       if (nPiOrbitalTorsions > 0 && logger.isLoggable(Level.FINEST)) {
3352         logger.finest(format("  Pi-Orbital Torsions:               %10d", nPiOrbitalTorsions));
3353       }
3354     } else {
3355       nPiOrbitalTorsions = 0;
3356       piOrbitalTorsions = null;
3357     }
3358 
3359     // Collect, count, pack and sort torsion-torsions.
3360     if (torsionTorsionTerm) {
3361       List<TorsionTorsion> torsionTorsion = molecularAssembly.getTorsionTorsionList();
3362       nTorsionTorsions = 0;
3363       for (TorsionTorsion r : torsionTorsion) {
3364         if (keep(r)) {
3365           nTorsionTorsions++;
3366         }
3367       }
3368       if (nTorsionTorsions >= torsionTorsions.length) {
3369         torsionTorsions = new TorsionTorsion[nTorsionTorsions];
3370       }
3371       fill(torsionTorsions, null);
3372       nTorsionTorsions = 0;
3373       for (TorsionTorsion r : torsionTorsion) {
3374         if (keep(r)) {
3375           torsionTorsions[nTorsionTorsions++] = r;
3376         }
3377       }
3378       if (nTorsionTorsions > 0 && logger.isLoggable(Level.FINEST)) {
3379         logger.finest(format("  Torsion-Torsions:                  %10d", nTorsionTorsions));
3380       }
3381     } else {
3382       nTorsionTorsions = 0;
3383       torsionTorsions = null;
3384     }
3385 
3386     // Collect, count, pack and sort improper torsions.
3387     if (improperTorsionTerm) {
3388       List<ImproperTorsion> improperTorsion = molecularAssembly.getImproperTorsionList();
3389       nImproperTorsions = 0;
3390       for (ImproperTorsion r : improperTorsion) {
3391         if (keep(r)) {
3392           nImproperTorsions++;
3393         }
3394       }
3395       if (nImproperTorsions >= improperTorsions.length) {
3396         improperTorsions = new ImproperTorsion[nImproperTorsions];
3397       }
3398       fill(improperTorsions, null);
3399       nImproperTorsions = 0;
3400       for (ImproperTorsion r : improperTorsion) {
3401         if (keep(r)) {
3402           improperTorsions[nImproperTorsions++] = r;
3403         }
3404       }
3405       if (nImproperTorsions > 0 && logger.isLoggable(Level.FINEST)) {
3406         logger.finest(format("  Improper Torsions:                 %10d", nImproperTorsions));
3407       }
3408     } else {
3409       nImproperTorsions = 0;
3410       improperTorsions = null;
3411     }
3412 
3413     if (vanderWaalsTerm) {
3414       if (esvTerm) {
3415         vanderWaals.setAtoms(esvSystem.getExtendedAtoms(), esvSystem.getExtendedMolecule());
3416       } else {
3417         vanderWaals.setAtoms(atoms, molecule);
3418       }
3419     }
3420 
3421     if (multipoleTerm) {
3422       if (esvTerm) {
3423         particleMeshEwald.setAtoms(esvSystem.getExtendedAtoms(), esvSystem.getExtendedMolecule());
3424       } else {
3425         particleMeshEwald.setAtoms(atoms, molecule);
3426       }
3427     }
3428 
3429     if (ncsTerm) {
3430       logger.severe(" NCS energy term cannot be used with variable systems sizes.");
3431     }
3432 
3433     if (restrainTerm) {
3434       logger.severe(" Restrain energy term cannot be used with variable systems sizes.");
3435     }
3436 
3437     if (comTerm) {
3438       logger.severe(" COM restrain energy term cannot be used with variable systems sizes.");
3439     }
3440 
3441     bondedRegion = new BondedRegion();
3442   }
3443 
3444   /**
3445    * {@inheritDoc}
3446    *
3447    * <p>The acceleration array should only contain acceleration data for active atoms.
3448    */
3449   @Override
3450   public void setAcceleration(double[] acceleration) {
3451     if (acceleration == null) {
3452       return;
3453     }
3454     int index = 0;
3455     double[] accel = new double[3];
3456     for (int i = 0; i < nAtoms; i++) {
3457       if (atoms[i].isActive()) {
3458         accel[0] = acceleration[index++];
3459         accel[1] = acceleration[index++];
3460         accel[2] = acceleration[index++];
3461         atoms[i].setAcceleration(accel);
3462       }
3463     }
3464   }
3465 
3466   /**
3467    * The coordinate array should only contain active atoms.
3468    *
3469    * @param coords an array of {@link double} objects.
3470    */
3471   public void setCoordinates(double[] coords) {
3472     if (coords == null) {
3473       return;
3474     }
3475     int index = 0;
3476     for (int i = 0; i < nAtoms; i++) {
3477       Atom a = atoms[i];
3478       if (a.isActive() && !a.isBackground()) {
3479         double x = coords[index++];
3480         double y = coords[index++];
3481         double z = coords[index++];
3482         a.moveTo(x, y, z);
3483       }
3484     }
3485   }
3486 
3487   /**
3488    * Set the boundary conditions for this calculation.
3489    *
3490    * @param crystal Crystal to set.
3491    * @param checkReplicatesCell Check if a replicates cell must be created.
3492    */
3493   public void setCrystal(Crystal crystal, boolean checkReplicatesCell) {
3494     if (checkReplicatesCell) {
3495       this.crystal = ReplicatesCrystal.replicatesCrystalFactory(crystal.getUnitCell(), cutOff2);
3496     } else {
3497       this.crystal = crystal;
3498     }
3499     /*
3500      Update VanDerWaals first, in case the NeighborList needs to be
3501      re-allocated to include a larger number of replicated cells.
3502     */
3503     if (vanderWaalsTerm) {
3504       vanderWaals.setCrystal(this.crystal);
3505     }
3506     if (multipoleTerm) {
3507       particleMeshEwald.setCrystal(this.crystal);
3508     }
3509   }
3510 
3511   /**
3512    * setFixedCharges.
3513    *
3514    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
3515    */
3516   public void setFixedCharges(Atom[] atoms) {
3517     if (particleMeshEwald != null) {
3518       particleMeshEwald.setFixedCharges(atoms);
3519     }
3520   }
3521 
3522   /**
3523    * setLambdaMultipoleScale.
3524    *
3525    * @param scale a double.
3526    */
3527   public void setLambdaMultipoleScale(double scale) {
3528     if (particleMeshEwald != null) {
3529       particleMeshEwald.setLambdaMultipoleScale(scale);
3530     }
3531   }
3532 
3533   /**
3534    * {@inheritDoc}
3535    *
3536    * <p>The previousAcceleration array should only contain previous acceleration data for active
3537    * atoms.
3538    */
3539   @Override
3540   public void setPreviousAcceleration(double[] previousAcceleration) {
3541     if (previousAcceleration == null) {
3542       return;
3543     }
3544     int index = 0;
3545     double[] prev = new double[3];
3546     for (int i = 0; i < nAtoms; i++) {
3547       if (atoms[i].isActive()) {
3548         prev[0] = previousAcceleration[index++];
3549         prev[1] = previousAcceleration[index++];
3550         prev[2] = previousAcceleration[index++];
3551         atoms[i].setPreviousAcceleration(prev);
3552       }
3553     }
3554   }
3555 
3556   /**
3557    * Sets the printOnFailure flag; if override is true, over-rides any existing property. Essentially
3558    * sets the default value of printOnFailure for an algorithm. For example, rotamer optimization
3559    * will generally run into force field issues in the normal course of execution as it tries
3560    * unphysical self and pair configurations, so the algorithm should not print out a large number of
3561    * error PDBs.
3562    *
3563    * @param onFail To set
3564    * @param override Override properties
3565    */
3566   public void setPrintOnFailure(boolean onFail, boolean override) {
3567     if (override) {
3568       // Ignore any pre-existing value
3569       printOnFailure = onFail;
3570     } else {
3571       try {
3572         molecularAssembly.getForceField().getBoolean("PRINT_ON_FAILURE");
3573         /*
3574          If the call was successful, the property was explicitly set
3575          somewhere and should be kept. If an exception was thrown, the
3576          property was never set explicitly, so over-write.
3577         */
3578       } catch (Exception ex) {
3579         printOnFailure = onFail;
3580       }
3581     }
3582   }
3583 
3584   /**
3585    * Setter for the field <code>resolution</code>.
3586    *
3587    * @param resolution a {@link ffx.potential.bonded.Atom.Resolution} object.
3588    */
3589   public void setResolution(Resolution resolution) {
3590     this.resolution = resolution;
3591 
3592     if (vanderWaals != null) {
3593       vanderWaals.setResolution(resolution);
3594     }
3595 
3596     if (resolution == Resolution.FIXEDCHARGE) {
3597       multipoleTerm = false;
3598       polarizationTerm = false;
3599       generalizedKirkwoodTerm = false;
3600     }
3601   }
3602 
3603   /**
3604    * setRestraintBond
3605    *
3606    * @param a1 a {@link ffx.potential.bonded.Atom} object.
3607    * @param a2 a {@link ffx.potential.bonded.Atom} object.
3608    * @param distance a double.
3609    * @param forceConstant the force constant in kcal/mole
3610    */
3611   public void setRestraintBond(Atom a1, Atom a2, double distance, double forceConstant) {
3612     setRestraintBond(a1, a2, distance, forceConstant, 0);
3613   }
3614 
3615   /**
3616    * setRestraintBond
3617    *
3618    * @param a1 a {@link ffx.potential.bonded.Atom} object.
3619    * @param a2 a {@link ffx.potential.bonded.Atom} object.
3620    * @param distance a double.
3621    * @param forceConstant the force constant in kcal/mole.
3622    * @param flatBottom Radius of a flat-bottom potential in Angstroms.
3623    */
3624   public void setRestraintBond(
3625       Atom a1, Atom a2, double distance, double forceConstant, double flatBottom) {
3626     setRestraintBond(
3627         a1,
3628         a2,
3629         distance,
3630         forceConstant,
3631         flatBottom,
3632         RestraintBond.DEFAULT_RB_LAM_START,
3633         RestraintBond.DEFAULT_RB_LAM_END,
3634         new ConstantSwitch());
3635   }
3636 
3637   /**
3638    * {@inheritDoc}
3639    *
3640    * <p>The velocity array should only contain velocity data for active atoms.
3641    */
3642   @Override
3643   public void setVelocity(double[] velocity) {
3644     if (velocity == null) {
3645       return;
3646     }
3647     int index = 0;
3648     double[] vel = new double[3];
3649     for (int i = 0; i < nAtoms; i++) {
3650       if (atoms[i].isActive()) {
3651         vel[0] = velocity[index++];
3652         vel[1] = velocity[index++];
3653         vel[2] = velocity[index++];
3654         atoms[i].setVelocity(vel);
3655       }
3656     }
3657   }
3658 
3659   /** {@inheritDoc} */
3660   @Override
3661   public String toString() {
3662     StringBuilder sb = new StringBuilder();
3663 
3664     if (bondTerm && nBonds > 0) {
3665       sb.append(
3666           format(
3667               "  %s %16.8f %12d %12.3f (%8.5f)\n",
3668               "Bond Stretching   ", bondEnergy, nBonds, bondTime * toSeconds, bondRMSD));
3669     }
3670     if (angleTerm && nAngles > 0) {
3671       sb.append(
3672           format(
3673               "  %s %16.8f %12d %12.3f (%8.5f)\n",
3674               "Angle Bending     ", angleEnergy, nAngles, angleTime * toSeconds, angleRMSD));
3675     }
3676     if (stretchBendTerm && nStretchBends > 0) {
3677       sb.append(
3678           format(
3679               "  %s %16.8f %12d %12.3f\n",
3680               "Stretch-Bend      ", stretchBendEnergy, nStretchBends, stretchBendTime * toSeconds));
3681     }
3682     if (ureyBradleyTerm && nUreyBradleys > 0) {
3683       sb.append(
3684           format(
3685               "  %s %16.8f %12d %12.3f\n",
3686               "Urey-Bradley      ", ureyBradleyEnergy, nUreyBradleys, ureyBradleyTime * toSeconds));
3687     }
3688     if (outOfPlaneBendTerm && nOutOfPlaneBends > 0) {
3689       sb.append(
3690           format(
3691               "  %s %16.8f %12d %12.3f\n",
3692               "Out-of-Plane Bend ",
3693               outOfPlaneBendEnergy,
3694               nOutOfPlaneBends,
3695               outOfPlaneBendTime * toSeconds));
3696     }
3697     if (torsionTerm && nTorsions > 0) {
3698       sb.append(
3699           format(
3700               "  %s %16.8f %12d %12.3f\n",
3701               "Torsional Angle   ", torsionEnergy, nTorsions, torsionTime * toSeconds));
3702     }
3703     if (piOrbitalTorsionTerm && nPiOrbitalTorsions > 0) {
3704       sb.append(
3705           format(
3706               "  %s %16.8f %12d %12.3f\n",
3707               "Pi-Orbital Torsion",
3708               piOrbitalTorsionEnergy,
3709               nPiOrbitalTorsions,
3710               piOrbitalTorsionTime * toSeconds));
3711     }
3712     if (stretchTorsionTerm && nStretchTorsions > 0) {
3713       sb.append(
3714           format(
3715               "  %s %16.8f %12d %12.3f\n",
3716               "Stretch-Torsion   ",
3717               stretchTorsionEnergy,
3718               nStretchTorsions,
3719               stretchTorsionTime * toSeconds));
3720     }
3721     if (angleTorsionTerm && nAngleTorsions > 0) {
3722       sb.append(
3723           format(
3724               "  %s %16.8f %12d %12.3f\n",
3725               "Angle-Torsion     ",
3726               angleTorsionEnergy,
3727               nAngleTorsions,
3728               angleTorsionTime * toSeconds));
3729     }
3730     if (torsionTorsionTerm && nTorsionTorsions > 0) {
3731       sb.append(
3732           format(
3733               "  %s %16.8f %12d %12.3f\n",
3734               "Torsion-Torsion   ",
3735               torsionTorsionEnergy,
3736               nTorsionTorsions,
3737               torsionTorsionTime * toSeconds));
3738     }
3739     if (improperTorsionTerm && nImproperTorsions > 0) {
3740       sb.append(
3741           format(
3742               "  %s %16.8f %12d %12.3f\n",
3743               "Improper Torsion  ",
3744               improperTorsionEnergy,
3745               nImproperTorsions,
3746               improperTorsionTime * toSeconds));
3747     }
3748     if (restraintBondTerm && nRestraintBonds > 0) {
3749       sb.append(
3750           format(
3751               "  %s %16.8f %12d %12.3f\n",
3752               "Bond Restraint    ",
3753               restraintBondEnergy,
3754               nRestraintBonds,
3755               restraintBondTime * toSeconds));
3756     }
3757     if (ncsTerm) {
3758       sb.append(
3759           format(
3760               "  %s %16.8f %12d %12.3f\n",
3761               "NCS Restraint     ", ncsEnergy, nAtoms, ncsTime * toSeconds));
3762     }
3763     if (restrainTerm && !coordRestraints.isEmpty()) {
3764       int nRests = 0;
3765       for (CoordRestraint restraint : coordRestraints) {
3766         nRests += restraint.getNumAtoms();
3767       }
3768       sb.append(
3769           format(
3770               "  %s %16.8f %12d %12.3f\n",
3771               "Coord. Restraints ", restrainEnergy, nRests, coordRestraintTime * toSeconds));
3772     }
3773     if (comTerm) {
3774       sb.append(
3775           format(
3776               "  %s %16.8f %12d %12.3f\n",
3777               "COM Restraint     ", comRestraintEnergy, nAtoms, comRestraintTime * toSeconds));
3778     }
3779     if (restrainGroupTerm) {
3780       sb.append(
3781           format(
3782               "  %s %16.8f %12d %12.3f\n",
3783               "Restrain Groups   ",
3784               restrainGroupEnergy,
3785               nRestrainGroups,
3786               restrainGroupTime * toSeconds));
3787     }
3788     if (rTorsTerm) {
3789       sb.append(
3790           format("  %s %16.8f %12d %12.3f\n",
3791               "Dihedral Restraints",
3792               rTorsEnergy,
3793               nRestTors,
3794               rTorsTime * toSeconds)
3795       );
3796     }
3797     if (vanderWaalsTerm && nVanDerWaalInteractions > 0) {
3798       sb.append(
3799           format(
3800               "  %s %16.8f %12d %12.3f\n",
3801               "Van der Waals     ",
3802               vanDerWaalsEnergy,
3803               nVanDerWaalInteractions,
3804               vanDerWaalsTime * toSeconds));
3805     }
3806     if (multipoleTerm && nPermanentInteractions > 0) {
3807       String pmeTitle = "Atomic Multipoles ";
3808       if (polarizationTerm) {
3809         sb.append(format("  %s %16.8f %12d\n", pmeTitle, permanentMultipoleEnergy,
3810             nPermanentInteractions));
3811       } else {
3812         sb.append(
3813             format(
3814                 "  %s %16.8f %12d %12.3f\n",
3815                 pmeTitle,
3816                 permanentMultipoleEnergy,
3817                 nPermanentInteractions,
3818                 electrostaticTime * toSeconds));
3819       }
3820     }
3821     if (polarizationTerm && nPermanentInteractions > 0) {
3822       sb.append(
3823           format(
3824               "  %s %16.8f %12d %12.3f\n",
3825               "Polarization      ",
3826               polarizationEnergy,
3827               nPermanentInteractions,
3828               electrostaticTime * toSeconds));
3829     }
3830     if (generalizedKirkwoodTerm && nGKInteractions > 0) {
3831       sb.append(
3832           format("  %s %16.8f %12d\n", "Solvation         ", solvationEnergy, nGKInteractions));
3833     }
3834     if (relativeSolvationTerm) {
3835       sb.append(
3836           format(
3837               "  %s %16.8f %12d\n",
3838               "Relative Solvation", relativeSolvationEnergy, nRelativeSolvations));
3839     }
3840 
3841     if (esvTerm) {
3842       sb.append(
3843           format("  %s %16.8f  %s\n", "ExtendedSystemBias", esvBias, esvSystem.getLambdaList()));
3844       sb.append(esvSystem.getBiasDecomposition());
3845     }
3846 
3847     sb.append(
3848         format(
3849             "  %s %16.8f  %s %12.3f (sec)",
3850             "Total Potential   ", totalEnergy, "(Kcal/mole)", totalTime * toSeconds));
3851 
3852     int nsymm = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
3853     if (nsymm > 1) {
3854       sb.append(format("\n  %s %16.8f", "Unit Cell         ", totalEnergy * nsymm));
3855     }
3856     if (crystal.getUnitCell() != crystal) {
3857       nsymm = crystal.spaceGroup.getNumberOfSymOps();
3858       sb.append(format("\n  %s %16.8f", "Replicates Cell   ", totalEnergy * nsymm));
3859     }
3860     return sb.toString();
3861   }
3862 
3863   private boolean keep(BondedTerm term) {
3864     switch (resolution) {
3865       case AMOEBA:
3866         return term.isResolution(Resolution.AMOEBA);
3867       case FIXEDCHARGE:
3868         return term.containsResolution(Resolution.FIXEDCHARGE);
3869       default:
3870         return true;
3871     }
3872   }
3873 
3874   /**
3875    * Setter for the field <code>lambdaBondedTerms</code>.
3876    *
3877    * @param lambdaBondedTerms a boolean.
3878    * @param lambdaAllBondedTerms a boolean.
3879    */
3880   void setLambdaBondedTerms(boolean lambdaBondedTerms, boolean lambdaAllBondedTerms) {
3881     this.lambdaBondedTerms = lambdaBondedTerms;
3882     this.lambdaAllBondedTerms = lambdaAllBondedTerms;
3883   }
3884 
3885   /**
3886    * Return the non-bonded components of energy (vdW, electrostatics).
3887    *
3888    * @param includeSolv Include solvation energy
3889    * @return Nonbonded energy
3890    */
3891   private double getNonbondedEnergy(boolean includeSolv) {
3892     return (includeSolv ? (totalNonBondedEnergy + solvationEnergy) : totalNonBondedEnergy);
3893   }
3894 
3895   /**
3896    * Private method for internal use, so we don't have subclasses calling super.energy, and this
3897    * class delegating to the subclass's getGradient method.
3898    *
3899    * @param g Gradient array to fill.
3900    * @return Gradient array.
3901    */
3902   private double[] fillGradient(double[] g) {
3903     assert (g != null);
3904     double[] grad = new double[3];
3905     int n = getNumberOfVariables();
3906     if (g == null || g.length < n) {
3907       g = new double[n];
3908     }
3909     int index = 0;
3910     for (int i = 0; i < nAtoms; i++) {
3911       Atom a = atoms[i];
3912       if (a.isActive()) {
3913         a.getXYZGradient(grad);
3914         double gx = grad[0];
3915         double gy = grad[1];
3916         double gz = grad[2];
3917         if (isNaN(gx)
3918             || isInfinite(gx)
3919             || isNaN(gy)
3920             || isInfinite(gy)
3921             || isNaN(gz)
3922             || isInfinite(gz)) {
3923           StringBuilder sb =
3924               new StringBuilder(
3925                   format(
3926                       "The gradient of atom %s is (%8.3f,%8.3f,%8.3f).", a.toString(), gx, gy, gz));
3927           double[] vals = new double[3];
3928           a.getVelocity(vals);
3929           sb.append(format("\n Velocities: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
3930           a.getAcceleration(vals);
3931           sb.append(format("\n Accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
3932           a.getPreviousAcceleration(vals);
3933           sb.append(
3934               format("\n Previous accelerations: %8.3g %8.3g %8.3g", vals[0], vals[1], vals[2]));
3935 
3936           throw new EnergyException(sb.toString());
3937         }
3938         g[index++] = gx;
3939         g[index++] = gy;
3940         g[index++] = gz;
3941       }
3942     }
3943     return g;
3944   }
3945 
3946   /**
3947    * setRestraintBond
3948    *
3949    * @param a1 a {@link ffx.potential.bonded.Atom} object.
3950    * @param a2 a {@link ffx.potential.bonded.Atom} object.
3951    * @param distance a double.
3952    * @param forceConstant the force constant in kcal/mole.
3953    * @param flatBottom Radius of a flat-bottom potential in Angstroms.
3954    * @param lamStart At what lambda does the restraint begin to take effect?
3955    * @param lamEnd At what lambda does the restraint hit full strength?
3956    * @param switchingFunction Switching function to use as a lambda dependence.
3957    */
3958   private void setRestraintBond(
3959       Atom a1,
3960       Atom a2,
3961       double distance,
3962       double forceConstant,
3963       double flatBottom,
3964       double lamStart,
3965       double lamEnd,
3966       UnivariateSwitchingFunction switchingFunction) {
3967     restraintBondTerm = true;
3968     boolean rbLambda = !(switchingFunction instanceof ConstantSwitch) && lambdaTerm;
3969     RestraintBond rb =
3970         new RestraintBond(a1, a2, crystal, rbLambda, lamStart, lamEnd, switchingFunction);
3971     int[] classes = {a1.getAtomType().atomClass, a2.getAtomType().atomClass};
3972     if (flatBottom != 0) {
3973       rb.setBondType(
3974           new BondType(
3975               classes,
3976               forceConstant,
3977               distance,
3978               BondType.BondFunction.FLAT_BOTTOM_HARMONIC,
3979               flatBottom));
3980     } else {
3981       rb.setBondType(
3982           (new BondType(classes, forceConstant, distance, BondType.BondFunction.HARMONIC)));
3983     }
3984 
3985     // As long as we continue to add elements one-at-a-time to an array, this code will continue to
3986     // be ugly.
3987     RestraintBond[] newRbs = new RestraintBond[++nRestraintBonds];
3988     if (restraintBonds != null && restraintBonds.length != 0) {
3989       System.arraycopy(restraintBonds, 0, newRbs, 0, (nRestraintBonds - 1));
3990     }
3991     newRbs[nRestraintBonds - 1] = rb;
3992     restraintBonds = newRbs;
3993     rb.energy(false);
3994     rb.log();
3995   }
3996 
3997   private Crystal configureNCS(ForceField forceField, Crystal unitCell) {
3998     // MTRIXn to be permuted with standard space group in NCSCrystal.java for experimental
3999     // refinement.
4000     if (forceField.getProperties().containsKey("MTRIX1")
4001         && forceField.getProperties().containsKey("MTRIX2")
4002         && forceField.getProperties().containsKey("MTRIX3")) {
4003       Crystal unitCell2 =
4004           new Crystal(
4005               unitCell.a,
4006               unitCell.b,
4007               unitCell.c,
4008               unitCell.alpha,
4009               unitCell.beta,
4010               unitCell.gamma,
4011               unitCell.spaceGroup.pdbName);
4012       SpaceGroup spaceGroup = unitCell2.spaceGroup;
4013       // Separate string list MTRIXn into Double matricies then pass into symops
4014       CompositeConfiguration properties = forceField.getProperties();
4015       String[] MTRX1List = properties.getStringArray("MTRIX1");
4016       String[] MTRX2List = properties.getStringArray("MTRIX2");
4017       String[] MTRX3List = properties.getStringArray("MTRIX3");
4018       spaceGroup.symOps.clear();
4019       double number1;
4020       double number2;
4021       double number3;
4022       for (int i = 0; i < MTRX1List.length; i++) {
4023         double[][] Rot_MTRX = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
4024         double[] Tr_MTRX = {0, 0, 0};
4025         String[] tokens1 = MTRX1List[i].trim().split(" +"); // 4 items: rot [0][1-3] * trans[0]
4026         String[] tokens2 = MTRX2List[i].trim().split(" +"); // 4 items: rot [1][1-3] * trans[1]
4027         String[] tokens3 = MTRX3List[i].trim().split(" +"); // 4 items: rot [2][1-3] * trans[2]
4028         for (int k = 0; k < 4; k++) {
4029           number1 = Double.parseDouble(tokens1[k]);
4030           number2 = Double.parseDouble(tokens2[k]);
4031           number3 = Double.parseDouble(tokens3[k]);
4032           if (k != 3) {
4033             Rot_MTRX[0][k] = number1;
4034             Rot_MTRX[1][k] = number2;
4035             Rot_MTRX[2][k] = number3;
4036           } else {
4037             Tr_MTRX[0] = number1;
4038             Tr_MTRX[1] = number2;
4039             Tr_MTRX[2] = number3;
4040           }
4041         }
4042         SymOp symOp = new SymOp(Rot_MTRX, Tr_MTRX);
4043         if (logger.isLoggable(Level.FINEST)) {
4044           logger.info(
4045               format(" MTRIXn SymOp: %d of %d\n" + symOp.toString(), i + 1, MTRX1List.length));
4046         }
4047         spaceGroup.symOps.add(symOp);
4048       }
4049       unitCell = NCSCrystal.NCSCrystalFactory(unitCell, spaceGroup.symOps);
4050       unitCell.updateCrystal();
4051     }
4052     return unitCell;
4053   }
4054 
4055   /**
4056    * Getter for the field <code>restraintBonds</code>.
4057    *
4058    * @return a {@link java.util.List} object.
4059    */
4060   protected List<RestraintBond> getRestraintBonds() {
4061     if (restraintBonds != null && restraintBonds.length > 0) {
4062       return Arrays.asList(restraintBonds);
4063     } else {
4064       return null;
4065     }
4066   }
4067 
4068   /**
4069    * Getter for the RestrainGroup field.
4070    *
4071    * @return Returns RestrainGroup.
4072    */
4073   protected RestrainGroups getRestrainGroups() {
4074     return restrainGroups;
4075   }
4076 
4077   /**
4078    * Platform describes a set of force field implementations that include a pure Java reference
4079    * implementation (FFX), and two OpenMM implementations (OMM_CUDA and OMM_REF are supported)
4080    *
4081    * <p>FFX: reference FFX implementation
4082    *
4083    * <p>OMM: Currently an alias for OMM_CUDA, may eventually become "try to find best OpenMM
4084    * implementation" OMM_CUDA:
4085    *
4086    * <p>OpenMM CUDA implementation OMM_REF: OpenMM reference implementation
4087    *
4088    * <p>OMM_OPTCPU: Optimized OpenMM CPU implementation (no AMOEBA)
4089    *
4090    * <p>OMM_OPENCL: OpenMM OpenCL implementation (no AMOEBA)
4091    */
4092   public enum Platform {
4093     FFX,
4094     OMM,
4095     OMM_CUDA,
4096     OMM_REF,
4097     OMM_OPTCPU,
4098     OMM_OPENCL
4099   }
4100 
4101   private class BondedRegion extends ParallelRegion {
4102 
4103     // Shared RMSD variables.
4104     private final SharedDouble sharedBondRMSD;
4105     private final SharedDouble sharedAngleRMSD;
4106     // Shared force field bonded energy accumulation variables.
4107     private final SharedDouble sharedBondEnergy;
4108     private final SharedDouble sharedAngleEnergy;
4109     private final SharedDouble sharedOutOfPlaneBendEnergy;
4110     private final SharedDouble sharedPiOrbitalTorsionEnergy;
4111     private final SharedDouble sharedStretchBendEnergy;
4112     private final SharedDouble sharedUreyBradleyEnergy;
4113     private final SharedDouble sharedImproperTorsionEnergy;
4114     private final SharedDouble sharedTorsionEnergy;
4115     private final SharedDouble sharedStretchTorsionEnergy;
4116     private final SharedDouble sharedAngleTorsionEnergy;
4117     private final SharedDouble sharedTorsionTorsionEnergy;
4118     // Shared restraint terms.
4119     private final SharedDouble sharedRestraintBondEnergy;
4120     private final SharedDouble sharedRestTorsEnergy;
4121     // Number of threads.
4122     private final int nThreads;
4123     // Gradient loops.
4124     private final GradInitLoop[] gradInitLoops;
4125     private final GradReduceLoop[] gradReduceLoops;
4126     // Force field bonded energy parallel loops.
4127     private final BondedTermLoop[] bondLoops;
4128     private final BondedTermLoop[] angleLoops;
4129     private final BondedTermLoop[] outOfPlaneBendLoops;
4130     private final BondedTermLoop[] improperTorsionLoops;
4131     private final BondedTermLoop[] piOrbitalTorsionLoops;
4132     private final BondedTermLoop[] stretchBendLoops;
4133     private final BondedTermLoop[] torsionLoops;
4134     private final BondedTermLoop[] stretchTorsionLoops;
4135     private final BondedTermLoop[] angleTorsionLoops;
4136     private final BondedTermLoop[] torsionTorsionLoops;
4137     private final BondedTermLoop[] ureyBradleyLoops;
4138     // Retraint energy parallel loops.
4139     private final BondedTermLoop[] restraintBondLoops;
4140     private final BondedTermLoop[] rTorsLoops;
4141     private final AtomicDoubleArray3D grad;
4142     // Flag to indicate gradient computation.
4143     private boolean gradient = false;
4144     private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
4145     private AtomicDoubleArray3D lambdaGrad;
4146 
4147     BondedRegion() {
4148 
4149       // Allocate shared RMSD variables.
4150       sharedAngleRMSD = new SharedDouble();
4151       sharedBondRMSD = new SharedDouble();
4152 
4153       // Allocate shared force field bonded energy variables.
4154       sharedBondEnergy = new SharedDouble();
4155       sharedAngleEnergy = new SharedDouble();
4156       sharedImproperTorsionEnergy = new SharedDouble();
4157       sharedOutOfPlaneBendEnergy = new SharedDouble();
4158       sharedPiOrbitalTorsionEnergy = new SharedDouble();
4159       sharedStretchBendEnergy = new SharedDouble();
4160       sharedTorsionEnergy = new SharedDouble();
4161       sharedStretchTorsionEnergy = new SharedDouble();
4162       sharedAngleTorsionEnergy = new SharedDouble();
4163       sharedTorsionTorsionEnergy = new SharedDouble();
4164       sharedUreyBradleyEnergy = new SharedDouble();
4165       sharedRestTorsEnergy = new SharedDouble();
4166 
4167       // Allocate shared restraint variables.
4168       sharedRestraintBondEnergy = new SharedDouble();
4169 
4170       nThreads = parallelTeam.getThreadCount();
4171 
4172       // Gradient initialization loops.
4173       gradInitLoops = new GradInitLoop[nThreads];
4174       gradReduceLoops = new GradReduceLoop[nThreads];
4175 
4176       // Allocate memory for force field bonded energy loop arrays.
4177       angleLoops = new BondedTermLoop[nThreads];
4178       bondLoops = new BondedTermLoop[nThreads];
4179       improperTorsionLoops = new BondedTermLoop[nThreads];
4180       outOfPlaneBendLoops = new BondedTermLoop[nThreads];
4181       piOrbitalTorsionLoops = new BondedTermLoop[nThreads];
4182       stretchBendLoops = new BondedTermLoop[nThreads];
4183       torsionLoops = new BondedTermLoop[nThreads];
4184       stretchTorsionLoops = new BondedTermLoop[nThreads];
4185       angleTorsionLoops = new BondedTermLoop[nThreads];
4186       torsionTorsionLoops = new BondedTermLoop[nThreads];
4187       ureyBradleyLoops = new BondedTermLoop[nThreads];
4188       rTorsLoops = new BondedTermLoop[nThreads];
4189 
4190       // Allocate memory for restrain energy terms.
4191       restraintBondLoops = new BondedTermLoop[nThreads];
4192 
4193       // Define how the gradient will be accumulated.
4194       atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
4195       ForceField forceField = molecularAssembly.getForceField();
4196 
4197       String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
4198       try {
4199         atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
4200       } catch (Exception e) {
4201         logger.info(
4202             format(
4203                 " Unrecognized ARRAY-REDUCTION %s; defaulting to %s",
4204                 value, atomicDoubleArrayImpl));
4205       }
4206       logger.fine(format("  Bonded using %s arrays.", atomicDoubleArrayImpl.toString()));
4207 
4208       grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, nThreads);
4209       if (lambdaTerm) {
4210         lambdaGrad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, nThreads);
4211       }
4212     }
4213 
4214     @Override
4215     public void finish() {
4216       // Finalize bond and angle RMSD values.
4217       if (bondTerm) {
4218         bondRMSD = sqrt(sharedBondRMSD.get() / bonds.length);
4219       }
4220       if (angleTerm) {
4221         angleRMSD = sqrt(sharedAngleRMSD.get() / angles.length);
4222       }
4223 
4224       // Load shared energy values into their respective fields.
4225       angleEnergy = sharedAngleEnergy.get();
4226       bondEnergy = sharedBondEnergy.get();
4227       improperTorsionEnergy = sharedImproperTorsionEnergy.get();
4228       outOfPlaneBendEnergy = sharedOutOfPlaneBendEnergy.get();
4229       piOrbitalTorsionEnergy = sharedPiOrbitalTorsionEnergy.get();
4230       stretchBendEnergy = sharedStretchBendEnergy.get();
4231       ureyBradleyEnergy = sharedUreyBradleyEnergy.get();
4232       torsionEnergy = sharedTorsionEnergy.get();
4233       stretchTorsionEnergy = sharedStretchTorsionEnergy.get();
4234       angleTorsionEnergy = sharedAngleTorsionEnergy.get();
4235       torsionTorsionEnergy = sharedTorsionTorsionEnergy.get();
4236       ureyBradleyEnergy = sharedUreyBradleyEnergy.get();
4237 
4238       // Load shared restraint energy values.
4239       restraintBondEnergy = sharedRestraintBondEnergy.get();
4240       rTorsEnergy = sharedRestTorsEnergy.get();
4241     }
4242 
4243     @Override
4244     public void run() throws Exception {
4245       int threadID = getThreadIndex();
4246 
4247       // Initialize the Gradient to zero.
4248       if (gradient) {
4249         if (gradInitLoops[threadID] == null) {
4250           gradInitLoops[threadID] = new GradInitLoop();
4251         }
4252         execute(0, nAtoms - 1, gradInitLoops[threadID]);
4253       }
4254 
4255       // Evaluate force field bonded energy terms in parallel.
4256       if (angleTerm) {
4257         if (angleLoops[threadID] == null) {
4258           angleLoops[threadID] = new BondedTermLoop(angles, sharedAngleEnergy, sharedAngleRMSD);
4259         }
4260         if (threadID == 0) {
4261           angleTime = -System.nanoTime();
4262         }
4263         execute(0, nAngles - 1, angleLoops[threadID]);
4264         if (threadID == 0) {
4265           angleTime += System.nanoTime();
4266         }
4267       }
4268 
4269       if (bondTerm) {
4270         if (bondLoops[threadID] == null) {
4271           bondLoops[threadID] = new BondedTermLoop(bonds, sharedBondEnergy, sharedBondRMSD);
4272         }
4273         if (threadID == 0) {
4274           bondTime = -System.nanoTime();
4275         }
4276         execute(0, nBonds - 1, bondLoops[threadID]);
4277         if (threadID == 0) {
4278           bondTime += System.nanoTime();
4279         }
4280       }
4281 
4282       if (improperTorsionTerm) {
4283         if (improperTorsionLoops[threadID] == null) {
4284           improperTorsionLoops[threadID] =
4285               new BondedTermLoop(improperTorsions, sharedImproperTorsionEnergy);
4286         }
4287         if (threadID == 0) {
4288           improperTorsionTime = -System.nanoTime();
4289         }
4290         execute(0, nImproperTorsions - 1, improperTorsionLoops[threadID]);
4291         if (threadID == 0) {
4292           improperTorsionTime += System.nanoTime();
4293         }
4294       }
4295 
4296       if (outOfPlaneBendTerm) {
4297         if (outOfPlaneBendLoops[threadID] == null) {
4298           outOfPlaneBendLoops[threadID] =
4299               new BondedTermLoop(outOfPlaneBends, sharedOutOfPlaneBendEnergy);
4300         }
4301         if (threadID == 0) {
4302           outOfPlaneBendTime = -System.nanoTime();
4303         }
4304         execute(0, nOutOfPlaneBends - 1, outOfPlaneBendLoops[threadID]);
4305         if (threadID == 0) {
4306           outOfPlaneBendTime += System.nanoTime();
4307         }
4308       }
4309 
4310       if (piOrbitalTorsionTerm) {
4311         if (piOrbitalTorsionLoops[threadID] == null) {
4312           piOrbitalTorsionLoops[threadID] =
4313               new BondedTermLoop(piOrbitalTorsions, sharedPiOrbitalTorsionEnergy);
4314         }
4315         if (threadID == 0) {
4316           piOrbitalTorsionTime = -System.nanoTime();
4317         }
4318         execute(0, nPiOrbitalTorsions - 1, piOrbitalTorsionLoops[threadID]);
4319         if (threadID == 0) {
4320           piOrbitalTorsionTime += System.nanoTime();
4321         }
4322       }
4323 
4324       if (stretchBendTerm) {
4325         if (stretchBendLoops[threadID] == null) {
4326           stretchBendLoops[threadID] = new BondedTermLoop(stretchBends, sharedStretchBendEnergy);
4327         }
4328         if (threadID == 0) {
4329           stretchBendTime = -System.nanoTime();
4330         }
4331         execute(0, nStretchBends - 1, stretchBendLoops[threadID]);
4332         if (threadID == 0) {
4333           stretchBendTime += System.nanoTime();
4334         }
4335       }
4336 
4337       if (torsionTerm) {
4338         if (torsionLoops[threadID] == null) {
4339           torsionLoops[threadID] = new BondedTermLoop(torsions, sharedTorsionEnergy);
4340         }
4341         if (threadID == 0) {
4342           torsionTime = -System.nanoTime();
4343         }
4344         execute(0, nTorsions - 1, torsionLoops[threadID]);
4345         if (threadID == 0) {
4346           torsionTime += System.nanoTime();
4347         }
4348       }
4349 
4350       if (stretchTorsionTerm) {
4351         if (stretchTorsionLoops[threadID] == null) {
4352           stretchTorsionLoops[threadID] =
4353               new BondedTermLoop(stretchTorsions, sharedStretchTorsionEnergy);
4354         }
4355         if (threadID == 0) {
4356           stretchTorsionTime = -System.nanoTime();
4357         }
4358         execute(0, nStretchTorsions - 1, stretchTorsionLoops[threadID]);
4359         if (threadID == 0) {
4360           stretchTorsionTime += System.nanoTime();
4361         }
4362       }
4363 
4364       if (angleTorsionTerm) {
4365         if (angleTorsionLoops[threadID] == null) {
4366           angleTorsionLoops[threadID] = new BondedTermLoop(angleTorsions, sharedAngleTorsionEnergy);
4367         }
4368         if (threadID == 0) {
4369           angleTorsionTime = -System.nanoTime();
4370         }
4371         execute(0, nAngleTorsions - 1, angleTorsionLoops[threadID]);
4372         if (threadID == 0) {
4373           angleTorsionTime += System.nanoTime();
4374         }
4375       }
4376 
4377       if (torsionTorsionTerm) {
4378         if (torsionTorsionLoops[threadID] == null) {
4379           torsionTorsionLoops[threadID] =
4380               new BondedTermLoop(torsionTorsions, sharedTorsionTorsionEnergy);
4381         }
4382         if (threadID == 0) {
4383           torsionTorsionTime = -System.nanoTime();
4384         }
4385         execute(0, nTorsionTorsions - 1, torsionTorsionLoops[threadID]);
4386         if (threadID == 0) {
4387           torsionTorsionTime += System.nanoTime();
4388         }
4389       }
4390 
4391       if (ureyBradleyTerm) {
4392         if (ureyBradleyLoops[threadID] == null) {
4393           ureyBradleyLoops[threadID] = new BondedTermLoop(ureyBradleys, sharedUreyBradleyEnergy);
4394         }
4395         if (threadID == 0) {
4396           ureyBradleyTime = -System.nanoTime();
4397         }
4398         execute(0, nUreyBradleys - 1, ureyBradleyLoops[threadID]);
4399         if (threadID == 0) {
4400           ureyBradleyTime += System.nanoTime();
4401         }
4402       }
4403 
4404       // Evaluate restraint terms in parallel.
4405       if (restraintBondTerm) {
4406         if (restraintBondLoops[threadID] == null) {
4407           restraintBondLoops[threadID] =
4408               new BondedTermLoop(restraintBonds, sharedRestraintBondEnergy);
4409         }
4410         if (threadID == 0) {
4411           restraintBondTime = -System.nanoTime();
4412         }
4413         execute(0, nRestraintBonds - 1, restraintBondLoops[threadID]);
4414         if (threadID == 0) {
4415           restraintBondTime += System.nanoTime();
4416         }
4417       }
4418 
4419       if (rTorsTerm) {
4420         if (rTorsLoops[threadID] == null) {
4421           rTorsLoops[threadID] = new BondedTermLoop(rTors, sharedRestTorsEnergy);
4422         }
4423         if (threadID == 0) {
4424           rTorsTime = -System.nanoTime();
4425         }
4426         execute(0, nRestTors - 1, rTorsLoops[threadID]);
4427         if (threadID == 0) {
4428           rTorsTime += System.nanoTime();
4429         }
4430       }
4431 
4432       // Reduce the Gradient and load it into Atom instances.
4433       if (gradient) {
4434         if (gradReduceLoops[threadID] == null) {
4435           gradReduceLoops[threadID] = new GradReduceLoop();
4436         }
4437         execute(0, nAtoms - 1, gradReduceLoops[threadID]);
4438       }
4439     }
4440 
4441     public void setGradient(boolean gradient) {
4442       this.gradient = gradient;
4443     }
4444 
4445     @Override
4446     public void start() {
4447       // Zero out shared RMSD values.
4448       sharedAngleRMSD.set(0.0);
4449       sharedBondRMSD.set(0.0);
4450 
4451       // Zero out shared energy values.
4452       sharedAngleEnergy.set(0.0);
4453       sharedBondEnergy.set(0.0);
4454       sharedImproperTorsionEnergy.set(0.0);
4455       sharedOutOfPlaneBendEnergy.set(0.0);
4456       sharedPiOrbitalTorsionEnergy.set(0.0);
4457       sharedStretchBendEnergy.set(0.0);
4458       sharedTorsionEnergy.set(0.0);
4459       sharedStretchTorsionEnergy.set(0.0);
4460       sharedAngleTorsionEnergy.set(0.0);
4461       sharedTorsionTorsionEnergy.set(0.0);
4462       sharedUreyBradleyEnergy.set(0.0);
4463 
4464       // Zero out shared restraint energy values.
4465       sharedRestraintBondEnergy.set(0.0);
4466       sharedRestTorsEnergy.set(0.0);
4467 
4468       // Assure capacity of the gradient arrays.
4469       if (gradient) {
4470         grad.alloc(nAtoms);
4471       }
4472       if (lambdaTerm) {
4473         lambdaGrad.alloc(nAtoms);
4474       }
4475     }
4476 
4477     private class GradInitLoop extends IntegerForLoop {
4478 
4479       @Override
4480       public void run(int first, int last) throws Exception {
4481         int threadID = getThreadIndex();
4482         if (gradient) {
4483           grad.reset(threadID, first, last);
4484           for (int i = first; i <= last; i++) {
4485             atoms[i].setXYZGradient(0.0, 0.0, 0.0);
4486           }
4487         }
4488         if (lambdaTerm) {
4489           lambdaGrad.reset(threadID, first, last);
4490           for (int i = first; i <= last; i++) {
4491             atoms[i].setLambdaXYZGradient(0.0, 0.0, 0.0);
4492           }
4493         }
4494       }
4495 
4496       @Override
4497       public IntegerSchedule schedule() {
4498         return IntegerSchedule.fixed();
4499       }
4500     }
4501 
4502     private class GradReduceLoop extends IntegerForLoop {
4503 
4504       @Override
4505       public void run(int first, int last) throws Exception {
4506         if (gradient) {
4507           grad.reduce(first, last);
4508           for (int i = first; i <= last; i++) {
4509             Atom a = atoms[i];
4510             a.setXYZGradient(grad.getX(i), grad.getY(i), grad.getZ(i));
4511           }
4512         }
4513         if (lambdaTerm) {
4514           lambdaGrad.reduce(first, last);
4515           for (int i = first; i <= last; i++) {
4516             Atom a = atoms[i];
4517             a.setLambdaXYZGradient(lambdaGrad.getX(i), lambdaGrad.getY(i), lambdaGrad.getZ(i));
4518           }
4519         }
4520       }
4521 
4522       @Override
4523       public IntegerSchedule schedule() {
4524         return IntegerSchedule.fixed();
4525       }
4526     }
4527 
4528     private class BondedTermLoop extends IntegerForLoop {
4529 
4530       private final BondedTerm[] terms;
4531       private final SharedDouble sharedEnergy;
4532       private final SharedDouble sharedRMSD;
4533       private final boolean computeRMSD;
4534       private double localEnergy;
4535       private double localRMSD;
4536       private int threadID;
4537 
4538       BondedTermLoop(BondedTerm[] terms, SharedDouble sharedEnergy) {
4539         this(terms, sharedEnergy, null);
4540       }
4541 
4542       BondedTermLoop(BondedTerm[] terms, SharedDouble sharedEnergy, SharedDouble sharedRMSD) {
4543         this.terms = terms;
4544         this.sharedEnergy = sharedEnergy;
4545         this.sharedRMSD = sharedRMSD;
4546         computeRMSD = (sharedRMSD != null);
4547       }
4548 
4549       @Override
4550       public void finish() {
4551         sharedEnergy.addAndGet(localEnergy);
4552         if (computeRMSD) {
4553           sharedRMSD.addAndGet(localRMSD);
4554         }
4555       }
4556 
4557       @Override
4558       public void run(int first, int last) throws Exception {
4559         for (int i = first; i <= last; i++) {
4560           BondedTerm term = terms[i];
4561           /*
4562            * The logic here deals with how DualTopologyEnergy (DTE) works.
4563            * If this ForceFieldEnergy (FFE) is not part of a DTE, lambdaBondedTerms is always false, and the term is always evaluated.
4564            *
4565            * Most bonded terms should be at full-strength, regardless of lambda, "complemented" to 1.0*strength.
4566            * Outside FFE, DTE scales the initial, !lambdaBondedTerms evaluation by f(lambda).
4567            *
4568            * In the case of a bonded term lacking any softcore atoms, this will be externally complemented by the other FFE topology.
4569            * If there is internal lambda-scaling, that will separately apply at both ends.
4570            *
4571            * In the case of a bonded term with a softcore atom, it's a bit trickier.
4572            * If it's unscaled by lambda, it needs to be internally complemented; we re-evaluate it with lambdaBondedTerms true.
4573            * This second evaluation is scaled by DTE by a factor of f(1-lambda), and becomes "restraintEnergy".
4574            * If it is scaled internally by lambda, we assume that the energy term is not meant to be internally complemented.
4575            * In that case, we skip evaluation into restraintEnergy.
4576            */
4577           boolean used = !lambdaBondedTerms || lambdaAllBondedTerms || (term.applyLambda()
4578               && !term.isLambdaScaled());
4579           if (used) {
4580             localEnergy += term.energy(gradient, threadID, grad, lambdaGrad);
4581             if (computeRMSD) {
4582               double value = term.getValue();
4583               localRMSD += value * value;
4584             }
4585           }
4586         }
4587       }
4588 
4589       @Override
4590       public void start() {
4591         localEnergy = 0.0;
4592         localRMSD = 0.0;
4593         threadID = getThreadIndex();
4594       }
4595     }
4596   }
4597 }