1 #!/bin/env python 2 # 3 # File: OpenFECalculatePartialCharges.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 37 # OpenFE imports... 38 try: 39 import openfe 40 except ImportError as ErrMsg: 41 sys.stderr.write("\nFailed to import OpenFE related module/package: %s\n" % ErrMsg) 42 sys.stderr.write("Check/update your OpenFE environment and try again.\n\n") 43 sys.exit(1) 44 45 # RDKit imports... 46 try: 47 from rdkit import rdBase 48 except ImportError as ErrMsg: 49 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg) 50 sys.stderr.write("Check/update your RDKit environment and try again.\n\n") 51 sys.exit(1) 52 53 # MayaChemTools imports... 54 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python")) 55 try: 56 from docopt import docopt 57 import MiscUtil 58 import RDKitUtil 59 import OpenFEUtil 60 except ImportError as ErrMsg: 61 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg) 62 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n") 63 sys.exit(1) 64 65 ScriptName = os.path.basename(sys.argv[0]) 66 Options = {} 67 OptionsInfo = {} 68 69 70 def main(): 71 """Start execution of the script.""" 72 73 MiscUtil.PrintInfo( 74 "\n%s (OpenFE v%s; OpenMM v%s; RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 75 % ( 76 ScriptName, 77 openfe.version("openfe"), 78 openfe.version("openmm"), 79 rdBase.rdkitVersion, 80 MiscUtil.GetMayaChemToolsVersion(), 81 time.asctime(), 82 ) 83 ) 84 85 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 86 87 # Retrieve command line arguments and options... 88 RetrieveOptions() 89 90 # Process and validate command line arguments and options... 91 ProcessOptions() 92 93 # Perform actions required by the script... 94 CalculatePartialCharges() 95 96 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 97 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 98 99 100 def CalculatePartialCharges(): 101 """Calculate partial charges and write them out to a SD file.""" 102 103 # Process molecules... 104 Mols = ProcessMolecules() 105 106 # Calculate charges... 107 ChargedMols = CalculateCharges(Mols) 108 109 # Write charges... 110 WriteCharges(ChargedMols) 111 112 113 def CalculateCharges(Mols): 114 """Calculate partial charges.""" 115 116 MiscUtil.PrintInfo("\nCalculating partial atomic charges (%s)..." % OptionsInfo["Charge"]) 117 118 ChargedMols, Status = OpenFEUtil.CalculatePartialCharges(Mols, OptionsInfo["Charge"], OptionsInfo["ChargeParams"]) 119 120 if not Status: 121 MiscUtil.PrintError("Failed to calculate partial atomic charges.") 122 123 return ChargedMols 124 125 126 def WriteCharges(ChargedMols): 127 """Write charges.""" 128 129 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"]) 130 if Writer is None: 131 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"]) 132 MiscUtil.PrintInfo("\nGenerating file %s..." % OptionsInfo["Outfile"]) 133 134 RDKitChargedMols = [openfe.SmallMoleculeComponent.to_rdkit(Mol) for Mol in ChargedMols] 135 136 for Mol in RDKitChargedMols: 137 FormatPartialCharges(Mol) 138 Writer.write(Mol) 139 Writer.close() 140 141 142 def FormatPartialCharges(Mol): 143 """Format partial charges.""" 144 145 Precision = OptionsInfo["ChargeParams"]["Precision"] 146 LineSize = OptionsInfo["ChargeParams"]["LineSize"] 147 148 # RDkit uses Chem.CreateAtomDoublePropertyList(rdmol, "PartialCharge") to set 149 # up 'atom.dprop.PartialCharge' employing "\n" as new line delimiter. 150 PropName = OpenFEUtil.GetPartialChargePropName() 151 LineDelim = "\n" 152 153 ChargesString = Mol.GetProp(PropName) 154 155 FormattedChargesLines = [] 156 CurrentLine = None 157 CurrentLineSize = 0 158 159 for ChargesLine in ChargesString.split(LineDelim): 160 for Value in ChargesLine.split(): 161 FormattedValue = "%.*f" % (Precision, float(Value)) 162 FormattedValueSize = len(FormattedValue) 163 164 if (FormattedValueSize + CurrentLineSize + 1) >= LineSize: 165 if CurrentLine is not None: 166 FormattedChargesLines.append(" ".join(CurrentLine)) 167 168 CurrentLine = [FormattedValue] 169 CurrentLineSize = FormattedValueSize 170 else: 171 CurrentLineSize += FormattedValueSize 172 if CurrentLine is None: 173 CurrentLine = [FormattedValue] 174 else: 175 # Increment line size to account for space delimiter.... 176 CurrentLineSize += 1 177 CurrentLine.append(FormattedValue) 178 179 if CurrentLine is not None: 180 FormattedChargesLines.append(" ".join(CurrentLine)) 181 182 FormattedChargesString = LineDelim.join(FormattedChargesLines) 183 Mol.SetProp(PropName, FormattedChargesString) 184 185 186 def ProcessMolecules(): 187 """Process molecules.""" 188 189 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"]) 190 Mols, MolCount, ValidMolCount = OpenFEUtil.ReadAndValidateMolecules( 191 OptionsInfo["InfilePath"], **OptionsInfo["InfileParams"] 192 ) 193 194 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount) 195 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount) 196 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount)) 197 198 if ValidMolCount == 0: 199 MiscUtil.PrintInfo("") 200 MiscUtil.PrintError("No valid molecules found in input file.\n") 201 202 return Mols 203 204 205 def ConfigureLogging(): 206 """Configure logging.""" 207 208 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"] 209 210 if re.match("^Warning$", OptionsInfo["LoggingLevel"], re.I): 211 LoggingLevel = logging.WARNING 212 else: 213 LoggingLevel = logging.INFO 214 215 logging.basicConfig(format="%(levelname)s: %(message)s", level=LoggingLevel) 216 217 218 def ProcessOptions(): 219 """Process and validate command line arguments and options.""" 220 221 MiscUtil.PrintInfo("Processing options...") 222 223 # Validate options... 224 ValidateOptions() 225 226 # Configure logging... 227 ConfigureLogging() 228 229 OptionsInfo["Infile"] = Options["--infile"] 230 OptionsInfo["InfilePath"] = os.path.abspath(OptionsInfo["Infile"]) 231 232 ParamsDefaultInfoOverride = {"RemoveHydrogens": False} 233 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters( 234 "--infileParams", 235 Options["--infileParams"], 236 InfileName=Options["--infile"], 237 ParamsDefaultInfo=ParamsDefaultInfoOverride, 238 ) 239 240 OptionsInfo["Outfile"] = Options["--outfile"] 241 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters( 242 "--outfileParams", Options["--outfileParams"] 243 ) 244 245 OptionsInfo["Charge"] = OpenFEUtil.ProcessOptionOpenFECharge("-c, --charge", Options["--charge"]) 246 OptionsInfo["ChargeParams"] = OpenFEUtil.ProcessOptionOpenFEChargeParameters( 247 "--chargeParams", Options["--chargeParams"], OptionsInfo["Charge"] 248 ) 249 250 OptionsInfo["LoggingLevel"] = Options["--loggingLevel"] 251 252 OptionsInfo["Overwrite"] = Options["--overwrite"] 253 254 255 def RetrieveOptions(): 256 """Retrieve command line arguments and options.""" 257 258 # Get options... 259 global Options 260 Options = docopt(_docoptUsage_) 261 262 # Set current working directory to the specified directory... 263 WorkingDir = Options["--workingdir"] 264 if WorkingDir: 265 os.chdir(WorkingDir) 266 267 # Handle examples option... 268 if "--examples" in Options and Options["--examples"]: 269 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 270 sys.exit(0) 271 272 273 def ValidateOptions(): 274 """Validate option values.""" 275 276 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"]) 277 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol") 278 279 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd") 280 MiscUtil.ValidateOptionsOutputFileOverwrite( 281 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 282 ) 283 MiscUtil.ValidateOptionsDistinctFileNames( 284 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"] 285 ) 286 287 MiscUtil.ValidateOptionTextValue( 288 "-c, --charge", Options["--charge"], "AM1BCC AM1-Mulliken Espaloma Gasteiger MMFF94 NAGL" 289 ) 290 291 MiscUtil.ValidateOptionTextValue("--loggingLevel", Options["--loggingLevel"], "Info Warning") 292 293 294 # Setup a usage string for docopt... 295 _docoptUsage_ = """ 296 OpenFECalculatePartialCharges.py - Generate ligand network 297 298 Usage: 299 OpenFECalculatePartialCharges.py [--charge <text>] [--chargeParams <Name,Value,..>] 300 [--infileParams <Name,Value,...>] [--loggingLevel <Info or Warning>] 301 [--outfileParams <Name,Value,...>] [--overwrite] [-w <dir>] -i <infile> -o <outfile> 302 OpenFECalculatePartialCharges.py -h | --help | -e | --examples 303 304 Description: 305 Calculate partial atomic charges for molecules in an input file and write 306 them out to a SD file. 307 308 The partial charges are written to SD file as values of the data field label 309 'atom.dprop.PartialCharge'. These values are automatically processed by 310 RDKit during the loading of a SD file and assigned to atoms. You may retrieve 311 these value using RDKit method Atom.GetDoubleProp('PartialCharge') 312 313 You must specify a valid input file containing 3D coordinates for all 314 molecules. In addition, the hydrogens must be present for all molecules 315 in the input file. 316 317 The supported input file formats are: Mol (.mol), SD (.sdf, .sd) 318 319 The supported output file formats is: SD (.sdf, .sd) 320 321 Options: 322 -e, --examples 323 Print examples. 324 -h, --help 325 Print this help message. 326 -i, --infile <infile> 327 Input file name. 328 -c, --charge <text> [default: AM1BCC] 329 Type of partial atomic charges to calculate. Possible values: AM1BCC, 330 AM1-Mulliken, Espaloma, Gasteiger, MMFF94, or NAGL. 331 --chargeParams <Name,Value,..> [default: auto] 332 A comma delimited list of parameter name and value pairs for calculating 333 partial aromic charges. 334 335 The supported parameter names along with their default values are 336 337 naglModel, auto [ Possible value: A valid NAGL model name. By 338 default, it corresponds to the latest AM1BCC production model ] 339 toolkit, auto [ Possible values: RDKit or AmberTools. Default value: 340 RDKit for Gasteiger and MMFF94; AmberTools for AM1BCC and 341 AM1-Mulliken; Not used for Espaloma and NAGL. ] 342 343 numProcessors, 1 [ Only used for AM1BCC, AM1-Mulliken, Espaloma, 344 and NAGL ] 345 346 precision, 4 347 lineSize, 90 348 349 useConformer, auto [ Use current conformer. Possible values: yes or 350 no. Default value: no for Gasteiger using AmberToolkit; 351 otherwise, yes. ] 352 353 A brief description of parameters is provided below: 354 355 naglModel: NAGL model name. The latest AM1BCC NAGL production 356 model is used by default. You must specify it explicitly in case no 357 production model is available. 358 toolkit: Toolkit name. RDKit for Gasteiger and MMFF94; AmberTools 359 for AM1BCC, AM1-Mulliken, and Gasteiger. 360 361 numProcessors: Number of processors. This is only used during the 362 calculation of AM1BCC, AM1-Mulliken, Espaloma, and NAGL 363 employing OpenFE method bulk_assign_partial_charges(). 364 365 precision: Floating point precision for writing the calculated 366 partial atomic charges. 367 lineSize: Line size for writing the calculated partial aromic 368 charges to SD file as a string value for data field label 369 'atom.dprop.PartialCharge'. 370 371 useConformer: Use current conformer. The current conformer is 372 always used to calculate AM1BCC, Espaloma abd NAGL charges 373 using OpenFE method bulk_assign_partial_charges() and this 374 option is ignored. In addition, the option value is passed to 375 OpenFF method assign_partial_charges() during the calculation 376 of AM1-Mulliken, Gasteiger and MMFF94 charges employing 377 AmberTools or RDKit. The RDKit functions, however, ignore the 378 conformer during the calculation of Gasteiger and MMFF94 379 charges. The current conformer appears not used to calculate 380 Gasteiger charges employing AmberTools. 381 382 --infileParams <Name,Value,...> [default: auto] 383 A comma delimited list of parameter name and value pairs for reading 384 molecules from files. The supported parameter names for different file 385 formats, along with their default values, are shown below: 386 387 SD,MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes 388 389 --loggingLevel <Info or Warning> [default: Warning] 390 Logging level to configure the 'root logger' via logging.basicConfig() 391 function. The default logging level is changed from 'logging.INFO' to 392 'logging.WARNING'. Otherwise, OpenFE and its associated modules 393 may generate a lot of informational messages. 394 -o, --outfile <outfile> 395 Output file name. 396 --outfileParams <Name,Value,...> [default: auto] 397 A comma delimited list of parameter name and value pairs for writing 398 molecules to files. The supported parameter names for different file 399 formats, along with their default values, are shown below: 400 401 SD: kekulize,yes,forceV3000,no 402 403 --overwrite 404 Overwrite existing files. 405 -w, --workingdir <dir> 406 Location of working directory which defaults to the current directory. 407 408 Examples: 409 To calculate AM1BCC partial atomic charges for molecules in a SD file and 410 file and write them out to a file, type: 411 412 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 413 -o SampleTyk2LigandsAM1BCCOut.sdf 414 415 To run the first example for calculating AM1-Mulliken charges using 416 AmberTools, type: 417 418 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 419 -o SampleTyk2LigandsAM1MullikenOut.sdf -c AM1-Mulliken 420 421 To run the first example for calculating Gasteiger charges using RDKit, 422 type: 423 424 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 425 -o SampleTyk2LigandsGasteigerRDKit.sdf -c Gasteiger 426 427 To run the first example for calculating Gasteiger charges using 428 AmberTools, type: 429 430 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 431 -o SampleTyk2LigandsasteigerAmberTools.sdf -c Gasteiger 432 --chargeParams "toolkit,AmberTools" 433 434 To run the first example for calculating NAGL charges using the default 435 production AM1BCC model, type: 436 437 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 438 -o SampleTyk2LigandsNAGL.sdf -c NAGL 439 440 To run the first example for calculating NAGL charges using a specific 441 model, type: 442 443 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 444 -o SampleTyk2LigandsNAGL.sdf -c NAGL 445 --chargeParams "naglmodel, openff-gnn-am1bcc-0.1.0-rc.3.pt" 446 447 To run the first example by specifying explicit values for various parameters, 448 type: 449 450 % OpenFECalculatePartialCharges.py -i SampleTyk2Ligands.sdf 451 -o SampleTyk2LigandsAM1BCCOut.sdf -c AM1BCC 452 --chargeParams "numProcessors, 4, precision, 4, lineSize, 90" 453 --loggingLevel Warning 454 455 Author: 456 Manish Sud(msud@san.rr.com) 457 458 See also: 459 OpenFECalculateAbsoluteHydrationFreeEnergy.py, OpenFEGenerateLigandNetwork.py, 460 OpenFECalculateRelativeBindingFreeEnergy.py, 461 OpenFECalculateRelativeHydrationFreeEnergy.py 462 463 Copyright: 464 Copyright (C) 2026 Manish Sud. All rights reserved. 465 466 The functionality available in this script is implemented using OpenFE, an 467 open source molecuar for alchemical free energy calculations. 468 469 This file is part of MayaChemTools. 470 471 MayaChemTools is free software; you can redistribute it and/or modify it under 472 the terms of the GNU Lesser General Public License as published by the Free 473 Software Foundation; either version 3 of the License, or (at your option) any 474 later version. 475 476 """ 477 478 if __name__ == "__main__": 479 main()
