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