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