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