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-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.String.format;
42  import static org.apache.commons.math3.util.FastMath.exp;
43  import static org.apache.commons.math3.util.FastMath.min;
44  
45  import ffx.algorithms.dynamics.thermostats.Thermostat;
46  import ffx.potential.ForceFieldEnergy;
47  import ffx.potential.MolecularAssembly;
48  import ffx.potential.bonded.AminoAcidUtils;
49  import ffx.potential.bonded.Atom;
50  import ffx.potential.bonded.Residue;
51  import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
52  import ffx.potential.bonded.ResidueState;
53  import ffx.potential.bonded.Torsion;
54  import ffx.potential.parsers.PDBFilter;
55  import java.io.File;
56  import java.util.ArrayList;
57  import java.util.List;
58  import java.util.concurrent.ThreadLocalRandom;
59  import java.util.logging.Logger;
60  import org.apache.commons.io.FilenameUtils;
61  
62  /**
63   * Orientational Biased Monte Carlo (as applied to chi0 torsion of peptide side-chains).
64   *
65   * <p>As described by Frenkel/Smit in "Understanding Molecular Simulation" Chapter 13.1.2. This uses
66   * the "orientational biasing" method to select chi[0] moves that are frequently accepted.
67   *
68   * @author Stephen D. LuCore
69   */
70  public class RosenbluthOBMC implements MonteCarloListener {
71  
72    private static final Logger logger = Logger.getLogger(RosenbluthOBMC.class.getName());
73  
74    /** MolecularAssembly to operate on. */
75    private final MolecularAssembly molecularAssembly;
76    /** Force field energy to use. */
77    private final ForceFieldEnergy forceFieldEnergy;
78  
79    private final Thermostat thermostat;
80    /** At each move, one of these residues will be chosen as the target. */
81    private final List<Residue> targets;
82    /** Number of mcUpdate() calls (e.g. MD steps) between move proposals. */
83    private final int mcFrequency;
84    /** Size of the trial sets, k. */
85    private final int trialSetSize;
86    /** Keeps track of calls to mcUpdate (e.g. MD steps). */
87    private int steps = 0;
88    /** Rosenbluth factor for the forward move. */
89    private double Wn;
90    /** Rosenbluth factor for the backward move. */
91    private double Wo;
92    /** Counters for proposed and accepted moves. */
93    private int numMovesProposed = 0;
94    /** Creates verbose output. */
95    private StringBuilder report = new StringBuilder();
96    /** Writes PDBs of each trial set and original/proposed configurations. */
97    private boolean writeSnapshots = false;
98  
99    /**
100    * RRMC constructor.
101    *
102    * @param targets Residues to undergo RRMC.
103    * @param mcFrequency Number of MD steps between RRMC proposals.
104    * @param trialSetSize Larger values cost more but increase acceptance.
105    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
106    * @param forceFieldEnergy a {@link ffx.potential.ForceFieldEnergy} object.
107    * @param thermostat a {@link ffx.algorithms.dynamics.thermostats.Thermostat} object.
108    */
109   public RosenbluthOBMC(MolecularAssembly molecularAssembly, ForceFieldEnergy forceFieldEnergy,
110       Thermostat thermostat, List<Residue> targets, int mcFrequency, int trialSetSize) {
111     this.targets = targets;
112     this.mcFrequency = mcFrequency;
113     this.trialSetSize = trialSetSize;
114     this.molecularAssembly = molecularAssembly;
115     this.forceFieldEnergy = forceFieldEnergy;
116     this.thermostat = thermostat;
117   }
118 
119   /**
120    * Constructor for RosenbluthOBMC.
121    *
122    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
123    * @param forceFieldEnergy a {@link ffx.potential.ForceFieldEnergy} object.
124    * @param thermostat a {@link ffx.algorithms.dynamics.thermostats.Thermostat} object.
125    * @param targets a {@link java.util.List} object.
126    * @param mcFrequency a int.
127    * @param trialSetSize a int.
128    * @param writeSnapshots a boolean.
129    */
130   public RosenbluthOBMC(MolecularAssembly molecularAssembly, ForceFieldEnergy forceFieldEnergy,
131       Thermostat thermostat, List<Residue> targets, int mcFrequency, int trialSetSize,
132       boolean writeSnapshots) {
133     this(molecularAssembly, forceFieldEnergy, thermostat, targets, mcFrequency, trialSetSize);
134     this.writeSnapshots = writeSnapshots;
135   }
136 
137   /** {@inheritDoc} */
138   @Override
139   public boolean mcUpdate(double temperature) {
140     steps++;
141     if (steps % mcFrequency == 0) {
142       return mcStep();
143     }
144     return false;
145   }
146 
147   /**
148    * Does all the work for a move. Moveset is a continuous 360 degree spin of the chi[0] torsion.
149    * U_or in Frenkel's notation (uDep here) is the associated torsion energy. Evaluation criterion:
150    * P_accept = Min( 1, (Wn/Wo)*exp(-beta(U[n]-U[o]) )
151    */
152   private boolean mcStep() {
153     numMovesProposed++;
154     boolean accepted;
155 
156     // Select a target residue.
157     int index = ThreadLocalRandom.current().nextInt(targets.size());
158     Residue target = targets.get(index);
159     ResidueState origState = target.storeState();
160     Torsion chi0 = getChiZeroTorsion(target);
161     writeSnapshot("orig");
162 
163     /*
164        Create old and new trial sets, calculate Wn and Wo, and choose a move bn.
165        When doing strictly chi[0] moves, Frenkel/Smit's 'old' and 'new' configurations
166        are the same state. The distinction is made here only to aid in future generalization.
167     */
168     List<MCMove> oldTrialSet = createTrialSet(target, origState, trialSetSize - 1);
169     List<MCMove> newTrialSet = createTrialSet(target, origState, trialSetSize);
170     report = new StringBuilder();
171     report.append(format(" Rosenbluth Rotamer MC Move: %4d\n", numMovesProposed));
172     report.append(format("    residue:   %s\n", target));
173     report.append(format("    chi0:      %s\n", chi0.toString()));
174     MCMove proposal = calculateRosenbluthFactors(target, chi0, origState, oldTrialSet, origState,
175         newTrialSet);
176 
177     /*
178        Calculate the independent portion of the total old-conf energy.
179        Then apply the move and calculate the independent total new-conf energy.
180     */
181     setState(target, origState);
182     writeSnapshot("uIndO");
183     double uIndO = getTotalEnergy() - getTorsionEnergy(chi0);
184     proposal.move();
185     writeSnapshot("uIndN");
186     double uIndN = getTotalEnergy() - getTorsionEnergy(chi0);
187 
188     // Apply acceptance criterion.
189     double temperature = thermostat.getCurrentTemperature();
190     double beta = 1.0 / (R * temperature);
191     double dInd = uIndN - uIndO;
192     double dIndE = exp(-beta * dInd);
193     double criterion = (Wn / Wo) * exp(-beta * (uIndN - uIndO));
194     double metropolis = min(1, criterion);
195     double rng = ThreadLocalRandom.current().nextDouble();
196 
197     report.append(format("    theta:     %3.2f\n", ((RosenbluthChi0Move) proposal).theta));
198     report.append(format("    criterion: %1.4f\n", criterion));
199     report.append(format("       Wn/Wo:     %.2f\n", Wn / Wo));
200     report.append(format("       uIndN,O:  %7.2f\t%7.2f\n", uIndN, uIndO));
201     report.append(format("       dInd(E):  %7.2f\t%7.2f\n", dInd, dIndE));
202     report.append(format("    rng:       %1.4f\n", rng));
203     if (rng < metropolis) {
204       report.append(" Accepted.\n");
205       accepted = true;
206     } else {
207       proposal.revertMove();
208       report.append(" Denied.\n");
209       accepted = false;
210     }
211     logger.info(report.toString());
212 
213     // Cleanup.
214     Wn = 0.0;
215     Wo = 0.0;
216     return accepted;
217   }
218 
219   /**
220    * Generates trial movesets around new and old configurations. This involves loading the
221    * move-dependent energy component into each member of the trial sets.
222    */
223   private List<MCMove> createTrialSet(Residue target, ResidueState state, int setSize) {
224     List<MCMove> moves = new ArrayList<>();
225     // Trial set around old configuration is of size (k-1).
226     setState(target, state);
227     for (int i = 0; i < setSize; i++) {
228       moves.add(new RosenbluthChi0Move(target));
229     }
230     return moves;
231   }
232 
233   /**
234    * Chooses a move, bn, from amongst the new trial set, {b}k, based on the Boltzmann-weighted
235    * orientational energy, U_or[n]. Also calculates Rosenbluth factors for both sets, Wn and Wo. Wn =
236    * sum(i=1:k, exp(-beta * uDep[i])) Wo = exp(-beta * uDep[current]) + sum(i=2:k, exp(-beta *
237    * uDep[i]))
238    */
239   private MCMove calculateRosenbluthFactors(Residue target, Torsion chi0, ResidueState oldConf,
240       List<MCMove> oldTrialSet, ResidueState newConf, List<MCMove> newTrialSet) {
241     double temperature = thermostat.getCurrentTemperature();
242     double beta = 1.0 / (R * temperature);
243 
244     // Initialize and add up Wo.
245     Wo = exp(-beta * getTorsionEnergy(chi0));
246     report.append(format("    TestSet (Old): %5s\t%7s\t\t%7s\n", "uDepO", "uDepOe", "Sum(Wo)"));
247     report.append(format("       Orig %d:   %7.4f\t%7.4f\t\t%7.4f\n", 0, getTorsionEnergy(chi0),
248         exp(-beta * getTorsionEnergy(chi0)), Wo));
249     for (int i = 0; i < oldTrialSet.size(); i++) {
250       setState(target, oldConf);
251       MCMove move = oldTrialSet.get(i);
252       move.move();
253       double uDepO = getTorsionEnergy(chi0);
254       double uDepOe = exp(-beta * uDepO);
255       Wo += uDepOe;
256       if (i < 5 || i >= oldTrialSet.size() - 5) {
257         report.append(format("       Prop %d:   %7.4f\t%7.4f\t\t%7.4f\n", i + 1, uDepO, uDepOe, Wo));
258         writeSnapshot("ots");
259       } else if (i == 5) {
260         report.append("        ... \n");
261       }
262     }
263 
264     // Initialize and add up Wn.  Record dependent energy of each trial set member.
265     Wn = 0.0;
266     double[] uDepN = new double[newTrialSet.size()];
267     double[] uDepNe = new double[newTrialSet.size()];
268     report.append(format("    TestSet (New): %5s\t%7s\t\t%7s\n", "uDepN", "uDepNe", "Sum(Wn)"));
269     for (int i = 0; i < newTrialSet.size(); i++) {
270       setState(target, newConf);
271       MCMove move = newTrialSet.get(i);
272       move.move();
273       uDepN[i] = getTorsionEnergy(chi0);
274       uDepNe[i] = exp(-beta * uDepN[i]);
275       Wn += uDepNe[i];
276       if (i < 5 || i >= newTrialSet.size() - 5) {
277         report.append(
278             format("       Prop %d:   %7.4f\t%7.4f\t\t%7.4f\n", i, uDepN[i], uDepNe[i], Wn));
279         writeSnapshot("nts");
280       } else if (i == 5) {
281         report.append("        ... \n");
282       }
283     }
284     setState(target, oldConf);
285 
286     // Choose a proposal move from the new trial set.
287     MCMove proposal = null;
288     double rng = ThreadLocalRandom.current().nextDouble(Wn);
289     double running = 0.0;
290     for (int i = 0; i < newTrialSet.size(); i++) {
291       running += uDepNe[i];
292       if (rng < running) {
293         proposal = newTrialSet.get(i);
294         double prob = uDepNe[i] / Wn * 100;
295         report.append(
296             format("       Chose %d   %7.4f\t%7.4f\t  %4.1f%%\n", i, uDepN[i], uDepNe[i], prob));
297         break;
298       }
299     }
300     if (proposal == null) {
301       logger.severe("Programming error.");
302     }
303     return proposal;
304   }
305 
306   private double getTotalEnergy() {
307     double[] x = new double[forceFieldEnergy.getNumberOfVariables() * 3];
308     forceFieldEnergy.getCoordinates(x);
309     return forceFieldEnergy.energy(x);
310   }
311 
312   private double getTorsionEnergy(Torsion torsion) {
313     return torsion.energy(false);
314   }
315 
316   private Torsion getChiZeroTorsion(Residue residue) {
317     AminoAcid3 name = AminoAcidUtils.AminoAcid3.valueOf(residue.getName());
318     List<Torsion> torsions = residue.getTorsionList();
319     switch (name) {
320       case VAL -> {
321         Atom N = (Atom) residue.getAtomNode("N");
322         Atom CA = (Atom) residue.getAtomNode("CA");
323         Atom CB = (Atom) residue.getAtomNode("CB");
324         Atom CG1 = (Atom) residue.getAtomNode("CG1");
325         for (Torsion torsion : torsions) {
326           if (torsion.compare(N, CA, CB, CG1)) {
327             return torsion;
328           }
329         }
330       }
331       case ILE -> {
332         Atom N = (Atom) residue.getAtomNode("N");
333         Atom CA = (Atom) residue.getAtomNode("CA");
334         Atom CB = (Atom) residue.getAtomNode("CB");
335         Atom CG1 = (Atom) residue.getAtomNode("CG1");
336         for (Torsion torsion : torsions) {
337           if (torsion.compare(N, CA, CB, CG1)) {
338             return torsion;
339           }
340         }
341       }
342       case SER -> {
343         Atom N = (Atom) residue.getAtomNode("N");
344         Atom CA = (Atom) residue.getAtomNode("CA");
345         Atom CB = (Atom) residue.getAtomNode("CB");
346         Atom OG = (Atom) residue.getAtomNode("OG");
347         for (Torsion torsion : torsions) {
348           if (torsion.compare(N, CA, CB, OG)) {
349             return torsion;
350           }
351         }
352       }
353       case THR -> {
354         Atom N = (Atom) residue.getAtomNode("N");
355         Atom CA = (Atom) residue.getAtomNode("CA");
356         Atom CB = (Atom) residue.getAtomNode("CB");
357         Atom OG1 = (Atom) residue.getAtomNode("OG1");
358         for (Torsion torsion : torsions) {
359           if (torsion.compare(N, CA, CB, OG1)) {
360             return torsion;
361           }
362         }
363       }
364       case CYX -> {
365         Atom N = (Atom) residue.getAtomNode("N");
366         Atom CA = (Atom) residue.getAtomNode("CA");
367         Atom CB = (Atom) residue.getAtomNode("CB");
368         Atom SG = (Atom) residue.getAtomNode("SG");
369         for (Torsion torsion : torsions) {
370           if (torsion.compare(N, CA, CB, SG)) {
371             return torsion;
372           }
373         }
374       }
375       case CYD -> {
376         Atom N = (Atom) residue.getAtomNode("N");
377         Atom CA = (Atom) residue.getAtomNode("CA");
378         Atom CB = (Atom) residue.getAtomNode("CB");
379         Atom SG = (Atom) residue.getAtomNode("SG");
380         for (Torsion torsion : torsions) {
381           if (torsion.compare(N, CA, CB, SG)) {
382             return torsion;
383           }
384         }
385       }
386       default -> { // All other residues' chi[0] are defined by N,CA,CB,CG.
387         Atom N = (Atom) residue.getAtomNode("N");
388         Atom CA = (Atom) residue.getAtomNode("CA");
389         Atom CB = (Atom) residue.getAtomNode("CB");
390         Atom CG = (Atom) residue.getAtomNode("CG");
391         for (Torsion torsion : torsions) {
392           if (torsion.compare(N, CA, CB, CG)) {
393             return torsion;
394           }
395         }
396         logger.info("Couldn't find chi[0] for residue " + residue);
397         return null;
398       }
399     }
400     logger.info("Couldn't find chi[0] for residue " + residue);
401     return null;
402   }
403 
404   /**
405    * Calls through to residue.revertState() but also updates the Torsion objects associated with that
406    * residue (so they contain appropriate chi values).
407    */
408   private void setState(Residue target, ResidueState state) {
409     target.revertState(state);
410     for (Torsion torsion : target.getTorsionList()) {
411       torsion.update();
412     }
413   }
414 
415   private void writeSnapshot(String suffix) {
416     if (!writeSnapshots) {
417       return;
418     }
419     String filename =
420         FilenameUtils.removeExtension(molecularAssembly.getFile().toString()) + "." + suffix + "-"
421             + numMovesProposed;
422     File file = new File(filename);
423     PDBFilter writer = new PDBFilter(file, molecularAssembly, null, null);
424     writer.writeFile(file, false);
425   }
426 }