1 #!/bin/env python
2 #
3 # File: OpenFECalculateRelativeBindingFreeEnergySepTop.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 OpenFE, an
9 # open source package for alchemical free energy calculations.
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 import logging
36 import pathlib
37 import pandas as pd
38
39 # OpenFE imports...
40 try:
41 import openfe
42 from openfe.protocols.openmm_septop import SepTopProtocol
43 except ImportError as ErrMsg:
44 sys.stderr.write("\nFailed to import OpenFE related module/package: %s\n" % ErrMsg)
45 sys.stderr.write("Check/update your OpenFE environment and try again.\n\n")
46 sys.exit(1)
47
48 # RDKit imports...
49 try:
50 from rdkit import rdBase
51 except ImportError as ErrMsg:
52 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
53 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
54 sys.exit(1)
55
56 # MayaChemTools imports...
57 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
58 try:
59 from docopt import docopt
60 import MiscUtil
61 import OpenFEUtil
62 except ImportError as ErrMsg:
63 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
64 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
65 sys.exit(1)
66
67 ScriptName = os.path.basename(sys.argv[0])
68 Options = {}
69 OptionsInfo = {}
70
71
72 def main():
73 """Start execution of the script."""
74
75 MiscUtil.PrintInfo(
76 "\n%s (OpenFE v%s; OpenMM v%s; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
77 % (
78 ScriptName,
79 openfe.version("openfe"),
80 openfe.version("openmm"),
81 rdBase.rdkitVersion,
82 MiscUtil.GetMayaChemToolsVersion(),
83 time.asctime(),
84 )
85 )
86
87 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
88
89 # Retrieve command line arguments and options...
90 RetrieveOptions()
91
92 if Options["--list"]:
93 ProcessListOption()
94 else:
95 # Process and validate command line arguments and options...
96 ProcessOptions()
97
98 # Perform actions required by the script...
99 CalculateRelativeBindingFreeEnergy()
100
101 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
102 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
103
104
105 def CalculateRelativeBindingFreeEnergy():
106 """Calculate relative binding free energy using separated topologies."""
107
108 # Process input files...
109 MacroMol, Mols = ProcessInputFiles()
110
111 # Process molecule names...
112 MolPairs = ProcessMoleculeNames(Mols)
113
114 # Check for missing partial charges...
115 CheckMissingPartialCharges(MolPairs)
116
117 # Initialize RBFE protocols...
118 RBFEProtocol = InitializeSeparatedTopologyProtocol()
119
120 # Initialize solvent...
121 Solvent = InitializeSolventComponent()
122
123 # Setup transformations...
124 MolAToMolBTransformations = SetupTransformations(MolPairs, MacroMol, Solvent, RBFEProtocol)
125
126 # Setup protocol DAGs...
127 MolAToMolBProtocolDAGs = SetupProtocolDAGs(MolAToMolBTransformations)
128
129 # Execute protocol DAGs and gather results...
130 MolAToMolBProtocolResults = ExecuteProtocolDAGsAndGatherResults(MolAToMolBTransformations, MolAToMolBProtocolDAGs)
131
132 # Process protocol results...
133 ProcessProtocolResults(MolAToMolBTransformations, MolAToMolBProtocolResults)
134
135
136 def InitializeSeparatedTopologyProtocol():
137 """Initialize separated toplology protocol."""
138
139 MiscUtil.PrintInfo("\nInitializing separated topologies protocol...")
140
141 RBFESettings = OpenFEUtil.SetupRelativeFreeEnergySeparatedTopologySettings(
142 "-r, --rbfeParams", OptionsInfo["RBFEParams"]
143 )
144 RBFEProtocol = OpenFEUtil.InitializeRelativeFreeEngerySeparatedTopologyProtocol(RBFESettings)
145
146 return RBFEProtocol
147
148
149 def InitializeSolventComponent():
150 """Initialize solvent component."""
151
152 SolventParams = OptionsInfo["SolventParams"]
153 MiscUtil.PrintInfo(
154 "\nInitializing solvent component (PositiveIon: %s; NegativeIon: %s; Neutralize: %s; IonConcentration: %s)..."
155 % (
156 SolventParams["PositiveIon"],
157 SolventParams["NegativeIon"],
158 SolventParams["Neutralize"],
159 SolventParams["IonConcentration"],
160 )
161 )
162
163 Solvent = OpenFEUtil.InitializeSolventComponent(SolventParams)
164
165 return Solvent
166
167
168 def SetupTransformations(MolPairs, MacroMol, Solvent, RBFEProtocol):
169 """Set up transformations for moleule pairs."""
170
171 MiscUtil.PrintInfo("\nSetting up transformations (Count: %s)..." % (int(len(MolPairs) / 2)))
172
173 MolAToMolBTransformations = []
174 for Index in range(0, len(MolPairs), 2):
175 MolA = MolPairs[Index]
176 MolB = MolPairs[Index + 1]
177
178 # Setup chemical systems...
179 MolAComplexSolventSystem = OpenFEUtil.InitializeChemicalSystem(
180 SmallMol=MolA, MacroMol=MacroMol, Solvent=Solvent, Name="%s_Complex_Solvent" % MolA.name
181 )
182 MolBComplexSolventSystem = OpenFEUtil.InitializeChemicalSystem(
183 SmallMol=MolB, MacroMol=MacroMol, Solvent=Solvent, Name="%s_Complex_Solvent" % MolB.name
184 )
185
186 # Setup MolAComplexSolvent to MolBComplexSolvent transformation...
187 TransformationName = "%s_To_%s_Complex_Solvent" % (MolA.name, MolB.name)
188 MolAToMolBComplexSolventTransformation = OpenFEUtil.InitializeTransformation(
189 StateA=MolAComplexSolventSystem,
190 StateB=MolBComplexSolventSystem,
191 Mapping=None,
192 Protocol=RBFEProtocol,
193 Name=TransformationName,
194 Validate=False,
195 )
196 MolAToMolBTransformations.append(MolAToMolBComplexSolventTransformation)
197
198 # Write out transformatios...
199 WriteTransformations(MolAToMolBTransformations)
200
201 return MolAToMolBTransformations
202
203
204 def WriteTransformations(MolAToMolBTransformations):
205 """Write out transformations."""
206
207 TransformationsOutDirPath = pathlib.Path(OptionsInfo["TransformationsOutDirPath"])
208
209 MiscUtil.PrintInfo(
210 "Writing transformations files (Files: *.json; Count: %s; Subdirectory: %s)..."
211 % (len(MolAToMolBTransformations), OptionsInfo["TransformationsOutDir"])
212 )
213
214 for Transformation in MolAToMolBTransformations:
215 TransformationFilePath = TransformationsOutDirPath.joinpath("%s.json" % Transformation.name)
216 Transformation.dump(TransformationFilePath)
217
218
219 def SetupProtocolDAGs(MolAToMolBTransformations):
220 """Setup protocol Directed Acyclic Graphs (DAGs) for each transformation to
221 to perform calculations.
222 """
223
224 MiscUtil.PrintInfo("\nSetting up protocol DAGs (Count: %s)..." % len(MolAToMolBTransformations))
225
226 MolAToMolBProtocolDAGs = []
227 for Transformation in MolAToMolBTransformations:
228 ProtocolDAG = OpenFEUtil.InitializeProtocolDAG(Transformation, Name=Transformation.name)
229 MolAToMolBProtocolDAGs.append(ProtocolDAG)
230
231 return MolAToMolBProtocolDAGs
232
233
234 def ExecuteProtocolDAGsAndGatherResults(MolTransformations, MolProtocolDAGs):
235 """Execute protocol DAGs and gather results."""
236
237 ResultsSharedOutDirPath = OptionsInfo["ResultsOutDirPath"]
238 ResultsScratchOutDirPath = OptionsInfo["ResultsScratchOutDirPath"]
239 ExecuteDAGParams = OptionsInfo["ExecuteDAGParams"]
240
241 MolProtocolResults = OpenFEUtil.ExecuteProtocolDAGsAndGatherResults(
242 MolTransformations,
243 MolProtocolDAGs,
244 ResultsSharedOutDirPath,
245 ResultsScratchOutDirPath,
246 KeepShared=ExecuteDAGParams["KeepShared"],
247 KeepScratch=ExecuteDAGParams["KeepScratch"],
248 NRetries=ExecuteDAGParams["NRetries"],
249 WriteResults=True,
250 )
251
252 return MolProtocolResults
253
254
255 def ProcessProtocolResults(MolAToMolBTransformations, MolAToMolBProtocolResults):
256 """Process protocol results."""
257
258 ResultFileParams = OptionsInfo["ResultFileParams"]
259
260 ResultFile = "%s_RBFE_Results.%s" % (OptionsInfo["OutfilePrefix"], ResultFileParams["Ext"])
261 ResultFilePath = os.path.join(OptionsInfo["OutfileDirPath"], ResultFile)
262 MiscUtil.PrintInfo("\nWriting %s..." % ResultFile)
263
264 Precision = ResultFileParams["Precision"]
265
266 ResultData = []
267 for Index in range(0, len(MolAToMolBProtocolResults), 1):
268 MolAToMolBProtocolResult = MolAToMolBProtocolResults[Index]
269
270 # Setup mol names using transformation...
271 MolAToMolBTransformation = MolAToMolBTransformations[Index]
272 MolA = MolAToMolBTransformation.stateA.components["ligand"]
273 MolB = MolAToMolBTransformation.stateB.components["ligand"]
274 MolAName = MolA.name
275 MolBName = MolB.name
276
277 if MolAToMolBProtocolResult is None:
278 DeltaDeltaGBinding = "NA"
279 DeltaDeltaGBindingUncertainty = "NA"
280 else:
281 # Setup binding value without the units...
282 DeltaDeltaGBinding = MolAToMolBProtocolResult.get_estimate()
283 DeltaDeltaGBinding = "%.*f" % (Precision, DeltaDeltaGBinding.m)
284
285 # Setup uncertainty value without the units...
286 DeltaDeltaGBindingUncertainty = MolAToMolBProtocolResult.get_uncertainty()
287 DeltaDeltaGBindingUncertainty = "%.*f" % (Precision, DeltaDeltaGBindingUncertainty.m)
288
289 ResultData.append([MolAName, MolBName, DeltaDeltaGBinding, DeltaDeltaGBindingUncertainty])
290
291 ResultDF = pd.DataFrame(
292 ResultData,
293 columns=["MolAName", "MolBName", "DeltaDeltaG (MolA->MolB; RBFE) (kcal/mol)", "Uncertainty (kcal/mol)"],
294 )
295 ResultDF.to_csv(ResultFilePath, sep=ResultFileParams["Delim"], lineterminator="\n", index=False)
296
297
298 def ProcessInputFiles():
299 """Process input files."""
300
301 # Read PDB file...
302 MiscUtil.PrintInfo("\nReading PDB file %s..." % OptionsInfo["Infile"])
303 MacroMol = OpenFEUtil.ReadPDBFile(OptionsInfo["InfilePath"], Name=OptionsInfo["InfileRoot"])
304
305 # Read small molecule input file...
306 MiscUtil.PrintInfo("\nReading small molecule file %s..." % OptionsInfo["SmallMolFile"])
307 Mols, MolCount, ValidMolCount = OpenFEUtil.ReadAndValidateMolecules(
308 OptionsInfo["SmallMolFilePath"], **OptionsInfo["SmallMolFileParams"]
309 )
310
311 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
312 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
313 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
314
315 if ValidMolCount == 0:
316 MiscUtil.PrintInfo("")
317 MiscUtil.PrintError("No valid molecules found in small molecule input file.\n")
318
319 if ValidMolCount < 2:
320 MiscUtil.PrintInfo("")
321 MiscUtil.PrintError("Small molecule Input file must contain at least 2 molecules.\n")
322
323 return (MacroMol, Mols)
324
325
326 def ProcessMoleculeNames(Mols):
327 """Process molecule names."""
328
329 MolPairs = OpenFEUtil.ProcessMoleculePairs(Mols, OptionsInfo["MoleculePairsList"])
330
331 return MolPairs
332
333
334 def CheckMissingPartialCharges(Mols):
335 """Check missing partial charges for small molecules."""
336
337 MiscUtil.PrintInfo("\nChecking missing partial charges for specidfed small molecule pairs...")
338
339 MissingChargesMolCount = OpenFEUtil.GetMissingPartialChargesMolCount(Mols)
340 MiscUtil.PrintInfo("Number of molecules with missing partial charges: %s" % MissingChargesMolCount)
341
342 if MissingChargesMolCount == 0:
343 return
344
345 if re.match("^Stop$", OptionsInfo["MissingChargeMode"], re.I):
346 MiscUtil.PrintInfo("")
347 MiscUtil.PrintError(
348 'The small molecule input file contains molecules with missing partial charges. The execution of the script has been terminated for "Stop" value of "--missingChargedMode" option. You may continue the execution of the script by specifying "Calculate" value for "--missingChargedMode" option.\n\nThe missing charges will be automatically calculated by OpenFE RelativeHybridTopologyProtocol module during the calculation of RBFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--rbfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate RBFE.\n'
349 )
350 else:
351 MiscUtil.PrintInfo("")
352 MiscUtil.PrintWarning(
353 'The small molecule input file contains molecules with missing partial charges. The missing charges will be automatically calculated by OpenFE RelativeHybridTopologyProtocol module during the calculation of RBFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--rbfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate RBFE.\n'
354 )
355
356
357 def ProcessOutfilePrefixOption():
358 """Process outfile prefix option."""
359
360 OutfilePrefix = Options["--outfilePrefix"]
361
362 if re.match("^auto$", OutfilePrefix, re.I):
363 OutfilePrefix = OptionsInfo["SmallMolFileRoot"]
364
365 OptionsInfo["OutfilePrefix"] = OutfilePrefix
366
367
368 def ProcessOutfileDirOption():
369 """Process outfile directory Option."""
370
371 # Setup output directory...
372 OutfileDir = Options["--outfileDir"]
373 OutfileDirPath = os.path.abspath(OutfileDir)
374 if not os.path.exists(OutfileDir):
375 MiscUtil.PrintInfo("\nCreating output directory %s..." % (OutfileDir))
376 os.mkdir(OutfileDirPath)
377 OptionsInfo["OutfileDir"] = OutfileDir
378 OptionsInfo["OutfileDirPath"] = OutfileDirPath
379
380 # Setup a transformations subdirectory...
381 TransformationsOutDir = "Transformations"
382 TransformationsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], TransformationsOutDir)
383 if not os.path.exists(TransformationsOutDirPath):
384 os.mkdir(TransformationsOutDirPath)
385 OptionsInfo["TransformationsOutDir"] = TransformationsOutDir
386 OptionsInfo["TransformationsOutDirPath"] = TransformationsOutDirPath
387
388 # Setup a results subdirectory...
389 ResultsOutDir = "Results"
390 ResultsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], ResultsOutDir)
391 if not os.path.exists(ResultsOutDirPath):
392 os.mkdir(ResultsOutDirPath)
393 OptionsInfo["ResultsOutDir"] = ResultsOutDir
394 OptionsInfo["ResultsOutDirPath"] = ResultsOutDirPath
395
396 # Use results subdirectory for scratch results...
397 OptionsInfo["ResultsScratchOutDir"] = ResultsOutDir
398 OptionsInfo["ResultsScratchOutDirPath"] = ResultsOutDirPath
399
400
401 def ProcessListOption():
402 """Process list protocol settings option."""
403
404 RBFESettings = SepTopProtocol.default_settings()
405
406 MiscUtil.PrintInfo("\nListing RBFE seperated topology settings...")
407 OpenFEUtil.ListOpenFESettings(RBFESettings)
408
409
410 def ConfigureLogging():
411 """Configure logging."""
412
413 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
414
415 if re.match("^Error$", OptionsInfo["LoggingLevel"], re.I):
416 LoggingLevel = logging.ERROR
417 elif re.match("^Warning$", OptionsInfo["LoggingLevel"], re.I):
418 LoggingLevel = logging.WARNING
419 else:
420 LoggingLevel = logging.INFO
421
422 logging.basicConfig(format="%(levelname)s: %(message)s", level=LoggingLevel)
423
424 # Turn warnings issued by warnings.warn() into log message to avoid display
425 # of a stack trace...
426 logging.captureWarnings(True)
427
428
429 def ProcessOptions():
430 """Process and validate command line arguments and options."""
431
432 MiscUtil.PrintInfo("Processing options...")
433
434 # Validate options...
435 ValidateOptions()
436
437 # Configure logging...
438 ConfigureLogging()
439
440 OptionsInfo["Infile"] = Options["--infile"]
441 OptionsInfo["InfilePath"] = os.path.abspath(OptionsInfo["Infile"])
442 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
443 OptionsInfo["InfileRoot"] = FileName
444
445 OptionsInfo["SmallMolFile"] = Options["--smallMolFile"]
446 OptionsInfo["SmallMolFilePath"] = os.path.abspath(OptionsInfo["SmallMolFile"])
447 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["SmallMolFile"])
448 OptionsInfo["SmallMolFileRoot"] = FileName
449
450 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
451 OptionsInfo["SmallMolFileParams"] = MiscUtil.ProcessOptionInfileParameters(
452 "--smallMolFileParams",
453 Options["--smallMolFileParams"],
454 InfileName=Options["--smallMolFile"],
455 ParamsDefaultInfo=ParamsDefaultInfoOverride,
456 )
457
458 OptionsInfo["ExecuteDAGParams"] = OpenFEUtil.ProcessOptionOpenFEExecuteDAGParameters(
459 "--executeDAGParams", Options["--executeDAGParams"]
460 )
461
462 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
463
464 OptionsInfo["MissingChargeMode"] = OpenFEUtil.ProcessOptionOpenFEMissingChargeMode(
465 "--missingChargeMode", Options["--missingChargeMode"]
466 )
467
468 OptionsInfo["MoleculePairs"] = Options["--moleculePairs"]
469 OptionsInfo["MoleculePairsList"] = OpenFEUtil.ProcessOptionOpenFEMoleculePairs(
470 "--moleculePairs", Options["--moleculePairs"]
471 )
472 OptionsInfo["MoleculePairsMolList"] = None
473
474 OptionsInfo["ResultFileParams"] = OpenFEUtil.ProcessOptionOpenFEResultFileParameters(
475 "--resultFileParams", Options["--resultFileParams"]
476 )
477
478 ParamsDefaultInfoOverride = {"EngineComputePlatform": "CPU"}
479 OptionsInfo["RBFEParams"] = OpenFEUtil.ProcessOptionOpenFERelativeFreeEnergySeparatedTopologyParameters(
480 "--rbfeParams", Options["--rbfeParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
481 )
482
483 OptionsInfo["SolventParams"] = OpenFEUtil.ProcessOptionOpenFESolventParameters(
484 "--solventParams", Options["--solventParams"]
485 )
486
487 ProcessOutfilePrefixOption()
488 ProcessOutfileDirOption()
489
490 OptionsInfo["Overwrite"] = Options["--overwrite"]
491
492 # Track top level working directory...
493 OptionsInfo["TopWorkingDir"] = os.getcwd()
494
495
496 def RetrieveOptions():
497 """Retrieve command line arguments and options."""
498
499 # Get options...
500 global Options
501 Options = docopt(_docoptUsage_)
502
503 # Set current working directory to the specified directory...
504 WorkingDir = Options["--workingdir"]
505 if WorkingDir:
506 os.chdir(WorkingDir)
507
508 # Handle examples option...
509 if "--examples" in Options and Options["--examples"]:
510 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
511 sys.exit(0)
512
513
514 def ValidateOptions():
515 """Validate option values."""
516
517 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
518 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
519
520 MiscUtil.ValidateOptionFilePath("-s, --smallMolFile", Options["--smallMolFile"])
521 MiscUtil.ValidateOptionFileExt("-s, --smallMolFile", Options["--smallMolFile"], "sdf sd")
522
523 MiscUtil.ValidateOptionDirPath("-o, --outfileDir", Options["--outfileDir"])
524 MiscUtil.ValidateOptionsOutputDirOverwrite(
525 "-o, --outfileDir", Options["--outfileDir"], "--overwrite", Options["--overwrite"]
526 )
527
528 MiscUtil.ValidateOptionTextValue("--loggingLevel", Options["--loggingLevel"], "Info Warning Error")
529
530 MiscUtil.ValidateOptionTextValue("--missingChargeMode", Options["--missingChargeMode"], "Calculate Stop")
531
532 MoleculePairs = Options["--moleculePairs"]
533 if not re.match("^auto$", MoleculePairs, re.I):
534 MoleculePairsList = MoleculePairs.split(",")
535 if len(MoleculePairsList) % 2:
536 MiscUtil.PrintError(
537 'The number of comma delimited values, %d, specified using "--moleculePairs" option must be an even number.'
538 % (len(MoleculePairsList))
539 )
540
541
542 # Setup a usage string for docopt...
543 _docoptUsage_ = """
544 OpenFECalculateRelativeBindingFreeEnergySepTop.py - Calculate relative binding free energy
545
546 Usage:
547 OpenFECalculateRelativeBindingFreeEnergySepTop.py [--executeDAGParams <Name,Value,..>] [--loggingLevel <Info, Warning or Error>]
548 [--missingChargeMode <Calculate or Stop>] [--moleculePairs <MolName1,MolName2,..>]
549 [--outfilePrefix <text>] [--overwrite] [--rbfeParams <Name,Value,...>] [--resultFileParams <Name,Value,..>]
550 [--smallMolFileParams <Name,Value,...> ] [--solventParams <Name,Value,...>]
551 [-w <dir>] -i <infile> -s <smallmolfile> -o <outifiledir>
552 OpenFECalculateRelativeBindingFreeEnergySepTop.py -l | --list
553 OpenFECalculateRelativeBindingFreeEnergySepTop.py -h | --help | -e | --examples
554
555 Description:
556 Calculate Relative Binding Free Energy (RBFE) for a pair of molecules in a
557 small molecule input file using a Separated Topologies (SepTop) [ Ref 183 ]
558 approach.
559
560 You may calculate RBFEs for any arbitrary pairs of molecules. The SepTop
561 protocol calculates RBFE by performing two ABFE calculations simultaneously
562 in opposite directions. No atom mapping is required between pairs of molecules
563 to create hybrid topology. Consequently, the SepTop based RBFEs calculations
564 may be used for scaffold hopping.
565
566 A brief description of the SepTop protocol, taken from OpenFE documenation is
567 provided below:
568
569 The relative binding free energy is obtained through a thermodynamic
570 cycle. The molecules are transformed into each other both in solvent,
571 giving ddG (solvent), and in the complex, giving ddG (complex), which
572 allows the calculation of the relative binding free energy, ddG. Each
573 molecule is represented with its own set of coordinates, meaning that
574 the interactions of all atoms of one molecule are turned off while
575 simultaneously turning on the interactions of all atoms of the other
576 molecule. Therefore, restraints are required: Molecules are restrained
577 to the protein in the complex states using orientational (Boresch-style)
578 restraints; in the solvent states molecules are restrained to remain
579 apart from each other using a single harmonic distance restraint between
580 the molecules.
581
582 The coulombic interactions of the molecule are fully turned off
583 (annihilated), while the Lennard-Jones interactions are decoupled,
584 meaning the intermolecular interactions are turned off, while keeping
585 the intramolecular Lennard-Jones interactions.
586
587 The input file must contain a macromolecule already prepared for simulation.
588 The preparation of the macromolecule for a simulation generally involves the
589 following tasks: identification and replacement of non-standard residues;
590 addition of missing residues; addition of missing heavy atoms; addition of
591 missing hydrogens.
592
593 In addition, the small molecule input file must contain molecules already
594 prepared for simulation. It must contain appropriate 3D coordinates relative
595 to the macromolecule along with no missing hydrogens.
596
597 Protocol repeats, 3
598
599 Time step size: 4.0 femtosecond
600
601 Complex equilibration phase:
602
603 Max minimization steps: 5,000
604 NVT equilibration length: 0.1 nanosecond
605 NPT equilibration length: 0.1 nanosecond
606 NPT length: 2.0 nanosecond
607
608 Complex production phase:
609
610 Max minimization steps: 5,000
611 NPT equilibration length: 1.0 nanosecond
612 NPT length: 10.0 nanosecond
613
614 Solvent equilibration phase:
615
616 Max minimization steps: 5,000
617 NVT equilibration length: 0.1 nanosecond
618 NPT equilibration length: 0.1 nanosecond
619 NPT length: 2.0 nanosecond
620
621 Solvent production phase:
622
623 Max minimization steps: 5,000
624 NPT equilibration length: 1.0 nanosecond
625 NPT length: nanosecond 10.0 nanosecond
626
627 Each complex and solvent simulation, by default, may run for 13.2 and 13.2
628 nanosecond respectively, for a total of 26.4 nanoseconds. The total MD
629 simulation time for correspond to 79.2 nanosecond to repeat the protocol
630 3 times for the complex and solvent simulations.
631
632 The supported macromolecule input file formats are: PDB (.pdb) and
633 CIF (.cif)
634
635 The supported small molecule input file format are : SD (.sdf, .sd)
636
637 Possible outfile prefix:
638
639 <OutfilePrefix> or <SmallMolFileRoot>
640
641 Possible output directories:
642
643 <OutfileDir>
644
645 <OutfileDir>/Transformations
646 <OutfileDir>/Results
647
648 Possible output files and directories under <OutfileDir>:
649
650 <OutfilePrefix>_RBFE_Results.<csv or tsv>
651
652 ... ... ...
653
654 Transformations/<MolAName>_To_<MolBName>_Complex_Solvent.json
655 ... ... ...
656
657 Results/<MolAName>_To_<MolBName>_Complex_Solvent_Results.json
658 Results/shared_SepTopComplexRunUnit-*/
659 Results/shared_SepTopComplexSetupUnit-*/
660 Results/shared_SepTopSolventRunUnit-*/
661 Results/shared_SepTopSolventSetupUnit-*/
662 ... ... ...
663
664 Options:
665 -e, --examples
666 Print examples.
667 --executeDAGParams <Name,Value,..> [default: auto]
668 A comma delimited list of parameter name and value pairs for executing
669 protocol DAGs (Directed Acyclic Graph) to run RBFE calculations.
670
671 The supported parameter names along with their default values are
672 are shown below:
673
674 keepShared, yes [ Possible values: yes or no ]
675 keepScratch, no [ Possible values: yes or no ]
676 nRetries, 2 [ Possible values: >= 0. A value of 0 implies only
677 1 try. ]
678
679 A brief description of parameters is provided below:
680
681 keepShared: Keep shared directories after the execution of DAG.
682 keepScratch: Keep scratch directories after the execution of DAG.
683 nRetries: Number of times to attempt the execution.
684
685 -h, --help
686 Print this help message.
687 -i, --infile <infile>
688 Input file name containing a macromolecule.
689 -l, --list
690 List default RBFE protocol settings provided by OpenFE module
691 RelativeHybridTopologyProtocol.
692 --loggingLevel <Info, Warning or Error> [default: Error]
693 Logging level to configure the 'root logger' via logging.basicConfig()
694 function. The default logging level is changed from 'logging.INFO' to
695 'logging.ERROR'. Otherwise, OpenFE and its associated modules
696 may generate a lot of informational messages.
697 --missingChargeMode <Calculate or Stop> [default: Stop]
698 Calculate missing partial charges for molecules before running RBFE
699 calculations or terminate the execution of the script. The missing
700 partial charges will be automatically calculated by OpenFE module
701 RelativeHybridTopologyProtocol during the calculation of RBFE. You
702 may control the calculation of partial charges by specifying values for
703 partialCharge* parameters using '--rbfeParams' option.
704 -m, --moleculePairs <MolName1,MolName2,..> [default: auto]
705 A comma delimited list of molecule name pairs for calculating RBFEs.
706 Default: the names of the first and second molecule in small molecule
707 input file.
708 -o, --outfileDir <outfiledir>
709 Output directory.
710 --outfilePrefix <text> [default: auto]
711 Prefix for generating output files under output directory.
712 --overwrite
713 Overwrite existing files.
714 --resultFileParams <Name,Value,..> [default: auto]
715 A comma delimited list of parameter name and value pairs for writing
716 calculated RBFEs values to a results file.
717
718 The supported parameter names along with their default values are
719 are shown below:
720
721 precision, 4 [ Possible values: > 0 ]
722 delimiter, comma [ Possible values: comma or tab ]
723
724 -r, --rbfeParams <Name,Value,...> [default: auto]
725 A comma delimited list of parameter name and value pairs for RBFE separated
726 topology protocol settings employed during the calculation of RBFEs.
727
728 The default values are automatically updated to match settings provided by
729 OpenFE module SepTopProtocol.
730
731 You must specify valid OpenFE values for these parameters. An extensive
732 validation is not performed.
733
734 The supported parameter names along with their default values are
735 are shown below:
736
737 protocolRepeats, 3
738
739 Complex equil output settings:
740
741 complexEquilOutputCheckpointInterval, 1 [ Units: nanosecond ]
742 complexEquilOutputCheckpointStorageFilename, checkpoint.chk
743 complexEquilOutputEquilNPTStructure, equil_npt.pdb
744 complexEquilOutputEquilNVTstructure, None
745 complexEquilOutputForcefieldCache, db.json
746 complexEquilOutputLogOutput, equil_simulation.log
747 complexEquilOutputMinimizedStructure, minimized.pdb
748 complexEquilOutputIndices, all [ Possible value: Any valid
749 selection. ]
750 complexEquilOutputPreminimizedStructure, system.pdb
751 complexEquilOutputProductionTrajectoryFilename, production_equil.xtc
752 complexEquilOutputTrajectoryWriteInterval, 20.0 [ Units:
753 picosecond ]
754
755 Complex equil simulation settings:
756
757 complexEquilSimulationEquilibrationLength, 0.1 [ Units: nanosecond ]
758 complexEquilSimulationEquilibrationLengthNVT, 0.1 [ Units:
759 nanosecond ]
760 complexEquilSimulationMinimizationSteps, 5000
761 complexEquilSimulationProductionLength, 2.0 [ Units: nanosecond ]
762
763 Complex lambda settings:
764
765 complexLambdaElecA, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.5 0.75
766 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 [ Possible values: A space
767 delimited list of values between 0.0 and 1.0 ]
768 complexLambdaElecB, 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.75 0.5 0.25
769 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 [ Possible values: A space
770 delimited list of values between 0.0 and 1.0 ]
771 complexLambdaRestraintsA, 0.0 0.05 0.1 0.3 0.5 0.75 1.0 1.0 1.0
772 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 [ Possible values: A
773 space delimited list of values between 0.0 and 1.0 ]
774 complexLambdaRestraintsB, 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
775 1.0 1.0 1.0 0.75 0.5 0.3 0.1 0.05 0.0 [ Possible values: A space
776 delimited list of values between 0.0 and 1.0 ]
777 complexLambdaVdwA, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
778 0.143 0.286 0.429 0.572 0.715 0.857 1.0 [ Possible values: A
779 delimited list of values between 0.0 and 1.0 ]
780 complexLambdaVdwB, 1.0 0.857 0.715 0.572 0.429 0.286 0.143 0.0
781 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 [ Possible values: A
782 delimited list of values between 0.0 and 1.0 ]
783
784 Complex output settings:
785
786 complexOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
787 complexOutputCheckpointStorageFilename, complex_checkpoint.nc
788 complexOutputForcefieldCache, db.json
789 complexOutputFilename, complex.nc
790 complexOutputIndices, not water [ Possible value: Any valid
791 selection. ]
792 complexOutputStructure, alchemical_system.pdb
793 complexOutputPositionsWriteFrequency, 100.0 [ Units: picosecond ]
794 complexOutputVelocitiesWriteFrequency, None [ Possible
795 values: > 0; Units: picosecond ]
796
797 Complex restraint settings:
798
799 complexRestraintKPhiA, 334.72 [ Units: kilojoule_per_mole/radian**2
800 The default value is equivalent to 80 kcal/mol/radian**2 ]
801 complexRestraintKPhiB, 334.72 [ Units: kilojoule_per_mole/radian**2
802 The default value is equivalent to 80 kcal/mol/radian**2 ]
803 complexRestraintKPhiC, 334.72 [ Units: kilojoule_per_mole/radian**2
804 The default value is equivalent to 80 kcal/mol/radian**2 ]
805 complexRestraintKR, 4184.0 [ Units: kilojoule_per_mole/nanometer**2
806 The default value is equivalent to 10 kcal/mol/angstrom**2
807 complexRestraintKThetaA, 334.72 [ Units:kilojoule_per_mole/radian**2
808 The default value is equivalent to 80 kcal/mol/radian**2 ]
809 complexRestraintKThetaB, 334.72 [ Units:kilojoule_per_mole/radian**2
810 The default value is equivalent to 80 kcal/mol/radian**2 ]
811 complexRestraintAnchorFindingStrategy, bonded [ Possible values:
812 multi-residue or bonded ]
813 complexRestraintDsspFilter, yes [ Possible values: yes or no ]
814 complexRestraintHostMaxDistance, 1.5 [ Units: nanometer ]
815 complexRestraintHostMinDistance, 0.5 [ Units: nanometer ]
816 complexRestraintHostSelection, backbone [ Possible value: Any valid
817 selection. ]
818 complexRestraintRmsfCutoff, 0.1 [ Units: nanometer ]
819
820 Complex simulation settings:
821
822 complexSimulationEarlyTerminationTargetError, 0.0 [ Units:
823 kilocalorie_per_mole ]
824 complexSimulationEquilibrationLength, 1.0 [ Units: nanosecond ]
825 complexSimulationMinimizationSteps, 5000
826 complexSimulationNReplicas, 19
827 complexSimulationProductionLength, 10.0 [ Units: nanosecond ]
828 complexSimulationRealTimeAnalysisInterval, 250.0 [ Units:
829 picosecond ]
830 complexSimulationRealTimeAnalysisMinimumTime, 500.0 [ Units:
831 picosecond ]
832 complexSimulationSamplerMethod, repex [ Possible values: repex,
833 sams, or independent ]
834 complexSimulationSamsFlatnessCriteria, logZ-flatness [ Possible
835 values: logZ-flatness, minimum-visits or histogram-flatness ]
836 complexSimulationSamsGamma0, 1.0
837 complexSimulationTimePerIteration, 2.5 [ Units: picosecond ]
838
839 Complex solvation settings:
840
841 complexSolvationBoxShape, dodecahedron [ Possible values: cube,
842 dodecahedron, or octahedron ]
843 complexSolvationBoxSize, None [ Possible value: A triplet of space
844 X Y Z values; Units: nanometer ]
845 complexSolvationSolventModel, tip3p [ Possible values: tip3p, spce,
846 tip4pew, or tip5p ]
847 complexSolvationSolventPadding, 1.0 [ Units: nanometer ]
848
849 Engine settings:
850
851 engineComputePlatform, CPU [ Possible values: CPU, CUDA,
852 OpenCL, or Reference ]
853 engineGpuDeviceIndex, None [ Possible values: 0, 0 1, etc. ]
854
855 Forcefield settings:
856
857 forcefieldConstraints, HBonds [ Possible values: HBonds,
858 AllBonds, or HAngles ]
859 forcefields, amber/ff14SB.xml amber/tip3p_standard.xml
860 amber/tip3p_HFE_multivalent.xml amber/phosaa10.xml
861 [ Possible values: A space delimited list of valid names. ]
862 forcefieldHydrogenMass, 3.0 [ Units: amu ]
863 forcefieldNonbondedCutoff, 0.9 [ Units: nanometer ]
864 forcefieldNonbondedMethod, PME [ Possible values: PME or
865 NoCutoff ]
866 forcefieldRigidWater, yes [ Possible values: yes or no ]
867 forcefieldSmallMoleculeForcefield, openff-2.1.1 [ Possible
868 value: A valid forcefield name. ]
869
870 Integrator settings:
871
872 integratorBarostatFrequency, 25.0 * timestep [ The specified value
873 is a multiple of integratorTimestep. ]
874 integratorConstraintTolerance, 1e-06
875 integratorLangevinCollisionRate, 1.0 [ Units: 1 / picosecond ]
876 integratorNRestartAttempts, 20
877 integratorReassignVelocities, no [ Possible values: yes or no ]
878 integratorRemoveCom, no [ Possible values: yes or no ]
879 integratorTimestep, 4.0 [ Units: femtosecond ]
880
881 Partial charge settings:
882
883 partialChargeNaglModel, None [ Default: Production AM1BCC model for
884 NAGL; Possible value: Any valid name. ]
885 partialChargeNumberOfConformers, None [ Possible value: > 0 ]
886 partialChargeOffToolkitBackend, AmberTools [ Possible values:
887 AmberTools or RDKit ]
888 partialChargeMethod, AM1BCC [ Possble values: AM1BCC, Espaloma,
889 or NAGL ]
890
891 Solvent equil output settings:
892
893 solventEquilOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
894 solventEquilOutputCheckpointStorageFilename, checkpoint.chk
895 solventEquilOutputEquilNPTStructure, equil_npt.pdb
896 solventEquilOutputEquilNVTstructure, None
897 solventEquilOutputForcefieldCache, db.json
898 solventEquilOutputLogOutput, equil_simulation.log
899 solventEquilOutputMinimizedStructure, minimized.pdb
900 solventEquilOutputIndices, all [ Possible value: Any valid
901 selection. ]
902 solventEquilOutputPreminimizedStructure, system.pdb
903 solventEquilOutputProductionTrajectoryFilename, equil_npt.xtc
904 solventEquilOutputTrajectoryWriteInterval, 20.0 [ Units:
905 picosecond ]
906
907 Solvent_equil_simulation_settings:
908
909 solventEquilSimulationEquilibrationLength, 0.1 [ Units: nanosecond ]
910 solventEquilSimulationEquilibrationLengthNVT, 0.1 [ Units:
911 nanosecond ]
912 solventEquilSimulationMinimizationSteps, 5000
913 solventEquilSimulationProductionLength, 2.0 [ Units: nanosecond ]
914
915 Solvent lambda settings:
916
917 solventLambdaElecA, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125
918 0.25 0.375 0.5 0.625 0.75 0.875 1.0 1.0 1.0 1.0 1.0 1.0 1.0
919 1.0 1.0 1.0 [ Possible values: A space delimited list of values
920 between 0.0 and 1.0 ]
921 solventLambdaElecB, 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 0.875
922 0.75 0.625 0.5 0.375 0.25 0.125 0.0 0.0 0.0 0.0 0.0 0.0 0.0
923 0.0 0.0 0.0 [ Possible values: A space delimited list of values
924 between 0.0 and 1.0 ]
925 solventLambdaRestraintsA, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
926 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
927 0.0 [ Possible values: A space delimited list of values between
928 0.0 and 1.0 ]
929 solventLambdaRestraintsB, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
930 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
931 0.0 [ Possible values: A space delimited list of values between
932 0.0 and 1.0 ]
933 solventLambdaVdwA, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
934 0.0 0.0 0.0 0.0 0.0 0.0 0.15 0.23 0.3 0.4 0.52 0.64 0.76 0.88
935 1.0 [ Possible values: A space delimited list of values between
936 0.0 and 1.0 ]
937 solventLambdaVdwB, 1.0 0.85 0.77 0.7 0.6 0.48 0.36 0.24 0.12 0.0
938 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
939 0.0 [ Possible values: A space delimited list of values between
940 0.0 and 1.0 ]
941
942 Solvent output settings:
943
944 solventOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
945 solventOutputCheckpointStorageFilename, solvent_checkpoint.nc
946 solventOutputForcefieldCache, db.json
947 solventOutputFilename, solvent.nc
948 solventOutputIndices, not water [ Possible value: Any valid
949 selection. ]
950 solventOutputStructure, alchemical_system.pdb
951 solventOutputPositionsWriteFrequency, 100.0 [ Units: picosecond ]
952 solventOutputVelocitiesWriteFrequency, None [ Possible
953 values: > 0; Units: picosecond ]
954
955 Solvent restraint settings:
956
957 solventRestraintCentralAtomsOnly, No [ Possible values: yes or no ]
958 solventRestraintSpringConstant, 1000.0 [ Units: kilojoule_per_mole /
959 nanometer ** 2. The default value is equivalent to 2.40
960 kilocalorie_per_mole / angstromg ** 2 ]
961
962 Solvent simulation settings:
963
964 solventSimulationEarlyTerminationTargetError, 0.0 [ Units:
965 kilocalorie_per_mole ]
966 solventSimulationEquilibrationLength, 1.0 [ Units: nanosecond ]
967 solventSimulationMinimizationSteps, 5000
968 solventSimulationNReplicas, 27
969 solventSimulationProductionLength, 10.0 [ Units: nanosecond ]
970 solventSimulationRealTimeAnalysisInterval, 250.0 [ Unit: picosecond ]
971 solventSimulationRealTimeAnalysisMinimumTime, 500.0 [ Units:
972 picosecond ]
973 solventSimulationSamplerMethod, repex [ Possible values: repex,
974 sams, or independent ]
975 solventSimulationSamsFlatnessCriteria, logZ-flatness [ Possible
976 values: logZ-flatness, minimum-visits or histogram-flatness ]
977 solventSimulationSamsGamma0, 1.0
978 solventSimulationTimePerIteration, 2.5 [ Units: picosecond ]
979
980 Solvent solvation settings:
981
982 solventSolvationBoxShape, dodecahedron [ Possible values: cube,
983 dodecahedron, or octahedron ]
984 solventSolvationBoxSize, None [ Possible value: A triplet of space
985 X Y Z values; Units: nanometer ]
986 solventSolvationSolventModel, tip3p [ Possible values: tip3p, spce,
987 tip4pew, or tip5p ]
988 solventSolvationSolventPadding, 1.5 [ Units: nanometer ]
989
990 Thermo settings:
991
992 thermoPh, None [ Possible values: > 0 ]
993 thermoPressure, 1.0 [ Units: bar ]
994 thermoRedoxPotential, None [ Possible values: A valid float.
995 Units: millivolts (mV) ]
996 thermoTemperature, 298.15 [ Units: kelvin ]
997
998 A brief description of parameters, taken from OpenFE documentation, is
999 provided below:
1000
1001 protocolRepeats: Number of completely independent repeats of the
1002 entire sampling process.
1003
1004 Complex settings:
1005
1006 Complex parameters for the system, including the solvent model and
1007 the solvent padding.
1008
1009 Complex equil output settings:
1010
1011 Parameters controlling simulation output during equilibration
1012 phase of complex transformation.
1013
1014 complexEquilOutputCheckpointInterval: Frequency to write the
1015 checkpoint file.
1016 complexEquilOutputCheckpointStorageFilename: Checkpoint filename.
1017 complexEquilOutputEquilNPTStructure: NPT structure filename.
1018 complexEquilOutputEquilNVTstructure: NVT strucure filename.
1019 complexEquilOutputForcefieldCache: Filename for caching small
1020 molecule residue templates.
1021 complexEquilOutputLogOutput: Simulation log filename.
1022 complexEquilOutputMinimizedStructure: Minimized structure filename.
1023 complexEquilOutputIndices: Selection string for selecting
1024 coordinates to write.
1025 complexEquilOutputPremnimizedStructure: Initial structure filename.
1026 complexEquilOutputProductionTrajectoryFilename: Trajectory filename.
1027 complexEquilOutputTrajectoryWriteInterval: Frequency for writing
1028 velocities to trajectory file.
1029
1030 Complex equil simulation settings:
1031
1032 Parameters controlling simulation during equilibration phase of
1033 complex transformation.
1034
1035 complexEquilSimulationEquilibrationLength: Length of the NPT
1036 equilibration phase.
1037 complexEquilSimulationEquilibrationLengthNVT: Length of the NVT
1038 equilibration phase.
1039 complexEquilSimulationMinimizationSteps: Maximum number of
1040 minimization steps to perform.
1041 complexEquilSimulationProductionLength: Length of the NPT
1042 production phase.
1043
1044 Complex lambda settings:
1045
1046 Lambda protocol parameters for complex transformation.
1047
1048 complexLambdaElecA: List of lambda values for electrostatics. The
1049 values of 0 and 1 imply state A and state B respectively.
1050 complexLambdaElecB: List of lambda values for electrostatics. The
1051 values of 0 and 1 imply state A and state B respectively.
1052 complexLambdaRestraintsA: List of lambda values for restraints. The
1053 values of 0 and 1 imply state A and state B respectively.
1054 complexLambdaRestraintsB: List of lambda values for restraints. The
1055 values of 0 and 1 imply state A and state B respectively.
1056 complexLambdaVdwA: List of lamda values for van der Waals. The
1057 values of of 0 and 1 imply state A and state B respectively.
1058 complexLambdaVdwB: List of lamda values for van der Waals. The
1059 values of of 0 and 1 imply state A and state B respectively.
1060
1061 Complex output settings:
1062
1063 Parameters controlling simulation output during final phase of
1064 complex transformation.
1065
1066 complexOutputCheckpointInterval: Frequency to write the checkpoint
1067 file.
1068 complexOutputCheckpointStorageFilename: Checkpoint filename.
1069 complexOutputForcefieldCache: Filename for caching small molecule
1070 residue templates.
1071 complexOutputFilename: Trajectory filename.
1072 complexOutputIndices: Selection string for selecting coordinates to
1073 write.
1074 complexOutputStructure: Topology structure filename.
1075 complexOutputPositionsWriteFrequency: Frequency for writing
1076 positions to trajectory file.
1077 complexOutputVelocitiesWriteFrequency: Frequency for writing
1078 velocities to trajectory file.
1079
1080 Complex restraint settings:
1081
1082 Parameters to configure Boresch-style restraint between two groups
1083 of atoms named host (Hx) and guest (Gx).
1084
1085 complexRestraintKPhiA: Equilibrium force constant for the dihedral
1086 formed by H2-H1-H0-G0.
1087 complexRestraintKPhiB: Equilibrium force constant for the dihedral
1088 formed by H1-H0-G0-G1.
1089 complexRestraintKPhiC: Equilibrium force constant for the dihedral
1090 formed by H0-G0-G1-G2.
1091 complexRestraintKR: Bond spring constant between H0 and G0.
1092 restraintKThetaA: Spring constant for the angle formed by H1-H0-G0.
1093 complexRestraintKThetaA: Spring constant for the angle formed by
1094 H1-H0-G0.
1095 complexRestraintKThetaB: Spring constant for the angle formed by
1096 H0-G0-G1.
1097 complexRestraintAnchorFindingStrategy: Boresch atom picking strategy
1098 to use. bonded: pick host atoms that are bonded to each other.
1099 multi-residue: pick host atoms which can span multiple residues.
1100 complexRestraintDsspFilter: Apply DSSP filter to the host atoms.
1101 complexRestraintHostMaxDistance: Maximum distance between any
1102 host atom and the guest G0 atom.
1103 complexRestraintHostMinDistance: Minimum distance between any
1104 host atom and the guest G0 atom
1105 complexRestraintHostSelection: A valid selection string to
1106 sub-select the host atoms which will be involved in the
1107 restraint.
1108 complexRestraintRmsfCutoff: Cutoff value for filtering atoms by their
1109 root mean square fluctuation. Atoms with values above this
1110 cutoff are ignored.
1111
1112 Complex simulation settings:
1113
1114 Parameters controlling simulation during final phase of complex
1115 transformation.
1116
1117 complexSimulationEarlyTerminationTargetError: Target error for the
1118 real time analysis measured in kcal/mol. Once the MBAR error of
1119 the free energy is at or below this value, the simulation will
1120 be considered complete. The suggested value of 0.12 has shown to
1121 be effective in both hydration and binding free energy
1122 benchmarks.
1123 complexSimulationEquilibrationLength: Length of the equilibration
1124 phase. The specified value must be divisible by
1125 'integratorTimestep'.
1126 complexSimulationMinimizationSteps: Maximum number of minimization
1127 steps to perform.
1128 complexSimulationNReplicas: Number of replicas to use.
1129 complexSimulationProductionLength: Length of the production phase.
1130 The specified value must be divisible by 'integratorTimestep'.
1131 complexSimulationRealTimeAnalysisMinimumTime: Time interval for
1132 performing analysis of the free energies. At each interval, real
1133 time analysis data will be written to a yaml file named
1134 <outputFileName>_real_time_analysis.yaml. The current error
1135 in the estimate will also be assessed and the simulation will
1136 be terminated when it drops below
1137 'complexSimulationEarlyTerminationTargetError'.
1138 complexSimulationSamplerMethod: Alchemical sampling method to use:
1139 REPEX (Hamiltonian REPlica EXchange), SAMS (Self-Adjusted
1140 Mixture Sampling), or Independent (Independently sampled lambda
1141 windows).
1142 complexSimulationSamsFlatnessCriteria:Method for assessing when to
1143 switch to asymptomatically optimal scheme for SAMS.
1144 complexSimulationsamsGamma0: Initial weight adaptation rate for
1145 SAMS.
1146 complexSimulationTimePerIteration: Simulation time between each
1147 MCMC move attempt
1148
1149 Complex solvation settings:
1150
1151 Solvation parameters for the system, including the solvent model and
1152 the solvent padding.
1153
1154 complexSolvationBoxShape: Shape of the periodic solvent box.
1155 complexSolvationBoxSize: Lengths of the unit cell for a solvent box.
1156 complexSolvationSolventModel: Forcefield water model to use during
1157 solvation and defining the model properties.
1158 complexSolvationSolventPadding: Minimum distance from any solute
1159 bounding sphere to the edge of the box.
1160
1161 Engine settings:
1162
1163 Parameters configuring the compute platform used by the OpenMM to
1164 perform the simulation.
1165
1166 engineComputePlatform: Platform to use for running OpenMM MD
1167 calculations.
1168 engineGpuDeviceIndex: Space delimited list of device indices
1169 to use for running OpenMM MD calculations.
1170
1171 Forcefield settings:
1172
1173 forcefieldConstraints:Constraints to use.
1174 forcefields: List of valid forcefield paths for all components
1175 except small molecules.
1176 forcefieldHydrogenMass: Mass to be repartitioned to hydrogens
1177 from neighboring heavy atoms.
1178 forcefieldNonbondedCutoff: Cutoff for short range nonbonded
1179 interactions.
1180 forcefieldNonbondedMethod: Method for treating nonbonded
1181 interactions.
1182 forcefieldRigidWater: Use a rigid water model.
1183 forcefieldSmallMoleculeForcefield: A valid forcefield name to use
1184 for small molecules.
1185
1186 Integrator settings:
1187
1188 Parameters controlling the LangevinSplittingDynamicsMove integrator
1189 used for simulation.
1190
1191 integratorBarostatFrequency: Frequency at which volume scaling
1192 changes should be attempted.
1193 integratorConstraintTolerance: Tolerance for constraint solver.
1194 integratorLangevinCollisionRate: Collision frequency.
1195 integratorNRestartAttempts: Number of attempts to restart from
1196 Context in case there are NaNs in the energies after
1197 integration.
1198 integratorReassignVelocities: Reassign velocities from the
1199 Maxwell-Boltzmann distribution at the beginning of each
1200 Monte Carlo move.
1201 integratorRemoveCom: Remove the center of mass motion.
1202 integratorTimestep: Size of the simulation timestep.
1203
1204 Partial charge settings:
1205
1206 Parameters for automatically assigning missing partial charges to
1207 small molecules, including the partial charge method.
1208
1209 partialChargeNaglModel: Model to use for partial charge assignment.
1210 A value of None implies the use of the latest available
1211 production AM1BCC model.
1212 partialChargeNumberOfConformers: Number of conformers to generate
1213 as part of the partial charge assignment. A value of None
1214 implies the use of the existing conformer.
1215 partialChargeOffToolkitBackend: OpenFF toolkit registry backend to
1216 use for calculating partial charges.
1217 partialChargeMethod: Method to use for calculating partial charges.
1218
1219 Solvent equil output settings:
1220 Solvent equil simulation settings:
1221 Solvent lambda settings:
1222 Solvent output settings:
1223 Solvent restraint settings:
1224 Solvent simulation settings:
1225 Solvent solvation settings:
1226
1227 The solvent settings are similar to the complex settings already
1228 described under various sections for complex. The prefix 'solvent'
1229 is used for the names of the pramaters instead of the prefix
1230 'complex.'
1231
1232 Thermo settings:
1233
1234 Thermodynamic parameters, including the temperature and the pressure
1235 of the system.
1236
1237 thermoPh: Simulation pH
1238 thermoPressure: Simulation pressure.
1239 thermoRedoxPotential:Simulation redox potential.
1240 thermoTemperature: Simulation temperature.
1241
1242 -s, --smallMolFile <SmallMolFile>
1243 Input file containing small molecules.
1244 --smallMolFileParams <Name,Value,...> [default: auto]
1245 A comma delimited list of parameter name and value pairs for reading
1246 molecules from files. The supported parameter names for different file
1247 formats, along with their default values, are shown below:
1248
1249 SD: removeHydrogens,no,sanitize,yes,strictParsing,yes
1250
1251 --solventParams <Name,Value,...> [default: auto]
1252 A comma delimited list of parameter name and value pairs for solvent
1253 component. You must specify valid OpenFE values. No extensive validation
1254 is performed. These parameters are used in conjunction with solvation*
1255 parameters available through '--rbfeParams' to perform solvation.
1256
1257 The supported parameter names along with their default values are
1258 are shown below:
1259
1260 positiveIon, Na+ [ Possible value: Li+, Na+, K+, Rb+, or Cs+ ]
1261 negativeIon, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
1262 neutralize, yes [ Possible values: yes or no ]
1263 ionConcentration, 0.15 [ Units: molar ]
1264
1265 A brief description of parameters is provided below:
1266
1267 positiveIon, negativeion: Pair of ions used to neutralize and bring
1268 the solvent to required ionic concentration.
1269 neutralize: Neutralize the net charge on the chemical state by the
1270 ions in the solvent component.
1271 ionConcentration: Ionic concentration.
1272
1273 -w, --workingdir <dir>
1274 Location of working directory which defaults to the current directory.
1275
1276 Examples:
1277 The sample protein and ligand files for tyrosine kinase 2 (Tyk2) are
1278 distributed with MayaChemTools and are available in data directory. These
1279 files have been taken from OpenFE distribution for example notebooks. The
1280 AM1BCC partial charges have been calculated for the ligands in SD file to
1281 facilitate calculations. You may review OpenFE tutorial notebooks for the
1282 expected results.
1283
1284 To calculate RBFE for a pair molecules corresponding to the fist and second
1285 molecules in a SD file, performing 3 independent repeats of the entire MD
1286 sampling process to calculate RBFE for a pair of molecules, each complex
1287 and solvent MD repeat consisting of equilibration phase (Complex:
1288 Minimization - 5,000; NVT - 0.1 ns; NPT - 0.1; NPT prod - 2.0 ns; Solvent:
1289 Minimization - 5,000; NVT - 0.1; NPT - 0.1 ns; NPT prod - 2.0 ns) and
1290 production phase ( Complex: Minimization - 5,000; NPT equil - 1.0 ns; NPT prod:
1291 10.0 ns; Solvent: Minimization - 5,000; NPT equil - 1.0 ns; NPT prod - 10 ns)
1292 using a step size of of 4 fs, writing out appropriate trajectory and PDB files for
1293 each MD repeat in Results subdirectory under output directory, type:
1294
1295 % OpenFECalculateRelativeBindingFreeEnergySepTop.py -i SampleTyk2.pdb
1296 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFESepTop
1297
1298 To run the first example for calculating RBFE for a specific pair molecules
1299 using CUDA platform on your machine to perform MD simulations and generate
1300 various output files, type:
1301
1302 % OpenFECalculateRelativeBindingFreeEnergySepTop.py -i SampleTyk2.pdb
1303 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFESepTop
1304 --moleculePairs "lig_ejm_31, lig_ejm_47"
1305 --rbfeParams "engineComputePlatform,CUDA"
1306
1307 To run the second example to see all warning messages produced by OpenFE
1308 modules and write various output files, type;
1309
1310 % OpenFECalculateRelativeBindingFreeEnergySepTop.py -i SampleTyk2.pdb
1311 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFESepTop
1312 --moleculePairs "lig_ejm_31, lig_ejm_47"
1313 --rbfeParams "engineComputePlatform,CUDA"
1314 --loggingLevel Warning
1315
1316 To run the first example for calculating RBFE for a specific pair molecules
1317 using CUDA platform on your machine to perform MD simulations, automatically
1318 calculate missing partial charges for molecules, and generate various output
1319 files, type:
1320
1321 % OpenFECalculateRelativeBindingFreeEnergySepTop.py -i SampleTyk2.pdb
1322 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFESepTop
1323 --moleculePairs "lig_ejm_31, lig_ejm_47"
1324 --rbfeParams "engineComputePlatform,CUDA"
1325 --missingChargeMode Calculate
1326
1327 To run the second example by specifying explict values for various parametres
1328 and generate various output files, type:
1329
1330 % OpenFECalculateRelativeBindingFreeEnergySepTop.py -i SampleTyk2.pdb
1331 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFESepTop
1332 --moleculePairs "lig_ejm_31, lig_ejm_47"
1333 --loggingLevel Warning
1334 --executeDAGParams "keepShared, yes, nRetries, 2"
1335 --missingChargeMode Stop
1336 --rbfeParams "engineComputePlatform,CUDA,
1337 protocolRepeats,3,
1338 engineComputePlatform,CUDA,
1339 forcefieldConstraints, HBonds, forcefieldHydrogenMass, 3.0,
1340 forcefieldNonbondedMethod, PME, integratorTimestep, 4.0,
1341 complexEquilSimulationMinimizationSteps, 5000,
1342 complexEquilSimulationEquilibrationLength, 0.1,
1343 complexEquilSimulationEquilibrationLengthNVT, 0.1,
1344 complexEquilSimulationProductionLength, 2.0,
1345 complexSimulationMinimizationSteps, 5000,
1346 complexSimulationEquilibrationLength, 1.0,
1347 complexSimulationProductionLength, 10.0,
1348 solventEquilSimulationMinimizationSteps, 5000,
1349 solventEquilSimulationEquilibrationLength, 0.1,
1350 solventEquilSimulationEquilibrationLengthNVT, 0.1,
1351 solventEquilSimulationProductionLength, 2.0,
1352 solventSimulationMinimizationSteps, 5000,
1353 solventSimulationEquilibrationLength, 1.0,
1354 solventSimulationProductionLength, 10.0"
1355 --solventParams "positiveIon, Na+, negativeIon, Cl-"
1356
1357 Author:
1358 Manish Sud(msud@san.rr.com)
1359
1360 See also:
1361 OpenFECalculateAbsoluteBindingFreeEnergy.py,
1362 OpenFECalculateAbsoluteHydrationFreeEnergy.py, OpenFECalculatePartialCharges.py,
1363 OpenFECalculateRelativeHydrationFreeEnergy.py, OpenFEGenerateLigandNetwork.py
1364
1365 Copyright:
1366 Copyright (C) 2026 Manish Sud. All rights reserved.
1367
1368 The functionality available in this script is implemented using OpenFE, an
1369 open source molecuar for alchemical free energy calculations.
1370
1371 This file is part of MayaChemTools.
1372
1373 MayaChemTools is free software; you can redistribute it and/or modify it under
1374 the terms of the GNU Lesser General Public License as published by the Free
1375 Software Foundation; either version 3 of the License, or (at your option) any
1376 later version.
1377
1378 """
1379
1380 if __name__ == "__main__":
1381 main()
