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-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.potential.nonbonded;
39  
40  import edu.rit.pj.IntegerSchedule;
41  import edu.rit.pj.ParallelTeam;
42  import edu.rit.pj.reduction.SharedDouble;
43  import edu.rit.util.Range;
44  import ffx.crystal.Crystal;
45  import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
46  import ffx.numerics.atomic.AtomicDoubleArray3D;
47  import ffx.numerics.multipole.MultipoleTensor;
48  import ffx.potential.Platform;
49  import ffx.potential.bonded.Atom;
50  import ffx.potential.bonded.Bond;
51  import ffx.potential.bonded.LambdaInterface;
52  import ffx.potential.extended.ExtendedSystem;
53  import ffx.potential.nonbonded.pme.*;
54  import ffx.potential.parameters.AtomType;
55  import ffx.potential.parameters.ForceField;
56  import ffx.potential.parameters.ForceField.ELEC_FORM;
57  import ffx.potential.parameters.ForceField.ForceFieldType;
58  import ffx.potential.parameters.MultipoleType;
59  import ffx.potential.parameters.PolarizeType;
60  import ffx.potential.utils.EnergyException;
61  import ffx.utilities.Constants;
62  import ffx.utilities.FFXProperty;
63  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
64  import org.apache.commons.math3.linear.EigenDecomposition;
65  
66  import java.util.ArrayList;
67  import java.util.Arrays;
68  import java.util.List;
69  import java.util.logging.Level;
70  import java.util.logging.Logger;
71  
72  import static ffx.potential.nonbonded.pme.EwaldParameters.DEFAULT_EWALD_COEFFICIENT;
73  import static ffx.potential.parameters.ForceField.ELEC_FORM.PAM;
74  import static ffx.potential.parameters.ForceField.toEnumForm;
75  import static ffx.potential.parameters.MultipoleType.*;
76  import static ffx.utilities.Constants.ELEC_ANG_TO_DEBYE;
77  import static ffx.utilities.Constants.NS2SEC;
78  import static ffx.utilities.PropertyGroup.ElectrostaticsFunctionalForm;
79  import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
80  import static java.lang.String.format;
81  import static java.lang.System.arraycopy;
82  import static java.util.Arrays.fill;
83  import static java.util.Collections.sort;
84  import static org.apache.commons.math3.util.FastMath.*;
85  
86  /**
87   * This Particle Mesh Ewald class implements PME for the AMOEBA polarizable mutlipole force field in
88   * parallel using a {@link NeighborList} for any {@link Crystal} space group. The real space
89   * contribution is contained within this class and the reciprocal space contribution is delegated to
90   * the {@link ReciprocalSpace} class.
91   *
92   * @author Michael J. Schnieders<br> derived from:<br> TINKER code by Jay Ponder, Pengyu Ren and Tom
93   * Darden.<br>
94   * @see <br>
95   * <a href="http://dx.doi.org/10.1021/ct300035u" target="_blank"> M. J. Schnieders, J.
96   * Baltrusaitis, Y. Shi, G. Chattree, L. Zheng, W. Yang and P. Ren, The Structure,
97   * Thermodynamics, and Solubility of Organic Crystals from Simulation with a Polarizable Force
98   * Field, Journal of Chemical Theory and Computation 8 (5), 1721-36 (2012)</a>
99   * @see <br>
100  * <a href="http://dx.doi.org/10.1021/ct100506d" target="_blank"> M. J. Schnieders, T. D. Fenn
101  * and V. S. Pande, Polarizable atomic multipole X-ray refinement: Particle-mesh Ewald
102  * electrostatics for macromolecular crystals. Journal of Chemical Theory and Computation 7 (4),
103  * 1141-56 (2011)</a>
104  * @see <br>
105  * <a href="http://dx.doi.org/10.1063/1.1630791" target="_blank"> C. Sagui, L. G. Pedersen, and
106  * T. A. Darden, Journal of Chemical Physics 120 (1), 73 (2004)</a>
107  * @see <br>
108  * <a href="http://link.aip.org/link/?JCPSA6/98/10089/1" target="_©lank"> T. Darden, D. York,
109  * and L. Pedersen, Journal of Chemical Physics 98 (12), 10089 (1993)</a>
110  * @see <br>
111  * <a href="http://www.ccp5.org" target="_blank"> W. Smith, "Point Multipoles in the Ewald
112  * Summation (Revisited)", CCP5 Newsletter, 46, 18-30 (1998)</a>
113  */
114 @SuppressWarnings("deprecation")
115 public class ParticleMeshEwald implements LambdaInterface {
116 
117   private static final Logger logger = Logger.getLogger(ParticleMeshEwald.class.getName());
118   /**
119    * Number of unique tensors for given order.
120    */
121   private static final int tensorCount = MultipoleTensor.tensorCount(3);
122   /**
123    * If lambdaTerm is true, some ligand atom interactions with the environment are being turned
124    * on/off.
125    */
126   private final boolean lambdaTerm;
127   /**
128    * If nnTerm is true, some atomic interactions are treated with a neural network.
129    */
130   private final boolean nnTerm;
131   private final boolean reciprocalSpaceTerm;
132   /**
133    * Reference to the force field being used.
134    */
135   private final ForceField forceField;
136   private static final double DEFAULT_POLAR_EPS = 1.0e-6;
137   @FFXProperty(name = "polar-eps", propertyGroup = ElectrostaticsFunctionalForm, defaultValue = "1.0e-6",
138       description = """
139           Sets the convergence criterion applied during computation of self-consistent induced dipoles.
140           The calculation is deemed to have converged when the rms change in Debyes in
141           the induced dipole at all polarizable sites is less than the value specified with this property.
142           The default value in the absence of the keyword is 1.0e-6 Debyes.
143           """)
144   private final double poleps;
145   private final boolean directFallback;
146   /**
147    * Specify an SCF predictor algorithm.
148    */
149   private final SCFPredictorParameters scfPredictorParameters;
150   private final EwaldParameters ewaldParameters;
151   private final ScaleParameters scaleParameters;
152   private final AlchemicalParameters alchemicalParameters;
153   private final RealSpaceNeighborParameters realSpaceNeighborParameters;
154   private final PMETimings pmeTimings;
155   /**
156    * By default, maxThreads is set to the number of available SMP cores.
157    */
158   private final int maxThreads;
159   /**
160    * The default ParallelTeam encapsulates the maximum number of threads used to parallelize the
161    * electrostatics calculation.
162    */
163   private final ParallelTeam parallelTeam;
164   /**
165    * The sectionTeam encapsulates 1 or 2 threads.
166    *
167    * <p>If it contains 1 thread, the real and reciprocal space calculations are done sequentially.
168    *
169    * <p>If it contains 2 threads, the real and reciprocal space calculations will be done
170    * concurrently.
171    */
172   private final ParallelTeam sectionTeam;
173   /**
174    * If the real and reciprocal space parts of PME are done sequentially, then the realSpaceTeam is
175    * equal parallalTeam.
176    *
177    * <p>If the real and reciprocal space parts of PME are done concurrently, then the realSpaceTeam
178    * will have fewer threads than the default parallelTeam.
179    */
180   private final ParallelTeam realSpaceTeam;
181   /**
182    * If the real and reciprocal space parts of PME are done sequentially, then the
183    * reciprocalSpaceTeam is equal parallalTeam.
184    *
185    * <p>If the real and reciprocal space parts of PME are done concurrently, then the
186    * reciprocalSpaceTeam will have fewer threads than the default parallelTeam.
187    */
188   private final ParallelTeam fftTeam;
189   private final NeighborList neighborList;
190   private final InitializationRegion initializationRegion;
191   private final PermanentFieldRegion permanentFieldRegion;
192   private final InducedDipoleFieldRegion inducedDipoleFieldRegion;
193   private final InducedDipoleFieldReduceRegion inducedDipoleFieldReduceRegion;
194   private final ExpandInducedDipolesRegion expandInducedDipolesRegion;
195   private final DirectRegion directRegion;
196   private final SORRegion sorRegion;
197   private final OPTRegion optRegion;
198   private final PCGSolver pcgSolver;
199   private final ReciprocalSpace reciprocalSpace;
200   private final ReciprocalEnergyRegion reciprocalEnergyRegion;
201   private final RealSpaceEnergyRegion realSpaceEnergyRegion;
202   private final ReduceRegion reduceRegion;
203   private final GeneralizedKirkwood generalizedKirkwood;
204   /**
205    * Polarization modes include "direct", in which induced dipoles do not interact, and "mutual" that
206    * converges the self-consistent field to a tolerance specified by the "polar-eps" keyword.
207    */
208   public Polarization polarization;
209   /**
210    * Dimensions of [nsymm][xyz][nAtoms].
211    */
212   public double[][][] coordinates;
213   /**
214    * Neighbor lists, including atoms beyond the real space cutoff. [nsymm][nAtoms][nAllNeighbors]
215    */
216   public int[][][] neighborLists;
217   /**
218    * Cartesian multipoles in the global frame with dimensions of [nsymm][nAtoms][10]
219    */
220   public double[][][] globalMultipole;
221   /**
222    * Fractional multipoles in the global frame with dimensions of [nsymm][nAtoms][10]
223    */
224   public double[][][] fractionalMultipole;
225   /**
226    * Dimensions of [nsymm][nAtoms][3]
227    */
228   public double[][][] inducedDipole;
229   public double[][][] inducedDipoleCR;
230   /**
231    * Direct induced dipoles.
232    */
233   public double[][] directDipole;
234   public double[][] directDipoleCR;
235   public double[][] directField;
236   public double[][] directFieldCR;
237   /**
238    * Vacuum induced dipoles
239    */
240   public double[][][] vacuumInducedDipole;
241   public double[][][] vacuumInducedDipoleCR;
242   /**
243    * Vacuum induced dipoles
244    */
245   public double[][] vacuumDirectDipole;
246   public double[][] vacuumDirectDipoleCR;
247   /**
248    * Coulomb constant in units of kcal*Ang/(mol*electron^2)
249    */
250   @FFXProperty(name = "electric", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "332.063713",
251       description = """
252           Specifies a value for the so-called "electric constant" allowing conversion unit of electrostatic
253           potential energy values from electrons^2/Angstrom to kcal/mol. Internally, FFX stores a default value
254           for this constant as 332.063713 based on CODATA reference values. Since different force fields are
255           intended for use with slightly different values, this keyword allows overriding the default value.
256           """)
257   public double electric;
258   /**
259    * The requested permittivity for the solute.
260    */
261   public double soluteDielectric;
262   /**
263    * An ordered array of atoms in the system.
264    */
265   protected Atom[] atoms;
266   /**
267    * The number of atoms in the system.
268    */
269   protected int nAtoms;
270   /**
271    * Polarization groups.
272    */
273   protected int[][] ip11;
274   protected int[][] ip12;
275   protected int[][] ip13;
276   /**
277    * Total multipole energy = permanentMultipoleEnergy + polarizationEnergy. <br> This does not
278    * include GK.
279    */
280   protected double totalMultipoleEnergy;
281   /**
282    * Permanent multipole energy = permanentRealSpaceEnergy + permanentSelfEnergy +
283    * permanentReciprocalEnergy.
284    */
285   protected double permanentMultipoleEnergy;
286   protected double permanentRealSpaceEnergy;
287   protected double permanentSelfEnergy;
288   protected double permanentReciprocalEnergy;
289   /**
290    * Polarization energy = inducedRealSpaceEnergy + inducedSelfEnergy + inducedReciprocalEnergy.
291    */
292   protected double polarizationEnergy;
293   protected double inducedRealSpaceEnergy;
294   protected double inducedSelfEnergy;
295   protected double inducedReciprocalEnergy;
296   protected SCFAlgorithm scfAlgorithm;
297   /**
298    * Partial derivative with respect to Lambda.
299    */
300   private final SharedDouble shareddEdLambda;
301   /**
302    * Second partial derivative with respect to Lambda.
303    */
304   private final SharedDouble sharedd2EdLambda2;
305   /**
306    * The electrostatics functional form in use.
307    */
308   private final ELEC_FORM elecForm;
309   /**
310    * Unit cell and spacegroup information.
311    */
312   private Crystal crystal;
313   /**
314    * Number of symmetry operators.
315    */
316   private int nSymm;
317   /**
318    * Flag to indicate use of generalized Kirkwood.
319    */
320   private boolean generalizedKirkwoodTerm;
321   /**
322    * If true, compute coordinate gradient.
323    */
324   private boolean gradient = false;
325   /**
326    * Number of PME multipole interactions.
327    */
328   private int interactions;
329   /**
330    * Number of generalized Kirkwood interactions.
331    */
332   private int gkInteractions;
333   /**
334    * Generalized Kirkwood energy.
335    */
336   private double solvationEnergy;
337   /**
338    * The current LambdaMode of this PME instance (or OFF for no lambda dependence).
339    */
340   private LambdaMode lambdaMode = LambdaMode.OFF;
341   /**
342    * Current state.
343    */
344   private double lambda = 1.0;
345   /**
346    * Permanent multipoles in their local frame.
347    */
348   private double[][] localMultipole;
349   private MultipoleFrameDefinition[] frame;
350   private int[][] axisAtom;
351   private double[] ipdamp;
352   private double[] thole;
353   private double[] polarizability;
354   /**
355    * 1-2, 1-3, 1-4 and 1-5 connectivity lists.
356    */
357   private int[][] mask12;
358   private int[][] mask13;
359   private int[][] mask14;
360   private int[][] mask15;
361   /**
362    * Flag for ligand atoms.
363    */
364   private boolean[] isSoft;
365   /**
366    * Molecule number for each atom.
367    */
368   private int[] molecule;
369   /**
370    * When computing the polarization energy at Lambda there are 3 pieces.
371    *
372    * <p>1.) Upol(1) = The polarization energy computed normally (i.e. system with ligand).
373    *
374    * <p>2.) Uenv = The polarization energy of the system without the ligand.
375    *
376    * <p>3.) Uligand = The polarization energy of the ligand by itself.
377    *
378    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
379    *
380    * <p>Set the "use" array to true for all atoms for part 1. Set the "use" array to true for all
381    * atoms except the ligand for part 2. Set the "use" array to true only for the ligand atoms for
382    * part 3.
383    *
384    * <p>The "use" array can also be employed to turn off atoms for computing the electrostatic
385    * energy of sub-structures.
386    */
387   private boolean[] use;
388   private IntegerSchedule permanentSchedule;
389   private double[][] cartesianMultipolePhi;
390   private double[][] fracMultipolePhi;
391   private double[][] cartesianInducedDipolePhi;
392   private double[][] cartesianInducedDipolePhiCR;
393   private double[][] fractionalInducedDipolePhi;
394   private double[][] fractionalInducedDipolePhiCR;
395   private double[][] cartesianVacuumDipolePhi;
396   private double[][] cartesianVacuumDipolePhiCR;
397   private double[][] fractionalVacuumDipolePhi;
398   private double[][] fractionalVacuumDipolePhiCR;
399   /**
400    * AtomicDoubleArray implementation to use.
401    */
402   private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
403   /**
404    * Field array for each thread. [threadID][X/Y/Z][atomID]
405    */
406   private AtomicDoubleArray3D field;
407   /**
408    * Chain rule field array for each thread. [threadID][X/Y/Z][atomID]
409    */
410   private AtomicDoubleArray3D fieldCR;
411   /**
412    * Gradient array for each thread. [threadID][X/Y/Z][atomID]
413    */
414   private AtomicDoubleArray3D grad;
415   /**
416    * Torque array for each thread. [threadID][X/Y/Z][atomID]
417    */
418   private AtomicDoubleArray3D torque;
419   /**
420    * Partial derivative of the gradient with respect to Lambda. [threadID][X/Y/Z][atomID]
421    */
422   private AtomicDoubleArray3D lambdaGrad;
423   /**
424    * Partial derivative of the torque with respect to Lambda. [threadID][X/Y/Z][atomID]
425    */
426   private AtomicDoubleArray3D lambdaTorque;
427 
428   /**
429    * ParticleMeshEwald constructor.
430    *
431    * @param atoms        the Atom array to do electrostatics on.
432    * @param molecule     the molecule number for each atom.
433    * @param forceField   the ForceField the defines the electrostatics parameters.
434    * @param crystal      The boundary conditions.
435    * @param elecForm     The electrostatics functional form.
436    * @param neighborList The NeighborList for both van der Waals and PME.
437    * @param ewaldCutoff  The Ewald real space cutoff.
438    * @param gkCutoff     The generalized Kirkwood cutoff.
439    * @param parallelTeam A ParallelTeam that delegates parallelization.
440    */
441   public ParticleMeshEwald(
442       Atom[] atoms,
443       int[] molecule,
444       ForceField forceField,
445       Crystal crystal,
446       NeighborList neighborList,
447       ELEC_FORM elecForm,
448       double ewaldCutoff,
449       double gkCutoff,
450       ParallelTeam parallelTeam) {
451     this.atoms = atoms;
452     this.molecule = molecule;
453     this.forceField = forceField;
454     this.crystal = crystal;
455     this.parallelTeam = parallelTeam;
456     this.neighborList = neighborList;
457     this.elecForm = elecForm;
458     neighborLists = neighborList.getNeighborList();
459     permanentSchedule = neighborList.getPairwiseSchedule();
460     nAtoms = atoms.length;
461     nSymm = crystal.spaceGroup.getNumberOfSymOps();
462     maxThreads = parallelTeam.getThreadCount();
463 
464     electric = forceField.getDouble("ELECTRIC", Constants.DEFAULT_ELECTRIC);
465 
466     // Solute dielectric is ignored if it's less than 1.0.
467     soluteDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
468     if (soluteDielectric > 1.0) {
469       electric = electric / soluteDielectric;
470     } else {
471       soluteDielectric = 1.0;
472     }
473 
474     poleps = forceField.getDouble("POLAR_EPS", DEFAULT_POLAR_EPS);
475 
476     // If PME-specific lambda term not set, default to force field-wide lambda term.
477     lambdaTerm = forceField.getBoolean("ELEC_LAMBDATERM",
478         forceField.getBoolean("LAMBDATERM", false));
479 
480     nnTerm = forceField.getBoolean("NNTERM", false);
481 
482     if (nnTerm && lambdaTerm) {
483       logger.severe(" Use a neural network potential with alchemical simulations is not yet supported.");
484     }
485 
486     double aewald = forceField.getDouble("EWALD_ALPHA", DEFAULT_EWALD_COEFFICIENT);
487     ewaldParameters = new EwaldParameters(ewaldCutoff, aewald);
488     scaleParameters = new ScaleParameters(elecForm, forceField);
489     reciprocalSpaceTerm = forceField.getBoolean("RECIPTERM", true);
490 
491     SCFPredictor scfPredictor;
492     try {
493       String predictor = forceField.getString("SCF_PREDICTOR", "NONE");
494       predictor = predictor.replaceAll("-", "_").toUpperCase();
495       scfPredictor = SCFPredictor.valueOf(predictor);
496     } catch (Exception e) {
497       scfPredictor = SCFPredictor.NONE;
498     }
499 
500     scfPredictorParameters = new SCFPredictorParameters(scfPredictor, nAtoms);
501     if (scfPredictor != SCFPredictor.NONE) {
502       scfPredictorParameters.init(forceField);
503     }
504 
505     String algorithm = forceField.getString("SCF_ALGORITHM", "CG");
506     try {
507       algorithm = algorithm.replaceAll("-", "_").toUpperCase();
508       scfAlgorithm = SCFAlgorithm.valueOf(algorithm);
509     } catch (Exception e) {
510       scfAlgorithm = SCFAlgorithm.CG;
511     }
512 
513     // Fall back to the direct
514     directFallback = forceField.getBoolean("DIRECT_SCF_FALLBACK", true);
515 
516     // Define how force arrays will be accumulated.
517     String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
518     try {
519       atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
520     } catch (Exception e) {
521       atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
522       logger.info(
523           format(
524               " Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value, atomicDoubleArrayImpl));
525     }
526     logger.fine(format("  PME using %s arrays.", atomicDoubleArrayImpl.toString()));
527 
528     if (!scfAlgorithm.isSupported(Platform.FFX)) {
529       // Can't know a-priori whether this is being constructed under an FFX or OpenMM
530       // ForceFieldEnergy, so fine logging.
531       logger.fine(format(
532               " SCF algorithm %s is not supported by FFX reference implementation; falling back to CG!",
533               scfAlgorithm));
534       scfAlgorithm = SCFAlgorithm.CG;
535     }
536 
537     pcgSolver = new PCGSolver(maxThreads, poleps, forceField, nAtoms);
538 
539     alchemicalParameters = new AlchemicalParameters(forceField, lambdaTerm, nnTerm, polarization);
540     if (nnTerm) {
541       initNN();
542     }
543 
544     String polar = forceField.getString("POLARIZATION", "MUTUAL");
545     if (elecForm == ELEC_FORM.FIXED_CHARGE) {
546       polar = "NONE";
547     }
548     boolean polarizationTerm = forceField.getBoolean("POLARIZETERM", true);
549     if (!polarizationTerm || polar.equalsIgnoreCase("NONE")) {
550       polarization = Polarization.NONE;
551     } else if (polar.equalsIgnoreCase("DIRECT")) {
552       polarization = Polarization.DIRECT;
553     } else {
554       polarization = Polarization.MUTUAL;
555     }
556 
557     if (lambdaTerm) {
558       shareddEdLambda = new SharedDouble();
559       sharedd2EdLambda2 = new SharedDouble();
560     } else {
561       shareddEdLambda = null;
562       sharedd2EdLambda2 = null;
563       lambdaGrad = null;
564       lambdaTorque = null;
565     }
566 
567     directRegion = new DirectRegion(maxThreads);
568     sorRegion = new SORRegion(maxThreads, forceField);
569     optRegion = new OPTRegion(maxThreads);
570 
571     if (logger.isLoggable(Level.INFO)) {
572       StringBuilder sb = new StringBuilder();
573       sb.append("\n Electrostatics\n");
574       sb.append(format("   Polarization:                       %8s\n", polarization.toString()));
575       if (polarization == Polarization.MUTUAL) {
576         sb.append(format("    SCF Convergence Criteria:         %8.3e\n", poleps));
577         sb.append(format("    SCF Predictor:                     %8s\n",
578             scfPredictorParameters.scfPredictor));
579         sb.append(format("    SCF Algorithm:                     %8s\n", scfAlgorithm));
580         if (scfAlgorithm == SCFAlgorithm.SOR) {
581           sb.append(format("    SOR Parameter:                     %8.3f\n", sorRegion.getSOR()));
582         } else {
583           sb.append(format("    CG Preconditioner Cut-Off:         %8.3f\n",
584               pcgSolver.getPreconditionerCutoff()));
585           sb.append(format("    CG Preconditioner Ewald Coeff.:    %8.3f\n",
586               pcgSolver.getPreconditionerEwald()));
587           sb.append(format("    CG Preconditioner Scale:           %8.3f\n",
588               pcgSolver.getPreconditionerScale()));
589           sb.append(
590               format("    CG Preconditioner Mode:     %15s\n", pcgSolver.getPreconditionerMode()));
591         }
592       }
593       if (ewaldParameters.aewald > 0.0) {
594         sb.append("   Particle-mesh Ewald\n");
595         sb.append(format("    Ewald Coefficient:                 %8.3f\n", ewaldParameters.aewald));
596         sb.append(format("    Particle Cut-Off:                  %8.3f (A)", ewaldParameters.off));
597       } else if (ewaldParameters.off != Double.POSITIVE_INFINITY) {
598         sb.append(format("    Electrostatics Cut-Off:            %8.3f (A)\n",
599             ewaldParameters.off));
600       } else {
601         sb.append("    Electrostatics Cut-Off:                NONE\n");
602       }
603 
604       logger.info(sb.toString());
605     }
606 
607     // Either 1 or 2; see description below.
608     int sectionThreads;
609     /*
610      If real and reciprocal space are done sequentially or OpenCL is used,
611      then realSpaceThreads == maxThreads. Otherwise, the number of
612      realSpaceThreads is set to ffx.realSpaceThreads.
613     */
614     int realSpaceThreads;
615     /*
616      If real and reciprocal space are done sequentially then reciprocalThreads == maxThreads.
617      If CUDA is used, reciprocalThreads == 1 Otherwise,
618      reciprocalThreads = maxThreads - realSpaceThreads
619     */
620     int reciprocalThreads;
621 
622     boolean concurrent;
623     int realThreads = 1;
624     try {
625       realThreads = forceField.getInteger("PME_REAL_THREADS");
626       if (realThreads >= maxThreads || realThreads < 1) {
627         throw new Exception("pme-real-threads must be < ffx.nt and greater than 0");
628       }
629       concurrent = true;
630     } catch (Exception e) {
631       concurrent = false;
632     }
633     if (concurrent) {
634       sectionThreads = 2;
635       realSpaceThreads = realThreads;
636       reciprocalThreads = maxThreads - realThreads;
637       sectionTeam = new ParallelTeam(sectionThreads);
638       realSpaceTeam = new ParallelTeam(realSpaceThreads);
639       fftTeam = new ParallelTeam(reciprocalThreads);
640     } else {
641       // If pme-real-threads is not defined, then do real and reciprocal space parts sequentially.
642       sectionThreads = 1;
643       sectionTeam = new ParallelTeam(sectionThreads);
644       realSpaceTeam = parallelTeam;
645       fftTeam = parallelTeam;
646     }
647 
648     realSpaceNeighborParameters = new RealSpaceNeighborParameters(maxThreads);
649     initializationRegion = new InitializationRegion(this, maxThreads, forceField);
650     expandInducedDipolesRegion = new ExpandInducedDipolesRegion(maxThreads);
651     initAtomArrays();
652 
653     /*
654      Note that we always pass on the unit cell crystal to ReciprocalSpace
655      instance even if the real space calculations require a
656      ReplicatesCrystal.
657     */
658     if (ewaldParameters.aewald > 0.0 && reciprocalSpaceTerm) {
659       reciprocalSpace = new ReciprocalSpace(this,
660           crystal.getUnitCell(), forceField, atoms,
661           ewaldParameters.aewald, fftTeam, parallelTeam);
662       reciprocalEnergyRegion = new ReciprocalEnergyRegion(maxThreads, ewaldParameters.aewald, electric);
663     } else {
664       reciprocalSpace = null;
665       reciprocalEnergyRegion = null;
666     }
667     permanentFieldRegion = new PermanentFieldRegion(realSpaceTeam, forceField, lambdaTerm);
668     inducedDipoleFieldRegion = new InducedDipoleFieldRegion(realSpaceTeam, forceField, lambdaTerm);
669     inducedDipoleFieldReduceRegion = new InducedDipoleFieldReduceRegion(maxThreads);
670     realSpaceEnergyRegion = new RealSpaceEnergyRegion(elecForm, maxThreads, forceField, lambdaTerm, electric);
671     reduceRegion = new ReduceRegion(maxThreads, forceField);
672 
673     pmeTimings = new PMETimings(maxThreads);
674 
675     if (lambdaTerm) {
676       logger.info(alchemicalParameters.toString());
677     }
678 
679     // The GK reaction field is added to the intra-molecular field to give the self-consistent
680     // reaction field.
681     generalizedKirkwoodTerm = forceField.getBoolean("GKTERM", false);
682     if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
683       generalizedKirkwood = new GeneralizedKirkwood(forceField, atoms, this, crystal, parallelTeam, gkCutoff);
684     } else {
685       generalizedKirkwood = null;
686     }
687   }
688 
689   /**
690    * Returns the ELEC_FORM.
691    */
692   public ELEC_FORM getElecForm() {
693     return elecForm;
694   }
695 
696   public void computeInduceDipoleField() {
697     expandInducedDipoles();
698 
699     if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
700       reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
701     }
702     field.reset(parallelTeam);
703     fieldCR.reset(parallelTeam);
704     inducedDipoleFieldRegion.init(
705         atoms, crystal, use, molecule, ipdamp, thole, coordinates, realSpaceNeighborParameters,
706         inducedDipole, inducedDipoleCR, reciprocalSpaceTerm, reciprocalSpace, lambdaMode,
707         ewaldParameters, field, fieldCR, pmeTimings);
708     inducedDipoleFieldRegion.executeWith(sectionTeam);
709     pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
710 
711     if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
712       reciprocalSpace.computeInducedPhi(
713           cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
714           fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
715     }
716 
717     if (generalizedKirkwoodTerm) {
718       pmeTimings.gkEnergyTotal = -System.nanoTime();
719       generalizedKirkwood.computeInducedGKField();
720       pmeTimings.gkEnergyTotal += System.nanoTime();
721       logger.fine(
722           format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
723     }
724 
725     inducedDipoleFieldReduceRegion.init(
726         atoms,
727         inducedDipole,
728         inducedDipoleCR,
729         generalizedKirkwoodTerm,
730         generalizedKirkwood,
731         ewaldParameters,
732         soluteDielectric,
733         cartesianInducedDipolePhi,
734         cartesianInducedDipolePhiCR,
735         field,
736         fieldCR);
737     inducedDipoleFieldReduceRegion.executeWith(parallelTeam);
738   }
739 
740   public void destroy() {
741     if (fftTeam != null) {
742       try {
743         fftTeam.shutdown();
744       } catch (Exception ex) {
745         logger.warning(" Exception in shutting down fftTeam");
746       }
747     }
748     if (sectionTeam != null) {
749       try {
750         sectionTeam.shutdown();
751       } catch (Exception ex) {
752         logger.warning(" Exception in shutting down sectionTeam");
753       }
754     }
755     if (realSpaceTeam != null) {
756       try {
757         realSpaceTeam.shutdown();
758       } catch (Exception ex) {
759         logger.warning(" Exception in shutting down realSpaceTeam");
760       }
761     }
762   }
763 
764   /**
765    * Calculate the PME electrostatic energy.
766    *
767    * @param gradient If <code>true</code>, the gradient will be calculated.
768    * @param print    If <code>true</code>, extra logging is enabled.
769    * @return return the total electrostatic energy (permanent + polarization).
770    */
771   public double energy(boolean gradient, boolean print) {
772 
773     this.gradient = gradient;
774 
775     // Initialize energy variables.
776     totalMultipoleEnergy = 0.0;
777     permanentMultipoleEnergy = 0.0;
778     permanentRealSpaceEnergy = 0.0;
779     permanentSelfEnergy = 0.0;
780     permanentReciprocalEnergy = 0.0;
781     polarizationEnergy = 0.0;
782     inducedRealSpaceEnergy = 0.0;
783     inducedSelfEnergy = 0.0;
784     inducedReciprocalEnergy = 0.0;
785     solvationEnergy = 0.0;
786 
787     // Initialize number of interactions.
788     interactions = 0;
789     gkInteractions = 0;
790 
791     // Initialize timing variables.
792     pmeTimings.init();
793     if (reciprocalSpace != null) {
794       reciprocalSpace.initTimings();
795     }
796 
797     // Initialize Lambda variables.
798     if (lambdaTerm) {
799       shareddEdLambda.set(0.0);
800       sharedd2EdLambda2.set(0.0);
801     }
802 
803     if (esvTerm) {
804       extendedSystem.initEsvPermElec();
805       extendedSystem.initEsvIndElec();
806     }
807 
808     alchemicalParameters.doPermanentRealSpace = true;
809     alchemicalParameters.permanentScale = 1.0;
810     alchemicalParameters.doPolarization = true;
811     alchemicalParameters.polarizationScale = 1.0;
812 
813     // Expand coordinates and rotate multipoles into the global frame.
814     initializationRegion.init(lambdaTerm, extendedSystem,
815         atoms, coordinates, crystal, frame, axisAtom, globalMultipole,
816         dMultipoledTirationESV, dMultipoledTautomerESV,
817         polarizability, dPolardTitrationESV, dPolardTautomerESV,
818         thole, ipdamp, use, neighborLists,
819         realSpaceNeighborParameters.realSpaceLists,
820         alchemicalParameters.vaporLists,
821         grad, torque, lambdaGrad, lambdaTorque);
822     initializationRegion.executeWith(parallelTeam);
823 
824     // Initialize GeneralizedKirkwood.
825     if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
826       generalizedKirkwood.init();
827     }
828 
829     // Total permanent + polarization energy.
830     double energy;
831     if (!lambdaTerm) {
832       lambdaMode = LambdaMode.OFF;
833       energy = computeEnergy(print);
834       if (nnTerm) {
835         lambdaMode = LambdaMode.VAPOR;
836         double temp = energy;
837         energy = nnElec();
838         logger.fine(format(" Vacuum energy: %20.8f", energy - temp));
839       }
840     } else {
841       // Condensed phase with all atoms.
842       lambdaMode = LambdaMode.CONDENSED;
843       energy = condensedEnergy();
844       if (logger.isLoggable(Level.FINE)) {
845         logger.fine(format(" Condensed energy: %20.8f", energy));
846       }
847 
848       // Condensed phase SCF without ligand atoms.
849       lambdaMode = LambdaMode.CONDENSED_NO_LIGAND;
850       double temp = energy;
851       energy = condensedNoLigandSCF();
852       if (logger.isLoggable(Level.FINE)) {
853         logger.fine(format(" Condensed no ligand energy: %20.8f", energy - temp));
854       }
855 
856       // Vapor ligand electrostatics.
857       if (alchemicalParameters.doLigandVaporElec) {
858         lambdaMode = LambdaMode.VAPOR;
859         temp = energy;
860         energy = ligandElec();
861         if (logger.isLoggable(Level.FINE)) {
862           logger.fine(format(" Vacuum energy: %20.8f", energy - temp));
863         }
864       }
865     }
866 
867     /*
868      Convert torques to gradients on multipole frame defining atoms. Add
869      to electrostatic gradient to the total XYZ gradient.
870     */
871     if (gradient || lambdaTerm) {
872       reduceRegion.init(lambdaTerm, gradient,
873           atoms, coordinates, frame, axisAtom,
874           grad, torque, lambdaGrad, lambdaTorque);
875       reduceRegion.excuteWith(parallelTeam);
876     }
877 
878     // Log some timings.
879     if (logger.isLoggable(Level.FINE)) {
880       pmeTimings.printRealSpaceTimings(maxThreads, realSpaceEnergyRegion);
881       if (ewaldParameters.aewald > 0.0 && reciprocalSpaceTerm) {
882         reciprocalSpace.printTimings();
883       }
884     }
885 
886     return permanentMultipoleEnergy + polarizationEnergy;
887   }
888 
889   public void expandInducedDipoles() {
890     if (nSymm > 1) {
891       expandInducedDipolesRegion.init(atoms, crystal, inducedDipole, inducedDipoleCR);
892       expandInducedDipolesRegion.executeWith(parallelTeam);
893     }
894   }
895 
896   public int[][] getAxisAtoms() {
897     return axisAtom;
898   }
899 
900   public double getCavitationEnergy() {
901     return generalizedKirkwood.getCavitationEnergy();
902   }
903 
904   public double[][][] getCoordinates() {
905     return coordinates;
906   }
907 
908   public double getDispersionEnergy() {
909     return generalizedKirkwood.getDispersionEnergy();
910   }
911 
912   public double getEwaldCoefficient() {
913     return ewaldParameters.aewald;
914   }
915 
916   public double getEwaldCutoff() {
917     return ewaldParameters.off;
918   }
919 
920   public GeneralizedKirkwood getGK() {
921     return generalizedKirkwood;
922   }
923 
924   /**
925    * getGeneralizedKirkwoodEnergy.
926    *
927    * @return a double.
928    */
929   public double getGKEnergy() {
930     return generalizedKirkwood.getGeneralizedKirkwoordEnergy();
931   }
932 
933   /**
934    * getGKInteractions
935    *
936    * @return The number of GK interactions.
937    */
938   public int getGKInteractions() {
939     return gkInteractions;
940   }
941 
942   /**
943    * Getter for the field <code>interactions</code>.
944    *
945    * @return The number of PME interactions.
946    */
947   public int getInteractions() {
948     return interactions;
949   }
950 
951   /**
952    * {@inheritDoc}
953    *
954    * <p>Get the current lambda scale value.
955    */
956   @Override
957   public double getLambda() {
958     return lambda;
959   }
960 
961   /**
962    * {@inheritDoc}
963    *
964    * <p>Set the electrostatic lambda scaling factor.
965    */
966   @Override
967   public void setLambda(double lambda) {
968     assert (lambda >= 0.0 && lambda <= 1.0);
969     if (!lambdaTerm) {
970       return;
971     }
972     this.lambda = lambda;
973 
974     initSoftCore();
975     alchemicalParameters.update(lambda);
976     if (generalizedKirkwoodTerm) {
977       generalizedKirkwood.setLambda(alchemicalParameters.polLambda);
978       generalizedKirkwood.setLambdaFunction(
979           alchemicalParameters.lPowPol,
980           alchemicalParameters.dlPowPol,
981           alchemicalParameters.d2lPowPol);
982     }
983   }
984 
985   public double getPermanentEnergy() {
986     return permanentMultipoleEnergy;
987   }
988 
989   /**
990    * getIndRealEnergy.
991    *
992    * @return a double.
993    */
994   public double getIndRealEnergy() {
995     return inducedRealSpaceEnergy;
996   }
997 
998   /**
999    * getIndRecipEnergy.
1000    *
1001    * @return a double.
1002    */
1003   public double getIndRecipEnergy() {
1004     return inducedReciprocalEnergy;
1005   }
1006 
1007   /**
1008    * getIndSelfEnergy.
1009    *
1010    * @return a double.
1011    */
1012   public double getIndSelfEnergy() {
1013     return inducedSelfEnergy;
1014   }
1015 
1016   /**
1017    * getPermSelfEnergy.
1018    *
1019    * @return a double.
1020    */
1021   public double getPermSelfEnergy() {
1022     return permanentSelfEnergy;
1023   }
1024 
1025   public double getPermRealEnergy() {
1026     return permanentRealSpaceEnergy;
1027   }
1028 
1029   public double getPermRecipEnergy() {
1030     return permanentReciprocalEnergy;
1031   }
1032 
1033   public double getPolarEps() {
1034     return poleps;
1035   }
1036 
1037   public int[][] getPolarization11() {
1038     return ip11;
1039   }
1040 
1041   public int[][] getPolarization12() {
1042     return ip12;
1043   }
1044 
1045   public int[][] getPolarization13() {
1046     return ip13;
1047   }
1048 
1049   /**
1050    * Getter for the field <code>polarizationEnergy</code>.
1051    *
1052    * @return a double.
1053    */
1054   public double getPolarizationEnergy() {
1055     return polarizationEnergy;
1056   }
1057 
1058   public Polarization getPolarizationType() {
1059     return polarization;
1060   }
1061 
1062   public ReciprocalSpace getReciprocalSpace() {
1063     return reciprocalSpace;
1064   }
1065 
1066   public double getScale14() {
1067     return scaleParameters.m14scale;
1068   }
1069 
1070   /**
1071    * getGKEnergy
1072    *
1073    * @return a double.
1074    */
1075   public double getSolvationEnergy() {
1076     return solvationEnergy;
1077   }
1078 
1079   /**
1080    * Get the MultipoleType for Atom i.
1081    *
1082    * @param i The atom index.
1083    * @return The MultipoleType.
1084    */
1085   public MultipoleType getMultipoleType(int i) {
1086     Atom atom = atoms[i];
1087     MultipoleType multipoleType = atom.getMultipoleType();
1088     if (esvTerm && extendedSystem.isTitrating(i)) {
1089       double[] multipole = multipoleType.getMultipole();
1090       double[] esvMultipole = new double[10];
1091       arraycopy(multipole, 0, esvMultipole, 0, multipole.length);
1092       double titrationLambda = extendedSystem.getTitrationLambda(i);
1093       double tautomerLambda = extendedSystem.getTautomerLambda(i);
1094       esvMultipole = extendedSystem.getTitrationUtils()
1095           .getMultipole(atom, titrationLambda, tautomerLambda, esvMultipole);
1096       // Create a new MultipoleType for the tritrating site.
1097       multipoleType = new MultipoleType(multipoleType, esvMultipole);
1098     }
1099     return multipoleType;
1100   }
1101 
1102   /**
1103    * Get the PolarizeType for Atom i.
1104    *
1105    * @param i The atom index.
1106    * @return The PolarizeType.
1107    */
1108   public PolarizeType getPolarizeType(int i) {
1109     Atom atom = atoms[i];
1110     PolarizeType polarizeType = atom.getPolarizeType();
1111     if (polarizeType != null) {
1112       if (esvTerm && extendedSystem.isTitrating(i) && (extendedSystem.isTitratingHydrogen(i)
1113           || extendedSystem.isTitratingHeavy(i))) {
1114         double titrationLambda = extendedSystem.getTitrationLambda(i);
1115         double tautomerLambda = extendedSystem.getTautomerLambda(i);
1116         double esvPolarizability = extendedSystem.getTitrationUtils()
1117             .getPolarizability(atom, titrationLambda, tautomerLambda, polarizeType.polarizability);
1118         polarizeType = new PolarizeType(polarizeType, esvPolarizability);
1119       }
1120     }
1121     return polarizeType;
1122   }
1123 
1124   /**
1125    * {@inheritDoc}
1126    */
1127   @Override
1128   public double getd2EdL2() {
1129     if (sharedd2EdLambda2 == null || !lambdaTerm) {
1130       return 0.0;
1131     }
1132     double d2EdL2 = sharedd2EdLambda2.get();
1133     if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
1134       d2EdL2 += generalizedKirkwood.getd2EdL2();
1135     }
1136     return d2EdL2;
1137   }
1138 
1139   /**
1140    * {@inheritDoc}
1141    */
1142   @Override
1143   public double getdEdL() {
1144     if (shareddEdLambda == null || !lambdaTerm) {
1145       return 0.0;
1146     }
1147     double dEdL = shareddEdLambda.get();
1148     if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
1149       dEdL += generalizedKirkwood.getdEdL();
1150     }
1151     return dEdL;
1152   }
1153 
1154   /**
1155    * {@inheritDoc}
1156    */
1157   @Override
1158   public void getdEdXdL(double[] gradient) {
1159     if (lambdaGrad == null || !lambdaTerm) {
1160       return;
1161     }
1162     // Note that the Generalized Kirkwood contributions are already in the lambdaGrad array.
1163     int index = 0;
1164     for (int i = 0; i < nAtoms; i++) {
1165       if (atoms[i].isActive()) {
1166         gradient[index++] += lambdaGrad.getX(i);
1167         gradient[index++] += lambdaGrad.getY(i);
1168         gradient[index++] += lambdaGrad.getZ(i);
1169       }
1170     }
1171   }
1172 
1173   public void setAtoms(Atom[] atoms, int[] molecule) {
1174     if (lambdaTerm && atoms.length != nAtoms) {
1175       logger.severe(" Changing the number of atoms is not compatible with use of Lambda.");
1176     }
1177     this.atoms = atoms;
1178     this.molecule = molecule;
1179     nAtoms = atoms.length;
1180     initAtomArrays();
1181 
1182     if (reciprocalSpace != null) {
1183       reciprocalSpace.setAtoms(atoms);
1184     }
1185 
1186     if (generalizedKirkwood != null) {
1187       generalizedKirkwood.setAtoms(atoms);
1188     }
1189   }
1190 
1191   public void setCrystal(Crystal crystal) {
1192     // Check if memory allocation is required.
1193     int nSymmNew = crystal.spaceGroup.getNumberOfSymOps();
1194     if (nSymm < nSymmNew) {
1195       coordinates = new double[nSymmNew][3][nAtoms];
1196       globalMultipole = new double[nSymmNew][nAtoms][10];
1197       fractionalMultipole = new double[nSymmNew][nAtoms][10];
1198       inducedDipole = new double[nSymmNew][nAtoms][3];
1199       inducedDipoleCR = new double[nSymmNew][nAtoms][3];
1200       if (generalizedKirkwood != null) {
1201         generalizedKirkwood.setAtoms(atoms);
1202       }
1203       realSpaceNeighborParameters.allocate(nAtoms, nSymmNew);
1204       pcgSolver.allocateLists(nSymmNew, nAtoms);
1205     }
1206     nSymm = nSymmNew;
1207     neighborLists = neighborList.getNeighborList();
1208     this.crystal = crystal;
1209     if (reciprocalSpace != null) {
1210       reciprocalSpace.setCrystal(crystal.getUnitCell());
1211     }
1212   }
1213 
1214   private void initAtomArrays() {
1215     if (localMultipole == null || localMultipole.length < nAtoms) {
1216       localMultipole = new double[nAtoms][10];
1217       frame = new MultipoleFrameDefinition[nAtoms];
1218       axisAtom = new int[nAtoms][];
1219       cartesianMultipolePhi = new double[nAtoms][tensorCount];
1220       fracMultipolePhi = new double[nAtoms][tensorCount];
1221       directDipole = new double[nAtoms][3];
1222       directDipoleCR = new double[nAtoms][3];
1223       directField = new double[nAtoms][3];
1224       directFieldCR = new double[nAtoms][3];
1225       vacuumDirectDipole = new double[nAtoms][3];
1226       vacuumDirectDipoleCR = new double[nAtoms][3];
1227       if (optRegion != null) {
1228         int optOrder = optRegion.optOrder;
1229         optRegion.optDipole = new double[optOrder + 1][nAtoms][3];
1230         optRegion.optDipoleCR = new double[optOrder + 1][nAtoms][3];
1231       }
1232       // Total Induced Dipole Phi (i.e. with GK)
1233       cartesianInducedDipolePhi = new double[nAtoms][tensorCount];
1234       cartesianInducedDipolePhiCR = new double[nAtoms][tensorCount];
1235       fractionalInducedDipolePhi = new double[nAtoms][tensorCount];
1236       fractionalInducedDipolePhiCR = new double[nAtoms][tensorCount];
1237       // Vacuum Induced Dipole Phi (i.e. no GK)
1238       cartesianVacuumDipolePhi = new double[nAtoms][tensorCount];
1239       cartesianVacuumDipolePhiCR = new double[nAtoms][tensorCount];
1240       fractionalVacuumDipolePhi = new double[nAtoms][tensorCount];
1241       fractionalVacuumDipolePhiCR = new double[nAtoms][tensorCount];
1242 
1243       mask12 = new int[nAtoms][];
1244       mask13 = new int[nAtoms][];
1245       mask14 = new int[nAtoms][];
1246       mask15 = new int[nAtoms][];
1247       ip11 = new int[nAtoms][];
1248       ip12 = new int[nAtoms][];
1249       ip13 = new int[nAtoms][];
1250       thole = new double[nAtoms];
1251       ipdamp = new double[nAtoms];
1252       polarizability = new double[nAtoms];
1253 
1254       if (scfAlgorithm == SCFAlgorithm.CG) {
1255         pcgSolver.allocateVectors(nAtoms);
1256       }
1257       pcgSolver.allocateLists(nSymm, nAtoms);
1258 
1259       if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
1260         int predictorOrder = scfPredictorParameters.predictorOrder;
1261         if (lambdaTerm) {
1262           scfPredictorParameters.predictorInducedDipole = new double[3][predictorOrder][nAtoms][3];
1263           scfPredictorParameters.predictorInducedDipoleCR =
1264               new double[3][predictorOrder][nAtoms][3];
1265         } else {
1266           scfPredictorParameters.predictorInducedDipole = new double[1][predictorOrder][nAtoms][3];
1267           scfPredictorParameters.predictorInducedDipoleCR =
1268               new double[1][predictorOrder][nAtoms][3];
1269         }
1270       }
1271 
1272       // Initialize per-thread memory for collecting the gradient, torque, field and chain-rule
1273       // field.
1274       grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1275       torque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1276       field = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1277       fieldCR = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1278       lambdaGrad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1279       lambdaTorque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1280       isSoft = new boolean[nAtoms];
1281       use = new boolean[nAtoms];
1282 
1283       coordinates = new double[nSymm][3][nAtoms];
1284       globalMultipole = new double[nSymm][nAtoms][10];
1285       fractionalMultipole = new double[nSymm][nAtoms][10];
1286       inducedDipole = new double[nSymm][nAtoms][3];
1287       inducedDipoleCR = new double[nSymm][nAtoms][3];
1288       vacuumInducedDipole = new double[nSymm][nAtoms][3];
1289       vacuumInducedDipoleCR = new double[nSymm][nAtoms][3];
1290 
1291       // The size of reduced neighbor list depends on the size of the real space cutoff.
1292       realSpaceNeighborParameters.allocate(nAtoms, nSymm);
1293 
1294       // Lambda factors are different for OST and ESV interactions.
1295       lambdaFactors = new LambdaFactors[maxThreads];
1296       for (int i = 0; i < maxThreads; i++) {
1297         /*if (esvTerm) {
1298           // Invoked every time through inner loops.
1299           lambdaFactors[i] = new LambdaFactorsESV();
1300         } else */
1301         if (lambdaTerm) {
1302           // Invoked on calls to setLambda().
1303           lambdaFactors[i] = new LambdaFactorsOST();
1304         } else {
1305           // Invoked never; inoperative defaults.
1306           lambdaFactors[i] = LambdaDefaults;
1307         }
1308       }
1309     }
1310 
1311     // Initialize the soft core lambda mask to false for all atoms.
1312     fill(isSoft, false);
1313 
1314     // Initialize the use mask to true for all atoms.
1315     fill(use, true);
1316 
1317     // Assign multipole parameters.
1318     assignMultipoles();
1319 
1320     // Assign polarization groups.
1321     if (elecForm == PAM) {
1322       assignPolarizationGroups();
1323     }
1324 
1325     // Fill the thole, inverse polarization damping and polarizability arrays.
1326     for (int i = 0; i < nAtoms; i++) {
1327       Atom ai = atoms[i];
1328 
1329       if (elecForm == PAM) {
1330         PolarizeType polarizeType = ai.getPolarizeType();
1331         int index = ai.getIndex() - 1;
1332         thole[index] = polarizeType.thole;
1333         ipdamp[index] = polarizeType.pdamp;
1334         if (!(ipdamp[index] > 0.0)) {
1335           ipdamp[index] = Double.POSITIVE_INFINITY;
1336         } else {
1337           ipdamp[index] = 1.0 / ipdamp[index];
1338         }
1339         polarizability[index] = polarizeType.polarizability;
1340       }
1341 
1342       // Collect 1-2 interactions.
1343       List<Atom> n12 = ai.get12List();
1344       mask12[i] = new int[n12.size()];
1345       int j = 0;
1346       for (Atom a12 : n12) {
1347         mask12[i][j++] = a12.getIndex() - 1;
1348       }
1349 
1350       // Collect 1-3 interactions.
1351       List<Atom> n13 = ai.get13List();
1352       mask13[i] = new int[n13.size()];
1353       j = 0;
1354       for (Atom a13 : n13) {
1355         mask13[i][j++] = a13.getIndex() - 1;
1356       }
1357 
1358       // Collect 1-4 interactions.
1359       List<Atom> n14 = ai.get14List();
1360       mask14[i] = new int[n14.size()];
1361       j = 0;
1362       for (Atom a14 : n14) {
1363         mask14[i][j++] = a14.getIndex() - 1;
1364       }
1365 
1366       // Collect 1-5 interactions.
1367       List<Atom> n15 = ai.get15List();
1368       mask15[i] = new int[n15.size()];
1369       j = 0;
1370       for (Atom a15 : n15) {
1371         mask15[i][j++] = a15.getIndex() - 1;
1372       }
1373     }
1374   }
1375 
1376   /**
1377    * Initialize a boolean array of soft atoms and, if requested, ligand vapor electrostatics.
1378    */
1379   private void initSoftCore() {
1380 
1381     boolean rebuild = false;
1382 
1383     // Initialize a boolean array that marks soft atoms.
1384     StringBuilder sb = new StringBuilder("\n Softcore Atoms:\n");
1385     int count = 0;
1386     for (int i = 0; i < nAtoms; i++) {
1387       Atom ai = atoms[i];
1388       boolean soft = ai.applyLambda();
1389       if (soft != isSoft[i]) {
1390         rebuild = true;
1391       }
1392       isSoft[i] = soft;
1393       if (soft) {
1394         count++;
1395         sb.append(ai).append("\n");
1396       }
1397     }
1398 
1399     // Force rebuild ligand vapor electrostatics are being computed and vaporCrystal is null.
1400     if (alchemicalParameters.doLigandVaporElec
1401         && alchemicalParameters.vaporCrystal == null) {
1402       rebuild = true;
1403     }
1404 
1405     if (!rebuild) {
1406       return;
1407     }
1408 
1409     if (count > 0 && logger.isLoggable(Level.FINE)) {
1410       logger.fine(format(" Softcore atom count: %d", count));
1411       logger.fine(sb.toString());
1412     }
1413 
1414     // Initialize boundary conditions,
1415     // the n^2 neighbor list and parallel scheduling for ligand vapor electrostatics.
1416     if (alchemicalParameters.doLigandVaporElec) {
1417       double vacuumOff = getVacuumOff();
1418       alchemicalParameters.vaporCrystal =
1419           new Crystal(3 * vacuumOff, 3 * vacuumOff, 3 * vacuumOff, 90.0, 90.0, 90.0, "P1");
1420       alchemicalParameters.vaporCrystal.setAperiodic(true);
1421       NeighborList vacuumNeighborList = new NeighborList(null,
1422           alchemicalParameters.vaporCrystal, atoms, vacuumOff, 2.0, parallelTeam);
1423       vacuumNeighborList.setIntermolecular(false);
1424 
1425       alchemicalParameters.vaporLists = new int[1][nAtoms][];
1426       double[][] coords = new double[1][nAtoms * 3];
1427       for (int i = 0; i < nAtoms; i++) {
1428         coords[0][i * 3] = atoms[i].getX();
1429         coords[0][i * 3 + 1] = atoms[i].getY();
1430         coords[0][i * 3 + 2] = atoms[i].getZ();
1431       }
1432       boolean print = logger.isLoggable(Level.FINE);
1433       vacuumNeighborList.buildList(coords, alchemicalParameters.vaporLists, isSoft,
1434           true, print);
1435       alchemicalParameters.vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
1436       alchemicalParameters.vaporEwaldSchedule = alchemicalParameters.vaporPermanentSchedule;
1437       alchemicalParameters.vacuumRanges = new Range[maxThreads];
1438       vacuumNeighborList.setDisableUpdates(
1439           forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
1440     } else {
1441       alchemicalParameters.vaporCrystal = null;
1442       alchemicalParameters.vaporLists = null;
1443       alchemicalParameters.vaporPermanentSchedule = null;
1444       alchemicalParameters.vaporEwaldSchedule = null;
1445       alchemicalParameters.vacuumRanges = null;
1446     }
1447   }
1448 
1449   private double getVacuumOff() {
1450     double maxr = 10.0;
1451     for (int i = 0; i < nAtoms; i++) {
1452       Atom ai = atoms[i];
1453       if (ai.applyLambda()) {
1454 
1455         // Determine ligand size.
1456         for (int j = i + 1; j < nAtoms; j++) {
1457           Atom aj = atoms[j];
1458           if (aj.applyLambda()) {
1459             double dx = ai.getX() - aj.getX();
1460             double dy = ai.getY() - aj.getY();
1461             double dz = ai.getZ() - aj.getZ();
1462             double r = sqrt(dx * dx + dy * dy + dz * dz);
1463             maxr = max(r, maxr);
1464           }
1465         }
1466       }
1467     }
1468 
1469     double vacuumOff = 2 * maxr;
1470     return vacuumOff;
1471   }
1472 
1473   /**
1474    * Initialize a boolean array of soft atoms and, if requested, ligand vapor electrostatics.
1475    */
1476   private void initNN() {
1477 
1478     // Initialize a boolean array that marks soft atoms.
1479     StringBuilder sb = new StringBuilder("\n NN Atoms:\n");
1480     boolean[] nnAtoms = new boolean[nAtoms];
1481 
1482     int count = 0;
1483     for (int i = 0; i < nAtoms; i++) {
1484       Atom ai = atoms[i];
1485       nnAtoms[i] = ai.isNeuralNetwork();
1486       if (nnAtoms[i]) {
1487         count++;
1488         sb.append(ai).append("\n");
1489       }
1490     }
1491 
1492     if (count > 0 && logger.isLoggable(Level.FINE)) {
1493       logger.fine(format(" Neural network atom count: %d", count));
1494       logger.fine(sb.toString());
1495     }
1496 
1497     // Initialize boundary conditions,
1498     // An n^2 neighbor list and parallel scheduling for ligand vapor electrostatics.
1499     if (alchemicalParameters.doLigandVaporElec) {
1500       alchemicalParameters.vaporCrystal = new Crystal(10.0, 10.0, 10.0,
1501           90.0, 90.0, 90.0, "P1");
1502       alchemicalParameters.vaporCrystal.setAperiodic(true);
1503       NeighborList vacuumNeighborList = new NeighborList(null,
1504           alchemicalParameters.vaporCrystal, atoms, Double.POSITIVE_INFINITY, 2.0, parallelTeam);
1505       vacuumNeighborList.setIntermolecular(false);
1506 
1507       alchemicalParameters.vaporLists = new int[1][nAtoms][];
1508       double[][] coords = new double[1][nAtoms * 3];
1509       for (int i = 0; i < nAtoms; i++) {
1510         coords[0][i * 3] = atoms[i].getX();
1511         coords[0][i * 3 + 1] = atoms[i].getY();
1512         coords[0][i * 3 + 2] = atoms[i].getZ();
1513       }
1514       boolean print = logger.isLoggable(Level.FINE);
1515       vacuumNeighborList.buildList(coords, alchemicalParameters.vaporLists, nnAtoms, true, print);
1516       alchemicalParameters.vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
1517       alchemicalParameters.vaporEwaldSchedule = alchemicalParameters.vaporPermanentSchedule;
1518       alchemicalParameters.vacuumRanges = new Range[maxThreads];
1519       vacuumNeighborList.setDisableUpdates(forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
1520     } else {
1521       alchemicalParameters.vaporCrystal = null;
1522       alchemicalParameters.vaporLists = null;
1523       alchemicalParameters.vaporPermanentSchedule = null;
1524       alchemicalParameters.vaporEwaldSchedule = null;
1525       alchemicalParameters.vacuumRanges = null;
1526     }
1527   }
1528 
1529   /**
1530    * Total system under PBC.
1531    * A. Softcore real space for Ligand-Protein and Ligand-Ligand.
1532    * B. Reciprocal space scaled by lambda.
1533    * C. Polarization scaled by lambda.
1534    */
1535   private double condensedEnergy() {
1536     if (lambda < alchemicalParameters.polLambdaStart) {
1537       /*
1538        If the polarization has been completely decoupled, the
1539        contribution of the complete system is zero.
1540        We can skip the SCF for part 1 for efficiency.
1541       */
1542       alchemicalParameters.polarizationScale = 0.0;
1543       alchemicalParameters.doPolarization = false;
1544     } else if (lambda <= alchemicalParameters.polLambdaEnd) {
1545       alchemicalParameters.polarizationScale = alchemicalParameters.lPowPol;
1546       alchemicalParameters.doPolarization = true;
1547     } else {
1548       alchemicalParameters.polarizationScale = 1.0;
1549       alchemicalParameters.doPolarization = true;
1550     }
1551     alchemicalParameters.doPermanentRealSpace = true;
1552     alchemicalParameters.permanentScale = alchemicalParameters.lPowPerm;
1553     alchemicalParameters.dEdLSign = 1.0;
1554 
1555     return computeEnergy(false);
1556   }
1557 
1558   /**
1559    * Condensed phase system without the ligand.
1560    * A. No permanent real space electrostatics needs to be calculated because this was handled analytically in step 1.
1561    * B. Permanent reciprocal space scaled by (1 - lambda).
1562    * C. Polarization scaled by (1 - lambda).
1563    */
1564   private double condensedNoLigandSCF() {
1565     // Turn off the ligand.
1566     boolean skip = true;
1567     for (int i = 0; i < nAtoms; i++) {
1568       if (atoms[i].applyLambda()) {
1569         use[i] = false;
1570       } else {
1571         use[i] = true;
1572         skip = false;
1573       }
1574     }
1575 
1576     // Permanent real space is done for the condensed phase. Scale the reciprocal space part.
1577     alchemicalParameters.doPermanentRealSpace = false;
1578     alchemicalParameters.permanentScale =
1579         1.0 - alchemicalParameters.lPowPerm;
1580     alchemicalParameters.dEdLSign = -1.0;
1581 
1582     // If we are past the end of the polarization lambda window, then only the condensed phase is
1583     // necessary.
1584     if (lambda <= alchemicalParameters.polLambdaEnd
1585         && alchemicalParameters.doNoLigandCondensedSCF) {
1586       alchemicalParameters.doPolarization = true;
1587       alchemicalParameters.polarizationScale =
1588           1.0 - alchemicalParameters.lPowPol;
1589     } else {
1590       alchemicalParameters.doPolarization = false;
1591       alchemicalParameters.polarizationScale = 0.0;
1592     }
1593 
1594     // Turn off GK.
1595     boolean gkBack = generalizedKirkwoodTerm;
1596     generalizedKirkwoodTerm = false;
1597 
1598     // If we are disappearing the entire system (i.e. a small crystal) then
1599     // the energy of this step is 0, and we can skip it.
1600     double energy;
1601     if (skip) {
1602       energy = permanentMultipoleEnergy + polarizationEnergy + solvationEnergy;
1603     } else {
1604       energy = computeEnergy(false);
1605       Arrays.fill(use, true);
1606     }
1607 
1608     generalizedKirkwoodTerm = gkBack;
1609 
1610     return energy;
1611   }
1612 
1613   /**
1614    * Aperiodic ligand electrostatics.
1615    * A. Real space with an Ewald coefficient of 0.0 (no reciprocal space).
1616    * B. Polarization scaled as in Step 2 by (1 - lambda).
1617    */
1618   private double ligandElec() {
1619     for (int i = 0; i < nAtoms; i++) {
1620       use[i] = atoms[i].applyLambda();
1621     }
1622 
1623     // Scale the permanent vacuum electrostatics. The softcore alpha is not
1624     // necessary (nothing in vacuum to collide with).
1625     alchemicalParameters.doPermanentRealSpace = true;
1626     alchemicalParameters.permanentScale = 1.0 - alchemicalParameters.lPowPerm;
1627     alchemicalParameters.dEdLSign = -1.0;
1628     double lAlphaBack = alchemicalParameters.lAlpha;
1629     double dlAlphaBack = alchemicalParameters.dlAlpha;
1630     double d2lAlphaBack = alchemicalParameters.d2lAlpha;
1631     alchemicalParameters.lAlpha = 0.0;
1632     alchemicalParameters.dlAlpha = 0.0;
1633     alchemicalParameters.d2lAlpha = 0.0;
1634 
1635     // If we are past the end of the polarization lambda window, then only
1636     // the condensed phase is necessary.
1637     if (lambda <= alchemicalParameters.polLambdaEnd) {
1638       alchemicalParameters.doPolarization = true;
1639       alchemicalParameters.polarizationScale = 1.0 - alchemicalParameters.lPowPol;
1640     } else {
1641       alchemicalParameters.doPolarization = false;
1642       alchemicalParameters.polarizationScale = 0.0;
1643     }
1644 
1645     // Save the current real space PME parameters.
1646     double offBack = ewaldParameters.off;
1647     double aewaldBack = ewaldParameters.aewald;
1648     ewaldParameters.setEwaldParameters(Double.MAX_VALUE, 0.0);
1649     // Save the current parallelization schedule.
1650     IntegerSchedule permanentScheduleBack = permanentSchedule;
1651     IntegerSchedule ewaldScheduleBack = realSpaceNeighborParameters.realSpaceSchedule;
1652     Range[] rangesBack = realSpaceNeighborParameters.realSpaceRanges;
1653     permanentSchedule = alchemicalParameters.vaporPermanentSchedule;
1654     realSpaceNeighborParameters.realSpaceSchedule = alchemicalParameters.vaporEwaldSchedule;
1655     realSpaceNeighborParameters.realSpaceRanges = alchemicalParameters.vacuumRanges;
1656 
1657     // Use vacuum crystal / vacuum neighborLists.
1658     Crystal crystalBack = crystal;
1659     int nSymmBack = nSymm;
1660     int[][][] listsBack = neighborLists;
1661     neighborLists = alchemicalParameters.vaporLists;
1662     crystal = alchemicalParameters.vaporCrystal;
1663     nSymm = 1;
1664 
1665     // Turn off GK if in use.
1666     boolean gkBack = generalizedKirkwoodTerm;
1667 
1668     // Turn off Pre-conditioned conjugate gradient SCF solver.
1669     SCFAlgorithm scfBack = scfAlgorithm;
1670     scfAlgorithm = SCFAlgorithm.SOR;
1671 
1672     if (alchemicalParameters.doLigandGKElec) {
1673       generalizedKirkwoodTerm = true;
1674       generalizedKirkwood.setNeighborList(alchemicalParameters.vaporLists);
1675       generalizedKirkwood.setLambda(lambda);
1676       generalizedKirkwood.setCutoff(ewaldParameters.off);
1677       generalizedKirkwood.setCrystal(alchemicalParameters.vaporCrystal);
1678       generalizedKirkwood.setLambdaFunction(alchemicalParameters.polarizationScale,
1679           alchemicalParameters.dEdLSign * alchemicalParameters.dlPowPol,
1680           alchemicalParameters.dEdLSign * alchemicalParameters.d2lPowPol);
1681     } else {
1682       generalizedKirkwoodTerm = false;
1683     }
1684 
1685     double energy = computeEnergy(false);
1686 
1687     // Revert to the saved parameters.
1688     ewaldParameters.setEwaldParameters(offBack, aewaldBack);
1689     neighborLists = listsBack;
1690     crystal = crystalBack;
1691     nSymm = nSymmBack;
1692     permanentSchedule = permanentScheduleBack;
1693     realSpaceNeighborParameters.realSpaceSchedule = ewaldScheduleBack;
1694     realSpaceNeighborParameters.realSpaceRanges = rangesBack;
1695     alchemicalParameters.lAlpha = lAlphaBack;
1696     alchemicalParameters.dlAlpha = dlAlphaBack;
1697     alchemicalParameters.d2lAlpha = d2lAlphaBack;
1698     generalizedKirkwoodTerm = gkBack;
1699     scfAlgorithm = scfBack;
1700 
1701     fill(use, true);
1702 
1703     return energy;
1704   }
1705 
1706   /**
1707    * Aperiodic electrostatics when using neural networks.
1708    * A. Real space with an Ewald coefficient of 0.0 (no reciprocal space).
1709    * B. Permanent and Polarization scaled by -1.
1710    */
1711   private double nnElec() {
1712     for (int i = 0; i < nAtoms; i++) {
1713       use[i] = atoms[i].isNeuralNetwork();
1714     }
1715 
1716     // Scale the permanent vacuum electrostatics.
1717     // The softcore alpha is not necessary.
1718     alchemicalParameters.doPermanentRealSpace = true;
1719     alchemicalParameters.permanentScale = -1.0;
1720     alchemicalParameters.dEdLSign = 1.0;
1721     double lAlphaBack = alchemicalParameters.lAlpha;
1722     double dlAlphaBack = alchemicalParameters.dlAlpha;
1723     double d2lAlphaBack = alchemicalParameters.d2lAlpha;
1724     alchemicalParameters.lAlpha = 0.0;
1725     alchemicalParameters.dlAlpha = 0.0;
1726     alchemicalParameters.d2lAlpha = 0.0;
1727     alchemicalParameters.doPolarization = true;
1728     alchemicalParameters.polarizationScale = -1.0;
1729 
1730     // Save the current real space PME parameters.
1731     double offBack = ewaldParameters.off;
1732     double aewaldBack = ewaldParameters.aewald;
1733     ewaldParameters.setEwaldParameters(Double.POSITIVE_INFINITY, 0.0);
1734     // Save the current parallelization schedule.
1735     IntegerSchedule permanentScheduleBack = permanentSchedule;
1736     IntegerSchedule ewaldScheduleBack = realSpaceNeighborParameters.realSpaceSchedule;
1737     Range[] rangesBack = realSpaceNeighborParameters.realSpaceRanges;
1738     permanentSchedule = alchemicalParameters.vaporPermanentSchedule;
1739     realSpaceNeighborParameters.realSpaceSchedule = alchemicalParameters.vaporEwaldSchedule;
1740     realSpaceNeighborParameters.realSpaceRanges = alchemicalParameters.vacuumRanges;
1741 
1742     // Use vacuum crystal / vacuum neighborLists.
1743     Crystal crystalBack = crystal;
1744     int nSymmBack = nSymm;
1745     int[][][] listsBack = neighborLists;
1746     neighborLists = alchemicalParameters.vaporLists;
1747     crystal = alchemicalParameters.vaporCrystal;
1748     nSymm = 1;
1749 
1750     // Turn off GK if in use.
1751     boolean gkBack = generalizedKirkwoodTerm;
1752 
1753     // Turn off Pre-conditioned conjugate gradient SCF solver.
1754     SCFAlgorithm scfBack = scfAlgorithm;
1755     scfAlgorithm = SCFAlgorithm.SOR;
1756     generalizedKirkwoodTerm = false;
1757 
1758     double energy = computeEnergy(false);
1759 
1760     // Revert to the saved parameters.
1761     ewaldParameters.setEwaldParameters(offBack, aewaldBack);
1762     neighborLists = listsBack;
1763     crystal = crystalBack;
1764     nSymm = nSymmBack;
1765     permanentSchedule = permanentScheduleBack;
1766     realSpaceNeighborParameters.realSpaceSchedule = ewaldScheduleBack;
1767     realSpaceNeighborParameters.realSpaceRanges = rangesBack;
1768     alchemicalParameters.lAlpha = lAlphaBack;
1769     alchemicalParameters.dlAlpha = dlAlphaBack;
1770     alchemicalParameters.d2lAlpha = d2lAlphaBack;
1771     generalizedKirkwoodTerm = gkBack;
1772     scfAlgorithm = scfBack;
1773 
1774     fill(use, true);
1775 
1776     return energy;
1777   }
1778 
1779   /**
1780    * Calculate the PME electrostatic energy for a Lambda state.
1781    *
1782    * @param print If <code>true</code>, extra logging is enabled.
1783    * @return return the total electrostatic energy (permanent + polarization).
1784    */
1785   private double computeEnergy(boolean print) {
1786     // Find the permanent multipole potential, field, etc.
1787     permanentMultipoleField();
1788 
1789     // Compute Born radii if necessary.
1790     if (generalizedKirkwoodTerm) {
1791       pmeTimings.bornRadiiTotal -= System.nanoTime();
1792       generalizedKirkwood.setUse(use);
1793       generalizedKirkwood.computeBornRadii();
1794       pmeTimings.bornRadiiTotal += System.nanoTime();
1795     }
1796 
1797     // Do the self-consistent field calculation.
1798     if (polarization != Polarization.NONE && alchemicalParameters.doPolarization) {
1799       // Compute vacuum dipole moments.
1800       if (generalizedKirkwoodTerm) {
1801         generalizedKirkwoodTerm = false;
1802         // Run the vacuum SCF.
1803         selfConsistentField(logger.isLoggable(Level.FINE));
1804         // Store vacuum dipole moments
1805         saveInducedDipolesToVacuumDipoles();
1806 
1807         generalizedKirkwoodTerm = true;
1808         // Load the permanent multipole field.
1809         permanentMultipoleField();
1810       }
1811 
1812       // Compute induced dipoles.
1813       selfConsistentField(logger.isLoggable(Level.FINE));
1814 
1815       if (esvTerm && polarization != Polarization.NONE) {
1816         for (int i = 0; i < nAtoms; i++) {
1817           if (extendedSystem.isTitrating(i)
1818               && (extendedSystem.isTitratingHydrogen(i) || extendedSystem.isTitratingHeavy(i))) {
1819             double dx = field.getX(i);
1820             double dy = field.getY(i);
1821             double dz = field.getZ(i);
1822             double dxCR = fieldCR.getX(i);
1823             double dyCR = fieldCR.getY(i);
1824             double dzCR = fieldCR.getZ(i);
1825             // Add back permanent multipole field to total field for extended system derivatives
1826             // if mutual polarization is used
1827             if (polarization == Polarization.MUTUAL) {
1828               dx += directField[i][0];
1829               dy += directField[i][1];
1830               dz += directField[i][2];
1831               dxCR += directFieldCR[i][0];
1832               dyCR += directFieldCR[i][1];
1833               dzCR += directFieldCR[i][2];
1834             }
1835             double fix = dx * dPolardTitrationESV[i] * dxCR;
1836             double fiy = dy * dPolardTitrationESV[i] * dyCR;
1837             double fiz = dz * dPolardTitrationESV[i] * dzCR;
1838             double titrdUdL = fix + fiy + fiz;
1839             double tautdUdL = 0.0;
1840             if (extendedSystem.isTautomerizing(i)) {
1841               fix = dx * dPolardTautomerESV[i] * dxCR;
1842               fiy = dy * dPolardTautomerESV[i] * dyCR;
1843               fiz = dz * dPolardTautomerESV[i] * dzCR;
1844               tautdUdL = fix + fiy + fiz;
1845             }
1846             extendedSystem.addIndElecDeriv(i, titrdUdL * electric * -0.5, tautdUdL * electric * -0.5);
1847           }
1848         }
1849       }
1850 
1851       if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1852         if (gradient && polarization == Polarization.DIRECT) {
1853           reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
1854           field.reset(parallelTeam);
1855           fieldCR.reset(parallelTeam);
1856           inducedDipoleFieldRegion.init(
1857               atoms, crystal, use, molecule,
1858               ipdamp, thole, coordinates, realSpaceNeighborParameters,
1859               inducedDipole, inducedDipoleCR,
1860               reciprocalSpaceTerm, reciprocalSpace,
1861               lambdaMode, ewaldParameters,
1862               field, fieldCR, pmeTimings);
1863           inducedDipoleFieldRegion.executeWith(sectionTeam);
1864           reciprocalSpace.computeInducedPhi(
1865               cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
1866               fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
1867         }
1868       }
1869 
1870       if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
1871         scfPredictorParameters.saveMutualInducedDipoles(lambdaMode,
1872             inducedDipole, inducedDipoleCR, directDipole, directDipoleCR);
1873       }
1874     }
1875 
1876     /*
1877      Find the total real space energy. This includes:
1878       1) the permanent multipoles in their own real space potential
1879       2) the permenant multipoles in their own reciprocal space potential
1880       3) the permanent multipoles interacting with the induced dipole real space potential
1881       4) the permanent multipoles interacting with the induced dipole reciprocal space potential
1882     */
1883 
1884     // With GK, we need to compute the polarization energy with vacuum induced dipoles.
1885     double[][][] inputDipole = inducedDipole;
1886     double[][][] inputDipoleCR = inducedDipoleCR;
1887     double[][] inputPhi = cartesianInducedDipolePhi;
1888     double[][] inputPhiCR = cartesianInducedDipolePhiCR;
1889     double[][] fracInputPhi = fractionalInducedDipolePhi;
1890     double[][] fracInputPhiCR = fractionalInducedDipolePhiCR;
1891     if (generalizedKirkwoodTerm) {
1892       inputDipole = vacuumInducedDipole;
1893       inputDipoleCR = vacuumInducedDipoleCR;
1894       inputPhi = cartesianVacuumDipolePhi;
1895       inputPhiCR = cartesianVacuumDipolePhiCR;
1896       fracInputPhi = fractionalVacuumDipolePhi;
1897       fracInputPhiCR = fractionalVacuumDipolePhiCR;
1898     }
1899 
1900     double eself = 0.0;
1901     double erecip = 0.0;
1902     double eselfi = 0.0;
1903     double erecipi = 0.0;
1904     if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1905       reciprocalEnergyRegion.init(atoms, crystal, gradient, lambdaTerm, esvTerm, use,
1906           globalMultipole, fractionalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1907           cartesianMultipolePhi, fracMultipolePhi,
1908           polarization, inputDipole, inputDipoleCR, inputPhi, inputPhiCR, fracInputPhi,
1909           fracInputPhiCR, reciprocalSpace, alchemicalParameters, extendedSystem,
1910           // Output
1911           grad, torque, lambdaGrad, lambdaTorque, shareddEdLambda, sharedd2EdLambda2);
1912       reciprocalEnergyRegion.executeWith(parallelTeam);
1913       eself = reciprocalEnergyRegion.getPermanentSelfEnergy();
1914       erecip = reciprocalEnergyRegion.getPermanentReciprocalEnergy();
1915       eselfi = reciprocalEnergyRegion.getInducedDipoleSelfEnergy();
1916       erecipi = reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
1917       interactions += nAtoms;
1918     }
1919 
1920     pmeTimings.realSpaceEnergyTotal -= System.nanoTime();
1921     realSpaceEnergyRegion.init(atoms, crystal, extendedSystem, esvTerm, coordinates, frame, axisAtom,
1922         globalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1923         inputDipole, inputDipoleCR, use, molecule,
1924         ip11, mask12, mask13, mask14, mask15, isSoft, ipdamp, thole, realSpaceNeighborParameters,
1925         gradient, lambdaTerm, nnTerm, lambdaMode, polarization,
1926         ewaldParameters, scaleParameters, alchemicalParameters,
1927         pmeTimings.realSpaceEnergyTime,
1928         // Output
1929         grad, torque, lambdaGrad, lambdaTorque, shareddEdLambda, sharedd2EdLambda2);
1930     realSpaceEnergyRegion.executeWith(parallelTeam);
1931     double ereal = realSpaceEnergyRegion.getPermanentEnergy();
1932     double ereali = realSpaceEnergyRegion.getPolarizationEnergy();
1933     interactions += realSpaceEnergyRegion.getInteractions();
1934     pmeTimings.realSpaceEnergyTotal += System.nanoTime();
1935 
1936     if (generalizedKirkwoodTerm) {
1937       // Compute the polarization energy cost to polarize the induced dipoles
1938       // from vacuum to the SCRF values.
1939       double eGK = 0.0;
1940 
1941       // Create default alchemical parameters.
1942       AlchemicalParameters alchemicalParametersGK = new AlchemicalParameters(
1943           forceField, false, false, polarization);
1944       // No permanent multipole contribution.
1945       alchemicalParametersGK.permanentScale = 0.0;
1946       alchemicalParametersGK.doPermanentRealSpace = false;
1947       // Flip the sign on the vacuum polarization energy and derivatives.
1948       alchemicalParametersGK.polarizationScale = -1.0;
1949 
1950       // Store the derivative and torque contributions in the GK arrays.
1951       AtomicDoubleArray3D gradGK = generalizedKirkwood.getGrad();
1952       AtomicDoubleArray3D torqueGK = generalizedKirkwood.getTorque();
1953 
1954       if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1955         reciprocalEnergyRegion.init(
1956             atoms, crystal, gradient, false, esvTerm,
1957             use, globalMultipole, fractionalMultipole,
1958             dMultipoledTirationESV, dMultipoledTautomerESV,
1959             cartesianMultipolePhi, fracMultipolePhi,
1960             polarization,
1961             vacuumInducedDipole, vacuumInducedDipoleCR,
1962             cartesianVacuumDipolePhi, cartesianVacuumDipolePhiCR,
1963             fractionalVacuumDipolePhi, fractionalVacuumDipolePhiCR,
1964             reciprocalSpace, alchemicalParametersGK,
1965             extendedSystem,
1966             gradGK, torqueGK, null, null,
1967             shareddEdLambda, sharedd2EdLambda2);
1968         reciprocalEnergyRegion.executeWith(parallelTeam);
1969         eGK = reciprocalEnergyRegion.getInducedDipoleSelfEnergy()
1970             + reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
1971       }
1972 
1973       pmeTimings.gkEnergyTotal -= System.nanoTime();
1974       realSpaceEnergyRegion.init(atoms, crystal,
1975           extendedSystem, esvTerm,
1976           coordinates, frame, axisAtom,
1977           globalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1978           vacuumInducedDipole, vacuumInducedDipoleCR,
1979           use, molecule,
1980           ip11, mask12, mask13, mask14, mask15, isSoft, ipdamp, thole,
1981           realSpaceNeighborParameters,
1982           gradient, false, nnTerm, lambdaMode,
1983           polarization, ewaldParameters, scaleParameters, alchemicalParametersGK,
1984           pmeTimings.realSpaceEnergyTime,
1985           // Output
1986           gradGK,
1987           torqueGK,
1988           null,
1989           null,
1990           shareddEdLambda,
1991           sharedd2EdLambda2);
1992       realSpaceEnergyRegion.executeWith(parallelTeam);
1993       eGK += realSpaceEnergyRegion.getPolarizationEnergy();
1994 
1995       // Normal sign of +1 for the SCRF polarization energy and derivatives.
1996       alchemicalParametersGK.polarizationScale = 1.0;
1997 
1998       if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1999         reciprocalEnergyRegion.init(
2000             atoms, crystal, gradient, false, esvTerm, use,
2001             globalMultipole, fractionalMultipole,
2002             dMultipoledTirationESV, dMultipoledTautomerESV,
2003             cartesianMultipolePhi, fracMultipolePhi,
2004             polarization, inducedDipole, inducedDipoleCR,
2005             cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2006             fractionalInducedDipolePhi, fractionalInducedDipolePhiCR,
2007             reciprocalSpace, alchemicalParametersGK, extendedSystem,
2008             gradGK, torqueGK, null, null,
2009             shareddEdLambda, sharedd2EdLambda2);
2010         reciprocalEnergyRegion.executeWith(parallelTeam);
2011         eGK += reciprocalEnergyRegion.getInducedDipoleSelfEnergy()
2012             + reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
2013       }
2014 
2015       pmeTimings.gkEnergyTotal -= System.nanoTime();
2016       realSpaceEnergyRegion.init(
2017           atoms,
2018           crystal,
2019           extendedSystem,
2020           esvTerm,
2021           coordinates,
2022           frame,
2023           axisAtom,
2024           globalMultipole,
2025           dMultipoledTirationESV,
2026           dMultipoledTautomerESV,
2027           inducedDipole,
2028           inducedDipoleCR,
2029           use,
2030           molecule,
2031           ip11,
2032           mask12,
2033           mask13,
2034           mask14,
2035           mask15,
2036           isSoft,
2037           ipdamp,
2038           thole,
2039           realSpaceNeighborParameters,
2040           gradient,
2041           false,
2042           nnTerm,
2043           lambdaMode,
2044           polarization,
2045           ewaldParameters,
2046           scaleParameters,
2047           alchemicalParametersGK,
2048           pmeTimings.realSpaceEnergyTime,
2049           // Output
2050           gradGK,
2051           torqueGK,
2052           null,
2053           null,
2054           shareddEdLambda,
2055           sharedd2EdLambda2);
2056       realSpaceEnergyRegion.executeWith(parallelTeam);
2057       eGK += realSpaceEnergyRegion.getPolarizationEnergy();
2058 
2059       // Compute the solvation free energy.
2060       solvationEnergy += generalizedKirkwood.solvationEnergy(eGK, gradient, print);
2061       if (gradient) {
2062         // Add the GK derivative contributions into the overall derivatives.
2063         generalizedKirkwood.reduce(grad, torque, lambdaGrad, lambdaTorque);
2064       }
2065       gkInteractions += generalizedKirkwood.getInteractions();
2066       pmeTimings.gkEnergyTotal += System.nanoTime();
2067     }
2068 
2069     // Collect energy terms.
2070     permanentRealSpaceEnergy += ereal;
2071     permanentSelfEnergy += eself;
2072     permanentReciprocalEnergy += erecip;
2073     inducedRealSpaceEnergy += ereali;
2074     inducedSelfEnergy += eselfi;
2075     inducedReciprocalEnergy += erecipi;
2076     permanentMultipoleEnergy += eself + erecip + ereal;
2077     polarizationEnergy += eselfi + erecipi + ereali;
2078     totalMultipoleEnergy += ereal + eself + erecip + ereali + eselfi + erecipi;
2079 
2080     // Log some info.
2081     if (logger.isLoggable(Level.FINE)) {
2082       StringBuilder sb = new StringBuilder();
2083       sb.append(format("\n Global Cartesian PME, lambdaMode=%s\n", lambdaMode.toString()));
2084       sb.append(format(" Multipole Self-Energy:   %16.8f\n", eself));
2085       sb.append(format(" Multipole Reciprocal:    %16.8f\n", erecip));
2086       sb.append(format(" Multipole Real Space:    %16.8f\n", ereal));
2087       sb.append(format(" Polarization Self-Energy:%16.8f\n", eselfi));
2088       sb.append(format(" Polarization Reciprocal: %16.8f\n", erecipi));
2089       sb.append(format(" Polarization Real Space: %16.8f\n", ereali));
2090       if (generalizedKirkwoodTerm) {
2091         sb.append(format(" Generalized Kirkwood:    %16.8f\n", solvationEnergy));
2092       }
2093       logger.info(sb.toString());
2094     }
2095 
2096     return permanentMultipoleEnergy + polarizationEnergy + solvationEnergy;
2097   }
2098 
2099   /**
2100    * Find the permanent multipole potential, field, etc.
2101    */
2102   private void permanentMultipoleField() {
2103     try {
2104       // Compute b-Splines and permanent density.
2105       if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2106         reciprocalSpace.computeBSplines();
2107         reciprocalSpace.splinePermanentMultipoles(globalMultipole, fractionalMultipole, use);
2108       }
2109 
2110       field.reset(parallelTeam);
2111       fieldCR.reset(parallelTeam);
2112       permanentFieldRegion.init(atoms, crystal, coordinates,
2113           globalMultipole, inducedDipole, inducedDipoleCR, neighborLists,
2114           scaleParameters, use, molecule, ipdamp, thole, ip11,
2115           mask12, mask13, mask14, lambdaMode, reciprocalSpaceTerm, reciprocalSpace,
2116           ewaldParameters, pcgSolver, permanentSchedule, realSpaceNeighborParameters,
2117           field, fieldCR, pmeTimings);
2118       // The real space contribution can be calculated at the same time
2119       // the reciprocal space convolution is being done.
2120       sectionTeam.execute(permanentFieldRegion);
2121 
2122       pmeTimings.realSpacePermTotal = permanentFieldRegion.getRealSpacePermTotal();
2123 
2124       // Collect the reciprocal space field.
2125       if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2126         reciprocalSpace.computePermanentPhi(cartesianMultipolePhi, fracMultipolePhi);
2127       }
2128     } catch (RuntimeException e) {
2129       String message = "Fatal exception computing the permanent multipole field.\n";
2130       logger.log(Level.WARNING, message, e);
2131       throw e;
2132     } catch (Exception e) {
2133       String message = "Fatal exception computing the permanent multipole field.\n";
2134       logger.log(Level.SEVERE, message, e);
2135     }
2136   }
2137 
2138   private void saveInducedDipolesToVacuumDipoles() {
2139     for (int i = 0; i < nAtoms; i++) {
2140       arraycopy(directDipole[i], 0, vacuumDirectDipole[i], 0, 3);
2141       arraycopy(directDipoleCR[i], 0, vacuumDirectDipoleCR[i], 0, 3);
2142       arraycopy(inducedDipole[0][i], 0, vacuumInducedDipole[0][i], 0, 3);
2143       arraycopy(inducedDipoleCR[0][i], 0, vacuumInducedDipoleCR[0][i], 0, 3);
2144       arraycopy(cartesianInducedDipolePhi[i], 0, cartesianVacuumDipolePhi[i], 0, tensorCount);
2145       arraycopy(cartesianInducedDipolePhiCR[i], 0, cartesianVacuumDipolePhiCR[i], 0,
2146           tensorCount);
2147       arraycopy(fractionalInducedDipolePhi[i], 0, fractionalVacuumDipolePhi[i], 0,
2148           tensorCount);
2149       arraycopy(fractionalInducedDipolePhiCR[i], 0, fractionalVacuumDipolePhiCR[i], 0,
2150           tensorCount);
2151       if (nSymm > 1) {
2152         for (int s = 1; s < nSymm; s++) {
2153           arraycopy(inducedDipole[s][i], 0, vacuumInducedDipole[s][i], 0, 3);
2154           arraycopy(inducedDipoleCR[s][i], 0, vacuumInducedDipoleCR[s][i], 0, 3);
2155         }
2156       }
2157     }
2158   }
2159 
2160   /**
2161    * Apply the selected polarization model (NONE, Direct or Mutual).
2162    */
2163   private int selfConsistentField(boolean print) {
2164     if (polarization == Polarization.NONE) {
2165       return -1;
2166     }
2167     long startTime = System.nanoTime();
2168 
2169     // Compute the direct induced dipoles.
2170     if (generalizedKirkwoodTerm) {
2171       pmeTimings.gkEnergyTotal = -System.nanoTime();
2172       generalizedKirkwood.computePermanentGKField();
2173       pmeTimings.gkEnergyTotal += System.nanoTime();
2174       logger.fine(
2175           format(" Computed GK permanent field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2176     }
2177     directRegion.init(
2178         atoms,
2179         polarizability,
2180         globalMultipole,
2181         cartesianMultipolePhi,
2182         field,
2183         fieldCR,
2184         generalizedKirkwoodTerm,
2185         generalizedKirkwood,
2186         ewaldParameters,
2187         soluteDielectric,
2188         inducedDipole,
2189         inducedDipoleCR,
2190         directDipole,
2191         directDipoleCR,
2192         directField,
2193         directFieldCR
2194     );
2195     directRegion.executeWith(parallelTeam);
2196 
2197     // Return unless mutual polarization is selected.
2198     if (polarization != Polarization.MUTUAL) {
2199       expandInducedDipoles();
2200       return 0;
2201     }
2202 
2203     // Predict the current self-consistent induced dipoles using information from previous steps.
2204     if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
2205       switch (scfPredictorParameters.scfPredictor) {
2206         case ASPC -> scfPredictorParameters.aspcPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2207         case LS -> scfPredictorParameters.leastSquaresPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2208         case POLY -> scfPredictorParameters.polynomialPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2209         default -> {
2210         }
2211       }
2212     }
2213 
2214     // Expand the initial induced dipoles to P1 symmetry, if necessary.
2215     expandInducedDipoles();
2216 
2217     // Converge the self-consistent field.
2218     try {
2219       return switch (scfAlgorithm) {
2220         case SOR -> scfBySOR(print, startTime);
2221         case EPT -> scfByEPT(print, startTime);
2222         default -> {
2223           // PCG
2224           pcgSolver.init(atoms, coordinates, polarizability, ipdamp, thole,
2225               use, crystal, inducedDipole, inducedDipoleCR, directDipole, directDipoleCR,
2226               field, fieldCR, ewaldParameters, soluteDielectric, parallelTeam,
2227               realSpaceNeighborParameters.realSpaceSchedule,
2228               pmeTimings.realSpaceSCFTime);
2229           yield pcgSolver.scfByPCG(print, startTime, this);
2230         }
2231       };
2232     } catch (EnergyException ex) {
2233       if (directFallback) {
2234         // SCF Failure: warn and revert to direct polarization.
2235         logger.warning(ex.toString());
2236         // Compute the direct induced dipoles.
2237         if (generalizedKirkwoodTerm) {
2238           pmeTimings.gkEnergyTotal = -System.nanoTime();
2239           generalizedKirkwood.computePermanentGKField();
2240           pmeTimings.gkEnergyTotal += System.nanoTime();
2241           logger.fine(
2242               format(" Computed GK permanent field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2243         }
2244         directRegion.init(
2245             atoms,
2246             polarizability,
2247             globalMultipole,
2248             cartesianMultipolePhi,
2249             field,
2250             fieldCR,
2251             generalizedKirkwoodTerm,
2252             generalizedKirkwood,
2253             ewaldParameters,
2254             soluteDielectric,
2255             inducedDipole,
2256             inducedDipoleCR,
2257             directDipole,
2258             directDipoleCR,
2259             directField,
2260             directFieldCR
2261         );
2262         directRegion.executeWith(parallelTeam);
2263         expandInducedDipoles();
2264         logger.info(" Direct induced dipoles computed due to SCF failure.");
2265         return 0;
2266       } else {
2267         throw ex;
2268       }
2269     }
2270   }
2271 
2272   /**
2273    * Converge the SCF using Successive Over-Relaxation (SOR).
2274    */
2275   private int scfBySOR(boolean print, long startTime) {
2276     long directTime = System.nanoTime() - startTime;
2277 
2278     // A request of 0 SCF cycles simplifies mutual polarization to direct polarization.
2279     StringBuilder sb = null;
2280     if (print) {
2281       sb = new StringBuilder("\n Self-Consistent Field\n Iter  RMS Change (Debye)  Time\n");
2282     }
2283     int completedSCFCycles = 0;
2284     int maxSCFCycles = 1000;
2285     double eps = 100.0;
2286     double previousEps;
2287     boolean done = false;
2288     while (!done) {
2289       long cycleTime = -System.nanoTime();
2290       try {
2291         if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2292           reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
2293         }
2294         field.reset(parallelTeam);
2295         fieldCR.reset(parallelTeam);
2296         inducedDipoleFieldRegion.init(
2297             atoms, crystal, use, molecule, ipdamp, thole, coordinates,
2298             realSpaceNeighborParameters, inducedDipole, inducedDipoleCR,
2299             reciprocalSpaceTerm, reciprocalSpace, lambdaMode, ewaldParameters,
2300             field, fieldCR, pmeTimings);
2301         inducedDipoleFieldRegion.executeWith(sectionTeam);
2302         pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
2303         if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2304           reciprocalSpace.computeInducedPhi(
2305               cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2306               fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
2307         }
2308 
2309         if (generalizedKirkwoodTerm) {
2310           // GK field.
2311           pmeTimings.gkEnergyTotal = -System.nanoTime();
2312           generalizedKirkwood.computeInducedGKField();
2313           pmeTimings.gkEnergyTotal += System.nanoTime();
2314           logger.fine(
2315               format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2316         }
2317 
2318         sorRegion.init(
2319             atoms,
2320             polarizability,
2321             inducedDipole,
2322             inducedDipoleCR,
2323             directDipole,
2324             directDipoleCR,
2325             cartesianInducedDipolePhi,
2326             cartesianInducedDipolePhiCR,
2327             field,
2328             fieldCR,
2329             generalizedKirkwoodTerm,
2330             generalizedKirkwood,
2331             ewaldParameters);
2332         parallelTeam.execute(sorRegion);
2333 
2334         expandInducedDipoles();
2335       } catch (Exception e) {
2336         String message = "Exception computing mutual induced dipoles.";
2337         logger.log(Level.SEVERE, message, e);
2338       }
2339       completedSCFCycles++;
2340       previousEps = eps;
2341       eps = sorRegion.getEps();
2342       eps = Constants.ELEC_ANG_TO_DEBYE * sqrt(eps / (double) nAtoms);
2343       cycleTime += System.nanoTime();
2344       if (print) {
2345         sb.append(format(" %4d     %15.10f %7.4f\n", completedSCFCycles, eps, cycleTime * NS2SEC));
2346       }
2347 
2348       // If the RMS Debye change increases, fail the SCF process.
2349       if (eps > previousEps) {
2350         if (sb != null) {
2351           logger.warning(sb.toString());
2352         }
2353         String message = format(" SCF convergence failure: (%10.5f > %10.5f)\n", eps, previousEps);
2354         throw new EnergyException(message);
2355       }
2356 
2357       // The SCF should converge well before the max iteration check. Otherwise, fail the SCF
2358       // process.
2359       if (completedSCFCycles >= maxSCFCycles) {
2360         if (sb != null) {
2361           logger.warning(sb.toString());
2362         }
2363         String message = format(" Maximum SCF iterations reached: (%d)\n", completedSCFCycles);
2364         throw new EnergyException(message);
2365       }
2366 
2367       // Check if the convergence criteria has been achieved.
2368       if (eps < poleps) {
2369         done = true;
2370       }
2371     }
2372     if (print) {
2373       sb.append(format(" Direct:                  %7.4f\n", NS2SEC * directTime));
2374       startTime = System.nanoTime() - startTime;
2375       sb.append(format(" Total:                   %7.4f", startTime * NS2SEC));
2376       logger.info(sb.toString());
2377     }
2378     return completedSCFCycles;
2379   }
2380 
2381   /**
2382    * Set the induced dipoles using extrapolated perturbation theory.
2383    */
2384   private int scfByEPT(boolean print, long startTime) {
2385 
2386     // Zeroth order OPT dipoles are the direct dipoles.
2387     for (int i = 0; i < nAtoms; i++) {
2388       for (int j = 0; j < 3; j++) {
2389         optRegion.optDipole[0][i][j] = directDipole[i][j];
2390         optRegion.optDipoleCR[0][i][j] = directDipoleCR[i][j];
2391       }
2392     }
2393 
2394     int optOrder = optRegion.optOrder;
2395     // Collect OPT dipole contributions from 1st to Nth.
2396     for (int currentOptOrder = 1; currentOptOrder <= optOrder; currentOptOrder++) {
2397       try {
2398         if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2399           reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
2400         }
2401         field.reset(parallelTeam);
2402         fieldCR.reset(parallelTeam);
2403         inducedDipoleFieldRegion.init(
2404             atoms, crystal, use, molecule, ipdamp, thole, coordinates, realSpaceNeighborParameters,
2405             inducedDipole, inducedDipoleCR, reciprocalSpaceTerm, reciprocalSpace, lambdaMode,
2406             ewaldParameters, field, fieldCR, pmeTimings);
2407         inducedDipoleFieldRegion.executeWith(sectionTeam);
2408         pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
2409 
2410         if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2411           reciprocalSpace.computeInducedPhi(
2412               cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2413               fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
2414         }
2415 
2416         if (generalizedKirkwoodTerm) {
2417           // GK field.
2418           pmeTimings.gkEnergyTotal = -System.nanoTime();
2419           generalizedKirkwood.computeInducedGKField();
2420           pmeTimings.gkEnergyTotal += System.nanoTime();
2421           logger.fine(
2422               format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2423         }
2424 
2425         optRegion.init(
2426             currentOptOrder,
2427             atoms,
2428             polarizability,
2429             inducedDipole,
2430             inducedDipoleCR,
2431             cartesianInducedDipolePhi,
2432             cartesianInducedDipolePhiCR,
2433             field,
2434             fieldCR,
2435             generalizedKirkwoodTerm,
2436             generalizedKirkwood,
2437             ewaldParameters,
2438             soluteDielectric);
2439         parallelTeam.execute(optRegion);
2440 
2441         expandInducedDipoles();
2442 
2443       } catch (Exception e) {
2444         String message = "Exception computing opt induced dipoles.";
2445         logger.log(Level.SEVERE, message, e);
2446       }
2447     }
2448 
2449     // Use OPT dipole components to construct Induced Dipoles for the asymmetric unit.
2450     for (int i = 0; i < nAtoms; i++) {
2451       for (int j = 0; j < 3; j++) {
2452         inducedDipole[0][i][j] = 0.0;
2453         inducedDipoleCR[0][i][j] = 0.0;
2454         double sum = 0.0;
2455         double sump = 0.0;
2456         for (int k = 0; k <= optOrder; k++) {
2457           sum += optRegion.optDipole[k][i][j];
2458           sump += optRegion.optDipoleCR[k][i][j];
2459           inducedDipole[0][i][j] += optRegion.optCoefficients[k] * sum;
2460           inducedDipoleCR[0][i][j] += optRegion.optCoefficients[k] * sump;
2461         }
2462       }
2463     }
2464 
2465     // Expand asymmetric OPT dipoles for crystal symmetry.
2466     expandInducedDipoles();
2467 
2468     return optOrder;
2469   }
2470 
2471   /**
2472    * Getter for the field <code>totalMultipoleEnergy</code>.
2473    *
2474    * @return a double.
2475    */
2476   public double getTotalMultipoleEnergy() {
2477     return permanentMultipoleEnergy + polarizationEnergy;
2478   }
2479 
2480   /**
2481    * Setter for the field <code>polarization</code>.
2482    *
2483    * @param set a {@link Polarization} object.
2484    */
2485   public void setPolarization(Polarization set) {
2486     this.polarization = set;
2487   }
2488 
2489   /**
2490    * Given an array of atoms (with atom types), assign multipole types and reference sites.
2491    */
2492   private void assignMultipoles() {
2493     if (forceField == null) {
2494       String message = "No force field is defined.\n";
2495       logger.log(Level.SEVERE, message);
2496       return;
2497     }
2498 
2499     if (forceField.getForceFieldTypeCount(ForceFieldType.MULTIPOLE) < 1
2500         && forceField.getForceFieldTypeCount(ForceFieldType.CHARGE) < 1) {
2501       String message = "Force field has no permanent electrostatic types.\n";
2502       logger.log(Level.SEVERE, message);
2503       return;
2504     }
2505     if (nAtoms < 1) {
2506       String message = "No atoms are defined.\n";
2507       logger.log(Level.SEVERE, message);
2508       return;
2509     }
2510     for (int i = 0; i < nAtoms; i++) {
2511       Atom atom = atoms[i];
2512 
2513       if (!assignMultipole(elecForm, atom, forceField, localMultipole[i], i, axisAtom, frame)) {
2514         logger.info(
2515             format("No MultipoleType could be assigned:\n %s --> %s", atom, atom.getAtomType()));
2516         StringBuilder sb = new StringBuilder();
2517         List<Bond> bonds = atom.getBonds();
2518         for (Bond bond12 : bonds) {
2519           Atom a2 = bond12.get1_2(atom);
2520           AtomType aType2 = a2.getAtomType();
2521           sb.append(format("\n  1-2 %s --> %s", a2, aType2));
2522         }
2523         for (Bond bond12 : bonds) {
2524           Atom atom2 = bond12.get1_2(atom);
2525           bonds = atom2.getBonds();
2526           for (Bond bond23 : bonds) {
2527             Atom a2 = bond23.get1_2(atom2);
2528             AtomType aType2 = a2.getAtomType();
2529             sb.append(format("\n  1-3 %s --> %s", a2, aType2));
2530           }
2531         }
2532 
2533         List<MultipoleType> multipoleTypes = forceField.getMultipoleTypes(atom.getAtomType().getKey());
2534         if (multipoleTypes != null && !multipoleTypes.isEmpty()) {
2535           sb.append("\n Similar Multipole types:");
2536           for (MultipoleType multipoleType : multipoleTypes) {
2537             sb.append(format("\n %s", multipoleType));
2538           }
2539         }
2540 
2541         logger.log(Level.SEVERE, sb.toString());
2542       }
2543     }
2544 
2545     // Check for multipoles that were not assigned correctly.
2546     StringBuilder sb = new StringBuilder();
2547     for (int i = 0; i < nAtoms; i++) {
2548       boolean flag = false;
2549       for (int j = 0; j < 10; j++) {
2550         if (Double.isNaN(localMultipole[i][j])) {
2551           flag = true;
2552           break;
2553         }
2554       }
2555       if (flag) {
2556         sb.append("\n").append(atoms[i].toString()).append("\n");
2557         sb.append(format("%d", i + 1));
2558         for (int j = 0; j < 10; j++) {
2559           sb.append(format(" %8.3f", localMultipole[i][j]));
2560         }
2561         sb.append("\n");
2562       }
2563     }
2564     if (!sb.isEmpty()) {
2565       String message = "Fatal exception: Error assigning multipoles. " + sb;
2566       logger.log(Level.SEVERE, message);
2567     }
2568   }
2569 
2570   /**
2571    * AssignPolarizationGroups.
2572    */
2573   protected void assignPolarizationGroups() {
2574     // Find directly connected group members for each atom.
2575     List<Integer> group = new ArrayList<>();
2576     for (int i = 0; i < nAtoms; i++) {
2577       Atom a = atoms[i];
2578       if (a.getIndex() - 1 != i) {
2579         logger.info(format(" PME Index: %d Atom Index: %d\n Atom: %s",
2580             i, a.getIndex() - 1, a));
2581         logger.severe(" Atom indexing is not consistent in PME.");
2582       }
2583     }
2584     for (Atom ai : atoms) {
2585       group.clear();
2586       int index = ai.getIndex() - 1;
2587       group.add(index);
2588       PolarizeType polarizeType = ai.getPolarizeType();
2589       if (polarizeType != null) {
2590         if (polarizeType.polarizationGroup != null) {
2591           PolarizeType.growGroup(group, ai);
2592           sort(group);
2593           ip11[index] = new int[group.size()];
2594           int j = 0;
2595           for (int k : group) {
2596             ip11[index][j++] = k;
2597           }
2598         } else {
2599           ip11[index] = new int[group.size()];
2600           int j = 0;
2601           for (int k : group) {
2602             ip11[index][j++] = k;
2603           }
2604         }
2605       } else {
2606         String message = "The polarize keyword was not found for atom "
2607             + (index + 1) + " with type " + ai.getType();
2608         logger.severe(message);
2609       }
2610     }
2611     // Find 1-2 group relationships.
2612     int[] mask = new int[nAtoms];
2613     List<Integer> list = new ArrayList<>();
2614     List<Integer> keep = new ArrayList<>();
2615     for (int i = 0; i < nAtoms; i++) {
2616       mask[i] = -1;
2617     }
2618     for (int i = 0; i < nAtoms; i++) {
2619       list.clear();
2620       for (int j : ip11[i]) {
2621         list.add(j);
2622         mask[j] = i;
2623       }
2624       keep.clear();
2625       for (int j : list) {
2626         Atom aj = atoms[j];
2627         List<Bond> bonds = aj.getBonds();
2628         for (Bond b : bonds) {
2629           Atom ak = b.get1_2(aj);
2630           int k = ak.getIndex() - 1;
2631           if (mask[k] != i) {
2632             keep.add(k);
2633           }
2634         }
2635       }
2636       list.clear();
2637       for (int j : keep) {
2638         for (int k : ip11[j]) {
2639           list.add(k);
2640         }
2641       }
2642       sort(list);
2643       ip12[i] = new int[list.size()];
2644       int j = 0;
2645       for (int k : list) {
2646         ip12[i][j++] = k;
2647       }
2648     }
2649 
2650     // Find 1-3 group relationships.
2651     for (int i = 0; i < nAtoms; i++) {
2652       mask[i] = -1;
2653     }
2654     for (int i = 0; i < nAtoms; i++) {
2655       for (int j : ip11[i]) {
2656         mask[j] = i;
2657       }
2658       for (int j : ip12[i]) {
2659         mask[j] = i;
2660       }
2661       list.clear();
2662       for (int j : ip12[i]) {
2663         for (int k : ip12[j]) {
2664           if (mask[k] != i) {
2665             if (!list.contains(k)) {
2666               list.add(k);
2667             }
2668           }
2669         }
2670       }
2671       ip13[i] = new int[list.size()];
2672       sort(list);
2673       int j = 0;
2674       for (int k : list) {
2675         ip13[i][j++] = k;
2676       }
2677     }
2678   }
2679 
2680   /**
2681    * Compute multipole moments for an array of atoms.
2682    *
2683    * @param activeAtoms Atom array to consider.
2684    * @param forceEnergy Force calculation of the electrostatic energy (rotate multipoles, perform
2685    *                    SCF).
2686    */
2687   public void computeMoments(Atom[] activeAtoms, boolean forceEnergy) {
2688     // Zero out total charge, dipole and quadrupole components.
2689     var netchg = 0.0;
2690     var netdpl = 0.0;
2691     var xdpl = 0.0;
2692     var ydpl = 0.0;
2693     var zdpl = 0.0;
2694     var xxqdp = 0.0;
2695     var xyqdp = 0.0;
2696     var xzqdp = 0.0;
2697     var yxqdp = 0.0;
2698     var yyqdp = 0.0;
2699     var yzqdp = 0.0;
2700     var zxqdp = 0.0;
2701     var zyqdp = 0.0;
2702     var zzqdp = 0.0;
2703 
2704     // Find the center of mass of the set of active atoms.
2705     double xmid = 0.0;
2706     double ymid = 0.0;
2707     double zmid = 0.0;
2708     double totalMass = 0;
2709     for (Atom atom : activeAtoms) {
2710       var m = atom.getMass();
2711       totalMass += m;
2712       xmid = xmid + atom.getX() * m;
2713       ymid = ymid + atom.getY() * m;
2714       zmid = zmid + atom.getZ() * m;
2715     }
2716     if (totalMass > 0) {
2717       xmid /= totalMass;
2718       ymid /= totalMass;
2719       zmid /= totalMass;
2720     }
2721     int n = activeAtoms.length;
2722     double[] xcm = new double[n];
2723     double[] ycm = new double[n];
2724     double[] zcm = new double[n];
2725     int k = 0;
2726     for (Atom atom : activeAtoms) {
2727       xcm[k] = atom.getX() - xmid;
2728       ycm[k] = atom.getY() - ymid;
2729       zcm[k] = atom.getZ() - zmid;
2730       k++;
2731     }
2732 
2733     if (forceEnergy) {
2734       energy(false, false);
2735     }
2736 
2737     // Account for charge, dipoles and induced dipoles.
2738     k = 0;
2739     for (Atom atom : activeAtoms) {
2740       int i = atom.getIndex() - 1;
2741       double[] globalMultipolei = globalMultipole[0][i];
2742       double[] inducedDipolei = inducedDipole[0][i];
2743 
2744       var ci = globalMultipolei[t000];
2745       var dix = globalMultipolei[t100];
2746       var diy = globalMultipolei[t010];
2747       var diz = globalMultipolei[t001];
2748       var uix = inducedDipolei[0];
2749       var uiy = inducedDipolei[1];
2750       var uiz = inducedDipolei[2];
2751 
2752       netchg += ci;
2753       xdpl += xcm[k] * ci + dix + uix;
2754       ydpl += ycm[k] * ci + diy + uiy;
2755       zdpl += zcm[k] * ci + diz + uiz;
2756       xxqdp += xcm[k] * xcm[k] * ci + 2.0 * xcm[k] * (dix + uix);
2757       xyqdp += xcm[k] * ycm[k] * ci + xcm[k] * (diy + uiy) + ycm[k] * (dix + uix);
2758       xzqdp += xcm[k] * zcm[k] * ci + xcm[k] * (diz + uiz) + zcm[k] * (dix + uix);
2759       yxqdp += ycm[k] * xcm[k] * ci + ycm[k] * (dix + uix) + xcm[k] * (diy + uiy);
2760       yyqdp += ycm[k] * ycm[k] * ci + 2.0 * ycm[k] * (diy + uiy);
2761       yzqdp += ycm[k] * zcm[k] * ci + ycm[k] * (diz + uiz) + zcm[k] * (diy + uiy);
2762       zxqdp += zcm[k] * xcm[k] * ci + zcm[k] * (dix + uix) + xcm[k] * (diz + uiz);
2763       zyqdp += zcm[k] * ycm[k] * ci + zcm[k] * (diy + uiy) + ycm[k] * (diz + uiz);
2764       zzqdp += zcm[k] * zcm[k] * ci + 2.0 * zcm[k] * (diz + uiz);
2765       k++;
2766     }
2767 
2768     // Convert the quadrupole from traced to traceless form.
2769     var qave = (xxqdp + yyqdp + zzqdp) / 3.0;
2770     xxqdp = 1.5 * (xxqdp - qave);
2771     xyqdp = 1.5 * xyqdp;
2772     xzqdp = 1.5 * xzqdp;
2773     yxqdp = 1.5 * yxqdp;
2774     yyqdp = 1.5 * (yyqdp - qave);
2775     yzqdp = 1.5 * yzqdp;
2776     zxqdp = 1.5 * zxqdp;
2777     zyqdp = 1.5 * zyqdp;
2778     zzqdp = 1.5 * (zzqdp - qave);
2779 
2780     // Add the traceless atomic quadrupoles to total quadrupole.
2781     for (Atom atom : activeAtoms) {
2782       int i = atom.getIndex() - 1;
2783       double[] globalMultipolei = globalMultipole[0][i];
2784       var qixx = globalMultipolei[t200];
2785       var qiyy = globalMultipolei[t020];
2786       var qizz = globalMultipolei[t002];
2787       var qixy = globalMultipolei[t110];
2788       var qixz = globalMultipolei[t101];
2789       var qiyz = globalMultipolei[t011];
2790       xxqdp += qixx;
2791       xyqdp += qixy;
2792       xzqdp += qixz;
2793       yxqdp += qixy;
2794       yyqdp += qiyy;
2795       yzqdp += qiyz;
2796       zxqdp += qixz;
2797       zyqdp += qiyz;
2798       zzqdp += qizz;
2799     }
2800 
2801     // Convert dipole to Debye and quadrupole to Buckingham.
2802     xdpl = xdpl * ELEC_ANG_TO_DEBYE;
2803     ydpl = ydpl * ELEC_ANG_TO_DEBYE;
2804     zdpl = zdpl * ELEC_ANG_TO_DEBYE;
2805     xxqdp = xxqdp * ELEC_ANG_TO_DEBYE;
2806     xyqdp = xyqdp * ELEC_ANG_TO_DEBYE;
2807     xzqdp = xzqdp * ELEC_ANG_TO_DEBYE;
2808     yxqdp = yxqdp * ELEC_ANG_TO_DEBYE;
2809     yyqdp = yyqdp * ELEC_ANG_TO_DEBYE;
2810     yzqdp = yzqdp * ELEC_ANG_TO_DEBYE;
2811     zxqdp = zxqdp * ELEC_ANG_TO_DEBYE;
2812     zyqdp = zyqdp * ELEC_ANG_TO_DEBYE;
2813     zzqdp = zzqdp * ELEC_ANG_TO_DEBYE;
2814 
2815     // Get dipole magnitude and diagonalize quadrupole tensor.
2816     netdpl = sqrt(xdpl * xdpl + ydpl * ydpl + zdpl * zdpl);
2817     double[][] a = new double[3][3];
2818     a[0][0] = xxqdp;
2819     a[0][1] = xyqdp;
2820     a[0][2] = xzqdp;
2821     a[1][0] = yxqdp;
2822     a[1][1] = yyqdp;
2823     a[1][2] = yzqdp;
2824     a[2][0] = zxqdp;
2825     a[2][1] = zyqdp;
2826     a[2][2] = zzqdp;
2827     EigenDecomposition e = new EigenDecomposition(new Array2DRowRealMatrix(a));
2828     // Eigenvalues are returned in descending order, but logged below in ascending order.
2829     var netqdp = e.getRealEigenvalues();
2830 
2831     logger.info("\n Electric Moments\n");
2832     logger.info(format("  Total Electric Charge:    %13.5f Electrons\n", netchg));
2833     logger.info(format("  Dipole Moment Magnitude:  %13.5f Debye\n", netdpl));
2834     logger.info(format("  Dipole X,Y,Z-Components:  %13.5f %13.5f %13.5f\n", xdpl, ydpl, zdpl));
2835     logger.info(format("  Quadrupole Moment Tensor: %13.5f %13.5f %13.5f", xxqdp, xyqdp, xzqdp));
2836     logger.info(format("       (Buckinghams)        %13.5f %13.5f %13.5f", yxqdp, yyqdp, yzqdp));
2837     logger.info(format("                            %13.5f %13.5f %13.5f\n", zxqdp, zyqdp, zzqdp));
2838     logger.info(
2839         format(
2840             "  Principal Axes Quadrupole %13.5f %13.5f %13.5f\n", netqdp[2], netqdp[1], netqdp[0]));
2841   }
2842 
2843   /**
2844    * Log the real space electrostatics interaction.
2845    *
2846    * @param i   Atom i.
2847    * @param k   Atom j.
2848    * @param r   The distance rij.
2849    * @param eij The interaction energy.
2850    * @since 1.0
2851    */
2852   private void log(int i, int k, double r, double eij) {
2853     logger.info(
2854         format(
2855             "%s %6d-%s %6d-%s %10.4f  %10.4f",
2856             "ELEC",
2857             atoms[i].getIndex(),
2858             atoms[i].getAtomType().name,
2859             atoms[k].getIndex(),
2860             atoms[k].getAtomType().name,
2861             r,
2862             eij));
2863   }
2864 
2865   /**
2866    * Flag to indicate use of extended variables that interpolate permanent multipoles and/or atomic
2867    * polarizability between states (e.g., for protonation and/or tautamers).
2868    */
2869   private boolean esvTerm = false;
2870 
2871   /**
2872    * The ExtendedSystem instance controls use of extended variables.
2873    */
2874   private ExtendedSystem extendedSystem = null;
2875 
2876   /**
2877    * Attach system with extended variable such as titrations.
2878    *
2879    * @param system a {@link ffx.potential.extended.ExtendedSystem} object.
2880    */
2881   public void attachExtendedSystem(ExtendedSystem system) {
2882     // Set object handles.
2883     esvTerm = true;
2884     extendedSystem = system;
2885     // Update atoms and reinitialize arrays for consistency with the ExtendedSystem.
2886     setAtoms(extendedSystem.getExtendedAtoms(), extendedSystem.getExtendedMolecule());
2887     // Allocate space for dM/dTitratonESV
2888     if (dMultipoledTirationESV == null || dMultipoledTirationESV.length != nSymm
2889         || dMultipoledTirationESV[0].length != nAtoms) {
2890       dMultipoledTirationESV = new double[nSymm][nAtoms][10];
2891       dMultipoledTautomerESV = new double[nSymm][nAtoms][10];
2892     }
2893 
2894     if (dPolardTitrationESV == null || dPolardTitrationESV.length != nAtoms) {
2895       dPolardTitrationESV = new double[nAtoms];
2896       dPolardTautomerESV = new double[nAtoms];
2897     }
2898 
2899     logger.info(format(" Attached extended system (%d variables) to PME.\n", extendedSystem.getNumberOfVariables()));
2900   }
2901 
2902   /**
2903    * The scalar ESV that is operating on each atom (or 1.0 if the atom is not under ESV control).
2904    */
2905   double[] perAtomTitrationESV = null;
2906 
2907   /**
2908    * The partial derivative of each multipole with respect to its titration/tautomer ESV (or 0.0 if
2909    * the atom is not under titration ESV control).
2910    */
2911   double[][][] dMultipoledTirationESV = null;
2912   double[][][] dMultipoledTautomerESV = null;
2913 
2914   /**
2915    * The partial derivative of each polarizability with respect to its titration/tautomer ESV (or 0.0
2916    * if the atom is not under titration ESV control).
2917    */
2918   double[] dPolardTitrationESV = null;
2919   double[] dPolardTautomerESV = null;
2920   /**
2921    * OST and ESV specific factors that effect real space interactions.
2922    */
2923   LambdaFactors[] lambdaFactors = null;
2924 
2925   /**
2926    * The setFactors(i,k,lambdaMode) method is called every time through the inner PME loops, avoiding
2927    * an "if (esv)" branch statement.
2928    * <p>
2929    * A plain OST run will have an object of type LambdaFactorsOST instead, which contains an empty
2930    * version of setFactors(i,k,lambdaMode).
2931    * <p>
2932    * The OST version instead sets new factors only on lambda updates, in setLambda(i,k).
2933    * <p>
2934    * Interactions involving neither lambda receive the (inoperative) defaults below.
2935    */
2936   public class LambdaFactors {
2937 
2938     /**
2939      * lambda * esvLambda[i] * esvLambda[k]
2940      */
2941     protected double lambdaProd = 1.0;
2942     /**
2943      * Interatomic buffer distance: alpha*(1-lambda)*(1-lambda).
2944      */
2945     protected double lfAlpha = 0.0;
2946     /**
2947      * First lambda derivative of buffer distance.
2948      */
2949     protected double dlfAlpha = 0.0;
2950     /**
2951      * Second lambda derivative of buffer distance.
2952      */
2953     protected double d2lfAlpha = 0.0;
2954     /**
2955      * Lambda to its permanent exponent.
2956      */
2957     protected double lfPowPerm = 1.0;
2958     /**
2959      * First lambda derivative of lPowPerm.
2960      */
2961     protected double dlfPowPerm = 0.0;
2962     /**
2963      * Second lambda derivative of lPowPerm.
2964      */
2965     protected double d2lfPowPerm = 0.0;
2966     /**
2967      * Lambda to its polarization exponent.
2968      */
2969     protected double lfPowPol = 1.0;
2970     /**
2971      * First lambda derivative of lPowPol.
2972      */
2973     protected double dlfPowPol = 0.0;
2974     /**
2975      * Second lambda derivative of lPowPol.
2976      */
2977     protected double d2lfPowPol = 0.0;
2978     /**
2979      * Derivative of lambdaProduct w.r.t. lambda.
2980      */
2981     protected double dLpdL = 1.0;
2982     /**
2983      * Derivative of lambdaProduct w.r.t. esvLambda[i].
2984      */
2985     protected double dLpdLi = 1.0;
2986     /**
2987      * Derivative of lambdaProduct w.r.t. esvLambda[k].
2988      */
2989     protected double dLpdLk = 1.0;
2990     /**
2991      * Atom indices used only by LambdaFactorsESV::print.
2992      */
2993     protected int[] ik = new int[2];
2994 
2995     public void print() {
2996       StringBuilder sb = new StringBuilder();
2997       if (this instanceof LambdaFactorsESV) {
2998         sb.append(format("  (QI-ESV)  i,k:%d,%d", ik[0], ik[1]));
2999       } else {
3000         sb.append("  (QI-OST)");
3001       }
3002       sb.append(format(
3003           "  lambda:%.2f  lAlpha:%.2f,%.2f,%.2f  lPowPerm:%.2f,%.2f,%.2f  lPowPol:%.2f,%.2f,%.2f",
3004           lambdaProd,
3005           lfAlpha,
3006           dlfAlpha,
3007           d2lfAlpha,
3008           lfPowPerm,
3009           dlfPowPerm,
3010           d2lfPowPerm,
3011           lfPowPol,
3012           dlfPowPol,
3013           d2lfPowPol));
3014 
3015       sb.append(format(
3016           "\n    permExp:%.2f  permAlpha:%.2f  permWindow:%.2f,%.2f  polExp:%.2f  polWindow:%.2f,%.2f",
3017           alchemicalParameters.permLambdaExponent,
3018           alchemicalParameters.permLambdaAlpha,
3019           alchemicalParameters.permLambdaStart,
3020           alchemicalParameters.permLambdaEnd,
3021           alchemicalParameters.polLambdaExponent,
3022           alchemicalParameters.polLambdaStart,
3023           alchemicalParameters.polLambdaEnd));
3024       logger.info(sb.toString());
3025     }
3026 
3027     /**
3028      * Overriden by the ESV version which updates with every softcore interaction.
3029      *
3030      * @param i    Atom i index.
3031      * @param k    Atom k index.
3032      * @param mode the LambdaMode
3033      */
3034     public void setFactors(int i, int k, LambdaMode mode) {
3035       /* no-op */
3036     }
3037 
3038     /**
3039      * Overriden by the OST version which updates only during setLambda().
3040      */
3041     public void setFactors() {
3042       /* no-op */
3043     }
3044   }
3045 
3046   public class LambdaFactorsOST extends LambdaFactors {
3047 
3048     @Override
3049     public void setFactors() {
3050       lambdaProd = lambda;
3051       lfAlpha = alchemicalParameters.lAlpha;
3052       dlfAlpha = alchemicalParameters.dlAlpha;
3053       d2lfAlpha = alchemicalParameters.d2lAlpha;
3054       lfPowPerm = alchemicalParameters.permanentScale;
3055       dlfPowPerm =
3056           alchemicalParameters.dlPowPerm * alchemicalParameters.dEdLSign;
3057       d2lfPowPerm = alchemicalParameters.d2lPowPerm
3058           * alchemicalParameters.dEdLSign;
3059       lfPowPol = alchemicalParameters.polarizationScale;
3060       dlfPowPol =
3061           alchemicalParameters.dlPowPol * alchemicalParameters.dEdLSign;
3062       d2lfPowPol =
3063           alchemicalParameters.d2lPowPol * alchemicalParameters.dEdLSign;
3064     }
3065   }
3066 
3067   public class LambdaFactorsESV extends LambdaFactors {
3068 
3069     @Override
3070     public void setFactors(int i, int k, LambdaMode mode) {
3071       logger.info(format("Invoked Qi setFactors() method with i,k=%d,%d", i, k));
3072       ik[0] = i;
3073       ik[1] = k;
3074       final double L = lambda;
3075       lambdaProd = L * perAtomTitrationESV[i] * perAtomTitrationESV[k];
3076       final double esvli = perAtomTitrationESV[i];
3077       final double esvlk = perAtomTitrationESV[k];
3078       dLpdL = esvli * esvlk;
3079       dLpdLi = L * esvlk;
3080       dLpdLk = L * esvli;
3081 
3082       // Permanent Lambda Window (e.g., 0 .. 1).
3083       double permLambdaExponent = alchemicalParameters.permLambdaExponent;
3084       double permLambdaStart = alchemicalParameters.permLambdaStart;
3085       double permLambdaEnd = alchemicalParameters.permLambdaEnd;
3086 
3087       double permWindow = 1.0 / (permLambdaEnd - permLambdaStart);
3088       double permLambda = (lambdaProd - permLambdaStart) * permWindow;
3089       lfPowPerm = pow(permLambda, alchemicalParameters.permLambdaExponent);
3090       dlfPowPerm = (permLambdaExponent < 1)
3091           ? 0.0 : permLambdaExponent * pow(permLambda, permLambdaExponent - 1) * permWindow;
3092       d2lfPowPerm = (permLambdaExponent < 2)
3093           ? 0.0 : permLambdaExponent
3094           * (permLambdaExponent - 1)
3095           * pow(permLambda, permLambdaExponent - 2)
3096           * permWindow
3097           * permWindow;
3098 
3099       // Polarization Lambda Window (e.g., 0 .. 1).
3100       double polLambdaExponent = alchemicalParameters.polLambdaExponent;
3101       double polLambdaStart = alchemicalParameters.polLambdaStart;
3102       double polLambdaEnd = alchemicalParameters.polLambdaEnd;
3103 
3104       double polWindow = 1.0 / (polLambdaEnd - polLambdaStart);
3105       double polLambda = (lambdaProd - polLambdaStart) * polWindow;
3106       lfPowPol = pow(polLambda, polLambdaExponent);
3107       dlfPowPol = (polLambdaExponent < 1)
3108           ? 0.0 : polLambdaExponent * pow(polLambda, polLambdaExponent - 1) * polWindow;
3109       d2lfPowPol = (polLambdaExponent < 2)
3110           ? 0.0 : polLambdaExponent
3111           * (polLambdaExponent - 1)
3112           * pow(polLambda, polLambdaExponent - 2)
3113           * polWindow
3114           * polWindow;
3115 
3116       // Permanent Lambda Softcore Alpha.
3117       double permLambdaAlpha = alchemicalParameters.permLambdaAlpha;
3118       lfAlpha = permLambdaAlpha * (1.0 - permLambda) * (1.0 - permLambda);
3119       dlfAlpha = permLambdaAlpha * (1.0 - permLambda) * permWindow;
3120       d2lfAlpha = -permLambdaAlpha * permWindow * permWindow;
3121 
3122       /*
3123         Follow the logic of
3124         1) condensedEnergy
3125         2) noLigand
3126         and
3127         3) ligandElec
3128 
3129         to set permanentScale (lfPowPerm) and polarizationScale (lfPowPol),
3130         substituting lambdaProduct for lambda.
3131        */
3132       switch (mode) {
3133         case CONDENSED -> {
3134           lfPowPerm = 1.0 - lfPowPerm;
3135           dlfPowPerm = -dlfPowPerm; // handles dEdLSign
3136           d2lfPowPerm = -d2lfPowPerm; // handles dEdLSign
3137           if (polarization == Polarization.NONE || lambda < polLambdaStart) {
3138             lfPowPol = 0.0;
3139             dlfPowPol = 0.0;
3140             d2lfPowPol = 0.0;
3141           } else if (lambda <= polLambdaEnd) {
3142             // No change necessary.
3143           } else {
3144             lfPowPol = 1.0;
3145             dlfPowPol = 0.0;
3146             d2lfPowPol = 0.0;
3147           }
3148         }
3149         case CONDENSED_NO_LIGAND -> {
3150           // There's no treatment of polLambdaStart in ::condensedNoLigandSCF?
3151           lfPowPerm = 1.0 - lfPowPerm;
3152           dlfPowPerm = -dlfPowPerm;
3153           d2lfPowPerm = -d2lfPowPerm;
3154           if (polarization == Polarization.NONE) {
3155             lfPowPol = 0.0;
3156             dlfPowPol = 0.0;
3157             d2lfPowPol = 0.0;
3158           } else if (lambda <= polLambdaEnd) {
3159             lfPowPol = 1.0 - lfPowPol;
3160             dlfPowPol = -dlfPowPol;
3161             d2lfPowPol = -d2lfPowPol;
3162           } else {
3163             lfPowPol = 0.0;
3164             dlfPowPol = 0.0;
3165             d2lfPowPol = 0.0;
3166           }
3167         }
3168         case VAPOR -> {
3169           // There's no treatment of polLambdaStart in ::ligandElec?
3170           lfPowPerm = 1.0 - lfPowPerm;
3171           dlfPowPerm = -dlfPowPerm;
3172           d2lfPowPerm = -d2lfPowPerm;
3173           lfAlpha = 0.0;
3174           dlfAlpha = 0.0;
3175           d2lfAlpha = 0.0;
3176           if (polarization == Polarization.NONE || lambdaProd > polLambdaEnd) {
3177             lfPowPol = 0.0;
3178             dlfPowPol = 0.0;
3179             d2lfPowPol = 0.0;
3180           } else if (lambdaProd <= polLambdaEnd) {
3181             lfPowPol = 1.0 - lfPowPol;
3182             dlfPowPol = -dlfPowPol;
3183             d2lfPowPol = -d2lfPowPol;
3184           }
3185         }
3186         default -> {
3187         }
3188       }
3189     }
3190   }
3191 
3192   /**
3193    * The defaults are effectively final, as the implementation of setFactors in the base class is
3194    * always a no-op.
3195    */
3196   public final LambdaFactors LambdaDefaults = new LambdaFactors();
3197 }