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