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.ParallelTeam;
41  import ffx.crystal.Crystal;
42  import ffx.numerics.atomic.AtomicDoubleArray;
43  import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
44  import ffx.numerics.atomic.AtomicDoubleArray3D;
45  import ffx.potential.bonded.Atom;
46  import ffx.potential.bonded.LambdaInterface;
47  import ffx.potential.nonbonded.implicit.BornGradRegion;
48  import ffx.potential.nonbonded.implicit.BornRadiiRegion;
49  import ffx.potential.nonbonded.implicit.BornTanhRescaling;
50  import ffx.potential.nonbonded.implicit.ChandlerCavitation;
51  import ffx.potential.nonbonded.implicit.ConnollyRegion;
52  import ffx.potential.nonbonded.implicit.DispersionRegion;
53  import ffx.potential.nonbonded.implicit.GKEnergyRegion;
54  import ffx.potential.nonbonded.implicit.GaussVol;
55  import ffx.potential.nonbonded.implicit.InducedGKFieldRegion;
56  import ffx.potential.nonbonded.implicit.InitializationRegion;
57  import ffx.potential.nonbonded.implicit.PermanentGKFieldRegion;
58  import ffx.potential.nonbonded.implicit.SurfaceAreaRegion;
59  import ffx.potential.nonbonded.pme.Polarization;
60  import ffx.potential.parameters.AtomType;
61  import ffx.potential.parameters.ForceField;
62  import ffx.potential.parameters.SoluteType;
63  import ffx.potential.parameters.SoluteType.SOLUTE_RADII_TYPE;
64  import ffx.utilities.Constants;
65  import ffx.utilities.FFXProperty;
66  
67  import java.util.Arrays;
68  import java.util.HashMap;
69  import java.util.logging.Level;
70  import java.util.logging.Logger;
71  
72  import static ffx.potential.nonbonded.implicit.DispersionRegion.DEFAULT_DISPERSION_OFFSET;
73  import static ffx.potential.parameters.ForceField.toEnumForm;
74  import static ffx.potential.parameters.SoluteType.setSoluteRadii;
75  import static ffx.utilities.Constants.dWater;
76  import static ffx.utilities.PropertyGroup.ImplicitSolvent;
77  import static java.lang.String.format;
78  import static java.util.Arrays.fill;
79  import static org.apache.commons.math3.util.FastMath.max;
80  
81  /**
82   * This Generalized Kirkwood class implements GK for the AMOEBA polarizable atomic multipole force
83   * field in parallel using a {@link ffx.potential.nonbonded.NeighborList}.
84   *
85   * @author Michael J. Schnieders<br> derived from:<br> TINKER code by Michael J. Schnieders and Jay
86   * W. Ponder<br>
87   * @see <a href="http://dx.doi.org/10.1021/ct7001336" target="_blank">M. J. Schnieders and J. W.
88   * Ponder, Polarizable atomic multipole solutes in a generalized Kirkwood continuum, Journal of
89   * Chemical Theory and Computation 2007, 3, (6), 2083-2097.</a><br>
90   * @since 1.0
91   */
92  public class GeneralizedKirkwood implements LambdaInterface {
93  
94    /**
95     * Default dielectric offset
96     */
97    public static final double DEFAULT_DIELECTRIC_OFFSET = 0.09;
98  
99    private static final Logger logger = Logger.getLogger(GeneralizedKirkwood.class.getName());
100   /**
101    * Default Bondi scale factor.
102    */
103   private static final double DEFAULT_SOLUTE_SCALE = 1.0;
104 
105   /**
106    * Dielectric offset from:
107    *
108    * <p>W. C. Still, A. Tempczyk, R. C. Hawley and T. Hendrickson, "A Semianalytical Treatment of
109    * Solvation for Molecular Mechanics and Dynamics", J. Amer. Chem. Soc., 112, 6127-6129 (1990)
110    */
111   private final double dOffset = DEFAULT_DIELECTRIC_OFFSET;
112   /**
113    * Force field in use.
114    */
115   private final ForceField forceField;
116   /**
117    * Treatment of polarization.
118    */
119   private final Polarization polarization;
120   /**
121    * Particle mesh Ewald instance, which contains variables such as expanded coordinates and
122    * multipoles in the global frame that GK uses.
123    */
124   private final ParticleMeshEwald particleMeshEwald;
125   /**
126    * Parallel team object for shared memory parallelization.
127    */
128   private final ParallelTeam parallelTeam;
129   /**
130    * Initialize GK output variables.
131    */
132   private final InitializationRegion initializationRegion;
133   /**
134    * Parallel computation of Born Radii.
135    */
136   private final BornRadiiRegion bornRadiiRegion;
137   /**
138    * Parallel computation of the Permanent GK Field.
139    */
140   private final PermanentGKFieldRegion permanentGKFieldRegion;
141   /**
142    * Parallel computation of the Induced GK Field.
143    */
144   private final InducedGKFieldRegion inducedGKFieldRegion;
145   /**
146    * Parallel computation of the GK continuum electrostatics energy.
147    */
148   private final GKEnergyRegion gkEnergyRegion;
149   /**
150    * Parallel computation of Born radii chain rule term.
151    */
152   private final BornGradRegion bornGradRegion;
153   /**
154    * Parallel computation of Dispersion energy.
155    */
156   private final DispersionRegion dispersionRegion;
157   /**
158    * Parallel computation of Cavitation.
159    */
160   private final SurfaceAreaRegion surfaceAreaRegion;
161   /**
162    * Volume to Surface Area Cavitation Dependence.
163    */
164   private final ChandlerCavitation chandlerCavitation;
165   /**
166    * Maps radii overrides (by AtomType) specified from the command line. e.g.
167    * -DradiiOverride=134r1.20,135r1.20 sets atom types 134,135 to Bondi=1.20
168    */
169   private final HashMap<Integer, Double> radiiOverride = new HashMap<>();
170   /**
171    * Conversion from electron**2/Ang to kcal/mole.
172    */
173   public double electric;
174   /**
175    * Empirical scaling of the Bondi radii.
176    */
177   private final double bondiScale;
178 
179   /**
180    * Array of Atoms being considered.
181    */
182   private Atom[] atoms;
183   /**
184    * Number of atoms.
185    */
186   private int nAtoms;
187   /**
188    * Cartesian coordinates of each atom.
189    */
190   private double[] x, y, z;
191   /**
192    * Base radius of each atom.
193    */
194   private double[] baseRadius;
195   /**
196    * Descreen radius of each atom.
197    */
198   private double[] descreenRadius;
199   /**
200    * Overlap scale factor for each atom, when using the Hawkins, Cramer & Truhlar pairwise
201    * descreening algorithm.
202    *
203    * <p>G. D. Hawkins, C. J. Cramer and D. G. Truhlar, "Parametrized Models of Aqueous Free Energies
204    * of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium",
205    * J. Phys. Chem., 100, 19824-19839 (1996).
206    */
207   private double[] overlapScale;
208   /**
209    * Sneck scaling parameter for each atom. Set based on maximum Sneck scaling parameter and number
210    * of bound non-hydrogen atoms
211    */
212   private double[] neckScale;
213   /**
214    * If true, the descreening integral includes overlaps with the volume of the descreened atom
215    */
216   private final boolean perfectHCTScale;
217 
218   /**
219    * The requested permittivity for the solvent.
220    */
221   @FFXProperty(name = "solvent-dielectric", propertyGroup = ImplicitSolvent, defaultValue = "78.3", description = """
222       The dielectric constant used for the solvent in generalized Kirkwood calculations.
223       The default of 78.3 corresponds to water.
224       """)
225   private final double solventDielectric;
226   /**
227    * The requested permittivity for the solute.
228    */
229   @FFXProperty(name = "solute-dielectric", propertyGroup = ImplicitSolvent, defaultValue = "1.0", description = """
230       The dielectric constant used for the solute(s) in generalized Kirkwood calculations.
231       The default of 1.0 is consistent with all solute dielectric response arising from either
232       polarization via induced dipoles and/or permanent dipole realignment.
233       """)
234   private final double soluteDielectric;
235 
236   /**
237    * Default overlap scale factor for the Hawkins, Cramer & Truhlar pairwise descreening algorithm:
238    * 0.69 New default overlap scale factor set during implicit solvent model optimization: 0.72
239    */
240   private static final double DEFAULT_HCT_SCALE = 0.72;
241 
242   /**
243    * Base overlap HCT overlap scale factor.
244    */
245   @FFXProperty(name = "hct-scale", propertyGroup = ImplicitSolvent, defaultValue = "0.72", description =
246       "The default overlap scale factor for Hawkins-Cramer-Truhlar pairwise descreening.")
247   private final double hctScale;
248 
249   /**
250    * If true, HCT overlap scale factors are element-specific
251    */
252   @FFXProperty(name = "element-hct-scale", clazz = Boolean.class, propertyGroup = ImplicitSolvent, defaultValue = "true",
253       description =
254           "Flag to turn on element specific overlap scale factors for Hawkins-Cramer-Truhlar pairwise descreening.")
255   private final boolean elementHCTScale;
256 
257   /**
258    * Element-specific HCT overlap scale factors
259    */
260   private final HashMap<Integer, Double> elementHCTScaleFactors;
261 
262   /**
263    * If true, the descreening size of atoms is based on their force field vdW radius.
264    */
265   @FFXProperty(name = "descreen-vdw", propertyGroup = ImplicitSolvent, defaultValue = "true", description =
266       "If true, the descreening size of each atom is based on its force field van der Waals radius.")
267   private final boolean descreenVDW;
268 
269   /**
270    * If true, hydrogen atoms displace solvent during the pairwise descreening integral.
271    */
272   @FFXProperty(name = "descreen-hydrogen", propertyGroup = ImplicitSolvent, defaultValue = "false", description =
273       "If true, hydrogen atoms are contribute to the pairwise descreening integrals.")
274   private final boolean descreenHydrogen;
275 
276   private static final double DEFAULT_DESCREEN_OFFSET = 0.3;
277   /**
278    * Offset applied to the pairwise descreening integral to improve stability at small separation.
279    */
280   @FFXProperty(name = "descreen-offset", propertyGroup = ImplicitSolvent, defaultValue = "0.3", description =
281       "Offset applied to the pairwise descreening integral to improve stability at small separation.")
282   private final double descreenOffset;
283 
284   /**
285    * Apply a neck correction during descreening.
286    */
287   @FFXProperty(name = "neck-correction", propertyGroup = ImplicitSolvent, defaultValue = "true", description =
288       "Apply a neck correction during descreening.")
289   private final boolean neckCorrection;
290 
291   /**
292    * Default Sneck scaling factor (for proteins)
293    */
294   private static final double DEFAULT_NECK_SCALE = 0.1350;
295   /**
296    * Maximum Sneck scaling parameter value
297    */
298   @FFXProperty(name = "neck-scale", propertyGroup = ImplicitSolvent, defaultValue = "0.1350", description =
299       "The overlap scale factor to use during the descreening neck correction.")
300   private double sneck;
301 
302   /**
303    * Use the Corrigan et al. chemically aware neck correction; atoms with more heavy atom bonds are
304    * less capable of forming interstitial necks.
305    */
306   @FFXProperty(name = "chemically-aware-neck-scale", propertyGroup = ImplicitSolvent, defaultValue = "true",
307       description = """
308           If the neck descreening correction is being used, apply a smaller overlap scale
309           factors as the number of bonded heavy atoms increases.
310           """)
311   private final boolean chemicallyAwareSneck;
312 
313   /**
314    * If true, the descreening integral includes the tanh correction to better approximate molecular
315    * surface
316    */
317   @FFXProperty(name = "tanh-correction", propertyGroup = ImplicitSolvent, defaultValue = "true", description = """
318       If the neck descreening correction is being used, apply a smaller overlap scale
319       factors as the number of bonded heavy atoms increases.
320       """)
321   private final boolean tanhCorrection;
322 
323   /**
324    * Default value of beta0 for tanh scaling
325    */
326   public static final double DEFAULT_TANH_BETA0 = 0.9563;
327   /**
328    * Default value of beta1 for tanh scaling
329    */
330   public static final double DEFAULT_TANH_BETA1 = 0.2578;
331   /**
332    * Default value of beta2 for tanh scaling
333    */
334   public static final double DEFAULT_TANH_BETA2 = 0.0810;
335 
336   /**
337    * The coefficient beta0 for tanh rescaling of descreening integrals.
338    */
339   @FFXProperty(name = "tanh-beta0", propertyGroup = ImplicitSolvent, defaultValue = "0.9563", description =
340       "The coefficient beta0 for tanh rescaling of descreening integrals.")
341   private double beta0;
342 
343   /**
344    * The coefficient beta1 for tanh rescaling of descreening integrals.
345    */
346   @FFXProperty(name = "tanh-beta1", propertyGroup = ImplicitSolvent, defaultValue = "0.2578", description =
347       "The coefficient beta1 for tanh rescaling of descreening integrals.")
348   private double beta1;
349 
350   /**
351    * The coefficient beta2 for tanh rescaling of descreening integrals.
352    */
353   @FFXProperty(name = "tanh-beta2", propertyGroup = ImplicitSolvent, defaultValue = "0.0810", description =
354       "The coefficient beta2 for tanh rescaling of descreening integrals.")
355   private double beta2;
356 
357   /**
358    * Default constant for the Generalized Kirkwood cross-term.
359    */
360   public static final double DEFAULT_GKC = 2.455;
361   /**
362    * The Generalized Kirkwood cross-term parameter.
363    */
364   @FFXProperty(name = "gkc", propertyGroup = ImplicitSolvent, defaultValue = "2.455", description =
365       "The Generalized Kirkwood cross-term parameter.")
366   public final double gkc;
367 
368   private static final NonPolarModel DEFAULT_NONPOLAR_MODEL = NonPolarModel.GAUSS_DISP;
369 
370   /**
371    * Treatment of non-polar interactions.
372    */
373   @FFXProperty(name = "nonpolar-model", clazz = String.class, propertyGroup = ImplicitSolvent,
374       defaultValue = "gauss-disp", description = """ 
375       [CAV / CAV-DISP / DISP / GAUSS-DISP / SEV-DISP / NONE ]
376       The non-polar contribution to the implicit solvent.
377       """)
378   private final NonPolarModel nonPolarModel;
379 
380   /**
381    * Default surface tension for non-polar models without an explicit dispersion term. This is lower
382    * than CAVDISP, since the favorable dispersion term is implicitly included.
383    */
384   private static final double DEFAULT_CAVONLY_SURFACE_TENSION = 0.0049;
385 
386   /**
387    * Water probe radius.
388    */
389   public final double probe;
390 
391   /**
392    * Default solvent pressure for apolar models with an explicit volume term.
393    *
394    * <p>From work by Chandler et al., the following relationship for cavitation free energy as a
395    * function of spherical particle size was found: Cross-Over = 3 * S.T. / S.P.
396    *
397    * <p>A S.P. of 0.0334 kcal/mol/A^3 was obtained using explicit AMOEBA water simulations and
398    * solvent excluded volumes.
399    *
400    * <p>A S.P. of 0.0343 kcal/mol/A^3 is obtained assuming a macroscopic surface tension of 0.103
401    * kcal/mol/A^3 and a cross-over of 9.0 (i.e. S.P. = 3 * S.T. / Cross-Over)
402    *
403    * <p>Both values are in reasonably good agreement, and 0.0334 is chosen as our default.
404    */
405   public static final double DEFAULT_SOLVENT_PRESSURE = 0.0334;
406 
407   /**
408    * Default surface tension for apolar models with an explicit dispersion term.
409    *
410    * <p>Experimental value: 0.103 kcal/mol/Ang^2
411    * <p>More physical value, used for simulations: 0.080 kcal/mol/Ang^2
412    */
413   public static final double DEFAULT_CAVDISP_SURFACE_TENSION = 0.080;
414 
415   /**
416    * Using a S.P. of 0.0334 kcal/mol/A^3, and a limiting surface tension of 0.103 kcal/mol/A^2, the
417    * cross-over point is 9.2515 A.
418    *
419    * <p>Using a S.P. of 0.0334 kcal/mol/A^3, and a limiting surface tension of 0.08 (i.e. 80% of the
420    * experimentally observed surface tension of 0.103 kcal/mol/A^2), we derive a cross-over of:
421    *
422    * <p>9.251 A = 3 * 0.103 kcal/mol/A^2 / 0.0334 kcal/mol/A^3.
423    */
424   public static final double DEFAULT_CROSSOVER =
425       3.0 * DEFAULT_CAVDISP_SURFACE_TENSION / DEFAULT_SOLVENT_PRESSURE;
426 
427   /**
428    * Cavitation surface tension coefficient (kcal/mol/A^2).
429    */
430   @FFXProperty(name = "surface-tension", propertyGroup = ImplicitSolvent, defaultValue = "0.080", description =
431       "The cavitation surface tension coefficient (kcal/mol/A^2).")
432   private final double surfaceTension;
433 
434   /**
435    * Cavitation solvent pressure coefficient (kcal/mol/A^3).
436    */
437   @FFXProperty(name = "solvent-pressure", propertyGroup = ImplicitSolvent, defaultValue = "0.0334", description =
438       "The solvent pressure for nonpolar models with an explicit volume term (kcal/mol/A^3).")
439   private final double solventPressue;
440 
441   /**
442    * The base radii to use for GK.
443    */
444   @FFXProperty(name = "gk-radius", clazz = String.class, propertyGroup = ImplicitSolvent, defaultValue = "solute",
445       description = """
446           [SOLUTE / VDW / CONSENSUS]
447           The base atomic radii to use for generalized Kirkwood calculations. The default is to use solute radii,
448           which were fit to experimental solvation free energy differences. Alternatively, force field
449           specific van der Waals radii (vdw) or consensus Bondi radii (consensus) can be chosen.
450           """)
451   private SOLUTE_RADII_TYPE soluteRadiiType;
452 
453   /**
454    * Born radius of each atom.
455    */
456   private double[] born;
457   /**
458    * Flag to indicate if an atom should be included.
459    */
460   private boolean[] use = null;
461   /**
462    * Periodic boundary conditions and symmetry.
463    */
464   private Crystal crystal;
465   /**
466    * Atomic coordinates for each symmetry operator.
467    */
468   private double[][][] sXYZ;
469   /**
470    * Multipole moments for each symmetry operator.
471    */
472   private double[][][] globalMultipole;
473   /**
474    * Induced dipoles for each symmetry operator.
475    */
476   private double[][][] inducedDipole;
477   /**
478    * Induced dipole chain rule terms for each symmetry operator.
479    */
480   private double[][][] inducedDipoleCR;
481   /**
482    * AtomicDoubleArray implementation to use.
483    */
484   private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
485   /**
486    * Atomic Gradient array.
487    */
488   private AtomicDoubleArray3D grad;
489   /**
490    * Atomic Torque array.
491    */
492   private AtomicDoubleArray3D torque;
493   /**
494    * Atomic Born radii chain-rule array.
495    */
496   private AtomicDoubleArray bornRadiiChainRule;
497   /**
498    * Atomic GK field array.
499    */
500   private final AtomicDoubleArray3D fieldGK;
501   /**
502    * Atomic GK field chain-rule array.
503    */
504   private final AtomicDoubleArray3D fieldGKCR;
505   /**
506    * Neighbor lists for each atom and symmetry operator.
507    */
508   private int[][][] neighborLists;
509   /**
510    * This field is because re-initializing the force field resizes some arrays but not others; that
511    * second category must, when called on, be resized not to the current number of atoms but to the
512    * maximum number of atoms (and thus to the size of the first category of arrays).
513    */
514   private int maxNumAtoms;
515   /**
516    * GK cut-off distance.
517    */
518   private double cutoff;
519   /**
520    * GK cut-off distance squared.
521    */
522   private double cut2;
523   /**
524    * Boolean flag to indicate GK will be scaled by the lambda state variable.
525    */
526   private boolean lambdaTerm;
527   /**
528    * The current value of the lambda state variable.
529    */
530   private double lambda = 1.0;
531   /**
532    * lPow equals lambda^polarizationLambdaExponent, where polarizationLambdaExponent is also used by
533    * PME.
534    */
535   private double lPow = 1.0;
536   /**
537    * First derivative of lPow with respect to l.
538    */
539   private double dlPow = 0.0;
540   /**
541    * Second derivative of lPow with respect to l.
542    */
543   private double dl2Pow = 0.0;
544   /**
545    * Total Solvation Energy.
546    */
547   private double solvationEnergy = 0.0;
548   /**
549    * Electrostatic Solvation Energy.
550    */
551   private double gkEnergy = 0.0;
552   /**
553    * Electrostatic Solvation Energy from Permanent Multipoles.
554    */
555   private double gkPermanentEnergy = 0.0;
556   /**
557    * Electrostatic Solvation Energy from Induced Dipoles.
558    */
559   private double gkPolarizationEnergy = 0.0;
560   /**
561    * Dispersion Solvation Energy.
562    */
563   private double dispersionEnergy = 0.0;
564   /**
565    * Cavitation Solvation Energy.
566    */
567   private double cavitationEnergy = 0.0;
568   /**
569    * Time to compute GK electrostatics.
570    */
571   private long gkTime = 0;
572   /**
573    * Time to compute Dispersion energy.
574    */
575   private long dispersionTime = 0;
576   /**
577    * Time to compute Cavitation energy.
578    */
579   private long cavitationTime = 0;
580   /**
581    * Forces all atoms to be considered during Born radius updates.
582    */
583   private final boolean nativeEnvironmentApproximation;
584   /**
585    * Flag to turn on use of perfect Born radii.
586    */
587   private final boolean perfectRadii;
588 
589   /**
590    * Constructor for GeneralizedKirkwood.
591    *
592    * @param forceField        a {@link ffx.potential.parameters.ForceField} object.
593    * @param atoms             an array of {@link ffx.potential.bonded.Atom} objects.
594    * @param particleMeshEwald a {@link ParticleMeshEwald} object.
595    * @param crystal           a {@link ffx.crystal.Crystal} object.
596    * @param parallelTeam      a {@link edu.rit.pj.ParallelTeam} object.
597    */
598   public GeneralizedKirkwood(ForceField forceField, Atom[] atoms, ParticleMeshEwald particleMeshEwald,
599                              Crystal crystal, ParallelTeam parallelTeam, double gkCutoff) {
600     this.forceField = forceField;
601     this.atoms = atoms;
602     this.particleMeshEwald = particleMeshEwald;
603     this.crystal = crystal;
604     this.parallelTeam = parallelTeam;
605     this.electric = electric;
606     this.cutoff = gkCutoff;
607     nAtoms = atoms.length;
608     maxNumAtoms = nAtoms;
609     polarization = particleMeshEwald.polarization;
610 
611     // Read in the value of electric conversion factor.
612     electric = forceField.getDouble("ELECTRIC", Constants.DEFAULT_ELECTRIC);
613     // Set the Kirkwood multipolar reaction field constants for solvent.
614     solventDielectric = forceField.getDouble("SOLVENT_DIELECTRIC", dWater);
615     // Set the Kirkwood multipolar reaction field constants for solute.
616     soluteDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
617 
618     // Define how force arrays will be accumulated.
619     String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
620     try {
621       atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
622     } catch (Exception e) {
623       atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
624       logger.info(
625           format(" Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value,
626               atomicDoubleArrayImpl));
627     }
628 
629     String gkRadius = forceField.getString("GK_RADIUS", "SOLUTE");
630     try {
631       soluteRadiiType = SOLUTE_RADII_TYPE.valueOf(gkRadius.trim().toUpperCase());
632     } catch (Exception e) {
633       soluteRadiiType = SOLUTE_RADII_TYPE.SOLUTE;
634     }
635 
636     // Define default Bondi scale factor, and HCT overlap scale factors.
637     if (soluteRadiiType != SOLUTE_RADII_TYPE.SOLUTE) {
638       bondiScale = forceField.getDouble("SOLUTE_SCALE", DEFAULT_SOLUTE_SCALE);
639     } else {
640       // Default Bondi scale factor for Solute Radii is 1.0.
641       bondiScale = forceField.getDouble("SOLUTE_SCALE", 1.0);
642     }
643 
644     hctScale = forceField.getDouble("HCT_SCALE", DEFAULT_HCT_SCALE);
645     elementHCTScale = forceField.getBoolean("ELEMENT_HCT_SCALE", true);
646     descreenVDW = forceField.getBoolean("DESCREEN_VDW", true);
647     descreenHydrogen = forceField.getBoolean("DESCREEN_HYDROGEN", false);
648     if (descreenVDW && !descreenHydrogen) {
649       perfectHCTScale = forceField.getBoolean("PERFECT_HCT_SCALE", false);
650     } else {
651       perfectHCTScale = false;
652     }
653     descreenOffset = forceField.getDouble("DESCREEN_OFFSET", DEFAULT_DESCREEN_OFFSET);
654     // If true, the descreening integral includes the neck correction to better approximate molecular surface.
655     neckCorrection = forceField.getBoolean("NECK_CORRECTION", true);
656     double sn = forceField.getDouble("SNECK", DEFAULT_NECK_SCALE);
657     sneck = forceField.getDouble("NECK_SCALE", sn);
658     chemicallyAwareSneck = forceField.getBoolean("CHEMICALLY_AWARE_SNECK", true);
659     tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
660     double b0 = forceField.getDouble("BETA0", DEFAULT_TANH_BETA0);
661     double b1 = forceField.getDouble("BETA1", DEFAULT_TANH_BETA1);
662     double b2 = forceField.getDouble("BETA2", DEFAULT_TANH_BETA2);
663     beta0 = forceField.getDouble("TANH_BETA0", b0);
664     beta1 = forceField.getDouble("TANH_BETA1", b1);
665     beta2 = forceField.getDouble("TANH_BETA2", b2);
666 
667     BornTanhRescaling.setBeta0(beta0);
668     BornTanhRescaling.setBeta1(beta1);
669     BornTanhRescaling.setBeta2(beta2);
670 
671     // Default overlap element specific scale factors for the Hawkins, Cramer & Truhlar pairwise descreening algorithm.
672     HashMap<Integer, Double> DEFAULT_HCT_ELEMENTS = new HashMap<>();
673     // Fit default values from Corrigan et. al. interstitial spaces work
674     DEFAULT_HCT_ELEMENTS.put(1, 0.7200);
675     DEFAULT_HCT_ELEMENTS.put(6, 0.6950);
676     DEFAULT_HCT_ELEMENTS.put(7, 0.7673);
677     DEFAULT_HCT_ELEMENTS.put(8, 0.7965);
678     DEFAULT_HCT_ELEMENTS.put(15, 0.6117);
679     DEFAULT_HCT_ELEMENTS.put(16, 0.7204);
680 
681     // Add default values for all elements
682     elementHCTScaleFactors = new HashMap<>();
683     elementHCTScaleFactors.put(1, DEFAULT_HCT_ELEMENTS.get(1));
684     elementHCTScaleFactors.put(6, DEFAULT_HCT_ELEMENTS.get(6));
685     elementHCTScaleFactors.put(7, DEFAULT_HCT_ELEMENTS.get(7));
686     elementHCTScaleFactors.put(8, DEFAULT_HCT_ELEMENTS.get(8));
687     elementHCTScaleFactors.put(15, DEFAULT_HCT_ELEMENTS.get(15));
688     elementHCTScaleFactors.put(16, DEFAULT_HCT_ELEMENTS.get(16));
689 
690     // Fill elementHCTScaleFactors array based on input hct-element property flags
691     // Overwrite defaults with input values, add new values if input elements aren't represented
692     // Input lines read in as "hct-element atomicNumber hctValue"
693     // So "hct-element 1 0.7200" would set the HCT scaling factor value for hydrogen (atomic number 1) to 0.7200
694     String[] tmpElemScaleFactors = forceField.getProperties().getStringArray("hct-element");
695     for (String elemScale : tmpElemScaleFactors) {
696       String[] singleScale = elemScale.trim().split(" +");
697       if (elementHCTScaleFactors.containsKey(Integer.parseInt(singleScale[0]))) {
698         elementHCTScaleFactors.replace(Integer.parseInt(singleScale[0]), Double.parseDouble(singleScale[1]));
699       } else {
700         elementHCTScaleFactors.put(Integer.parseInt(singleScale[0]), Double.parseDouble(singleScale[1]));
701       }
702     }
703 
704     perfectRadii = forceField.getBoolean("PERFECT_RADII", false);
705 
706     // Process any radii override values.
707     String radiiProp = forceField.getString("GK_RADIIOVERRIDE", null);
708     if (radiiProp != null) {
709       String[] tokens = radiiProp.split("A");
710       for (String token : tokens) {
711         if (!token.contains("r")) {
712           logger.severe("Invalid radius override.");
713         }
714         int separator = token.indexOf("r");
715         int type = Integer.parseInt(token.substring(0, separator));
716         double factor = Double.parseDouble(token.substring(separator + 1));
717         logger.info(format(" (GK) Scaling AtomType %d with bondi factor %.2f", type, factor));
718         radiiOverride.put(type, factor);
719       }
720     }
721 
722     NonPolarModel npModel;
723     try {
724       String cavModel = forceField.getString("CAVMODEL", DEFAULT_NONPOLAR_MODEL.name());
725       cavModel = forceField.getString("NONPOLAR_MODEL", cavModel);
726       npModel = getNonPolarModel(cavModel.toUpperCase());
727     } catch (Exception ex) {
728       npModel = NonPolarModel.NONE;
729       logger.warning(format(" Error parsing non-polar model (set to NONE) %s", ex));
730     }
731     nonPolarModel = npModel;
732 
733     int threadCount = parallelTeam.getThreadCount();
734     fieldGK = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
735     fieldGKCR = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
736 
737     nativeEnvironmentApproximation =
738         forceField.getBoolean("NATIVE_ENVIRONMENT_APPROXIMATION", false);
739     probe = forceField.getDouble("PROBE_RADIUS", 1.4);
740     cut2 = cutoff * cutoff;
741     lambdaTerm = forceField.getBoolean("GK_LAMBDATERM", forceField.getBoolean("LAMBDATERM", false));
742 
743     // If polarization lambda exponent is set to 0.0, then we're running
744     // Dual-Topology and the GK energy will be scaled with the overall system lambda value.
745     double polLambdaExp = forceField.getDouble("POLARIZATION_LAMBDA_EXPONENT", 3.0);
746     if (polLambdaExp == 0.0) {
747       lambdaTerm = false;
748       logger.info(" GK lambda term set to false.");
749     }
750 
751     // If PME includes polarization and is a function of lambda, GK must also.
752     if (!lambdaTerm
753         && particleMeshEwald.getPolarizationType() != Polarization.NONE) {
754       if (forceField.getBoolean("ELEC_LAMBDATERM",
755           forceField.getBoolean("LAMBDATERM", false))) {
756         logger.info(" If PME includes polarization and is a function of lambda, GK must also.");
757         lambdaTerm = true;
758       }
759     }
760 
761     initAtomArrays();
762 
763     double tensionDefault;
764     switch (nonPolarModel) {
765       case CAV:
766         tensionDefault = DEFAULT_CAVONLY_SURFACE_TENSION;
767         surfaceAreaRegion = new SurfaceAreaRegion(atoms, x, y, z, use,
768             neighborLists, grad, threadCount, probe, tensionDefault);
769         dispersionRegion = null;
770         chandlerCavitation = null;
771         break;
772       case CAV_DISP:
773         tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
774         surfaceAreaRegion = new SurfaceAreaRegion(atoms, x, y, z, use,
775             neighborLists, grad, threadCount, probe, tensionDefault);
776         dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
777         chandlerCavitation = null;
778         break;
779       case DISP:
780         tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
781         dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
782         surfaceAreaRegion = null;
783         chandlerCavitation = null;
784         break;
785       case SEV_DISP:
786         tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
787         double[] radii = new double[nAtoms];
788         int index = 0;
789         for (Atom atom : atoms) {
790           radii[index] = atom.getVDWType().radius / 2.0;
791           index++;
792         }
793         ConnollyRegion connollyRegion = new ConnollyRegion(atoms, radii, threadCount);
794         // connollyRegion.setProbe(probe);
795         // connollyRegion.setExclude(0.0);
796         double wiggle = forceField.getDouble("WIGGLE", ConnollyRegion.DEFAULT_WIGGLE);
797         connollyRegion.setWiggle(wiggle);
798         chandlerCavitation = new ChandlerCavitation(atoms, connollyRegion, forceField);
799         dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
800         surfaceAreaRegion = null;
801         break;
802       case GAUSS_DISP:
803         dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
804         surfaceAreaRegion = null;
805         tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
806         GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
807         chandlerCavitation = new ChandlerCavitation(atoms, gaussVol, forceField);
808         break;
809       case BORN_CAV_DISP:
810         tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
811         surfaceAreaRegion = null;
812         chandlerCavitation = null;
813         dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
814         break;
815       case HYDROPHOBIC_PMF:
816       case BORN_SOLV:
817       case NONE:
818       default:
819         tensionDefault = DEFAULT_CAVONLY_SURFACE_TENSION;
820         surfaceAreaRegion = null;
821         dispersionRegion = null;
822         chandlerCavitation = null;
823         break;
824     }
825 
826     gkc = forceField.getDouble("GKC", DEFAULT_GKC);
827 
828     surfaceTension = forceField.getDouble("SURFACE_TENSION", tensionDefault);
829     solventPressue = forceField.getDouble("SOLVENT_PRESSURE", DEFAULT_SOLVENT_PRESSURE);
830     // Volume to surface area cross-over point (A).
831     double crossOver = forceField.getDouble("CROSS_OVER", DEFAULT_CROSSOVER);
832     if (chandlerCavitation != null) {
833       // Set the cross-over first.
834       chandlerCavitation.setCrossOver(crossOver);
835       // Surface tension and solvent pressure will over-write cross-over if it's not appropriate.
836       chandlerCavitation.setSurfaceTension(surfaceTension);
837       chandlerCavitation.setSolventPressure(solventPressue);
838       crossOver = chandlerCavitation.getCrossOver();
839     }
840     if (surfaceAreaRegion != null) {
841       surfaceAreaRegion.setSurfaceTension(surfaceTension);
842     }
843     double dispersionOffset = forceField.getDouble("DISPERSION_OFFSET", DEFAULT_DISPERSION_OFFSET);
844     if (dispersionRegion != null) {
845       dispersionRegion.setDispersionOffest(dispersionOffset);
846     }
847 
848     initializationRegion = new InitializationRegion(threadCount);
849     bornRadiiRegion = new BornRadiiRegion(threadCount, nAtoms, forceField, neckCorrection,
850         tanhCorrection, perfectHCTScale);
851     permanentGKFieldRegion = new PermanentGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
852     inducedGKFieldRegion = new InducedGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
853     if (!perfectRadii) {
854       bornGradRegion = new BornGradRegion(threadCount, neckCorrection, tanhCorrection, perfectHCTScale);
855     } else {
856       // No Born chain-rule terms when using Perfect Born Radii.
857       bornGradRegion = null;
858     }
859     boolean gkQI = forceField.getBoolean("GK_QI", false);
860     gkEnergyRegion = new GKEnergyRegion(threadCount, polarization, nonPolarModel, surfaceTension,
861         probe, electric, soluteDielectric, solventDielectric, gkc, gkQI);
862 
863     if (gkQI) {
864       logger.info("  Continuum Solvation (QI)");
865     } else {
866       logger.info("  Continuum Solvation");
867     }
868     logger.info(format("   Radii:                              %8s", soluteRadiiType));
869     if (cutoff != Double.POSITIVE_INFINITY) {
870       logger.info(format("   Generalized Kirkwood Cut-Off:       %8.4f (A)", cutoff));
871     } else {
872       logger.info("   Generalized Kirkwood Cut-Off:           NONE");
873     }
874     logger.info(format("   GKC:                                %8.4f", gkc));
875     logger.info(format("   Solvent Dielectric:                 %8.4f", solventDielectric));
876     logger.info(format("   Solute Dielectric:                  %8.4f", soluteDielectric));
877 
878     if (perfectRadii) {
879       logger.info(format("   Use Perfect Born Radii:             %8B", perfectRadii));
880     } else {
881       logger.info(format("   Descreen with vdW Radii:            %8B", descreenVDW));
882       logger.info(format("   Descreen with Hydrogen Atoms:       %8B", descreenHydrogen));
883       logger.info(format("   Descreen Offset:                    %8.4f", descreenOffset));
884       if (neckCorrection) {
885         logger.info(format("   Use Neck Correction:                %8B", neckCorrection));
886         logger.info(format("   Sneck Scale Factor:                 %8.4f", sneck));
887         logger.info(format("   Chemically Aware Sneck:             %8B", chemicallyAwareSneck));
888       }
889       if (tanhCorrection) {
890         logger.info(format("   Use Tanh Correction                 %8B", tanhCorrection));
891         logger.info(format("    Beta0:                             %8.4f", beta0));
892         logger.info(format("    Beta1:                             %8.4f", beta1));
893         logger.info(format("    Beta2:                             %8.4f", beta2));
894       }
895       if (perfectHCTScale) {
896         logger.info(format("   GaussVol HCT Scale Factors:         %8B", perfectHCTScale));
897       }
898       logger.info(format("   General HCT Scale Factor:           %8.4f",
899           forceField.getDouble("HCT-SCALE", DEFAULT_HCT_SCALE)));
900       if (elementHCTScale) {
901         logger.info(format("   Element-Specific HCT Scale Factors: %8B", elementHCTScale));
902         Integer[] elementHCTkeyset = elementHCTScaleFactors.keySet().toArray(new Integer[0]);
903         Arrays.sort(elementHCTkeyset);
904         for (Integer key : elementHCTkeyset) {
905           logger.info(format("    HCT-Element # %d:                   %8.4f", key, elementHCTScaleFactors.get(key)));
906         }
907       }
908     }
909 
910     logger.info(
911         format("   Non-Polar Model:                  %10s",
912             nonPolarModel.toString().replace('_', '-')));
913 
914     if (nonPolarModel.equals(NonPolarModel.GAUSS_DISP)) {
915       logger.info(format("    GaussVol Radii Offset:               %2.4f",
916           forceField.getDouble("GAUSSVOL_RADII_OFFSET", 0.0)));
917       logger.info(format("    GaussVol Radii Scale:                %2.4f",
918           forceField.getDouble("GAUSSVOL_RADII_SCALE", 1.15)));
919     }
920 
921     if (dispersionRegion != null) {
922       logger.info(format("   Dispersion Integral Offset:         %8.4f (A)",
923           dispersionRegion.getDispersionOffset()));
924     }
925 
926     if (surfaceAreaRegion != null) {
927       logger.info(format("   Cavitation Probe Radius:            %8.4f (A)", probe));
928       logger.info(
929           format("   Cavitation Surface Tension:         %8.4f (Kcal/mol/A^2)", surfaceTension));
930     } else if (chandlerCavitation != null) {
931       logger.info(
932           format("   Cavitation Solvent Pressure:        %8.4f (Kcal/mol/A^3)", solventPressue));
933       logger.info(
934           format("   Cavitation Surface Tension:         %8.4f (Kcal/mol/A^2)", surfaceTension));
935       logger.info(format("   Cavitation Cross-Over Radius:       %8.4f (A)", crossOver));
936     }
937 
938     // Print out all Base Radii
939     if (logger.isLoggable(Level.FINE)) {
940       logger.fine("      Base Radii  Descreen Radius  Overlap Scale  Neck Scale");
941       for (int i = 0; i < nAtoms; i++) {
942         logger.info(
943             format("   %s %8.6f %8.6f %5.3f %5.3f",
944                 atoms[i].toString(), baseRadius[i], descreenRadius[i], overlapScale[i],
945                 neckScale[i]));
946       }
947     }
948   }
949 
950   /**
951    * GK is using perfect radii where available.
952    *
953    * @return True if using perfect radii.
954    */
955   public boolean getUsePerfectRadii() {
956     return perfectRadii;
957   }
958 
959   /**
960    * Return perfect Born radii read in as keywords, or base radii if perfect radii are not
961    * available.
962    *
963    * @return Array of perfect Born radii.
964    */
965   public double[] getPerfectRadii() {
966     bornRadiiRegion.init(atoms, crystal, sXYZ, neighborLists, baseRadius, descreenRadius,
967         overlapScale, neckScale, descreenOffset, use, cut2, nativeEnvironmentApproximation, born);
968     return bornRadiiRegion.getPerfectRadii();
969   }
970 
971   /**
972    * getNonPolarModel.
973    *
974    * @param nonpolarModel a {@link java.lang.String} object.
975    * @return a {@link NonPolarModel} object.
976    */
977   public static NonPolarModel getNonPolarModel(String nonpolarModel) {
978     try {
979       return NonPolarModel.valueOf(toEnumForm(nonpolarModel));
980     } catch (IllegalArgumentException ex) {
981       logger.warning(" Unrecognized nonpolar model requested; defaulting to NONE.");
982       return NonPolarModel.NONE;
983     }
984   }
985 
986   /**
987    * computeBornRadii
988    */
989   public void computeBornRadii() {
990     // The solute radii can change during titration based rotamer optimization.
991     applySoluteRadii();
992 
993     // Update tanh correction parameters.
994     if (tanhCorrection) {
995       BornTanhRescaling.setBeta0(beta0);
996       BornTanhRescaling.setBeta1(beta1);
997       BornTanhRescaling.setBeta2(beta2);
998     }
999 
1000     try {
1001       bornRadiiRegion.init(
1002           atoms,
1003           crystal,
1004           sXYZ,
1005           neighborLists,
1006           baseRadius,
1007           descreenRadius,
1008           overlapScale,
1009           neckScale,
1010           descreenOffset,
1011           use,
1012           cut2,
1013           nativeEnvironmentApproximation,
1014           born);
1015       parallelTeam.execute(bornRadiiRegion);
1016     } catch (Exception e) {
1017       String message = "Fatal exception computing Born radii.";
1018       logger.log(Level.SEVERE, message, e);
1019     }
1020   }
1021 
1022   /**
1023    * computeInducedGKField
1024    */
1025   public void computeInducedGKField() {
1026     try {
1027       fieldGK.reset(parallelTeam);
1028       fieldGKCR.reset(parallelTeam);
1029       inducedGKFieldRegion.init(atoms, inducedDipole, inducedDipoleCR, crystal, sXYZ,
1030           neighborLists, use, cut2, born, fieldGK, fieldGKCR);
1031       parallelTeam.execute(inducedGKFieldRegion);
1032       fieldGK.reduce(parallelTeam);
1033       fieldGKCR.reduce(parallelTeam);
1034     } catch (Exception e) {
1035       String message = "Fatal exception computing induced GK field.";
1036       logger.log(Level.SEVERE, message, e);
1037     }
1038   }
1039 
1040   /**
1041    * computePermanentGKField
1042    */
1043   public void computePermanentGKField() {
1044     try {
1045       fieldGK.reset(parallelTeam);
1046       permanentGKFieldRegion.init(
1047           atoms, globalMultipole, crystal, sXYZ, neighborLists, use, cut2, born, fieldGK);
1048       parallelTeam.execute(permanentGKFieldRegion);
1049       fieldGK.reduce(parallelTeam);
1050     } catch (Exception e) {
1051       String message = "Fatal exception computing permanent GK field.";
1052       logger.log(Level.SEVERE, message, e);
1053     }
1054   }
1055 
1056   /**
1057    * Returns the base radii used for the GK calculation.
1058    *
1059    * <p>These are the radii used to compute the Born radii, and are not the same as the Born radii.
1060    *
1061    * @return Base radii for GK calculations.
1062    */
1063   public double[] getBaseRadii() {
1064     return baseRadius;
1065   }
1066 
1067   /**
1068    * Returns the descreening radii used for the GK calculation.
1069    *
1070    * <p>These are the radii used to during pairwise descreening.
1071    *
1072    * @return Descreening radii for GK calculations.
1073    */
1074   public double[] getDescreenRadii() {
1075     return descreenRadius;
1076   }
1077 
1078   /**
1079    * Returns the cavitation component of the solvation energy.
1080    *
1081    * @return Cavitation energy
1082    */
1083   public double getCavitationEnergy() {
1084     return cavitationEnergy;
1085   }
1086 
1087   /**
1088    * Return the Chandler Cavitation instance.
1089    *
1090    * @return ChandlerCavitation instance.
1091    */
1092   public ChandlerCavitation getChandlerCavitation() {
1093     return chandlerCavitation;
1094   }
1095 
1096   /**
1097    * Getter for the field <code>cutoff</code>.
1098    *
1099    * @return a double.
1100    */
1101   public double getCutoff() {
1102     return cutoff;
1103   }
1104 
1105   /**
1106    * Setter for the field <code>cutoff</code>.
1107    *
1108    * @param cutoff a double.
1109    */
1110   public void setCutoff(double cutoff) {
1111     this.cutoff = cutoff;
1112     this.cut2 = cutoff * cutoff;
1113   }
1114 
1115   /**
1116    * Returns the dielectric offset (in Angstroms).
1117    *
1118    * @return Currently: 0.09 Angstroms.
1119    */
1120   public double getDielecOffset() {
1121     return dOffset;
1122   }
1123 
1124   /**
1125    * Return the descreening dielectric offset.
1126    *
1127    * @return The offset (A).
1128    */
1129   public double getDescreenOffset() {
1130     return descreenOffset;
1131   }
1132 
1133   /**
1134    * Returns the dispersion component of the solvation energy.
1135    *
1136    * @return Dispersion energy
1137    */
1138   public double getDispersionEnergy() {
1139     return dispersionEnergy;
1140   }
1141 
1142   public DispersionRegion getDispersionRegion() {
1143     return dispersionRegion;
1144   }
1145 
1146   public AtomicDoubleArray3D getFieldGK() {
1147     return fieldGK;
1148   }
1149 
1150   public AtomicDoubleArray3D getFieldGKCR() {
1151     return fieldGKCR;
1152   }
1153 
1154   public SurfaceAreaRegion getSurfaceAreaRegion() {
1155     return surfaceAreaRegion;
1156   }
1157 
1158   /**
1159    * Returns the GK component of the solvation energy.
1160    *
1161    * @return GK electrostatic energy
1162    */
1163   public double getGeneralizedKirkwoordEnergy() {
1164     return gkEnergy;
1165   }
1166 
1167   /**
1168    * Returns the GK component of the solvation energy.
1169    *
1170    * @return GK electrostatic energy
1171    */
1172   public double getGeneralizedKirkwoordPermanentEnergy() {
1173     return gkPermanentEnergy;
1174   }
1175 
1176   /**
1177    * Returns the GK component of the solvation energy.
1178    *
1179    * @return GK electrostatic energy
1180    */
1181   public double getGeneralizedKirkwoordPolariztionEnergy() {
1182     return gkPolarizationEnergy;
1183   }
1184 
1185   public AtomicDoubleArray3D getGrad() {
1186     return grad;
1187   }
1188 
1189   /**
1190    * getInteractions
1191    *
1192    * @return a int.
1193    */
1194   public int getInteractions() {
1195     return gkEnergyRegion.getInteractions();
1196   }
1197 
1198   /**
1199    * {@inheritDoc}
1200    */
1201   @Override
1202   public double getLambda() {
1203     return lambda;
1204   }
1205 
1206   /**
1207    * {@inheritDoc}
1208    *
1209    * <p>Updates the value of lPow.
1210    */
1211   @Override
1212   public void setLambda(double lambda) {
1213     if (lambdaTerm) {
1214       this.lambda = lambda;
1215     } else {
1216       // If the lambdaTerm flag is false, lambda is set to one.
1217       this.lambda = 1.0;
1218       lPow = 1.0;
1219       dlPow = 0.0;
1220       dl2Pow = 0.0;
1221     }
1222   }
1223 
1224   /**
1225    * Checks whether GK uses the Native Environment Approximation.
1226    *
1227    * <p>This (previously known as born-use-all) is useful for rotamer optimization under continuum
1228    * solvent. If a large number of sidechain atoms are completely removed from the GK/GB calculation,
1229    * the remaining sidechains are overly solvated. The NEA says "we will keep all sidechains not
1230    * under optimization in some default state and let them contribute to Born radii calculations, but
1231    * still exclude their solvation energy components."
1232    *
1233    * @return Whether the NEA is in use.
1234    */
1235   public boolean getNativeEnvironmentApproximation() {
1236     return nativeEnvironmentApproximation;
1237   }
1238 
1239   /**
1240    * getNonPolarModel.
1241    *
1242    * @return a {@link NonPolarModel} object.
1243    */
1244   public NonPolarModel getNonPolarModel() {
1245     return nonPolarModel;
1246   }
1247 
1248   /**
1249    * Returns the overlap scale factors used for the GK calculation.
1250    *
1251    * @return Overlap scale factors for GK calculations.
1252    */
1253   public double[] getOverlapScale() {
1254     return overlapScale;
1255   }
1256 
1257   /**
1258    * Returns the neck scale factors used for the GK calculation.
1259    *
1260    * @return Neck scale factors for GK calculations.
1261    */
1262   public double[] getNeckScale() {
1263     return neckScale;
1264   }
1265 
1266   /**
1267    * Returns the tanh correction factor.
1268    *
1269    * <p>When true, the tanh correction is applied to the Born radii.
1270    *
1271    * @return True if tanh correction is applied.
1272    */
1273   public boolean getTanhCorrection() {
1274     return tanhCorrection;
1275   }
1276 
1277   /**
1278    * Returns the probe radius (typically 1.4 Angstroms).
1279    *
1280    * @return Radius of the solvent probe.
1281    */
1282   public double getProbeRadius() {
1283     return probe;
1284   }
1285 
1286   /**
1287    * Returns the solvent relative permittivity (typically 78.3).
1288    *
1289    * @return Relative permittivity of the solvent.
1290    */
1291   public double getSolventPermittivity() {
1292     return solventDielectric;
1293   }
1294 
1295   /**
1296    * Returns the solvent relative permittivity (typically 1.0).
1297    *
1298    * @return Relative permittivity of the solute.
1299    */
1300   public double getSolutePermittivity() {
1301     return soluteDielectric;
1302   }
1303 
1304   /**
1305    * Getter for the field <code>surfaceTension</code>.
1306    *
1307    * @return a double.
1308    */
1309   public double getSurfaceTension() {
1310     return surfaceTension;
1311   }
1312 
1313   public AtomicDoubleArray3D getTorque() {
1314     return torque;
1315   }
1316 
1317   /**
1318    * {@inheritDoc}
1319    *
1320    * <p>The 2nd derivative is 0.0. (U=Lambda*Egk, dU/dL=Egk, d2U/dL2=0.0)
1321    */
1322   @Override
1323   public double getd2EdL2() {
1324     if (lambdaTerm) {
1325       return dl2Pow * solvationEnergy;
1326     } else {
1327       return 0.0;
1328     }
1329   }
1330 
1331   /**
1332    * {@inheritDoc}
1333    */
1334   @Override
1335   public double getdEdL() {
1336     if (lambdaTerm) {
1337       return dlPow * solvationEnergy;
1338     }
1339     return 0.0;
1340   }
1341 
1342   /**
1343    * {@inheritDoc}
1344    *
1345    * <p>These contributions are already aggregated into the arrays used by PME.
1346    */
1347   @Override
1348   public void getdEdXdL(double[] gradient) {
1349     // This contribution is collected by GeneralizedKirkwood.reduce
1350   }
1351 
1352   public void init() {
1353     initializationRegion.init(this, atoms, lambdaTerm, grad, torque, bornRadiiChainRule);
1354     initializationRegion.executeWith(parallelTeam);
1355   }
1356 
1357   public void reduce(
1358       AtomicDoubleArray3D g,
1359       AtomicDoubleArray3D t,
1360       AtomicDoubleArray3D lg,
1361       AtomicDoubleArray3D lt) {
1362     grad.reduce(parallelTeam);
1363     torque.reduce(parallelTeam);
1364     for (int i = 0; i < nAtoms; i++) {
1365       g.add(0, i, lPow * grad.getX(i), lPow * grad.getY(i), lPow * grad.getZ(i));
1366       t.add(0, i, lPow * torque.getX(i), lPow * torque.getY(i), lPow * torque.getZ(i));
1367       if (lambdaTerm) {
1368         lg.add(0, i, dlPow * grad.getX(i), dlPow * grad.getY(i), dlPow * grad.getZ(i));
1369         lt.add(0, i, dlPow * torque.getX(i), dlPow * torque.getY(i), dlPow * torque.getZ(i));
1370       }
1371     }
1372   }
1373 
1374   public AtomicDoubleArray getSelfEnergy() {
1375     // Init and execute gkEnergyRegion (?)
1376     return gkEnergyRegion.getSelfEnergy();
1377   }
1378 
1379   public double[] getBorn() {
1380     return bornRadiiRegion.getBorn();
1381   }
1382 
1383   /**
1384    * Setter for element-specific HCT overlap scale factors
1385    *
1386    * @param elementHCT HashMap containing element name keys and scale factor values
1387    */
1388   public void setElementHCTScaleFactors(HashMap<Integer, Double> elementHCT) {
1389     for (Integer element : elementHCT.keySet()) {
1390       if (elementHCTScaleFactors.containsKey(element)) {
1391         elementHCTScaleFactors.replace(element, elementHCT.get(element));
1392       } else {
1393         logger.info("No HCT scale factor set for element: " + element);
1394       }
1395     }
1396     initAtomArrays();
1397   }
1398 
1399   /**
1400    * Setter for the field <code>atoms</code>.
1401    *
1402    * @param atoms an array of {@link ffx.potential.bonded.Atom} objects.
1403    */
1404   public void setAtoms(Atom[] atoms) {
1405     this.atoms = atoms;
1406     nAtoms = atoms.length;
1407     maxNumAtoms = max(nAtoms, maxNumAtoms);
1408     initAtomArrays();
1409   }
1410 
1411   /**
1412    * Setter for the field <code>crystal</code>.
1413    *
1414    * @param crystal a {@link ffx.crystal.Crystal} object.
1415    */
1416   public void setCrystal(Crystal crystal) {
1417     this.crystal = crystal;
1418   }
1419 
1420   /**
1421    * Setter for the field <code>neighborLists</code>.
1422    *
1423    * @param neighbors The neighbor list for the GK calculation.
1424    */
1425   public void setNeighborList(int[][][] neighbors) {
1426     this.neighborLists = neighbors;
1427   }
1428 
1429   /**
1430    * Setter for the field <code>use</code>.
1431    *
1432    * @param use the use array indicating which atoms are used in the GK calculation.
1433    */
1434   public void setUse(boolean[] use) {
1435     this.use = use;
1436   }
1437 
1438   /**
1439    * solvationEnergy
1440    *
1441    * @param gradient a boolean.
1442    * @param print    a boolean.
1443    * @return a double.
1444    */
1445   public double solvationEnergy(boolean gradient, boolean print) {
1446     return solvationEnergy(0.0, gradient, print);
1447   }
1448 
1449   /**
1450    * solvationEnergy
1451    *
1452    * @param gkInducedCorrectionEnergy GK vacuum to SCRF polarization energy cost.
1453    * @param gradient                  a boolean.
1454    * @param print                     a boolean.
1455    * @return a double.
1456    */
1457   public double solvationEnergy(double gkInducedCorrectionEnergy, boolean gradient, boolean print) {
1458     cavitationEnergy = 0.0;
1459     dispersionEnergy = 0.0;
1460     gkEnergy = 0.0;
1461 
1462     try {
1463       // Find the GK energy.
1464       gkTime = -System.nanoTime();
1465       gkEnergyRegion.init(
1466           atoms,
1467           globalMultipole,
1468           inducedDipole,
1469           inducedDipoleCR,
1470           crystal,
1471           sXYZ,
1472           neighborLists,
1473           use,
1474           cut2,
1475           baseRadius,
1476           born,
1477           gradient,
1478           parallelTeam,
1479           grad,
1480           torque,
1481           bornRadiiChainRule);
1482       parallelTeam.execute(gkEnergyRegion);
1483       gkTime += System.nanoTime();
1484 
1485       // Find the nonpolar energy.
1486       switch (nonPolarModel) {
1487         case CAV:
1488           cavitationTime = -System.nanoTime();
1489           parallelTeam.execute(surfaceAreaRegion);
1490           cavitationEnergy = surfaceAreaRegion.getEnergy();
1491           cavitationTime += System.nanoTime();
1492           break;
1493         case CAV_DISP:
1494           dispersionTime = -System.nanoTime();
1495           dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1496           parallelTeam.execute(dispersionRegion);
1497           dispersionEnergy = dispersionRegion.getEnergy();
1498           dispersionTime += System.nanoTime();
1499           cavitationTime = -System.nanoTime();
1500           parallelTeam.execute(surfaceAreaRegion);
1501           cavitationEnergy = surfaceAreaRegion.getEnergy();
1502           cavitationTime += System.nanoTime();
1503           break;
1504         case DISP:
1505           dispersionTime = -System.nanoTime();
1506           dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1507           parallelTeam.execute(dispersionRegion);
1508           dispersionEnergy = dispersionRegion.getEnergy();
1509           dispersionTime += System.nanoTime();
1510           break;
1511         case SEV_DISP:
1512           dispersionTime = -System.nanoTime();
1513           dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1514           parallelTeam.execute(dispersionRegion);
1515           dispersionEnergy = dispersionRegion.getEnergy();
1516           dispersionTime += System.nanoTime();
1517           cavitationTime = -System.nanoTime();
1518           cavitationTime = -System.nanoTime();
1519           double[][] positions = new double[nAtoms][3];
1520           for (int i = 0; i < nAtoms; i++) {
1521             positions[i][0] = atoms[i].getX();
1522             positions[i][1] = atoms[i].getY();
1523             positions[i][2] = atoms[i].getZ();
1524           }
1525           chandlerCavitation.energyAndGradient(positions, grad);
1526           cavitationEnergy = chandlerCavitation.getEnergy();
1527           cavitationTime += System.nanoTime();
1528           break;
1529         case GAUSS_DISP:
1530           dispersionTime = -System.nanoTime();
1531           dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1532           parallelTeam.execute(dispersionRegion);
1533           dispersionEnergy = dispersionRegion.getEnergy();
1534           dispersionTime += System.nanoTime();
1535           cavitationTime = -System.nanoTime();
1536           positions = new double[nAtoms][3];
1537           for (int i = 0; i < nAtoms; i++) {
1538             positions[i][0] = atoms[i].getX();
1539             positions[i][1] = atoms[i].getY();
1540             positions[i][2] = atoms[i].getZ();
1541           }
1542           chandlerCavitation.energyAndGradient(positions, grad);
1543           cavitationEnergy = chandlerCavitation.getEnergy();
1544           cavitationTime += System.nanoTime();
1545           break;
1546         case BORN_CAV_DISP:
1547           dispersionTime = -System.nanoTime();
1548           dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1549           parallelTeam.execute(dispersionRegion);
1550           dispersionEnergy = dispersionRegion.getEnergy();
1551           dispersionTime += System.nanoTime();
1552           break;
1553         case HYDROPHOBIC_PMF:
1554         case BORN_SOLV:
1555         case NONE:
1556         default:
1557           break;
1558       }
1559     } catch (Exception e) {
1560       String message = "Fatal exception computing the continuum solvation energy.";
1561       logger.log(Level.SEVERE, message, e);
1562     }
1563 
1564     // Compute the Born radii chain rule term.
1565     if (gradient && !perfectRadii) {
1566       try {
1567         gkTime -= System.nanoTime();
1568         bornGradRegion.init(
1569             atoms, crystal, sXYZ, neighborLists,
1570             baseRadius, descreenRadius, overlapScale, neckScale, descreenOffset,
1571             bornRadiiRegion.getUnscaledBornIntegral(), use, cut2,
1572             nativeEnvironmentApproximation, born, grad, bornRadiiChainRule);
1573         bornGradRegion.executeWith(parallelTeam);
1574         gkTime += System.nanoTime();
1575       } catch (Exception e) {
1576         String message = "Fatal exception computing Born radii chain rule term.";
1577         logger.log(Level.SEVERE, message, e);
1578       }
1579     }
1580 
1581     gkEnergy = gkEnergyRegion.getEnergy() + gkInducedCorrectionEnergy;
1582     gkPermanentEnergy = gkEnergyRegion.getPermanentEnergy();
1583     gkPolarizationEnergy = gkEnergy - gkPermanentEnergy;
1584 
1585     // The following expression is equivalent to the former.
1586     // gkPolarizationEnergy = gkEnergyRegion.getPolarizationEnergy() + gkInducedCorrectionEnergy;
1587 
1588     // Solvation energy is the sum of cavitation, dispersion and GK
1589     solvationEnergy = cavitationEnergy + dispersionEnergy + gkEnergy;
1590 
1591     if (print) {
1592       logger.info(format(" Generalized Kirkwood%16.8f %10.3f", gkEnergy, gkTime * 1e-9));
1593       switch (nonPolarModel) {
1594         case CAV:
1595           logger.info(format(" Cavitation          %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1596           break;
1597         case DISP:
1598           logger.info(format(" Dispersion          %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1599           break;
1600         case CAV_DISP:
1601         case SEV_DISP:
1602         case GAUSS_DISP:
1603           logger.info(format(" Cavitation          %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1604           logger.info(format(" Dispersion          %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1605           break;
1606         case BORN_CAV_DISP:
1607           logger.info(format(" Dispersion          %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1608           break;
1609         case HYDROPHOBIC_PMF:
1610         case BORN_SOLV:
1611         case NONE:
1612         default:
1613           break;
1614       }
1615     }
1616 
1617     if (lambdaTerm) {
1618       return lPow * solvationEnergy;
1619     } else {
1620       return solvationEnergy;
1621     }
1622   }
1623 
1624   private void initAtomArrays() {
1625     sXYZ = particleMeshEwald.coordinates;
1626 
1627     x = sXYZ[0][0];
1628     y = sXYZ[0][1];
1629     z = sXYZ[0][2];
1630 
1631     globalMultipole = particleMeshEwald.globalMultipole;
1632     inducedDipole = particleMeshEwald.inducedDipole;
1633     inducedDipoleCR = particleMeshEwald.inducedDipoleCR;
1634     neighborLists = particleMeshEwald.neighborLists;
1635 
1636     if (grad == null) {
1637       int threadCount = parallelTeam.getThreadCount();
1638       grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1639       torque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1640       bornRadiiChainRule = atomicDoubleArrayImpl.createInstance(threadCount, nAtoms);
1641     } else {
1642       grad.alloc(nAtoms);
1643       torque.alloc(nAtoms);
1644       bornRadiiChainRule.alloc(nAtoms);
1645     }
1646 
1647     fieldGK.alloc(nAtoms);
1648     fieldGKCR.alloc(nAtoms);
1649 
1650     if (baseRadius == null || baseRadius.length < nAtoms) {
1651       baseRadius = new double[nAtoms];
1652       descreenRadius = new double[nAtoms];
1653       overlapScale = new double[nAtoms];
1654       neckScale = new double[nAtoms];
1655       born = new double[nAtoms];
1656       use = new boolean[nAtoms];
1657     }
1658 
1659     fill(use, true);
1660 
1661     setSoluteRadii(forceField, atoms, soluteRadiiType);
1662     applySoluteRadii();
1663 
1664     if (dispersionRegion != null) {
1665       dispersionRegion.allocate(atoms);
1666     }
1667 
1668     if (surfaceAreaRegion != null) {
1669       surfaceAreaRegion.init();
1670     }
1671 
1672     if (chandlerCavitation != null) {
1673       GaussVol gaussVol = chandlerCavitation.getGaussVol();
1674       try {
1675         gaussVol.allocate(atoms);
1676       } catch (Exception e) {
1677         logger.severe(" " + e);
1678       }
1679     }
1680   }
1681 
1682   /**
1683    * Update GK solute parameters for a given atom. This should only be called after each atom is
1684    * assigned a "SoluteType".
1685    *
1686    * @param i The atom to update.
1687    */
1688   public void udpateSoluteParameters(int i) {
1689     Atom atom = atoms[i];
1690     SoluteType soluteType = atom.getSoluteType();
1691     if (soluteType == null) {
1692       logger.severe(format(" No SoluteType for atom %s.", atom));
1693       return;
1694     }
1695 
1696     // Assign the base radius.
1697     baseRadius[i] = soluteType.gkDiameter * 0.5 * bondiScale;
1698     // Assign a default overlap scale factor.
1699     overlapScale[i] = hctScale;
1700     // Use element specific HCT scaling factors.
1701     if (elementHCTScale) {
1702       int atomicNumber = atom.getAtomicNumber();
1703       if (elementHCTScaleFactors.containsKey(atomicNumber)) {
1704         overlapScale[i] = elementHCTScaleFactors.get(atomicNumber);
1705       } else {
1706         logger.fine(format(" No element-specific HCT scale factor for %s", atom));
1707         overlapScale[i] = hctScale;
1708       }
1709     }
1710     // Assign the default descreen radius to equal the base radius.
1711     descreenRadius[i] = baseRadius[i];
1712 
1713     // Handle radii override values.
1714     AtomType atomType = atom.getAtomType();
1715     if (radiiOverride.containsKey(atomType.type)) {
1716       double override = radiiOverride.get(atomType.type);
1717       // Remove default bondiFactor, and apply override.
1718       baseRadius[i] = baseRadius[i] * override / bondiScale;
1719       logger.fine(format(
1720           " Scaling %s (atom type %d) to %7.4f (Bondi factor %7.4f)",
1721           atom, atomType.type, baseRadius[i], override));
1722       descreenRadius[i] = baseRadius[i];
1723     }
1724 
1725     // Apply the descreenWithVDW flag.
1726     if (descreenVDW) {
1727       descreenRadius[i] = atom.getVDWType().radius / 2.0;
1728     }
1729 
1730     // Apply the descreenWithHydrogen flag.
1731     if (!descreenHydrogen && atom.getAtomicNumber() == 1) {
1732       overlapScale[i] = 0.0;
1733     }
1734 
1735     // Set Sneck scaling parameters based on atom type and number of bound non-hydrogen atoms.
1736     // The Sneck values for hydrogen atoms are controlled by their heavy atom.
1737     if (!atom.isHydrogen()) {
1738       // If the overlap scale factor is zero, then so is the neck overlap.
1739       if (!neckCorrection || overlapScale[i] == 0.0) {
1740         neckScale[i] = 0.0;
1741       } else {
1742         if (chemicallyAwareSneck) {
1743           // Determine number of bound heavy atoms for each non-hydrogen atom if chemically aware Sneck is being used
1744           int numBoundHeavyAtoms = 0;
1745           for (Atom boundAtom : atom.get12List()) {
1746             if (!boundAtom.isHydrogen()) {
1747               numBoundHeavyAtoms++;
1748             }
1749           }
1750           // Use this number to determine Sneck scaling parameter
1751           if (numBoundHeavyAtoms == 0) {
1752             // Sneck for lone ions or molecules like methane, which are not descreened by any other atoms
1753             neckScale[i] = 1.0;
1754           } else {
1755             neckScale[i] = atom.getSoluteType().sneck * (5.0 - numBoundHeavyAtoms) / 4.0;
1756           }
1757         } else {
1758           // Non-chemically aware Sneck - set neckScale to the max (input) Sneck value for all non-hydrogen atoms
1759           neckScale[i] = atom.getSoluteType().sneck;
1760         }
1761       }
1762 
1763       // Set hydrogen atom neck scaling factors to match those of their heavy atom binding partner.
1764       // By default, hydrogen atoms don't descreen other atoms. However, when they're descreened,
1765       // their contribution to the mixed neck value should matches their heavy atom.
1766       for (Atom boundAtom : atom.get12List()) {
1767         if (boundAtom.isHydrogen()) {
1768           int hydrogenIndex = boundAtom.getIndex();
1769           neckScale[hydrogenIndex - 1] = neckScale[i];
1770         }
1771       }
1772     }
1773   }
1774 
1775   /**
1776    * Apply solute radii definitions used to calculate Born radii.
1777    */
1778   public void applySoluteRadii() {
1779     // Set base radii, descreen radii, HCT overlap scale factor and neck scale factor.
1780     for (int i = 0; i < nAtoms; i++) {
1781       udpateSoluteParameters(i);
1782     }
1783 
1784     // Compute "perfect" HCT scale factors.
1785     if (perfectHCTScale) {
1786       double[][] coords = new double[nAtoms][3];
1787       int index = 0;
1788       for (Atom atom : atoms) {
1789         coords[index][0] = atom.getX();
1790         coords[index][1] = atom.getY();
1791         coords[index][2] = atom.getZ();
1792         index++;
1793       }
1794       GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
1795       gaussVol.computeVolumeAndSA(coords);
1796       double[] selfVolumesFractions = gaussVol.getSelfVolumeFractions();
1797       for (int i = 0; i < nAtoms; i++) {
1798         // Use the self volume fractions, plus add the GK overlap scale.
1799         overlapScale[i] = selfVolumesFractions[i] * hctScale;
1800 
1801         // Apply the descreenWithHydrogen flag.
1802         if (!descreenHydrogen && atoms[i].getAtomicNumber() == 1) {
1803           overlapScale[i] = 0.0;
1804         }
1805       }
1806     }
1807   }
1808 
1809   public void setSneck(double sneck_input) {
1810     this.sneck = sneck_input;
1811     initAtomArrays();
1812   }
1813 
1814   public double[] getTanhBetas() {
1815     double[] betas = {beta0, beta1, beta2};
1816     return betas;
1817   }
1818 
1819   public void setTanhBetas(double[] betas) {
1820     this.beta0 = betas[0];
1821     this.beta1 = betas[1];
1822     this.beta2 = betas[2];
1823   }
1824 
1825   void setLambdaFunction(double lPow, double dlPow, double dl2Pow) {
1826     if (lambdaTerm) {
1827       this.lPow = lPow;
1828       this.dlPow = dlPow;
1829       this.dl2Pow = dl2Pow;
1830     } else {
1831       // If the lambdaTerm flag is false, lambda must be set to one.
1832       this.lambda = 1.0;
1833       this.lPow = 1.0;
1834       this.dlPow = 0.0;
1835       this.dl2Pow = 0.0;
1836     }
1837   }
1838 
1839   public enum NonPolarModel {
1840     CAV,
1841     CAV_DISP,
1842     DISP,
1843     SEV_DISP,
1844     GAUSS_DISP,
1845     HYDROPHOBIC_PMF,
1846     BORN_CAV_DISP,
1847     BORN_SOLV,
1848     NONE
1849   }
1850 }