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 ffx.algorithms.dynamics.thermostats.Thermostat;
41 import ffx.potential.ForceFieldEnergy;
42 import ffx.potential.MolecularAssembly;
43 import ffx.potential.bonded.AminoAcidUtils;
44 import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
45 import ffx.potential.bonded.Residue;
46 import ffx.potential.parsers.PDBFilter;
47 import org.apache.commons.io.FilenameUtils;
48
49 import java.io.File;
50 import java.util.List;
51 import java.util.concurrent.ThreadLocalRandom;
52 import java.util.logging.Logger;
53
54 import static java.lang.String.format;
55 import static org.apache.commons.math3.util.FastMath.min;
56
57
58
59
60
61
62
63
64
65
66 public class RosenbluthCBMC implements MonteCarloListener {
67
68 private static final Logger logger = Logger.getLogger(RosenbluthCBMC.class.getName());
69
70
71 private final MolecularAssembly molecularAssembly;
72
73 private final ForceFieldEnergy forceFieldEnergy;
74
75 private final double temperature;
76
77 private final List<Residue> targets;
78
79 private final int mcFrequency;
80
81 private final int trialSetSize;
82
83 private int steps = 0;
84
85 private int numMovesProposed = 0;
86
87 private final boolean writeSnapshots;
88
89 private PDBFilter writer = null;
90
91
92
93
94
95
96
97
98
99
100
101
102 public RosenbluthCBMC(MolecularAssembly molecularAssembly, ForceFieldEnergy ffe,
103 Thermostat thermostat, List<Residue> targets, int mcFrequency, int trialSetSize,
104 boolean writeSnapshots) {
105 this.targets = targets;
106 this.mcFrequency = mcFrequency;
107 this.trialSetSize = trialSetSize;
108 this.molecularAssembly = molecularAssembly;
109 this.forceFieldEnergy = ffe;
110 this.writeSnapshots = writeSnapshots;
111
112 if (thermostat != null) {
113 temperature = thermostat.getTargetTemperature();
114 } else {
115 temperature = 298.15;
116 }
117
118 for (int i = targets.size() - 1; i >= 0; i--) {
119 AminoAcid3 name = AminoAcidUtils.AminoAcid3.valueOf(targets.get(i).getName());
120 if (name == AminoAcidUtils.AminoAcid3.GLY || name == AminoAcidUtils.AminoAcid3.PRO
121 || name == AminoAcidUtils.AminoAcid3.ALA) {
122 targets.remove(i);
123 }
124 }
125 if (targets.isEmpty()) {
126 logger.severe(" Empty target list for CMBC.");
127 }
128 }
129
130
131
132
133
134
135 public boolean cbmcStep() {
136 numMovesProposed++;
137 boolean accepted;
138
139
140 int index = ThreadLocalRandom.current().nextInt(targets.size());
141 Residue target = targets.get(index);
142 RosenbluthChiAllMove cbmcMove = new RosenbluthChiAllMove(molecularAssembly, target, trialSetSize,
143 forceFieldEnergy, temperature, writeSnapshots, numMovesProposed, true);
144 if (cbmcMove.getMode() == RosenbluthChiAllMove.MODE.CHEAP) {
145 return cbmcMove.wasAccepted();
146 }
147 double Wn = cbmcMove.getWn();
148 double Wo = cbmcMove.getWo();
149 double criterion = min(1, Wn / Wo);
150 double rng = ThreadLocalRandom.current().nextDouble();
151 logger.info(format(" rng: %5.2f", rng));
152 if (rng < criterion) {
153 cbmcMove.move();
154 logger.info(format(" Accepted! Energy: %.4f\n", cbmcMove.finalEnergy));
155 accepted = true;
156 write();
157 } else {
158 logger.info(" Denied.\n");
159 accepted = false;
160 }
161
162 return accepted;
163 }
164
165
166
167
168
169
170 public boolean controlStep() {
171 int index = ThreadLocalRandom.current().nextInt(targets.size());
172 Residue target = targets.get(index);
173 RosenbluthChiAllMove cbmcMove = new RosenbluthChiAllMove(molecularAssembly, target, -1,
174 forceFieldEnergy, temperature, false, numMovesProposed, true);
175 return cbmcMove.wasAccepted();
176 }
177
178
179 @Override
180 public boolean mcUpdate(double temperature) {
181 steps++;
182 if (steps % mcFrequency == 0) {
183 return cbmcStep();
184 }
185 return false;
186 }
187
188
189 private void write() {
190 if (writer == null) {
191 writer = new PDBFilter(molecularAssembly.getFile(), molecularAssembly, null, null);
192 }
193 String filename = molecularAssembly.getFile().getAbsolutePath();
194 if (!filename.contains("_mc")) {
195 filename = FilenameUtils.removeExtension(filename) + "_mc.pdb";
196 }
197 File file = new File(filename);
198 writer.writeFile(file, false);
199 }
200 }