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 PerformMinimization()
83
84 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
85 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
86
87 def PerformMinimization():
88 """Perform minimization."""
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 # Perform energy minimization...
103 PerformEnergyMinimization(Simulation)
104
105 def PrepareSystem():
106 """Prepare system for simulation."""
107
108 System, Topology, Positions = OpenMMUtil.InitializeSystem(OptionsInfo["Infile"], OptionsInfo["ForcefieldParams"], OptionsInfo["SystemParams"], OptionsInfo["WaterBox"], OptionsInfo["WaterBoxParams"], OptionsInfo["SmallMolFile"], OptionsInfo["SmallMolID"])
109
110 # Write out a PDB file for the system...
111 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["InitialPDBOutfile"])
112 OpenMMUtil.WritePDBFile(OptionsInfo["InitialPDBOutfile"], Topology, Positions, OptionsInfo["OutputParams"]["PDBOutKeepIDs"])
113
114 return (System, Topology, Positions)
115
116 def SetupIntegrator():
117 """Setup integrator. """
118
119 Integrator = OpenMMUtil.InitializeIntegrator(OptionsInfo["IntegratorParams"], OptionsInfo["SystemParams"]["ConstraintErrorTolerance"])
120
121 return Integrator
122
123 def SetupSimulation(System, Integrator, Topology, Positions):
124 """Setup simulation. """
125
126 Simulation = OpenMMUtil.InitializeSimulation(System, Integrator, Topology, Positions, OptionsInfo["PlatformParams"])
127
128 return Simulation
129
130 def PerformEnergyMinimization(Simulation):
131 """Perform energy minimization."""
132
133 SimulationParams = OpenMMUtil.SetupSimulationParameters(OptionsInfo["SimulationParams"])
134
135 OutputParams = OptionsInfo["OutputParams"]
136
137 # Setup a local minimization reporter...
138 MinimizeReporter = None
139 if OutputParams["MinimizationDataStdout"] or OutputParams["MinimizationDataLog"]:
140 MinimizeReporter = LocalMinimizationReporter()
141
142 if MinimizeReporter is not None:
143 MiscUtil.PrintInfo("\nAdding minimization reporters...")
144 if OutputParams["MinimizationDataLog"]:
145 MiscUtil.PrintInfo("Adding data log minimization reporter (Steps: %s; File: %s)..." % (OutputParams["MinimizationDataSteps"], OutputParams["MinimizationDataLogFile"]))
146 if OutputParams["MinimizationDataStdout"]:
147 MiscUtil.PrintInfo("Adding data stdout minimization reporter (Steps: %s)..." % (OutputParams["MinimizationDataSteps"]))
148 else:
149 MiscUtil.PrintInfo("\nSkipping addition of minimization reporters...")
150
151 MaxSteps = SimulationParams["MinimizationMaxSteps"]
152
153 MaxStepsMsg = "MaxSteps: %s" % ("UntilConverged" if MaxSteps == 0 else MaxSteps)
154 ToleranceMsg = "Tolerance: %.2f kcal/mol/A (%.2f kjoules/mol/nm)" % (SimulationParams["MinimizationToleranceInKcal"], SimulationParams["MinimizationToleranceInJoules"])
155
156 MiscUtil.PrintInfo("\nPerforming energy minimization (%s; %s)..." % (MaxStepsMsg, ToleranceMsg))
157
158 if OutputParams["MinimizationDataStdout"]:
159 HeaderLine = SetupMinimizationDataOutHeaderLine()
160 print("\n%s" % HeaderLine)
161
162 Simulation.minimizeEnergy(tolerance = SimulationParams["MinimizationTolerance"], maxIterations = MaxSteps, reporter = MinimizeReporter)
163
164 if OutputParams["MinimizationDataLog"]:
165 WriteMinimizationDataLogFile(MinimizeReporter.DataOutValues)
166
167 # Write out minimized structure...
168 MiscUtil.PrintInfo("\nWriting PDB file %s..." % OptionsInfo["MinimizedPDBOutfile"])
169 OpenMMUtil.WriteSimulationStatePDBFile(Simulation, OptionsInfo["MinimizedPDBOutfile"], OptionsInfo["OutputParams"]["PDBOutKeepIDs"])
170
171 class LocalMinimizationReporter(mm.MinimizationReporter):
172 """Setup a local minimization reporter. """
173
174 (DataSteps, DataOutTypeList, DataOutDelimiter, StdoutStatus) = [None] * 4
175
176 DataOutValues = []
177 First = True
178
179 def report(self, Iteration, PositonsList, GradientsList, DataStatisticsMap):
180 """Report and track minimization."""
181
182 if self.First:
183 # Initialize...
184 self.DataSteps = OptionsInfo["OutputParams"]["MinimizationDataSteps"]
185 self.DataOutTypeList = OptionsInfo["OutputParams"]["MinimizationDataOutTypeOpenMMNameList"]
186 self.DataOutDelimiter = OptionsInfo["OutputParams"]["DataOutDelimiter"]
187 self.StdoutStatus = True if OptionsInfo["OutputParams"]["MinimizationDataStdout"] else False
188
189 self.First = False
190
191 if Iteration % self.DataSteps == 0:
192 # Setup data values...
193 DataValues = []
194 DataValues.append("%s" % Iteration)
195 for DataType in self.DataOutTypeList:
196 DataValue = "%.4f" % DataStatisticsMap[DataType]
197 DataValues.append(DataValue)
198
199 # Track data...
200 self.DataOutValues.append(DataValues)
201
202 # Print data values...
203 if self.StdoutStatus:
204 print("%s" % self.DataOutDelimiter.join(DataValues))
205
206 # This method must return a bool. You may return true for early termination.
207 return False
208
209 def WriteMinimizationDataLogFile(DataOutValues):
210 """Write minimization data log file. """
211
212 OutputParams = OptionsInfo["OutputParams"]
213
214 Outfile = OutputParams["MinimizationDataLogFile"]
215 OutDelimiter = OutputParams["DataOutDelimiter"]
216
217 MiscUtil.PrintInfo("\nWriting minimization log file %s..." % Outfile)
218 OutFH = open(Outfile, "w")
219
220 HeaderLine = SetupMinimizationDataOutHeaderLine()
221 OutFH.write("%s\n" % HeaderLine)
222
223 for LineWords in DataOutValues:
224 Line = OutDelimiter.join(LineWords)
225 OutFH.write("%s\n" % Line)
226
227 OutFH.close()
228
229 def SetupMinimizationDataOutHeaderLine():
230 """Setup minimization data output header line. """
231
232 LineWords = ["Iteration"]
233 for Label in OptionsInfo["OutputParams"]["MinimizationDataOutTypeList"]:
234 if re.match("^(SystemEnergy|RestraintEnergy)$", Label, re.I):
235 LineWords.append("%s(kjoules/mol)" % Label)
236 elif re.match("^RestraintStrength$", Label, re.I):
237 LineWords.append("%s(kjoules/mol/nm^2)" % Label)
238 else:
239 LineWords.append(Label)
240
241 Line = OptionsInfo["OutputParams"]["DataOutDelimiter"].join(LineWords)
242
243 return Line
244
245 def FreezeRestraintAtoms(System, Topology, Positions):
246 """Handle freezing and restraining of atoms."""
247
248 # Get atoms for freezing...
249 FreezeAtomList = None
250 if OptionsInfo["FreezeAtoms"]:
251 FreezeAtomList = OpenMMUtil.GetAtoms(Topology, OptionsInfo["FreezeAtomsParams"]["CAlphaProteinStatus"], OptionsInfo["FreezeAtomsParams"]["ResidueNames"], OptionsInfo["FreezeAtomsParams"]["Negate"])
252 if FreezeAtomList is None:
253 MiscUtil.PrintError("The freeze atoms parameters specified, \"selection, %s, selectionSpec, %s, negate, %s\", - using \"--freezeAtomsParams\" option didn't match any atoms in the system. You must specify a valid set of parameters for freezing atoms or disable freezing using \"No\" value for \"--freezeAtoms\" option.." % (OptionsInfo["FreezeAtomsParams"]["Selection"], OptionsInfo["FreezeAtomsParams"]["SelectionSpec"], OptionsInfo["FreezeAtomsParams"]["Negate"]))
254
255 # Get atoms for restraining...
256 RestraintAtomList = None
257 if OptionsInfo["RestraintAtoms"]:
258 RestraintAtomList = OpenMMUtil.GetAtoms(Topology, OptionsInfo["RestraintAtomsParams"]["CAlphaProteinStatus"], OptionsInfo["RestraintAtomsParams"]["ResidueNames"], OptionsInfo["RestraintAtomsParams"]["Negate"])
259 if RestraintAtomList is None:
260 MiscUtil.PrintError("The restraint atoms parameters specified, \"selection, %s, selectionSpec, %s, negate, %s\", - using \"--restraintAtomsParams\" option didn't match any atoms in the system. You must specify a valid set of parameters for restraining atoms or disable restraining using \"No\" value for \"--restraintAtoms\" option." % (OptionsInfo["RestraintAtomsParams"]["Selection"], OptionsInfo["RestraintAtomsParams"]["SelectionSpec"], OptionsInfo["RestraintAtomsParams"]["Negate"]))
261
262 # Check for atoms to freeze or restraint...
263 if FreezeAtomList is None and RestraintAtomList is None:
264 return
265
266 # Check for overlap between freeze and restraint atoms...
267 if OpenMMUtil.DoAtomListsOverlap(FreezeAtomList, RestraintAtomList):
268 MiscUtil.PrintError("The atoms specified using \"--freezeAtomsParams\" and \"--restraintAtomsParams\" options appear to overlap. You must specify unique sets of atoms to freeze and restraint.")
269
270 # Check overlap of freeze atoms with system constraints...
271 if OpenMMUtil.DoesAtomListOverlapWithSystemConstraints(System, FreezeAtomList):
272 MiscUtil.PrintError("The atoms specified using \"--freezeAtomsParams\" appear to overlap with atoms being constrained corresponding to the value specified, %s, for paramater name \"constraints\" using \"--systemParams\" option.\n\nYou must specify a unique set of atoms to freeze or turn off system constaints by specifying value, None, for \"constraints\" parameter using option \"--systemsParams\".\n\nIn addtion, you may want specify, no, value for \"rigidWater\" option.\n\nThe atoms are frozen by setting their particle mass to zero. OpenMM doesn't allow to both constraint atoms and set their mass to zero." % (OptionsInfo["SystemParams"]["Constraints"]))
273
274 # Check overlap of restraint atoms with system constraints...
275 if OpenMMUtil.DoesAtomListOverlapWithSystemConstraints(System, RestraintAtomList):
276 MiscUtil.PrintInfo("")
277 MiscUtil.PrintWarning("The atoms specified using \"--restraintAtomsParams\" appear to overlap with atoms being constrained corresponding to the value specified, %s, for paramater name \"constraints\" using \"--systemParams\" option. You may want to specify a unique set of atoms to restraints or turn off system constaints by specifying value, None, for \"constraints\" parameter using option \"--systemsParams\"." % (OptionsInfo["SystemParams"]["Constraints"]))
278
279 # Freeze atoms...
280 if FreezeAtomList is None:
281 MiscUtil.PrintInfo("\nSkipping freezing of atoms...")
282 else:
283 MiscUtil.PrintInfo("\nFreezing atoms (Selection: %s)..." % OptionsInfo["FreezeAtomsParams"]["Selection"])
284 OpenMMUtil.FreezeAtoms(System, FreezeAtomList)
285
286 # Restraint atoms...
287 if RestraintAtomList is None:
288 MiscUtil.PrintInfo("\nSkipping restraining of atoms...")
289 else:
290 MiscUtil.PrintInfo("\nRestraining atoms (Selection: %s)..." % OptionsInfo["RestraintAtomsParams"]["Selection"])
291
292 SprintConstantInKcal = OptionsInfo["RestraintSpringConstant"]
293 OpenMMUtil.RestraintAtoms(System, Positions, RestraintAtomList, SprintConstantInKcal)
294
295 def ProcessOutfilePrefixParameter():
296 """Process outfile prefix paramater."""
297
298 OutfilePrefix = Options["--outfilePrefix"]
299
300 if not re.match("^auto$", OutfilePrefix, re.I):
301 OptionsInfo["OutfilePrefix"] = OutfilePrefix
302 return
303
304 if OptionsInfo["SmallMolFileMode"]:
305 OutfilePrefix = "%s_%s_Complex" % (OptionsInfo["InfileRoot"], OptionsInfo["SmallMolFileRoot"])
306 else:
307 OutfilePrefix = "%s" % (OptionsInfo["InfileRoot"])
308
309 if re.match("^yes$", Options["--waterBox"], re.I):
310 OutfilePrefix = "%s_Solvated" % (OutfilePrefix)
311
312 OptionsInfo["OutfilePrefix"] = OutfilePrefix
313
314 def ProcessOutfileNames():
315 """Process outfile names."""
316
317 OutputParams = OptionsInfo["OutputParams"]
318
319 OutfileParamName = "MinimizationDataLogFile"
320 OutfileParamValue = OutputParams[OutfileParamName]
321 if not Options["--overwrite"]:
322 if os.path.exists(OutfileParamValue):
323 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))
324
325 InitialPDBOutfile = "%s_Initial.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
326 MinimizedPDBOutfile = "%s_Minimized.%s" % (OptionsInfo["OutfilePrefix"], OutputParams["PDBOutfileExt"])
327 for Outfile in [InitialPDBOutfile, MinimizedPDBOutfile]:
328 if not Options["--overwrite"]:
329 if os.path.exists(Outfile):
330 MiscUtil.PrintError("The file name, %s, generated using option \"--outfilePrefix\" already exist. Use option \"--ov\" or \"--overwrite\" and try again. " % (Outfile))
331 OptionsInfo["InitialPDBOutfile"] = InitialPDBOutfile
332 OptionsInfo["MinimizedPDBOutfile"] = MinimizedPDBOutfile
333
334 def ProcessWaterBoxParameters():
335 """Process water box parameters. """
336
337 OptionsInfo["WaterBox"] = True if re.match("^yes$", Options["--waterBox"], re.I) else False
338 OptionsInfo["WaterBoxParams"] = OpenMMUtil.ProcessOptionOpenMMWaterBoxParameters("--waterBoxParams", Options["--waterBoxParams"])
339
340 if OptionsInfo["WaterBox"]:
341 if OptionsInfo["ForcefieldParams"]["ImplicitWater"]:
342 MiscUtil.PrintInfo("")
343 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"]))
344
345 def ProcessOptions():
346 """Process and validate command line arguments and options."""
347
348 MiscUtil.PrintInfo("Processing options...")
349
350 ValidateOptions()
351
352 OptionsInfo["Infile"] = Options["--infile"]
353 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
354 OptionsInfo["InfileRoot"] = FileName
355
356 SmallMolFile = Options["--smallMolFile"]
357 SmallMolID = Options["--smallMolID"]
358 SmallMolFileMode = False
359 SmallMolFileRoot = None
360 if SmallMolFile is not None:
361 FileDir, FileName, FileExt = MiscUtil.ParseFileName(SmallMolFile)
362 SmallMolFileRoot = FileName
363 SmallMolFileMode = True
364
365 OptionsInfo["SmallMolFile"] = SmallMolFile
366 OptionsInfo["SmallMolFileRoot"] = SmallMolFileRoot
367 OptionsInfo["SmallMolFileMode"] = SmallMolFileMode
368 OptionsInfo["SmallMolID"] = SmallMolID.upper()
369
370 ProcessOutfilePrefixParameter()
371
372 ParamsDefaultInfoOverride = {"MinimizationDataSteps": 100, "MinimizationDataStdout": False, "MinimizationDataLog": True}
373 OptionsInfo["OutputParams"] = OpenMMUtil.ProcessOptionOpenMMOutputParameters("--outputParams", Options["--outputParams"], OptionsInfo["OutfilePrefix"], ParamsDefaultInfoOverride)
374 ProcessOutfileNames()
375
376 OptionsInfo["ForcefieldParams"] = OpenMMUtil.ProcessOptionOpenMMForcefieldParameters("--forcefieldParams", Options["--forcefieldParams"])
377
378 OptionsInfo["FreezeAtoms"] = True if re.match("^yes$", Options["--freezeAtoms"], re.I) else False
379 if OptionsInfo["FreezeAtoms"]:
380 OptionsInfo["FreezeAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters("--freezeAtomsParams", Options["--freezeAtomsParams"])
381 else:
382 OptionsInfo["FreezeAtomsParams"] = None
383
384 ParamsDefaultInfoOverride = {"Name": Options["--platform"], "Threads": 1}
385 OptionsInfo["PlatformParams"] = OpenMMUtil.ProcessOptionOpenMMPlatformParameters("--platformParams", Options["--platformParams"], ParamsDefaultInfoOverride)
386
387 OptionsInfo["RestraintAtoms"] = True if re.match("^yes$", Options["--restraintAtoms"], re.I) else False
388 if OptionsInfo["RestraintAtoms"]:
389 OptionsInfo["RestraintAtomsParams"] = OpenMMUtil.ProcessOptionOpenMMAtomsSelectionParameters("--restraintAtomsParams", Options["--restraintAtomsParams"])
390 else:
391 OptionsInfo["RestraintAtomsParams"] = None
392 OptionsInfo["RestraintSpringConstant"] = float(Options["--restraintSpringConstant"])
393
394 OptionsInfo["SystemParams"] = OpenMMUtil.ProcessOptionOpenMMSystemParameters("--systemParams", Options["--systemParams"])
395
396 OptionsInfo["IntegratorParams"] = OpenMMUtil.ProcessOptionOpenMMIntegratorParameters("--integratorParams", Options["--integratorParams"], HydrogenMassRepartioningStatus = OptionsInfo["SystemParams"]["HydrogenMassRepartioning"])
397
398 OptionsInfo["SimulationParams"] = OpenMMUtil.ProcessOptionOpenMMSimulationParameters("--simulationParams", Options["--simulationParams"])
399
400 ProcessWaterBoxParameters()
401
402 OptionsInfo["Overwrite"] = Options["--overwrite"]
403
404 def RetrieveOptions():
405 """Retrieve command line arguments and options."""
406
407 # Get options...
408 global Options
409 Options = docopt(_docoptUsage_)
410
411 # Set current working directory to the specified directory...
412 WorkingDir = Options["--workingdir"]
413 if WorkingDir:
414 os.chdir(WorkingDir)
415
416 # Handle examples option...
417 if "--examples" in Options and Options["--examples"]:
418 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
419 sys.exit(0)
420
421 def ValidateOptions():
422 """Validate option values."""
423
424 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
425 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
426
427 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--infile"])
428 OutfilePrefix = Options["--outfilePrefix"]
429 if not re.match("^auto$", OutfilePrefix, re.I):
430 if re.match("^(%s)$" % OutfilePrefix, FileName, re.I):
431 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))
432
433 if Options["--smallMolFile"] is not None:
434 MiscUtil.ValidateOptionFilePath("-l, --smallMolFile", Options["--smallMolFile"])
435 MiscUtil.ValidateOptionFileExt("-l, --smallMolFile", Options["--smallMolFile"], "sd sdf")
436
437 SmallMolID = Options["--smallMolID"]
438 if len(SmallMolID) != 3:
439 MiscUtil.PrintError("The value specified, %s, for option \"--smallMolID\" is not valid. You must specify a three letter small molecule ID." % (SmallMolID))
440
441 MiscUtil.ValidateOptionTextValue("--freezeAtoms", Options["--freezeAtoms"], "yes no")
442 if re.match("^yes$", Options["--freezeAtoms"], re.I):
443 if Options["--freezeAtomsParams"] is None:
444 MiscUtil.PrintError("No value specified for option \"--freezeAtomsParams\". You must specify valid values during, yes, value for \"--freezeAtoms\" option.")
445
446 MiscUtil.ValidateOptionTextValue("-p, --platform", Options["--platform"], "CPU CUDA OpenCL Reference")
447
448 MiscUtil.ValidateOptionTextValue("--restraintAtoms", Options["--restraintAtoms"], "yes no")
449 if re.match("^yes$", Options["--restraintAtoms"], re.I):
450 if Options["--restraintAtomsParams"] is None:
451 MiscUtil.PrintError("No value specified for option \"--restraintAtomsParams\". You must specify valid values during, yes, value for \"--restraintAtoms\" option.")
452
453 MiscUtil.ValidateOptionFloatValue("--restraintSpringConstant", Options["--restraintSpringConstant"], {">": 0})
454
455 MiscUtil.ValidateOptionTextValue("--waterBox", Options["--waterBox"], "yes no")
456
457 # Setup a usage string for docopt...
458 _docoptUsage_ = """
459 OpenMMPerformMinimization.py - Perform an energy minimization.
460
461 Usage:
462 OpenMMPerformMDSimulation.py [--forcefieldParams <Name,Value,..>] [--freezeAtoms <yes or no>]
463 [--freezeAtomsParams <Name,Value,..>] [--integratorParams <Name,Value,..>]
464 [--outputParams <Name,Value,..>] [--outfilePrefix <text>]
465 [--overwrite] [--platform <text>] [--platformParams <Name,Value,..>]
466 [--restraintAtoms <yes or no>]
467 [--restraintAtomsParams <Name,Value,..>] [--restraintSpringConstant <number>]
468 [--simulationParams <Name,Value,..>] [--smallMolFile <SmallMolFile>] [--smallMolID <text>]
469 [--systemParams <Name,Value,..>] [--waterBox <yes or no>]
470 [--waterBoxParams <Name,Value,..>] [-w <dir>] -i <infile>
471 OpenMMPerformMDSimulation.py -h | --help | -e | --examples
472
473 Description:
474 Perform energy minimization for a macromolecule or a macromolecule in a
475 complex with small molecule. You may optionally add a water box and
476 freeze/restraint atoms to your system before minimization.
477
478 The input file must contain a macromolecule already prepared for simulation.
479 The preparation of the macromolecule for a simulation generally involves the
480 following: tasks: Identification and replacement non-standard residues;
481 Addition of missing residues; Addition of missing heavy atoms; Addition of
482 missing hydrogens; Addition of a water box which is optional.
483
484 In addition, the small molecule input file must contain a molecule already
485 prepared for simulation. It must contain appropriate 3D coordinates relative
486 to the macromolecule along with no missing hydrogens.
487
488 The supported macromolecule input file formats are: PDB (.pdb) and
489 CIF (.cif)
490
491 The supported small molecule input file format are : SD (.sdf, .sd)
492
493 Possible outfile prefixes:
494
495 <InfieRoot>
496 <InfieRoot>_Solvated
497 <InfieRoot>_<SmallMolFileRoot>_Complex
498 <InfieRoot>_<SmallMolFileRoot>_Complex_Solvated
499
500 Possible output files:
501
502 <OutfilePrefix>_Initial.<pdb or cif>
503 <OutfilePrefix>_Minimized.<pdb or cif>
504 <OutfilePrefix>_Minimization.csv
505
506 Options:
507 -e, --examples
508 Print examples.
509 -f, --forcefieldParams <Name,Value,..> [default: auto]
510 A comma delimited list of parameter name and value pairs for biopolymer,
511 water, and small molecule forcefields.
512
513 The supported parameter names along with their default values are
514 are shown below:
515
516 biopolymer, amber14-all.xml [ Possible values: Any Valid value ]
517 smallMolecule, openff-2.2.1 [ Possible values: Any Valid value ]
518 water, auto [ Possible values: Any Valid value ]
519 additional, none [ Possible values: Space delimited list of any
520 valid value ]
521
522 Possible biopolymer forcefield values:
523
524 amber14-all.xml, amber99sb.xml, amber99sbildn.xml, amber03.xml,
525 amber10.xml
526 charmm36.xml, charmm_polar_2019.xml
527 amoeba2018.xml
528
529 Possible small molecule forcefield values:
530
531 openff-2.2.1, openff-2.0.0, openff-1.3.1, openff-1.2.1,
532 openff-1.1.1, openff-1.1.0,...
533 smirnoff99Frosst-1.1.0, smirnoff99Frosst-1.0.9,...
534 gaff-2.11, gaff-2.1, gaff-1.81, gaff-1.8, gaff-1.4,...
535
536 The default water forcefield valus is dependent on the type of the
537 biopolymer forcefield as shown below:
538
539 Amber: amber14/tip3pfb.xml
540 CHARMM: charmm36/water.xml or None for charmm_polar_2019.xml
541 Amoeba: None (Explicit)
542
543 Possible water forcefield values:
544
545 amber14/tip3p.xml, amber14/tip3pfb.xml, amber14/spce.xml,
546 amber14/tip4pew.xml, amber14/tip4pfb.xml,
547 charmm36/water.xml, charmm36/tip3p-pme-b.xml,
548 charmm36/tip3p-pme-f.xml, charmm36/spce.xml,
549 charmm36/tip4pew.xml, charmm36/tip4p2005.xml,
550 charmm36/tip5p.xml, charmm36/tip5pew.xml,
551 implicit/obc2.xml, implicit/GBn.xml, implicit/GBn2.xml,
552 amoeba2018_gk.xml (Implict water)
553 None (Explicit water for amoeba)
554
555 The additional forcefield value is a space delimited list of any valid
556 forcefield values and is passed on to the OpenMMForcefields
557 SystemGenerator along with the specified forcefield values for
558 biopolymer, water, and mall molecule. Possible additional forcefield
559 values are:
560
561 amber14/DNA.OL15.xml amber14/RNA.OL3.xml
562 amber14/lipid17.xml amber14/GLYCAM_06j-1.xml
563 ... ... ...
564
565 You may specify any valid forcefield names supported by OpenMM. No
566 explicit validation is performed.
567 --freezeAtoms <yes or no> [default: no]
568 Freeze atoms during energy minimization. The specified atoms are kept
569 completely fixed by setting their masses to zero. Their positions do not
570 change during energy minimization.
571 --freezeAtomsParams <Name,Value,..>
572 A comma delimited list of parameter name and value pairs for freezing
573 atoms during energy minimization. You must specify these parameters for 'yes'
574 value of '--freezeAtoms' option.
575
576 The supported parameter names along with their default values are
577 are shown below:
578
579 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
580 Protein, Residues, or Water ]
581 selectionSpec, auto [ Possible values: A space delimited list of
582 residue names ]
583 negate, no [ Possible values: yes or no ]
584
585 A brief description of parameters is provided below:
586
587 selection: Atom selection to freeze.
588 selectionSpec: A space delimited list of residue names for
589 selecting atoms to freeze. You must specify its value during
590 'Ligand' and 'Protein' value for 'selection'. The default values
591 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
592 and 'Water' values of 'selection' as shown below:
593
594 CAlphaProtein: List of stadard protein residues from pdbfixer
595 for selecting CAlpha atoms.
596 Ions: Li Na K Rb Cs Cl Br F I
597 Water: HOH
598 Protein: List of standard protein residues from pdbfixer.
599
600 negate: Negate atom selection match to select atoms for freezing.
601
602 In addition, you may specify an explicit space delimited list of residue
603 names using 'selectionSpec' for any 'selection". The specified residue
604 names are appended to the appropriate default values during the
605 selection of atoms for freezing.
606 -h, --help
607 Print this help message.
608 -i, --infile <infile>
609 Input file name containing a macromolecule.
610 --integratorParams <Name,Value,..> [default: auto]
611 A comma delimited list of parameter name and value pairs for integrator
612 to setup the system for local energy minimization. No MD simulation is
613 performed.
614
615 The supported parameter names along with their default values are
616 are shown below:
617
618 temperature, 300.0 [ Units: kelvin ]
619
620 --outfilePrefix <text> [default: auto]
621 File prefix for generating the names of output files. The default value
622 depends on the names of input files for macromolecule and small molecule
623 along with the type of statistical ensemble and the nature of the solvation.
624
625 The possible values for outfile prefix are shown below:
626
627 <InfieRoot>
628 <InfieRoot>_Solvated
629 <InfieRoot>_<SmallMolFileRoot>_Complex
630 <InfieRoot>_<SmallMolFileRoot>_Complex_Solvated
631
632 --outputParams <Name,Value,..> [default: auto]
633 A comma delimited list of parameter name and value pairs for generating
634 output during energy minimization..
635
636 The supported parameter names along with their default values are
637 are shown below:
638
639 minimizationDataSteps, 100
640 minimizationDataStdout, no [ Possible values: yes or no ]
641 minimizationDataLog, yes [ Possible values: yes or no ]
642 minimizationDataLogFile, auto [ Default:
643 <OutfilePrefix>_MinimizationOut.csv ]
644 minimizationDataOutType, auto [ Possible values: A space delimited
645 list of valid parameter names. Default: SystemEnergy
646 RestraintEnergy MaxConstraintError.
647 Other valid names: RestraintStrength ]
648
649 pdbOutFormat, PDB [ Possible values: PDB or CIF ]
650 pdbOutKeepIDs, yes [ Possible values: yes or no ]
651
652 A brief description of parameters is provided below:
653
654 minimizationDataSteps: Frequency of writing data to stdout
655 and log file.
656 minimizationDataStdout: Write data to stdout.
657 minimizationDataLog: Write data to log file.
658 minimizationDataLogFile: Data log fie name.
659 minimizationDataOutType: Type of data to write to stdout
660 and log file.
661
662 pdbOutFormat: Format of output PDB files.
663 pdbOutKeepIDs: Keep existing chain and residue IDs.
664
665 --overwrite
666 Overwrite existing files.
667 -p, --platform <text> [default: CPU]
668 Platform to use for running MD simulation. Possible values: CPU, CUDA,
669 OpenCL, or Reference.
670 --platformParams <Name,Value,..> [default: auto]
671 A comma delimited list of parameter name and value pairs to configure
672 platform for running energy minimization calculations..
673
674 The supported parameter names along with their default values for
675 different platforms are shown below:
676
677 CPU:
678
679 threads, 1 [ Possible value: >= 0 or auto. The value of 'auto'
680 or zero implies the use of all available CPUs for threading. ]
681
682 CUDA:
683
684 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
685 deterministicForces, auto [ Possible values: yes or no ]
686 precision, single [ Possible values: single, double, or mix ]
687 tempDirectory, auto [ Possible value: DirName ]
688 useBlockingSync, auto [ Possible values: yes or no ]
689 useCpuPme, auto [ Possible values: yes or no ]
690
691 OpenCL:
692
693 deviceIndex, auto [ Possible values: 0, '0 1' etc. ]
694 openCLPlatformIndex, auto [ Possible value: Number]
695 precision, single [ Possible values: single, double, or mix ]
696 useCpuPme, auto [ Possible values: yes or no ]
697
698 A brief description of parameters is provided below:
699
700 CPU:
701
702 threads: Number of threads to use for simulation.
703
704 CUDA:
705
706 deviceIndex: Space delimited list of device indices to use for
707 calculations.
708 deterministicForces: Generate reproducible results at the cost of a
709 small decrease in performance.
710 precision: Number precision to use for calculations.
711 tempDirectory: Directory name for storing temporary files.
712 useBlockingSync: Control run-time synchronization between CPU and
713 GPU.
714 useCpuPme: Use CPU-based PME implementation.
715
716 OpenCL:
717
718 deviceIndex: Space delimited list of device indices to use for
719 simulation.
720 openCLPlatformIndex: Platform index to use for calculations.
721 precision: Number precision to use for calculations.
722 useCpuPme: Use CPU-based PME implementation.
723
724 --restraintAtoms <yes or no> [default: no]
725 Restraint atoms during energy minimization. The motion of specified atoms is
726 restricted by adding a harmonic force that binds them to their starting
727 positions. The atoms are not completely fixed unlike freezing of atoms.
728 Their motion, however, is restricted and they are not able to move far away
729 from their starting positions during energy minimization.
730 --restraintAtomsParams <Name,Value,..>
731 A comma delimited list of parameter name and value pairs for restraining
732 atoms during energy minimization. You must specify these parameters for 'yes'
733 value of '--restraintAtoms' option.
734
735 The supported parameter names along with their default values are
736 are shown below:
737
738 selection, none [ Possible values: CAlphaProtein, Ions, Ligand,
739 Protein, Residues, or Water ]
740 selectionSpec, auto [ Possible values: A space delimited list of
741 residue names ]
742 negate, no [ Possible values: yes or no ]
743
744 A brief description of parameters is provided below:
745
746 selection: Atom selection to restraint.
747 selectionSpec: A space delimited list of residue names for
748 selecting atoms to restraint. You must specify its value during
749 'Ligand' and 'Protein' value for 'selection'. The default values
750 are automatically set for 'CAlphaProtein', 'Ions', 'Protein',
751 and 'Water' values of 'selection' as shown below:
752
753 CAlphaProtein: List of stadard protein residues from pdbfixer
754 for selecting CAlpha atoms.
755 Ions: Li Na K Rb Cs Cl Br F I
756 Water: HOH
757 Protein: List of standard protein residues from pdbfixer.
758
759 negate: Negate atom selection match to select atoms for freezing.
760
761 In addition, you may specify an explicit space delimited list of residue
762 names using 'selectionSpec' for any 'selection". The specified residue
763 names are appended to the appropriate default values during the
764 selection of atoms for restraining.
765 --restraintSpringConstant <number> [default: 2.5]
766 Restraint spring constant for applying external restraint force to restraint
767 atoms relative to their initial positions during 'yes' value of '--restraintAtoms'
768 option. Default units: kcal/mol/A**2. The default value, 2.5, corresponds to
769 1046.0 kjoules/mol/nm**2. The default value is automatically converted into
770 units of kjoules/mol/nm**2 before its usage.
771 --simulationParams <Name,Value,..> [default: auto]
772 A comma delimited list of parameter name and value pairs for simulation.
773
774 The supported parameter names along with their default values are
775 are shown below:
776
777 minimizationMaxSteps, auto [ Possible values: >= 0. The value of
778 zero implies until the minimization is converged. ]
779 minimizationTolerance, 0.24 [ Units: kcal/mol/A. The default value
780 0.25, corresponds to OpenMM default of value of 10.04
781 kjoules/mol/nm. It is automatically converted into OpenMM
782 default units before its usage. ]
783
784 A brief description of parameters is provided below:
785
786 minimizationMaxSteps: Maximum number of minimization steps. The
787 value of zero implies until the minimization is converged.
788 minimizationTolerance: Energy convergence tolerance during
789 minimization.
790
791 -s, --smallMolFile <SmallMolFile>
792 Small molecule input file name. The macromolecue and small molecule are
793 merged for simulation and the complex is written out to a PDB file.
794 --smallMolID <text> [default: LIG]
795 Three letter small molecule residue ID. The small molecule ID corresponds
796 to the residue name of the small molecule and is written out to a PDB file
797 containing the complex.
798 --systemParams <Name,Value,..> [default: auto]
799 A comma delimited list of parameter name and value pairs to configure
800 a system for energy minimization. No MD simulation is performed.
801
802 The supported parameter names along with their default values are
803 are shown below:
804
805 constraints, BondsInvolvingHydrogens [ Possible values: None,
806 WaterOnly, BondsInvolvingHydrogens, AllBonds, or
807 AnglesInvolvingHydrogens ]
808 constraintErrorTolerance, 0.000001
809 ewaldErrorTolerance, 0.0005
810
811 nonbondedMethodPeriodic, PME [ Possible values: NoCutoff,
812 CutoffNonPeriodic, or PME ]
813 nonbondedMethodNonPeriodic, NoCutoff [ Possible values:
814 NoCutoff or CutoffNonPeriodic]
815 nonbondedCutoff, 1.0 [ Units: nm ]
816
817 hydrogenMassRepartioning, yes [ Possible values: yes or no ]
818 hydrogenMass, 1.5 [ Units: amu]
819
820 removeCMMotion, yes [ Possible values: yes or no ]
821 rigidWater, auto [ Possible values: yes or no. Default: 'No' for
822 'None' value of constraints; Otherwise, yes ]
823
824 A brief description of parameters is provided below:
825
826 constraints: Type of system constraints to use for simulation. These
827 constraints are different from freezing and restraining of any
828 atoms in the system.
829 constraintErrorTolerance: Distance tolerance for constraints as a
830 fraction of the constrained distance.
831 ewaldErrorTolerance: Ewald error tolerance for a periodic system.
832
833 nonbondedMethodPeriodic: Nonbonded method to use during the
834 calculation of long range interactions for a periodic system.
835 nonbondedMethodNonPeriodic: Nonbonded method to use during the
836 calculation of long range interactions for a non-periodic system.
837 nonbondedCutoff: Cutoff distance to use for long range interactions
838 in both perioidic non-periodic systems.
839
840 hydrogenMassRepartioning: Use hydrogen mass repartioning. It
841 increases the mass of the hydrogen atoms attached to the heavy
842 atoms and decreasing the mass of the bonded heavy atom to
843 maintain constant system mass. This allows the use of larger
844 integration step size (4 fs) during a simulation.
845 hydrogenMass: Hydrogen mass to use during repartioning.
846
847 removeCMMotion: Remove all center of mass motion at every time step.
848 rigidWater: Keep water rigid during a simulation. This is determined
849 automatically based on the value of 'constraints' parameter.
850
851 --waterBox <yes or no> [default: no]
852 Add water box.
853 --waterBoxParams <Name,Value,..> [default: auto]
854 A comma delimited list of parameter name and value pairs for adding
855 a water box.
856
857 The supported parameter names along with their default values are
858 are shown below:
859
860 model, tip3p [ Possible values: tip3p, spce, tip4pew, tip5p or
861 swm4ndp ]
862 mode, Padding [ Possible values: Size or Padding ]
863 padding, 1.0
864 size, None [ Possible value: xsize ysize zsize ]
865 shape, cube [ Possible values: cube, dodecahedron, or octahedron ]
866 ionPositive, Na+ [ Possible values: Li+, Na+, K+, Rb+, or Cs+ ]
867 ionNegative, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
868 ionicStrength, 0.0
869
870 A brief description of parameters is provided below:
871
872 model: Water model to use for adding water box. The van der
873 Waals radii and atomic charges are determined using the
874 specified water forcefield. You must specify an appropriate
875 water forcefield. No validation is performed.
876 mode: Specify the size of the waterbox explicitly or calculate it
877 automatically for a macromolecule along with adding padding
878 around ther macromolecule.
879 padding: Padding around a macromolecule in nanometers for filling
880 box with water. It must be specified during 'Padding' value of
881 'mode' parameter.
882 size: A space delimited triplet of values corresponding to water
883 size in nanometers. It must be specified during 'Size' value of
884 'mode' parameter.
885 ionPositive: Type of positive ion to add during the addition of a
886 water box.
887 ionNegative: Type of negative ion to add during the addition of a
888 water box.
889 ionicStrength: Total concentration of both positive and negative
890 ions to add excluding the ions added to neutralize the system
891 during the addition of a water box.
892
893 -w, --workingdir <dir>
894 Location of working directory which defaults to the current directory.
895
896 Examples:
897 To perform energy minimization for a macromolecule in a PDB file until the
898 energy is converged, writing information to log file every 100 steps and
899 generate a PDB file for the minimized system, type:
900
901 % OpenMMPerformMinimization.py -i Sample13.pdb
902
903 To run the first example for writing information to both stdout and log file
904 every 100 steps and generate various output files, type:
905
906 % OpenMMPerformMinimization.py -i Sample13.pdb --outputParams
907 "minimizationDataStdout, yes"
908
909 To run the first example for performing OpenMM simulation using multi-
910 threading employing all available CPUs on your machine and generate various
911 output files, type:
912
913 % OpenMMPerformMinimization.py -i Sample13.pdb
914 --platformParams "threads,0"
915
916 To run the first example for performing OpenMM simulation using CUDA platform
917 on your machine and generate various output files, type:
918
919 % OpenMMPerformMinimization.py -i Sample13.pdb -p CUDA
920
921 To run the second example for a marcomolecule in a complex with a small
922 molecule and generate various output files, type:
923
924 % OpenMMPerformMinimization.py -i Sample13.pdb -s Sample13Ligand.sdf
925 --platformParams "threads,0"
926
927 To run the second example by adding a water box to the system and generate
928 various output files, type:
929
930 % OpenMMPerformMinimization.py -i Sample13.pdb --waterBox yes
931 --platformParams "threads,0"
932
933 To run the second example for a marcomolecule in a complex with a small
934 molecule by adding a water box to the system and generate various output
935 files, type:
936
937 % OpenMMPerformMinimization.py -i Sample13.pdb -s Sample13Ligand.sdf
938 --waterBox yes --platformParams "threads,0"
939
940 To run the second example by freezing CAlpha atoms in a macromolecule without
941 using any system constraints to avoid any issues with the freezing of the same atoms
942 and generate various output files, type:
943
944 % OpenMMPerformMinimization.py -i Sample13.pdb
945 --freezeAtoms yes --freezeAtomsParams "selection,CAlphaProtein"
946 --systemParams "constraints, None"
947 --platformParams "threads,0"
948
949 % OpenMMPerformMinimization.py -i Sample13.pdb
950 --freezeAtoms yes --freezeAtomsParams "selection,CAlphaProtein"
951 --systemParams "constraints, None"
952 --platformParams "threads,0" --waterBox yes
953
954 To run the second example by restrainting CAlpha atoms in a macromolecule and
955 and generate various output files, type:
956
957 % OpenMMPerformMinimization.py -i Sample13.pdb
958 --restraintAtomsParams "selection,CAlphaProtein"
959 --platformParams "threads,0"
960
961 % OpenMMPerformMinimization.py -i Sample13.pdb
962 --restraintAtomsParams "selection,CAlphaProtein"
963 --platformParams "threads,0"
964 --waterBox yes
965
966 To run the second example by specifying explict values for various parametres
967 and generate various output files, type:
968
969 % OpenMMPerformMinimization.py -i Sample13.pdb
970 -f " biopolymer,amber14-all.xml,smallMolecule, openff-2.2.1,
971 water,amber14/tip3pfb.xml" --waterBox yes
972 --outputParams "minimizationDataSteps, 100, minimizationDataStdout,
973 yes,minimizationDataLog,yes"
974 -p CPU --platformParams "threads,0"
975 --simulationParams "minimizationMaxSteps,500,
976 minimizationTolerance, 0.24"
977 --systemParams "constraints,BondsInvolvingHydrogens,
978 nonbondedMethodPeriodic,PME,nonbondedMethodNonPeriodic,NoCutoff,
979 hydrogenMassRepartioning, yes"
980
981 Author:
982 Manish Sud(msud@san.rr.com)
983
984 See also:
985 OpenMMPrepareMacromolecule.py, OpenMMPerformMDSimulation.py,
986 OpenMMPerformSimulatedAnnealing.py
987
988 Copyright:
989 Copyright (C) 2025 Manish Sud. All rights reserved.
990
991 The functionality available in this script is implemented using OpenMM, an
992 open source molecuar simulation package.
993
994 This file is part of MayaChemTools.
995
996 MayaChemTools is free software; you can redistribute it and/or modify it under
997 the terms of the GNU Lesser General Public License as published by the Free
998 Software Foundation; either version 3 of the License, or (at your option) any
999 later version.
1000
1001 """
1002
1003 if __name__ == "__main__":
1004 main()
