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.cli;
39
40 import ffx.algorithms.AlgorithmListener;
41 import ffx.algorithms.dynamics.MDEngine;
42 import ffx.algorithms.dynamics.MolecularDynamics;
43 import ffx.algorithms.thermodynamics.HistogramData;
44 import ffx.algorithms.thermodynamics.LambdaData;
45 import ffx.algorithms.thermodynamics.MonteCarloOST;
46 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
47 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering.OptimizationParameters;
48 import ffx.crystal.CrystalPotential;
49 import ffx.potential.MolecularAssembly;
50 import ffx.potential.bonded.LambdaInterface;
51 import ffx.potential.cli.WriteoutOptions;
52 import org.apache.commons.configuration2.CompositeConfiguration;
53 import org.apache.commons.configuration2.Configuration;
54 import picocli.CommandLine.ArgGroup;
55 import picocli.CommandLine.Option;
56
57 import javax.annotation.Nullable;
58 import java.io.File;
59 import java.util.logging.Logger;
60
61
62
63
64
65
66
67
68
69
70 public class OSTOptions {
71
72 private static final Logger logger = Logger.getLogger(OSTOptions.class.getName());
73
74
75
76
77 @ArgGroup(heading = "%n Orthogonal Space Tempering Options%n", validate = false)
78 private final OSTOptionGroup group = new OSTOptionGroup();
79
80
81
82
83 @ArgGroup(heading = "%n Monte Carlo Orthogonal Space Tempering Options%n", validate = false)
84 private final MCOSTOptionGroup mcGroup = new MCOSTOptionGroup();
85
86
87
88
89
90
91
92
93
94
95
96 public CrystalPotential applyAllOSTOptions(OrthogonalSpaceTempering orthogonalSpaceTempering,
97 MolecularAssembly firstAssembly, DynamicsOptions dynamicsOptions,
98 BarostatOptions barostatOptions) {
99 if (dynamicsOptions.getOptimize()) {
100 OptimizationParameters opt = orthogonalSpaceTempering.getOptimizationParameters();
101 opt.setOptimization(true, firstAssembly);
102 }
103 return barostatOptions.checkNPT(firstAssembly, orthogonalSpaceTempering);
104 }
105
106
107
108
109
110
111
112
113 public void beginMCOST(MonteCarloOST monteCarloOST, DynamicsOptions dynamicsOptions,
114 ThermodynamicsOptions thermodynamicsOptions) {
115 long nEquil = thermodynamicsOptions.getEquilSteps();
116
117 if (nEquil > 0) {
118 logger.info("\n Beginning MC-OST equilibration.");
119 monteCarloOST.setTotalSteps(nEquil);
120 if (mcGroup.twoStep) {
121 monteCarloOST.sampleTwoStep();
122 } else {
123 monteCarloOST.sampleOneStep();
124 }
125 monteCarloOST.setEquilibration(false);
126 logger.info("\n Finished MC-OST equilibration.");
127 }
128
129 logger.info("\n Beginning MC-OST sampling.");
130 monteCarloOST.setLambdaStdDev(mcGroup.mcLambdaStdDev);
131 monteCarloOST.setTotalSteps(dynamicsOptions.getSteps());
132
133 if (mcGroup.twoStep) {
134 monteCarloOST.sampleTwoStep();
135 } else {
136 monteCarloOST.sampleOneStep();
137 }
138 }
139
140
141
142
143
144
145
146
147
148
149 public MolecularDynamics assembleMolecularDynamics(MolecularAssembly[] molecularAssemblies,
150 CrystalPotential crystalPotential, DynamicsOptions dynamicsOptions,
151 AlgorithmListener algorithmListener) {
152
153 MolecularAssembly firstTop = molecularAssemblies[0];
154
155 dynamicsOptions.init();
156
157 MolecularDynamics molDyn = MolecularDynamics.dynamicsFactory(firstTop, crystalPotential,
158 algorithmListener, dynamicsOptions.thermostat, dynamicsOptions.integrator, MDEngine.FFX);
159 for (int i = 1; i < molecularAssemblies.length; i++) {
160 molDyn.addAssembly(molecularAssemblies[i]);
161 }
162 molDyn.setRestartFrequency(dynamicsOptions.getCheckpoint());
163
164 return molDyn;
165 }
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182 public OrthogonalSpaceTempering constructOST(CrystalPotential crystalPotential,
183 @Nullable File lambdaRestartFile,
184 File histogramRestartFile, MolecularAssembly firstAssembly,
185 @Nullable Configuration addedProperties,
186 DynamicsOptions dynamicsOptions, ThermodynamicsOptions thermodynamicsOptions,
187 LambdaParticleOptions lambdaParticleOptions,
188 @Nullable AlgorithmListener algorithmListener, boolean async) {
189
190 LambdaInterface lambdaInterface = (LambdaInterface) crystalPotential;
191 CompositeConfiguration compositeConfiguration = new CompositeConfiguration(firstAssembly.getProperties());
192 if (addedProperties != null) {
193 compositeConfiguration.addConfiguration(addedProperties);
194 }
195
196
197 HistogramData histogramData = HistogramData.readHistogram(histogramRestartFile);
198
199 int histogramIndex = 0;
200 if (histogramData.wasHistogramRead()) {
201 logger.info("\n Read histogram restart from: " + histogramRestartFile);
202 logger.info(histogramData.toString());
203 } else {
204
205 histogramData.applyProperties(compositeConfiguration);
206
207 histogramData.setIndependentWalkers(group.independentWalkers);
208 histogramData.setWriteIndependent(group.independentWalkers);
209 histogramData.setAsynchronous(async);
210 histogramData.setTemperingFactor(getTemperingParameter(histogramIndex));
211 if (thresholdsSet()) {
212 histogramData.setTemperingOffset(getTemperingThreshold(histogramIndex));
213 }
214 histogramData.setMetaDynamics(group.metaDynamics);
215 histogramData.setBiasMag(getBiasMag(histogramIndex));
216 histogramData.setCountInterval(group.countInterval);
217 logger.info("\n OST Histogram Settings");
218 logger.info(histogramData.toString());
219 }
220
221
222 if (lambdaRestartFile == null) {
223 String filename = histogramRestartFile.toString().replaceFirst("\\.his$", ".lam");
224 lambdaRestartFile = new File(filename);
225 }
226 LambdaData lambdaData = LambdaData.readLambdaData(lambdaRestartFile);
227
228 if (lambdaData.wasLambdaRead()) {
229 logger.info(" Read lambda restart from: " + lambdaRestartFile);
230 } else {
231 lambdaData.setHistogramIndex(histogramIndex);
232 if (thermodynamicsOptions.getResetNumSteps()) {
233 lambdaData.setStepsTaken(0);
234 }
235 }
236
237 OrthogonalSpaceTempering orthogonalSpaceTempering = new OrthogonalSpaceTempering(lambdaInterface,
238 crystalPotential, histogramData, lambdaData, compositeConfiguration,
239 dynamicsOptions, lambdaParticleOptions, algorithmListener);
240 orthogonalSpaceTempering.setHardWallConstraint(mcGroup.mcHardWall);
241
242
243 return orthogonalSpaceTempering;
244 }
245
246
247
248
249
250
251
252
253
254
255
256
257
258 public void beginMDOST(OrthogonalSpaceTempering orthogonalSpaceTempering,
259 MolecularAssembly[] molecularAssemblies, CrystalPotential crystalPotential,
260 DynamicsOptions dynamicsOptions, WriteoutOptions writeoutOptions,
261 ThermodynamicsOptions thermodynamicsOptions, File dynFile,
262 AlgorithmListener algorithmListener) {
263
264 dynamicsOptions.init();
265
266 MolecularDynamics molDyn = assembleMolecularDynamics(molecularAssemblies, crystalPotential,
267 dynamicsOptions, algorithmListener);
268
269 boolean initVelocities = true;
270 long nSteps = dynamicsOptions.getSteps();
271
272 long nEquil = thermodynamicsOptions.getEquilSteps();
273 if (nEquil > 0) {
274 logger.info("\n Beginning equilibration");
275 orthogonalSpaceTempering.setPropagateLambda(false);
276 runDynamics(molDyn, nEquil, dynamicsOptions, writeoutOptions, true, dynFile);
277 logger.info(" Beginning OST sampling");
278 orthogonalSpaceTempering.setPropagateLambda(true);
279 } else {
280 logger.info(" Beginning OST sampling without equilibration");
281 if (!thermodynamicsOptions.getResetNumSteps()) {
282 long nEnergyCount = orthogonalSpaceTempering.getEnergyCount();
283 if (nEnergyCount > 0) {
284 nSteps -= nEnergyCount;
285 logger.info(String.format(" Lambda file: %12d steps picked up, now sampling %12d steps",
286 nEnergyCount, nSteps));
287 initVelocities = false;
288 }
289 }
290 }
291 if (nSteps > 0) {
292 runDynamics(molDyn, nSteps, dynamicsOptions, writeoutOptions, initVelocities, dynFile);
293 } else {
294 logger.info(" No steps remaining for this process!");
295 }
296 }
297
298
299
300
301
302
303 public boolean getIndependentWalkers() {
304 return group.independentWalkers;
305 }
306
307
308
309
310
311
312 public boolean isMonteCarlo() {
313 return mcGroup.monteCarlo;
314 }
315
316 public void setMonteCarlo(boolean monteCarlo) {
317 mcGroup.monteCarlo = monteCarlo;
318 }
319
320
321
322
323
324
325 public boolean isTwoStep() {
326 return mcGroup.twoStep;
327 }
328
329 public void setTwoStep(boolean twoStep) {
330 mcGroup.twoStep = twoStep;
331 }
332
333
334
335
336
337
338
339
340
341
342
343
344
345 public MonteCarloOST setupMCOST(OrthogonalSpaceTempering orthogonalSpaceTempering,
346 MolecularAssembly[] molecularAssemblies, CrystalPotential crystalPotential,
347 DynamicsOptions dynamicsOptions,
348 ThermodynamicsOptions thermodynamicsOptions, boolean verbose,
349 File dynRestart, AlgorithmListener algorithmListener) {
350 dynamicsOptions.init();
351
352 MonteCarloOST monteCarloOST = new MonteCarloOST(crystalPotential, orthogonalSpaceTempering,
353 molecularAssemblies[0], molecularAssemblies[0].getProperties(),
354 algorithmListener, dynamicsOptions, verbose, mcGroup.mcMDSteps, dynRestart);
355
356 MolecularDynamics md = monteCarloOST.getMD();
357 for (int i = 1; i < molecularAssemblies.length; i++) {
358 md.addAssembly(molecularAssemblies[i]);
359 }
360
361 long nEquil = thermodynamicsOptions.getEquilSteps();
362 if (nEquil > 0) {
363 monteCarloOST.setEquilibration(true);
364 }
365 return monteCarloOST;
366 }
367
368 private boolean thresholdsSet() {
369 return group.temperingThreshold.length != 1 || group.temperingThreshold[0] >= 0;
370 }
371
372
373
374
375
376
377
378 private double getBiasMag(int i) {
379 return group.biasMag.length > 1 ? group.biasMag[i] : group.biasMag[0];
380 }
381
382 private void runDynamics(MolecularDynamics molecularDynamics, long numSteps,
383 DynamicsOptions dynamicsOptions, WriteoutOptions writeoutOptions, boolean initVelocities,
384 File dyn) {
385 molecularDynamics.dynamic(numSteps, dynamicsOptions.getDt(), dynamicsOptions.getReport(),
386 dynamicsOptions.getWrite(), dynamicsOptions.getTemperature(), initVelocities,
387 writeoutOptions.getFileType(), dynamicsOptions.getCheckpoint(), dyn);
388 }
389
390
391
392
393
394
395
396 private double getTemperingThreshold(int i) {
397 return group.temperingThreshold.length > 1 ? group.temperingThreshold[i]
398 : group.temperingThreshold[0];
399 }
400
401
402
403
404
405
406
407 private double getTemperingParameter(int i) {
408 return group.temperingRate.length > 1 ? group.temperingRate[i] : group.temperingRate[0];
409 }
410
411
412
413
414
415
416 public int getCountInterval() {
417 return group.countInterval;
418 }
419
420 public void setCountInterval(int countInterval) {
421 group.countInterval = countInterval;
422 }
423
424
425
426
427
428
429 public double[] getBiasMag() {
430 return group.biasMag;
431 }
432
433 public void setBiasMag(double[] biasMag) {
434 group.biasMag = biasMag;
435 }
436
437
438
439
440
441
442 public boolean isIndependentWalkers() {
443 return group.independentWalkers;
444 }
445
446 public void setIndependentWalkers(boolean independentWalkers) {
447 group.independentWalkers = independentWalkers;
448 }
449
450
451
452
453
454
455 public boolean isMetaDynamics() {
456 return group.metaDynamics;
457 }
458
459 public void setMetaDynamics(boolean metaDynamics) {
460 group.metaDynamics = metaDynamics;
461 }
462
463
464
465
466
467
468 public double[] getTemperingRate() {
469 return group.temperingRate;
470 }
471
472 public void setTemperingRate(double[] temperingRate) {
473 group.temperingRate = temperingRate;
474 }
475
476
477
478
479
480
481 public double[] getTemperingThreshold() {
482 return group.temperingThreshold;
483 }
484
485 public void setTemperingThreshold(double[] temperingThreshold) {
486 group.temperingThreshold = temperingThreshold;
487 }
488
489
490
491
492
493
494
495 public boolean isMcHardWall() {
496 return mcGroup.mcHardWall;
497 }
498
499 public void setMcHardWall(boolean mcHardWall) {
500 mcGroup.mcHardWall = mcHardWall;
501 }
502
503
504
505
506
507
508 public int getMcMDSteps() {
509 return mcGroup.mcMDSteps;
510 }
511
512 public void setMcMDSteps(int mcMDSteps) {
513 mcGroup.mcMDSteps = mcMDSteps;
514 }
515
516
517
518
519
520
521 public double getMcLambdaStdDev() {
522 return mcGroup.mcLambdaStdDev;
523 }
524
525 public void setMcLambdaStdDev(double mcLambdaStdDev) {
526 mcGroup.mcLambdaStdDev = mcLambdaStdDev;
527 }
528
529
530
531
532
533
534 public double getLambdaWriteOut() {
535 return group.lambdaWriteOut;
536 }
537
538 public void setLambdaWriteOut(double lambdaWriteOut) {
539 group.lambdaWriteOut = lambdaWriteOut;
540 }
541
542
543
544
545 private static class OSTOptionGroup {
546
547
548
549
550 @Option(names = {"-C", "--count"}, paramLabel = "10", defaultValue = "10",
551 description = "Time steps between MD Orthogonal Space counts.")
552 private int countInterval = 10;
553
554
555
556
557 @Option(names = {"--bM", "--biasMag"}, paramLabel = "0.05", defaultValue = "0.05", split = ",",
558 description = "Orthogonal Space Gaussian bias magnitude (kcal/mol); RepEx OST uses a comma-separated list.")
559 private double[] biasMag = {0.05};
560
561
562
563
564 @Option(names = {"--iW", "--independentWalkers"}, defaultValue = "false",
565 description = "Enforces that each walker maintains their own histogram.")
566 private boolean independentWalkers = false;
567
568
569
570
571 @Option(names = {"--meta", "--metaDynamics"}, defaultValue = "false",
572 description = "Use a 1D metadynamics style bias.")
573 private boolean metaDynamics = false;
574
575
576
577
578 @Option(names = {"--tp", "--temperingRate"}, paramLabel = "4.0", defaultValue = "4.0", split = ",",
579 description = "Tempering rate parameter in multiples of kBT; RepEx OST uses a comma-separated list.")
580 private double[] temperingRate = {4.0};
581
582
583
584
585 @Option(names = {"--tth", "--temperingThreshold"}, paramLabel = "20*bias", defaultValue = "-1", split = ",",
586 description = "Tempering threshold in kcal/mol; RepEx OST uses a comma-separated list.")
587 private double[] temperingThreshold = {-1.0};
588
589
590
591
592
593 @Option(names = {"--lw", "--lambdaWriteOut"}, paramLabel = "0.0", defaultValue = "0.0",
594 description = "Only write out snapshots if lambda is greater than the value specified.")
595 private double lambdaWriteOut = 0.0;
596 }
597
598
599
600
601 private static class MCOSTOptionGroup {
602
603
604
605
606 @Option(names = {"--mc", "--monteCarlo"}, defaultValue = "false",
607 description = "Specify use of Monte Carlo OST")
608 private boolean monteCarlo = false;
609
610
611
612
613
614 @Option(names = {"--mcHW", "--mcHardWall"}, defaultValue = "false",
615 description = "Monte Carlo OST hard wall constraint.")
616 private boolean mcHardWall = false;
617
618
619
620
621 @Option(names = {"--mcMD", "--mcMDSteps"}, paramLabel = "100", defaultValue = "100",
622 description = "Number of dynamics steps to take for each MD trajectory for Monte Carlo OST")
623 private int mcMDSteps = 100;
624
625
626
627
628 @Option(names = {"--mcL", "--mcLambdaStdDev"}, paramLabel = "0.01", defaultValue = "0.01",
629 description = "Standard deviation for lambda move.")
630 private double mcLambdaStdDev = 0.01;
631
632
633
634
635 @Option(names = {"--ts", "--twoStep"}, defaultValue = "false",
636 description = "MC Orthogonal Space sampling using separate lambda and MD moves.")
637 private boolean twoStep = false;
638 }
639
640 }