1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024.
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.numerics.estimator;
39
40 import static ffx.numerics.estimator.EstimateBootstrapper.getBootstrapIndices;
41 import static java.util.Arrays.copyOf;
42 import static org.apache.commons.math3.util.FastMath.exp;
43 import static org.apache.commons.math3.util.FastMath.log;
44
45 import ffx.utilities.Constants;
46
47 import java.util.Random;
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50 import java.util.stream.IntStream;
51
52 /**
53 * The Zwanzig class implements exponential averaging/free energy perturbation using the Zwanzig
54 * relationship, in either the forwards or backwards direction (not both).
55 *
56 * @author Michael J. Schnieders
57 * @author Jacob M. Litman
58 * @since 1.0
59 */
60 public class Zwanzig extends SequentialEstimator implements BootstrappableEstimator {
61
62 private static final Logger logger = Logger.getLogger(SequentialEstimator.class.getName());
63 public final Directionality directionality;
64 private final boolean forwards;
65 /**
66 * Number of windows.
67 */
68 private final int nWindows;
69 /**
70 * Free energy difference for each window.
71 */
72 private final double[] windowFreeEnergyDifferences;
73 private final double[] windowEnthalpies;
74 /**
75 * Free energy difference uncertainty for each window.
76 */
77 private final double[] windowFreeEnergyUncertainties;
78 private final Random random;
79 /**
80 * Total free energy difference as a sum over windows.
81 */
82 private double totalFreeEnergyDifference;
83 /**
84 * Total free energy difference uncertainty: totalFreeEnergyUncertainty = Sqrt [ Sum over Windows [
85 * Window Uncertainty Squared ] ]
86 */
87 private double totalFreeEnergyUncertainty;
88
89 /**
90 * Total Enthalpy
91 */
92 private double totalEnthalpy;
93
94 /**
95 * Estimates a free energy using the Zwanzig relationship. The temperature array can be of length 1
96 * if all elements are meant to be the same temperature.
97 *
98 * <p>The first dimension of the energies arrays corresponds to the lambda values/windows. The
99 * second dimension (can be of uneven length) corresponds to potential energies of snapshots
100 * sampled from that lambda value, calculated either at that lambda value, the lambda value below,
101 * or the lambda value above. The arrays energiesLow[0] and energiesHigh[n-1] is expected to be all
102 * NaN.
103 *
104 * @param lambdaValues Values of lambda dynamics was run at.
105 * @param energiesLow Potential energies of trajectory L at lambda L-dL. Ignored for forwards
106 * FEP.
107 * @param energiesAt Potential energies of trajectory L at lambda L.
108 * @param energiesHigh Potential energies of trajectory L at lambda L+dL. Ignored for backwards
109 * FEP.
110 * @param temperature Temperature each lambda window was run at (single-element indicates
111 * identical temperatures).
112 * @param directionality Forwards vs. backwards FEP.
113 */
114 public Zwanzig(double[] lambdaValues, double[][] energiesLow, double[][] energiesAt,
115 double[][] energiesHigh,
116 double[] temperature, Directionality directionality) {
117 super(lambdaValues, energiesLow, energiesAt, energiesHigh, temperature);
118 this.directionality = directionality;
119 nWindows = nTrajectories - 1;
120
121 windowFreeEnergyDifferences = new double[nWindows];
122 windowFreeEnergyUncertainties = new double[nWindows];
123 windowEnthalpies = new double[nWindows];
124
125 forwards = directionality.equals(Directionality.FORWARDS);
126 random = new Random();
127
128 estimateDG();
129 }
130
131 /**
132 * {@inheritDoc}
133 */
134 @Override
135 public Zwanzig copyEstimator() {
136 return new Zwanzig(lamValues, eLow, eAt, eHigh, temperatures, directionality);
137 }
138
139 /**
140 * {@inheritDoc}
141 */
142 @Override
143 public final void estimateDG(final boolean randomSamples) {
144 double cumDG = 0;
145
146 Level warningLevel = randomSamples ? Level.FINE : Level.WARNING;
147
148 for (int i = 0; i < nWindows; i++) {
149 int windowIndex = forwards ? 0 : 1;
150 windowIndex += i;
151 double[] e1 = eAt[windowIndex];
152 double[] e2 = forwards ? eHigh[windowIndex] : eLow[windowIndex];
153 double sign = forwards ? 1.0 : -1.0;
154 double beta = -temperatures[windowIndex] * Constants.R;
155 double invBeta = 1.0 / beta;
156
157 int len = e1.length;
158 if (len == 0) {
159 logger.log(warningLevel, " Skipping frame " + i + " due to lack of snapshots!");
160 continue;
161 }
162
163 // With no iteration-to-convergence, generating a fresh random index is OK.
164 int[] samples =
165 randomSamples ? getBootstrapIndices(len, random) : IntStream.range(0, len).toArray();
166
167 double sum = 0.0;
168 double num = 0.0;
169 double uAve = 0.0;
170 for (int indJ = 0; indJ < len; indJ++) {
171 int j = samples[indJ];
172 double dE = e2[j] - e1[j];
173 if (sign > 0) {
174 uAve += e1[j];
175 var term = exp(invBeta * dE);
176 sum += term;
177 num += e2[j] * term;
178 } else {
179 uAve += e2[j];
180 var term = exp(invBeta * dE);
181 sum += term;
182 num += e1[j] * term;
183 }
184 }
185
186 uAve = uAve / len;
187
188 double dG = sign * beta * log(sum / len);
189 windowFreeEnergyDifferences[i] = dG;
190 windowEnthalpies[i] = sign * ((num / sum) - uAve);
191 windowFreeEnergyUncertainties[i] = 0.0;
192 totalEnthalpy += windowEnthalpies[i];
193 cumDG += dG;
194 }
195
196 totalFreeEnergyDifference = cumDG;
197 totalFreeEnergyUncertainty = 0.0;
198 }
199
200 /**
201 * {@inheritDoc}
202 */
203 @Override
204 public final void estimateDG() {
205 estimateDG(false);
206 }
207
208 /**
209 * {@inheritDoc}
210 */
211 @Override
212 public double[] getBinEnergies() {
213 return copyOf(windowFreeEnergyDifferences, nWindows);
214 }
215
216 /**
217 * {@inheritDoc}
218 */
219 @Override
220 public double[] getBinUncertainties() {
221 return copyOf(windowFreeEnergyUncertainties, nWindows);
222 }
223
224 /**
225 * {@inheritDoc}
226 */
227 @Override
228 public double getFreeEnergy() {
229 return totalFreeEnergyDifference;
230 }
231
232 /**
233 * {@inheritDoc}
234 */
235 @Override
236 public double getUncertainty() {
237 return totalFreeEnergyUncertainty;
238 }
239
240 /**
241 * {@inheritDoc}
242 */
243 @Override
244 public int numberOfBins() {
245 return nWindows;
246 }
247
248 /**
249 * {@inheritDoc}
250 */
251 @Override
252 public double[] getBinEnthalpies() {
253 return copyOf(windowEnthalpies, nWindows);
254 }
255
256 /**
257 * Directionality of the Zwanzig estimation (forwards perturbation or backwards perturbation).
258 * TODO: Implement bidirectional Zwanzig with simple estimation (i.e. 0.5*(forwards + backward)).
259 */
260 public enum Directionality {
261 FORWARDS,
262 BACKWARDS
263 }
264 }