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.commands.test;
39
40 import edu.rit.pj.Comm;
41 import ffx.algorithms.cli.AlgorithmsCommand;
42 import ffx.algorithms.cli.DynamicsOptions;
43 import ffx.algorithms.cli.RepExOptions;
44 import ffx.algorithms.dynamics.MDEngine;
45 import ffx.algorithms.dynamics.MolecularDynamics;
46 import ffx.algorithms.dynamics.MolecularDynamicsOpenMM;
47 import ffx.algorithms.dynamics.PhReplicaExchange;
48 import ffx.numerics.Potential;
49 import ffx.potential.ForceFieldEnergy;
50 import ffx.potential.cli.AtomSelectionOptions;
51 import ffx.potential.cli.WriteoutOptions;
52 import ffx.potential.extended.ExtendedSystem;
53 import ffx.potential.parsers.SystemFilter;
54 import ffx.potential.parsers.XPHFilter;
55 import ffx.potential.parsers.XYZFilter;
56 import ffx.utilities.FFXBinding;
57 import org.apache.commons.configuration2.CompositeConfiguration;
58 import org.apache.commons.io.FilenameUtils;
59 import picocli.CommandLine.Command;
60 import picocli.CommandLine.Mixin;
61 import picocli.CommandLine.Option;
62 import picocli.CommandLine.Parameters;
63
64 import java.io.BufferedReader;
65 import java.io.BufferedWriter;
66 import java.io.File;
67 import java.io.FileReader;
68 import java.io.FileWriter;
69 import java.io.IOException;
70 import java.util.Arrays;
71 import java.util.Collections;
72 import java.util.List;
73
74 import static java.lang.String.format;
75
76
77
78
79
80
81
82
83 @Command(description = " Run constant pH dynamics on a system.", name = "PhDynamics")
84 public class PhDynamics extends AlgorithmsCommand {
85
86 @Mixin
87 private AtomSelectionOptions atomSelectionOptions;
88
89 @Mixin
90 private DynamicsOptions dynamicsOptions;
91
92 @Mixin
93 private WriteoutOptions writeOutOptions;
94
95 @Mixin
96 private RepExOptions repEx;
97
98
99
100
101 @Option(names = {"--pH", "--constantPH"}, paramLabel = "7.4", defaultValue = "7.4",
102 description = "Constant pH value for molecular dynamics")
103 private double pH;
104
105 @Option(names = {"--titrationSteps"}, paramLabel = "10", defaultValue = "10",
106 description = "Number of steps done titrating protons on CPU in one cycle")
107 private int titrSteps;
108
109 @Option(names = {"--coordinateSteps"}, paramLabel = "100", defaultValue = "100",
110 description = "Number of steps done propagating coordinates only on GPU in one cycle")
111 private int coordSteps;
112
113 @Option(names = {"--cycles", "--OMMcycles"}, paramLabel = "5", defaultValue = "5",
114 description = "Number of times to cycle between titrating protons on CPU and propagating coordinates only on GPU")
115 private int cycles;
116
117 @Option(names = {"--titrationReport", "--esvLog"}, paramLabel = "0.25 (psec)", defaultValue = "0.25",
118 description = "Interval in psec to report ESV energy and lambdas when cycling between GPU and CPU.")
119 private double titrReport;
120
121 @Option(names = {"--pHGaps"}, paramLabel = "1", defaultValue = "1",
122 description = "pH gap between replica exchange windows.")
123 private double pHGap;
124
125 @Option(names = {"--initDynamics"}, paramLabel = "10000", defaultValue = "10000",
126 description = "Number of initialization steps to take before replica exchange windows start.")
127 private int initDynamics;
128
129 @Option(names = "--sort", paramLabel = "false", defaultValue = "false",
130 description = "Sort archive files by pH")
131 private boolean sort;
132
133 @Option(names = "--printRatioData", paramLabel = "false", defaultValue = "false",
134 description = "Print out the protonation ratios from throughout the simulation at the end")
135 private boolean printRatio;
136
137
138
139
140 @Parameters(arity = "1..*", paramLabel = "files",
141 description = "XYZ or PDB input files.")
142 private String filename;
143
144
145
146
147
148 private Potential potential;
149 public MolecularDynamics molecularDynamics = null;
150
151
152
153
154 public PhDynamics() {
155 super();
156 }
157
158
159
160
161
162 public PhDynamics(FFXBinding binding) {
163 super(binding);
164 }
165
166
167
168
169
170 public PhDynamics(String[] args) {
171 super(args);
172 }
173
174 public MolecularDynamics getMolecularDynamics() {
175 return molecularDynamics;
176 }
177
178 public Potential getPotentialObject() {
179 return potential;
180 }
181
182
183
184
185 @Override
186 public PhDynamics run() {
187
188 if (!init()) {
189 return this;
190 }
191
192 dynamicsOptions.init();
193
194 activeAssembly = getActiveAssembly(filename);
195 if (activeAssembly == null) {
196 logger.info(helpString());
197 return this;
198 }
199
200 String filename = activeAssembly.getFile().getAbsolutePath();
201
202 atomSelectionOptions.setActiveAtoms(activeAssembly);
203
204 potential = activeAssembly.getPotentialEnergy();
205
206
207 File esv = new File(FilenameUtils.removeExtension(filename) + ".esv");
208 if (!esv.exists()) {
209 esv = null;
210 }
211
212
213 ExtendedSystem esvSystem = new ExtendedSystem(activeAssembly, pH, esv);
214 esvSystem.setConstantPh(pH);
215 if (potential instanceof ForceFieldEnergy) {
216 ((ForceFieldEnergy) potential).attachExtendedSystem(esvSystem);
217 }
218
219 int numESVs = esvSystem.getExtendedResidueList().size();
220 logger.info(format(" Attached extended system with %d residues.", numESVs));
221
222 double[] x = new double[potential.getNumberOfVariables()];
223 potential.getCoordinates(x);
224 potential.energy(x, true);
225 SystemFilter systemFilter = algorithmFunctions.getFilter();
226 if (systemFilter instanceof XYZFilter) {
227 XPHFilter xphFilter = new XPHFilter(activeAssembly.getFile(), activeAssembly, activeAssembly.getForceField(),
228 activeAssembly.getProperties(), esvSystem);
229 xphFilter.readFile();
230 logger.info("Reading ESV lambdas from XPH file");
231 potential.getCoordinates(x);
232 potential.energy(x, true);
233 }
234
235 logger.info("\n Running molecular dynamics on " + filename);
236
237 molecularDynamics =
238 dynamicsOptions.getDynamics(writeOutOptions, potential, activeAssembly, algorithmListener);
239
240 molecularDynamics.attachExtendedSystem(esvSystem, titrReport);
241
242
243 File dyn = new File(FilenameUtils.removeExtension(filename) + ".dyn");
244 if (!dyn.exists()) {
245 dyn = null;
246 }
247
248 if (!(molecularDynamics instanceof MolecularDynamicsOpenMM)) {
249 if (repEx.getRepEx()) {
250 Comm world = Comm.world();
251 int size = world.size();
252 if (size == 1) {
253 System.exit(0);
254 }
255
256 String pHWindows = Arrays.toString(PhReplicaExchange.setEvenSpacePhLadder(pHGap, pH, size));
257 pHWindows = pHWindows.replace("[", "")
258 .replace("]", "")
259 .replace(",", " ");
260 CompositeConfiguration properties = activeAssembly.getProperties();
261 pHWindows = properties.getString("pH.Windows", pHWindows);
262 String[] temp = pHWindows.split(" +");
263 if (temp.length != size) {
264 logger.severe("pHLadder specified in properties/key file has " +
265 "incorrect number of windows given world.size()");
266 }
267 double[] pHLadder = new double[size];
268 for (int i = 0; i < temp.length; i++) {
269 pHLadder[i] = Double.parseDouble(temp[i]);
270 }
271
272 logger.info("\n Running replica exchange molecular dynamics on " + filename);
273 int rank = world.rank();
274
275 File structureFile = new File(filename);
276 final String newMolAssemblyFile = structureFile.getParent() + File.separator + rank +
277 File.separator + structureFile.getName();
278 logger.info(" Set activeAssembly filename: " + newMolAssemblyFile);
279 activeAssembly.setFile(new File(newMolAssemblyFile));
280
281 try {
282 PhReplicaExchange pHReplicaExchange = new PhReplicaExchange(molecularDynamics, structureFile, pH, pHLadder,
283 dynamicsOptions.getTemperature(), esvSystem, x, world.size());
284
285 long totalSteps = dynamicsOptions.getNumSteps();
286 int nSteps = repEx.getReplicaSteps();
287 int exchangeCycles = (int) (totalSteps / nSteps);
288
289 if (exchangeCycles <= 0) {
290 exchangeCycles = 1;
291 }
292
293 pHReplicaExchange.sample(exchangeCycles, nSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(),
294 dynamicsOptions.getDt() * titrSteps - 1, initDynamics);
295
296 if (sort) {
297 sortMyArc(structureFile, size, pHReplicaExchange.getPhScale()[world.rank()], world.rank());
298 }
299 } catch (Exception e) {
300 logger.severe(" Error during pH replica exchange: " + e.getMessage());
301 return this;
302 }
303
304 } else {
305
306 molecularDynamics.dynamic(dynamicsOptions.getSteps(), dynamicsOptions.getDt(),
307 dynamicsOptions.getReport(), dynamicsOptions.getWrite(), dynamicsOptions.getTemperature(), true, dyn);
308 esvSystem.writeLambdaHistogram(true);
309 }
310 } else {
311
312
313
314 MolecularDynamicsOpenMM molecularDynamicsOpenMM = (MolecularDynamicsOpenMM) molecularDynamics;
315
316
317 molecularDynamics = dynamicsOptions.getDynamics(writeOutOptions, potential,
318 activeAssembly, algorithmListener, MDEngine.FFX);
319 molecularDynamics.attachExtendedSystem(esvSystem, titrReport);
320
321 if (repEx.getRepEx()) {
322 Comm world = Comm.world();
323 int size = world.size();
324
325 logger.info("\n Running replica exchange molecular dynamics on " + filename);
326 int rank = (size > 1) ? world.rank() : 0;
327 logger.info(" Rank: " + rank);
328
329 String pHWindows = Arrays.toString(PhReplicaExchange.setEvenSpacePhLadder(pHGap, pH, size));
330 pHWindows = pHWindows.replace("[", "")
331 .replace("]", "")
332 .replace(",", " ");
333 CompositeConfiguration properties = activeAssembly.getProperties();
334 pHWindows = properties.getString("pH.Windows", pHWindows);
335 String[] temp = pHWindows.split(" +");
336 if (temp.length != size) {
337 logger.severe("pHLadder specified in properties/key file has " +
338 "incorrect number of windows given world.size()");
339 }
340 double[] pHLadder = new double[size];
341 for (int i = 0; i < temp.length; i++) {
342 pHLadder[i] = Double.parseDouble(temp[i]);
343 }
344
345 File structureFile = new File(filename);
346 File rankDirectory = new File(structureFile.getParent() + File.separator + Integer.toString(rank));
347
348 final String newMolAssemblyFile = rankDirectory.getPath() + File.separator + structureFile.getName();
349 logger.info(" Set activeAssembly filename: " + newMolAssemblyFile);
350 activeAssembly.setFile(new File(newMolAssemblyFile));
351
352 try {
353 PhReplicaExchange pHReplicaExchange = new PhReplicaExchange(molecularDynamics, structureFile, pH, pHLadder,
354 dynamicsOptions.getTemperature(), esvSystem, x, molecularDynamicsOpenMM, potential, world.size());
355
356 pHReplicaExchange.
357 sample(cycles, titrSteps, coordSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(),
358 dynamicsOptions.getDt() * titrSteps, initDynamics);
359
360 if (sort) {
361 sortMyArc(structureFile, world.size(), pHReplicaExchange.getPhScale()[world.rank()], world.rank());
362 }
363 } catch (Exception e) {
364 logger.severe("Error during pH replica exchange with OpenMM: " + e.getMessage());
365 return this;
366 }
367
368 } else {
369 for (int i = 0; i < cycles; i++) {
370
371 molecularDynamics.setCoordinates(x);
372 double forceWriteInterval = titrSteps * 0.001;
373 molecularDynamics.dynamic(titrSteps, dynamicsOptions.getDt(), titrReport, forceWriteInterval,
374 dynamicsOptions.getTemperature(), true, dyn);
375 x = molecularDynamics.getCoordinates();
376 esvSystem.writeLambdaHistogram(true);
377
378
379 potential.energy(x);
380 molecularDynamicsOpenMM.setCoordinates(x);
381 if (coordSteps != 0) {
382 molecularDynamicsOpenMM.dynamic(coordSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(),
383 dynamicsOptions.getWrite(), dynamicsOptions.getTemperature(), true, dyn);
384 x = molecularDynamicsOpenMM.getCoordinates();
385 }
386 }
387 }
388 }
389
390 if (printRatio) {
391 esvSystem.printProtonationRatios();
392 }
393
394 return this;
395 }
396
397
398
399
400 @Override
401 public List<Potential> getPotentials() {
402 List<Potential> potentials;
403 if (potential == null) {
404 potentials = Collections.emptyList();
405 } else {
406 potentials = Collections.singletonList(potential);
407 }
408 return potentials;
409 }
410
411
412
413
414 public static void sortMyArc(File structureFile, int nReplicas, double pH, int myRank) {
415 logger.info(" Sorting archive for rank " + myRank);
416 String parent = structureFile.getParent();
417 String arcName = FilenameUtils.removeExtension(structureFile.getName()) + ".arc";
418 BufferedReader[] bufferedReaders = new BufferedReader[nReplicas];
419 File output = new File(parent + File.separator + myRank + File.separator + arcName + "_sorted");
420
421 try (BufferedWriter out = new BufferedWriter(new FileWriter(output))) {
422
423 File temp = new File(parent + File.separator + 0 + File.separator + arcName);
424 int snapLength = 0;
425 int totalLines = 0;
426 try (BufferedReader brTemp = new BufferedReader(new FileReader(temp))) {
427 String data = brTemp.readLine();
428 boolean startSnap = false;
429 while (data != null) {
430 if (data.contains("pH:")) {
431 startSnap = !startSnap;
432 if (!startSnap) {
433 break;
434 }
435 }
436 data = brTemp.readLine();
437 snapLength++;
438 }
439 totalLines = snapLength;
440 while (data != null) {
441 totalLines++;
442 data = brTemp.readLine();
443 }
444 }
445 int numSnaps = totalLines / snapLength;
446
447
448 for (int i = 0; i < nReplicas; i++) {
449 File file = new File(parent + File.separator + i + File.separator + arcName);
450 bufferedReaders[i] = new BufferedReader(new FileReader(file));
451 }
452
453 try {
454
455 for (int i = 0; i < numSnaps; i++) {
456 for (int j = 0; j < nReplicas; j++) {
457
458 String data = bufferedReaders[j].readLine();
459 while (data != null) {
460 if (data.contains("pH:")) {
461 break;
462 }
463 data = bufferedReaders[j].readLine();
464 }
465
466 String[] tokens = data.split(" +");
467 double snapPh = Double.parseDouble(tokens[tokens.length - 3]);
468
469
470 for (int k = 0; k < snapLength - 1; k++) {
471 if (snapPh == pH) {
472 out.write(data + "\n");
473 }
474 data = bufferedReaders[j].readLine();
475 }
476 if (snapPh == pH) {
477 out.write("\n");
478 }
479 }
480 }
481 } catch (IOException e) {
482 e.printStackTrace();
483 } finally {
484
485 for (int i = 0; i < nReplicas; i++) {
486 if (bufferedReaders[i] != null) {
487 try {
488 bufferedReaders[i].close();
489 } catch (IOException e) {
490 e.printStackTrace();
491 }
492 }
493 }
494 }
495 } catch (IOException e) {
496 e.printStackTrace();
497 }
498 }
499 }