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.algorithms.mc;
39
40 import static ffx.utilities.Constants.R;
41 import static java.lang.Double.isFinite;
42 import static java.lang.String.format;
43 import static org.apache.commons.math3.util.FastMath.exp;
44 import static org.apache.commons.math3.util.FastMath.min;
45
46 import java.util.ArrayList;
47 import java.util.List;
48 import java.util.Random;
49 import java.util.logging.Logger;
50
51 /**
52 * The BoltzmannMC abstract class is a skeleton for Boltzmann-weighted Metropolis Monte Carlo
53 * simulations.
54 *
55 * @author Michael J. Schnieders
56 * @author Jacob M. Litman
57 */
58 public abstract class BoltzmannMC implements MetropolisMC {
59
60 /** Constant <code>logger</code> */
61 private static final Logger logger = Logger.getLogger(BoltzmannMC.class.getName());
62
63 protected final Random random = new Random();
64 /**
65 * Room temperature in Kelvin.
66 */
67 private double temperature = 298.15;
68 /**
69 * Constant factor for Monte Carlo moves: 1.0 / (kB * T)
70 */
71 private double invKBT = 1.0 / (R * temperature);
72
73 private boolean print = true;
74 private double e1 = 0.0;
75 private double e2 = 0.0;
76 private double lastE = 0.0;
77 private boolean lastAccept = false;
78
79 /**
80 * Boltzmann-weighted acceptance probability
81 *
82 * @param invKT 1.0 / (Boltzmann constant * temperature)
83 * @param e1 Energy before move
84 * @param e2 Proposed energy
85 * @return Chance of accepting this move
86 */
87 public static double acceptChance(double invKT, double e1, double e2) {
88 double dE = e2 - e1;
89 return min(exp(-invKT * dE), 1.0);
90 }
91
92 /**
93 * Boltzmann-weighted acceptance probability
94 *
95 * @param random Source of randomness
96 * @param invKT 1.0 / (Boltzmann constant * temperature)
97 * @param e1 Energy before move
98 * @param e2 Proposed energy
99 * @return Whether to accept the move.
100 */
101 public static boolean evaluateMove(Random random, double invKT, double e1, double e2) {
102 boolean e1Finite = isFinite(e1);
103 boolean e2Finite = isFinite(e2);
104 if (!e1Finite && !e2Finite) {
105 throw new IllegalArgumentException(
106 format(" Attempting to evaluate " + "move with non-finite energies %f and %f!", e1, e2));
107 } else if (!e1Finite) {
108 throw new IllegalArgumentException(
109 format(" Attempting to evaluate move " + "from non-finite %f to finite %.4f!", e1, e2));
110 } else if (!e2Finite) {
111 throw new IllegalArgumentException(
112 format(" Attempting to evaluate move " + "from finite %.4f to non-finite %.4f!", e1, e2));
113 }
114 if (e2 <= e1) {
115 return true;
116 } else {
117 // p = exp(-dE / kB*T)
118 double prob = acceptChance(invKT, e1, e2);
119 assert (prob >= 0.0
120 && prob <= 1.0) : "The probability of an MC move up in energy should be between [0 ..1].";
121 double trial = random.nextDouble();
122 return (trial <= prob);
123 }
124 }
125
126 /**
127 * {@inheritDoc}
128 *
129 * <p>Criterion for accept/reject a move; intended to be used mostly internally.
130 */
131 @Override
132 public boolean evaluateMove(double e1, double e2) {
133 return evaluateMove(random, invKBT, e1, e2);
134 }
135
136 /** {@inheritDoc} */
137 @Override
138 public boolean getAccept() {
139 return lastAccept;
140 }
141
142 /** {@inheritDoc} */
143 @Override
144 public double getE1() {
145 return e1;
146 }
147
148 /** {@inheritDoc} */
149 @Override
150 public double getE2() {
151 return e2;
152 }
153
154 /** {@inheritDoc} */
155 @Override
156 public double getTemperature() {
157 return temperature;
158 }
159
160 /** {@inheritDoc} */
161 @Override
162 public void setTemperature(double temp) {
163 temperature = temp;
164 invKBT = 1.0 / (R * temperature);
165 }
166
167 /** {@inheritDoc} */
168 @Override
169 public double lastEnergy() {
170 return lastE;
171 }
172
173 /** {@inheritDoc} */
174 @Override
175 public boolean mcStep(MCMove move) {
176 return mcStep(move, currentEnergy());
177 }
178
179 /** {@inheritDoc} */
180 @Override
181 public boolean mcStep(MCMove move, double en1) {
182 List<MCMove> moveList = new ArrayList<>(1);
183 moveList.add(move);
184 return mcStep(moveList, en1);
185 }
186
187 /** {@inheritDoc} */
188 @Override
189 public boolean mcStep(List<MCMove> moves) {
190 return mcStep(moves, currentEnergy());
191 }
192
193 /**
194 * {@inheritDoc}
195 *
196 * <p>Performs a Boltzmann-weighted Monte Carlo step with an arbitrary list of moves and a defined
197 * starting energy. The list of MCMoves should be of a type with O(1) element access, as the
198 * current implementation utilizes an indexed for loop.
199 */
200 @Override
201 public boolean mcStep(List<MCMove> moves, double en1) {
202 storeState();
203 e1 = en1;
204
205 int nMoves = moves.size();
206 for (MCMove move : moves) {
207 move.move();
208 }
209
210 lastE = currentEnergy(); // Is reset to e1 if move rejected.
211 e2 = lastE;
212 if (evaluateMove(e1, e2)) {
213 lastAccept = true;
214 if (print) {
215 logger.info(format(" Monte Carlo step accepted: e1 -> e2 %10.6f -> %10.6f", e1, e2));
216 }
217 return true;
218 } else {
219 lastAccept = false;
220 for (int i = nMoves - 1; i >= 0; i--) {
221 moves.get(i).revertMove();
222 }
223 lastE = e1;
224 if (print) {
225 logger.info(format(" Monte Carlo step rejected: e1 -> e2 %10.6f -> %10.6f", e1, e2));
226 }
227 return false;
228 }
229 }
230
231 /** {@inheritDoc} */
232 @Override
233 public void setPrint(boolean print) {
234 this.print = print;
235 }
236
237 /**
238 * Set the random seed.
239 *
240 * @param randomSeed The seed.
241 */
242 protected void setRandomSeed(int randomSeed) {
243 random.setSeed(randomSeed);
244 }
245
246 /**
247 * Must return the current energy of the system.
248 *
249 * @return Current system energy
250 */
251 protected abstract double currentEnergy();
252
253 /**
254 * Store the state for reverting a move. Must be properly implemented for revertStep() to function
255 * properly; otherwise, the implementation of revertStep() should throw an
256 * OperationNotSupportedException.
257 */
258 protected abstract void storeState();
259 }