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