1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
64
65
66
67
68
69
70 public class RosenbluthOBMC implements MonteCarloListener {
71
72 private static final Logger logger = Logger.getLogger(RosenbluthOBMC.class.getName());
73
74
75 private final MolecularAssembly molecularAssembly;
76
77 private final ForceFieldEnergy forceFieldEnergy;
78
79 private final Thermostat thermostat;
80
81 private final List<Residue> targets;
82
83 private final int mcFrequency;
84
85 private final int trialSetSize;
86
87 private int steps = 0;
88
89 private double Wn;
90
91 private double Wo;
92
93 private int numMovesProposed = 0;
94
95 private StringBuilder report = new StringBuilder();
96
97 private boolean writeSnapshots = false;
98
99
100
101
102
103
104
105
106
107
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
121
122
123
124
125
126
127
128
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
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
149
150
151
152 private boolean mcStep() {
153 numMovesProposed++;
154 boolean accepted;
155
156
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
165
166
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
179
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
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
214 Wn = 0.0;
215 Wo = 0.0;
216 return accepted;
217 }
218
219
220
221
222
223 private List<MCMove> createTrialSet(Residue target, ResidueState state, int setSize) {
224 List<MCMove> moves = new ArrayList<>();
225
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
235
236
237
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
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
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
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 -> {
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
406
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 }