1 #!/bin/env python
2 #
3 # File: OpenMMPerformMDSimulation.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Copyright (C) 2025 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 # Add local python path to the global path and import standard library modules...
32 import os
33 import sys; sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
34 import time
35 import re
36
37 # OpenMM imports...
38 try:
39 import openmm as mm
40 import openmm.app
41 except ImportError as ErrMsg:
42 sys.stderr.write("\nFailed to import OpenMM related module/package: %s\n" % ErrMsg)
43 sys.stderr.write("Check/update your OpenMM environment and try again.\n\n")
44 sys.exit(1)
45
46 # MDTraj import...
47 try:
48 import mdtraj
49 except ImportError as ErrMsg:
50 sys.stderr.write("\nFailed to import MDTraj: %s\n" % ErrMsg)
51 sys.stderr.write("Check/update your MDTraj environment and try again.\n\n")
52 sys.exit(1)
53
54 # MayaChemTools imports...
55 try:
56 from docopt import docopt
57 import MiscUtil
58 import OpenMMUtil
59 except ImportError as ErrMsg:
60 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
61 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
62 sys.exit(1)
63
64 ScriptName = os.path.basename(sys.argv[0])
65 Options = {}
66 OptionsInfo = {}
67
68 def main():
69 """Start execution of the script."""
70
71 MiscUtil.PrintInfo("\n%s (OpenMM v%s; MayaChemTools v%s; %s): Starting...\n" % (ScriptName, mm.Platform.getOpenMMVersion(), MiscUtil.GetMayaChemToolsVersion(), time.asctime()))
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 def PerformMDSimulation():
88 """Perform MD simulation."""
89
90 # Prepare system for simulation...
91 System, Topology, Positions = PrepareSystem()
92
93 # Freeze and restraint atoms...
94 FreezeRestraintAtoms(System, Topology, Positions)
95
96 # Setup integrator...
97 Integrator = SetupIntegrator()
98
99 # Setup simulation...
100 Simulation = SetupSimulation(System, Integrator, Topology, Positions)
101
102 # Write setup files...
103 WriteSimulationSetupFiles(System, Integrator)
104
105 # Perform minimization...
106 PerformMinimization(Simulation)
107
108 # Set up intial velocities...
109 SetupInitialVelocities(Simulation)
110
111 # Perform equilibration...
112 PerformEquilibration(Simulation)
113
114 # Setup reporters for production run...
115 SetupReporters(Simulation)
116
117 # Perform or restart production run...
118 PerformOrRestartProductionRun(Simulation)
119
120 # Save final state files...
121 WriteFinalStateFiles(Simulation)
122
123 # Reimage and realign trajectory for periodic systems...
124 ProcessTrajectory(System, Topology)
125
126 def PrepareSystem():
127 """Prepare system for simulation."""
128
129 System, Topology, Positions = OpenMMUtil.InitializeSystem(OptionsInfo["Infile"], OptionsInfo["ForcefieldParams"], OptionsInfo["SystemParams"], OptionsInfo["WaterBox"], OptionsInfo["WaterBoxParams"], OptionsInfo["SmallMolFile"], OptionsInfo["SmallMolID"])
130
131 if OptionsInfo["NPTMode"]:
132 if not OpenMMUtil.DoesSystemUsesPeriodicBoundaryConditions(System):
133 MiscUtil.PrintInfo("")
134 MiscUtil.PrintWarning("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. ")
135
136 BarostatHandle = OpenMMUtil.InitializeBarostat(OptionsInfo["IntegratorParams"])
137 MiscUtil.PrintInfo("Adding barostat for NPT simulation...")
138 try:
139 System.addForce(BarostatHandle)
140 except Exception as ErrMsg:
141 MiscUtil.PrintInfo("")
142 MiscUtil.PrintError("Failed to add barostat:\n%s\n" % (ErrMsg))
143
144 # Write out a PDB file for the system...
145 PDBFile = OptionsInfo["PDBOutfile"]
146 if OptionsInfo["RestartMode"]:
147 MiscUtil.PrintInfo("\nSkipping wrriting of PDB file %s during restart..." % PDBFile)
148 else:
149 MiscUtil.PrintInfo("\nWriting PDB file %s..." % PDBFile)
150 OpenMMUtil.WritePDBFile(PDBFile, Topology, Positions, OptionsInfo["OutputParams"]["PDBOutKeepIDs"])
151
152 return (System, Topology, Positions)
153
154 def SetupIntegrator():
155 """Setup integrator. """
156
157 Integrator = OpenMMUtil.InitializeIntegrator(OptionsInfo["IntegratorParams"], OptionsInfo["SystemParams"]["ConstraintErrorTolerance"])
158
159 return Integrator
160
161 def SetupSimulation(System, Integrator, Topology, Positions):
162 """Setup simulation. """
163
164 Simulation = OpenMMUtil.InitializeSimulation(System, Integrator, Topology, Positions, OptionsInfo["PlatformParams"])
165
166 return Simulation
167
168 def SetupInitialVelocities(Simulation):
169 """Setup initial velocities."""
170
171 # Set velocities to random values choosen from a Boltzman distribution at a given
172 # temperature...
173 if OptionsInfo["RestartMode"]:
174 MiscUtil.PrintInfo("\nSkipping setting of intial velocities to temperature during restart...")
175 else:
176 MiscUtil.PrintInfo("\nSetting intial velocities to temperature...")
177 IntegratorParamsInfo = OpenMMUtil.SetupIntegratorParameters(OptionsInfo["IntegratorParams"])
178 Simulation.context.setVelocitiesToTemperature(IntegratorParamsInfo["Temperature"])
179
180 def PerformMinimization(Simulation):
181 """Perform minimization."""
182
183 SimulationParams = OpenMMUtil.SetupSimulationParameters(OptionsInfo["SimulationParams"])
184
185 if OptionsInfo["RestartMode"]:
186 MiscUtil.PrintInfo("\nSkipping energy minimization during restart...")
187 return
188 else:
189 if not SimulationParams["Minimization"]:
190 MiscUtil.PrintInfo("\nSkipping energy minimization...")
191 return
192
193 OutputParams = OptionsInfo["OutputParams"]
194
195 # Setup a local minimization reporter...
196 MinimizeReporter = None
197 if OutputParams["MinimizationDataStdout"] or OutputParams["MinimizationDataLog"]:
198 MinimizeReporter = LocalMinimizationReporter()
199
200 if MinimizeReporter is not None:
201 MiscUtil.PrintInfo("\nAdding minimization reporters...")
202 if OutputParams["MinimizationDataLog"]:
203 MiscUtil.PrintInfo("Adding data log minimization reporter (Steps: %s; File: %s)..." % (OutputParams["MinimizationDataSteps"], OutputParams["MinimizationDataLogFile"]))
204 if OutputParams["MinimizationDataStdout"]:
205 MiscUtil.PrintInfo("Adding data stdout minimization reporter (Steps: %s)..." % (OutputParams["MinimizationDataSteps"]))
206 else:
207 MiscUtil.PrintInfo("\nSkipping addition of minimization reporters...")
208
209 MaxSteps = SimulationParams["MinimizationMaxSteps"]
210
211 MaxStepsMsg = "MaxSteps: %s" % ("UntilConverged" if MaxSteps == 0 else MaxSteps)
212 ToleranceMsg = "Tolerance: %.2f kcal/mol/A (%.2f kjoules/mol/nm)" % (SimulationParams["MinimizationToleranceInKcal"], SimulationParams["MinimizationToleranceInJoules"])
213
214 MiscUtil.PrintInfo("\nPerforming energy minimization (%s; %s)..." % (MaxStepsMsg, ToleranceMsg))
215
216 if OutputParams["MinimizationDataStdout"]:
217 HeaderLine = SetupMinimizationDataOutHeaderLine()
218 print("\n%s" % HeaderLine)
219
220 Simulation.minimizeEnergy(tolerance = SimulationParams["MinimizationTolerance"], maxIterations = MaxSteps, reporter = MinimizeReporter)
221
222 if OutputParams["MinimizationDataLog"]:
223 WriteMinimizationDataLogFile(MinimizeReporter.DataOutValues)
224
225 if OutputParams["PDBOutMinimized"]:
226 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["MinimizedPDBOutfile"])
227 OpenMMUtil.WriteSimulationStatePDBFile(Simulation, OptionsInfo["MinimizedPDBOutfile"], OutputParams["PDBOutKeepIDs"])
228
229 def PerformEquilibration(Simulation):
230 """Perform equilibration."""
231
232 Mode = OptionsInfo["Mode"]
233 SimulationParams = OptionsInfo["SimulationParams"]
234 OutputParams = OptionsInfo["OutputParams"]
235
236 if OptionsInfo["RestartMode"]:
237 MiscUtil.PrintInfo("\nSkipping equilibration during restart...")
238 return
239 else:
240 if not SimulationParams["Equilibration"]:
241 MiscUtil.PrintInfo("\nSkipping equilibration...")
242 return
243
244 EquilibrationSteps = SimulationParams["EquilibrationSteps"]
245
246 IntegratorParamsInfo = OpenMMUtil.SetupIntegratorParameters(OptionsInfo["IntegratorParams"])
247 StepSize = IntegratorParamsInfo["StepSize"]
248
249 TotalTime = OpenMMUtil.GetFormattedTotalSimulationTime(StepSize, EquilibrationSteps)
250 MiscUtil.PrintInfo("\nPerforming equilibration (Mode: %s; Steps: %s; StepSize: %s; TotalTime: %s)..." % (Mode, EquilibrationSteps, StepSize, TotalTime))
251
252 # Equilibrate...
253 Simulation.step(EquilibrationSteps)
254
255 if OutputParams["PDBOutEquilibrated"]:
256 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["EquilibratedPDBOutfile"])
257 OpenMMUtil.WriteSimulationStatePDBFile(Simulation, OptionsInfo["EquilibratedPDBOutfile"], OutputParams["PDBOutKeepIDs"])
258
259 def PerformOrRestartProductionRun(Simulation):
260 """Perform or restart production run. """
261
262 if OptionsInfo["RestartMode"]:
263 RestartProductionRun(Simulation)
264 else:
265 PerformProductionRun(Simulation)
266
267 def PerformProductionRun(Simulation):
268 """Perform production run."""
269
270 Mode = OptionsInfo["Mode"]
271 SimulationParams = OptionsInfo["SimulationParams"]
272 Steps = SimulationParams["Steps"]
273
274 IntegratorParamsInfo = OpenMMUtil.SetupIntegratorParameters(OptionsInfo["IntegratorParams"])
275 StepSize = IntegratorParamsInfo["StepSize"]
276
277 if SimulationParams["Equilibration"]:
278 Simulation.currentStep = 0
279 Simulation.context.setTime(0.0 * mm.unit.picoseconds)
280 MiscUtil.PrintInfo("\nSetting current step and current time to 0 for starting production run after equilibration...")
281
282 TotalTime = OpenMMUtil.GetFormattedTotalSimulationTime(StepSize, Steps)
283 MiscUtil.PrintInfo("\nPerforming production run (Mode: %s; Steps: %s; StepSize: %s; TotalTime: %s)..." % (Mode, Steps, StepSize, TotalTime))
284
285 Simulation.step(Steps)
286
287 def RestartProductionRun(Simulation):
288 """Restart production run."""
289
290 SimulationParams = OptionsInfo["SimulationParams"]
291 RestartParams = OptionsInfo["RestartParams"]
292
293 RestartFile = RestartParams["FinalStateFile"]
294
295 Steps = SimulationParams["Steps"]
296
297 IntegratorParamsInfo = OpenMMUtil.SetupIntegratorParameters(OptionsInfo["IntegratorParams"])
298 StepSize = IntegratorParamsInfo["StepSize"]
299
300 TotalTime = OpenMMUtil.GetFormattedTotalSimulationTime(StepSize, Steps)
301 MiscUtil.PrintInfo("\nRestarting production run (Steps: %s; StepSize: %s; TotalTime: %s)..." % (Steps, StepSize, TotalTime))
302
303 if RestartParams["FinalStateFileCheckpointMode"]:
304 MiscUtil.PrintInfo("Loading final state checkpoint file %s...\n" % RestartFile)
305 Simulation.loadCheckpoint(RestartFile)
306 elif RestartParams["FinalStateFileXMLMode"]:
307 MiscUtil.PrintInfo("Loading final state XML f ile %s...\n" % RestartFile)
308 Simulation.loadState(RestartFile)
309 else:
310 MiscUtil.PrintError("The specified final state restart file, %s, is not a valid. Supported file formats: chk or xml" % (RestartFile))
311
312 Simulation.step(Steps)
313
314 def SetupReporters(Simulation):
315 """Setup reporters. """
316
317 (TrajReporter, DataLogReporter, DataStdoutReporter, CheckpointReporter) = OpenMMUtil.InitializeReporters(OptionsInfo["OutputParams"], OptionsInfo["SimulationParams"]["Steps"], OptionsInfo["DataOutAppendMode"])
318
319 if TrajReporter is None and DataLogReporter is None and DataStdoutReporter is None and CheckpointReporter is None:
320 MiscUtil.PrintInfo("\nSkip adding reporters...")
321 return
322
323 MiscUtil.PrintInfo("\nAdding reporters...")
324
325 OutputParams = OptionsInfo["OutputParams"]
326 AppendMsg = ""
327 if OptionsInfo["RestartMode"]:
328 AppendMsg = "; Append: Yes" if OptionsInfo["DataOutAppendMode"] else "; Append: No"
329 if TrajReporter is not None:
330 MiscUtil.PrintInfo("Adding trajectory reporter (Steps: %s; File: %s%s)..." % (OutputParams["TrajSteps"], OutputParams["TrajFile"], AppendMsg))
331 Simulation.reporters.append(TrajReporter)
332
333 if CheckpointReporter is not None:
334 MiscUtil.PrintInfo("Adding checkpoint reporter (Steps: %s; File: %s)..." % (OutputParams["CheckpointSteps"], OutputParams["CheckpointFile"]))
335 Simulation.reporters.append(CheckpointReporter)
336
337 if DataLogReporter is not None:
338 MiscUtil.PrintInfo("Adding data log reporter (Steps: %s; File: %s%s)..." % (OutputParams["DataLogSteps"], OutputParams["DataLogFile"], AppendMsg))
339 Simulation.reporters.append(DataLogReporter)
340
341 if DataStdoutReporter is not None:
342 MiscUtil.PrintInfo("Adding data stdout reporter (Steps: %s)..." % (OutputParams["DataStdoutSteps"]))
343 Simulation.reporters.append(DataStdoutReporter)
344
345 class LocalMinimizationReporter(mm.MinimizationReporter):
346 """Setup a local minimization reporter. """
347
348 (DataSteps, DataOutTypeList, DataOutDelimiter, StdoutStatus) = [None] * 4
349
350 DataOutValues = []
351 First = True
352
353 def report(self, Iteration, PositonsList, GradientsList, DataStatisticsMap):
354 """Report and track minimization."""
355
356 if self.First:
357 # Initialize...
358 self.DataSteps = OptionsInfo["OutputParams"]["MinimizationDataSteps"]
359 self.DataOutTypeList = OptionsInfo["OutputParams"]["MinimizationDataOutTypeOpenMMNameList"]
360 self.DataOutDelimiter = OptionsInfo["OutputParams"]["DataOutDelimiter"]
361 self.StdoutStatus = True if OptionsInfo["OutputParams"]["MinimizationDataStdout"] else False
362
363 self.First = False
364
365 if Iteration % self.DataSteps == 0:
366 # Setup data values...
367 DataValues = []
368 DataValues.append("%s" % Iteration)
369 for DataType in self.DataOutTypeList:
370 DataValue = "%.4f" % DataStatisticsMap[DataType]
371 DataValues.append(DataValue)
372
373 # Track data...
374 self.DataOutValues.append(DataValues)
375
376 # Print data values...
377 if self.StdoutStatus:
378 print("%s" % self.DataOutDelimiter.join(DataValues))
379
380 # This method must return a bool. You may return true for early termination.
381 return False
382
383 def WriteMinimizationDataLogFile(DataOutValues):
384 """Write minimization data log file. """
385
386 OutputParams = OptionsInfo["OutputParams"]
387
388 Outfile = OutputParams["MinimizationDataLogFile"]
389 OutDelimiter = OutputParams["DataOutDelimiter"]
390
391 MiscUtil.PrintInfo("\nWriting minimization log file %s..." % Outfile)
392 OutFH = open(Outfile, "w")
393
394 HeaderLine = SetupMinimizationDataOutHeaderLine()
395 OutFH.write("%s\n" % HeaderLine)
396
397 for LineWords in DataOutValues:
398 Line = OutDelimiter.join(LineWords)
399 OutFH.write("%s\n" % Line)
400
401 OutFH.close()
402
403 def SetupMinimizationDataOutHeaderLine():
404 """Setup minimization data output header line. """
405
406 LineWords = ["Iteration"]
407 for Label in OptionsInfo["OutputParams"]["MinimizationDataOutTypeList"]:
408 if re.match("^(SystemEnergy|RestraintEnergy)$", Label, re.I):
409 LineWords.append("%s(kjoules/mol)" % Label)
410 elif re.match("^RestraintStrength$", Label, re.I):
411 LineWords.append("%s(kjoules/mol/nm^2)" % Label)
412 else:
413 LineWords.append(Label)
414
415 Line = OptionsInfo["OutputParams"]["DataOutDelimiter"].join(LineWords)
416
417 return Line
418
419 def FreezeRestraintAtoms(System, Topology, Positions):
420 """Handle freezing and restraining of atoms."""
421
422 FreezeAtomList, RestraintAtomList = OpenMMUtil.ValidateAndFreezeRestraintAtoms(OptionsInfo["FreezeAtoms"], OptionsInfo["FreezeAtomsParams"], OptionsInfo["RestraintAtoms"], OptionsInfo["RestraintAtomsParams"], OptionsInfo["RestraintSpringConstant"], OptionsInfo["SystemParams"], System, Topology, Positions)
423
424 # Check and adjust step size...
425 if FreezeAtomList is not None or RestraintAtomList is not None:
426 if re.match("^auto$", OptionsInfo["IntegratorParams"]["StepSizeSpecified"], re.I):
427 # Automatically set stepSize to 2.0 fs..
428 OptionsInfo["IntegratorParams"]["StepSize"] = 2.0
429 MiscUtil.PrintInfo("")
430 MiscUtil.PrintWarning("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." % (OptionsInfo["IntegratorParams"]["StepSize"]))
431 elif OptionsInfo["IntegratorParams"]["StepSize"] > 2:
432 MiscUtil.PrintInfo("")
433 MiscUtil.PrintWarning("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." % (OptionsInfo["IntegratorParams"]["StepSize"] ))
434 MiscUtil.PrintInfo("")
435
436 def WriteSimulationSetupFiles(System, Integrator):
437 """Write simulation setup files for system and integrator."""
438
439 OutputParams = OptionsInfo["OutputParams"]
440
441 if OutputParams["XmlSystemOut"] or OutputParams["XmlIntegratorOut"]:
442 MiscUtil.PrintInfo("")
443
444 if OutputParams["XmlSystemOut"]:
445 Outfile = OutputParams["XmlSystemFile"]
446 MiscUtil.PrintInfo("Writing system setup XML file %s..." % Outfile)
447 with open(Outfile, mode = "w") as OutFH:
448 OutFH.write(mm.XmlSerializer.serialize(System))
449
450 if OutputParams["XmlIntegratorOut"]:
451 Outfile = OutputParams["XmlIntegratorFile"]
452 MiscUtil.PrintInfo("Writing integrator setup XML file %s..." % Outfile)
453 with open(Outfile, mode = "w") as OutFH:
454 OutFH.write(mm.XmlSerializer.serialize(Integrator))
455
456 def WriteFinalStateFiles(Simulation):
457 """Write final state files. """
458
459 OutputParams = OptionsInfo["OutputParams"]
460
461 if OutputParams["SaveFinalStateCheckpoint"] or OutputParams["SaveFinalStateXML"] or OutputParams["PDBOutFinal"]:
462 MiscUtil.PrintInfo("")
463
464 if OutputParams["SaveFinalStateCheckpoint"]:
465 Outfile = OutputParams["SaveFinalStateCheckpointFile"]
466 MiscUtil.PrintInfo("Writing final state checkpoint file %s..." % Outfile)
467 Simulation.saveCheckpoint(Outfile)
468
469 if OutputParams["SaveFinalStateXML"]:
470 Outfile = OutputParams["SaveFinalStateXMLFile"]
471 MiscUtil.PrintInfo("Writing final state XML file %s..." % Outfile)
472 Simulation.saveState(Outfile)
473
474 if OutputParams["PDBOutFinal"]:
475 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["FinalPDBOutfile"])
476 OpenMMUtil.WriteSimulationStatePDBFile(Simulation, OptionsInfo["FinalPDBOutfile"], OutputParams["PDBOutKeepIDs"])
477
478 def ProcessTrajectory(System, Topology):
479 """Reimage and realign trajectory for periodic systems. """
480
481 ReimageFrames = True if OpenMMUtil.DoesSystemUsesPeriodicBoundaryConditions(System) else False
482 if not ReimageFrames:
483 MiscUtil.PrintInfo("\nSkipping reimaging and realigning of trajectory for a system not using periodic boundary conditions...")
484 return
485
486 MiscUtil.PrintInfo("\nReimaging and realigning trajectory for a system using periodic boundary conditions...")
487
488 # Reimage and realign trajectory file...
489 RealignFrames = True
490 TrajTopologyFile = OptionsInfo["PDBOutfile"]
491 TrajFile = OptionsInfo["OutputParams"]["TrajFile"]
492 Traj, ReimagedStatus, RealignedStatus = OpenMMUtil.ReimageRealignTrajectory(Topology, TrajFile, ReimageFrames, RealignFrames)
493
494 if (Traj is None) or (not ReimagedStatus and not RealignedStatus):
495 MiscUtil.PrintInfo("Skipping writing first frame to PDB file %s..." % PDBOutfile)
496 MiscUtil.PrintInfo("Skippig writing trajectory file %s..." % TrajOutfile)
497 return
498
499 PDBOutFormat = OptionsInfo["OutputParams"]["PDBOutFormat"]
500 PDBOutfile = OptionsInfo["ReimagedPDBOutfile"]
501 TrajOutfile = OptionsInfo["ReimagedTrajOutfile"]
502
503 # Write out first frame...
504 MiscUtil.PrintInfo("Writing first frame to PDB file %s..." % PDBOutfile)
505 if re.match("^CIF$", PDBOutFormat, re.I):
506 # MDTraj doesn't appear to support CIF format. Write it out as a temporary
507 # PDB file and convert it using OpenMM...
508 FileDir, FileRoot, FileExt = MiscUtil.ParseFileName(PDBOutfile)
509 TmpPDBOutfile = "%s_PID%s.pdb" % (FileRoot, os.getpid())
510 Traj[0].save(TmpPDBOutfile)
511
512 PDBHandle = OpenMMUtil.ReadPDBFile(TmpPDBOutfile)
513 OpenMMUtil.WritePDBFile(PDBOutfile, PDBHandle.topology, PDBHandle.positions, OptionsInfo["OutputParams"]["PDBOutKeepIDs"])
514
515 os.remove(TmpPDBOutfile)
516 else:
517 Traj[0].save(PDBOutfile)
518
519 # Write out reimaged and realinged trajectory...
520 MiscUtil.PrintInfo("Writing trajectory file %s..." % TrajOutfile)
521 Traj.save(TrajOutfile)
522
523 def ProcessOutfilePrefixParameter():
524 """Process outfile prefix paramater."""
525
526 OutfilePrefix = Options["--outfilePrefix"]
527
528 if not re.match("^auto$", OutfilePrefix, re.I):
529 OptionsInfo["OutfilePrefix"] = OutfilePrefix
530 return
531
532 if OptionsInfo["SmallMolFileMode"]:
533 OutfilePrefix = "%s_%s_Complex" % (OptionsInfo["InfileRoot"], OptionsInfo["SmallMolFileRoot"])
534 else:
535 OutfilePrefix = "%s" % (OptionsInfo["InfileRoot"])
536
537 if re.match("^yes$", Options["--waterBox"], re.I):
538 OutfilePrefix = "%s_Solvated" % (OutfilePrefix)
539
540 OutfilePrefix = "%s_%s" % (OutfilePrefix, OptionsInfo["Mode"])
541
542 OptionsInfo["OutfilePrefix"] = OutfilePrefix
543
544 def ProcessOutfileNames():
545 """Process outfile names."""
546
547 OutputParams = OptionsInfo["OutputParams"]
548 OutfileParamNames = ["CheckpointFile", "DataLogFile", "MinimizationDataLogFile", "SaveFinalStateCheckpointFile", "SaveFinalStateXMLFile", "TrajFile", "XmlSystemFile", "XmlIntegratorFile"]
549 for OutfileParamName in OutfileParamNames:
550 OutfileParamValue = OutputParams[OutfileParamName]
551 if not Options["--overwrite"]:
552 if os.path.exists(OutfileParamValue):
553 MiscUtil.PrintError("The file specified, %s, for parameter name, %s, using option \"--outfileParams\" already exist. Use option \"--ov\" or \"--overwrite\" and try again. " % (OutfileParamValue, OutfileParamName))
554
555 PDBOutfile = "%s.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
556 ReimagedPDBOutfile = "%s_Reimaged.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
557 ReimagedTrajOutfile = "%s_Reimaged.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["TrajFileExt"])
558
559 MinimizedPDBOutfile = "%s_Minimized.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
560 EquilibratedPDBOutfile = "%s_Equilibrated.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
561 FinalPDBOutfile = "%s_Final.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
562
563 for Outfile in [PDBOutfile, ReimagedPDBOutfile, ReimagedTrajOutfile, MinimizedPDBOutfile, EquilibratedPDBOutfile, FinalPDBOutfile]:
564 if not Options["--overwrite"]:
565 if os.path.exists(Outfile):
566 MiscUtil.PrintError("The file name, %s, generated using option \"--outfilePrefix\" already exist. Use option \"--ov\" or \"--overwrite\" and try again. " % (Outfile))
567 OptionsInfo["PDBOutfile"] = PDBOutfile
568 OptionsInfo["ReimagedPDBOutfile"] = ReimagedPDBOutfile
569 OptionsInfo["ReimagedTrajOutfile"] = ReimagedTrajOutfile
570
571 OptionsInfo["MinimizedPDBOutfile"] = MinimizedPDBOutfile
572 OptionsInfo["EquilibratedPDBOutfile"] = EquilibratedPDBOutfile
573 OptionsInfo["FinalPDBOutfile"] = FinalPDBOutfile
574
575 def ProcessRestartParameters():
576 """Process restart parameters. """
577
578 OptionsInfo["RestartMode"] = True if re.match("^yes$", Options["--restart"], re.I) else False
579 OptionsInfo["RestartParams"] = OpenMMUtil.ProcessOptionOpenMMRestartParameters("--restartParams", Options["--restartParams"], OptionsInfo["OutfilePrefix"])
580 if OptionsInfo["RestartMode"]:
581 RestartFile = OptionsInfo["RestartParams"]["FinalStateFile"]
582 if not os.path.exists(RestartFile):
583 MiscUtil.PrintError("The file specified, %s, for parameter name, finalStateFile, using option \"--restartParams\" doesn't exist." % (RestartFile))
584
585 DataOutAppendMode = False
586 if OptionsInfo["RestartMode"]:
587 DataOutAppendMode = True if OptionsInfo["RestartParams"]["DataAppend"] else False
588 OptionsInfo["DataOutAppendMode"] = DataOutAppendMode
589
590 def ProcessWaterBoxParameters():
591 """Process water box parameters. """
592
593 OptionsInfo["WaterBox"] = True if re.match("^yes$", Options["--waterBox"], re.I) else False
594 OptionsInfo["WaterBoxParams"] = OpenMMUtil.ProcessOptionOpenMMWaterBoxParameters("--waterBoxParams", Options["--waterBoxParams"])
595
596 if OptionsInfo["WaterBox"]:
597 if OptionsInfo["ForcefieldParams"]["ImplicitWater"]:
598 MiscUtil.PrintInfo("")
599 MiscUtil.PrintWarning("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." % (Options["--waterBox"], OptionsInfo["ForcefieldParams"]["Biopolymer"], OptionsInfo["ForcefieldParams"]["Water"]))
600
601 def ProcessOptions():
602 """Process and validate command line arguments and options."""
603
604 MiscUtil.PrintInfo("Processing options...")
605
606 ValidateOptions()
607
608 OptionsInfo["Infile"] = Options["--infile"]
609 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
610 OptionsInfo["InfileRoot"] = FileName
611
612 SmallMolFile = Options["--smallMolFile"]
613 SmallMolID = Options["--smallMolID"]
614 SmallMolFileMode = False
615 SmallMolFileRoot = None
616 if SmallMolFile is not None:
617 FileDir, FileName, FileExt = MiscUtil.ParseFileName(SmallMolFile)
618 SmallMolFileRoot = FileName
619 SmallMolFileMode = True
620
621 OptionsInfo["SmallMolFile"] = SmallMolFile
622 OptionsInfo["SmallMolFileRoot"] = SmallMolFileRoot
623 OptionsInfo["SmallMolFileMode"] = SmallMolFileMode
624 OptionsInfo["SmallMolID"] = SmallMolID.upper()
625
626 OptionsInfo["Mode"] = Options["--mode"].upper()
627 OptionsInfo["NPTMode"] = True if re.match("^NPT$", OptionsInfo["Mode"]) else False
628 OptionsInfo["NVTMode"] = True if re.match("^NVT$", OptionsInfo["Mode"]) else False
629
630 ProcessOutfilePrefixParameter()
631
632 if OptionsInfo["NVTMode"]:
633 ParamsDefaultInfoOverride = {"DataOutType": "Step Speed Progress PotentialEnergy Temperature Time Volume"}
634 else:
635 ParamsDefaultInfoOverride = {"DataOutType": "Step Speed Progress PotentialEnergy Temperature Time Density"}
636 OptionsInfo["OutputParams"] = OpenMMUtil.ProcessOptionOpenMMOutputParameters("--outputParams", Options["--outputParams"], OptionsInfo["OutfilePrefix"], ParamsDefaultInfoOverride)
637 ProcessOutfileNames()
638
639 OptionsInfo["ForcefieldParams"] = OpenMMUtil.ProcessOptionOpenMMForcefieldParameters("--forcefieldParams", Options["--forcefieldParams"])
640
641 OptionsInfo["FreezeAtoms"] = True if re.match("^yes$", Options["--freezeAtoms"], re.I) else False
642 if OptionsInfo["FreezeAtoms"]:
643 OptionsInfo["FreezeAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters("--freezeAtomsParams", Options["--freezeAtomsParams"])
644 else:
645 OptionsInfo["FreezeAtomsParams"] = None
646
647 ParamsDefaultInfoOverride = {"Name": Options["--platform"], "Threads": 1}
648 OptionsInfo["PlatformParams"] = OpenMMUtil.ProcessOptionOpenMMPlatformParameters("--platformParams", Options["--platformParams"], ParamsDefaultInfoOverride)
649
650 OptionsInfo["RestraintAtoms"] = True if re.match("^yes$", Options["--restraintAtoms"], re.I) else False
651 if OptionsInfo["RestraintAtoms"]:
652 OptionsInfo["RestraintAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters("--restraintAtomsParams", Options["--restraintAtomsParams"])
653 else:
654 OptionsInfo["RestraintAtomsParams"] = None
655 OptionsInfo["RestraintSpringConstant"] = float(Options["--restraintSpringConstant"])
656
657 ProcessRestartParameters()
658
659 OptionsInfo["SystemParams"] = OpenMMUtil.ProcessOptionOpenMMSystemParameters("--systemParams", Options["--systemParams"])
660
661 OptionsInfo["IntegratorParams"] = OpenMMUtil.ProcessOptionOpenMMIntegratorParameters("--integratorParams", Options["--integratorParams"], HydrogenMassRepartioningStatus = OptionsInfo["SystemParams"]["HydrogenMassRepartioning"])
662
663 OptionsInfo["SimulationParams"] = OpenMMUtil.ProcessOptionOpenMMSimulationParameters("--simulationParams", Options["--simulationParams"])
664
665 ProcessWaterBoxParameters()
666
667 OptionsInfo["Overwrite"] = Options["--overwrite"]
668
669 def RetrieveOptions():
670 """Retrieve command line arguments and options."""
671
672 # Get options...
673 global Options
674 Options = docopt(_docoptUsage_)
675
676 # Set current working directory to the specified directory...
677 WorkingDir = Options["--workingdir"]
678 if WorkingDir:
679 os.chdir(WorkingDir)
680
681 # Handle examples option...
682 if "--examples" in Options and Options["--examples"]:
683 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
684 sys.exit(0)
685
686 def ValidateOptions():
687 """Validate option values."""
688
689 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
690 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
691
692 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--infile"])
693 OutfilePrefix = Options["--outfilePrefix"]
694 if not re.match("^auto$", OutfilePrefix, re.I):
695 if re.match("^(%s)$" % OutfilePrefix, FileName, re.I):
696 MiscUtil.PrintError("The value specified, %s, for option \"--outfilePrefix\" is not valid. You must specify a value different from, %s, the root of infile name." % (OutfilePrefix, FileName))
697
698 if Options["--smallMolFile"] is not None:
699 MiscUtil.ValidateOptionFilePath("-l, --smallMolFile", Options["--smallMolFile"])
700 MiscUtil.ValidateOptionFileExt("-l, --smallMolFile", Options["--smallMolFile"], "sd sdf")
701
702 SmallMolID = Options["--smallMolID"]
703 if len(SmallMolID) != 3:
704 MiscUtil.PrintError("The value specified, %s, for option \"--smallMolID\" is not valid. You must specify a three letter small molecule ID." % (SmallMolID))
705
706 MiscUtil.ValidateOptionTextValue("--freezeAtoms", Options["--freezeAtoms"], "yes no")
707 if re.match("^yes$", Options["--freezeAtoms"], re.I):
708 if Options["--freezeAtomsParams"] is None:
709 MiscUtil.PrintError("No value specified for option \"--freezeAtomsParams\". You must specify valid values during, yes, value for \"--freezeAtoms\" option.")
710
711 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "NPT NVT")
712
713 MiscUtil.ValidateOptionTextValue("-p, --platform", Options["--platform"], "CPU CUDA OpenCL Reference")
714
715 MiscUtil.ValidateOptionTextValue("-r, --restart ", Options["--restart"], "yes no")
716
717 MiscUtil.ValidateOptionTextValue("--restraintAtoms", Options["--restraintAtoms"], "yes no")
718 if re.match("^yes$", Options["--restraintAtoms"], re.I):
719 if Options["--restraintAtomsParams"] is None:
720 MiscUtil.PrintError("No value specified for option \"--restraintAtomsParams\". You must specify valid values during, yes, value for \"--restraintAtoms\" option.")
721
722 MiscUtil.ValidateOptionFloatValue("--restraintSpringConstant", Options["--restraintSpringConstant"], {">": 0})
723
724 MiscUtil.ValidateOptionTextValue("--waterBox", Options["--waterBox"], "yes no")
725
726 # Setup a usage string for docopt...
727 _docoptUsage_ = """
728 OpenMMPerformMDSimulation.py - Perform a MD simulation
729
730 Usage:
731 OpenMMPerformMDSimulation.py [--forcefieldParams <Name,Value,..>] [--freezeAtoms <yes or no>]
732 [--freezeAtomsParams <Name,Value,..>] [--integratorParams <Name,Value,..>]
733 [--mode <NVT or NPT>] [--outputParams <Name,Value,..>] [--outfilePrefix <text>]
734 [--overwrite] [--platform <text>] [--platformParams <Name,Value,..>] [--restart <yes or no>]
735 [--restartParams <Name,Value,..>] [--restraintAtoms <yes or no>]
736 [--restraintAtomsParams <Name,Value,..>] [--restraintSpringConstant <number>]
737 [--simulationParams <Name,Value,..>] [--smallMolFile <SmallMolFile>] [--smallMolID <text>]
738 [--systemParams <Name,Value,..>] [--waterBox <yes or no>]
739 [--waterBoxParams <Name,Value,..>] [-w <dir>] -i <infile>
740 OpenMMPerformMDSimulation.py -h | --help | -e | --examples
741
742 Description:
743 Perform a MD simulation using an NPT or NVT statistical ensemble. You may
744 run a simulation using a macromolecule or a macromolecule in a complex with
745 small molecule. By default, the system is minimized and equilibrated before the
746 production run.
747
748 The input file must contain a macromolecule already prepared for simulation.
749 The preparation of the macromolecule for a simulation generally involves the
750 following: tasks: Identification and replacement non-standard residues;
751 Addition of missing residues; Addition of missing heavy atoms; Addition of
752 missing hydrogens; Addition of a water box which is optional.
753
754 In addition, the small molecule input file must contain a molecule already
755 prepared for simulation. It must contain appropriate 3D coordinates relative
756 to the macromolecule along with no missing hydrogens.
757
758 You may optionally add a water box and freeze/restraint atoms for the
759 simulation.
760
761 The supported macromolecule input file formats are: PDB (.pdb) and
762 CIF (.cif)
763
764 The supported small molecule input file format are : SD (.sdf, .sd)
765
766 Possible outfile prefixes:
767
768 <InfieRoot>_<Mode>
769 <InfieRoot>_Solvated_<Mode>
770 <InfieRoot>_<SmallMolFileRoot>_Complex_<Mode>,
771 <InfieRoot>_<SmallMolFileRoot>_Complex_Solvated_<Mode>
772
773 Possible output files:
774
775 <OutfilePrefix>.<pdb or cif> [ Initial sytem ]
776 <OutfilePrefix>.<dcd or xtc>
777
778 <OutfilePrefix>_Reimaged.<pdb or cif> [ First frame ]
779 <OutfilePrefix>_Reimaged.<dcd or xtc>
780
781 <OutfilePrefix>_Minimized.<pdb or cif>
782 <OutfilePrefix>_Equilibrated.<pdb or cif>
783 <OutfilePrefix>_Final.<pdb or cif> [ Final system ]
784
785 <OutfilePrefix>.chk
786 <OutfilePrefix>.csv
787 <OutfilePrefix>_Minimization.csv
788 <OutfilePrefix>_FinalState.chk
789 <OutfilePrefix>_FinalState.xml
790
791 <OutfilePrefix>_System.xml
792 <OutfilePrefix>_Integrator.xml
793
794 The reimaged PDB file, <OutfilePrefix>_Reimaged.pdb, corresponds to the first
795 frame in the trajectory. The reimaged trajectory file contains all the frames
796 aligned to the first frame after reimaging of the frames for periodic systems.
797
798 Options:
799 -e, --examples
800 Print examples.
801 -f, --forcefieldParams <Name,Value,..> [default: auto]
802 A comma delimited list of parameter name and value pairs for biopolymer,
803 water, and small molecule forcefields.
804
805 The supported parameter names along with their default values are
806 are shown below:
807
808 biopolymer, amber14-all.xml [ Possible values: Any Valid value ]
809 smallMolecule, openff-2.2.1 [ Possible values: Any Valid value ]
810 water, auto [ Possible values: Any Valid value ]
811 additional, none [ Possible values: Space delimited list of any
812 valid value ]
813
814 Possible biopolymer forcefield values:
815
816 amber14-all.xml, amber99sb.xml, amber99sbildn.xml, amber03.xml,
817 amber10.xml
818 charmm36.xml, charmm_polar_2019.xml
819 amoeba2018.xml
820
821 Possible small molecule forcefield values:
822
823 openff-2.2.1, openff-2.0.0, openff-1.3.1, openff-1.2.1,
824 openff-1.1.1, openff-1.1.0,...
825 smirnoff99Frosst-1.1.0, smirnoff99Frosst-1.0.9,...
826 gaff-2.11, gaff-2.1, gaff-1.81, gaff-1.8, gaff-1.4,...
827
828 The default water forcefield valus is dependent on the type of the
829 biopolymer forcefield as shown below:
830
831 Amber: amber14/tip3pfb.xml
832 CHARMM: charmm36/water.xml or None for charmm_polar_2019.xml
833 Amoeba: None (Explicit)
834
835 Possible water forcefield values:
836
837 amber14/tip3p.xml, amber14/tip3pfb.xml, amber14/spce.xml,
838 amber14/tip4pew.xml, amber14/tip4pfb.xml,
839 charmm36/water.xml, charmm36/tip3p-pme-b.xml,
840 charmm36/tip3p-pme-f.xml, charmm36/spce.xml,
841 charmm36/tip4pew.xml, charmm36/tip4p2005.xml,
842 charmm36/tip5p.xml, charmm36/tip5pew.xml,
843 implicit/obc2.xml, implicit/GBn.xml, implicit/GBn2.xml,
844 amoeba2018_gk.xml (Implict water)
845 None (Explicit water for amoeba)
846
847 The additional forcefield value is a space delimited list of any valid
848 forcefield values and is passed on to the OpenMMForcefields
849 SystemGenerator along with the specified forcefield values for
850 biopolymer, water, and mall molecule. Possible additional forcefield
851 values are:
852
853 amber14/DNA.OL15.xml amber14/RNA.OL3.xml
854 amber14/lipid17.xml amber14/GLYCAM_06j-1.xml
855 ... ... ...
856
857 You may specify any valid forcefield names supported by OpenMM. No
858 explicit validation is performed.
859 --freezeAtoms <yes or no> [default: no]
860 Freeze atoms during a simulation. The specified atoms are kept completely
861 fixed by setting their masses to zero. Their positions do not change during
862 local energy minimization and MD simulation, and they do not contribute
863 to the kinetic energy of the system.
864 --freezeAtomsParams <Name,Value,..>
865 A comma delimited list of parameter name and value pairs for freezing
866 atoms during a simulation. You must specify these parameters for 'yes'
867 value of '--freezeAtoms' option.
868
869 The supported parameter names along with their default values are
870 are shown below:
871
872 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
873 Protein, Residues, or Water ]
874 selectionSpec, auto [ Possible values: A space delimited list of
875 residue names ]
876 negate, no [ Possible values: yes or no ]
877
878 A brief description of parameters is provided below:
879
880 selection: Atom selection to freeze.
881 selectionSpec: A space delimited list of residue names for
882 selecting atoms to freeze. You must specify its value during
883 'Ligand' and 'Protein' value for 'selection'. The default values
884 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
885 and 'Water' values of 'selection' as shown below:
886
887 CAlphaProtein: List of stadard protein residues from pdbfixer
888 for selecting CAlpha atoms.
889 Ions: Li Na K Rb Cs Cl Br F I
890 Water: HOH
891 Protein: List of standard protein residues from pdbfixer.
892
893 negate: Negate atom selection match to select atoms for freezing.
894
895 In addition, you may specify an explicit space delimited list of residue
896 names using 'selectionSpec' for any 'selection". The specified residue
897 names are appended to the appropriate default values during the
898 selection of atoms for freezing.
899 -h, --help
900 Print this help message.
901 -i, --infile <infile>
902 Input file name containing a macromolecule.
903 --integratorParams <Name,Value,..> [default: auto]
904 A comma delimited list of parameter name and value pairs for integrator
905 during a simulation.
906
907 The supported parameter names along with their default values are
908 are shown below:
909
910 integrator, LangevinMiddle [ Possible values: LangevinMiddle,
911 Langevin, NoseHoover, Brownian ]
912
913 randomSeed, auto [ Possible values: > 0 ]
914
915 frictionCoefficient, 1.0 [ Units: 1/ps ]
916 stepSize, auto [ Units: fs; Default value: 4 fs during yes value of
917 hydrogen mass repartioning with no freezing/restraining of atoms;
918 otherwsie, 2 fs ]
919 temperature, 300.0 [ Units: kelvin ]
920
921 barostat, MonteCarlo [ Possible values: MonteCarlo or
922 MonteCarloMembrane ]
923 barostatInterval, 25
924 pressure, 1.0 [ Units: atm ]
925
926 Parameters used only for MonteCarloMembraneBarostat with default
927 values corresponding to Amber forcefields:
928
929 surfaceTension, 0.0 [ Units: atm*A. It is automatically converted
930 into OpenMM default units of atm*nm before its usage. ]
931 xymode, Isotropic [ Possible values: Anisotropic or Isotropic ]
932 zmode, Free [ Possible values: Free or Fixed ]
933
934 A brief description of parameters is provided below:
935
936 integrator: Type of integrator
937
938 randomSeed: Random number seed for barostat and integrator. Not
939 supported for NoseHoover integrator.
940
941 frictionCoefficient: Friction coefficient for coupling the system to
942 the heat bath..
943 stepSize: Simulation time step size.
944 temperature: Simulation temperature.
945
946 barostat: Barostat type.
947 barostatInterval: Barostat interval step size during NPT
948 simulation for applying Monte Carlo pressure changes.
949 pressure: Pressure during NPT simulation.
950
951 Parameters used only for MonteCarloMembraneBarostat:
952
953 surfaceTension: Surface tension acting on the system.
954 xymode: Behavior along X and Y axes. You may allow the X and Y axes
955 to vary independently of each other or always scale them by the same
956 amount to keep the ratio of their lengths constant.
957 zmode: Beahvior along Z axis. You may allow the Z axis to vary
958 independently of the other axes or keep it fixed.
959
960 -m, --mode <NPT or NVT> [default: NPT]
961 Type of statistical ensemble to use for simulation. Possible values:
962 NPT (constant Number of particles, Pressure, and Temperature) or
963 NVT ((constant Number of particles, Volume and Temperature)
964 --outfilePrefix <text> [default: auto]
965 File prefix for generating the names of output files. The default value
966 depends on the names of input files for macromolecule and small molecule
967 along with the type of statistical ensemble and the nature of the solvation.
968
969 The possible values for outfile prefix are shown below:
970
971 <InfieRoot>_<Mode>
972 <InfieRoot>_Solvated_<Mode>
973 <InfieRoot>_<SmallMolFileRoot>_Complex_<Mode>,
974 <InfieRoot>_<SmallMolFileRoot>_Complex_Solvated_<Mode>
975
976 --outputParams <Name,Value,..> [default: auto]
977 A comma delimited list of parameter name and value pairs for generating
978 output during a simulation.
979
980 The supported parameter names along with their default values are
981 are shown below:
982
983 checkpoint, no [ Possible values: yes or no ]
984 checkpointFile, auto [ Default: <OutfilePrefix>.chk ]
985 checkpointSteps, 10000
986
987 dataOutType, auto [ Possible values: A space delimited list of valid
988 parameter names.
989 NPT simulation default: Density Step Speed Progress
990 PotentialEnergy Temperature Time.
991 NVT simulation default: Step Speed Progress PotentialEnergy
992 Temperature Time Volume
993 Other valid names: ElapsedTime RemainingTime KineticEnergy
994 TotalEnergy ]
995
996 dataLog, yes [ Possible values: yes or no ]
997 dataLogFile, auto [ Default: <OutfilePrefix>.csv ]
998 dataLogSteps, 1000
999
1000 dataStdout, no [ Possible values: yes or no ]
1001 dataStdoutSteps, 1000
1002
1003 minimizationDataSteps, 100
1004 minimizationDataStdout, no [ Possible values: yes or no ]
1005 minimizationDataLog, no [ Possible values: yes or no ]
1006 minimizationDataLogFile, auto [ Default:
1007 <OutfilePrefix>_MinimizationOut.csv ]
1008 minimizationDataOutType, auto [ Possible values: A space delimited
1009 list of valid parameter names. Default: SystemEnergy
1010 RestraintEnergy MaxConstraintError.
1011 Other valid names: RestraintStrength ]
1012
1013 pdbOutFormat, PDB [ Possible values: PDB or CIF ]
1014 pdbOutKeepIDs, yes [ Possible values: yes or no ]
1015
1016 pdbOutMinimized, no [ Possible values: yes or no ]
1017 pdbOutEquilibrated, no [ Possible values: yes or no ]
1018 pdbOutFinal, no [ Possible values: yes or no ]
1019
1020 saveFinalStateCheckpoint, yes [ Possible values: yes or no ]
1021 saveFinalStateCheckpointFile, auto [ Default:
1022 <OutfilePrefix>_FinalState.chk ]
1023 saveFinalStateXML, no [ Possible values: yes or no ]
1024 saveFinalStateXMLFile, auto [ Default:
1025 <OutfilePrefix>_FinalState.xml]
1026
1027 traj, yes [ Possible values: yes or no ]
1028 trajFile, auto [ Default: <OutfilePrefix>.<TrajFormat> ]
1029 trajFormat, DCD [ Possible values: DCD or XTC ]
1030 trajSteps, 10000 [ The default value corresponds to 40 ps for step
1031 size of 4 fs. ]
1032
1033 xmlSystemOut, no [ Possible values: yes or no ]
1034 xmlSystemFile, auto [ Default: <OutfilePrefix>_System.xml ]
1035 xmlIntegratorOut, no [ Possible values: yes or no ]
1036 xmlIntegratorFile, auto [ Default: <OutfilePrefix>_Integrator.xml ]
1037
1038 A brief description of parameters is provided below:
1039
1040 checkpoint: Write intermediate checkpoint file.
1041 checkpointFile: Intermediate checkpoint file name.
1042 checkpointSteps: Frequency of writing intermediate checkpoint file.
1043
1044 dataOutType: Type of data to write to stdout and log file.
1045
1046 dataLog: Write data to log file.
1047 dataLogFile: Data log file name.
1048 dataLogSteps: Frequency of writing data to log file.
1049
1050 dataStdout: Write data to stdout.
1051 dataStdoutSteps: Frequency of writing data to stdout.
1052
1053 minimizationDataSteps: Frequency of writing data to stdout
1054 and log file.
1055 minimizationDataStdout: Write data to stdout.
1056 minimizationDataLog: Write data to log file.
1057 minimizationDataLogFile: Data log fie name.
1058 minimizationDataOutType: Type of data to write to stdout
1059 and log file.
1060
1061 saveFinalStateCheckpoint: Save final state checkpoint file.
1062 saveFinalStateCheckpointFile: Name of final state checkpoint file.
1063 saveFinalStateXML: Save final state XML file.
1064 saveFinalStateXMLFile: Name of final state XML file.
1065
1066 pdbOutFormat: Format of output PDB files.
1067 pdbOutKeepIDs: Keep existing chain and residue IDs.
1068
1069 pdbOutMinimized: Write PDB file after minimization.
1070 pdbOutEquilibrated: Write PDB file after equilibration.
1071 pdbOutFinal: Write final PDB file after production run.
1072
1073 traj: Write out trajectory file.
1074 trajFile: Trajectory file name.
1075 trajFormat: Trajectory file format.
1076 trajSteps: Frequency of writing trajectory file.
1077
1078 xmlSystemOut: Write system XML file.
1079 xmlSystemFile: System XML file name.
1080 xmlIntegratorOut: Write integrator XML file.
1081 xmlIntegratorFile: Integrator XML file name.
1082
1083 --overwrite
1084 Overwrite existing files.
1085 -p, --platform <text> [default: CPU]
1086 Platform to use for running MD simulation. Possible values: CPU, CUDA,
1087 OpenCL, or Reference.
1088 --platformParams <Name,Value,..> [default: auto]
1089 A comma delimited list of parameter name and value pairs to configure
1090 platform for running MD simulation.
1091
1092 The supported parameter names along with their default values for
1093 different platforms are shown below:
1094
1095 CPU:
1096
1097 threads, 1 [ Possible value: >= 0 or auto. The value of 'auto'
1098 or zero implies the use of all available CPUs for threading. ]
1099
1100 CUDA:
1101
1102 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
1103 deterministicForces, auto [ Possible values: yes or no ]
1104 precision, single [ Possible values: single, double, or mix ]
1105 tempDirectory, auto [ Possible value: DirName ]
1106 useBlockingSync, auto [ Possible values: yes or no ]
1107 useCpuPme, auto [ Possible values: yes or no ]
1108
1109 OpenCL:
1110
1111 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
1112 openCLPlatformIndex, auto [ Possible value: Number]
1113 precision, single [ Possible values: single, double, or mix ]
1114 useCpuPme, auto [ Possible values: yes or no ]
1115
1116 A brief description of parameters is provided below:
1117
1118 CPU:
1119
1120 threads: Number of threads to use for simulation.
1121
1122 CUDA:
1123
1124 deviceIndex: Space delimited list of device indices to use for
1125 calculations.
1126 deterministicForces: Generate reproducible results at the cost of a
1127 small decrease in performance.
1128 precision: Number precision to use for calculations.
1129 tempDirectory: Directory name for storing temporary files.
1130 useBlockingSync: Control run-time synchronization between CPU and
1131 GPU.
1132 useCpuPme: Use CPU-based PME implementation.
1133
1134 OpenCL:
1135
1136 deviceIndex: Space delimited list of device indices to use for
1137 simulation.
1138 openCLPlatformIndex: Platform index to use for calculations.
1139 precision: Number precision to use for calculations.
1140 useCpuPme: Use CPU-based PME implementation.
1141
1142 -r, --restart <yes or no> [default: no]
1143 Restart simulation using a previously saved final state checkpoint or
1144 XML file.
1145 --restartParams <Name,Value,..> [default: auto]
1146 A comma delimited list of parameter name and value pairs for restarting
1147 a simulation.
1148
1149 The supported parameter names along with their default values are
1150 are shown below:
1151
1152 finalStateFile, <OutfilePrefix>_FinalState.<chk> [ Possible values:
1153 Valid final state checkpoint or XML filename ]
1154 dataAppend, yes [ Possible values: yes or no]
1155
1156 A brief description of parameters is provided below:
1157
1158 finalStateFile: Final state checkpoint or XML file
1159 dataAppend: Append data to existing trajectory and data log files
1160 during the restart of a simulation using a previously saved final
1161 state checkpoint or XML file.
1162
1163 --restraintAtoms <yes or no> [default: no]
1164 Restraint atoms during a simulation. The motion of specified atoms is
1165 restricted by adding a harmonic force that binds them to their starting
1166 positions. The atoms are not completely fixed unlike freezing of atoms.
1167 Their motion, however, is restricted and they are not able to move far away
1168 from their starting positions during local energy minimization and MD
1169 simulation.
1170 --restraintAtomsParams <Name,Value,..>
1171 A comma delimited list of parameter name and value pairs for restraining
1172 atoms during a simulation. You must specify these parameters for 'yes'
1173 value of '--restraintAtoms' option.
1174
1175 The supported parameter names along with their default values are
1176 are shown below:
1177
1178 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
1179 Protein, Residues, or Water ]
1180 selectionSpec, auto [ Possible values: A space delimited list of
1181 residue names ]
1182 negate, no [ Possible values: yes or no ]
1183
1184 A brief description of parameters is provided below:
1185
1186 selection: Atom selection to restraint.
1187 selectionSpec: A space delimited list of residue names for
1188 selecting atoms to restraint. You must specify its value during
1189 'Ligand' and 'Protein' value for 'selection'. The default values
1190 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
1191 and 'Water' values of 'selection' as shown below:
1192
1193 CAlphaProtein: List of stadard protein residues from pdbfixer
1194 for selecting CAlpha atoms.
1195 Ions: Li Na K Rb Cs Cl Br F I
1196 Water: HOH
1197 Protein: List of standard protein residues from pdbfixer.
1198
1199 negate: Negate atom selection match to select atoms for freezing.
1200
1201 In addition, you may specify an explicit space delimited list of residue
1202 names using 'selectionSpec' for any 'selection". The specified residue
1203 names are appended to the appropriate default values during the
1204 selection of atoms for restraining.
1205 --restraintSpringConstant <number> [default: 2.5]
1206 Restraint spring constant for applying external restraint force to restraint
1207 atoms relative to their initial positions during 'yes' value of '--restraintAtoms'
1208 option. Default units: kcal/mol/A**2. The default value, 2.5, corresponds to
1209 1046.0 kjoules/mol/nm**2. The default value is automatically converted into
1210 units of kjoules/mol/nm**2 before its usage.
1211 --simulationParams <Name,Value,..> [default: auto]
1212 A comma delimited list of parameter name and value pairs for simulation.
1213
1214 The supported parameter names along with their default values are
1215 are shown below:
1216
1217 steps, 1000000 [ Possible values: > 0. The default value
1218 corresponds to 4 ns for step size of 4 fs. ]
1219
1220 minimization, yes [ Possible values: yes or no ]
1221 minimizationMaxSteps, auto [ Possible values: >= 0. The value of
1222 zero implies until the minimization is converged. ]
1223 minimizationTolerance, 0.24 [ Units: kcal/mol/A. The default value
1224 0.25, corresponds to OpenMM default of value of 10.04
1225 kjoules/mol/nm. It is automatically converted into OpenMM
1226 default units before its usage. ]
1227
1228 equilibration, yes [ Possible values: yes or no ]
1229 equilibrationSteps, 1000 [ Possible values: > 0 ]
1230
1231 A brief description of parameters is provided below:
1232
1233 steps: Number of steps for production run.
1234
1235 minimization: Perform minimization before equilibration and
1236 production run.
1237 minimizationMaxSteps: Maximum number of minimization steps. The
1238 value of zero implies until the minimization is converged.
1239 minimizationTolerance: Energy convergence tolerance during
1240 minimization.
1241
1242 equilibration: Perform equilibration before the production run.
1243 equilibrationSteps: Number of steps for equilibration.
1244
1245 -s, --smallMolFile <SmallMolFile>
1246 Small molecule input file name. The macromolecue and small molecule are
1247 merged for simulation and the complex is written out to a PDB file.
1248 --smallMolID <text> [default: LIG]
1249 Three letter small molecule residue ID. The small molecule ID corresponds
1250 to the residue name of the small molecule and is written out to a PDB file
1251 containing the complex.
1252 --systemParams <Name,Value,..> [default: auto]
1253 A comma delimited list of parameter name and value pairs to configure
1254 a system for simulation.
1255
1256 The supported parameter names along with their default values are
1257 are shown below:
1258
1259 constraints, BondsInvolvingHydrogens [ Possible values: None,
1260 WaterOnly, BondsInvolvingHydrogens, AllBonds, or
1261 AnglesInvolvingHydrogens ]
1262 constraintErrorTolerance, 0.000001
1263 ewaldErrorTolerance, 0.0005
1264
1265 nonbondedMethodPeriodic, PME [ Possible values: NoCutoff,
1266 CutoffNonPeriodic, or PME ]
1267 nonbondedMethodNonPeriodic, NoCutoff [ Possible values:
1268 NoCutoff or CutoffNonPeriodic]
1269 nonbondedCutoff, 1.0 [ Units: nm ]
1270
1271 hydrogenMassRepartioning, yes [ Possible values: yes or no ]
1272 hydrogenMass, 1.5 [ Units: amu]
1273
1274 removeCMMotion, yes [ Possible values: yes or no ]
1275 rigidWater, auto [ Possible values: yes or no. Default: 'No' for
1276 'None' value of constraints; Otherwise, yes ]
1277
1278 A brief description of parameters is provided below:
1279
1280 constraints: Type of system constraints to use for simulation. These
1281 constraints are different from freezing and restraining of any
1282 atoms in the system.
1283 constraintErrorTolerance: Distance tolerance for constraints as a
1284 fraction of the constrained distance.
1285 ewaldErrorTolerance: Ewald error tolerance for a periodic system.
1286
1287 nonbondedMethodPeriodic: Nonbonded method to use during the
1288 calculation of long range interactions for a periodic system.
1289 nonbondedMethodNonPeriodic: Nonbonded method to use during the
1290 calculation of long range interactions for a non-periodic system.
1291 nonbondedCutoff: Cutoff distance to use for long range interactions
1292 in both perioidic non-periodic systems.
1293
1294 hydrogenMassRepartioning: Use hydrogen mass repartioning. It
1295 increases the mass of the hydrogen atoms attached to the heavy
1296 atoms and decreasing the mass of the bonded heavy atom to
1297 maintain constant system mass. This allows the use of larger
1298 integration step size (4 fs) during a simulation.
1299 hydrogenMass: Hydrogen mass to use during repartioning.
1300
1301 removeCMMotion: Remove all center of mass motion at every time step.
1302 rigidWater: Keep water rigid during a simulation. This is determined
1303 automatically based on the value of 'constraints' parameter.
1304
1305 --waterBox <yes or no> [default: no]
1306 Add water box.
1307 --waterBoxParams <Name,Value,..> [default: auto]
1308 A comma delimited list of parameter name and value pairs for adding
1309 a water box.
1310
1311 The supported parameter names along with their default values are
1312 are shown below:
1313
1314 model, tip3p [ Possible values: tip3p, spce, tip4pew, tip5p or
1315 swm4ndp ]
1316 mode, Padding [ Possible values: Size or Padding ]
1317 padding, 1.0
1318 size, None [ Possible value: xsize ysize zsize ]
1319 shape, cube [ Possible values: cube, dodecahedron, or octahedron ]
1320 ionPositive, Na+ [ Possible values: Li+, Na+, K+, Rb+, or Cs+ ]
1321 ionNegative, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
1322 ionicStrength, 0.0
1323
1324 A brief description of parameters is provided below:
1325
1326 model: Water model to use for adding water box. The van der
1327 Waals radii and atomic charges are determined using the
1328 specified water forcefield. You must specify an appropriate
1329 water forcefield. No validation is performed.
1330 mode: Specify the size of the waterbox explicitly or calculate it
1331 automatically for a macromolecule along with adding padding
1332 around ther macromolecule.
1333 padding: Padding around a macromolecule in nanometers for filling
1334 box with water. It must be specified during 'Padding' value of
1335 'mode' parameter.
1336 size: A space delimited triplet of values corresponding to water
1337 size in nanometers. It must be specified during 'Size' value of
1338 'mode' parameter.
1339 ionPositive: Type of positive ion to add during the addition of a
1340 water box.
1341 ionNegative: Type of negative ion to add during the addition of a
1342 water box.
1343 ionicStrength: Total concentration of both positive and negative
1344 ions to add excluding the ions added to neutralize the system
1345 during the addition of a water box.
1346
1347 -w, --workingdir <dir>
1348 Location of working directory which defaults to the current directory.
1349
1350 Examples:
1351 To perform a MD simulation for a macromolecule in a PDB file by using an NPT
1352 ensemble, applying system constraints for bonds involving hydrogens along
1353 with hydrogen mass repartioning, using a step size of 4 fs, performing minimization
1354 until it's converged along with equilibration for 1,000 steps ( 4 ps), performing
1355 production run for 1,000,000 steps (4 ns), writing trajectory file every 10,000
1356 steps (40 ps), writing data log file every 1,000 steps (4 ps), generating a checkpoint
1357 file after the completion of the, and generating a PDB for the final system, type:
1358
1359 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes
1360
1361 To run the first example for performing OpenMM simulation using multi-
1362 threading employing all available CPUs on your machine and generate various
1363 output files, type:
1364
1365 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes
1366 --platformParams "threads,0"
1367
1368 To run the first example for performing OpenMM simulation using CUDA platform
1369 on your machine and generate various output files, type:
1370
1371 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes -p CUDA
1372
1373 To run the second example for performing NVT simulation minimizing for a
1374 maximum of 2,000 steps, performing production run of 10,00 steps (40 ps),
1375 writing trajectory file every 1,000 steps (4 ps), and generate various output
1376 files, type:
1377
1378 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes
1379 --platformParams "threads,0" --simulationParams "steps,10000,
1380 minimizationMaxSteps, 1000" --outputParams "trajSteps,1000"
1381
1382 To run the second example for a marcomolecule in a complex with a small
1383 molecule and generate various output files, type:
1384
1385 % OpenMMPerformMDSimulation.py -i Sample13.pdb -s Sample13Ligand.sdf
1386 --waterBox yes --platformParams "threads,0"
1387
1388 To run the second example for performing NVT simulation and generate various
1389 output files, type:
1390
1391 % OpenMMPerformMDSimulation.py -i Sample13.pdb -s Sample13Ligand.sdf
1392 --mode NVT --platformParams "threads,0"
1393
1394 % OpenMMPerformMDSimulation.py -i Sample13.pdb -s Sample13Ligand.sdf
1395 --mode NVT --waterBox yes --platformParams "threads,0"
1396
1397 To run the second example for a macromolecule in a lipid bilayer, write out
1398 reimaged and realigned trajectory file for the periodic system, along with a
1399 PDB file corresponding to the first frame, and generate various output files,
1400 type:
1401
1402 % OpenMMPerformMDSimulation.py -i Sample12LipidBilayer.pdb
1403 --platformParams "threads,0" --integratorParams
1404 "barostat,MonteCarloMembrane"
1405
1406 To run the second example by freezing CAlpha atoms in a macromolecule without
1407 using any system constraints to avoid any issues with the freezing of the same atoms,
1408 using a step size of 2 fs, and generate various output files, type:
1409
1410 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes
1411 --freezeAtoms yes --freezeAtomsParams "selection,CAlphaProtein"
1412 --systemParams "constraints, None"
1413 --platformParams "threads,0" --integratorParams "stepSize,2"
1414
1415 To run the second example by restrainting CAlpha atoms in a macromolecule and
1416 and generate various output files, type:
1417
1418 % OpenMMPerformMDSimulation.py -i Sample13.pdb --waterBox yes
1419 --restraintAtomsParams "selection,CAlphaProtein"
1420 --platformParams "threads,0" --integratorParams "stepSize,2"
1421
1422 To run the second example for a marcomolecule in a complex with a small
1423 molecule by using implicit water and generate various output files, type:
1424
1425 % OpenMMPerformMDSimulation.py -i Sample13.pdb -s Sample13Ligand.sdf
1426 --mode NVT --platformParams "threads,0"
1427 --forcefieldParams "biopolymer, amber14-all.xml, water,
1428 implicit/obc2.xml"
1429
1430 To run the second example by specifying explict values for various parametres
1431 and generate various output files, type:
1432
1433 % OpenMMPerformMDSimulation.py -m NPT -i Sample13.pdb
1434 -f " biopolymer,amber14-all.xml,smallMolecule, openff-2.2.1,
1435 water,amber14/tip3pfb.xml" --waterBox yes
1436 --integratorParams "integrator,LangevinMiddle,randomSeed,42,
1437 stepSize,2, temperature, 300.0,pressure, 1.0"
1438 --outputParams "checkpoint,yes,dataLog,yes,dataStdout,yes,
1439 minimizationDataStdout,yes,minimizationDataLog,yes,
1440 pdbOutFormat,CIF,pdbOutKeepIDs,yes,saveFinalStateCheckpoint, yes,
1441 traj,yes,xmlSystemOut,yes,xmlIntegratorOut,yes"
1442 -p CPU --platformParams "threads,0"
1443 --simulationParams "steps,10000,minimization,yes,
1444 minimizationMaxSteps,500,equilibration,yes,equilibrationSteps,1000"
1445 --systemParams "constraints,BondsInvolvingHydrogens,
1446 nonbondedMethodPeriodic,PME,nonbondedMethodNonPeriodic,NoCutoff,
1447 hydrogenMassRepartioning, yes"
1448
1449 Author:
1450 Manish Sud(msud@san.rr.com)
1451
1452 See also:
1453 OpenMMPrepareMacromolecule.py, OpenMMPerformMinimization.py,
1454 OpenMMPerformSimulatedAnnealing.py
1455
1456 Copyright:
1457 Copyright (C) 2025 Manish Sud. All rights reserved.
1458
1459 The functionality available in this script is implemented using OpenMM, an
1460 open source molecuar simulation package.
1461
1462 This file is part of MayaChemTools.
1463
1464 MayaChemTools is free software; you can redistribute it and/or modify it under
1465 the terms of the GNU Lesser General Public License as published by the Free
1466 Software Foundation; either version 3 of the License, or (at your option) any
1467 later version.
1468
1469 """
1470
1471 if __name__ == "__main__":
1472 main()
