1 #!/bin/env python
2 #
3 # File: OpenFECalculateRelativeBindingFreeEnergy.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 numpy as np
38 import pandas as pd
39
40 # OpenFE imports...
41 try:
42 import openfe
43 except ImportError as ErrMsg:
44 sys.stderr.write("\nFailed to import OpenFE related module/package: %s\n" % ErrMsg)
45 sys.stderr.write("Check/update your OpenFE environment and try again.\n\n")
46 sys.exit(1)
47
48 # RDKit imports...
49 try:
50 from rdkit import rdBase
51 except ImportError as ErrMsg:
52 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
53 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
54 sys.exit(1)
55
56 # MayaChemTools imports...
57 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
58 try:
59 from docopt import docopt
60 import MiscUtil
61 import OpenFEUtil
62 except ImportError as ErrMsg:
63 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
64 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
65 sys.exit(1)
66
67 ScriptName = os.path.basename(sys.argv[0])
68 Options = {}
69 OptionsInfo = {}
70
71
72 def main():
73 """Start execution of the script."""
74
75 MiscUtil.PrintInfo(
76 "\n%s (OpenFE v%s; OpenMM v%s; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
77 % (
78 ScriptName,
79 openfe.version("openfe"),
80 openfe.version("openmm"),
81 rdBase.rdkitVersion,
82 MiscUtil.GetMayaChemToolsVersion(),
83 time.asctime(),
84 )
85 )
86
87 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
88
89 # Retrieve command line arguments and options...
90 RetrieveOptions()
91
92 if Options["--list"]:
93 ProcessListOption()
94 else:
95 # Process and validate command line arguments and options...
96 ProcessOptions()
97
98 # Perform actions required by the script...
99 CalculateRelativeBindingFreeEnergy()
100
101 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
102 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
103
104
105 def CalculateRelativeBindingFreeEnergy():
106 """Calculate relative binding free energy."""
107
108 # Process input files...
109 MacroMol, Mols = ProcessInputFiles()
110
111 # Validate molecule names...
112 ValidateMoleculeNames(Mols)
113
114 # Check for miising partial charges...
115 CheckMissingPartialCharges(Mols)
116
117 # Setup atom mapping...
118 MolAToMolBMappings = GenerateAtomMappings(Mols)
119
120 # Initialize RBFE protocols...
121 RBFEProtocol, RBFEProtocolChargeCorrection = InitializeRelativeHybridTopologyProtocol()
122
123 # Initialize solvent...
124 Solvent = InitializeSolventComponent()
125
126 # Setup transformations...
127 MolAToMolBTransformations = SetupTransformations(
128 MolAToMolBMappings, MacroMol, Solvent, RBFEProtocol, RBFEProtocolChargeCorrection
129 )
130
131 # Setup protocol DAGs...
132 MolAToMolBProtocolDAGs = SetupProtocolDAGs(MolAToMolBTransformations)
133
134 # Execute protocol DAGs and gather results...
135 MolAToMolBProtocolResults = ExecuteProtocolDAGsAndGatherResults(MolAToMolBTransformations, MolAToMolBProtocolDAGs)
136
137 # Process protocol results...
138 ProcessProtocolResults(MolAToMolBTransformations, MolAToMolBProtocolResults)
139
140
141 def InitializeRelativeHybridTopologyProtocol():
142 """Initialize relative hybrid toplology protocol."""
143
144 MiscUtil.PrintInfo("\nInitializing relative hybrid topology protocol...")
145
146 RBFESettings = OpenFEUtil.SetupRelativeFreeEnergySettings("-r, --rbfeParams", OptionsInfo["RBFEParams"])
147 RBFEProtocol = OpenFEUtil.InitializeRelativeFreeEngeryHybridTopologyProtocol(RBFESettings)
148
149 RBFESettingsChargeCorrection = OpenFEUtil.SetupRelativeFreeEnergySettings(
150 "-r, --rbfeParams", OptionsInfo["RBFEParams"]
151 )
152 OpenFEUtil.UpdateRelativeFreeEnergySettingsForChargeCorrection(
153 "--rbfeChargeCorrectionParams", OptionsInfo["RBFEChargeCorrectionParams"], RBFESettingsChargeCorrection
154 )
155 RBFEProtocolChargeCorrection = OpenFEUtil.InitializeRelativeFreeEngeryHybridTopologyProtocol(
156 RBFESettingsChargeCorrection
157 )
158
159 return (RBFEProtocol, RBFEProtocolChargeCorrection)
160
161
162 def InitializeSolventComponent():
163 """Initialize solvent component."""
164
165 SolventParams = OptionsInfo["SolventParams"]
166 MiscUtil.PrintInfo(
167 "\nInitializing solvent component (PositiveIon: %s; NegativeIon: %s; Neutralize: %s; IonConcentration: %s)..."
168 % (
169 SolventParams["PositiveIon"],
170 SolventParams["NegativeIon"],
171 SolventParams["Neutralize"],
172 SolventParams["IonConcentration"],
173 )
174 )
175
176 Solvent = OpenFEUtil.InitializeSolventComponent(SolventParams)
177
178 return Solvent
179
180
181 def SetupTransformations(MolAToMolBMappings, MacroMol, Solvent, RBFEProtocol, RBFEProtocolChargeCorrection):
182 """Set up a transformation pair for each mapping."""
183
184 MiscUtil.PrintInfo("\nSetting up alchemical transformations (Count: %s)..." % (len(MolAToMolBMappings) * 2))
185
186 MolAToMolBTransformations = []
187
188 for MolAToMolBMapping in MolAToMolBMappings:
189 MolA = MolAToMolBMapping.componentA
190 MolB = MolAToMolBMapping.componentB
191
192 # Setup chemical systems...
193 MolAComplexSolventSystem = OpenFEUtil.InitializeChemicalSystem(
194 SmallMol=MolA, MacroMol=MacroMol, Solvent=Solvent, Name="%s_Complex_Solvent" % MolA.name
195 )
196 MolASolventSystem = OpenFEUtil.InitializeChemicalSystem(
197 SmallMol=MolA, MacroMol=None, Solvent=Solvent, Name="%s_Solvent" % MolA.name
198 )
199
200 MolBComplexSolventSystem = OpenFEUtil.InitializeChemicalSystem(
201 SmallMol=MolB, MacroMol=MacroMol, Solvent=Solvent, Name="%s_Complex_Solvent" % MolB.name
202 )
203 MolBSolventSystem = OpenFEUtil.InitializeChemicalSystem(
204 SmallMol=MolB, MacroMol=None, Solvent=Solvent, Name="%s_Solvent" % MolB.name
205 )
206
207 # Setup transformation protocol...
208 TransformationProtocol = SetupTransformationProtocol(
209 MolAToMolBMapping, RBFEProtocol, RBFEProtocolChargeCorrection
210 )
211
212 # Setup MolAComplexSolvent to MolBComplexSolvent transformation...
213 TransformationName = "%s_To_%s_Complex_Solvent" % (MolA.name, MolB.name)
214 MolAToMolBComplexSolventTransformation = OpenFEUtil.InitializeTransformation(
215 StateA=MolAComplexSolventSystem,
216 StateB=MolBComplexSolventSystem,
217 Mapping=MolAToMolBMapping,
218 Protocol=TransformationProtocol,
219 Name=TransformationName,
220 Validate=False,
221 )
222 MolAToMolBTransformations.append(MolAToMolBComplexSolventTransformation)
223
224 # Setup MolASolvent to MolBSolvent transformation...
225 TransformationName = "%s_To_%s_Solvent" % (MolA.name, MolB.name)
226 MolAToMolBSolventTransformation = OpenFEUtil.InitializeTransformation(
227 StateA=MolASolventSystem,
228 StateB=MolBSolventSystem,
229 Mapping=MolAToMolBMapping,
230 Protocol=RBFEProtocol,
231 Name=TransformationName,
232 Validate=False,
233 )
234 MolAToMolBTransformations.append(MolAToMolBSolventTransformation)
235
236 # Write out transformatios...
237 WriteTransformations(MolAToMolBTransformations)
238
239 return MolAToMolBTransformations
240
241
242 def SetupTransformationProtocol(MolAToMolBMapping, RBFEProtocol, RBFEProtocolChargeCorrection):
243 """Setup transformation protocol."""
244
245 from openfe.utils import ligand_utils
246
247 ChargeDifference = ligand_utils.get_alchemical_charge_difference(MolAToMolBMapping)
248
249 if ChargeDifference != 0:
250 MolA = MolAToMolBMapping.componentA
251 MolB = MolAToMolBMapping.componentB
252 if OptionsInfo["RBFEChargeCorrection"]:
253 TransformationProtocol = RBFEProtocolChargeCorrection
254 MiscUtil.PrintInfo("")
255 MiscUtil.PrintWarning(
256 'The transformation between molecules %s and %s involves a charge change of %s. The RBFE setting parameters have been automatically updated for "Yes" value of option "--rbfeChargeCorrection" to employ a more expensive set of parameters specified by option "--rbfeChargeCorrectionParams". '
257 % (MolA.name, MolB.name, ChargeDifference)
258 )
259 else:
260 TransformationProtocol = RBFEProtocol
261 MiscUtil.PrintInfo("")
262 MiscUtil.PrintWarning(
263 'The transformation between molecules %s and %s involves a charge change of %s. The RBFE setting parameters have not been automatically updated for "No" value of option "--rbfeChargeCorrection" to employ more expensive set of parameters specified by option "--rbfeChargeCorrectionParams". A word to the wise: You may want to consider sepecifying "Yes" value for option "--rbfeChargeCorrection".'
264 % (MolA.name, MolB.name, ChargeDifference)
265 )
266 else:
267 TransformationProtocol = RBFEProtocol
268
269 return TransformationProtocol
270
271
272 def WriteTransformations(MolAToMolBTransformations):
273 """Write out transformations."""
274
275 TransformationsOutDirPath = pathlib.Path(OptionsInfo["TransformationsOutDirPath"])
276
277 MiscUtil.PrintInfo(
278 "Writing transformations files (Files: *.json; Count: %s; Subdirectory: %s)..."
279 % (len(MolAToMolBTransformations), OptionsInfo["TransformationsOutDir"])
280 )
281
282 for Transformation in MolAToMolBTransformations:
283 TransformationFilePath = TransformationsOutDirPath.joinpath("%s.json" % Transformation.name)
284 Transformation.dump(TransformationFilePath)
285
286
287 def SetupProtocolDAGs(MolAToMolBTransformations):
288 """Setup protocol Directed Acyclic Graphs (DAGs) for each transformation to
289 to perform calculations.
290 """
291
292 MiscUtil.PrintInfo("\nSetting up protocol DAGs (Count: %s)..." % len(MolAToMolBTransformations))
293
294 MolAToMolBProtocolDAGs = []
295 for Transformation in MolAToMolBTransformations:
296 ProtocolDAG = OpenFEUtil.InitializeProtocolDAG(Transformation, Name=Transformation.name)
297 MolAToMolBProtocolDAGs.append(ProtocolDAG)
298
299 return MolAToMolBProtocolDAGs
300
301
302 def ExecuteProtocolDAGsAndGatherResults(MolTransformations, MolProtocolDAGs):
303 """Execute protocol DAGs and gather results."""
304
305 ResultsSharedOutDirPath = OptionsInfo["ResultsOutDirPath"]
306 ResultsScratchOutDirPath = OptionsInfo["ResultsScratchOutDirPath"]
307 ExecuteDAGParams = OptionsInfo["ExecuteDAGParams"]
308
309 MolProtocolResults = OpenFEUtil.ExecuteProtocolDAGsAndGatherResults(
310 MolTransformations,
311 MolProtocolDAGs,
312 ResultsSharedOutDirPath,
313 ResultsScratchOutDirPath,
314 KeepShared=ExecuteDAGParams["KeepShared"],
315 KeepScratch=ExecuteDAGParams["KeepScratch"],
316 NRetries=ExecuteDAGParams["NRetries"],
317 WriteResults=True,
318 )
319
320 return MolProtocolResults
321
322
323 def ProcessProtocolResults(MolAToMolBTransformations, MolAToMolBProtocolResults):
324 """Process protocol results."""
325
326 ResultFileParams = OptionsInfo["ResultFileParams"]
327
328 ResultFile = "%s_RBFE_Results.%s" % (OptionsInfo["OutfilePrefix"], ResultFileParams["Ext"])
329 ResultFilePath = os.path.join(OptionsInfo["OutfileDirPath"], ResultFile)
330 MiscUtil.PrintInfo("\nWriting %s..." % ResultFile)
331
332 Precision = ResultFileParams["Precision"]
333
334 ResultData = []
335 for Index in range(0, len(MolAToMolBProtocolResults), 2):
336 MolAToMolBComplexSolventProtocolResult = MolAToMolBProtocolResults[Index]
337 MolAToMolBSolventProtocolResult = MolAToMolBProtocolResults[Index + 1]
338
339 # Setup mol names using complex solvent transformation. The solvent transformation
340 # also contains the same pair of moleules.
341 MolAToMolBComplexSolventTransformation = MolAToMolBTransformations[Index]
342
343 MolA = MolAToMolBComplexSolventTransformation.stateA.components["ligand"]
344 MolB = MolAToMolBComplexSolventTransformation.stateB.components["ligand"]
345 MolAName = MolA.name
346 MolBName = MolB.name
347
348 if MolAToMolBComplexSolventProtocolResult is None or MolAToMolBSolventProtocolResult is None:
349 DeltaDeltaGBinding = "NA"
350 DeltaDeltaGBindingUncertainty = "NA"
351 else:
352 # Setup binding value without the units...
353 MolAToMolBComplexSolventDeltaG = MolAToMolBComplexSolventProtocolResult.get_estimate()
354 MolAToMolBSolventDeltaG = MolAToMolBSolventProtocolResult.get_estimate()
355
356 DeltaDeltaGBinding = MolAToMolBComplexSolventDeltaG.m - MolAToMolBSolventDeltaG.m
357 DeltaDeltaGBinding = "%.*f" % (Precision, DeltaDeltaGBinding)
358
359 # Setup uncertainty value without the units...
360 MolAToMolBComplexSolventDeltaGUncertainty = MolAToMolBComplexSolventProtocolResult.get_uncertainty()
361 MolAToMolBSolventDeltaGUncertainty = MolAToMolBSolventProtocolResult.get_uncertainty()
362
363 DeltaDeltaGBindingUncertainty = np.sqrt(
364 np.sum(np.square([MolAToMolBComplexSolventDeltaGUncertainty.m, MolAToMolBSolventDeltaGUncertainty.m]))
365 )
366 DeltaDeltaGBindingUncertainty = "%.*f" % (Precision, DeltaDeltaGBindingUncertainty)
367
368 ResultData.append([MolAName, MolBName, DeltaDeltaGBinding, DeltaDeltaGBindingUncertainty])
369
370 ResultDF = pd.DataFrame(
371 ResultData,
372 columns=["MolAName", "MolBName", "DeltaDeltaG (MolA->MolB; RBFE) (kcal/mol)", "Uncertainty (kcal/mol)"],
373 )
374 ResultDF.to_csv(ResultFilePath, sep=ResultFileParams["Delim"], lineterminator="\n", index=False)
375
376
377 def GenerateAtomMappings(Mols):
378 """Generate atom mappings."""
379
380 MiscUtil.PrintInfo("\nChanging directory to %s..." % OptionsInfo["OutfileDir"])
381 os.chdir(OptionsInfo["OutfileDirPath"])
382
383 # Initialize atom mappers...
384 MiscUtil.PrintInfo("\nInitializing atom mappers (%s)..." % " ".join(OptionsInfo["MapperList"]))
385 Mappers = OpenFEUtil.InitializeAtomMappers(OptionsInfo["MapperList"], OptionsInfo["MapperParams"])
386
387 MiscUtil.PrintInfo("\nInitializing atom mapper scorer (%s)..." % OptionsInfo["MapperScorer"])
388 MapperScorer = OpenFEUtil.InitializeAtomMapperScorer(OptionsInfo["MapperScorer"])
389
390 MolAToMolBMappings = None
391 if OptionsInfo["MoleculePairsMode"]:
392 MolAToMolBMappings = GenerateAtomMappingsForMoleculePairs(Mols, Mappers, MapperScorer)
393 elif OptionsInfo["MoleculeNetworkMode"]:
394 MolAToMolBMappings = GenerateAtomMappingForMoleculeNetwork(Mols, Mappers, MapperScorer)
395
396 if MolAToMolBMappings is None or len(MolAToMolBMappings) == 0:
397 MiscUtil.PrintError("Failed to generate atom mappings for small molecules.")
398
399 return MolAToMolBMappings
400
401
402 def GenerateAtomMappingsForMoleculePairs(Mols, Mappers, MapperScorer):
403 """Generate atom mappings for molecule pairs."""
404
405 MiscUtil.PrintInfo(
406 "\nGenerating atom mappings (%s: %d)..." % (OptionsInfo["Mode"], (len(OptionsInfo["MoleculePairsMolList"]) / 2))
407 )
408
409 # Generate atom mappings...
410 MolAToMolBMappings = OpenFEUtil.SuggestAtomMappingsForMoleculePairs(
411 OptionsInfo["MoleculePairsMolList"], Mappers, MapperScorer
412 )
413
414 # Write out image files...
415 WriteMoleculePairsOutputFiles(MolAToMolBMappings)
416
417 return MolAToMolBMappings
418
419
420 def WriteMoleculePairsOutputFiles(MolAToMolBMappings):
421 """Write mapping image output files for molecule pairs."""
422
423 if len(MolAToMolBMappings):
424 MiscUtil.PrintInfo(
425 "Writing molecule pairs output files (Files: <MolName1>_To_<MolName2>_*.png; Count: %s; Subdirectory: %s)..."
426 % (len(MolAToMolBMappings), OptionsInfo["PairImagesOutfileDir"])
427 )
428
429 for Mapping in MolAToMolBMappings:
430 PairOutfilePath = SetupMappingImageFilePath("Molecule_Pair", Mapping, OptionsInfo["PairImagesOutfileDirPath"])
431 OpenFEUtil.WriteMappingImageFile(Mapping, PairOutfilePath)
432
433
434 def GenerateAtomMappingForMoleculeNetwork(Mols, Mappers, MapperScorer):
435 """Setup atom mapping for molecule network."""
436
437 MiscUtil.PrintInfo("\nGenerating atom mappings (%s)..." % OptionsInfo["Mode"])
438
439 # Generate network...
440 LigandNetwork = OpenFEUtil.GenerateLigandNetwork(
441 Mols, OptionsInfo["Network"], OptionsInfo["NetworkParams"], Mappers, MapperScorer
442 )
443
444 # Write out network output files...
445 WriteMoleculeNetworkOutputFiles(LigandNetwork)
446
447 # Setup mappings...
448 MolAToMolBMappings = [Edge for Edge in LigandNetwork.edges]
449
450 return MolAToMolBMappings
451
452
453 def WriteMoleculeNetworkOutputFiles(LigandNetwork):
454 """Write out network output files."""
455
456 NetworkName = OptionsInfo["Network"]
457 MiscUtil.PrintInfo("\nGenerating ligand network (%s)..." % NetworkName)
458
459 # Write out ligand network graphml and image files...
460 NetworkOutfilePrefix = "%s_Network_%s_Mapper_%s" % (OptionsInfo["OutfilePrefix"], NetworkName, SetupMapperLabel())
461 GraphMLOutfile = "%s.graphml" % NetworkOutfilePrefix
462 ImageOutfile = "%s.%s" % (NetworkOutfilePrefix, OptionsInfo["NetworkParams"]["OutputNetworkFormat"])
463
464 GraphMLOutfilePath = os.path.join(OptionsInfo["OutfileDirPath"], GraphMLOutfile)
465 ImageOutfilePath = os.path.join(OptionsInfo["OutfileDirPath"], ImageOutfile)
466
467 MiscUtil.PrintInfo("Writing %s..." % GraphMLOutfile)
468 OpenFEUtil.WriteLigandNetworkGraphMLFile(LigandNetwork, GraphMLOutfilePath)
469
470 MiscUtil.PrintInfo("Writing %s..." % ImageOutfile)
471 OpenFEUtil.WriteLigandNetworkImageFile(LigandNetwork, ImageOutfilePath)
472
473 # Write out image files for edges...
474 NetworkEdges = [Edge for Edge in LigandNetwork.edges]
475 if OptionsInfo["NetworkParams"]["OutputEdges"]:
476 if len(NetworkEdges):
477 MiscUtil.PrintInfo(
478 "Writing edge output files (Files: <MolName1>_To_<MolName2>_*.png; Count: %s; Subdirectory: %s)..."
479 % (len(NetworkEdges), OptionsInfo["EdgeImagesOutfileDir"])
480 )
481
482 for Edge in NetworkEdges:
483 EdgeOutfilePath = SetupMappingImageFilePath("Network", Edge, OptionsInfo["EdgeImagesOutfileDirPath"])
484 OpenFEUtil.WriteMappingImageFile(Edge, EdgeOutfilePath)
485
486
487 def SetupMappingImageFilePath(ModeLabel, Mapping, OutfileDirPath):
488 """Setup mapping image file path."""
489
490 Outfile = "%s_To_%s_%s_Mapper_%s.png" % (
491 Mapping.componentA.name,
492 Mapping.componentB.name,
493 ModeLabel,
494 SetupMapperLabel(),
495 )
496 Outfile = re.sub(" ", "_", Outfile)
497
498 OutfilePath = os.path.join(OutfileDirPath, Outfile)
499
500 return OutfilePath
501
502
503 def SetupMapperLabel():
504 """Setup mapper label."""
505
506 return "_".join(OptionsInfo["MapperList"])
507
508
509 def ProcessInputFiles():
510 """Process input files."""
511
512 # Read PDB file...
513 MiscUtil.PrintInfo("\nReading PDB file %s..." % OptionsInfo["Infile"])
514 MacroMol = OpenFEUtil.ReadPDBFile(OptionsInfo["InfilePath"], Name=OptionsInfo["InfileRoot"])
515
516 # Read small molecule input file...
517 MiscUtil.PrintInfo("\nReading small molecule file %s..." % OptionsInfo["SmallMolFile"])
518 Mols, MolCount, ValidMolCount = OpenFEUtil.ReadAndValidateMolecules(
519 OptionsInfo["SmallMolFilePath"], **OptionsInfo["SmallMolFileParams"]
520 )
521
522 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
523 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
524 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
525
526 if ValidMolCount == 0:
527 MiscUtil.PrintInfo("")
528 MiscUtil.PrintError("No valid molecules found in small molecule input file.\n")
529
530 if ValidMolCount < 2:
531 MiscUtil.PrintInfo("")
532 MiscUtil.PrintError("Small molecule Input file must contain at least 2 molecules.\n")
533
534 return (MacroMol, Mols)
535
536
537 def ValidateMoleculeNames(Mols):
538 """Validate molecule names."""
539
540 if OptionsInfo["MoleculePairsMode"]:
541 OptionsInfo["MoleculePairsMolList"] = OpenFEUtil.ProcessMoleculePairs(Mols, OptionsInfo["MoleculePairsList"])
542 elif OptionsInfo["MoleculeNetworkMode"]:
543 if OptionsInfo["RadialNetworkStatus"]:
544 OptionsInfo["NetworkParams"]["RadialCentralLigandMol"] = OpenFEUtil.ProcessRadialCentralLigandName(
545 Mols, OptionsInfo["NetworkParams"]["RadialCentralLigand"]
546 )
547
548
549 def CheckMissingPartialCharges(Mols):
550 """Check missing partial charges for small molecules."""
551
552 MiscUtil.PrintInfo("\nChecking missing partial charges for small molecules...")
553
554 MissingChargesMolCount = OpenFEUtil.GetMissingPartialChargesMolCount(Mols)
555 MiscUtil.PrintInfo("Number of molecules with missing partial charges: %s" % MissingChargesMolCount)
556
557 if MissingChargesMolCount == 0:
558 return
559
560 if re.match("^Stop$", OptionsInfo["MissingChargeMode"], re.I):
561 MiscUtil.PrintInfo("")
562 MiscUtil.PrintError(
563 'The small molecule input file contains molecules with missing partial charges. The execution of the script has been terminated for "Stop" value of "--missingChargedMode" option. You may continue the execution of the script by specifying "Calculate" value for "--missingChargedMode" option.\n\nThe missing charges will be automatically calculated by OpenFE RelativeHybridTopologyProtocol module during the calculation of RBFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--rbfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate RBFE.\n'
564 )
565 else:
566 MiscUtil.PrintInfo("")
567 MiscUtil.PrintWarning(
568 'The small molecule input file contains molecules with missing partial charges. The missing charges will be automatically calculated by OpenFE RelativeHybridTopologyProtocol module during the calculation of RBFE. You may control the calculation of partial charges by specifying values for partialCharge* parameters using "--rbfeParams" option. Alternatively, you may employ the OpenFECalculatePartialCharges.py script to calculate partial charges and use the small molecule input file containing charges to calculate RBFE.\n'
569 )
570
571
572 def ProcessOutfilePrefixOption():
573 """Process outfile prefix option."""
574
575 OutfilePrefix = Options["--outfilePrefix"]
576
577 if re.match("^auto$", OutfilePrefix, re.I):
578 OutfilePrefix = OptionsInfo["SmallMolFileRoot"]
579
580 OptionsInfo["OutfilePrefix"] = OutfilePrefix
581
582
583 def ProcessOutfileDirOption():
584 """Process outfile directory Option."""
585
586 # Setup output directory...
587 OutfileDir = Options["--outfileDir"]
588 OutfileDirPath = os.path.abspath(OutfileDir)
589 if not os.path.exists(OutfileDir):
590 MiscUtil.PrintInfo("\nCreating output directory %s..." % (OutfileDir))
591 os.mkdir(OutfileDirPath)
592 OptionsInfo["OutfileDir"] = OutfileDir
593 OptionsInfo["OutfileDirPath"] = OutfileDirPath
594
595 # Setup a images subdirectory for a network...
596 EdgeImagesOutfileDir = "NetworkEdgeImages"
597 EdgeImagesOutfileDirPath = os.path.join(OptionsInfo["OutfileDirPath"], EdgeImagesOutfileDir)
598 if OptionsInfo["MoleculeNetworkMode"] and OptionsInfo["NetworkParams"]["OutputEdges"]:
599 if not os.path.exists(EdgeImagesOutfileDirPath):
600 os.mkdir(EdgeImagesOutfileDirPath)
601 OptionsInfo["EdgeImagesOutfileDir"] = EdgeImagesOutfileDir
602 OptionsInfo["EdgeImagesOutfileDirPath"] = EdgeImagesOutfileDirPath
603
604 # Setup a images subdirectory for molecule pairs...
605 PairImagesOutfileDir = "MoleculePairImages"
606 PairImagesOutfileDirPath = os.path.join(OptionsInfo["OutfileDirPath"], PairImagesOutfileDir)
607 if OptionsInfo["MoleculePairsMode"]:
608 if not os.path.exists(PairImagesOutfileDirPath):
609 os.mkdir(PairImagesOutfileDirPath)
610 OptionsInfo["PairImagesOutfileDir"] = PairImagesOutfileDir
611 OptionsInfo["PairImagesOutfileDirPath"] = PairImagesOutfileDirPath
612
613 # Setup a transformations subdirectory...
614 TransformationsOutDir = "Transformations"
615 TransformationsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], TransformationsOutDir)
616 if not os.path.exists(TransformationsOutDirPath):
617 os.mkdir(TransformationsOutDirPath)
618 OptionsInfo["TransformationsOutDir"] = TransformationsOutDir
619 OptionsInfo["TransformationsOutDirPath"] = TransformationsOutDirPath
620
621 # Setup a results subdirectory...
622 ResultsOutDir = "Results"
623 ResultsOutDirPath = os.path.join(OptionsInfo["OutfileDirPath"], ResultsOutDir)
624 if not os.path.exists(ResultsOutDirPath):
625 os.mkdir(ResultsOutDirPath)
626 OptionsInfo["ResultsOutDir"] = ResultsOutDir
627 OptionsInfo["ResultsOutDirPath"] = ResultsOutDirPath
628
629 # Use results subdirectory for scratch results...
630 OptionsInfo["ResultsScratchOutDir"] = ResultsOutDir
631 OptionsInfo["ResultsScratchOutDirPath"] = ResultsOutDirPath
632
633
634 def ProcessListOption():
635 """Process list protocol settings option."""
636
637 RBFESettings = openfe.protocols.openmm_rfe.RelativeHybridTopologyProtocol.default_settings()
638
639 MiscUtil.PrintInfo("\nListing RBFE settings...")
640 OpenFEUtil.ListOpenFESettings(RBFESettings)
641
642
643 def ConfigureLogging():
644 """Configure logging."""
645
646 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
647
648 if re.match("^Error$", OptionsInfo["LoggingLevel"], re.I):
649 LoggingLevel = logging.ERROR
650 elif re.match("^Warning$", OptionsInfo["LoggingLevel"], re.I):
651 LoggingLevel = logging.WARNING
652 else:
653 LoggingLevel = logging.INFO
654
655 logging.basicConfig(format="%(levelname)s: %(message)s", level=LoggingLevel)
656
657 # Turn warnings issued by warnings.warn() into log message to avoid display
658 # of a stack trace...
659 logging.captureWarnings(True)
660
661
662 def ProcessOptions():
663 """Process and validate command line arguments and options."""
664
665 MiscUtil.PrintInfo("Processing options...")
666
667 # Validate options...
668 ValidateOptions()
669
670 # Configure logging...
671 ConfigureLogging()
672
673 OptionsInfo["Infile"] = Options["--infile"]
674 OptionsInfo["InfilePath"] = os.path.abspath(OptionsInfo["Infile"])
675 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
676 OptionsInfo["InfileRoot"] = FileName
677
678 OptionsInfo["SmallMolFile"] = Options["--smallMolFile"]
679 OptionsInfo["SmallMolFilePath"] = os.path.abspath(OptionsInfo["SmallMolFile"])
680 FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["SmallMolFile"])
681 OptionsInfo["SmallMolFileRoot"] = FileName
682
683 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
684 OptionsInfo["SmallMolFileParams"] = MiscUtil.ProcessOptionInfileParameters(
685 "--smallMolFileParams",
686 Options["--smallMolFileParams"],
687 InfileName=Options["--smallMolFile"],
688 ParamsDefaultInfo=ParamsDefaultInfoOverride,
689 )
690
691 OptionsInfo["ExecuteDAGParams"] = OpenFEUtil.ProcessOptionOpenFEExecuteDAGParameters(
692 "--executeDAGParams", Options["--executeDAGParams"]
693 )
694
695 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"]
696
697 OptionsInfo["MapperList"] = OpenFEUtil.ProcessOptionOpenFEMapper("-m, --mapper", Options["--mapper"])
698 OptionsInfo["MapperParams"] = OpenFEUtil.ProcessOptionOpenFEMapperParameters(
699 "-m, --mapperParams", Options["--mapperParams"]
700 )
701 OptionsInfo["MapperScorer"] = Options["--mapperScorer"]
702
703 OptionsInfo["Mode"] = OpenFEUtil.ProcessOptionOpenFERelativeFreeEnergyMode("-m, --mode", Options["--mode"])
704 OptionsInfo["MoleculePairsMode"] = True if re.match("^MoleculePairs$", OptionsInfo["Mode"], re.I) else False
705 OptionsInfo["MoleculeNetworkMode"] = True if re.match("^MoleculeNetwork$", OptionsInfo["Mode"], re.I) else False
706
707 OptionsInfo["MissingChargeMode"] = OpenFEUtil.ProcessOptionOpenFEMissingChargeMode(
708 "--missingChargeMode", Options["--missingChargeMode"]
709 )
710
711 OptionsInfo["Network"] = OpenFEUtil.ProcessOptionOpenFENetwork("-n, --network", Options["--network"])
712 OptionsInfo["RadialNetworkStatus"] = True if re.match("^Radial$", OptionsInfo["Network"], re.I) else False
713
714 ParamsDefaultInfoOverride = {"OutputEdges": True}
715 OptionsInfo["NetworkParams"] = OpenFEUtil.ProcessOptionOpenFENetworkParameters(
716 "--networkParams",
717 Options["--networkParams"],
718 RadialNetworkStatus=OptionsInfo["RadialNetworkStatus"],
719 ParamsDefaultInfo=ParamsDefaultInfoOverride,
720 )
721
722 OptionsInfo["MoleculePairs"] = Options["--moleculePairs"]
723 OptionsInfo["MoleculePairsList"] = OpenFEUtil.ProcessOptionOpenFEMoleculePairs(
724 "--moleculePairs", Options["--moleculePairs"]
725 )
726 OptionsInfo["MoleculePairsMolList"] = None
727
728 OptionsInfo["ResultFileParams"] = OpenFEUtil.ProcessOptionOpenFEResultFileParameters(
729 "--resultFileParams", Options["--resultFileParams"]
730 )
731
732 ParamsDefaultInfoOverride = {"EngineComputePlatform": "CPU"}
733 OptionsInfo["RBFEParams"] = OpenFEUtil.ProcessOptionOpenFERelativeFreeEnergyParameters(
734 "--rbfeParams", Options["--rbfeParams"], ParamsDefaultInfo=ParamsDefaultInfoOverride
735 )
736
737 OptionsInfo["RBFEChargeCorrection"] = True if re.match("^yes$", Options["--rbfeChargeCorrection"]) else False
738 OptionsInfo["RBFEChargeCorrectionParams"] = (
739 OpenFEUtil.ProcessOptionOpenFERelativeFreeEnergyChargeCorrectionParameters(
740 "--rbfeChargeCorrectionParams", Options["--rbfeChargeCorrectionParams"]
741 )
742 )
743
744 OptionsInfo["SolventParams"] = OpenFEUtil.ProcessOptionOpenFESolventParameters(
745 "--solventParams", Options["--solventParams"]
746 )
747
748 ProcessOutfilePrefixOption()
749 ProcessOutfileDirOption()
750
751 OptionsInfo["Overwrite"] = Options["--overwrite"]
752
753 # Track top level working directory...
754 OptionsInfo["TopWorkingDir"] = os.getcwd()
755
756
757 def RetrieveOptions():
758 """Retrieve command line arguments and options."""
759
760 # Get options...
761 global Options
762 Options = docopt(_docoptUsage_)
763
764 # Set current working directory to the specified directory...
765 WorkingDir = Options["--workingdir"]
766 if WorkingDir:
767 os.chdir(WorkingDir)
768
769 # Handle examples option...
770 if "--examples" in Options and Options["--examples"]:
771 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
772 sys.exit(0)
773
774
775 def ValidateOptions():
776 """Validate option values."""
777
778 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
779 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
780
781 MiscUtil.ValidateOptionFilePath("-s, --smallMolFile", Options["--smallMolFile"])
782 MiscUtil.ValidateOptionFileExt("-s, --smallMolFile", Options["--smallMolFile"], "sdf sd")
783
784 MiscUtil.ValidateOptionDirPath("-o, --outfileDir", Options["--outfileDir"])
785 MiscUtil.ValidateOptionsOutputDirOverwrite(
786 "-o, --outfileDir", Options["--outfileDir"], "--overwrite", Options["--overwrite"]
787 )
788
789 MiscUtil.ValidateOptionTextValue("--loggingLevel", Options["--loggingLevel"], "Info Warning Error")
790
791 for Mapper in Options["--mapper"].split(","):
792 Mapper = Mapper.strip()
793 MiscUtil.ValidateOptionTextValue("--mapper", Mapper, "LOMAP Kartograf")
794
795 MiscUtil.ValidateOptionTextValue("--mapperScorer", Options["--mapperScorer"], "LOMAP")
796
797 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "MoleculePairs MoleculeNetwork")
798 MiscUtil.ValidateOptionTextValue("--missingChargeMode", Options["--missingChargeMode"], "Calculate Stop")
799
800 MoleculePairs = Options["--moleculePairs"]
801 if not re.match("^auto$", MoleculePairs, re.I):
802 MoleculePairsList = MoleculePairs.split(",")
803 if len(MoleculePairsList) % 2:
804 MiscUtil.PrintError(
805 'The number of comma delimited values, %d, specified using "--moleculePairs" option must be an even number.'
806 % (len(MoleculePairsList))
807 )
808
809 MiscUtil.ValidateOptionTextValue("-n, --network", Options["--network"], "LOMAP MinimalSpanning Radial")
810
811 MiscUtil.ValidateOptionTextValue("--rbfeChargeCorrection", Options["--rbfeChargeCorrection"], "yes no")
812
813
814 # Setup a usage string for docopt...
815 _docoptUsage_ = """
816 OpenFECalculateRelativeBindingFreeEnergy.py - Calculate relative binding free energy
817
818 Usage:
819 OpenFECalculateRelativeBindingFreeEnergy.py [--executeDAGParams <Name,Value,..>] [--loggingLevel <Info, Warning or Error>]
820 [--mapper <mapper1, mapper2,...>] [--mapperParams <Name,Value,..>] [--mapperScorer <LOMAP>]
821 [--mode <MoleculePairs or MoleculeNetwork>] [--missingChargeMode <Calculate or Stop>]
822 [--moleculePairs <MolName1,MolName2,..>] [--network <text>] [--networkParams <Name,Value,..>]
823 [--outfilePrefix <text>] [--overwrite] [--rbfeParams <Name,Value,...>] [--rbfeChargeCorrection <yes or no>]
824 [--rbfeChargeCorrectionParams <Name,Value,...>] [--resultFileParams <Name,Value,..>]
825 [--solventParams <Name,Value,...>] [--smallMolFileParams <Name,Value,...> ]
826 [-w <dir>] -i <infile> -s <smallmolfile> -o <outifiledir>
827 OpenFECalculateRelativeBindingFreeEnergy.py -l | --list
828 OpenFECalculateRelativeBindingFreeEnergy.py -h | --help | -e | --examples
829
830 Description:
831 Calculate Relative Binding Free Energy (RBFE) for a pair of molecules in a
832 small molecule input file. You may calculate RBFEs for specific pairs of
833 molecules or all molecule pairs corresponding to edges in a molecule network.
834
835 The input file must contain a macromolecule already prepared for simulation.
836 The preparation of the macromolecule for a simulation generally involves the
837 following tasks: identification and replacement of non-standard residues;
838 addition of missing residues; addition of missing heavy atoms; addition of
839 missing hydrogens.
840
841 In addition, the small molecule input file must contain molecules already
842 prepared for simulation. It must contain appropriate 3D coordinates relative
843 to the macromolecule along with no missing hydrogens.
844
845 The MD simulation workflow, employed for the calculation of RBFEs, involves the
846 following steps: initial minimization; NVT equilibration; NPT equilibration;
847 production NPT. The MD simulation protocol is repeated 3 times for each pair
848 pair of transformations, MolAToMolBComplexSolvent and MolAToMolBSolvent,
849 and the results are analyzed to estimate RBFEs. The default time and step size
850 settings for the MD protocol are shown below:
851
852 Protocol repeats, 3
853
854 Time step size: 4.0 femtosecond
855
856 Max minimization steps: 5,000
857
858 NVT equilibration length: 1.0 nanosecond
859 NPT equilibration length: 1.0 nanosecond
860 NPT production length: 5.0 nanosecond
861
862 Each MD simulation, by default, may run for 7 nanosecond, for a total of 21
863 nanosecond to repeat it 3 times. The total MD simulation time for each pair
864 of transformations, MolAToMolBComplexSolvent and MolAToMolBSolvent,
865 may correspond to more than 42 nanoseconds.
866
867 The supported macromolecule input file formats are: PDB (.pdb) and
868 CIF (.cif)
869
870 The supported small molecule input file format are : SD (.sdf, .sd)
871
872 Possible outfile prefix:
873
874 <OutfilePrefix> or <SmallMolFileRoot>
875
876 Possible output directories:
877
878 <OutfileDir>
879
880 <OutfileDir>/MoleculePairImages [ MoleculeNetwork mode ]
881 <OutfileDir>/NetworkEdgeImages [ MoleculePairs mode]
882
883 <OutfileDir>/Transformations
884 <OutfileDir>/Results
885
886 Possible output files and directories under <OutfileDir>:
887
888 <OutfilePrefix>_RBFE_Results.<csv or tsv>
889
890 MoleculePairImages/<MolAName>_To_<MolBName>_Molecule_Pair*.png
891 ... ... ...
892
893 <OutfilePrefix>_Network*.graphml
894 <OutfilePrefix>_Network*.svg
895 NetworkEdgeImages/<MolAName>_To_<MolBName>_Network*.png
896 ... ... ...
897
898 Transformations/<MolAName>_To_<MolBName>_Complex_Solvent.json
899 Transformations/<MolAName>_To_<MolBName>_Solvent.json
900 ... ... ...
901
902 Results/shared_RelativeHybridTopologyProtocolUnit-*/
903 Results/scratch_RelativeHybridTopologyProtocolUnit-*/
904 ... ... ...
905
906 Options:
907 -e, --examples
908 Print examples.
909 --executeDAGParams <Name,Value,..> [default: auto]
910 A comma delimited list of parameter name and value pairs for executing
911 protocol DAGs (Directed Acyclic Graph) to run RBFE calculations.
912
913 The supported parameter names along with their default values are
914 are shown below:
915
916 keepShared, yes [ Possible values: yes or no ]
917 keepScratch, no [ Possible values: yes or no ]
918 nRetries, 2 [ Possible values: >= 0. A value of 0 implies only
919 1 try. ]
920
921 A brief description of parameters is provided below:
922
923 keepShared: Keep shared directories after the execution of DAG.
924 keepScratch: Keep scratch directories after the execution of DAG.
925 nRetries: Number of times to attempt the execution.
926
927 -h, --help
928 Print this help message.
929 -i, --infile <infile>
930 Input file name containing a macromolecule.
931 -l, --list
932 List default RBFE protocol settings provided by OpenFE module
933 RelativeHybridTopologyProtocol.
934 --loggingLevel <Info, Warning or Error> [default: Error]
935 Logging level to configure the 'root logger' via logging.basicConfig()
936 function. The default logging level is changed from 'logging.INFO' to
937 'logging.ERROR'. Otherwise, OpenFE and its associated modules
938 may generate a lot of informational messages.
939 --mapper <mapper1, mapper2> [default: LOMAP]
940 A comma delimited names of atom mappers for generating atom mapping
941 corresponding to molecule pairs or edges in a molecule network. Possible
942 values: LOMAP [ Lead Optimization MAPer; Ref 176 ] or Kartograf [ Ref 177 ].
943 You may specify multiple mappers for generating mapping between pair of
944 molecules. All specified mappers are employed to identify the highest
945 scoring mapping between a pair of molecules.
946 --mapperParams <Name,Value,..> [default: auto]
947 A comma delimited list of parameter name and value pairs for atom mappers
948 employed to generate mapping between molecule pairs or edges in a molecule
949 network.
950
951 The supported parameter names along with their default values are
952 are shown below:
953
954 lomapTime, 20, [ Units: seconds ]
955 lomapThreeD, yes [ Possible values: yes or no ]
956 lomapMax3D, 1.0 [ Units: Angstrom ]
957 lomapElementChange, yes [ Possible values: yes or no]
958 lomapSeed, None [ Possible value: A string. An empty string causes
959 MCS search to start from scratch ]
960 lomapShift, no [ Possible values: yes or no]
961
962 kartografAtomMaxDistance, 0.95 [ Units: Angstrom ]
963 kartografAtomMapHydrogens, yes [ Possible values: yes or no ]
964 kartografMapHydrogensOnHydrogensOnly, No [ Possible values: yes or
965 no ]
966 kartografMapExactRingMatchesOnly, yes [ Possible values: yes or no ]
967 kartografAllowPartialFusedRings, yes [ Possible values: yes or no ]
968
969 A brief description of parameters is provided below:
970
971 lomapTime: Time out for MCS algorithm.
972 lomapThreeD: Use atom positions to prune symmetric mappings.
973 lomapMax3D: Forbid mapping between atoms with distance more than
974 specified value.
975 lomapElementChange: Allow mappings that change an atom element.
976 lomapSeed: An Empty SMARTS string causes MCS search to start from
977 scratch.
978 lomapShift: Keep pre-aligned atom positions for 3D position checks.
979
980 kartografAtomMaxDistance: Geometric criteria for two atoms
981 corresponding to maximum distance between them.
982 kartografAtomMapHydrogens: Map hydrogens.
983 kartografMapHydrogensOnHydrogensOnly: Map hydrogens only on
984 hydrogens.
985 kartografMapExactRingMatchesOnly: Map rings with only matching ring
986 size and bond orders. In addition, ring breaking is not
987 permitted.
988 kartografAllowPartialFusedRings: Allow mapping of partially fused
989 rings.
990
991 --mapperScorer <LOMAP> [default: LOMAP]
992 Atom mapper scorer to use for scoring mapping between molecule pairs or
993 edges in a molecule network. Possible value: LOMAP. The atom scorer is
994 not used during the generation of MinimalSpanning network.
995 -m, --mode <MoleculePairs or MoleculeNetwork> [default: MoleculePairs]
996 Calculate RBFEs for specified pairs of molecules or all molecule pairs
997 corresponding to edges in a molecule network.
998 --missingChargeMode <Calculate or Stop> [default: Stop]
999 Calculate missing partial charges for molecules before running RBFE
1000 calculations or terminate the execution of the script. The missing
1001 partial charges will be automatically calculated by OpenFE module
1002 RelativeHybridTopologyProtocol during the calculation of RBFE. You
1003 may control the calculation of partial charges by specifying values for
1004 partialCharge* parameters using '--rbfeParams' option.
1005 --moleculePairs <MolName1,MolName2,..> [default: auto]
1006 A comma delimited list of molecule name pairs for calculating RBFEs.
1007 Default: the names of the first and second molecule in small molecule
1008 input file. This option is only used during 'MoleculePairs' value for
1009 '-m, --mode' option.
1010 -n, --network <text> [default: MinimalSpanning]
1011 Name of a molecule network to generate for calculating RBFEs. Possible
1012 values: LOMAP, MinimalSpanning or Radial. This option is only used during
1013 'MoleculeNetwork' value for '-m, --mode' option.
1014 --networkParams <Name,Value,..> [default: auto]
1015 A comma delimited list of parameter name and value pairs for generating
1016 a molecule network.
1017
1018 The supported parameter names along with their default values are
1019 are shown below:
1020
1021 lomapDistanceCutoff, 0.4
1022 lomapMaxPathLength, 6
1023 lomapRequireCycleCovering, yes [ Possible values: yes or no ]
1024
1025 minimalSpanningProgress, no [ Possible values: yes or no ]
1026
1027 radialCentralLigand, None [ Possible values: Valid ligand name ]
1028
1029 outputEdges, no [ Possible values: yes or no ]
1030 outputNetworkFormat, svg [ Possible values: Any valid format. ]
1031
1032 A brief description of parameters is provided below:
1033
1034 lomapDistanceCutoff: Maximum distance/dissimilarity between two
1035 molecules for an edge to be accepted.
1036 lomapMaxPathLength: Maximum distance between any two molecules in
1037 the resulting network
1038 lomapRequireCycleCovering: Add cycles into the network
1039
1040 minimalSpanningProgress: Show progress using tqdm.
1041
1042 radialCentralLigand: Name of central ligand. A valid ligand name
1043 must be specified to generate a radial molecule network.
1044
1045 outputEdges: Generate PNG image files for all edges in a molecule
1046 network.
1047 outputNetworkFormat: Valid image file format for molecule network.
1048 You must specify a valid format supported by Python module
1049 Matplotlib. For example: PNG (.png), SVG (.svg), PDF (.pdf),
1050 etc. In addition, the graphml file is always generated.
1051
1052 -o, --outfileDir <outfiledir>
1053 Output directory.
1054 --outfilePrefix <text> [default: auto]
1055 Prefix for generating output files under output directory.
1056 --overwrite
1057 Overwrite existing files.
1058 --resultFileParams <Name,Value,..> [default: auto]
1059 A comma delimited list of parameter name and value pairs for writing
1060 calculated RBFEs values to a results file.
1061
1062 The supported parameter names along with their default values are
1063 are shown below:
1064
1065 precision, 4 [ Possible values: > 0 ]
1066 delimiter, comma [ Possible values: comma or tab ]
1067
1068 -r, --rbfeParams <Name,Value,...> [default: auto]
1069 A comma delimited list of parameter name and value pairs for RBFE protocol
1070 settings employed during the calculation of RBFEs.
1071
1072 The default values are automatically updated to match settings provided by
1073 OpenFE module RelativeHybridTopologyProtocol.
1074
1075 You must specify valid OpenFE values for these parameters. An extensive
1076 validation is not performed.
1077
1078 The supported parameter names along with their default values are
1079 are shown below:
1080
1081 protocolRepeats, 3
1082
1083 Alchemical settings:
1084
1085 alchemicalEndstateDispersionCorrection, no [ Possible values:
1086 yes or no ]
1087 alchemicalExplicitChargeCorrection, no [ Possible values:
1088 yes or no ]
1089 alchemicalExplicitChargeCorrectionCutoff, 0.8 [ Units: nanometer ]
1090 alchemicalSoftcoreLJ, Gapsys [ Possible values: Gapsys or Beutler ]
1091 alchemicalSoftcoreAlpha, 0.85
1092 alchemicalTurnOffCoreUniqueExceptions, no [ Possible values:
1093 yes or no ]
1094 alchemicalUseDispersionCorrection, no [ Possible values: yes or no ]
1095
1096 Engine settings:
1097
1098 engineComputePlatform, CPU [ Possible values: CPU, CUDA, OpenCL,
1099 or Reference ]
1100 engineGpuDeviceIndex, None [ Possible values: 0, 0 1, etc. ]
1101
1102 Forcefield settings:
1103
1104 forcefieldConstraints, HBonds [ Possible values: HBonds, AllBonds, or
1105 HAngles ]
1106 forcefields, amber/ff14SB.xml amber/tip3p_standard.xml
1107 amber/tip3p_HFE_multivalent.xml amber/phosaa10.xml
1108 [ Possible values: A space delimited list of valid names. ]
1109 forcefieldHydrogenMass, 3.0 [ Units: amu ]y
1110 forcefieldNonbondedCutoff, 0.9 [ Units: nanometer ]
1111 forcefieldNonbondedMethod, PME [ Possible values: PME or NoCutoff ]
1112 forcefieldRigidWater, yes [ Possible values: yes or no ]
1113 forcefieldSmallMoleculeForcefield, openff-2.1.1 [ Possible value:
1114 A valid forcefield name. ]
1115
1116 Integrator settings:
1117
1118 integratorBarostatFrequency, 25.0 * timestep [ The specified value
1119 is a multiple of integratorTimestep. ]
1120 integratorConstraintTolerance, 1e-06
1121 integratorLangevinCollisionRate, 1.0 [ Units: 1 / picosecond ]
1122 integratorNRestartAttempts, 20
1123 integratorReassignVelocities, no [ Possible values: yes or no ]
1124 integratorRemoveCom, no [ Possible values: yes or no ]
1125 integratorTimestep, 4.0 [ Units: femtosecond ]
1126
1127 Lambda settings:
1128
1129 lambdaFunctions, default [ Possible values: Default, namd, or
1130 quarters ]
1131 lambdaWindows, 11
1132
1133 Output settings:
1134
1135 outputCheckpointInterval, 1.0 [ Units: nanosecond ]
1136 outputCheckpointStorageFilename, checkpoint.chk
1137 outputForcefieldCache, db.json
1138 outputFilename, simulation.nc
1139 outputIndices, not water [ Possible value: Any valid selection. ]
1140 outputStructure, hybrid_system.pdb
1141 outputPositionsWriteFrequency, 100.0 [ Units: picosecond ]
1142 outputVelocitiesWriteFrequency, None [ Possible values: > 0;
1143 Units: picosecond ]
1144
1145 Partial charge settings:
1146
1147 partialChargeNaglModel, None [ Default: Production AM1BCC model for
1148 NAGL; Possible value: Any valid name. ]
1149 partialChargeNumberOfConformers, None [ Possible value: > 0 ]
1150 partialChargeOffToolkitBackend, AmberTools [ Possible values:
1151 AmberTools or RDKit ]
1152 partialChargeMethod, AM1BCC [ Possble values: AM1BCC, Espaloma,
1153 or NAGL ]
1154
1155 Simulation settings:
1156
1157 simulationEarlyTerminationTargetError, 0.0 [ Units:
1158 kilocalorie_per_mole ]
1159 simulationEquilibrationLength, 1.0 [ Units: nanosecond ]
1160 simulationMinimizationSteps, 5000
1161 simulationNReplicas, 11
1162 simulationProductionLength, 5.0 [ Units: nanosecond ]
1163 simulationRealTimeAnalysisInterval, 250.0 [ Units: picosecond ]
1164 simulationRealTimeAnalysisMinimumTime, 500.0 [ Units: picosecond ]
1165 simulationSamplerMethod, repex [ Possible values: repex, sams,
1166 or independent ]
1167 simulationSamsFlatnessCriteria, logZ-flatness [ Possible values:
1168 logZ-flatness, minimum-visits or histogram-flatness ]
1169 simulationSamsGamma0, 1.0
1170 simulationTimePerIteration, 2.5 [ Units: picosecond ]
1171
1172 Solvation settings:
1173
1174 solvationBoxShape, dodecahedron [ Possible values: cube,
1175 dodecahedron, or octahedron ]
1176 solvationBoxSize, None [ Possible value: A triplet of space
1177 X Y Z values; Units: nanometer ]
1178 solvationSolventModel, tip3p [ Possible values: tip3p, spce, tip4pew,
1179 or tip5p ]
1180 solvationSolventPadding, 1.5 [ Units: nanometer ]
1181
1182 Thermo settings:
1183
1184 thermoPh, None [ Possible values: > 0 ]
1185 thermoPressure, 1.0 [ Units: bar ]
1186 thermoRedoxPotential, None [ Possible values: A valid float.
1187 Units: millivolts (mV) ]
1188 thermoTemperature, 298.15 [ Units: kelvin ]
1189
1190 A brief description of parameters, taken from OpenFE documentation, is
1191 provided below:
1192
1193 protocolRepeats: Number of completely independent repeats of the
1194 entire sampling process.
1195
1196 Alchemical settings:
1197
1198 Parameters controlling the creation of the hybrid topology system,
1199 including various parameters ranging from softcore parameters to
1200 whether or not to apply an explicit charge correction for systems
1201 with net charge changes.
1202
1203 alchemicalEndstateDispersionCorrection: Employ extra unsampled
1204 endstate windows for long range correction.
1205 alchemicalExplicitChargeCorrection: Explicitly account for a charge
1206 difference during the alchemical transformation by transforming
1207 a water to a counterion of the opposite charge of the formal
1208 charge difference.
1209 alchemicalExplicitChargeCorrectionCutoff: Minimum distance from the
1210 system solutes from which an alchemical water can be chosen.
1211 alchemicalSoftcoreLJ: Use LJ softcore function as defined by Gapsys
1212 [ Ref 181 ] or Buetler [ Ref 182 ].
1213 alchemicalSoftcoreAlpha: Softcore alpha parameter.
1214 alchemicalTurnOffCoreUniqueExceptions: Turn off interactions for
1215 new exceptions (not just 1,4s) at lambda 0 and old exceptions at
1216 lambda 1 between unique atoms and core atoms.
1217 alchemicalUseDispersionCorrection: Use dispersion correction in the
1218 hybrid topology state.
1219
1220 Engine settings:
1221
1222 Parameters configuring the compute platform used by the OpenMM to
1223 perform the simulation.
1224
1225 engineComputePlatform: Platform to use for running OpenMM MD
1226 calculations.
1227 engineGpuDeviceIndex: Space delimited list of device indices to use
1228 for running OpenMM MD calculations.
1229
1230 Forcefield settings:
1231
1232 Parameters to set up the force field with OpenMM Force Fields,
1233 including the general force fields, the small molecule force field,
1234 the nonbonded method, and the nonbonded cutoff.
1235
1236 forcefieldConstraints: Constraints to use.
1237 forcefields: List of valid forcefield paths for all components
1238 except small molecules.
1239 forcefieldHydrogenMass: Mass to be repartitioned to hydrogens from
1240 neighboring heavy atoms.
1241 forcefieldNonbondedCutoff: Cutoff for short range nonbonded
1242 interactions.
1243 forcefieldNonbondedMethod: Method for treating nonbonded
1244 interactions.
1245 forcefieldRigidWater: Use a rigid water model.
1246 forcefieldSmallMoleculeForcefield: A valid forcefield name to use
1247 for small molecules.
1248
1249 Integrator settings
1250
1251 Parameters controlling the LangevinSplittingDynamicsMove integrator
1252 used for simulation.
1253
1254 integratorBarostatFrequency: Frequency at which volume scaling
1255 changes should be attempted.
1256 integratorConstraintTolerance: Tolerance for constraint solver.
1257 integratorLangevinCollisionRate: Collision frequency.
1258 integratorNRestartAttempts: Number of attempts to restart from
1259 Context in case there are NaNs in the energies after
1260 integration.
1261 integratorReassignVelocities: Reassign velocities from the
1262 Maxwell-Boltzmann distribution at the beginning of each
1263 Monte Carlo move.
1264 integratorRemoveCom: Remove the center of mass motion.
1265 integratorTimestep: Size of the simulation timestep.
1266
1267 Lambda settings:
1268
1269 Lambda protocol parameters, including number of lambda windows and
1270 lambda functions.
1271
1272 lambdaFunctions: Function name to use for alchemical mutation.
1273 lambdaWindows: Number of lambda windows to calculate.
1274
1275 Output settings:
1276
1277 Parameter controlling simulation output, including the frequency to
1278 write a checkpoint file, the selection string for writing selected
1279 coordinates, and the paths to the trajectory and output structure
1280 files.
1281
1282 outputCheckpointInterval: Frequency to write the checkpoint file.
1283 outputCheckpointStorageFilename: Checkpoint filename.
1284 outputForcefieldCache: Filename for caching small molecule residue
1285 templates.
1286 outputFilename: Trajectory filename.
1287 outputIndices: Selection string for selecting coordinates to write.
1288 outputStructure: Hybrid topology structure filename.
1289 outputPositionsWriteFrequency: Frequency for writing positions to
1290 trajectory file.
1291 outputVelocitiesWriteFrequency: Frequency for writing velocities to
1292 trajectory file.
1293
1294 Partial charge settings:
1295
1296 Parameters for automatically assigning missing partial charges to
1297 small molecules, including the partial charge method.
1298
1299 partialChargeNaglModel: Model to use for partial charge assignment.
1300 A value of None implies the use of the latest available
1301 production AM1BCC model.
1302 partialChargeNumberOfConformers: Number of conformers to generate
1303 as part of the partial charge assignment. A value of None
1304 implies the use of the existing conformer.
1305 partialChargeOffToolkitBackend: OpenFF toolkit registry backend to
1306 use for calculating partial charges.
1307 partialChargeMethod: Method to use for calculating partial charges.
1308
1309 Simulation settings:
1310
1311 Parameters controlling the simulation plan and the alchemical
1312 sampler, including the number of minimization steps, lengths of
1313 equilibration and production runs, the sampler method (e.g.
1314 Hamiltonian REPlica EXchange (repex), and the time interval at
1315 which to perform an analysis of the free energies.
1316
1317 simulationEarlyTerminationTargetError: Target error for the real
1318 time analysis measured in kcal/mol. Once the MBAR error of the
1319 free energy is at or below this value, the simulation will be
1320 considered complete. The suggested value of 0.12 has shown to
1321 be effective in both hydration and binding free energy
1322 benchmarks.
1323 simulationEquilibrationLength: Length of the equilibration phase.
1324 The specified value must be divisible by 'integratorTimestep'.
1325 simulationMinimizationSteps: Maximum number of minimization steps
1326 to perform.
1327 simulationNReplicas: Number of replicas to use.
1328 simulationProductionLength: Length of the production phase.
1329 The specified value must be divisible by 'integratorTimestep'.
1330 simulationRealTimeAnalysisInterval: Time interval for performing
1331 analysis of the free energies. At each interval, real time
1332 analysis data will be written to a yaml file named
1333 <outputFileName>_real_time_analysis.yaml. The current error
1334 in the estimate will also be assessed and the simulation will
1335 be terminated when it drops below
1336 'simulationEarlyTerminationTargetError'.
1337 simulationRealTimeAnalysisMinimumTime: Minimum simulation time
1338 after which the real time analysis is performed.
1339 simulationSamplerMethod: Alchemical sampling method to use:
1340 REPEX (Hamiltonian REPlica EXchange), SAMS (Self-Adjusted
1341 Mixture Sampling), or Independent (Independently sampled lambda
1342 windows).
1343 simulationSamsFlatnessCriteria:Method for assessing when to switch
1344 to asymptomatically optimal scheme for SAMS.
1345 simulationSamsGamma0: Initial weight adaptation rate for SAMS.
1346 simulationTimePerIteration: Simulation time between each MCMC move
1347 attempt
1348
1349 Solvation settings:
1350
1351 Solvation parameters for the system, including the solvent model and
1352 the solvent padding.
1353
1354 solvationBoxShape: Shape of the periodic solvent box to create.
1355 solvationBoxSize: Lengths of the unit cell for a solvent box.
1356 solvationSolventModel: Forcefield water model to use during
1357 solvation and defining the model properties.
1358 solvationSolventPadding: Minimum distance from any solute bounding
1359 sphere to the edge of the box.
1360
1361 Thermo settings:
1362
1363 Thermodynamic parameters, including the temperature and the pressure
1364 of the system.
1365
1366 thermoPh: Simulation pH
1367 thermoPressure: Simulation pressure.
1368 thermoRedoxPotential:Simulation redox potential.
1369 thermoTemperature: Simulation temperature.
1370
1371 --rbfeChargeCorrection <yes or no> [default: yes]
1372 Perform automatic charge correction for charge changing transformations
1373 during the calculation of RBFEs. The '--rbfeChargeCorrectionParams' are
1374 used during the automatic charge correction to override the corresponding
1375 values in '-r, --rbfeParams'.
1376 --rbfeChargeCorrectionParams <Name,Value,...> [default: auto]
1377 A comma delimited list of parameter name and value pairs to use for RBFE
1378 protocol settings during explicit charge correction for charge changing
1379 transformation between pair of molecules. These parameters override the
1380 corresponding values in '-r, --rbfeParams'.
1381
1382 The default parameter values for charge changing transformations are based
1383 on the industry benchmarking performed by OpenFE.
1384
1385 The supported parameter names along with their default values are
1386 are shown below:
1387
1388 alchemicalExplicitChargeCorrection, yes [ Possible values: yes or no ]
1389 simulationProductionLength, 20 [ Units: nanosecond ]
1390 simulationNReplicas, 22
1391 lambdaWindows, 22
1392
1393 A brief description of these parameters is available under the corresponding
1394 parameters in the section for '-r, --rbfeParams'.
1395 -s, --smallMolFile <SmallMolFile>
1396 Input file containing small molecules.
1397 --smallMolFileParams <Name,Value,...> [default: auto]
1398 A comma delimited list of parameter name and value pairs for reading
1399 molecules from files. The supported parameter names for different file
1400 formats, along with their default values, are shown below:
1401
1402 SD: removeHydrogens,no,sanitize,yes,strictParsing,yes
1403
1404 --solventParams <Name,Value,...> [default: auto]
1405 A comma delimited list of parameter name and value pairs for solvent
1406 component. You must specify valid OpenFE values. No extensive validation
1407 is performed. These parameters are used in conjunction with solvation*
1408 parameters available through '--rbfeParams' to perform solvation.
1409
1410 The supported parameter names along with their default values are
1411 are shown below:
1412
1413 positiveIon, Na+ [ Possible value: Li+, Na+, K+, Rb+, or Cs+ ]
1414 negativeIon, Cl- [ Possible values: Cl-, Br-, F-, or I- ]
1415 neutralize, yes [ Possible values: yes or no ]
1416 ionConcentration, 0.15 [ Units: molar ]
1417
1418 A brief description of parameters is provided below:
1419
1420 positiveIon, negativeion: Pair of ions used to neutralize and bring
1421 the solvent to required ionic concentration.
1422 neutralize: Neutralize the net charge on the chemical state by the
1423 ions in the solvent component.
1424 ionConcentration: Ionic concentration.
1425
1426 -w, --workingdir <dir>
1427 Location of working directory which defaults to the current directory.
1428
1429 Examples:
1430 The sample protein and ligand files for tyrosine kinase 2 (Tyk2) are
1431 distributed with MayaChemTools and are available in data directory. These
1432 files have been taken from OpenFE distribution for example notebooks. The
1433 AM1BCC partial charges have been calculated for the ligands in SD file to
1434 facilitate calculations. You may review OpenFE tutorial notebooks for the
1435 expected results.
1436
1437 To calculate RBFE for a pair molecules corresponding to the fist and second
1438 molecules in a SD file, performing 3 independent repeats of the entire MD sampling
1439 process to calculate RBFE for a pair of molecules, each MD repeat consisting of
1440 minimization (5,000 steps) followed by NVT and NPT equilibration (1 ns;
1441 250,000 steps) leading to NPT production (5ns; 1,250,000 steps) using a step size
1442 of 4 fs, writing out appropriate trajectory and PDB files for each MD repeat
1443 in Results subdirectory under output directory, generating final results file
1444 along with appropriate graph and image files under output directory, type:
1445
1446 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1447 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFE
1448
1449 To run the first example for calculating RBFE for a specific pair molecules
1450 using CUDA platform on your machine to perform MD simulations and generate
1451 various output files, type:
1452
1453 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1454 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFE -m MoleculePairs
1455 --moleculePairs "lig_ejm_31, lig_ejm_47"
1456 --rbfeParams "engineComputePlatform,CUDA"
1457
1458 To run the second example to see all warning messages produced by OpenFE
1459 modules and write various output files, type;
1460
1461 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1462 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFE -m MoleculePairs
1463 --moleculePairs "lig_ejm_31, lig_ejm_47"
1464 --rbfeParams "engineComputePlatform,CUDA"
1465 --loggingLevel Warning
1466
1467 To run the first example for calculating RBFE for all pairs of molecules
1468 corresponding to edges in a molecule network using CUDA platform on your
1469 machine to perform MD simulations and generate various output files, type:
1470
1471 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1472 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFE -m MoleculeNetwork
1473 --network MinimalSpanning
1474 --rbfeParams "engineComputePlatform,CUDA"
1475
1476 To run the first example for calculating RBFE for a specific pair molecules
1477 using CUDA platform on your machine to perform MD simulations, automatically
1478 calculate missing partial charges for molecules, and generate various output
1479 files, type:
1480
1481 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1482 -s SampleTyk2LigandsNoCharges.sdf -o SampleTyk2LigandsRBFE
1483 -m MoleculePairs --moleculePairs "lig_ejm_31, lig_ejm_47"
1484 --rbfeParams "engineComputePlatform,CUDA"
1485 --missingChargeMode Calculate
1486
1487 To run the second example by specifying explict values for various parametres
1488 and generate various output files, type:
1489
1490 % OpenFECalculateRelativeBindingFreeEnergy.py -i SampleTyk2.pdb
1491 -s SampleTyk2Ligands.sdf -o SampleTyk2LigandsRBFE -m MoleculePairs
1492 --moleculePairs "lig_ejm_31, lig_ejm_47"
1493 --loggingLevel Error
1494 --executeDAGParams "keepShared, yes, nRetries, 2" --mapper LOMAP
1495 --mapperParams "lomapTime, 20, lomapThreeD, yes"
1496 --missingChargeMode Stop --rbfeParams "protocolRepeats,3,
1497 alchemicalSoftcoreLJ, Gapsys, engineComputePlatform,CUDA,
1498 forcefieldConstraints, HBonds, forcefieldHydrogenMass, 3.0,
1499 forcefieldNonbondedMethod, PME, integratorTimestep, 4.0,
1500 lambdaWindows, 11, outputCheckpointInterval, 250.0,
1501 simulationMinimizationSteps, 5000, simulationEquilibrationLength, 1.0,
1502 simulationProductionLength, 5.0, solvationBoxShape, cube,
1503 solvationSolventPadding, 1.2, thermoPressure, 0.98692327,
1504 thermoTemperature, 298.15"
1505 --solventParams "positiveIon, Na+, negativeIon, Cl-"
1506
1507 Author:
1508 Manish Sud(msud@san.rr.com)
1509
1510 See also:
1511 OpenFECalculateAbsoluteBindingFreeEnergy.py,
1512 OpenFECalculateAbsoluteHydrationFreeEnergy.py, OpenFECalculatePartialCharges.py,
1513 OpenFECalculateRelativeHydrationFreeEnergy.py, OpenFEGenerateLigandNetwork.py
1514
1515 Copyright:
1516 Copyright (C) 2026 Manish Sud. All rights reserved.
1517
1518 The functionality available in this script is implemented using OpenFE, an
1519 open source molecuar for alchemical free energy calculations.
1520
1521 This file is part of MayaChemTools.
1522
1523 MayaChemTools is free software; you can redistribute it and/or modify it under
1524 the terms of the GNU Lesser General Public License as published by the Free
1525 Software Foundation; either version 3 of the License, or (at your option) any
1526 later version.
1527
1528 """
1529
1530 if __name__ == "__main__":
1531 main()
