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.dynamics;
39
40 import static ffx.utilities.Constants.FSEC_TO_PSEC;
41 import static ffx.utilities.Constants.KCAL_TO_GRAM_ANG2_PER_PS2;
42 import static ffx.utilities.Constants.NS2SEC;
43 import static java.lang.String.format;
44 import static java.util.Arrays.fill;
45
46 import edu.rit.pj.Comm;
47 import ffx.algorithms.AlgorithmListener;
48 import ffx.algorithms.Terminatable;
49 import ffx.algorithms.dynamics.integrators.BetterBeeman;
50 import ffx.algorithms.dynamics.integrators.Integrator;
51 import ffx.algorithms.dynamics.integrators.IntegratorEnum;
52 import ffx.algorithms.dynamics.integrators.Respa;
53 import ffx.algorithms.dynamics.integrators.Stochastic;
54 import ffx.algorithms.dynamics.integrators.VelocityVerlet;
55 import ffx.algorithms.dynamics.thermostats.Adiabatic;
56 import ffx.algorithms.dynamics.thermostats.Berendsen;
57 import ffx.algorithms.dynamics.thermostats.Bussi;
58 import ffx.algorithms.dynamics.thermostats.Thermostat;
59 import ffx.algorithms.dynamics.thermostats.ThermostatEnum;
60 import ffx.algorithms.thermodynamics.OrthogonalSpaceTempering;
61 import ffx.crystal.Crystal;
62 import ffx.numerics.Constraint;
63 import ffx.numerics.Potential;
64 import ffx.numerics.math.RunningStatistics;
65 import ffx.potential.openmm.OpenMMEnergy;
66 import ffx.potential.MolecularAssembly;
67 import ffx.potential.SystemState;
68 import ffx.potential.UnmodifiableState;
69 import ffx.potential.bonded.Atom;
70 import ffx.potential.bonded.LambdaInterface;
71 import ffx.potential.extended.ExtendedSystem;
72 import ffx.potential.parameters.ForceField;
73 import ffx.potential.parsers.DYNFilter;
74 import ffx.potential.parsers.PDBFilter;
75 import ffx.potential.parsers.XPHFilter;
76 import ffx.potential.parsers.XYZFilter;
77 import ffx.utilities.FileUtils;
78 import java.io.File;
79 import java.util.ArrayList;
80 import java.util.Arrays;
81 import java.util.EnumSet;
82 import java.util.List;
83 import java.util.logging.Level;
84 import java.util.logging.Logger;
85 import org.apache.commons.configuration2.CompositeConfiguration;
86 import org.apache.commons.io.FilenameUtils;
87
88
89
90
91
92
93
94 public class MolecularDynamics implements Runnable, Terminatable {
95
96 private static final Logger logger = Logger.getLogger(MolecularDynamics.class.getName());
97
98
99
100
101 private static final double DEFAULT_LOG_INTERVAL = 0.25;
102
103
104
105 private static final double DEFAULT_RESTART_INTERVAL = 1.0;
106
107
108
109 private static final double DEFAULT_TRAJECTORY_INTERVAL = 10.0;
110
111 private static final int DEFAULT_DYNAMICS_SLEEP_TIME = 25000;
112
113
114
115
116
117
118 protected MolecularAssembly[] molecularAssembly;
119
120 private final Potential potential;
121
122 protected final AlgorithmListener algorithmListener;
123
124 private final Integrator integrator;
125
126 private Thermostat thermostat;
127
128 boolean constantPressure = false;
129
130
131
132 Barostat barostat = null;
133
134
135
136 protected final SystemState state;
137
138
139
140 protected UnmodifiableState initialState;
141
142
143
144
145 private UnmodifiableState storedState;
146
147
148
149
150 private final RunningStatistics temperatureStats = new RunningStatistics();
151
152
153
154 private final RunningStatistics potentialEnergyStats = new RunningStatistics();
155
156
157
158 private final RunningStatistics kineticEnergyStats = new RunningStatistics();
159
160
161
162 private final RunningStatistics totalEnergyStats = new RunningStatistics();
163
164
165
166
167
168
169 protected double dt = 0.001;
170
171
172
173 private long nSteps = 1000;
174
175
176
177 double targetTemperature = 298.15;
178
179
180
181 protected boolean initVelocities = true;
182
183
184
185
186 protected boolean automaticWriteouts = true;
187
188
189
190 protected double totalSimTime = 0.0;
191
192
193 protected double restartInterval = DEFAULT_RESTART_INTERVAL;
194
195 protected int restartFrequency;
196
197 private double trajectoryInterval = DEFAULT_TRAJECTORY_INTERVAL;
198
199 protected int trajectoryFrequency;
200
201 private double logInterval = DEFAULT_LOG_INTERVAL;
202
203
204
205 protected int logFrequency;
206
207
208 protected String fileType = "XYZ";
209
210 boolean saveSnapshotAsPDB = true;
211
212 File restartFile = null;
213
214 boolean loadRestart = false;
215
216 DYNFilter dynFilter = null;
217
218 private File fallbackDynFile;
219
220
221
222
223 Level basicLogging = Level.INFO;
224
225 private MDVerbosity verbosityLevel = MDVerbosity.VERBOSE;
226
227
228 boolean initialized = false;
229
230 private boolean terminate = false;
231
232 protected boolean done;
233
234
235 private ExtendedSystem esvSystem;
236
237 private Stochastic esvIntegrator;
238
239 private Adiabatic esvThermostat;
240
241 private int printEsvFrequency = -1;
242
243
244
245
246
247
248
249
250
251
252
253 public MolecularDynamics(MolecularAssembly assembly, Potential potentialEnergy,
254 AlgorithmListener listener, ThermostatEnum requestedThermostat,
255 IntegratorEnum requestedIntegrator) {
256 this.molecularAssembly = new MolecularAssembly[1];
257 this.molecularAssembly[0] = assembly;
258 this.algorithmListener = listener;
259 this.potential = potentialEnergy;
260 if (potentialEnergy instanceof Barostat) {
261 constantPressure = true;
262 barostat = (Barostat) potentialEnergy;
263 }
264
265 int numberOfVariables = potentialEnergy.getNumberOfVariables();
266 state = new SystemState(numberOfVariables);
267 state.setMass(potential.getMass());
268
269 CompositeConfiguration properties = molecularAssembly[0].getProperties();
270
271 if (requestedIntegrator == null) {
272 String integrate = properties.getString("integrate", "VERLET").trim();
273 try {
274 integrate = integrate.toUpperCase().replace("-", "_");
275 requestedIntegrator = IntegratorEnum.valueOf(integrate.toUpperCase());
276 } catch (Exception e) {
277 requestedIntegrator = IntegratorEnum.VERLET;
278 }
279 }
280
281 boolean oMMLogging = potential instanceof OpenMMEnergy;
282 List<Constraint> constraints = potentialEnergy.getConstraints();
283
284
285 integrator = switch (requestedIntegrator) {
286 case RESPA, MTS -> {
287 Respa respa = new Respa(state);
288 int in = molecularAssembly[0].getProperties().getInt("respa-dt", 4);
289 if (in < 2) {
290 in = 2;
291 }
292 if (!oMMLogging) {
293 respa.setInnerTimeSteps(in);
294 }
295 logger.log(Level.FINE, format(" Created a RESPA integrator with %d inner time steps.", in));
296 yield respa;
297 }
298 case STOCHASTIC, LANGEVIN -> {
299 double friction = properties.getDouble("friction", 91.0);
300 logger.log(Level.FINE, format(" Friction set at %.3f collisions/picosecond", friction));
301 Stochastic stochastic = new Stochastic(friction, state);
302 if (properties.containsKey("randomseed")) {
303 stochastic.setRandomSeed(properties.getInt("randomseed", 0));
304 }
305
306
307
308 requestedThermostat = ThermostatEnum.ADIABATIC;
309 yield stochastic;
310 }
311 case BEEMAN -> new BetterBeeman(state);
312 default -> new VelocityVerlet(state);
313 };
314
315 integrator.addConstraints(constraints);
316
317
318 if (requestedThermostat == null) {
319 String thermo = properties.getString("thermostat", "Berendsen").trim();
320 try {
321 thermo = thermo.toUpperCase().replace("-", "_");
322 requestedThermostat = ThermostatEnum.valueOf(thermo.toUpperCase());
323 } catch (Exception e) {
324 requestedThermostat = ThermostatEnum.BERENDSEN;
325 }
326 }
327
328
329 thermostat = switch (requestedThermostat) {
330 case BERENDSEN -> {
331 double tau = properties.getDouble("tau-temperature", 0.2);
332 yield new Berendsen(state, potentialEnergy.getVariableTypes(), targetTemperature, tau,
333 constraints);
334 }
335 case BUSSI -> {
336 double tau = properties.getDouble("tau-temperature", 0.2);
337 Bussi bussi = new Bussi(state, potentialEnergy.getVariableTypes(), targetTemperature, tau,
338 constraints);
339 if (properties.containsKey("randomseed")) {
340 bussi.setRandomSeed(properties.getInt("randomseed", 0));
341 }
342 yield bussi;
343 }
344 default -> new Adiabatic(state, potentialEnergy.getVariableTypes(), constraints);
345 };
346
347 if (properties.containsKey("randomseed")) {
348 thermostat.setRandomSeed(properties.getInt("randomseed", 0));
349 }
350
351
352 if (integrator instanceof Stochastic) {
353 thermostat.setRemoveCenterOfMassMotion(false);
354 }
355
356 done = true;
357 fallbackDynFile = defaultFallbackDyn(assembly);
358 }
359
360
361
362
363
364
365
366
367
368
369
370
371
372 public static MolecularDynamics dynamicsFactory(MolecularAssembly assembly,
373 Potential potentialEnergy, AlgorithmListener listener, ThermostatEnum requestedThermostat,
374 IntegratorEnum requestedIntegrator) {
375
376 return dynamicsFactory(assembly, potentialEnergy, listener, requestedThermostat,
377 requestedIntegrator, defaultEngine(assembly, potentialEnergy));
378 }
379
380
381
382
383
384
385
386
387
388
389
390
391
392 public static MolecularDynamics dynamicsFactory(MolecularAssembly assembly,
393 Potential potentialEnergy, AlgorithmListener listener, ThermostatEnum requestedThermostat,
394 IntegratorEnum requestedIntegrator, MDEngine engine) {
395
396 switch (engine) {
397 case OPENMM, OMM -> {
398
399
400
401
402 boolean ommLeaves = potentialEnergy.getUnderlyingPotentials().stream()
403 .anyMatch((Potential p) -> p instanceof OpenMMEnergy);
404 ommLeaves = ommLeaves || potentialEnergy instanceof OpenMMEnergy;
405 if (ommLeaves) {
406 return new MolecularDynamicsOpenMM(assembly, potentialEnergy, listener,
407 requestedThermostat, requestedIntegrator);
408 } else {
409 throw new IllegalArgumentException(format(
410 " Requested OpenMM engine %s, but at least one leaf of the potential %s is not an OpenMM force field!",
411 engine, potentialEnergy));
412 }
413 }
414 default -> {
415 return new MolecularDynamics(assembly, potentialEnergy, listener, requestedThermostat,
416 requestedIntegrator);
417 }
418 }
419 }
420
421 private static MDEngine defaultEngine(MolecularAssembly molecularAssembly,
422 Potential potentialEnergy) {
423 CompositeConfiguration properties = molecularAssembly.getProperties();
424 String mdEngine = properties.getString("MD-engine");
425 if (mdEngine != null) {
426 if (mdEngine.equalsIgnoreCase("OMM")) {
427 logger.info(" Creating OpenMM Dynamics Object");
428 return MDEngine.OPENMM;
429 } else {
430 logger.info(" Creating FFX Dynamics Object");
431 return MDEngine.FFX;
432 }
433 } else {
434
435 boolean ommLeaves = potentialEnergy.getUnderlyingPotentials().stream()
436 .anyMatch((Potential p) -> p instanceof OpenMMEnergy);
437 ommLeaves = ommLeaves || potentialEnergy instanceof OpenMMEnergy;
438 if (ommLeaves) {
439 return MDEngine.OPENMM;
440 } else {
441 return MDEngine.FFX;
442 }
443 }
444 }
445
446
447
448
449
450
451
452 private static File defaultFallbackDyn(MolecularAssembly assembly) {
453 String firstFileName = FilenameUtils.removeExtension(assembly.getFile().getAbsolutePath());
454 return new File(firstFileName + ".dyn");
455 }
456
457
458
459
460
461
462
463 public void addAssembly(MolecularAssembly assembly) {
464 molecularAssembly = Arrays.copyOf(molecularAssembly, molecularAssembly.length + 1);
465 molecularAssembly[molecularAssembly.length - 1] = assembly;
466 }
467
468
469
470
471
472
473 public void attachExtendedSystem(ExtendedSystem system, double reportFreq) {
474 if (esvSystem != null) {
475 logger.warning("An ExtendedSystem is already attached to this MD!");
476 }
477 esvSystem = system;
478 SystemState esvState = esvSystem.getState();
479 this.esvIntegrator = new Stochastic(esvSystem.getThetaFriction(), esvState);
480 if(!esvSystem.getConstraints().isEmpty()){
481 esvIntegrator.addConstraints(esvSystem.getConstraints());
482 }
483 this.esvThermostat = new Adiabatic(esvState, potential.getVariableTypes());
484 printEsvFrequency = intervalToFreq(reportFreq, "Reporting (logging) interval");
485 logger.info(
486 format(" Attached extended system (%s) to molecular dynamics.", esvSystem.toString()));
487 logger.info(format(" Extended System Theta Friction: %f", esvSystem.getThetaFriction()));
488 logger.info(format(" Extended System Theta Mass: %f", esvSystem.getThetaMass()));
489 logger.info(format(" Extended System Lambda Print Frequency: %d (fsec)", printEsvFrequency));
490 }
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506 public void dynamic(final long nSteps, final double timeStep, final double loggingInterval,
507 final double trajectoryInterval, final double temperature, final boolean initVelocities,
508 String fileType, double restartInterval, final File dyn) {
509 this.fileType = fileType;
510 setRestartFrequency(restartInterval);
511 dynamic(nSteps, timeStep, loggingInterval, trajectoryInterval, temperature, initVelocities, dyn);
512 }
513
514
515
516
517
518
519
520
521
522
523
524
525
526 public void dynamic(final long nSteps, final double timeStep, final double loggingInterval,
527 final double trajectoryInterval, final double temperature, final boolean initVelocities,
528 final File dyn) {
529
530
531 if (!done) {
532 logger.warning(" Programming error - a thread invoked dynamic when it was already running.");
533 return;
534 }
535
536 init(nSteps, timeStep, loggingInterval, trajectoryInterval, fileType, restartInterval,
537 temperature, initVelocities, dyn);
538
539 Thread dynamicThread = new Thread(this);
540 dynamicThread.start();
541 synchronized (this) {
542 try {
543 while (dynamicThread.isAlive()) {
544 wait(0, DEFAULT_DYNAMICS_SLEEP_TIME);
545 }
546 } catch (InterruptedException e) {
547 String message = " Molecular dynamics interrupted.";
548 logger.log(Level.WARNING, message, e);
549 }
550 }
551 if (!verbosityLevel.isQuiet()) {
552 logger.info(" Done with an MD round.");
553 }
554 }
555
556
557
558
559
560
561 public MolecularAssembly[] getAssemblyArray() {
562 return molecularAssembly.clone();
563 }
564
565
566
567
568
569
570 public File getDynFile() {
571 return restartFile;
572 }
573
574
575
576
577
578
579 public double getInitialKineticEnergy() {
580 return initialState.kineticEnergy();
581 }
582
583
584
585
586
587
588 public double getInitialPotentialEnergy() {
589 return initialState.potentialEnergy();
590 }
591
592
593
594
595
596
597 public double getInitialTemperature() {
598 return initialState.temperature();
599 }
600
601
602
603
604
605
606 public double getInitialTotalEnergy() {
607 return initialState.getTotalEnergy();
608 }
609
610
611
612
613
614
615 public int getIntervalSteps() {
616 return 1;
617 }
618
619
620
621
622
623
624 public void setIntervalSteps(int intervalSteps) {
625
626 }
627
628
629
630
631
632
633 public double getKineticEnergy() {
634 return state.getKineticEnergy();
635 }
636
637
638
639
640
641
642 public double getPotentialEnergy() {
643 return state.getPotentialEnergy();
644 }
645
646
647
648
649
650
651 public double getTemperature() {
652 return state.getTemperature();
653 }
654
655
656
657
658
659
660 public Thermostat getThermostat() {
661 return thermostat;
662 }
663
664
665
666
667
668
669 public void setThermostat(Thermostat thermostat) {
670 this.thermostat = thermostat;
671 }
672
673
674
675
676
677
678 public double getTimeStep() {
679 return dt;
680 }
681
682
683
684
685
686
687 private void setTimeStep(double step) {
688 dt = step * FSEC_TO_PSEC;
689
690 setRestartFrequency(restartInterval);
691 setLoggingFrequency(logInterval);
692 setTrajectoryFrequency(trajectoryInterval);
693 }
694
695
696
697
698
699
700 public double getTotalEnergy() {
701 return state.getTotalEnergy();
702 }
703
704 public MDVerbosity getVerbosityLevel() {
705 return verbosityLevel;
706 }
707
708 public void setVerbosityLevel(MDVerbosity level) {
709 verbosityLevel = level;
710 if (level == MDVerbosity.SILENT) {
711 basicLogging = Level.FINE;
712 } else {
713 basicLogging = Level.INFO;
714 }
715 }
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731 public void init(final long nSteps, final double timeStep, final double loggingInterval,
732 final double trajectoryInterval, final String fileType, final double restartInterval,
733 final double temperature, final boolean initVelocities, final File dyn) {
734
735
736 if (!done) {
737 logger.warning(
738 " Programming error - attempt to modify parameters of a running MolecularDynamics instance.");
739 return;
740 }
741
742 if (integrator instanceof Stochastic) {
743 if (constantPressure) {
744 logger.log(basicLogging, "\n Stochastic dynamics in the NPT ensemble");
745 } else {
746 logger.log(basicLogging, "\n Stochastic dynamics in the NVT ensemble");
747 }
748 } else if (!(thermostat instanceof Adiabatic)) {
749 if (constantPressure) {
750 logger.log(basicLogging, "\n Molecular dynamics in the NPT ensemble");
751 } else {
752 logger.log(basicLogging, "\n Molecular dynamics in the NVT ensemble");
753 }
754 } else {
755 if (constantPressure) {
756 logger.severe("\n NPT Molecular dynamics requires a thermostat");
757 } else {
758 logger.log(basicLogging, "\n Molecular dynamics in the NVE ensemble");
759 }
760 }
761
762 this.nSteps = nSteps;
763 totalSimTime = 0.0;
764
765
766 setTimeStep(timeStep);
767
768
769 setLoggingFrequency(loggingInterval);
770 setTrajectoryFrequency(trajectoryInterval);
771 setRestartFrequency(restartInterval);
772
773 checkFrequency("Reporting ", logFrequency);
774 if (automaticWriteouts) {
775
776 checkFrequency("Trajectory Write-Out", trajectoryFrequency);
777 checkFrequency("Restart ", restartFrequency);
778 }
779
780
781 saveSnapshotAsPDB = true;
782 if (fileType.equalsIgnoreCase("XYZ") || fileType.equalsIgnoreCase("ARC")) {
783 saveSnapshotAsPDB = false;
784 } else if (!fileType.equalsIgnoreCase("PDB")) {
785 logger.warning("Snapshot file type unrecognized; saving snapshots as PDB.\n");
786 }
787
788 setArchiveFile();
789
790 this.restartFile = (dyn == null) ? fallbackDynFile : dyn;
791 loadRestart = restartFile.exists() && !initialized;
792
793 if (dynFilter == null) {
794 dynFilter = new DYNFilter(molecularAssembly[0].getName());
795 }
796
797 this.targetTemperature = temperature;
798 this.initVelocities = initVelocities;
799 done = false;
800
801 if (loadRestart) {
802 logger.info(" Continuing from " + restartFile.getAbsolutePath());
803 }
804
805 if (!verbosityLevel.isQuiet()) {
806 logger.info(format(" Integrator: %15s", integrator.toString()));
807 logger.info(format(" Thermostat: %15s", thermostat.toString()));
808 logger.info(format(" Number of steps: %8d", nSteps));
809 logger.info(format(" Time step: %8.3f (fsec)", timeStep));
810 logger.info(format(" Print interval: %8.3f (psec)", loggingInterval));
811 logger.info(format(" Save interval: %8.3f (psec)", trajectoryInterval));
812 if (molecularAssembly.length > 1) {
813 for (int i = 0; i < molecularAssembly.length; i++) {
814 File archiveFile = molecularAssembly[i].getArchiveFile();
815 logger.info(format(" Archive file %3d: %s", (i + 1), archiveFile.getAbsolutePath()));
816 }
817 } else {
818 logger.info(format(" Archive file: %s",
819 molecularAssembly[0].getArchiveFile().getAbsolutePath()));
820 }
821 logger.info(format(" Restart file: %s", restartFile.getAbsolutePath()));
822 }
823 }
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838 public void init(final long nSteps, final double timeStep, final double loggingInterval,
839 final double trajectoryInterval, final double temperature, final boolean initVelocities,
840 final File dyn) {
841 init(nSteps, timeStep, loggingInterval, trajectoryInterval, "XYZ", restartInterval, temperature,
842 initVelocities, dyn);
843 }
844
845
846 @Override
847 public void run() {
848 try {
849 preRunOps();
850 } catch (IllegalStateException ise) {
851 return;
852 }
853 initializeEnergies();
854 postInitEnergies();
855 mainLoop();
856 postRun();
857 }
858
859
860
861
862
863
864
865
866 public void setAutomaticWriteouts(boolean autoWriteout) {
867 this.automaticWriteouts = autoWriteout;
868 }
869
870
871
872
873
874
875 public void setFallbackDynFile(File fallback) {
876 this.fallbackDynFile = fallback;
877 }
878
879
880
881
882
883
884 public void setFileType(String fileType) {
885 this.fileType = fileType;
886 }
887
888
889
890
891
892
893
894 public void setObtainVelAcc(boolean obtainVA) {
895
896 }
897
898
899
900
901
902
903
904 public void setRestartFrequency(double restartInterval) throws IllegalArgumentException {
905 restartFrequency = intervalToFreq(restartInterval, "Restart interval");
906 this.restartInterval = restartInterval;
907 }
908
909
910
911
912
913
914 public void setArchiveFiles(File[] archiveFiles) {
915 logger.info(" Setting archive files:\n " + Arrays.toString(archiveFiles));
916 int n = molecularAssembly.length;
917 assert archiveFiles.length == n;
918 for (int i = 0; i < n; i++) {
919 molecularAssembly[i].setArchiveFile(archiveFiles[i]);
920 }
921 }
922
923
924 @Override
925 public void terminate() {
926 terminate = true;
927 while (!done) {
928 synchronized (this) {
929 try {
930 wait(1);
931 } catch (Exception e) {
932 logger.log(Level.WARNING, " Exception terminating dynamics.\n", e);
933 }
934 }
935 }
936 }
937
938
939
940
941
942
943
944
945
946
947 public EnumSet<MDWriteAction> writeFilesForStep(long step, boolean trySnapshot,
948 boolean tryRestart) {
949 return writeFilesForStep(step, trySnapshot, tryRestart, null);
950 }
951
952
953
954
955
956
957
958
959
960
961
962
963 public EnumSet<MDWriteAction> writeFilesForStep(long step, boolean trySnapshot, boolean tryRestart,
964 String[] extraLines) {
965 List<String> linesList =
966 (extraLines == null) ? new ArrayList<>() : new ArrayList<>(Arrays.asList(extraLines));
967
968 if (potential instanceof LambdaInterface) {
969 String lamString = format("Lambda: %.8f", ((LambdaInterface) potential).getLambda());
970 linesList.add(lamString);
971 }
972
973 String tempString = format("Temp: %.2f", thermostat.getTargetTemperature());
974 linesList.add(tempString);
975
976 String timeString = format(" Time: %7.3e", totalSimTime);
977 linesList.add(timeString);
978
979 if (esvSystem != null) {
980 String pHString = format("pH: %.2f", esvSystem.getConstantPh());
981 linesList.add(pHString);
982 }
983
984 Comm world = Comm.world();
985 if (world != null && world.size() > 1) {
986 String rankString = format("Rank: %d", world.rank());
987 linesList.add(rankString);
988 }
989
990 String[] allLines = new String[linesList.size()];
991 allLines = linesList.toArray(allLines);
992
993 EnumSet<MDWriteAction> written = EnumSet.noneOf(MDWriteAction.class);
994 if (step != 0) {
995
996 if (trySnapshot && trajectoryFrequency > 0 && step % trajectoryFrequency == 0) {
997
998
999 logger.info(format("\n Average Values for the Last %d Out of %d Dynamics Steps\n",
1000 trajectoryFrequency, step));
1001 logger.info(format(" Simulation Time %16.4f Picosecond", step * dt));
1002 logger.info(format(" Total Energy %16.4f Kcal/mole (+/-%9.4f)", totalEnergyStats.getMean(),
1003 totalEnergyStats.getStandardDeviation()));
1004 logger.info(format(" Potential Energy %16.4f Kcal/mole (+/-%9.4f)",
1005 potentialEnergyStats.getMean(), potentialEnergyStats.getStandardDeviation()));
1006 logger.info(format(" Kinetic Energy %16.4f Kcal/mole (+/-%9.4f)", kineticEnergyStats.getMean(),
1007 kineticEnergyStats.getStandardDeviation()));
1008 logger.info(format(" Temperature %16.4f Kelvin (+/-%9.4f)\n", temperatureStats.getMean(),
1009 temperatureStats.getStandardDeviation()));
1010 totalEnergyStats.reset();
1011 potentialEnergyStats.reset();
1012 kineticEnergyStats.reset();
1013 temperatureStats.reset();
1014
1015 if (esvSystem != null) {
1016 for (Atom atom : molecularAssembly[0].getAtomList()) {
1017 int atomIndex = atom.getIndex() - 1;
1018 atom.setOccupancy(esvSystem.getTitrationLambda(atomIndex));
1019 atom.setTempFactor(esvSystem.getTautomerLambda(atomIndex));
1020 }
1021 }
1022 appendSnapshot(allLines);
1023 written.add(MDWriteAction.SNAPSHOT);
1024 }
1025
1026
1027 if (tryRestart && restartFrequency > 0 && step % restartFrequency == 0) {
1028 writeRestart();
1029 written.add(MDWriteAction.RESTART);
1030 }
1031 }
1032 return written;
1033 }
1034
1035
1036 public void writeRestart() {
1037 String dynName = FileUtils.relativePathTo(restartFile).toString();
1038 double[] x = state.x();
1039 double[] v = state.v();
1040 double[] a = state.a();
1041 double[] aPrevious = state.aPrevious();
1042 if (dynFilter.writeDYN(restartFile, molecularAssembly[0].getCrystal(), x, v, a, aPrevious)) {
1043 logger.log(basicLogging, format(" Wrote dynamics restart to: %s.", dynName));
1044 } else {
1045 logger.log(basicLogging, format(" Writing dynamics restart failed: %s.", dynName));
1046 }
1047 if (esvSystem != null) {
1048 esvSystem.writeRestart();
1049 esvSystem.writeLambdaHistogram(false);
1050 }
1051 potential.writeAdditionalRestartInfo(true);
1052 }
1053
1054
1055
1056
1057 private void setArchiveFile() {
1058 for (MolecularAssembly assembly : molecularAssembly) {
1059 File file = assembly.getFile();
1060 String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
1061 File archiveFile = assembly.getArchiveFile();
1062 if (archiveFile == null) {
1063 archiveFile = new File(filename + ".arc");
1064 assembly.setArchiveFile(XYZFilter.version(archiveFile));
1065 }
1066 }
1067 }
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077 private int intervalToFreq(double interval, String describe) throws IllegalArgumentException {
1078 if (!Double.isFinite(interval) || interval <= 0) {
1079 throw new IllegalArgumentException(
1080 format(" %s must be " + "positive finite value in ps, was %10.4g", describe, interval));
1081 }
1082 if (interval >= dt) {
1083 return (int) (interval / dt);
1084 } else {
1085 logger.warning(format(" Specified %s of %.6f ps < time step %.6f ps; "
1086 + "interval is set to once per time step!", describe, interval, dt));
1087 return 1;
1088 }
1089 }
1090
1091
1092
1093
1094
1095
1096 private void setLoggingFrequency(double logInterval) {
1097 logFrequency = intervalToFreq(logInterval, "Reporting (logging) interval");
1098 this.logInterval = logInterval;
1099 }
1100
1101
1102
1103
1104
1105
1106 private void setTrajectoryFrequency(double snapshotInterval) {
1107 trajectoryFrequency = intervalToFreq(snapshotInterval, "Trajectory writeout interval");
1108 this.trajectoryInterval = snapshotInterval;
1109 }
1110
1111
1112 void preRunOps() throws IllegalStateException {
1113 done = false;
1114 terminate = false;
1115
1116
1117 thermostat.setTargetTemperature(targetTemperature);
1118 boolean quiet = verbosityLevel.isQuiet();
1119 thermostat.setQuiet(quiet);
1120 if (integrator instanceof Stochastic stochastic) {
1121 stochastic.setTemperature(targetTemperature);
1122 }
1123
1124
1125 integrator.setTimeStep(dt);
1126
1127 if (!initialized) {
1128
1129 if (loadRestart) {
1130 Crystal crystal = molecularAssembly[0].getCrystal();
1131 double[] x = state.x();
1132 double[] v = state.v();
1133 double[] a = state.a();
1134 double[] aPrevious = state.aPrevious();
1135 if (!dynFilter.readDYN(restartFile, crystal, x, v, a, aPrevious)) {
1136 String message = " Could not load the restart file - dynamics terminated.";
1137 logger.log(Level.WARNING, message);
1138 done = true;
1139 throw new IllegalStateException(message);
1140 } else {
1141 molecularAssembly[0].setCrystal(crystal);
1142 }
1143 } else {
1144
1145 potential.getCoordinates(state.x());
1146
1147 if (initVelocities) {
1148 thermostat.maxwell(targetTemperature);
1149 } else {
1150 fill(state.v(), 0.0);
1151 }
1152 }
1153 } else {
1154
1155
1156 if (initVelocities) {
1157 thermostat.maxwell(targetTemperature);
1158 }
1159 }
1160 }
1161
1162
1163 private void initializeEnergies() {
1164
1165 double[] x = state.x();
1166 double[] gradient = state.gradient();
1167
1168 boolean propagateLambda = true;
1169 if (potential instanceof OrthogonalSpaceTempering orthogonalSpaceTempering) {
1170 propagateLambda = orthogonalSpaceTempering.getPropagateLambda();
1171 orthogonalSpaceTempering.setPropagateLambda(false);
1172 }
1173
1174 if (esvSystem != null && potential instanceof OpenMMEnergy) {
1175 state.setPotentialEnergy(((OpenMMEnergy) potential).energyAndGradientFFX(x, gradient));
1176 } else {
1177 state.setPotentialEnergy(potential.energyAndGradient(x, gradient));
1178 }
1179
1180 if (potential instanceof OrthogonalSpaceTempering orthogonalSpaceTempering) {
1181 orthogonalSpaceTempering.setPropagateLambda(propagateLambda);
1182 }
1183
1184
1185 if (!loadRestart || initialized || integrator instanceof Respa) {
1186
1187 if (integrator instanceof Respa) {
1188 potential.setEnergyTermState(Potential.STATE.SLOW);
1189 potential.energyAndGradient(x, gradient);
1190 }
1191
1192 int numberOfVariables = state.getNumberOfVariables();
1193 double[] a = state.a();
1194 double[] mass = state.getMass();
1195 for (int i = 0; i < numberOfVariables; i++) {
1196 a[i] = -KCAL_TO_GRAM_ANG2_PER_PS2 * gradient[i] / mass[i];
1197 }
1198 state.copyAccelerationsToPrevious();
1199 }
1200
1201
1202 thermostat.computeKineticEnergy();
1203
1204
1205 initialState = new UnmodifiableState(state);
1206
1207
1208 temperatureStats.reset();
1209 potentialEnergyStats.reset();
1210 kineticEnergyStats.reset();
1211 totalEnergyStats.reset();
1212 }
1213
1214
1215 void postInitEnergies() {
1216 initialized = true;
1217
1218 logger.log(basicLogging,
1219 format("\n %8s %12s %12s %12s %8s %8s", "Time", "Kinetic", "Potential", "Total", "Temp", "CPU"));
1220 logger.log(basicLogging,
1221 format(" %8s %12s %12s %12s %8s %8s", "psec", "kcal/mol", "kcal/mol", "kcal/mol", "K", "sec"));
1222 logger.log(basicLogging, format(" %8s %12.4f %12.4f %12.4f %8.2f", "", state.getKineticEnergy(),
1223 state.getPotentialEnergy(), state.getTotalEnergy(), state.getTemperature()));
1224
1225
1226 storeState();
1227 }
1228
1229
1230 void postRun() {
1231
1232 if (integrator instanceof Respa) {
1233 potential.setEnergyTermState(Potential.STATE.BOTH);
1234 }
1235
1236
1237 if (!terminate) {
1238 logger.log(basicLogging, format(" Completed %8d time steps\n", nSteps));
1239 }
1240
1241
1242 done = true;
1243 terminate = false;
1244 }
1245
1246
1247
1248
1249
1250
1251 protected void appendSnapshot(String[] extraLines) {
1252
1253 for (MolecularAssembly assembly : molecularAssembly) {
1254 File archiveFile = assembly.getArchiveFile();
1255 ForceField forceField = assembly.getForceField();
1256 CompositeConfiguration properties = assembly.getProperties();
1257
1258
1259 if (archiveFile != null && !saveSnapshotAsPDB) {
1260 String aiName = FileUtils.relativePathTo(archiveFile).toString();
1261 if (esvSystem == null) {
1262 XYZFilter xyzFilter = new XYZFilter(archiveFile, assembly, forceField, properties);
1263 if (xyzFilter.writeFile(archiveFile, true, extraLines)) {
1264 logger.log(basicLogging, format(" Appended snapshot to: %s", aiName));
1265 } else {
1266 logger.warning( format(" Appending snapshot failed: %s", aiName));
1267 }
1268 } else {
1269 XPHFilter xphFilter = new XPHFilter(archiveFile, assembly, forceField, properties,
1270 esvSystem);
1271 if (xphFilter.writeFile(archiveFile, true, extraLines)) {
1272 logger.log(basicLogging, format(" Appended to XPH archive %s", aiName));
1273 } else {
1274 logger.warning(format(" Appending to XPH archive %s failed.", aiName));
1275 }
1276 }
1277 } else if (saveSnapshotAsPDB) {
1278 File file = assembly.getFile();
1279 String extName = FilenameUtils.getExtension(file.getName());
1280 File pdbFile;
1281 if (extName.toLowerCase().startsWith("pdb")) {
1282 pdbFile = file;
1283 } else {
1284 String filename = FilenameUtils.removeExtension(file.getAbsolutePath());
1285 pdbFile = new File(filename + ".pdb");
1286 }
1287 String aiName = FileUtils.relativePathTo(pdbFile).toString();
1288 PDBFilter pdbFilter = new PDBFilter(pdbFile, assembly, forceField, properties);
1289 if (pdbFilter.writeFile(pdbFile, true, extraLines)) {
1290 logger.log(basicLogging, format(" Appended to PDB file %s", aiName));
1291 } else {
1292 logger.warning(format(" Appending to PDB file to %s failed.", aiName));
1293 }
1294 }
1295 }
1296 }
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306 protected long logThermoForTime(long step, long time) {
1307 if (step % logFrequency == 0) {
1308 return logThermodynamics(time);
1309 } else {
1310 return time;
1311 }
1312 }
1313
1314
1315
1316
1317
1318
1319
1320 private long logThermodynamics(long time) {
1321 time = System.nanoTime() - time;
1322 logger.log(basicLogging,
1323 format(" %7.3e %12.4f %12.4f %12.4f %8.2f %8.3f", totalSimTime, state.getKineticEnergy(),
1324 state.getPotentialEnergy(), state.getTotalEnergy(), state.getTemperature(),
1325 time * NS2SEC));
1326 return System.nanoTime();
1327 }
1328
1329
1330 private void mainLoop() {
1331
1332
1333 long time = System.nanoTime();
1334 for (long step = 1; step <= nSteps; step++) {
1335 if (step > 1) {
1336 List<Constraint> constraints = potential.getConstraints();
1337
1338 long constraintFails = constraints.stream()
1339 .filter((Constraint c) -> !c.constraintSatisfied(state.x(), state.v(), 1E-7, 1E-7))
1340 .count();
1341 if (constraintFails > 0) {
1342 logger.info(format(" %d constraint failures in step %d", constraintFails, step));
1343 }
1344 }
1345
1346
1347 thermostat.halfStep(dt);
1348
1349
1350 integrator.preForce(potential);
1351 if (esvSystem != null) {
1352 esvIntegrator.preForce(potential);
1353
1354 esvSystem.preForce();
1355 }
1356
1357
1358 if (esvSystem != null && potential instanceof OpenMMEnergy) {
1359 state.setPotentialEnergy(((OpenMMEnergy) potential).energyAndGradientFFX(state.x(), state.gradient()));
1360 } else {
1361 state.setPotentialEnergy(potential.energyAndGradient(state.x(), state.gradient()));
1362 }
1363
1364
1365 if (integrator instanceof Respa r) {
1366 double potentialEnergy = state.getPotentialEnergy();
1367 state.setPotentialEnergy(potentialEnergy + r.getHalfStepEnergy());
1368 }
1369
1370
1371 integrator.postForce(state.gradient());
1372 if (esvSystem != null) {
1373 double[] dEdL = esvSystem.postForce();
1374 esvIntegrator.postForce(dEdL);
1375 }
1376
1377
1378 thermostat.computeKineticEnergy();
1379
1380
1381 thermostat.fullStep(dt);
1382
1383
1384 thermostat.computeKineticEnergy();
1385 if (esvSystem != null) {
1386
1387 esvThermostat.computeKineticEnergy();
1388 }
1389
1390
1391 int removeCOMMotionFrequency = 100;
1392 if (thermostat.getRemoveCenterOfMassMotion() && step % removeCOMMotionFrequency == 0) {
1393 thermostat.centerOfMassMotion(true, false);
1394 }
1395
1396
1397 if (esvSystem != null) {
1398 double kineticEnergy = thermostat.getKineticEnergy();
1399 double esvKineticEnergy = esvThermostat.getKineticEnergy();
1400 state.setKineticEnergy(kineticEnergy + esvKineticEnergy);
1401 }
1402
1403
1404 temperatureStats.addValue(state.getTemperature());
1405 potentialEnergyStats.addValue(state.getPotentialEnergy());
1406 kineticEnergyStats.addValue(state.getKineticEnergy());
1407 totalEnergyStats.addValue(state.getTotalEnergy());
1408
1409
1410 potential.setVelocity(state.v());
1411 potential.setAcceleration(state.a());
1412 potential.setPreviousAcceleration(state.aPrevious());
1413
1414
1415 totalSimTime += dt;
1416 time = logThermoForTime(step, time);
1417 if (step % printEsvFrequency == 0 && esvSystem != null) {
1418 logger.log(basicLogging, format(" %s", esvSystem.getLambdaList()));
1419 }
1420
1421 if (automaticWriteouts) {
1422 writeFilesForStep(step, true, true);
1423 }
1424
1425
1426 if (algorithmListener != null && step % logFrequency == 0) {
1427 for (MolecularAssembly assembly : molecularAssembly) {
1428 algorithmListener.algorithmUpdate(assembly);
1429 }
1430 }
1431
1432
1433 if (terminate) {
1434 logger.info(format("\n Terminating after %8d time steps\n", step));
1435 break;
1436 }
1437 }
1438 }
1439
1440
1441
1442
1443
1444
1445
1446
1447 private void checkFrequency(String describe, int frequency) {
1448 if (frequency > nSteps) {
1449 logger.fine(
1450 format(" Specified %s frequency of %d is greater than the number of steps %d", describe,
1451 frequency, nSteps));
1452 }
1453 }
1454
1455
1456
1457
1458
1459
1460 public void setCoordinates(double[] coords) {
1461 state.setCoordinates(coords);
1462 }
1463
1464
1465
1466
1467
1468
1469 public double[] getCoordinates() {
1470 return state.getCoordinatesCopy();
1471 }
1472
1473
1474
1475
1476 public void storeState() {
1477 storedState = state.getUnmodifiableState();
1478 }
1479
1480
1481
1482
1483
1484
1485 public void revertState() throws Exception {
1486 if (storedState == null) {
1487 throw new Exception();
1488 }
1489 revertState(storedState);
1490 }
1491
1492
1493
1494
1495
1496
1497 private void revertState(UnmodifiableState state) {
1498 this.state.revertState(state);
1499 potential.setVelocity(state.v());
1500 potential.setAcceleration(state.a());
1501 potential.setPreviousAcceleration(state.aPrevious());
1502
1503 int numberOfVariables = this.state.getNumberOfVariables();
1504 Atom[] atoms = molecularAssembly[0].getActiveAtomArray();
1505 if (atoms.length * 3 == numberOfVariables) {
1506 int index = 0;
1507 double[] x = state.x();
1508 double[] gradient = state.gradient();
1509 for (Atom atom : atoms) {
1510 atom.moveTo(x[index], x[index + 1], x[index + 2]);
1511 atom.setXYZGradient(gradient[index], gradient[index + 1], gradient[index + 2]);
1512 index += 3;
1513 }
1514 }
1515 }
1516
1517 }