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