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