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