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