1 #!/bin/env python
2 #
3 # File: OpenMMExecuteMDSimulationProtocol.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Acknowledgment: Paul Charifson
7 #
8 # Copyright (C) 2026 Manish Sud. All rights reserved.
9 #
10 # The functionality available in this script is implemented using OpenMM, an
11 # open source molecuar simulation package.
12 #
13 # This file is part of MayaChemTools.
14 #
15 # MayaChemTools is free software; you can redistribute it and/or modify it under
16 # the terms of the GNU Lesser General Public License as published by the Free
17 # Software Foundation; either version 3 of the License, or (at your option) any
18 # later version.
19 #
20 # MayaChemTools is distributed in the hope that it will be useful, but without
21 # any warranty; without even the implied warranty of merchantability of fitness
22 # for a particular purpose. See the GNU Lesser General Public License for more
23 # details.
24 #
25 # You should have received a copy of the GNU Lesser General Public License
26 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
27 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
28 # Boston, MA, 02111-1307, USA.
29 #
30
31
32 # Add local python path to the global path and import standard library modules...
33 import os
34 import sys
35 import time
36 import re
37
38 import pandas as pd
39 import matplotlib.pyplot as plt
40 import seaborn as sns
41
42 # OpenMM imports...
43 try:
44 import openmm as mm
45 import openmm.app
46 except ImportError as ErrMsg:
47 sys.stderr.write("\nFailed to import OpenMM related module/package: %s\n" % ErrMsg)
48 sys.stderr.write("Check/update your OpenMM environment and try again.\n\n")
49 sys.exit(1)
50
51 # MayaChemTools imports...
52 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
53 try:
54 from docopt import docopt
55 import MiscUtil
56 import OpenMMUtil
57 except ImportError as ErrMsg:
58 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
59 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
60 sys.exit(1)
61
62 ScriptName = os.path.basename(sys.argv[0])
63 Options = {}
64 OptionsInfo = {}
65
66
67 def main():
68 """Start execution of the script."""
69
70 MiscUtil.PrintInfo(
71 "\n%s (OpenMM v%s; MayaChemTools v%s; %s): Starting...\n"
72 % (ScriptName, mm.Platform.getOpenMMVersion(), MiscUtil.GetMayaChemToolsVersion(), time.asctime())
73 )
74
75 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
76
77 # Retrieve command line arguments and options...
78 RetrieveOptions()
79
80 # Process and validate command line arguments and options...
81 ProcessOptions()
82
83 # Perform actions required by the script...
84 ExecuteMDSimulationProtocol()
85
86 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
87 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
88
89
90 def ExecuteMDSimulationProtocol():
91 """Execute MD simulation protocol."""
92
93 # Prepare system for simulation...
94 System, Topology, Positions = PrepareSystem()
95
96 # Freeze and restraint atoms...
97 FreezeRestraintAtoms(System, Topology, Positions)
98
99 # Setup integrator...
100 Integrator = SetupIntegrator()
101
102 # Setup simulation...
103 Simulation = SetupSimulation(System, Integrator, Topology, Positions)
104
105 # Write setup files...
106 WriteSimulationSetupFiles(System, Integrator)
107
108 # Perform minimization...
109 PerformMinimization(Simulation)
110
111 # Set up intial velocities...
112 SetupInitialVelocities(Simulation)
113
114 # Execute MD protocol workflow...
115 ExecuteMDSimulationProtocolWorkflow(System, Simulation, Integrator)
116
117 # Save final state files...
118 WriteFinalStateFiles(Simulation)
119
120 # Reimage and realign trajectory for periodic systems...
121 ProcessTrajectory(System, Topology)
122
123 # Fix column name in data log file..
124 ProcessDataLogFile()
125
126 # Generate plots using data in log file...
127 GeneratePlots()
128
129
130 def PrepareSystem():
131 """Prepare system for simulation."""
132
133 System, Topology, Positions = OpenMMUtil.InitializeSystem(
134 OptionsInfo["Infile"],
135 OptionsInfo["ForcefieldParams"],
136 OptionsInfo["SystemParams"],
137 OptionsInfo["WaterBox"],
138 OptionsInfo["WaterBoxParams"],
139 OptionsInfo["SmallMolFile"],
140 OptionsInfo["SmallMolID"],
141 )
142
143 if OptionsInfo["MDProtocolParams"]["Phase4"] or OptionsInfo["MDProtocolParams"]["Phase5"]:
144 if not OpenMMUtil.DoesSystemUsesPeriodicBoundaryConditions(System):
145 MiscUtil.PrintInfo("")
146 MiscUtil.PrintWarning(
147 "A barostat is required for NPT equilibration and production simulations during phase 4 and 5. It appears that your system is a non-periodic system and OpenMM may fail during the addition of a barostat for phase 4 and 5. You must specify a periodic system or add water box to automatically set up a periodic system. "
148 )
149
150 MiscUtil.PrintInfo("\nChanging directory to %s..." % OptionsInfo["OutfileDir"])
151 os.chdir(OptionsInfo["OutfileDirPath"])
152
153 # Write out a PDB file for the system...
154 PDBFile = OptionsInfo["PDBOutfile"]
155 MiscUtil.PrintInfo("\nWriting PDB file %s..." % PDBFile)
156 OpenMMUtil.WritePDBFile(PDBFile, Topology, Positions, OptionsInfo["OutputParams"]["PDBOutKeepIDs"])
157
158 return (System, Topology, Positions)
159
160
161 def SetupIntegrator():
162 """Setup integrator."""
163
164 Integrator = OpenMMUtil.InitializeIntegrator(
165 OptionsInfo["IntegratorParams"], OptionsInfo["SystemParams"]["ConstraintErrorTolerance"]
166 )
167
168 return Integrator
169
170
171 def SetupSimulation(System, Integrator, Topology, Positions):
172 """Setup simulation."""
173
174 Simulation = OpenMMUtil.InitializeSimulation(System, Integrator, Topology, Positions, OptionsInfo["PlatformParams"])
175
176 return Simulation
177
178
179 def SetupInitialVelocities(Simulation):
180 """Setup initial velocities."""
181
182 # Set velocities to random values choosen from a Boltzman distribution at a given
183 # temperature...
184 MDProtocolParams = OptionsInfo["MDProtocolParams"]
185 MDProtocolParamsInfo = OpenMMUtil.SetupMDProtocolParameters(OptionsInfo["MDProtocolParams"])
186
187 TemperatureParamName = "Phase1InitialStart" if MDProtocolParams["Phase1"] else "Phase1InitialEnd"
188
189 MiscUtil.PrintInfo("\nSetting initial velocities to temperature (%s K)..." % MDProtocolParams[TemperatureParamName])
190 Simulation.context.setVelocitiesToTemperature(MDProtocolParamsInfo[TemperatureParamName])
191
192
193 def PerformMinimization(Simulation):
194 """Perform minimization."""
195
196 SimulationParams = OpenMMUtil.SetupSimulationParameters(OptionsInfo["SimulationParams"])
197
198 if not SimulationParams["Minimization"]:
199 MiscUtil.PrintInfo("\nSkipping energy minimization...")
200 return
201
202 OutputParams = OptionsInfo["OutputParams"]
203
204 # Setup a local minimization reporter...
205 MinimizeReporter = None
206 if OutputParams["MinimizationDataStdout"] or OutputParams["MinimizationDataLog"]:
207 MinimizeReporter = LocalMinimizationReporter()
208
209 if MinimizeReporter is not None:
210 MiscUtil.PrintInfo("\nAdding minimization reporters...")
211 if OutputParams["MinimizationDataLog"]:
212 MiscUtil.PrintInfo(
213 "Adding data log minimization reporter (Steps: %s; File: %s)..."
214 % (OutputParams["MinimizationDataSteps"], OutputParams["MinimizationDataLogFile"])
215 )
216 if OutputParams["MinimizationDataStdout"]:
217 MiscUtil.PrintInfo(
218 "Adding data stdout minimization reporter (Steps: %s)..." % (OutputParams["MinimizationDataSteps"])
219 )
220 else:
221 MiscUtil.PrintInfo("\nSkipping addition of minimization reporters...")
222
223 MaxSteps = SimulationParams["MinimizationMaxSteps"]
224
225 MaxStepsMsg = "MaxSteps: %s" % ("UntilConverged" if MaxSteps == 0 else MaxSteps)
226 ToleranceMsg = "Tolerance: %.2f kcal/mol/A (%.2f kjoules/mol/nm)" % (
227 SimulationParams["MinimizationToleranceInKcal"],
228 SimulationParams["MinimizationToleranceInJoules"],
229 )
230
231 MiscUtil.PrintInfo("\nPerforming energy minimization (%s; %s)..." % (MaxStepsMsg, ToleranceMsg))
232
233 if OutputParams["MinimizationDataStdout"]:
234 HeaderLine = SetupMinimizationDataOutHeaderLine()
235 print("\n%s" % HeaderLine)
236
237 Simulation.minimizeEnergy(
238 tolerance=SimulationParams["MinimizationTolerance"], maxIterations=MaxSteps, reporter=MinimizeReporter
239 )
240
241 if OutputParams["MinimizationDataLog"]:
242 WriteMinimizationDataLogFile(MinimizeReporter.DataOutValues)
243
244 if OutputParams["PDBOutMinimized"]:
245 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["MinimizedPDBOutfile"])
246 OpenMMUtil.WriteSimulationStatePDBFile(
247 Simulation, OptionsInfo["MinimizedPDBOutfile"], OutputParams["PDBOutKeepIDs"]
248 )
249
250
251 def ExecuteMDSimulationProtocolWorkflow(System, Simulation, Integrator):
252 """Execute MD simulation prorocol workflow."""
253
254 MiscUtil.PrintInfo("\nExecuting MD simulation protocol...")
255
256 if OptionsInfo["OutputReportersModeAllPhases"]:
257 SetupReporters(Simulation)
258
259 TotalSimulationSteps = 0
260
261 Phase1SimulationSteps = PerformPhase1InitialHeating(Simulation, Integrator)
262 TotalSimulationSteps += Phase1SimulationSteps
263
264 Phase2SimulationSteps = PerformPhase2HeatingAndCooling(Simulation, Integrator)
265 TotalSimulationSteps += Phase2SimulationSteps
266
267 Phase3SimulationSteps = PerformPhase3Equilibration(Simulation, Integrator)
268 TotalSimulationSteps += Phase3SimulationSteps
269
270 if OptionsInfo["MDProtocolParams"]["Phase4"] or OptionsInfo["MDProtocolParams"]["Phase5"]:
271 Barostat = OpenMMUtil.InitializeBarostat(OptionsInfo["IntegratorParams"])
272 MiscUtil.PrintInfo("Adding barostat for NPT simulation...")
273 try:
274 System.addForce(Barostat)
275 Simulation.context.reinitialize(preserveState=True)
276 except Exception as ErrMsg:
277 MiscUtil.PrintInfo("")
278 MiscUtil.PrintError("Failed to add barostat:\n%s\n" % (ErrMsg))
279
280 Phase4SimulationSteps = PerformPhase4Equilibration(Simulation, Integrator)
281 TotalSimulationSteps += Phase4SimulationSteps
282
283 if OptionsInfo["MDProtocolParams"]["Phase5"]:
284 if not OptionsInfo["OutputReportersModeAllPhases"]:
285 SetupReporters(Simulation)
286
287 Phase5SimulationSteps = PerformPhase5ProductionRun(Simulation, Integrator)
288 TotalSimulationSteps += Phase5SimulationSteps
289
290 MiscUtil.PrintInfo(
291 "\nFinishing executing MD protocol (TotalSteps: %s; TotalTime: %s)..."
292 % (TotalSimulationSteps, GetTotalSimulationTime(TotalSimulationSteps))
293 )
294
295
296 def PerformPhase1InitialHeating(Simulation, Integrator):
297 """Perform phase 1 NVT initial heating."""
298
299 if not OptionsInfo["MDProtocolParams"]["Phase1"]:
300 MiscUtil.PrintInfo("\nSkipping phase 1 initial heating...")
301 return 0
302
303 MDProtocolParams = OptionsInfo["MDProtocolParams"]
304 OutputParams = OptionsInfo["OutputParams"]
305
306 # Perform intial heating along with equilibration...
307 InitialStart = MDProtocolParams["Phase1InitialStart"]
308 InitialEnd = MDProtocolParams["Phase1InitialEnd"]
309 InitialChange = MDProtocolParams["Phase1InitialChange"]
310 InitialSteps = MDProtocolParams["Phase1InitialSteps"]
311
312 Barostat = None
313 TotalSimulationSteps = 0
314
315 MiscUtil.PrintInfo(
316 "\nPerforming phase 1 initial heating (Ensemble: NVT; Start: %.1f K; End: %.1f K; Change: %.1f K)..."
317 % (InitialStart, InitialEnd, InitialChange)
318 )
319 TotalInitialSimulationSteps = OpenMMUtil.PerformAnnealing(
320 Simulation, Integrator, Barostat, InitialStart, InitialEnd, InitialChange, InitialSteps
321 )
322 MiscUtil.PrintInfo(
323 "Finished initial heating (TotalSteps: %s; TotalTime: %s)..."
324 % (TotalInitialSimulationSteps, GetTotalSimulationTime(TotalInitialSimulationSteps))
325 )
326 TotalSimulationSteps += TotalInitialSimulationSteps
327
328 # Perform equilibration after intial heating...
329 InitialEquilibrationSteps = MDProtocolParams["Phase1InitialEquilibrationSteps"]
330 MiscUtil.PrintInfo(
331 "\nPerforming phase 1 equilibration after initial heating (Ensemble: NVT; Steps: %s; Time: %s)..."
332 % (InitialEquilibrationSteps, GetTotalSimulationTime(InitialEquilibrationSteps))
333 )
334 Simulation.step(InitialEquilibrationSteps)
335 TotalSimulationSteps += InitialEquilibrationSteps
336
337 if OutputParams["PDBOutPhase1HeatedNVT"]:
338 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["Phase1HeatedNVTPDBOutfile"])
339 OpenMMUtil.WriteSimulationStatePDBFile(
340 Simulation, OptionsInfo["Phase1HeatedNVTPDBOutfile"], OutputParams["PDBOutKeepIDs"]
341 )
342
343 return TotalSimulationSteps
344
345
346 def PerformPhase2HeatingAndCooling(Simulation, Integrator):
347 """Perform phase 2 NVT heating and cooling."""
348
349 if not OptionsInfo["MDProtocolParams"]["Phase2"]:
350 MiscUtil.PrintInfo("\nSkipping phase 2 heating and cooling...")
351 return 0
352
353 MDProtocolParams = OptionsInfo["MDProtocolParams"]
354 OutputParams = OptionsInfo["OutputParams"]
355
356 # Peform heating and coolling annealing cycles along with equilibration...
357 Cycles = MDProtocolParams["Phase2Cycles"]
358 CycleStart = MDProtocolParams["Phase2CycleStart"]
359 CycleEnd = MDProtocolParams["Phase2CycleEnd"]
360 CycleChange = MDProtocolParams["Phase2CycleChange"]
361 CycleSteps = MDProtocolParams["Phase2CycleSteps"]
362 CycleEquilibrationSteps = MDProtocolParams["Phase2CycleEquilibrationSteps"]
363
364 Barostat = None
365 TotalSimulationSteps = 0
366
367 MiscUtil.PrintInfo("\nPerforming phase 2 heating and cooling cycles (NumCycles: %s)..." % (Cycles))
368 for Cycle in range(Cycles):
369 MiscUtil.PrintInfo("\nPerforming heating and cooling cycle %s..." % (Cycle + 1))
370
371 MiscUtil.PrintInfo(
372 "\nPerforming heating (Start: %.1f K; End: %.1f K; Change: %.1f K)..." % (CycleStart, CycleEnd, CycleChange)
373 )
374 TotalCycleSimulationSteps = OpenMMUtil.PerformAnnealing(
375 Simulation, Integrator, Barostat, CycleStart, CycleEnd, CycleChange, CycleSteps
376 )
377 MiscUtil.PrintInfo(
378 "Finished heating cycle (TotalSteps: %s; TotalTime: %s)..."
379 % (TotalCycleSimulationSteps, GetTotalSimulationTime(TotalCycleSimulationSteps))
380 )
381 TotalSimulationSteps += TotalCycleSimulationSteps
382
383 MiscUtil.PrintInfo(
384 "\nPerforming equilibration (Steps: %s; Time: %s)..."
385 % (CycleEquilibrationSteps, GetTotalSimulationTime(CycleEquilibrationSteps))
386 )
387 Simulation.step(CycleEquilibrationSteps)
388 TotalSimulationSteps += CycleEquilibrationSteps
389
390 MiscUtil.PrintInfo(
391 "\nPerforming cooling (Start: %.1f K; End: %.1f K; Change: %.1f K)..." % (CycleEnd, CycleStart, CycleChange)
392 )
393 TotalCycleSimulationSteps = OpenMMUtil.PerformAnnealing(
394 Simulation, Integrator, Barostat, CycleEnd, CycleStart, CycleChange, CycleSteps
395 )
396 MiscUtil.PrintInfo(
397 "Finished cooling cycle (TotalSteps: %s; TotalTime: %s)..."
398 % (TotalCycleSimulationSteps, GetTotalSimulationTime(TotalCycleSimulationSteps))
399 )
400 TotalSimulationSteps += TotalCycleSimulationSteps
401
402 MiscUtil.PrintInfo(
403 "\nPerforming equilibration (Steps: %s; Time: %s)..."
404 % (CycleEquilibrationSteps, GetTotalSimulationTime(CycleEquilibrationSteps))
405 )
406 Simulation.step(CycleEquilibrationSteps)
407 TotalSimulationSteps += CycleEquilibrationSteps
408
409 MiscUtil.PrintInfo("\nFinished heating and cooling cycle %s..." % (Cycle + 1))
410
411 if OutputParams["PDBOutPhase2AnnealedNVT"]:
412 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["Phase2AnnealedNVTPDBOutfile"])
413 OpenMMUtil.WriteSimulationStatePDBFile(
414 Simulation, OptionsInfo["Phase2AnnealedNVTPDBOutfile"], OutputParams["PDBOutKeepIDs"]
415 )
416
417 return TotalSimulationSteps
418
419
420 def PerformPhase3Equilibration(Simulation, Integrator):
421 """Perform phase 3 NVT equilibration."""
422
423 if not OptionsInfo["MDProtocolParams"]["Phase3"]:
424 MiscUtil.PrintInfo("\nSkipping phase 3 equilibration...")
425 return 0
426
427 MDProtocolParams = OptionsInfo["MDProtocolParams"]
428 OutputParams = OptionsInfo["OutputParams"]
429
430 Phase3Steps = MDProtocolParams["Phase3Steps"]
431 MiscUtil.PrintInfo(
432 "\nPerforming phase 3 equilibration (Ensemble: NVT; Steps: %s; Time: %s)..."
433 % (Phase3Steps, GetTotalSimulationTime(Phase3Steps))
434 )
435 Simulation.step(Phase3Steps)
436
437 if OutputParams["PDBOutPhase3EquilibratedNVT"]:
438 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["Phase3EquilibratedNVTPDBOutfile"])
439 OpenMMUtil.WriteSimulationStatePDBFile(
440 Simulation, OptionsInfo["Phase3EquilibratedNVTPDBOutfile"], OutputParams["PDBOutKeepIDs"]
441 )
442
443 return Phase3Steps
444
445
446 def PerformPhase4Equilibration(Simulation, Integrator):
447 """Perform phase 4 NPT equilibration."""
448
449 if not OptionsInfo["MDProtocolParams"]["Phase4"]:
450 MiscUtil.PrintInfo("\nSkipping phase 4 equilibration...")
451 return 0
452
453 MDProtocolParams = OptionsInfo["MDProtocolParams"]
454 OutputParams = OptionsInfo["OutputParams"]
455
456 Phase4Steps = MDProtocolParams["Phase4Steps"]
457 MiscUtil.PrintInfo(
458 "\nPerforming phase 4 equilibration (Ensemble: NPT; Steps: %s; Time: %s)..."
459 % (Phase4Steps, GetTotalSimulationTime(Phase4Steps))
460 )
461 Simulation.step(Phase4Steps)
462
463 if OutputParams["PDBOutPhase4EquilibratedNPT"]:
464 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["Phase4EquilibratedNPTPDBOutfile"])
465 OpenMMUtil.WriteSimulationStatePDBFile(
466 Simulation, OptionsInfo["Phase4EquilibratedNPTPDBOutfile"], OutputParams["PDBOutKeepIDs"]
467 )
468
469 return Phase4Steps
470
471
472 def PerformPhase5ProductionRun(Simulation, Integrator):
473 """Perform phase 5 NPT production run.."""
474
475 if not OptionsInfo["MDProtocolParams"]["Phase5"]:
476 MiscUtil.PrintInfo("\nSkipping phase 5 equilibration...")
477 return 0
478
479 MDProtocolParamsInfo = OpenMMUtil.SetupMDProtocolParameters(OptionsInfo["MDProtocolParams"])
480 OutputParams = OptionsInfo["OutputParams"]
481
482 Phase5Steps = MDProtocolParamsInfo["Phase5Steps"]
483 Phase5StepSize = MDProtocolParamsInfo["Phase5StepSize"]
484 if Phase5StepSize is not None:
485 # Setup step size for phase 5 simulation..
486 MiscUtil.PrintInfo("\nModifying step size for phase 5 (StepSize: %s)..." % Phase5StepSize)
487 Integrator.setStepSize(Phase5StepSize)
488
489 MiscUtil.PrintInfo(
490 "\nPerforming phase 5 production run (Ensemble: NPT; Steps: %s; Time: %s)..."
491 % (Phase5Steps, GetTotalSimulationTime(Phase5Steps, Phase5StepSize))
492 )
493 Simulation.step(Phase5Steps)
494
495 if OutputParams["PDBOutPhase5ProductionNPT"]:
496 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["Phase5ProductionNPTPDBOutfile"])
497 OpenMMUtil.WriteSimulationStatePDBFile(
498 Simulation, OptionsInfo["Phase5ProductionNPTPDBOutfile"], OutputParams["PDBOutKeepIDs"]
499 )
500
501 return Phase5Steps
502
503
504 def GetTotalSimulationTime(SimulationSteps, StepSize=None):
505 """Get total simulation time."""
506
507 if StepSize is None:
508 IntegratorParamsInfo = OpenMMUtil.SetupIntegratorParameters(OptionsInfo["IntegratorParams"])
509 StepSize = IntegratorParamsInfo["StepSize"]
510
511 TotalTime = OpenMMUtil.GetFormattedTotalSimulationTime(StepSize, SimulationSteps)
512
513 return TotalTime
514
515
516 def SetupReporters(Simulation):
517 """Setup reporters."""
518
519 DataAppend = False
520 (TrajReporter, DataLogReporter, DataStdoutReporter, CheckpointReporter) = OpenMMUtil.InitializeReporters(
521 OptionsInfo["OutputParams"], OptionsInfo["SimulationParams"]["Steps"], DataAppend
522 )
523
524 if TrajReporter is None and DataLogReporter is None and DataStdoutReporter is None and CheckpointReporter is None:
525 MiscUtil.PrintInfo("\nSkip adding reporters...")
526 return
527
528 MiscUtil.PrintInfo("\nAdding reporters...")
529
530 OutputParams = OptionsInfo["OutputParams"]
531 AppendMsg = ""
532 if TrajReporter is not None:
533 MiscUtil.PrintInfo(
534 "Adding trajectory reporter (Steps: %s; File: %s%s)..."
535 % (OutputParams["TrajSteps"], OutputParams["TrajFile"], AppendMsg)
536 )
537 Simulation.reporters.append(TrajReporter)
538
539 if CheckpointReporter is not None:
540 MiscUtil.PrintInfo(
541 "Adding checkpoint reporter (Steps: %s; File: %s)..."
542 % (OutputParams["CheckpointSteps"], OutputParams["CheckpointFile"])
543 )
544 Simulation.reporters.append(CheckpointReporter)
545
546 if DataLogReporter is not None:
547 MiscUtil.PrintInfo(
548 "Adding data log reporter (Steps: %s; File: %s%s)..."
549 % (OutputParams["DataLogSteps"], OutputParams["DataLogFile"], AppendMsg)
550 )
551 Simulation.reporters.append(DataLogReporter)
552
553 if DataStdoutReporter is not None:
554 MiscUtil.PrintInfo("Adding data stdout reporter (Steps: %s)..." % (OutputParams["DataStdoutSteps"]))
555 Simulation.reporters.append(DataStdoutReporter)
556
557
558 class LocalMinimizationReporter(mm.MinimizationReporter):
559 """Setup a local minimization reporter."""
560
561 (DataSteps, DataOutTypeList, DataOutDelimiter, StdoutStatus) = [None] * 4
562
563 DataOutValues = []
564 First = True
565
566 def report(self, Iteration, PositonsList, GradientsList, DataStatisticsMap):
567 """Report and track minimization."""
568
569 if self.First:
570 # Initialize...
571 self.DataSteps = OptionsInfo["OutputParams"]["MinimizationDataSteps"]
572 self.DataOutTypeList = OptionsInfo["OutputParams"]["MinimizationDataOutTypeOpenMMNameList"]
573 self.DataOutDelimiter = OptionsInfo["OutputParams"]["DataOutDelimiter"]
574 self.StdoutStatus = True if OptionsInfo["OutputParams"]["MinimizationDataStdout"] else False
575
576 self.First = False
577
578 if Iteration % self.DataSteps == 0:
579 # Setup data values...
580 DataValues = []
581 DataValues.append("%s" % Iteration)
582 for DataType in self.DataOutTypeList:
583 DataValue = "%.4f" % DataStatisticsMap[DataType]
584 DataValues.append(DataValue)
585
586 # Track data...
587 self.DataOutValues.append(DataValues)
588
589 # Print data values...
590 if self.StdoutStatus:
591 print("%s" % self.DataOutDelimiter.join(DataValues))
592
593 # This method must return a bool. You may return true for early termination.
594 return False
595
596
597 def WriteMinimizationDataLogFile(DataOutValues):
598 """Write minimization data log file."""
599
600 OutputParams = OptionsInfo["OutputParams"]
601
602 Outfile = OutputParams["MinimizationDataLogFile"]
603 OutDelimiter = OutputParams["DataOutDelimiter"]
604
605 MiscUtil.PrintInfo("\nWriting minimization log file %s..." % Outfile)
606 OutFH = open(Outfile, "w")
607
608 HeaderLine = SetupMinimizationDataOutHeaderLine()
609 OutFH.write("%s\n" % HeaderLine)
610
611 for LineWords in DataOutValues:
612 Line = OutDelimiter.join(LineWords)
613 OutFH.write("%s\n" % Line)
614
615 OutFH.close()
616
617
618 def SetupMinimizationDataOutHeaderLine():
619 """Setup minimization data output header line."""
620
621 LineWords = ["Iteration"]
622 for Label in OptionsInfo["OutputParams"]["MinimizationDataOutTypeList"]:
623 if re.match("^(SystemEnergy|RestraintEnergy)$", Label, re.I):
624 LineWords.append("%s(kjoules/mol)" % Label)
625 elif re.match("^RestraintStrength$", Label, re.I):
626 LineWords.append("%s(kjoules/mol/nm^2)" % Label)
627 else:
628 LineWords.append(Label)
629
630 Line = OptionsInfo["OutputParams"]["DataOutDelimiter"].join(LineWords)
631
632 return Line
633
634
635 def FreezeRestraintAtoms(System, Topology, Positions):
636 """Handle freezing and restraining of atoms."""
637
638 FreezeAtomList, RestraintAtomList = OpenMMUtil.ValidateAndFreezeRestraintAtoms(
639 OptionsInfo["FreezeAtoms"],
640 OptionsInfo["FreezeAtomsParams"],
641 OptionsInfo["RestraintAtoms"],
642 OptionsInfo["RestraintAtomsParams"],
643 OptionsInfo["RestraintSpringConstant"],
644 OptionsInfo["SystemParams"],
645 System,
646 Topology,
647 Positions,
648 )
649
650 # Check and adjust step size...
651 if FreezeAtomList is not None or RestraintAtomList is not None:
652 if re.match("^auto$", OptionsInfo["IntegratorParams"]["StepSizeSpecified"], re.I):
653 # Automatically set stepSize to 2.0 fs..
654 OptionsInfo["IntegratorParams"]["StepSize"] = 2.0
655 MiscUtil.PrintInfo("")
656 MiscUtil.PrintWarning(
657 'The time step has been automatically set to %s fs during freezing or restraining of atoms. You may specify an explicit value for parameter name, stepSize, using "--integratorParams" option.'
658 % (OptionsInfo["IntegratorParams"]["StepSize"])
659 )
660 elif OptionsInfo["IntegratorParams"]["StepSize"] > 2:
661 MiscUtil.PrintInfo("")
662 MiscUtil.PrintWarning(
663 'A word to the wise: The parameter value specified, %s, for parameter name, stepSize, using "--integratorParams" option may be too large. You may want to consider using a smaller time step. Othwerwise, your simulation may blow up.'
664 % (OptionsInfo["IntegratorParams"]["StepSize"])
665 )
666 MiscUtil.PrintInfo("")
667
668
669 def WriteSimulationSetupFiles(System, Integrator):
670 """Write simulation setup files for system and integrator."""
671
672 OutputParams = OptionsInfo["OutputParams"]
673
674 if OutputParams["XmlSystemOut"] or OutputParams["XmlIntegratorOut"]:
675 MiscUtil.PrintInfo("")
676
677 if OutputParams["XmlSystemOut"]:
678 Outfile = OutputParams["XmlSystemFile"]
679 MiscUtil.PrintInfo("Writing system setup XML file %s..." % Outfile)
680 with open(Outfile, mode="w") as OutFH:
681 OutFH.write(mm.XmlSerializer.serialize(System))
682
683 if OutputParams["XmlIntegratorOut"]:
684 Outfile = OutputParams["XmlIntegratorFile"]
685 MiscUtil.PrintInfo("Writing integrator setup XML file %s..." % Outfile)
686 with open(Outfile, mode="w") as OutFH:
687 OutFH.write(mm.XmlSerializer.serialize(Integrator))
688
689
690 def WriteFinalStateFiles(Simulation):
691 """Write final state files."""
692
693 OutputParams = OptionsInfo["OutputParams"]
694
695 if OutputParams["SaveFinalStateCheckpoint"] or OutputParams["SaveFinalStateXML"] or OutputParams["PDBOutFinal"]:
696 MiscUtil.PrintInfo("")
697
698 if OutputParams["SaveFinalStateCheckpoint"]:
699 Outfile = OutputParams["SaveFinalStateCheckpointFile"]
700 MiscUtil.PrintInfo("Writing final state checkpoint file %s..." % Outfile)
701 Simulation.saveCheckpoint(Outfile)
702
703 if OutputParams["SaveFinalStateXML"]:
704 Outfile = OutputParams["SaveFinalStateXMLFile"]
705 MiscUtil.PrintInfo("Writing final state XML file %s..." % Outfile)
706 Simulation.saveState(Outfile)
707
708 if OutputParams["PDBOutFinal"]:
709 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["FinalPDBOutfile"])
710 OpenMMUtil.WriteSimulationStatePDBFile(
711 Simulation, OptionsInfo["FinalPDBOutfile"], OutputParams["PDBOutKeepIDs"]
712 )
713
714
715 def ProcessTrajectory(System, Topology):
716 """Reimage and realign trajectory for periodic systems."""
717
718 TrajTopologyFile = OptionsInfo["PDBOutfile"]
719
720 OpenMMUtil.GenerateReimagedRealignedTrajectoryFiles(
721 System,
722 Topology,
723 TrajTopologyFile,
724 OptionsInfo["ReimagedPDBOutfile"],
725 OptionsInfo["ReimagedTrajOutfile"],
726 OptionsInfo["OutputParams"],
727 RealignFrames=True,
728 )
729
730
731 def ProcessDataLogFile():
732 """Process data log file."""
733
734 OutputParams = OptionsInfo["OutputParams"]
735 if not OutputParams["DataLog"] or not os.path.exists(OutputParams["DataLogFile"]):
736 return
737
738 DataLogFile = OutputParams["DataLogFile"]
739 MiscUtil.PrintInfo("\nProcessing data log file %s..." % DataLogFile)
740
741 OpenMMUtil.FixColumNamesLineInDataLogFile(DataLogFile)
742
743
744 def GeneratePlots():
745 """Generate plots using data in log file."""
746
747 OutputParams = OptionsInfo["OutputParams"]
748 if (
749 not OutputParams["DataLog"]
750 or not OutputParams["DataOutTypePlot"]
751 or not os.path.exists(OutputParams["DataLogFile"])
752 ):
753 MiscUtil.PrintInfo("\nSkipping generation of plots...")
754 return
755
756 MiscUtil.PrintInfo("\nGenerating plots...")
757 InitializePlotParameters()
758
759 DataLogFile = OutputParams["DataLogFile"]
760
761 MiscUtil.PrintInfo("Processing file %s..." % DataLogFile)
762 DataLogDF = pd.read_csv(DataLogFile, sep=",")
763 DataLogColNames = DataLogDF.columns.tolist()
764
765 # Collect data types to plot...
766 DataOutTypePlotList = []
767 DataOutTypePlotList.append(OutputParams["DataOutTypePlotX"])
768 DataOutTypePlotList.extend(OutputParams["DataOutTypePlotYList"])
769 DataOutTypePlotColNames = OpenMMUtil.MapDataOutTypePlotToDataLogColumnNames(DataOutTypePlotList, DataLogColNames)
770
771 DataOutTypePlotFiles = OptionsInfo["DataOutTypePlotFiles"]
772
773 for PlotY in OutputParams["DataOutTypePlotYList"]:
774 PlotOutFile = DataOutTypePlotFiles[PlotY]
775
776 PlotX = OutputParams["DataOutTypePlotX"]
777 PlotXColName = DataOutTypePlotColNames[PlotX]
778 PlotYColName = DataOutTypePlotColNames[PlotY]
779
780 if PlotXColName is None or PlotYColName is None:
781 MiscUtil.PrintInfo(
782 "Skipping generation of plot file %s (Missing %s or %s data column in data log file)..."
783 % (PlotOutFile, PlotX, PlotY)
784 )
785 continue
786
787 PlotXLabel = PlotXColName
788 PlotYLabel = PlotYColName
789
790 if OptionsInfo["OutputReportersModeAllPhases"]:
791 PlotTitle = "MD Simulation Protocol"
792 else:
793 PlotTitle = "MD Production Simulation (NPT)"
794
795 GeneratePlotOutFile(PlotOutFile, DataLogDF, PlotXColName, PlotYColName, PlotXLabel, PlotYLabel, PlotTitle)
796
797
798 def GeneratePlotOutFile(PlotOutFile, DataLogDF, PlotXColName, PlotYColName, PlotXLabel, PlotYLabel, PlotTitle):
799 """Generate plot out file."""
800
801 OutPlotParams = OptionsInfo["OutPlotParams"]
802
803 MiscUtil.PrintInfo("Generating plot file %s..." % PlotOutFile)
804
805 # Create a new figure...
806 plt.figure()
807
808 # Draw plot...
809 PlotType = OutPlotParams["Type"]
810 if re.match("^line$", PlotType, re.I):
811 Axis = sns.lineplot(DataLogDF, x=PlotXColName, y=PlotYColName, legend=False)
812 elif re.match("^linepoint$", PlotType, re.I):
813 Axis = sns.lineplot(DataLogDF, x=PlotXColName, y=PlotYColName, marker="o", legend=False)
814 elif re.match("^scatter$", PlotType, re.I):
815 Axis = sns.scatterplot(DataLogDF, x=PlotXColName, y=PlotYColName, legend=False)
816 else:
817 MiscUtil.PrintError(
818 'The value, %s, specified for "type" using option "--outPlotParams" is not supported. Valid plot types: linepoint, scatter or line'
819 % (PlotType)
820 )
821
822 # Set labels and title...
823 Axis.set(xlabel=PlotXLabel, ylabel=PlotYLabel, title=PlotTitle)
824
825 # Save figure...
826 plt.savefig(PlotOutFile)
827
828 # Close the plot...
829 plt.close()
830
831
832 def InitializePlotParameters():
833 """Initialize plot parameters."""
834
835 if OptionsInfo["OutPlotInitialized"]:
836 return
837
838 # Initialize seaborn and matplotlib paramaters...
839 OptionsInfo["OutPlotInitialized"] = True
840
841 OutPlotParams = OptionsInfo["OutPlotParams"]
842 RCParams = {
843 "figure.figsize": (OutPlotParams["Width"], OutPlotParams["Height"]),
844 "axes.titleweight": OutPlotParams["TitleWeight"],
845 "axes.labelweight": OutPlotParams["LabelWeight"],
846 }
847 sns.set(
848 context=OutPlotParams["Context"],
849 style=OutPlotParams["Style"],
850 palette=OutPlotParams["Palette"],
851 font=OutPlotParams["Font"],
852 font_scale=OutPlotParams["FontScale"],
853 rc=RCParams,
854 )
855
856
857 def ProcessOutfilePrefixOption():
858 """Process outfile prefix option."""
859
860 OutfilePrefix = Options["--outfilePrefix"]
861
862 if not re.match("^auto$", OutfilePrefix, re.I):
863 OptionsInfo["OutfilePrefix"] = OutfilePrefix
864 return
865
866 if OptionsInfo["SmallMolFileMode"]:
867 OutfilePrefix = "%s_%s_Complex" % (OptionsInfo["InfileRoot"], OptionsInfo["SmallMolFileRoot"])
868 else:
869 OutfilePrefix = "%s" % (OptionsInfo["InfileRoot"])
870
871 if re.match("^yes$", Options["--waterBox"], re.I):
872 OutfilePrefix = "%s_Solvated" % (OutfilePrefix)
873
874 OptionsInfo["OutfilePrefix"] = OutfilePrefix
875
876
877 def ProcessOutfileDirOption():
878 """Process outfile directory Option."""
879
880 # Setup output directory...
881 OutfileDir = Options["--outfileDir"]
882 OutfileDirPath = os.path.abspath(OutfileDir)
883 if not os.path.exists(OutfileDir):
884 MiscUtil.PrintInfo("\nCreating output directory %s..." % (OutfileDir))
885 os.mkdir(OutfileDirPath)
886 OptionsInfo["OutfileDir"] = OutfileDir
887 OptionsInfo["OutfileDirPath"] = OutfileDirPath
888
889
890 def ProcessOutfileNames():
891 """Process outfile names."""
892
893 OutputParams = OptionsInfo["OutputParams"]
894
895 PDBOutfile = "%s.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
896 ReimagedPDBOutfile = "%s_Reimaged.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
897 ReimagedTrajOutfile = "%s_Reimaged.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["TrajFileExt"])
898
899 MinimizedPDBOutfile = "%s_Minimized.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
900 FinalPDBOutfile = "%s_Final.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
901
902 Phase1HeatedNVTPDBOutfile = "%s_NVT_Phase1_Heated_Equilibrated.%s" % (
903 OptionsInfo["OutfilePrefix"],
904 OutputParams["PDBOutfileExt"],
905 )
906 Phase2AnnealedNVTPDBOutfile = "%s_NVT_Phase2_Annealed_Equilibrated.%s" % (
907 OptionsInfo["OutfilePrefix"],
908 OutputParams["PDBOutfileExt"],
909 )
910 Phase3EquilibratedNVTPDBOutfile = "%s_NVT_Phase3_Equilibrated.%s" % (
911 OptionsInfo["OutfilePrefix"],
912 OutputParams["PDBOutfileExt"],
913 )
914 Phase4EquilibratedNPTPDBOutfile = "%s_NPT_Phase4_Equilibrated.%s" % (
915 OptionsInfo["OutfilePrefix"],
916 OutputParams["PDBOutfileExt"],
917 )
918 Phase5ProductionNPTPDBOutfile = "%s_NPT_Phase5_Production.%s" % (
919 OptionsInfo["OutfilePrefix"],
920 OutputParams["PDBOutfileExt"],
921 )
922
923 OptionsInfo["PDBOutfile"] = PDBOutfile
924 OptionsInfo["ReimagedPDBOutfile"] = ReimagedPDBOutfile
925 OptionsInfo["ReimagedTrajOutfile"] = ReimagedTrajOutfile
926
927 OptionsInfo["MinimizedPDBOutfile"] = MinimizedPDBOutfile
928 OptionsInfo["FinalPDBOutfile"] = FinalPDBOutfile
929
930 OptionsInfo["Phase1HeatedNVTPDBOutfile"] = Phase1HeatedNVTPDBOutfile
931 OptionsInfo["Phase2AnnealedNVTPDBOutfile"] = Phase2AnnealedNVTPDBOutfile
932 OptionsInfo["Phase3EquilibratedNVTPDBOutfile"] = Phase3EquilibratedNVTPDBOutfile
933 OptionsInfo["Phase4EquilibratedNPTPDBOutfile"] = Phase4EquilibratedNPTPDBOutfile
934 OptionsInfo["Phase5ProductionNPTPDBOutfile"] = Phase5ProductionNPTPDBOutfile
935
936 OutputParams = OptionsInfo["OutputParams"]
937 OutPlotParams = OptionsInfo["OutPlotParams"]
938
939 DataOutTypePlotYList = OutputParams["DataOutTypePlotYList"]
940 DataOutTypePlotFiles = {}
941 for PlotDataType in DataOutTypePlotYList:
942 Outfile = "%s_%sPlot.%s" % (OptionsInfo["OutfilePrefix"], PlotDataType, OutPlotParams["OutExt"])
943 if not Options["--overwrite"]:
944 if os.path.exists(Outfile):
945 MiscUtil.PrintError(
946 'The file name, %s, generated using option "--outfilePrefix" already exist. Use option "--ov" or "--overwrite" and try again. '
947 % (Outfile)
948 )
949 DataOutTypePlotFiles[PlotDataType] = Outfile
950 OptionsInfo["DataOutTypePlotFiles"] = DataOutTypePlotFiles
951
952
953 def ProcessWaterBoxParameters():
954 """Process water box parameters."""
955
956 OptionsInfo["WaterBox"] = True if re.match("^yes$", Options["--waterBox"], re.I) else False
957 OptionsInfo["WaterBoxParams"] = OpenMMUtil.ProcessOptionOpenMMWaterBoxParameters(
958 "--waterBoxParams", Options["--waterBoxParams"]
959 )
960
961 if OptionsInfo["WaterBox"]:
962 if OptionsInfo["ForcefieldParams"]["ImplicitWater"]:
963 MiscUtil.PrintInfo("")
964 MiscUtil.PrintWarning(
965 'The value, %s, specified using option "--waterBox" may not be valid for the combination of biopolymer and water forcefields, %s and %s, specified using "--forcefieldParams". You may consider using a valid combination of biopolymer and water forcefields for explicit water during the addition of a water box.'
966 % (
967 Options["--waterBox"],
968 OptionsInfo["ForcefieldParams"]["Biopolymer"],
969 OptionsInfo["ForcefieldParams"]["Water"],
970 )
971 )
972
973
974 def ProcessOutPlotParameters():
975 """Process out plot parameters."""
976
977 DefaultValues = {"Type": "line", "Width": 10.0, "Height": 5.6}
978 OptionsInfo["OutPlotParams"] = MiscUtil.ProcessOptionSeabornPlotParameters(
979 "--outPlotParams", Options["--outPlotParams"], DefaultValues
980 )
981 if not re.match("^(linepoint|scatter|Line)$", OptionsInfo["OutPlotParams"]["Type"], re.I):
982 MiscUtil.PrintError(
983 'The value, %s, specified for "type" using option "--outPlotParams" is not supported. Valid plot types: linepoint, scatter or line'
984 % (OptionsInfo["OutPlotParams"]["Type"])
985 )
986
987 for PlotParamName in ["XLabel", "YLabel", "Title"]:
988 if not re.match("^auto$", OptionsInfo["OutPlotParams"][PlotParamName], re.I):
989 MiscUtil.PrintError(
990 'The value, %s, specified for "%s" using option "--outPlotParams" is not supported. Valid value: auto'
991 % (PlotParamName, OptionsInfo["OutPlotParams"][PlotParamName])
992 )
993
994 OptionsInfo["OutPlotInitialized"] = False
995
996
997 def ProcessOptions():
998 """Process and validate command line arguments and options."""
999
1000 MiscUtil.PrintInfo("Processing options...")
1001
1002 ValidateOptions()
1003
1004 OptionsInfo["Infile"] = Options["--infile"]
1005 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
1006 OptionsInfo["InfileRoot"] = FileName
1007 OptionsInfo["InfilePath"] = os.path.abspath(OptionsInfo["Infile"])
1008
1009 SmallMolFile = Options["--smallMolFile"]
1010 SmallMolID = Options["--smallMolID"]
1011 SmallMolFilePath = None
1012 SmallMolFileMode = False
1013 SmallMolFileRoot = None
1014 if SmallMolFile is not None:
1015 FileDir, FileName, FileExt = MiscUtil.ParseFileName(SmallMolFile)
1016 SmallMolFileRoot = FileName
1017 SmallMolFileMode = True
1018 SmallMolFilePath = os.path.abspath(SmallMolFile)
1019
1020 OptionsInfo["SmallMolFile"] = SmallMolFile
1021 OptionsInfo["SmallMolFilePath"] = SmallMolFilePath
1022 OptionsInfo["SmallMolFileRoot"] = SmallMolFileRoot
1023 OptionsInfo["SmallMolFileMode"] = SmallMolFileMode
1024 OptionsInfo["SmallMolID"] = SmallMolID.upper()
1025
1026 ProcessOutfilePrefixOption()
1027 ProcessOutfileDirOption()
1028
1029 ParamsDefaultInfoOverride = {"DataOutType": "Step Speed PotentialEnergy Temperature Time Density Volume"}
1030 ParamsDefaultInfoOverride["DataOutTypePlotX"] = "Time"
1031 ParamsDefaultInfoOverride["DataOutTypePlotY"] = "PotentialEnergy Temperature Density Volume"
1032 for ParamName in [
1033 "PDBOutMinimized",
1034 "PDBOutFinal",
1035 "PDBOutPhase1HeatedNVT",
1036 "PDBOutPhase2AnnealedNVT",
1037 "PDBOutPhase3EquilibratedNVT",
1038 "PDBOutPhase4EquilibratedNPT",
1039 "PDBOutPhase5ProductionNPT",
1040 ]:
1041 ParamsDefaultInfoOverride[ParamName] = True
1042 OptionsInfo["OutputParams"] = OpenMMUtil.ProcessOptionOpenMMOutputParameters(
1043 "--outputParams", Options["--outputParams"], OptionsInfo["OutfilePrefix"], ParamsDefaultInfoOverride
1044 )
1045
1046 ProcessOutPlotParameters()
1047 ProcessOutfileNames()
1048
1049 OptionsInfo["MDProtocolParams"] = OpenMMUtil.ProcessOptionOpenMMMDProtocolParameters(
1050 "-m, --mdProtocolParams", Options["--mdProtocolParams"]
1051 )
1052
1053 OptionsInfo["ForcefieldParams"] = OpenMMUtil.ProcessOptionOpenMMForcefieldParameters(
1054 "--forcefieldParams", Options["--forcefieldParams"]
1055 )
1056
1057 OptionsInfo["FreezeAtoms"] = True if re.match("^yes$", Options["--freezeAtoms"], re.I) else False
1058 if OptionsInfo["FreezeAtoms"]:
1059 OptionsInfo["FreezeAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters(
1060 "--freezeAtomsParams", Options["--freezeAtomsParams"]
1061 )
1062 else:
1063 OptionsInfo["FreezeAtomsParams"] = None
1064
1065 OptionsInfo["OutputReportersMode"] = Options["--outputReportersMode"]
1066 OptionsInfo["OutputReportersModeAllPhases"] = (
1067 True if re.match("^AllPhases$", Options["--outputReportersMode"], re.I) else False
1068 )
1069
1070 ParamsDefaultInfoOverride = {"Name": Options["--platform"], "Threads": 1}
1071 OptionsInfo["PlatformParams"] = OpenMMUtil.ProcessOptionOpenMMPlatformParameters(
1072 "--platformParams", Options["--platformParams"], ParamsDefaultInfoOverride
1073 )
1074
1075 OptionsInfo["RestraintAtoms"] = True if re.match("^yes$", Options["--restraintAtoms"], re.I) else False
1076 if OptionsInfo["RestraintAtoms"]:
1077 OptionsInfo["RestraintAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters(
1078 "--restraintAtomsParams", Options["--restraintAtomsParams"]
1079 )
1080 else:
1081 OptionsInfo["RestraintAtomsParams"] = None
1082 OptionsInfo["RestraintSpringConstant"] = float(Options["--restraintSpringConstant"])
1083
1084 OptionsInfo["SystemParams"] = OpenMMUtil.ProcessOptionOpenMMSystemParameters(
1085 "--systemParams", Options["--systemParams"]
1086 )
1087
1088 OptionsInfo["IntegratorParams"] = OpenMMUtil.ProcessOptionOpenMMIntegratorParameters(
1089 "--integratorParams",
1090 Options["--integratorParams"],
1091 HydrogenMassRepartioningStatus=OptionsInfo["SystemParams"]["HydrogenMassRepartioning"],
1092 )
1093 if OptionsInfo["MDProtocolParams"]["Phase1"]:
1094 OptionsInfo["IntegratorParams"]["Temperature"] = OptionsInfo["MDProtocolParams"]["Phase1InitialStart"]
1095 else:
1096 OptionsInfo["IntegratorParams"]["Temperature"] = OptionsInfo["MDProtocolParams"]["Phase1InitialEnd"]
1097
1098 OptionsInfo["SimulationParams"] = OpenMMUtil.ProcessOptionOpenMMSimulationParameters(
1099 "--simulationParams", Options["--simulationParams"]
1100 )
1101
1102 ProcessWaterBoxParameters()
1103
1104 OptionsInfo["Overwrite"] = Options["--overwrite"]
1105
1106 # Track top level working directory...
1107 OptionsInfo["TopWorkingDir"] = os.getcwd()
1108
1109
1110 def RetrieveOptions():
1111 """Retrieve command line arguments and options."""
1112
1113 # Get options...
1114 global Options
1115 Options = docopt(_docoptUsage_)
1116
1117 # Set current working directory to the specified directory...
1118 WorkingDir = Options["--workingdir"]
1119 if WorkingDir:
1120 os.chdir(WorkingDir)
1121
1122 # Handle examples option...
1123 if "--examples" in Options and Options["--examples"]:
1124 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1125 sys.exit(0)
1126
1127
1128 def ValidateOptions():
1129 """Validate option values."""
1130
1131 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1132 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
1133
1134 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--infile"])
1135 OutfilePrefix = Options["--outfilePrefix"]
1136 if not re.match("^auto$", OutfilePrefix, re.I):
1137 if re.match("^(%s)$" % OutfilePrefix, FileName, re.I):
1138 MiscUtil.PrintError(
1139 'The value specified, %s, for option "--outfilePrefix" is not valid. You must specify a value different from, %s, the root of infile name.'
1140 % (OutfilePrefix, FileName)
1141 )
1142
1143 if Options["--smallMolFile"] is not None:
1144 MiscUtil.ValidateOptionFilePath("-l, --smallMolFile", Options["--smallMolFile"])
1145 MiscUtil.ValidateOptionFileExt("-l, --smallMolFile", Options["--smallMolFile"], "sd sdf")
1146
1147 SmallMolID = Options["--smallMolID"]
1148 if len(SmallMolID) != 3:
1149 MiscUtil.PrintError(
1150 'The value specified, %s, for option "--smallMolID" is not valid. You must specify a three letter small molecule ID.'
1151 % (SmallMolID)
1152 )
1153
1154 MiscUtil.ValidateOptionDirPath("-o, --outfileDir", Options["--outfileDir"])
1155 MiscUtil.ValidateOptionsOutputDirOverwrite(
1156 "-o, --outfileDir", Options["--outfileDir"], "--overwrite", Options["--overwrite"]
1157 )
1158
1159 MiscUtil.ValidateOptionTextValue("--freezeAtoms", Options["--freezeAtoms"], "yes no")
1160 if re.match("^yes$", Options["--freezeAtoms"], re.I):
1161 if Options["--freezeAtomsParams"] is None:
1162 MiscUtil.PrintError(
1163 'No value specified for option "--freezeAtomsParams". You must specify valid values during, yes, value for "--freezeAtoms" option.'
1164 )
1165
1166 MiscUtil.ValidateOptionTextValue(
1167 "--outputReportersMode", Options["--outputReportersMode"], "AllPhases ProductionPhaseOnly"
1168 )
1169
1170 MiscUtil.ValidateOptionTextValue("-p, --platform", Options["--platform"], "CPU CUDA OpenCL Reference")
1171
1172 MiscUtil.ValidateOptionTextValue("--restraintAtoms", Options["--restraintAtoms"], "yes no")
1173 if re.match("^yes$", Options["--restraintAtoms"], re.I):
1174 if Options["--restraintAtomsParams"] is None:
1175 MiscUtil.PrintError(
1176 'No value specified for option "--restraintAtomsParams". You must specify valid values during, yes, value for "--restraintAtoms" option.'
1177 )
1178
1179 MiscUtil.ValidateOptionFloatValue("--restraintSpringConstant", Options["--restraintSpringConstant"], {">": 0})
1180
1181 MiscUtil.ValidateOptionTextValue("--waterBox", Options["--waterBox"], "yes no")
1182
1183
1184 # Setup a usage string for docopt...
1185 _docoptUsage_ = """
1186 OpenMMExecuteMDSimulationProtocol.py - Execute MD simulation workflow
1187
1188 Usage:
1189 OpenMMExecuteMDSimulationProtocol.py [--forcefieldParams <Name,Value,..>] [--freezeAtoms <yes or no>]
1190 [--freezeAtomsParams <Name,Value,..>] [--integratorParams <Name,Value,..>]
1191 [--outfilePrefix <text>] [--outputParams <Name,Value,..>] [--outPlotParams <Name,Value,...>]
1192 [--outputReportersMode <text>] [--overwrite] [--platform <text>] [--mdProtocolParams <Name,Value,..>]
1193 [--platformParams <Name,Value,..>] [--restraintAtoms <yes or no>] [--restraintAtomsParams <Name,Value,..>]
1194 [--restraintSpringConstant <number>] [--simulationParams <Name,Value,..>] [--smallMolFile <SmallMolFile>]
1195 [--smallMolID <text>] [--systemParams <Name,Value,..>] [--waterBox <yes or no>]
1196 [--waterBoxParams <Name,Value,..>] [-w <dir>] -i <infile> -o <outifiledir>
1197 OpenMMExecuteMDSimulationProtocol.py -h | --help | -e | --examples
1198
1199 Description:
1200 Perform a MD simulation using a simulation protocol. You may run a simulation
1201 using a macromolecule or a macromolecule in a complex with small molecule.
1202 By default, the system is minimized before executing the MD simulation protocol.
1203
1204 The MD protocol consists of the following steps:
1205
1206 . Initial heating (NVT simulation)
1207 . Heating and cooling cycles (NVT simulation)
1208 . NVT equilibration
1209 . NPT equilibration
1210 . NPT production
1211
1212 The input file must contain a macromolecule already prepared for simulation.
1213 The preparation of the macromolecule for a simulation generally involves the
1214 following: identification and replacement non-standard residues; addition of
1215 missing residues; addition of missing heavy atoms; addition of missing
1216 hydrogens; addition of a water box which is optional.
1217
1218 In addition, the small molecule input file must contain a molecule already
1219 prepared for simulation. It must contain appropriate 3D coordinates relative
1220 to the macromolecule along with no missing hydrogens.
1221
1222 You may optionally add a water box and freeze/restraint atoms for the
1223 simulation.
1224
1225 The restart option is not available in the current script. You may employ
1226 another script named OpenMMPerformMDSimulation.py to restart the
1227 simulation using the final checkpoint file generated by the current script.
1228
1229 By default, the MD protocol is executed for a total of 8.15 ns as shown
1230 below:
1231
1232 ... ... ...
1233 MD protocol annealing (StepSize: 4 fs)
1234
1235 Phase1 - Initial heating (NVT simulation):
1236
1237 Initial heating (Start: 0.0 K; End: 300.0 K; Change: 5.0 K)
1238 TotalSteps: 305,000; TotalTime: 1.22 ns
1239
1240 Equilibration after initial heating (Steps: 100,000; Time: 400.00 ps)
1241
1242 Phase2 - Heating and cooling cycles (NVT simulation):
1243
1244 Heating and cooling cycles (NumCycles: 1)
1245
1246 Heating and cooling cycle 1
1247 Heating (Start: 300.0 K; End: 315.0 K; Change: 1.0 K)
1248 TotalSteps: 16,000; TotalTime: 64.00 ps
1249
1250 Equilibration after heating (Steps: 100,000; Time: 400.00 ps)
1251
1252 Cooling (Start: 315.0 K; End: 300.0 K; Step: 1.0 K)
1253 TotalSteps: 16,000; TotalTime: 64.00 ps
1254
1255 Equilibration after cooling (Steps: 100,000; Time: 400.00 ps)
1256
1257 Phase3 - NVT equilibration: (Steps: 200,000; Time: 800.00 ps)
1258
1259 Phase4 - NPT equilibration: (Steps: 200,000; Time: 800.00 ps)
1260
1261 Phase5 - NPT production: (Steps: 1,000.000; Time: 4.00 ns)
1262
1263 MD protocol Summary: (TotalSteps: 2,037,000; TotalTime: 8.15 ns)
1264
1265 ... ... ...
1266
1267 The supported macromolecule input file formats are: PDB (.pdb) and
1268 CIF (.cif)
1269
1270 The supported small molecule input file format are : SD (.sdf, .sd)
1271
1272 Possible outfile prefixes:
1273
1274 <InfileRoot>
1275 <InfileRoot>_Solvated
1276 <InfileRoot>_<SmallMolFileRoot>
1277 <InfileRoot>_<SmallMolFileRoot>_Complex_Solvated
1278
1279 Possible output files:
1280
1281 <OutfilePrefix>.<pdb or cif> [ Initial sytem ]
1282 <OutfilePrefix>.<dcd or xtc>
1283
1284 <OutfilePrefix>_Reimaged.<pdb or cif> [ First frame ]
1285 <OutfilePrefix>_Reimaged.<dcd or xtc>
1286
1287 <OutfilePrefix>_Minimized.<pdb or cif>
1288 <OutfilePrefix>_Final.<pdb or cif>
1289
1290 <OutfilePrefix>_NVT_Phase1_Heated_Equilibated.<pdb or cif>
1291 <OutfilePrefix>_NVT_Phase2_Annealed_Equilibrated.<pdb or cif>
1292 <OutfilePrefix>_NVT_Phase3_Equilibrated.<pdb or cif>
1293 <OutfilePrefix>_NPT_Phase4_Equilibrated.<pdb or cif>
1294 <OutfilePrefix>_NPT_Phase5_Production.<pdb or cif>
1295
1296 <OutfilePrefix>.chk
1297 <OutfilePrefix>.csv
1298 <OutfilePrefix>_Minimization.csv
1299 <OutfilePrefix>_FinalState.chk
1300 <OutfilePrefix>_FinalState.xml
1301
1302 <OutfilePrefix>_System.xml
1303 <OutfilePrefix>_Integrator.xml
1304
1305 <OutfilePrefix>_<DataOutTypePlotY1>Plot.<outExt>
1306 <OutfilePrefix>_<DataOutTypePlotY2>Plot.<outExt>
1307 ... ... ...
1308
1309 The reimaged PDB file, <OutfilePrefix>_Reimaged.pdb, corresponds to the first
1310 frame in the trajectory. The reimaged trajectory file contains all the frames
1311 aligned to the first frame after reimaging of the frames for periodic systems.
1312
1313 Options:
1314 -e, --examples
1315 Print examples.
1316 -f, --forcefieldParams <Name,Value,..> [default: auto]
1317 A comma delimited list of parameter name and value pairs for biopolymer,
1318 water, and small molecule forcefields.
1319
1320 The supported parameter names along with their default values are
1321 are shown below:
1322
1323 biopolymer, amber14-all.xml [ Possible values: Any Valid value ]
1324 smallMolecule, openff-2.2.1 [ Possible values: Any Valid value ]
1325 water, auto [ Possible values: Any Valid value ]
1326 additional, none [ Possible values: Space delimited list of any
1327 valid value ]
1328
1329 Possible biopolymer forcefield values:
1330
1331 amber14-all.xml, amber99sb.xml, amber99sbildn.xml, amber03.xml,
1332 amber10.xml
1333 charmm36.xml, charmm_polar_2019.xml
1334 amoeba2018.xml
1335
1336 Possible small molecule forcefield values:
1337
1338 openff-2.2.1, openff-2.0.0, openff-1.3.1, openff-1.2.1,
1339 openff-1.1.1, openff-1.1.0,...
1340 smirnoff99Frosst-1.1.0, smirnoff99Frosst-1.0.9,...
1341 gaff-2.11, gaff-2.1, gaff-1.81, gaff-1.8, gaff-1.4,...
1342
1343 The default water forcefield valus is dependent on the type of the
1344 biopolymer forcefield as shown below:
1345
1346 Amber: amber14/tip3pfb.xml
1347 CHARMM: charmm36/water.xml or None for charmm_polar_2019.xml
1348 Amoeba: None (Explicit)
1349
1350 Possible water forcefield values:
1351
1352 amber14/tip3p.xml, amber14/tip3pfb.xml, amber14/spce.xml,
1353 amber14/tip4pew.xml, amber14/tip4pfb.xml,
1354 charmm36/water.xml, charmm36/tip3p-pme-b.xml,
1355 charmm36/tip3p-pme-f.xml, charmm36/spce.xml,
1356 charmm36/tip4pew.xml, charmm36/tip4p2005.xml,
1357 charmm36/tip5p.xml, charmm36/tip5pew.xml,
1358 implicit/obc2.xml, implicit/GBn.xml, implicit/GBn2.xml,
1359 amoeba2018_gk.xml (Implict water)
1360 None (Explicit water for amoeba)
1361
1362 The additional forcefield value is a space delimited list of any valid
1363 forcefield values and is passed on to the OpenMMForcefields
1364 SystemGenerator along with the specified forcefield values for
1365 biopolymer, water, and mall molecule. Possible additional forcefield
1366 values are:
1367
1368 amber14/DNA.OL15.xml amber14/RNA.OL3.xml
1369 amber14/lipid17.xml amber14/GLYCAM_06j-1.xml
1370 ... ... ...
1371
1372 You may specify any valid forcefield names supported by OpenMM. No
1373 explicit validation is performed.
1374 --freezeAtoms <yes or no> [default: no]
1375 Freeze atoms during a simulation. The specified atoms are kept completely
1376 fixed by setting their masses to zero. Their positions do not change during
1377 local energy minimization and MD simulation, and they do not contribute
1378 to the kinetic energy of the system.
1379 --freezeAtomsParams <Name,Value,..>
1380 A comma delimited list of parameter name and value pairs for freezing
1381 atoms during a simulation. You must specify these parameters for 'yes'
1382 value of '--freezeAtoms' option.
1383
1384 The supported parameter names along with their default values are
1385 are shown below:
1386
1387 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
1388 Protein, Residues, or Water ]
1389 selectionSpec, auto [ Possible values: A space delimited list of
1390 residue names ]
1391 negate, no [ Possible values: yes or no ]
1392
1393 A brief description of parameters is provided below:
1394
1395 selection: Atom selection to freeze.
1396 selectionSpec: A space delimited list of residue names for
1397 selecting atoms to freeze. You must specify its value during
1398 'Ligand' and 'Protein' value for 'selection'. The default values
1399 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
1400 and 'Water' values of 'selection' as shown below:
1401
1402 CAlphaProtein: List of stadard protein residues from pdbfixer
1403 for selecting CAlpha atoms.
1404 Ions: Li Na K Rb Cs Cl Br F I
1405 Water: HOH
1406 Protein: List of standard protein residues from pdbfixer.
1407
1408 negate: Negate atom selection match to select atoms for freezing.
1409
1410 In addition, you may specify an explicit space delimited list of residue
1411 names using 'selectionSpec' for any 'selection". The specified residue
1412 names are appended to the appropriate default values during the
1413 selection of atoms for freezing.
1414 -h, --help
1415 Print this help message.
1416 -i, --infile <infile>
1417 Input file name containing a macromolecule.
1418 --integratorParams <Name,Value,..> [default: auto]
1419 A comma delimited list of parameter name and value pairs for integrator
1420 during a simulation.
1421
1422 The supported parameter names along with their default values are
1423 are shown below:
1424
1425 integrator, LangevinMiddle [ Possible values: LangevinMiddle,
1426 Langevin, NoseHoover, Brownian ]
1427
1428 randomSeed, auto [ Possible values: > 0 ]
1429
1430 frictionCoefficient, 1.0 [ Units: 1/ps ]
1431 stepSize, auto [ Units: fs; Default value: 4 fs during yes value of
1432 hydrogen mass repartioning with no freezing/restraining of atoms;
1433 otherwsie, 2 fs ]
1434
1435 barostat, MonteCarlo [ Possible values: MonteCarlo or
1436 MonteCarloMembrane ]
1437 barostatInterval, 25
1438 pressure, 1.0 [ Units: atm ]
1439
1440 Parameters used only for MonteCarloMembraneBarostat with default
1441 values corresponding to Amber forcefields:
1442
1443 surfaceTension, 0.0 [ Units: atm*A. It is automatically converted
1444 into OpenMM default units of atm*nm before its usage. ]
1445 xymode, Isotropic [ Possible values: Anisotropic or Isotropic ]
1446 zmode, Free [ Possible values: Free or Fixed ]
1447
1448 A brief description of parameters is provided below:
1449
1450 integrator: Type of integrator
1451
1452 randomSeed: Random number seed for barostat and integrator. Not
1453 supported for NoseHoover integrator.
1454
1455 frictionCoefficient: Friction coefficient for coupling the system to
1456 the heat bath..
1457 stepSize: Simulation time step size.
1458
1459 barostat: Barostat type.
1460 barostatInterval: Barostat interval step size during NPT
1461 simulation for applying Monte Carlo pressure changes.
1462 pressure: Pressure during NPT simulation.
1463
1464 Parameters used only for MonteCarloMembraneBarostat:
1465
1466 surfaceTension: Surface tension acting on the system.
1467 xymode: Behavior along X and Y axes. You may allow the X and Y axes
1468 to vary independently of each other or always scale them by the same
1469 amount to keep the ratio of their lengths constant.
1470 zmode: Beahvior along Z axis. You may allow the Z axis to vary
1471 independently of the other axes or keep it fixed.
1472
1473 -m, --mdProtocolParams <Name,Value,..> [default: auto]
1474 A comma delimited list of parameter name and value pairs for executing
1475 MD protocol.
1476
1477 The supported parameter names along with their default values are
1478 are shown below:
1479
1480 Phase1 - Initial heating parameters (NVT simulation):
1481
1482 phase1, yes [ Possible values: yes or no ]
1483 phase1InitialStart, 0.0 [ Units: kelvin ]
1484 phase1InitialEnd, 300.0 [ Units: kelvin ]
1485 phase1InitialChange, 5.0 [ Units: kelvin ]
1486 phase1InitialSteps, 5000
1487
1488 phase1InitialEquilibrationSteps, 100000
1489
1490 Phase2 - Heating and cooling cycles parameters (NVT simulation):
1491
1492 phase2, yes [ Possible values: yes or no ]
1493 phase2Cycles, 1
1494 phase2CycleStart, auto [ Units: kelvin. The default value is set to
1495 initialEnd ]
1496 phase2CycleEnd, 315.0 [ Units: kelvin ]
1497 phase2CycleChange, 1.0 [ Units: kelvin ]
1498 phase2CycleSteps, 1000
1499
1500 phase2CycleEquilibrationSteps, 100000
1501
1502 Phase3 - NVT equilibration parameters:
1503
1504 phase3, yes [ Possible values: yes or no ]
1505 phase3Steps, 200000
1506
1507 Phase4 - NPT equilibration parameters:
1508
1509 phase4, yes [ Possible values: yes or no ]
1510 phase4Steps, 200000
1511
1512 Phase5 - NPT production parameters:
1513
1514 phase5, yes [ Possible values: yes or no ]
1515 phase5Steps, 1000000
1516 phase5StepSize, auto [ Units: fs; Default value: Same as stepSize
1517 parameter in integratorParams option. ]
1518
1519 A brief description of parameters is provided below:
1520
1521 Phase1 - Initial heating parameters (NVT simulation):
1522
1523 phase1: Execute phase1.
1524 phase1InitialStart: Start temperature for initial heating.
1525 phase1InitialEnd: End temperature for initial heating.
1526 phas1InitialChange: Temperature change for increasing temperature
1527 during initial heating.
1528 phase1InitialSteps: Number of simulation steps after each
1529 heating step during initial heating
1530
1531 phase1InitialEquilibrationSteps: Number of equilibration steps
1532 after the completion of initial heating.
1533
1534 Phase2 - Heating and cooling cycles parameters (NVT simulation):
1535
1536 phase2: Execute phase2.
1537 phase2Cycles: Number of annealing cycles to perform. Each cycle
1538 consists of a heating and a cooling phase. The heating phase
1539 consists of the following steps: Heat system from start to
1540 end temperature using change size and perform simulation for a
1541 number of steps after each increase in temperature; Perform
1542 equilibration after the completion of heating. The cooling
1543 phase is reverse of the heating phase and cools the system
1544 from end to start temperature.
1545
1546 phase2CycleStart: Start temperature for annealing cycle.
1547 phase2CycleEnd: End temperature for annealing cycle.
1548 phase2CycleChange: Temperature change for increasing or decreasing
1549 temperature during annealing cycle.
1550 phase2CycleSteps: Number of simulation steps after each heating and
1551 cooling step during annealing cycle.
1552
1553 phase2CycleEquilibrationSteps: Number of equilibration steps
1554 after the completion of heating and cooling phase during a
1555 annealing cycle.
1556
1557 Phase3 - NVT equilibration parameters:
1558
1559 phase3: Execute phase3.
1560 phase3Steps: Number of NVT equilibration steps.
1561
1562 Phase4 - NPT equilibration parameters:
1563
1564 phase4: Execute phase4.
1565 phase4Steps: Number of NPT equilibration steps.
1566
1567 Phase5 - NPT production parameters:
1568
1569 phase5: Execute phase5
1570 phase5Steps: Number of NPT production steps.
1571 phase5StepSize: Simulation time step size for NPT production.
1572
1573 -o, --outfileDir <outfiledir>
1574 Output files directory.
1575 --outfilePrefix <text> [default: auto]
1576 File prefix for generating the names of output files. The default value
1577 depends on the names of input files for macromolecule and small molecule
1578 along with the type of statistical ensemble and the nature of the solvation.
1579
1580 The possible values for outfile prefix are shown below:
1581
1582 <InfileRoot>_<Mode>
1583 <InfileRoot>_Solvated_<Mode>
1584 <InfileRoot>_<SmallMolFileRoot>_Complex
1585 <InfileRoot>_<SmallMolFileRoot>_Complex_Solvated
1586
1587 --outputParams <Name,Value,..> [default: auto]
1588 A comma delimited list of parameter name and value pairs for generating
1589 output during a simulation.
1590
1591 The supported parameter names along with their default values are
1592 are shown below:
1593
1594 checkpoint, no [ Possible values: yes or no ]
1595 checkpointFile, auto [ Default: <OutfilePrefix>.chk ]
1596 checkpointSteps, 10000
1597
1598 dataOutType, auto [ Possible values: A space delimited list of valid
1599 parameter names.
1600 Default: Density Step Speed PotentialEnergy Temperature Time
1601 Volume
1602 Other valid names: ElapsedTime Progress RemainingTime
1603 KineticEnergy TotalEnergy ]
1604
1605 dataLog, yes [ Possible values: yes or no ]
1606 dataLogFile, auto [ Default: <OutfilePrefix>.csv ]
1607 dataLogSteps, 1000
1608
1609 dataStdout, no [ Possible values: yes or no ]
1610 dataStdoutSteps, 1000
1611
1612 dataOutTypePlot, yes [ Possible values: yes or no ]
1613 dataOutTypePlotX, auto [ Default: Time; Possible values: Step or
1614 Time ]
1615 dataOutTypePlotY, auto [ Possible values: A space delimited list
1616 of valid parameter names specified for dataOutType.
1617 Default: Density PotentialEnergy Temperature Volume
1618 Other valid names: KineticEnergy TotalEnergy]
1619
1620 minimizationDataSteps, 100
1621 minimizationDataStdout, no [ Possible values: yes or no ]
1622 minimizationDataLog, no [ Possible values: yes or no ]
1623 minimizationDataLogFile, auto [ Default:
1624 <OutfilePrefix>_MinimizationOut.csv ]
1625 minimizationDataOutType, auto [ Possible values: A space delimited
1626 list of valid parameter names. Default: SystemEnergy
1627 RestraintEnergy MaxConstraintError.
1628 Other valid names: RestraintStrength ]
1629
1630 pdbOutFormat, PDB [ Possible values: PDB or CIF ]
1631 pdbOutKeepIDs, yes [ Possible values: yes or no ]
1632
1633 pdbOutMinimized, yes [ Possible values: yes or no ]
1634 pdbOutFinal, yes [ Possible values: yes or no ]
1635
1636 pdbOutPhase1HeatedNVT, yes [ Possible values: yes or no ]
1637 pdbOutPhase2AnnealedNVT, yes [ Possible values: yes or no ]
1638 pdbOutPhase3EquilibratedNVT, yes [ Possible values: yes or no ]
1639 pdbOutPhase4EquilibratedNPT, yes [ Possible values: yes or no ]
1640 pdbOutPhase5ProductionNPT, yes [ Possible values: yes or no ]
1641
1642 saveFinalStateCheckpoint, yes [ Possible values: yes or no ]
1643 saveFinalStateCheckpointFile, auto [ Default:
1644 <OutfilePrefix>_FinalState.chk ]
1645 saveFinalStateXML, no [ Possible values: yes or no ]
1646 saveFinalStateXMLFile, auto [ Default:
1647 <OutfilePrefix>_FinalState.xml]
1648
1649 traj, yes [ Possible values: yes or no ]
1650 trajFile, auto [ Default: <OutfilePrefix>.<TrajFormat> ]
1651 trajFormat, DCD [ Possible values: DCD or XTC ]
1652 trajSteps, 10000 [ The default value corresponds to 40 ps for step
1653 size of 4 fs. ]
1654
1655 xmlSystemOut, no [ Possible values: yes or no ]
1656 xmlSystemFile, auto [ Default: <OutfilePrefix>_System.xml ]
1657 xmlIntegratorOut, no [ Possible values: yes or no ]
1658 xmlIntegratorFile, auto [ Default: <OutfilePrefix>_Integrator.xml ]
1659
1660 A brief description of parameters is provided below:
1661
1662 checkpoint: Write intermediate checkpoint file.
1663 checkpointFile: Intermediate checkpoint file name.
1664 checkpointSteps: Frequency of writing intermediate checkpoint file.
1665
1666 dataOutType: Type of data to write to stdout and log file.
1667
1668 dataLog: Write data to log file.
1669 dataLogFile: Data log file name.
1670 dataLogSteps: Frequency of writing data to log file.
1671
1672 dataStdout: Write data to stdout.
1673 dataStdoutSteps: Frequency of writing data to stdout.
1674
1675 dataOutTypePlot: Generate plots using data written to log file.
1676 dataOutTypePlotX: Data out type to plot on X axis.
1677 dataOutTypePlotY: Data out types to plot on Y axis. An individual plot
1678 is generated for each pair of X and Y vaues to be plotted.
1679
1680 minimizationDataSteps: Frequency of writing data to stdout
1681 and log file.
1682 minimizationDataStdout: Write data to stdout.
1683 minimizationDataLog: Write data to log file.
1684 minimizationDataLogFile: Data log fie name.
1685 minimizationDataOutType: Type of data to write to stdout
1686 and log file.
1687
1688 saveFinalStateCheckpoint: Save final state checkpoint file.
1689 saveFinalStateCheckpointFile: Name of final state checkpoint file.
1690 saveFinalStateXML: Save final state XML file.
1691 saveFinalStateXMLFile: Name of final state XML file.
1692
1693 pdbOutFormat: Format of output PDB files.
1694 pdbOutKeepIDs: Keep existing chain and residue IDs.
1695
1696 pdbOutMinimized: Write PDB file after minimization.
1697 pdbOutFinal: Write final PDB file.
1698
1699 pdbOutPhase1HeatedNVT: Write out PDB file after initial heatin
1700 pdbOutPhase2AnnealedNVT: Write out PDB file after heating and
1701 cooling.
1702 pdbOutPhase3EquilibratedNVT: Write out PDB file after equilibration.
1703 pdbOutPhase4EquilibratedNPT: Write out PDB file after equilibration.
1704 pdbOutPhase5ProductionNPT: Write out PDB file after production run.
1705
1706 traj: Write out trajectory file.
1707 trajFile: Trajectory file name.
1708 trajFormat: Trajectory file format.
1709 trajSteps: Frequency of writing trajectory file.
1710
1711 xmlSystemOut: Write system XML file.
1712 xmlSystemFile: System XML file name.
1713 xmlIntegratorOut: Write integrator XML file.
1714 xmlIntegratorFile: Integrator XML file name.
1715
1716 --outPlotParams <Name,Value,...> [default: auto]
1717 A comma delimited list of parameter name and value pairs for generating
1718 plots using Seaborn module. The supported parameter names along with their
1719 default values are shown below:
1720
1721 type,linepoint,outExt,svg,width,10,height,5.6,
1722 titleWeight,bold,labelWeight,bold, style,darkgrid,
1723 palette,deep,font,sans-serif,fontScale,1,
1724 context,notebook
1725
1726 Possible values:
1727
1728 type: linepoint, scatter, or line. Both points and lines are drawn
1729 for linepoint plot type.
1730 outExt: Any valid format supported by Python module Matplotlib.
1731 For example: PDF (.pdf), PNG (.png), PS (.ps), SVG (.svg)
1732 titleWeight, labelWeight: Font weight for title and axes labels.
1733 Any valid value.
1734 style: darkgrid, whitegrid, dark, white, ticks
1735 palette: deep, muted, pastel, dark, bright, colorblind
1736 font: Any valid font name
1737 context: paper, notebook, talk, poster, or any valid name
1738
1739 --outputReportersMode <text> [default: ProductionPhaseOnly]
1740 Add output reporters for production phase only or for all phases of the MD
1741 protocol. Possible values: AllPhases or ProductionPhaseOnly. The following
1742 reporters may be added based on the values of correponding output
1743 parameters specified using '--outputParams' option: TrajReporter,
1744 DataLogReporter, DataStdoutReporter, and CheckpointReporter.
1745 --overwrite
1746 Overwrite existing files.
1747 -p, --platform <text> [default: CPU]
1748 Platform to use for running MD simulation. Possible values: CPU, CUDA,
1749 OpenCL, or Reference.
1750 --platformParams <Name,Value,..> [default: auto]
1751 A comma delimited list of parameter name and value pairs to configure
1752 platform for running MD simulation.
1753
1754 The supported parameter names along with their default values for
1755 different platforms are shown below:
1756
1757 CPU:
1758
1759 threads, 1 [ Possible value: >= 0 or auto. The value of 'auto'
1760 or zero implies the use of all available CPUs for threading. ]
1761
1762 CUDA:
1763
1764 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
1765 deterministicForces, auto [ Possible values: yes or no ]
1766 precision, single [ Possible values: single, double, or mix ]
1767 tempDirectory, auto [ Possible value: DirName ]
1768 useBlockingSync, auto [ Possible values: yes or no ]
1769 useCpuPme, auto [ Possible values: yes or no ]
1770
1771 OpenCL:
1772
1773 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
1774 openCLPlatformIndex, auto [ Possible value: Number]
1775 precision, single [ Possible values: single, double, or mix ]
1776 useCpuPme, auto [ Possible values: yes or no ]
1777
1778 A brief description of parameters is provided below:
1779
1780 CPU:
1781
1782 threads: Number of threads to use for simulation.
1783
1784 CUDA:
1785
1786 deviceIndex: Space delimited list of device indices to use for
1787 calculations.
1788 deterministicForces: Generate reproducible results at the cost of a
1789 small decrease in performance.
1790 precision: Number precision to use for calculations.
1791 tempDirectory: Directory name for storing temporary files.
1792 useBlockingSync: Control run-time synchronization between CPU and
1793 GPU.
1794 useCpuPme: Use CPU-based PME implementation.
1795
1796 OpenCL:
1797
1798 deviceIndex: Space delimited list of device indices to use for
1799 simulation.
1800 openCLPlatformIndex: Platform index to use for calculations.
1801 precision: Number precision to use for calculations.
1802 useCpuPme: Use CPU-based PME implementation.
1803
1804 --restraintAtoms <yes or no> [default: no]
1805 Restraint atoms during a simulation. The motion of specified atoms is
1806 restricted by adding a harmonic force that binds them to their starting
1807 positions. The atoms are not completely fixed unlike freezing of atoms.
1808 Their motion, however, is restricted and they are not able to move far away
1809 from their starting positions during local energy minimization and MD
1810 simulation.
1811 --restraintAtomsParams <Name,Value,..>
1812 A comma delimited list of parameter name and value pairs for restraining
1813 atoms during a simulation. You must specify these parameters for 'yes'
1814 value of '--restraintAtoms' option.
1815
1816 The supported parameter names along with their default values are
1817 are shown below:
1818
1819 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
1820 Protein, Residues, or Water ]
1821 selectionSpec, auto [ Possible values: A space delimited list of
1822 residue names ]
1823 negate, no [ Possible values: yes or no ]
1824
1825 A brief description of parameters is provided below:
1826
1827 selection: Atom selection to restraint.
1828 selectionSpec: A space delimited list of residue names for
1829 selecting atoms to restraint. You must specify its value during
1830 'Ligand' and 'Protein' value for 'selection'. The default values
1831 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
1832 and 'Water' values of 'selection' as shown below:
1833
1834 CAlphaProtein: List of stadard protein residues from pdbfixer
1835 for selecting CAlpha atoms.
1836 Ions: Li Na K Rb Cs Cl Br F I
1837 Water: HOH
1838 Protein: List of standard protein residues from pdbfixer.
1839
1840 negate: Negate atom selection match to select atoms for freezing.
1841
1842 In addition, you may specify an explicit space delimited list of residue
1843 names using 'selectionSpec' for any 'selection". The specified residue
1844 names are appended to the appropriate default values during the
1845 selection of atoms for restraining.
1846 --restraintSpringConstant <number> [default: 2.5]
1847 Restraint spring constant for applying external restraint force to restraint
1848 atoms relative to their initial positions during 'yes' value of '--restraintAtoms'
1849 option. Default units: kcal/mol/A**2. The default value, 2.5, corresponds to
1850 1046.0 kjoules/mol/nm**2. The default value is automatically converted into
1851 units of kjoules/mol/nm**2 before its usage.
1852 --simulationParams <Name,Value,..> [default: auto]
1853 A comma delimited list of parameter name and value pairs for simulation.
1854
1855 The supported parameter names along with their default values are
1856 are shown below:
1857
1858 minimization, yes [ Possible values: yes or no ]
1859 minimizationMaxSteps, auto [ Possible values: >= 0. The value of
1860 zero implies until the minimization is converged. ]
1861 minimizationTolerance, 0.24 [ Units: kcal/mol/A. The default value
1862 0.24, corresponds to OpenMM default of value of 10.04
1863 kjoules/mol/nm. It is automatically converted into OpenMM
1864 default units before its usage. ]
1865
1866 A brief description of parameters is provided below:
1867
1868 minimization: Perform minimization before equilibration and
1869 production run.
1870 minimizationMaxSteps: Maximum number of minimization steps. The
1871 value of zero implies until the minimization is converged.
1872 minimizationTolerance: Energy convergence tolerance during
1873 minimization.
1874
1875 -s, --smallMolFile <SmallMolFile>
1876 Small molecule input file name. The macromolecue and small molecule are
1877 merged for simulation and the complex is written out to a PDB file.
1878 --smallMolID <text> [default: LIG]
1879 Three letter small molecule residue ID. The small molecule ID corresponds
1880 to the residue name of the small molecule and is written out to a PDB file
1881 containing the complex.
1882 --systemParams <Name,Value,..> [default: auto]
1883 A comma delimited list of parameter name and value pairs to configure
1884 a system for simulation.
1885
1886 The supported parameter names along with their default values are
1887 are shown below:
1888
1889 constraints, BondsInvolvingHydrogens [ Possible values: None,
1890 WaterOnly, BondsInvolvingHydrogens, AllBonds, or
1891 AnglesInvolvingHydrogens ]
1892 constraintErrorTolerance, 0.000001
1893 ewaldErrorTolerance, 0.0005
1894
1895 nonbondedMethodPeriodic, PME [ Possible values: NoCutoff,
1896 CutoffNonPeriodic, or PME ]
1897 nonbondedMethodNonPeriodic, NoCutoff [ Possible values:
1898 NoCutoff or CutoffNonPeriodic]
1899 nonbondedCutoff, 1.0 [ Units: nm ]
1900
1901 hydrogenMassRepartioning, yes [ Possible values: yes or no ]
1902 hydrogenMass, 1.5 [ Units: amu]
1903
1904 removeCMMotion, yes [ Possible values: yes or no ]
1905 rigidWater, auto [ Possible values: yes or no. Default: 'No' for
1906 'None' value of constraints; Otherwise, yes ]
1907
1908 A brief description of parameters is provided below:
1909
1910 constraints: Type of system constraints to use for simulation. These
1911 constraints are different from freezing and restraining of any
1912 atoms in the system.
1913 constraintErrorTolerance: Distance tolerance for constraints as a
1914 fraction of the constrained distance.
1915 ewaldErrorTolerance: Ewald error tolerance for a periodic system.
1916
1917 nonbondedMethodPeriodic: Nonbonded method to use during the
1918 calculation of long range interactions for a periodic system.
1919 nonbondedMethodNonPeriodic: Nonbonded method to use during the
1920 calculation of long range interactions for a non-periodic system.
1921 nonbondedCutoff: Cutoff distance to use for long range interactions
1922 in both perioidic non-periodic systems.
1923
1924 hydrogenMassRepartioning: Use hydrogen mass repartioning. It
1925 increases the mass of the hydrogen atoms attached to the heavy
1926 atoms and decreasing the mass of the bonded heavy atom to
1927 maintain constant system mass. This allows the use of larger
1928 integration step size (4 fs) during a simulation.
1929 hydrogenMass: Hydrogen mass to use during repartioning.
1930
1931 removeCMMotion: Remove all center of mass motion at every time step.
1932 rigidWater: Keep water rigid during a simulation. This is determined
1933 automatically based on the value of 'constraints' parameter.
1934
1935 --waterBox <yes or no> [default: no]
1936 Add water box.
1937 --waterBoxParams <Name,Value,..> [default: auto]
1938 A comma delimited list of parameter name and value pairs for adding
1939 a water box.
1940
1941 The supported parameter names along with their default values are
1942 are shown below:
1943
1944 model, tip3p [ Possible values: tip3p, spce, tip4pew, tip5p or
1945 swm4ndp ]
1946 mode, Padding [ Possible values: Size or Padding ]
1947 padding, 1.0
1948 size, None [ Possible value: xsize ysize zsize ]
1949 shape, cube [ Possible values: cube, dodecahedron, or octahedron ]
1950 ionPositive, Na+ [ Possible values: Li+, Na+, K+, Rb+, or Cs+ ]
1951 ionNegative, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
1952 ionicStrength, 0.0
1953
1954 A brief description of parameters is provided below:
1955
1956 model: Water model to use for adding water box. The van der
1957 Waals radii and atomic charges are determined using the
1958 specified water forcefield. You must specify an appropriate
1959 water forcefield. No validation is performed.
1960 mode: Specify the size of the waterbox explicitly or calculate it
1961 automatically for a macromolecule along with adding padding
1962 around ther macromolecule.
1963 padding: Padding around a macromolecule in nanometers for filling
1964 box with water. It must be specified during 'Padding' value of
1965 'mode' parameter.
1966 size: A space delimited triplet of values corresponding to water
1967 size in nanometers. It must be specified during 'Size' value of
1968 'mode' parameter.
1969 ionPositive: Type of positive ion to add during the addition of a
1970 water box.
1971 ionNegative: Type of negative ion to add during the addition of a
1972 water box.
1973 ionicStrength: Total concentration of both positive and negative
1974 ions to add excluding the ions added to neutralize the system
1975 during the addition of a water box.
1976
1977 -w, --workingdir <dir>
1978 Location of working directory which defaults to the current directory.
1979
1980 Examples:
1981 To execute MD simulation protocol for a macromolecule in a PDB file, applying
1982 system constraints for bonds involving hydrogens along with hydrogen mass
1983 repartioning, using a step size of 4 fs, performing minimization until it's
1984 converged, performing phase1 initial heating along for 305,000 steps (1.22 ns)
1985 along with an equilibration for 100,000 steps (400.00 ps) after the completion
1986 of iniital heating, performing phase 2one heating and cooling cycle along with
1987 equilibration for 232,000 steps (928 ps), performing phase3 NVT equilibration
1988 for 200,000 steps (800.00 ps), performing phase 4NPT equilibration for 200,000
1989 steps (800.00 ps), performing phase NPT production run for 1,000.000 steps
1990 (4.00 ns), writing trajectory and data log files every 10,000 steps (40 ps) and
1991 1,000 steps (4 ps) only during the production run, generating a checkpoint file
1992 after the completion of the calculation, and generating various PDB files for the
1993 system during the calculation, type:
1994
1995 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
1996 -o Sample13OutMDProtocol --waterBox yes
1997
1998 To run the first example for performing OpenMM simulation using multi-
1999 threading employing all available CPUs on your machine and generate various
2000 output files, type:
2001
2002 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2003 -o Sample13OutMDProtocol --waterBox yes
2004 --platformParams "threads,0"
2005
2006 To run the first example for performing OpenMM simulation using CUDA platform
2007 on your machine and generate various output files, type:
2008
2009 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2010 -o Sample13OutMDProtocol --waterBox yes
2011 -p CUDA
2012
2013 To run the second example for a marcomolecule in a complex with a small
2014 molecule and generate various output files, type:
2015
2016 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2017 -o Sample13OutMDProtocol --waterBox yes
2018 -s Sample13Ligand.sdf
2019 --platformParams "threads,0"
2020
2021 To run the second example to reporters for writing trajectory, data log, and
2022 checkpoint files during all phases of the execution of MD protocol, and
2023 generate various output files, type:
2024
2025 To run the second example by skipping phase 2 heating and cooling cycle and
2026 generate various output files, type:
2027
2028 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2029 -o Sample13OutMDProtocol --waterBox yes
2030 -s Sample13Ligand.sdf
2031 --platformParams "threads,0"
2032 --outputReportersMode AllPhases
2033
2034 To run the second example by freezing CAlpha atoms in a macromolecule without
2035 using any system constraints to avoid any issues with the freezing of the same atoms,
2036 using a step size of 2 fs, and generate various output files, type:
2037
2038 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2039 -o Sample13OutMDProtocol --waterBox yes
2040 --freezeAtoms yes --freezeAtomsParams "selection,CAlphaProtein"
2041 --systemParams "constraints, None"
2042 --platformParams "threads,0" --integratorParams "stepSize,2"
2043
2044 To run the second example by restrainting CAlpha atoms in a macromolecule and
2045 and generate various output files, type:
2046
2047 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2048 -o Sample13OutMDProtocol --waterBox yes --restraintAtoms yes
2049 --restraintAtomsParams "selection,CAlphaProtein"
2050 --platformParams "threads,0" --integratorParams "stepSize,2"
2051
2052 To run the second example by specifying explict values for various parametres
2053 and generate various output files, type:
2054
2055 % OpenMMExecuteMDSimulationProtocol.py -i Sample13.pdb
2056 -o Sample13OutMDProtocol --waterBox yes
2057 --mdProtocolParams "phase1, yes, phase1InitialStart, 0.0,
2058 phase1InitialEnd, 300.0, phase1InitialChange, 5.0,
2059 phase1InitialSteps, 5000, phase1InitialEquilibrationSteps, 100000,
2060 phase2, yes,phase2Cycles, 1, phase2CycleStart, auto,
2061 phase2CycleEnd, 315.0, phase2CycleSteps, 1000,
2062 phase2CycleEquilibrationSteps, 100000, phase3, yes,
2063 phase3Steps, 200000, phase4, yes, phase4Steps, 200000,
2064 phase5, yes, phase5Steps, 1000000"
2065 -f " biopolymer,amber14-all.xml,smallMolecule, openff-2.2.1,
2066 water,amber14/tip3pfb.xml"
2067 --integratorParams "integrator,LangevinMiddle,randomSeed,42,
2068 stepSize,2,pressure, 1.0"
2069 --outputParams "checkpoint,yes,dataLog,yes,dataStdout,yes,
2070 minimizationDataStdout,yes,minimizationDataLog,yes,
2071 pdbOutFormat,CIF,pdbOutKeepIDs,yes,saveFinalStateCheckpoint, yes,
2072 traj,yes,xmlSystemOut,yes,xmlIntegratorOut,yes"
2073 -p CPU --platformParams "threads,0"
2074 --simulationParams "minimization,yes, minimizationMaxSteps,
2075 5000,equilibration,yes"
2076 --systemParams "constraints,BondsInvolvingHydrogens,
2077 nonbondedMethodPeriodic,PME,nonbondedMethodNonPeriodic,NoCutoff,
2078 hydrogenMassRepartioning, yes"
2079
2080 Author:
2081 Manish Sud(msud@san.rr.com)
2082
2083 Acknowledgment:
2084 Paul Charifson
2085
2086 See also:
2087 OpenMMPrepareMacromolecule.py, OpenMMPerformMDSimulation.py,
2088 OpenMMPerformSimulatedAnnealing.py, OpenMMPerformMinimization.py
2089
2090 Copyright:
2091 Copyright (C) 2026 Manish Sud. All rights reserved.
2092
2093 The functionality available in this script is implemented using OpenMM, an
2094 open source molecuar simulation package.
2095
2096 This file is part of MayaChemTools.
2097
2098 MayaChemTools is free software; you can redistribute it and/or modify it under
2099 the terms of the GNU Lesser General Public License as published by the Free
2100 Software Foundation; either version 3 of the License, or (at your option) any
2101 later version.
2102
2103 """
2104
2105 if __name__ == "__main__":
2106 main()
