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