1 #!/bin/env python
2 #
3 # File: OpenFECalculateAbsoluteHydrationFreeEnergy.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 AbsoluteSolvationProtocol
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 CalculateAbsoluteHydrationFreeEnergy()
100
101 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
102 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
103
104
105 def CalculateAbsoluteHydrationFreeEnergy():
106 """Calculate absolute hydration free energy."""
107
108 # Process input file...
109 InMols = ProcessInputFile()
110
111 # Process molecule names...
112 Mols = ProcessMoleculeNames(InMols)
113
114 # Check for miising partial charges...
115 CheckMissingPartialCharges(Mols)
116
117 # Initialize AHFE protocols...
118 AHFEProtocol = InitializeAbsoluteSolvationProtocol()
119
120 # Initialize solvent...
121 Solvent = InitializeSolventComponent()
122
123 # Setup transformations...
124 MolTransformations = SetupTransformations(Mols, Solvent, AHFEProtocol)
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 InitializeAbsoluteSolvationProtocol():
137 """Initialize absolute solvation protocol."""
138
139 MiscUtil.PrintInfo("\nInitializing absolute solvation protocol...")
140
141 AHFESettings = OpenFEUtil.SetupAbsoluteHydrationFreeEnergySettings("-a, --ahfeParams", OptionsInfo["AHFEParams"])
142 AHFEProtocol = OpenFEUtil.InitializeAbsoluteSolvationFreeEngeryProtocol(AHFESettings)
143
144 return AHFEProtocol
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(Mols, Solvent, AHFEProtocol):
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 in the solvent...
175 MolSolventSystem = OpenFEUtil.InitializeChemicalSystem(
176 SmallMol=Mol, MacroMol=None, Solvent=Solvent, Name="%s_Solvent" % Mol.name
177 )
178
179 # Setup a system for a molecule fully decoupled in the solvent: Only need to use the solvent....
180 SolventOnlySystem = OpenFEUtil.InitializeChemicalSystem(
181 SmallMol=None, MacroMol=None, Solvent=Solvent, Name="Solvent_Only"
182 )
183
184 # Setup a transformation for absolute solvation protocol from MolSolvent to SolventOnly
185 # as the vacuum modelling is automatically handled by AbsoluteSolvationProtocol...
186 TransformationName = "%s_AbsoluteSolvation" % (Mol.name)
187 MolTransformation = OpenFEUtil.InitializeTransformation(
188 StateA=MolSolventSystem,
189 StateB=SolventOnlySystem,
190 Mapping=None,
191 Protocol=AHFEProtocol,
192 Name=TransformationName,
193 Validate=False,
194 )
195
196 MolTransformations.append(MolTransformation)
197
198 # Write out transformatios...
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_AHFE_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 DeltaGHydration = "NA"
277 DeltaGHydrationUncertainty = "NA"
278 else:
279 # Setup hydration value without the units...
280 DeltaGHydration = MolProtocolResult.get_estimate()
281 DeltaGHydration = "%.*f" % (Precision, DeltaGHydration.m)
282
283 # Setup uncertainty value without the units...
284 DeltaGHydrationUncertainty = MolProtocolResult.get_uncertainty()
285 DeltaGHydrationUncertainty = "%.*f" % (Precision, DeltaGHydrationUncertainty.m)
286
287 ResultData.append([MolName, DeltaGHydration, DeltaGHydrationUncertainty])
288
289 ResultDF = pd.DataFrame(ResultData, columns=["MolName", "AHFE DeltaG (kcal/mol)", "Uncertainty (kcal/mol)"])
290 ResultDF.to_csv(ResultFilePath, sep=ResultFileParams["Delim"], lineterminator="\n", index=False)
291
292
293 def ProcessInputFile():
294 """Process input file."""
295
296 # Read small molecule input file...
297 MiscUtil.PrintInfo("\nReading small molecule file %s..." % OptionsInfo["Infile"])
298 Mols, MolCount, ValidMolCount = OpenFEUtil.ReadAndValidateMolecules(
299 OptionsInfo["InfilePath"], **OptionsInfo["InfileParams"]
300 )
301
302 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
303 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
304 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
305
306 if ValidMolCount == 0:
307 MiscUtil.PrintInfo("")
308 MiscUtil.PrintError("No valid molecules found in small molecule input file.\n")
309
310 return Mols
311
312
313 def ProcessMoleculeNames(Mols):
314 """Process molecule names."""
315
316 SpecifiedMols = []
317 if OptionsInfo["FirstMoleculeMode"]:
318 SpecifiedMols.append(Mols[0])
319 elif OptionsInfo["AllMoleculesMode"]:
320 SpecifiedMols = Mols
321 elif OptionsInfo["MoleculesNamesMode"]:
322 SpecifiedMols = OpenFEUtil.ProcessMoleculeNames(Mols, OptionsInfo["MoleculeNamesList"])
323
324 return SpecifiedMols
325
326
327 def CheckMissingPartialCharges(Mols):
328 """Check missing partial charges for small molecules."""
329
330 MiscUtil.PrintInfo("\nChecking missing partial charges for small molecules...")
331
332 MissingChargesMolCount = OpenFEUtil.GetMissingPartialChargesMolCount(Mols)
333 MiscUtil.PrintInfo("Number of molecules with missing partial charges: %s" % MissingChargesMolCount)
334
335 if MissingChargesMolCount == 0:
336 return
337
338 if re.match("^Stop$", OptionsInfo["MissingChargeMode"], re.I):
339 MiscUtil.PrintInfo("")
340 MiscUtil.PrintError(
341 '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 AbsoluteSolvationProtocol module during the calculation of AHFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--ahfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate AHFE.\n'
342 )
343 else:
344 MiscUtil.PrintInfo("")
345 MiscUtil.PrintWarning(
346 'The small molecule input file contains molecules with missing partial charges. The missing charges will be automatically calculated by OpenFE AbsoluteSolvationProtocol module during the calculation of AHFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--ahfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate AHFE.\n'
347 )
348
349
350 def ProcessMoleculeNamesOption():
351 """Process molecule names Option."""
352
353 OptionsInfo["MoleculeNames"] = Options["--moleculeNames"]
354 OptionsInfo["MoleculeNamesList"] = None
355
356 if OptionsInfo["MoleculeNames"] is None:
357 return
358
359 MoleculeNamesList = []
360 for MoleculeName in OptionsInfo["MoleculeNames"].split(","):
361 MoleculeNamesList.append(MoleculeName.strip())
362
363 OptionsInfo["MoleculeNamesList"] = MoleculeNamesList
364
365
366 def ProcessOutfilePrefixOption():
367 """Process outfile prefix option."""
368
369 OutfilePrefix = Options["--outfilePrefix"]
370
371 if re.match("^auto$", OutfilePrefix, re.I):
372 OutfilePrefix = OptionsInfo["InfileRoot"]
373
374 OptionsInfo["OutfilePrefix"] = OutfilePrefix
375
376
377 def ProcessOutfileDirOption():
378 """Process outfile directory Option."""
379
380 # Setup output directory...
381 OutfileDir = Options["--outfileDir"]
382 OutfileDirPath = os.path.abspath(OutfileDir)
383 if not os.path.exists(OutfileDir):
384 MiscUtil.PrintInfo("\nCreating output directory %s..." % (OutfileDir))
385 os.mkdir(OutfileDirPath)
386 OptionsInfo["OutfileDir"] = OutfileDir
387 OptionsInfo["OutfileDirPath"] = OutfileDirPath
388
389 # Setup a transformations subdirectory...
390 TransformationsOutDir = "Transformations"
391 TransformationsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], TransformationsOutDir)
392 if not os.path.exists(TransformationsOutDirPath):
393 os.mkdir(TransformationsOutDirPath)
394 OptionsInfo["TransformationsOutDir"] = TransformationsOutDir
395 OptionsInfo["TransformationsOutDirPath"] = TransformationsOutDirPath
396
397 # Setup a results subdirectory...
398 ResultsOutDir = "Results"
399 ResultsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], ResultsOutDir)
400 if not os.path.exists(ResultsOutDirPath):
401 os.mkdir(ResultsOutDirPath)
402 OptionsInfo["ResultsOutDir"] = ResultsOutDir
403 OptionsInfo["ResultsOutDirPath"] = ResultsOutDirPath
404
405 # Use results subdirectory for scratch results...
406 OptionsInfo["ResultsScratchOutDir"] = ResultsOutDir
407 OptionsInfo["ResultsScratchOutDirPath"] = ResultsOutDirPath
408
409
410 def ProcessListOption():
411 """Process list protocol settings option."""
412
413 AHFESettings = AbsoluteSolvationProtocol.default_settings()
414
415 MiscUtil.PrintInfo("\nListing AHFE settings...")
416 OpenFEUtil.ListOpenFESettings(AHFESettings)
417
418
419 def ConfigureLogging():
420 """Configure logging."""
421
422 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
423
424 if re.match("^Error$", OptionsInfo["LoggingLevel"], re.I):
425 LoggingLevel = logging.ERROR
426 elif re.match("^Warning$", OptionsInfo["LoggingLevel"], re.I):
427 LoggingLevel = logging.WARNING
428 else:
429 LoggingLevel = logging.INFO
430
431 logging.basicConfig(format="%(levelname)s: %(message)s", level=LoggingLevel)
432
433 # Turn warnings issued by warnings.warn() into log message to avoid display
434 # of a stack trace...
435 logging.captureWarnings(True)
436
437
438 def ProcessOptions():
439 """Process and validate command line arguments and options."""
440
441 MiscUtil.PrintInfo("Processing options...")
442
443 # Validate options...
444 ValidateOptions()
445
446 # Configure logging...
447 ConfigureLogging()
448
449 OptionsInfo["Infile"] = Options["--infile"]
450 OptionsInfo["InfilePath"] = os.path.abspath(OptionsInfo["Infile"])
451 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
452 OptionsInfo["InfileRoot"] = FileName
453
454 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
455 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
456 "--infileParams",
457 Options["--infileParams"],
458 InfileName=Options["--infile"],
459 ParamsDefaultInfo=ParamsDefaultInfoOverride,
460 )
461
462 ParamsDefaultInfoOverride = {"SolventEngineComputePlatform": "CPU", "VacuumEngineComputePlatform": "CPU"}
463 ParamsDefaultInfoOverride = None
464 OptionsInfo["AHFEParams"] = OpenFEUtil.ProcessOptionOpenFEAbsoluteHydrationFreeEnergyParameters(
465 "--ahfeParams", Options["--ahfeParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
466 )
467
468 OptionsInfo["ExecuteDAGParams"] = OpenFEUtil.ProcessOptionOpenFEExecuteDAGParameters(
469 "--executeDAGParams", Options["--executeDAGParams"]
470 )
471 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
472
473 OptionsInfo["Mode"] = OpenFEUtil.ProcessOptionOpenFEAbsoluteFreeEnergyMode("-m, --mode", Options["--mode"])
474 OptionsInfo["FirstMoleculeMode"] = True if re.match("^FirstMolecule$", OptionsInfo["Mode"], re.I) else False
475 OptionsInfo["AllMoleculesMode"] = True if re.match("^AllMolecules$", OptionsInfo["Mode"], re.I) else False
476 OptionsInfo["MoleculesNamesMode"] = True if re.match("^MoleculeNames$", OptionsInfo["Mode"], re.I) else False
477
478 OptionsInfo["MissingChargeMode"] = OpenFEUtil.ProcessOptionOpenFEMissingChargeMode(
479 "--missingChargeMode", Options["--missingChargeMode"]
480 )
481
482 ProcessMoleculeNamesOption()
483
484 OptionsInfo["ResultFileParams"] = OpenFEUtil.ProcessOptionOpenFEResultFileParameters(
485 "--resultFileParams", Options["--resultFileParams"]
486 )
487 OptionsInfo["SolventParams"] = OpenFEUtil.ProcessOptionOpenFESolventParameters(
488 "--solventParams", Options["--solventParams"]
489 )
490
491 ProcessOutfilePrefixOption()
492 ProcessOutfileDirOption()
493
494 OptionsInfo["Overwrite"] = Options["--overwrite"]
495
496 # Track top level working directory...
497 OptionsInfo["TopWorkingDir"] = os.getcwd()
498
499
500 def RetrieveOptions():
501 """Retrieve command line arguments and options."""
502
503 # Get options...
504 global Options
505 Options = docopt(_docoptUsage_)
506
507 # Set current working directory to the specified directory...
508 WorkingDir = Options["--workingdir"]
509 if WorkingDir:
510 os.chdir(WorkingDir)
511
512 # Handle examples option...
513 if "--examples" in Options and Options["--examples"]:
514 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
515 sys.exit(0)
516
517
518 def ValidateOptions():
519 """Validate option values."""
520
521 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
522 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd")
523
524 MiscUtil.ValidateOptionDirPath("-o, --outfileDir", Options["--outfileDir"])
525 MiscUtil.ValidateOptionsOutputDirOverwrite(
526 "-o, --outfileDir", Options["--outfileDir"], "--overwrite", Options["--overwrite"]
527 )
528
529 MiscUtil.ValidateOptionTextValue("--loggingLevel", Options["--loggingLevel"], "Info Warning Error")
530
531 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "FirstMolecule AllMolecules MoleculeNames")
532 MiscUtil.ValidateOptionTextValue("--missingChargeMode", Options["--missingChargeMode"], "Calculate Stop")
533
534 if re.match("^MoleculeNames$", Options["--mode"], re.I):
535 if MiscUtil.IsEmpty(Options["--moleculeNames"]):
536 MiscUtil.PrintError(
537 'You must specify a value for "--moleculeNames" option during, MoleculeNames, value for "-m, --mode" option.'
538 )
539
540
541 # Setup a usage string for docopt...
542 _docoptUsage_ = """
543 OpenFECalculateAbsoluteHydrationFreeEnergy.py - Calculate absolute hydration free energy
544
545 Usage:
546 OpenFECalculateAbsoluteHydrationFreeEnergy.py [--ahfeParams <Name,Value,...>] [--executeDAGParams <Name,Value,..>]
547 [--infileParams <Name,Value,...>] [--loggingLevel <Info, Warning or Error>] [--mode <FirstMolecule, AllMolecules, or ...>]
548 [--missingChargeMode <Calculate or Stop>] [--moleculeNames <MolName1,MolName2,..>] [--outfilePrefix <text>]
549 [--resultFileParams <Name,Value,..>] [--solventParams <Name,Value,...>] [--overwrite]
550 [-w <dir>] -i <infile> -o <outifiledir>
551 OpenFECalculateAbsoluteHydrationFreeEnergy.py -l | --list
552 OpenFECalculateAbsoluteHydrationFreeEnergy.py -h | --help | -e | --examples
553
554 Description:
555 Calculate Absolute Hydration Free Energy (AHFE) for molecules in a small
556 molecule input file. You may calculate AHFEs for specific molecules or all
557 molecules in the input file.
558
559 The small molecule input file must contain molecules already prepared for
560 simulation. It must contain appropriate 3D coordinates along with no missing
561 hydrogens.
562
563 The MD simulation workflow, employed for the calculation of AHFEs, involves
564 the following steps: initial minimization; NVT equilibration; NPT equilibration;
565 production NPT. The MD simulation protocol is repeated 3 times for solvent
566 transformation from MolSolventToSolventOnly and vacumm simulation, and
567 he results are analyzed to estimate AHFEs. The default time and step size
568 settings for the MD protocol are shown below:
569
570 Protocol repeats, 3
571
572 Time step size: 4.0 femtosecond
573
574 Solvent equilibration phase:
575
576 Max minimization steps: 5,000
577 NVT equilibration length: 0.1 nanosecond
578 NPT equilibration length: 0.2 nanosecond
579 NPT length: 0.5 nanosecond
580
581 Solvent production phase:
582
583 Max minimization steps: 5,000
584 NPT equilibration length: 1.0 nanosecond
585 NPT length: 10.0 nanosecond
586
587 Vacuum equilibration phase:
588
589 Max minimization steps: 5,000
590 NVT equilibration length: None
591 NPT equilibration length: 0.2 nanosecond
592 NPT length: 0.5 nanosecond
593
594 Vacuum production phase:
595
596 Max minimization steps: 5,000
597 NVT equilibration length: None
598 NPT equilibration length: 0.5 nanosecond
599 NPT length: 2.0 nanosecond
600
601 Each solvent and vacuum simulation, by default, may run for 11.8 and 3.2
602 nanosecond respectively, for a total of 15 nanoseconds. The total MD
603 simulation time for correspond to 45 nanosecond to repeat the protocol
604 3 times for the solvent and vacuum simulations.
605
606 Possible outfile prefix:
607
608 <OutfilePrefix> or <InfileRoot>
609
610 Possible output directories:
611
612 <OutfileDir>
613
614 <OutfileDir>/Transformations
615 <OutfileDir>/Results
616
617 Possible output files and directories under <OutfileDir>:
618
619 <OutfilePrefix>_AHFE_Results.<csv or tsv>
620
621 Transformations/<MolName>_Solvent_To_Solvent_Only.json
622 ... ... ...
623
624 Results/<MolName>_AbsoluteSolvation_Results.json
625
626 Results/shared_AbsoluteSolvationSolventUnit-*/
627 Results/scratch_AbsoluteSolvationVacuumUnit-*/
628 ... ... ...
629
630 Options:
631 -a, --ahfeParams <Name,Value,...> [default: auto]
632 A comma delimited list of parameter name and value pairs for AHFE protocol
633 settings employed during the calculation of AHFEs.
634
635 The default values are automatically updated to match settings provided by
636 OpenFE module AbsoluteSolvationProtocol.
637
638 You must specify valid OpenFE values for these parameters. An extensive
639 validation is not performed.
640
641 The supported parameter names along with their default values are
642 are shown below:
643
644 protocolRepeats, 3
645
646 Integrator settings:
647
648 integratorBarostatFrequency, 25.0 * timestep [ The specified value
649 is a multiple of integratorTimestep. ]
650 integratorConstraintTolerance, 1e-06
651 integratorLangevinCollisionRate, 1.0 [ Units: 1 / picosecond ]
652 integratorNRestartAttempts, 20
653 integratorReassignVelocities, no [ Possible values: yes or no ]
654 integratorRemoveCom, no [ Possible values: yes or no ]
655 integratorTimestep, 4.0 [ Units: femtosecond ]
656
657 Lambda settings:
658
659 lambdaElec, 0.0 0.25 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
660 1.0 [ Possible values: A space delimited list of values
661 between 0.0 and 1.0 ]
662 lambdaRestraints, 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
663 0.0 0.0 [ Possible values: A space delimited list of values
664 between 0.0 and 1.0 ]
665 lambdaVdw, [0.0 0.0 0.0 0.0 0.0 0.12 0.24 0.36 0.48 0.6 0.7 0.77
666 0.85 1.0 [ Possible values: A space delimited list of values
667 between 0.0 and 1.0 ]
668
669 Partial charge settings:
670
671 partialChargeNaglModel, None [ Default: Production AM1BCC model for
672 NAGL; Possible value: Any valid name. ]
673 partialChargeNumberOfConformers, None [ Possible value: > 0 ]
674 partialChargeOffToolkitBackend, AmberTools [ Possible values:
675 AmberTools or RDKit ]
676 partialChargeMethod, AM1BCC [ Possble values: AM1BCC, Espaloma,
677 or NAGL ]
678
679 Solvation settings:
680
681 solvationBoxShape, dodecahedron [ Possible values: cube,
682 dodecahedron, or octahedron ]
683 solvationBoxSize, None [ Possible value: A triplet of space
684 X Y Z values; Units: nanometer ]
685 solvationSolventModel, tip3p [ Possible values: tip3p, spce, tip4pew,
686 or tip5p ]
687 solvationSolventPadding, 1.5 [ Units: nanometer ]
688
689 Solvent engine settings:
690
691 solventEngineComputePlatform, CPU [ Possible values: CPU, CUDA,
692 OpenCL, or Reference ]
693 solventEngineGpuDeviceIndex, None [ Possible values: 0, 0 1, etc. ]
694
695 Solvent equil output settings:
696
697 solventEquilOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
698 solventEquilOutputCheckpointStorageFilename, checkpoint.chk
699 solventEquilOutputNPTStructure, equil_npt_structure.pdb
700 solventEquilOutputNVTStructure, equil_nvt_structure.pdb
701 solventEquilOutputForcefieldCache, db.json
702 solventEquilOutputLogOutput, equil_simulation.log
703 solventEquilOutputMinimizedStructure, minimized.pdb
704 solventEquilOutputIndices, not water [ Possible value: Any valid
705 selection. ]
706 solventEquilOutputPreminimizedStructure, system.pdb
707 solventEquilOutputProductionTrajectoryFilename, production_equil.xtc
708 solventEquilOutputTrajectoryWriteInterval, 20.0 [ Units: picosecond ]
709
710 Solvent equil simulation settings:
711
712 solventEquilSimulationEquilLength, 0.2 [ Units: nanosecond ]
713 solventEquilSimulationEquiLengthNVT, 0.1 [ Units: nanosecond ]
714 solventEquilSimulationMinimizationSteps,5000
715 solventEquilSimulationProductionLength,0.5 [ Units: nanosecond ]
716
717 Solvent forcefield settings:
718
719 solventForcefieldConstraints, HBonds [ Possible values: HBonds,
720 AllBonds, or HAngles ]
721 solventForcefields, amber/ff14SB.xml, amber/tip3p_standard.xml
722 amber/tip3p_HFE_multivalent.xml amber/phosaa10.xml
723 [ Possible values: A space delimited list of valid names. ]
724 solventForcefieldHydrogenMass, 3.0 [ Units: amu ]
725 solventForcefieldNonbondedCutoff,0.9 [ Units: nanometer ]
726 solventForcefieldNonbondedMethod, PME [ Possible values: PME or
727 NoCutoff ]
728 solventForcefieldRigidWater, yes, [ Possible values: yes or no ]
729 solventForcefieldSmallMoleculeForcefield, openff-2.1.1 [ Possible
730 value: A valid forcefield name. ]
731
732 Solvent output settings:
733
734 solventOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
735 solventOutputCheckpointStorageFilename, solvent_checkpoint.nc
736 solventOutputForcefieldCache, db.json
737 solventOutputFilename, solvent.nc
738 solventOutputIndices, not water [ Possible value: Any valid
739 selection. ]
740 solventOutputStructure, hybrid_system.pdb
741 solventOutputPositionsWriteFrequency, 100.0 [ Units: picosecond ]
742 solventOutputVelocitiesWriteFrequency, None [ Possible
743 values: > 0; Units: picosecond ]
744
745 Solvent simulation settings:
746
747 solventSimulationEarlyTerminationTargetError, 0.0 [ Units:
748 kilocalorie_per_mole ]
749 solventSimulationEquilibrationLength, 1.0 [ Units: nanosecond ]
750 solventSimulationMinimizationSteps, 5000
751 solventSimulationNReplicas, 14
752 solventSimulationProductionLength, 10.0 [ Units: nanosecond ]
753 solventSimulationRealTimeAnalysisInterval, 250.0 [ Units:
754 picosecond ]
755 solventSimulationRealTimeAnalysisMinimumTime, 500.0 [ Units:
756 picosecond
757 solventSimulationSamplerMethod, repex [ Possible values: repex,
758 sams, or independent ]
759 solventSimulationSamsFlatnessCriteria, logZ-flatness [ Possible
760 values: logZ-flatness, minimum-visits or histogram-flatness ]
761 solventSimulationsamsGamma0, 1.0
762 solventSimulationTimePerIteration,2.5 [ Units: picosecond ]
763
764 Thermo settings:
765
766 thermoPh, None [ Possible values: > 0 ]
767 thermoPressure, 1.0 [ Units: bar ]
768 thermoRedoxPotential, None [ Possible values: A valid float.
769 Units: millivolts (mV) ]
770 thermoTemperature, 298.15 [ Units: kelvin ]
771
772 Vacuum engine settings:
773
774 vacuumEngineComputePlatform, CPU [ Possible values: CPU, CUDA,
775 OpenCL, or Reference ]
776 vacummEngineGpuDeviceIndex, None [ Possible values: 0, 0 1, etc. ]
777
778 Vacuum equil output settings:
779
780 vacuumEquilOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
781 vacuumEquilOutputCheckpointStorageFilename, checkpoint.chk
782 vacuumEquilOutputNPTStructure, equil_structure.pdb
783 vacuumEquilOutputNVTStructure,None
784 vacuumEquilOutputForcefieldCache, db.json
785 vacuumEquilOutputLogOutput, equil_simulation.log
786 vacuumEquilOutputMinimizedStructure, minimized.pdb
787 vacuumEquilOutputIndices, not water [ Possible value: Any valid
788 selection. ]
789 vacuumEquilOutputPreminimizedStructure, system.pdb
790 vacuumEquilOutputProductionTrajectoryFilename, production_equil.xtc
791 vacuumEquilOutputTrajectoryWriteInterval, 20.0 [ Units: picosecond ]
792
793 Vacuum equil simulation settings:
794
795 vacuumEquilSimulationEquilLength, 0.2 [ Units: nanosecond ]
796 vacuumEquilSimulationEquilLengthNVT, None [ Units: nanosecond ]
797 vacuumEquilSimulationMinimizationSteps, 5000
798 vacuumEquilSimulationProductionLength, 0.5 [ Units: nanosecond ]
799
800 Vacuum forcefield settings:
801
802 vacuumForcefieldConstraints, HBonds [ Possible values: HBonds,
803 AllBonds, or HAngles ]
804 vacuumForcefields, amber/ff14SB.xml, amber/tip3p_standard.xml
805 amber/tip3p_HFE_multivalent.xml amber/phosaa10.xml
806 [ Possible values: A space delimited list of valid names. ]
807 vacuumForcefieldHydrogenMass, 3.0 [ Units: amu ]
808 vacuumForcefieldNonbondedCutoff, 0.9 [ Units: nanometer ]
809 vacuumForcefieldNonbondedMethod, nocutoff [ Possible values: PME
810 or NoCutoff ]
811 vacuumForcefieldRigidWater, yes, [ Possible values: yes or no ]
812 vacuumForcefieldSmallMoleculeForcefield, openff-2.1.1 [ Possible
813 value: A valid forcefield name. ]
814
815 Vacuum output settings:
816
817 vacuumOutputCheckpointInterval, 1.0 [ Units: nanosecond ]
818 vacuumOutputCheckpointStorageFilename, vacuum_checkpoint.nc
819 vacuumOutputForcefieldCache, db.json
820 vacuumOutputFilename, vacuum.nc
821 vacuumOutputIndices, not water [ Possible value: Any valid
822 selection. ]
823 vacuumOutputStructure, hybrid_system.pdb
824 vacuumOutputPositionsWriteFrequency, 100.0 [ Units: picosecond ]
825 vacuumOutputVelocitiesWriteFrequency, None [ Possible
826 values: > 0; Units: picosecond ]
827
828 Vacuum simulation settings:
829
830 vacuumSimulationEarlyTerminationTargetError, 0.0 [ Units:
831 0.0 kilocalorie_per_mole ]
832 vacuumSimulationEquilibrationLength, 0.5 [ Units: nanosecond ]
833 vacuumSimulationMinimizationSteps, 5000
834 vacuumSimulationNReplicas, 14
835 vacuumSimulationProductionLength, 2.0 [ Units: nanosecond ]
836 vacuumSimulationRealTimeAnalysisInterval, 250.0 [ Units: picosecond ]
837 vacuumSimulationRealTimeAnalysisMinimumTime, 500.0 [ Units: picosecond]
838 vacuumSimulationSamplerMethod, repex [ Possible values: repex,
839 sams, or independent ]
840 vacuumSimulationSamsFlatnessCriteria, logZ-flatness [ Possible
841 values: logZ-flatness, minimum-visits or histogram-flatness ]
842 vacuumSimulationSamsGamma0, 1.0
843 vacuumSimulationTimePerIteration,2.5 [ Units: picosecond ]
844
845 A brief description of parameters, taken from OpenFE documentation, is
846 provided below:
847
848 protocolRepeats: Number of completely independent repeats of the
849 entire sampling process.
850
851 Integrator settings:
852
853 Parameters controlling the LangevinSplittingDynamicsMove integrator
854 used for simulation.
855
856 integratorBarostatFrequency: Frequency at which volume scaling
857 changes should be attempted.
858 integratorConstraintTolerance: Tolerance for constraint solver.
859 integratorLangevinCollisionRate: Collision frequency.
860 integratorNRestartAttempts: Number of attempts to restart from
861 Context in case there are NaNs in the energies after
862 integration.
863 integratorReassignVelocities: Reassign velocities from the
864 Maxwell-Boltzmann distribution at the beginning of each
865 Monte Carlo move.
866 integratorRemoveCom: Remove the center of mass motion.
867 integratorTimestep: Size of the simulation timestep.
868
869 Lambda settings:
870
871 Lambda protocol parameters.
872
873 lambdaElec: List of lambda values for electrostatics. The values of
874 0 and 1 imply state A and state B respectively.
875 lambdaRestraints: List of lambda values for restraints. The values
876 of 0 and 1 imply state A and state B respectively.
877 lambdaVdw: List of lamda values for van der Waals. The values of
878 of 0 and 1 imply state A and state B respectively.
879
880 Partial charge settings:
881
882 Parameters for automatically assigning missing partial charges to
883 small molecules, including the partial charge method.
884
885 partialChargeNaglModel: Model to use for partial charge assignment.
886 A value of None implies the use of the latest available
887 production AM1BCC model.
888 partialChargeNumberOfConformers: Number of conformers to generate
889 as part of the partial charge assignment. A value of None
890 implies the use of the existing conformer.
891 partialChargeOffToolkitBackend: OpenFF toolkit registry backend to
892 use for calculating partial charges.
893 partialChargeMethod: Method to use for calculating partial charges.
894
895 Solvation settings:
896
897 Solvation parameters for the system, including the solvent model and
898 the solvent padding.
899
900 solvationBoxShape: Shape of the periodic solvent box to create.
901 solvationBoxSize: Lengths of the unit cell for a solvent box.
902 solvationSolventModel: Forcefield water model to use during
903 solvation and defining the model properties.
904 solvationSolventPadding: Minimum distance from any solute bounding
905 sphere to the edge of the box.
906
907 Solvent engine settings:
908
909 Parameters configuring the compute platform used by the OpenMM to
910 perform the simulation.
911
912 solventEngineComputePlatform: Platform to use for running OpenMM MD
913 calculations.
914 solventEngineGpuDeviceIndex: Space delimited list of device indices
915 to use for running OpenMM MD calculations.
916
917 Solvent equil output settings:
918
919 Parameters controlling simulation output during equilibration
920 phase of solvent transformation.
921
922 solventEquilOutputCheckpointInterval: Frequency to write the
923 checkpoint file.
924 solventEquilOutputCheckpointStorageFilename: Checkpoint filename.
925 solventEquilOutputNPTStructure: NPT structure filename.
926 solventEquilOutputNVTStructure: NVT strucure filename.
927 solventEquilOutputForcefieldCache: Filename for caching small
928 molecule residue templates.
929 solventEquilOutputLogOutput: Simulation log filename.
930 solventEquilOutputMinimizedStructure: Minimized structire filename.
931 solventEquilOutputIndices: Selection string for selecting
932 coordinates to write.
933 solventEquilOutputPreminimizedStructure: Initial structure filename.
934 solventEquilOutputProductionTrajectoryFilename: Trajectory filename.
935 solventEquilOutputTrajectoryWriteInterval: Frequency for writing
936 velocities to trajectory file.
937
938 Solvent equil simulation settings:
939
940 Parameters controlling simulation during equilibration phase of
941 solvent transformation.
942
943 solventEquilSimulationEquilLength: Length of the NPT equilibration
944 phase.
945 solventEquilSimulationEquiLengthNVT: Length of the NVT equilibration
946 phase.
947 solventEquilSimulationMinimizationSteps: Maximum number of
948 minimization steps to perform.
949 solventEquilSimulationProductionLength: Length of the NPT production
950 phase.
951
952 Solvent forcefield settings:
953
954 Parameters to set up the force field with OpenMM Force Fields
955 equilibration phase of solvent transformation.
956
957 solventForcefieldConstraints: Constraints to use.
958 solventForcefields: List of valid forcefield paths for all
959 components except small molecules.
960 solventForcefieldHydrogenMass: Mass to be repartitioned to
961 hydrogens from neighboring heavy atoms.
962 solventForcefieldNonbondedCutoff: Cutoff for short range nonbonded
963 interactions.
964 solventForcefieldNonbondedMethod: Method for treating nonbonded
965 interactions.
966 solventForcefieldRigidWater: Use a rigid water model.
967 solventForcefieldSmallMoleculeForcefield: A valid forcefield name
968 to use for small molecules.
969
970 Solvent output settings:
971
972 Parameters controlling simulation output during final phase of
973 solvent transformation.
974
975 solventOutputCheckpointInterval: Frequency to write the checkpoint
976 file.
977 solventOutputCheckpointStorageFilename: Checkpoint filename.
978 solventOutputForcefieldCache: Filename for caching small molecule
979 residue templates.
980 solventOutputFilename: Trajectory filename.
981 solventOutputIndices: Selection string for selecting coordinates to
982 write.
983 solventOutputStructure: Hybrid topology structure filename.
984 solventOutputPositionsWriteFrequency: Frequency for writing
985 positions to trajectory file.
986 solventOutputVelocitiesWriteFrequency: Frequency for writing
987 velocities to trajectory file.
988
989 Solvent simulation settings:
990
991 Parameters controlling simulation during final phase of solvent
992 transformation.
993
994 solventSimulationEarlyTerminationTargetError: Target error for the
995 real time analysis measured in kcal/mol. Once the MBAR error of
996 the free energy is at or below this value, the simulation will
997 be considered complete. The suggested value of 0.12 has shown to
998 be effective in both hydration and binding free energy
999 benchmarks.
1000 solventSimulationEquilibrationLength: Length of the equilibration
1001 phase. The specified value must be divisible by 'integratorTimestep'.
1002 solventSimulationMinimizationSteps: Maximum number of minimization
1003 steps to perform.
1004 solventSimulationNReplicas: Number of replicas to use.
1005 solventSimulationProductionLength: Length of the production phase.
1006 The specified value must be divisible by 'integratorTimestep'.
1007 solventSimulationRealTimeAnalysisMinimumTime: Time interval for
1008 performing analysis of the free energies. At each interval, real
1009 time analysis data will be written to a yaml file named
1010 <outputFileName>_real_time_analysis.yaml. The current error
1011 in the estimate will also be assessed and the simulation will
1012 be terminated when it drops below
1013 'solventSimulationEarlyTerminationTargetError'.
1014 solventSimulationSamplerMethod: Alchemical sampling method to use:
1015 REPEX (Hamiltonian REPlica EXchange), SAMS (Self-Adjusted
1016 Mixture Sampling), or Independent (Independently sampled lambda
1017 windows).
1018 solventSimulationSamsFlatnessCriteria:Method for assessing when to
1019 switch to asymptomatically optimal scheme for SAMS.
1020 solventSimulationsamsGamma0: Initial weight adaptation rate for
1021 SAMS.
1022 solventSimulationTimePerIteration: Simulation time between each
1023 MCMC move attempt
1024
1025 Vacuum engine settings:
1026
1027 Parameters configuring the compute platform used by the OpenMM to
1028 perform the simulation.
1029
1030 vacuumEngineComputePlatform: Platform to use for running OpenMM MD
1031 calculations.
1032 vacuumEngineGpuDeviceIndex: Space delimited list of device indices
1033 to use for running OpenMM MD calculations.
1034
1035 The rest of the vacuum settings are similar to the solvent settings
1036 already described under various sections for solvent. The prefix
1037 'vacuum' is used for the names of the pramaters instead of the
1038 prefix 'solvent.'
1039
1040 Thermo settings:
1041
1042 Thermodynamic parameters, including the temperature and the pressure
1043 of the system.
1044
1045 thermoPh: Simulation pH
1046 thermoPressure: Simulation pressure.
1047 thermoRedoxPotential:Simulation redox potential.
1048 thermoTemperature: Simulation temperature.
1049
1050 -e, --examples
1051 Print examples.
1052 --executeDAGParams <Name,Value,..> [default: auto]
1053 A comma delimited list of parameter name and value pairs for executing
1054 protocol DAGs (Directed Acyclic Graph) to run AHFE calculations.
1055
1056 The supported parameter names along with their default values are
1057 are shown below:
1058
1059 keepShared, yes [ Possible values: yes or no ]
1060 keepScratch, no [ Possible values: yes or no ]
1061 nRetries, 2 [ Possible values: >= 0. A value of 0 implies only
1062 1 try. ]
1063
1064 A brief description of parameters is provided below:
1065
1066 keepShared: Keep shared directories after the execution of DAG.
1067 keepScratch: Keep scratch directories after the execution of DAG.
1068 nRetries: Number of times to attempt the execution.
1069
1070 -h, --help
1071 Print this help message.
1072 -i, --infile <infile>
1073 Input file containing small molecules.
1074 --infileParams <Name,Value,...> [default: auto]
1075 A comma delimited list of parameter name and value pairs for reading
1076 molecules from files. The supported parameter names for different file
1077 formats, along with their default values, are shown below:
1078
1079 SD: removeHydrogens,no,sanitize,yes,strictParsing,yes
1080
1081 -l, --list
1082 List default AHFE protocol settings provided by OpenFE module
1083 AbsoluteSolvationProtocol.
1084 --loggingLevel <Info, Warning or Error> [default: Error]
1085 Logging level to configure the 'root logger' via logging.basicConfig()
1086 function. The default logging level is changed from 'logging.INFO' to
1087 'logging.ERROR'. Otherwise, OpenFE and its associated modules
1088 may generate a lot of informational messages.
1089 -m, --mode <FirstMolecule, AllMolecules, or ...> [default: FirstMolecule]
1090 Calculate AHFE for the first molecule, the specified molecule names,
1091 or all molecules in an input file. Possible values: FirstMolecule,
1092 AllMolecules, or MoleculeNames. You must specify a comma delimited list
1093 of molecule names using '--moleculeNames'option during 'MoleculeNames'
1094 value for '--mode' option.
1095 --missingChargeMode <Calculate or Stop> [default: Stop]
1096 Calculate missing partial charges for molecules before running AHFE
1097 calculations or terminate the execution of the script. The missing
1098 partial charges will be automatically calculated by OpenFE module
1099 AbsoluteSolvationProtocol during the calculation of AHFE. You
1100 may control the calculation of partial charges by specifying values for
1101 partialCharge* parameters using '--ahfeParams' option.
1102 --moleculeNames <MolName1,MolName2,..>
1103 A comma delimited list of molecule names for calculating AHFEs.
1104 This option is only used during 'MoleculeNames' value for
1105 '--mode' option.
1106 -o, --outfileDir <outfiledir>
1107 Output directory.
1108 --outfilePrefix <text> [default: auto]
1109 Prefix for generating output files under output directory.
1110 --overwrite
1111 Overwrite existing files.
1112 --resultFileParams <Name,Value,..> [default: auto]
1113 A comma delimited list of parameter name and value pairs for writing
1114 calculated RHFEs values to a results file.
1115
1116 The supported parameter names along with their default values are
1117 are shown below:
1118
1119 precision, 4 [ Possible values: > 0 ]
1120 delimiter, comma [ Possible values: comma or tab ]
1121
1122 --solventParams <Name,Value,...> [default: auto]
1123 A comma delimited list of parameter name and value pairs for solvent
1124 component. You must specify valid OpenFE values. No extensive validation
1125 is performed. These parameters are used in conjunction with solvation*
1126 parameters available through '--ahfeParams' to perform solvation.
1127
1128 The supported parameter names along with their default values are
1129 are shown below:
1130
1131 positiveIon, Na+ [ Possible value: Li+, Na+, K+, Rb+, or Cs+ ]
1132 negativeIon, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
1133 neutralize, yes [ Possible values: yes or no ]
1134 ionConcentration, 0.15 [ Units: molar ]
1135
1136 A brief description of parameters is provided below:
1137
1138 positiveIon, negativeion: Pair of ions used to neutralize and bring
1139 the solvent to required ionic concentration.
1140 neutralize: Neutralize the net charge on the chemical state by the
1141 ions in the solvent component.
1142 ionConcentration: Ionic concentration.
1143
1144 -w, --workingdir <dir>
1145 Location of working directory which defaults to the current directory.
1146
1147 Examples:
1148 The sample protein and ligand files for tyrosine kinase 2 (Tyk2) are
1149 distributed with MayaChemTools and are available in data directory. These
1150 files have been taken from OpenFE distribution for example notebooks. The
1151 AM1BCC partial charges have been calculated for the ligands in SD file to
1152 facilitate calculations. You may review OpenFE tutorial notebooks for the
1153 expected results.
1154
1155 To calcuate AHFE for the first molecule in a SD file, performing 3 independent
1156 repeats of the entire MD sampling process to estimate AHFE for the molecule,
1157 each solvent and vacuum MD repeat consisting of equilibration phase ( Solvent:
1158 Minimization - 5,000; NVT - 0.1 ns; NPT - 0.2; NPT prod - 0.5 ns; Vacuum:
1159 Minimization - 5,000; NVT - None; NPT - 0.2 ns; NPT prod - 0.5 ns) and
1160 production phase ( Solvent: Minimization - 5,000; NPT equil - 1.0 ns; NPT prod:
1161 10.0 ns; Vacuum: Minimization - 5,000; NVT equil - None; NPT equil - 0.5 ns;
1162 NPT prod - 2 ns) using a step size of of 4 fs, writing out appropriate trajectory
1163 and PDB files for each MD repeat in Results subdirectory under output directory,
1164 type:
1165
1166 % OpenFECalculateAbsoluteHydrationFreeEnergy.py
1167 -i SampleTyk2Ligands.sdf -o SampleTyk2LigandsAHFE
1168
1169 To run the first example for calculating AHFE for specific molecules using CUDA
1170 platform on your machine to perform solvent and vacuum MD simulations and
1171 generate various output files, type:
1172
1173 % OpenFECalculateAbsoluteHydrationFreeEnergy.py
1174 -i SampleTyk2Ligands.sdf -o SampleTyk2LigandsAHFE -m MoleculeNames
1175 --moleculeNames "lig_ejm_31, lig_ejm_47"
1176 --ahfeParams "solventEngineComputePlatform,CUDA,
1177 vacuumEngineComputePlatform,CUDA"
1178
1179 To run the second example to see all warning messages produced by OpenFE
1180 modules and write various output files, type;
1181
1182 % OpenFECalculateAbsoluteHydrationFreeEnergy.py
1183 -i SampleTyk2Ligands.sdf -o SampleTyk2LigandsAHFE -m MoleculeNames
1184 --moleculeNames "lig_ejm_31, lig_ejm_47"
1185 --ahfeParams "solventEngineComputePlatform,CUDA,
1186 vacuumEngineComputePlatform,CUDA"
1187 --loggingLevel Warning
1188
1189 To run the first example for calculating AHFE for all molecules using CUDA
1190 platform on your machine to perform solvent and vacuum MD simulations,
1191 automatically calculate missing partial charges for molecules, and generate
1192 various output files, type:
1193
1194 % OpenFECalculateAbsoluteHydrationFreeEnergy.py
1195 -i SampleTyk2LigandsNoCharges.sdf -o SampleTyk2LigandsAHFE
1196 -m AllMolecules --ahfeParams "solventEngineComputePlatform,CUDA,
1197 vacuumEngineComputePlatform,CUDA"
1198 --missingChargeMode Calculate
1199
1200 To run the second example by specifying explict values for various parametres
1201 and generate various output files, type:
1202
1203 % OpenFECalculateAbsoluteHydrationFreeEnergy.py
1204 -i SampleTyk2Ligands.sdf -o SampleTyk2LigandsAHFE -m MoleculeNames
1205 --moleculeNames "lig_ejm_31, lig_ejm_47"
1206 --loggingLevel Error
1207 --executeDAGParams "keepShared, yes, nRetries, 2"
1208 --missingChargeMode Stop --ahfeParams "protocolRepeats,3,
1209 solventEngineComputePlatform,CUDA, vacuumEngineComputePlatform,CUDA,
1210 integratorTimestep, 4.0, solvationBoxShape, cube,
1211 solvationSolventModel, tip3p,
1212 solventEquilSimulationEquiLengthNVT, 0.1,
1213 solventEquilSimulationEquilLength, 0.2,
1214 solventEquilSimulationProductionLength,0.5,
1215 solventSimulationEquilibrationLength, 1.0,
1216 solventSimulationProductionLength, 10.0,
1217 vacuumEquilSimulationEquilLengthNVT, None,
1218 vacuumEquilSimulationEquilLength, 0.2,
1219 vacuumEquilSimulationProductionLength, 0.5,
1220 vacuumSimulationEquilibrationLength, 0.5,
1221 vacuumSimulationProductionLength, 2.0,
1222 thermoPressure, 0.98692327, thermoTemperature, 298.15"
1223
1224 Author:
1225 Manish Sud(msud@san.rr.com)
1226
1227 See also:
1228 OpenFECalculateAbsoluteBindingFreeEnergy.py,
1229 OpenFECalculatePartialCharges.py, OpenFECalculateRelativeBindingFreeEnergy.py,
1230 OpenFECalculateRelativeHydrationFreeEnergy.py, OpenFEGenerateLigandNetwork.py
1231
1232 Copyright:
1233 Copyright (C) 2026 Manish Sud. All rights reserved.
1234
1235 The functionality available in this script is implemented using OpenFE, an
1236 open source molecuar for alchemical free energy calculations.
1237
1238 This file is part of MayaChemTools.
1239
1240 MayaChemTools is free software; you can redistribute it and/or modify it under
1241 the terms of the GNU Lesser General Public License as published by the Free
1242 Software Foundation; either version 3 of the License, or (at your option) any
1243 later version.
1244
1245 """
1246
1247 if __name__ == "__main__":
1248 main()
