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.pme;
39  
40  import edu.rit.pj.IntegerSchedule;
41  import edu.rit.util.Range;
42  import ffx.crystal.Crystal;
43  import ffx.potential.parameters.ForceField;
44  
45  import java.util.logging.Logger;
46  
47  import static java.lang.String.format;
48  import static org.apache.commons.math3.util.FastMath.pow;
49  
50  public class AlchemicalParameters {
51  
52    private static final Logger logger = Logger.getLogger(AlchemicalParameters.class.getName());
53  
54    /**
55     * The polarization model to use.
56     */
57    private final Polarization polarization;
58    /**
59     * The alchemical mode to use.
60     */
61    public final AlchemicalMode mode;
62    /**
63     * Constant α in: r' = sqrt(r^2 + α*(1 - L)^2)
64     */
65    public double permLambdaAlpha = 1.0;
66    /**
67     * Power on L in front of the pairwise multipole potential.
68     */
69    public double permLambdaExponent = 3.0;
70    /**
71     * Begin turning on permanent multipoles at Lambda = 0.4;
72     */
73    public double permLambdaStart = 0.4;
74    /**
75     * Finish turning on permanent multipoles at Lambda = 1.0;
76     */
77    public double permLambdaEnd = 1.0;
78    /**
79     * Start turning on polarization later in the Lambda path to prevent SCF convergence problems when
80     * atoms nearly overlap.
81     */
82    public double polLambdaStart = 0.75;
83    public double polLambdaEnd = 1.0;
84    /**
85     * Power on L in front of the polarization energy.
86     */
87    public double polLambdaExponent = 3.0;
88    /**
89     * Intramolecular electrostatics for the ligand in vapor is included by default.
90     */
91    public boolean doLigandVaporElec = true;
92    /**
93     * Intramolecular electrostatics for the ligand in done in GK implicit solvent.
94     */
95    public boolean doLigandGKElec = false;
96    /**
97     * Condensed phase SCF without the ligand present is included by default. For DualTopologyEnergy
98     * calculations it can be turned off.
99     */
100   public boolean doNoLigandCondensedSCF = true;
101   /**
102    * lAlpha = α*(1 - L)^2
103    */
104   public double lAlpha = 0.0;
105   public double dlAlpha = 0.0;
106   public double d2lAlpha = 0.0;
107   public double dEdLSign = 1.0;
108   /**
109    * lPowPerm = L^permanentLambdaExponent
110    */
111   public double lPowPerm = 1.0;
112   public double dlPowPerm = 0.0;
113   public double d2lPowPerm = 0.0;
114   public boolean doPermanentRealSpace = true;
115   public double permanentScale = 1.0;
116   /**
117    * lPowPol = L^polarizationLambdaExponent
118    */
119   public double lPowPol = 1.0;
120   public double dlPowPol = 0.0;
121   public double d2lPowPol = 0.0;
122   public boolean doPolarization = true;
123   /**
124    * When computing the polarization energy at L there are 3 pieces.
125    *
126    * <p>1.) Upol(1) = The polarization energy computed normally (ie. system with ligand).
127    *
128    * <p>2.) Uenv = The polarization energy of the system without the ligand.
129    *
130    * <p>3.) Uligand = The polarization energy of the ligand by itself.
131    *
132    * <p>Upol(L) = L*Upol(1) + (1-L)*(Uenv + Uligand)
133    *
134    * <p>Set polarizationScale to L for part 1. Set polarizationScale to (1-L) for parts 2 and 3.
135    */
136   public double polarizationScale = 1.0;
137   /**
138    * The polarization Lambda value goes from 0.0 .. 1.0 as the global lambda value varies between
139    * polLambdaStart .. polLambadEnd.
140    */
141   public double polLambda = 1.0;
142   /**
143    * The permanent Lambda value goes from 0.0 .. 1.0 as the global lambda value varies between
144    * permLambdaStart .. permLambdaEnd.
145    */
146   public double permLambda = 1.0;
147   /**
148    * Boundary conditions for the vapor end of the alchemical path.
149    */
150   public Crystal vaporCrystal = null;
151   public int[][][] vaporLists = null;
152   public Range[] vacuumRanges = null;
153   public IntegerSchedule vaporPermanentSchedule = null;
154   public IntegerSchedule vaporEwaldSchedule = null;
155 
156   public AlchemicalParameters(ForceField forceField, boolean lambdaTerm,
157                               boolean nnTerm, Polarization polarization) {
158     this.polarization = polarization;
159 
160     AlchemicalMode tempMode;
161     try {
162       tempMode = AlchemicalMode.valueOf(forceField.getString("ALCHEMICAL_MODE", "OST").toUpperCase());
163     } catch (IllegalArgumentException e) {
164       logger.info(" Invalid value for alchemical-mode; reverting to OST");
165       tempMode = AlchemicalMode.OST;
166     }
167     mode = tempMode;
168 
169     if (lambdaTerm) {
170       // Values of PERMANENT_LAMBDA_ALPHA below 2 can lead to unstable  trajectories.
171       permLambdaAlpha = forceField.getDouble("PERMANENT_LAMBDA_ALPHA", 2.0);
172       if (permLambdaAlpha < 0.0 || permLambdaAlpha > 3.0) {
173         logger.warning("Invalid value for permanent-lambda-alpha (<0.0 || >3.0); reverting to 2.0");
174         permLambdaAlpha = 2.0;
175       }
176 
177       /*
178        A PERMANENT_LAMBDA_EXPONENT of 2 gives a non-zero d2U/dL2 at the
179        beginning of the permanent schedule. Choosing a power of 3 or
180        greater ensures a smooth dU/dL and d2U/dL2 over the schedule.
181 
182        A value of 0.0 is also admissible for when ExtendedSystem is
183        scaling multipoles rather than softcoring them.
184       */
185       permLambdaExponent = forceField.getDouble("PERMANENT_LAMBDA_EXPONENT", 3.0);
186       if (permLambdaExponent < 0.0) {
187         logger.warning("Invalid value for permanent-lambda-exponent (<0.0); reverting to 3.0");
188         permLambdaExponent = 3.0;
189       }
190 
191       /*
192        A POLARIZATION_LAMBDA_EXPONENT of 2 gives a non-zero d2U/dL2 at
193        the beginning of the polarization schedule. Choosing a power of 3
194        or greater ensures a smooth dU/dL and d2U/dL2 over the schedule.
195 
196        A value of 0.0 is also admissible: when polarization is not being
197        softcored but instead scaled, as by ExtendedSystem.
198       */
199       polLambdaExponent = forceField.getDouble("POLARIZATION_LAMBDA_EXPONENT", 3.0);
200       if (polLambdaExponent < 0.0) {
201         logger.warning("Invalid value for polarization-lambda-exponent (<0.0); reverting to 3.0");
202         polLambdaExponent = 3.0;
203       }
204 
205       // Values of PERMANENT_LAMBDA_START below 0.5 can lead to unstable trajectories.
206       permLambdaStart = forceField.getDouble("PERMANENT_LAMBDA_START", 0.4);
207       if (permLambdaStart < 0.0 || permLambdaStart > 1.0) {
208         logger.warning("Invalid value for perm-lambda-start (<0.0 || >1.0); reverting to 0.4");
209         permLambdaStart = 0.4;
210       }
211 
212       // Values of PERMANENT_LAMBDA_END must be greater than permLambdaStart and <= 1.0.
213       permLambdaEnd = forceField.getDouble("PERMANENT_LAMBDA_END", 1.0);
214       if (permLambdaEnd < permLambdaStart || permLambdaEnd > 1.0) {
215         logger.warning("Invalid value for perm-lambda-end (<start || >1.0); reverting to 1.0");
216         permLambdaEnd = 1.0;
217       }
218 
219       /*
220          The POLARIZATION_LAMBDA_START defines the point in the lambda
221          schedule when the condensed phase polarization of the ligand
222          begins to be turned on. If the condensed phase polarization
223          is considered near lambda=0, then SCF convergence is slow,
224          even with Thole damping. In addition, 2 (instead of 1)
225          condensed phase SCF calculations are necessary from the
226          beginning of the window to lambda=1.
227         */
228       polLambdaStart = forceField.getDouble("POLARIZATION_LAMBDA_START", 0.75);
229       if (polLambdaStart < 0.0 || polLambdaStart > 1.0) {
230         logger.warning("Invalid value for polarization-lambda-start; reverting to 0.75");
231         polLambdaStart = 0.75;
232       }
233 
234       /*
235          The POLARIZATION_LAMBDA_END defines the point in the lambda
236          schedule when the condensed phase polarization of ligand has
237          been completely turned on. Values other than 1.0 have not been tested.
238         */
239       polLambdaEnd = forceField.getDouble("POLARIZATION_LAMBDA_END", 1.0);
240       if (polLambdaEnd < polLambdaStart || polLambdaEnd > 1.0) {
241         logger.warning(
242             "Invalid value for polarization-lambda-end (<start || >1.0); reverting to 1.0");
243         polLambdaEnd = 1.0;
244       }
245 
246       // The LAMBDA_VAPOR_ELEC defines if intra-molecular electrostatics of the ligand in vapor
247       // will be considered.
248       if (mode == AlchemicalMode.OST) {
249         doLigandVaporElec = forceField.getBoolean("LIGAND_VAPOR_ELEC", true);
250         doLigandGKElec = forceField.getBoolean("LIGAND_GK_ELEC", false);
251         doNoLigandCondensedSCF = forceField.getBoolean("NO_LIGAND_CONDENSED_SCF", true);
252       } else {
253         doLigandVaporElec = false;
254         doLigandGKElec = false;
255         doNoLigandCondensedSCF = false;
256       }
257     } else if (nnTerm) {
258       permLambdaAlpha = 0.0;
259       permLambdaExponent = 1.0;
260       polLambdaExponent = 1.0;
261       permLambdaStart = 0.0;
262       permLambdaEnd = 1.0;
263       polLambdaStart = 0.0;
264       polLambdaEnd = 1.0;
265       // The LAMBDA_VAPOR_ELEC defines if intramolecular electrostatics of the neural network
266       // atoms in vapor will be removed.
267       if (mode == AlchemicalMode.OST) {
268         doLigandVaporElec = forceField.getBoolean("LIGAND_VAPOR_ELEC", true);
269       } else {
270         doLigandVaporElec = false;
271       }
272       doLigandGKElec = false;
273       doNoLigandCondensedSCF = false;
274     }
275   }
276 
277   public String toString() {
278     StringBuilder sb = new StringBuilder("   Alchemical Parameters\n");
279     sb.append(
280         format(
281             "    Permanent Multipole Range:      %5.3f-%5.3f\n", permLambdaStart, permLambdaEnd));
282     sb.append(format("    Permanent Multipole Softcore Alpha:   %5.3f\n", permLambdaAlpha));
283     sb.append(format("    Permanent Multipole Lambda Exponent:  %5.3f\n", permLambdaExponent));
284     if (polarization != Polarization.NONE) {
285       sb.append(format("    Polarization Lambda Exponent:         %5.3f\n", polLambdaExponent));
286       sb.append(
287           format(
288               "    Polarization Range:             %5.3f-%5.3f\n", polLambdaStart, polLambdaEnd));
289       sb.append(format("    Condensed SCF Without Ligand:         %B\n", doNoLigandCondensedSCF));
290     }
291     if (!doLigandGKElec) {
292       sb.append(format("    Vapor Electrostatics:                 %B\n", doLigandVaporElec));
293     } else {
294       sb.append(format("    GK Electrostatics at L=0:             %B\n", doLigandGKElec));
295     }
296     return sb.toString();
297   }
298 
299   /*
300    * f = sqrt(r^2 + lAlpha)
301    *
302    * df/dL = -alpha * (1.0 - lambda) / f
303    *
304    * g = 1 / sqrt(r^2 + lAlpha)
305    *
306    * dg/dL = alpha * (1.0 - lambda) / (r^2 + lAlpha)^(3/2)
307    *
308    * define dlAlpha = alpha * 1.0 - lambda)
309    *
310    * then df/dL = -dlAlpha / f and dg/dL = dlAlpha * g^3
311    *
312    * Multipoles are turned on from permLambdaStart .. permLambdaEnd.
313    *
314    * @param lambda
315    */
316   public void update(double lambda) {
317     lPowPerm = 1.0;
318     permLambda = 1.0;
319     dlPowPerm = 0.0;
320     d2lPowPerm = 0.0;
321     lAlpha = 0.0;
322     dlAlpha = 0.0;
323     d2lAlpha = 0.0;
324     if (lambda < permLambdaStart) {
325       lPowPerm = 0.0;
326       permLambda = 0.0;
327     } else if (lambda <= permLambdaEnd) {
328       double permWindow = permLambdaEnd - permLambdaStart;
329       double permLambdaScale = 1.0 / permWindow;
330       permLambda = permLambdaScale * (lambda - permLambdaStart);
331       if (mode == AlchemicalMode.OST) {
332         lAlpha = permLambdaAlpha * (1.0 - permLambda) * (1.0 - permLambda);
333         dlAlpha = permLambdaAlpha * (1.0 - permLambda);
334         d2lAlpha = -permLambdaAlpha;
335 
336         lPowPerm = pow(permLambda, permLambdaExponent);
337         dlPowPerm = permLambdaExponent * pow(permLambda, permLambdaExponent - 1.0);
338         d2lPowPerm = 0.0;
339         if (permLambdaExponent >= 2.0) {
340           d2lPowPerm =
341               permLambdaExponent
342                   * (permLambdaExponent - 1.0)
343                   * pow(permLambda, permLambdaExponent - 2.0);
344         }
345 
346         dlAlpha *= permLambdaScale;
347         d2lAlpha *= (permLambdaScale * permLambdaScale);
348         dlPowPerm *= permLambdaScale;
349         d2lPowPerm *= (permLambdaScale * permLambdaScale);
350       }
351     }
352 
353     // Polarization is turned on from polarizationLambdaStart .. polarizationLambdaEnd.
354     lPowPol = 1.0;
355     polLambda = 1.0;
356     dlPowPol = 0.0;
357     d2lPowPol = 0.0;
358     if (lambda < polLambdaStart) {
359       lPowPol = 0.0;
360       polLambda = 0.0;
361     } else if (lambda <= polLambdaEnd) {
362       double polWindow = polLambdaEnd - polLambdaStart;
363       double polLambdaScale = 1.0 / polWindow;
364       polLambda = polLambdaScale * (lambda - polLambdaStart);
365       if (mode == AlchemicalMode.OST) {
366         if (polLambdaExponent > 0.0) {
367           lPowPol = pow(polLambda, polLambdaExponent);
368           if (polLambdaExponent >= 1.0) {
369             dlPowPol = polLambdaExponent * pow(polLambda, polLambdaExponent - 1.0);
370             if (polLambdaExponent >= 2.0) {
371               d2lPowPol =
372                   polLambdaExponent
373                       * (polLambdaExponent - 1.0)
374                       * pow(polLambda, polLambdaExponent - 2.0);
375             }
376           }
377         }
378         // Add the chain rule term due to shrinking the lambda range for the polarization energy.
379         dlPowPol *= polLambdaScale;
380         d2lPowPol *= (polLambdaScale * polLambdaScale);
381       }
382     }
383   }
384 
385   /**
386    * For OST mode, we are calculating analytic dU/dL, d2U/dL2 and d2U/dL/dX for the permanent and
387    * polarization energy terms.
388    * <p>
389    * For SCALE mode, the permanent multipoles and polarizabilities are scaled by lambda for
390    * alchemical atoms.
391    */
392   public enum AlchemicalMode {
393     OST, SCALE
394   }
395 
396 }