1 #!/bin/env python 2 # 3 # File: RDKitPerformMinimization.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 RDKit, an 9 # open source toolkit for cheminformatics developed by Greg Landrum. 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 multiprocessing as mp 36 37 # RDKit imports... 38 try: 39 from rdkit import rdBase 40 from rdkit import Chem 41 from rdkit.Chem import AllChem 42 except ImportError as ErrMsg: 43 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg) 44 sys.stderr.write("Check/update your RDKit environment and try again.\n\n") 45 sys.exit(1) 46 47 # MayaChemTools imports... 48 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python")) 49 try: 50 from docopt import docopt 51 import MiscUtil 52 import RDKitUtil 53 except ImportError as ErrMsg: 54 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg) 55 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n") 56 sys.exit(1) 57 58 ScriptName = os.path.basename(sys.argv[0]) 59 Options = {} 60 OptionsInfo = {} 61 62 63 def main(): 64 """Start execution of the script.""" 65 66 MiscUtil.PrintInfo( 67 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n" 68 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime()) 69 ) 70 71 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime() 72 73 # Retrieve command line arguments and options... 74 RetrieveOptions() 75 76 # Process and validate command line arguments and options... 77 ProcessOptions() 78 79 # Perform actions required by the script... 80 PerformMinimization() 81 82 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName) 83 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime)) 84 85 86 def PerformMinimization(): 87 """Perform minimization.""" 88 89 # Setup a molecule reader... 90 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"]) 91 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"]) 92 93 # Set up a molecule writer... 94 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"]) 95 if Writer is None: 96 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"]) 97 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"]) 98 99 MolCount, ValidMolCount, MinimizationFailedCount, WriteFailedCount = ProcessMolecules(Mols, Writer) 100 101 if Writer is not None: 102 Writer.close() 103 104 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount) 105 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount) 106 MiscUtil.PrintInfo( 107 "Number of molecules failed during conformation generation or minimization: %d" % MinimizationFailedCount 108 ) 109 MiscUtil.PrintInfo("Number of molecules failed during writing: %d" % WriteFailedCount) 110 MiscUtil.PrintInfo( 111 "Number of ignored molecules: %d" % (MolCount - ValidMolCount + MinimizationFailedCount + WriteFailedCount) 112 ) 113 114 115 def ProcessMolecules(Mols, Writer): 116 """Process and minimize molecules.""" 117 118 if OptionsInfo["MPMode"]: 119 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer) 120 else: 121 return ProcessMoleculesUsingSingleProcess(Mols, Writer) 122 123 124 def ProcessMoleculesUsingSingleProcess(Mols, Writer): 125 """Process and minimize molecules using a single process.""" 126 127 if OptionsInfo["SkipConformerGeneration"]: 128 MiscUtil.PrintInfo("\nPerforming minimization without generation of conformers...") 129 else: 130 MiscUtil.PrintInfo("\nPerforming minimization with generation of conformers...") 131 132 (MolCount, ValidMolCount, MinimizationFailedCount, WriteFailedCount) = [0] * 4 133 for Mol in Mols: 134 MolCount += 1 135 136 if Mol is None: 137 continue 138 139 if RDKitUtil.IsMolEmpty(Mol): 140 if not OptionsInfo["QuietMode"]: 141 MolName = RDKitUtil.GetMolName(Mol, MolCount) 142 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 143 continue 144 ValidMolCount += 1 145 146 Mol, CalcStatus, ConfID, Energy = MinimizeMoleculeOrConformers(Mol, MolCount) 147 148 if not CalcStatus: 149 MinimizationFailedCount += 1 150 continue 151 152 WriteStatus = WriteMolecule(Writer, Mol, MolCount, ConfID, Energy) 153 if not WriteStatus: 154 WriteFailedCount += 1 155 156 return (MolCount, ValidMolCount, MinimizationFailedCount, WriteFailedCount) 157 158 159 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer): 160 """Process and minimize molecules using multiprocessing.""" 161 162 if OptionsInfo["SkipConformerGeneration"]: 163 MiscUtil.PrintInfo("\nPerforming minimization without generation of conformers using multiprocessing...") 164 else: 165 MiscUtil.PrintInfo("\nPerforming minimization with generation of conformers using multiprocessing...") 166 167 MPParams = OptionsInfo["MPParams"] 168 169 # Setup data for initializing a worker process... 170 InitializeWorkerProcessArgs = ( 171 MiscUtil.ObjectToBase64EncodedString(Options), 172 MiscUtil.ObjectToBase64EncodedString(OptionsInfo), 173 ) 174 175 # Setup a encoded mols data iterable for a worker process... 176 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols) 177 178 # Setup process pool along with data initialization for each process... 179 MiscUtil.PrintInfo( 180 "\nConfiguring multiprocessing using %s method..." 181 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()") 182 ) 183 MiscUtil.PrintInfo( 184 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n" 185 % ( 186 MPParams["NumProcesses"], 187 MPParams["InputDataMode"], 188 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]), 189 ) 190 ) 191 192 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs) 193 194 # Start processing... 195 if re.match("^Lazy$", MPParams["InputDataMode"], re.I): 196 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 197 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I): 198 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"]) 199 else: 200 MiscUtil.PrintError( 201 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"]) 202 ) 203 204 (MolCount, ValidMolCount, MinimizationFailedCount, WriteFailedCount) = [0] * 4 205 for Result in Results: 206 MolCount += 1 207 MolIndex, EncodedMol, CalcStatus, ConfID, Energy = Result 208 209 if EncodedMol is None: 210 continue 211 ValidMolCount += 1 212 213 if not CalcStatus: 214 MinimizationFailedCount += 1 215 continue 216 217 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 218 WriteStatus = WriteMolecule(Writer, Mol, MolCount, ConfID, Energy) 219 if not WriteStatus: 220 WriteFailedCount += 1 221 222 return (MolCount, ValidMolCount, MinimizationFailedCount, WriteFailedCount) 223 224 225 def InitializeWorkerProcess(*EncodedArgs): 226 """Initialize data for a worker process.""" 227 228 global Options, OptionsInfo 229 230 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid()) 231 232 # Decode Options and OptionInfo... 233 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0]) 234 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1]) 235 236 237 def WorkerProcess(EncodedMolInfo): 238 """Process data for a worker process.""" 239 240 MolIndex, EncodedMol = EncodedMolInfo 241 242 CalcStatus = False 243 ConfID = None 244 Energy = None 245 246 if EncodedMol is None: 247 return [MolIndex, None, CalcStatus, ConfID, Energy] 248 249 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol) 250 if RDKitUtil.IsMolEmpty(Mol): 251 if not OptionsInfo["QuietMode"]: 252 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1)) 253 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName) 254 return [MolIndex, None, CalcStatus, ConfID, Energy] 255 256 Mol, CalcStatus, ConfID, Energy = MinimizeMoleculeOrConformers(Mol, (MolIndex + 1)) 257 258 return [ 259 MolIndex, 260 RDKitUtil.MolToBase64EncodedMolString( 261 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps 262 ), 263 CalcStatus, 264 ConfID, 265 Energy, 266 ] 267 268 269 def MinimizeMoleculeOrConformers(Mol, MolNum=None): 270 """Minimize molecule or conformers.""" 271 272 ConfID = None 273 if OptionsInfo["SkipConformerGeneration"]: 274 Mol, CalcStatus, Energy = MinimizeMolecule(Mol, MolNum) 275 else: 276 Mol, CalcStatus, ConfID, Energy = MinimizeConformers(Mol, MolNum) 277 278 return (Mol, CalcStatus, ConfID, Energy) 279 280 281 def MinimizeMolecule(Mol, MolNum=None): 282 """Minimize molecule.""" 283 284 if OptionsInfo["AddHydrogens"]: 285 Mol = Chem.AddHs(Mol, addCoords=True) 286 287 Status = 0 288 try: 289 if OptionsInfo["UseUFF"]: 290 Status = AllChem.UFFOptimizeMolecule(Mol, maxIters=OptionsInfo["MaxIters"]) 291 elif OptionsInfo["UseMMFF"]: 292 Status = AllChem.MMFFOptimizeMolecule( 293 Mol, maxIters=OptionsInfo["MaxIters"], mmffVariant=OptionsInfo["MMFFVariant"] 294 ) 295 else: 296 MiscUtil.PrintError( 297 "Minimization couldn't be performed: Specified forcefield, %s, is not supported" 298 % OptionsInfo["ForceField"] 299 ) 300 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg: 301 if not OptionsInfo["QuietMode"]: 302 MolName = RDKitUtil.GetMolName(Mol, MolNum) 303 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg)) 304 return (Mol, False, None) 305 306 if Status != 0: 307 if not OptionsInfo["QuietMode"]: 308 MolName = RDKitUtil.GetMolName(Mol, MolNum) 309 MiscUtil.PrintWarning( 310 'Minimization failed to converge for molecule %s in %d steps. Try using higher value for "--maxIters" option...\n' 311 % (MolName, OptionsInfo["MaxIters"]) 312 ) 313 314 Energy = None 315 if OptionsInfo["EnergyOut"]: 316 EnergyStatus, Energy = GetEnergy(Mol) 317 if EnergyStatus: 318 Energy = "%.2f" % Energy 319 else: 320 if not OptionsInfo["QuietMode"]: 321 MolName = RDKitUtil.GetMolName(Mol, MolNum) 322 MiscUtil.PrintWarning( 323 "Failed to retrieve calculated energy for molecule %s. Try again after removing any salts or cleaing up the molecule...\n" 324 % (MolName) 325 ) 326 327 if OptionsInfo["RemoveHydrogens"]: 328 Mol = Chem.RemoveHs(Mol) 329 330 return (Mol, True, Energy) 331 332 333 def MinimizeConformers(Mol, MolNum=None): 334 """Generate and minimize conformers for a molecule to get the lowest energy conformer.""" 335 336 if OptionsInfo["AddHydrogens"]: 337 Mol = Chem.AddHs(Mol) 338 339 ConfIDs = EmbedMolecule(Mol, MolNum) 340 if not len(ConfIDs): 341 if not OptionsInfo["QuietMode"]: 342 MolName = RDKitUtil.GetMolName(Mol, MolNum) 343 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName) 344 return (Mol, False, None, None) 345 346 CalcEnergyMap = {} 347 for ConfID in ConfIDs: 348 try: 349 if OptionsInfo["UseUFF"]: 350 Status = AllChem.UFFOptimizeMolecule(Mol, confId=ConfID, maxIters=OptionsInfo["MaxIters"]) 351 elif OptionsInfo["UseMMFF"]: 352 Status = AllChem.MMFFOptimizeMolecule( 353 Mol, confId=ConfID, maxIters=OptionsInfo["MaxIters"], mmffVariant=OptionsInfo["MMFFVariant"] 354 ) 355 else: 356 MiscUtil.PrintError( 357 "Minimization couldn't be performed: Specified forcefield, %s, is not supported" 358 % OptionsInfo["ForceField"] 359 ) 360 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg: 361 if not OptionsInfo["QuietMode"]: 362 MolName = RDKitUtil.GetMolName(Mol, MolNum) 363 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg)) 364 return (Mol, False, None, None) 365 366 EnergyStatus, Energy = GetEnergy(Mol, ConfID) 367 if not EnergyStatus: 368 if not OptionsInfo["QuietMode"]: 369 MolName = RDKitUtil.GetMolName(Mol, MolNum) 370 MiscUtil.PrintWarning( 371 "Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n" 372 % (ConfID, MolName) 373 ) 374 return (Mol, False, None, None) 375 376 if Status != 0: 377 if not OptionsInfo["QuietMode"]: 378 MolName = RDKitUtil.GetMolName(Mol, MolNum) 379 MiscUtil.PrintWarning( 380 'Minimization failed to converge for conformation number %d of molecule %s in %d steps. Try using higher value for "--maxIters" option...\n' 381 % (ConfID, MolName, OptionsInfo["MaxIters"]) 382 ) 383 384 CalcEnergyMap[ConfID] = Energy 385 386 SortedConfIDs = sorted(ConfIDs, key=lambda ConfID: CalcEnergyMap[ConfID]) 387 MinEnergyConfID = SortedConfIDs[0] 388 389 if OptionsInfo["RemoveHydrogens"]: 390 Mol = Chem.RemoveHs(Mol) 391 392 Energy = "%.2f" % CalcEnergyMap[MinEnergyConfID] if OptionsInfo["EnergyOut"] else None 393 394 return (Mol, True, MinEnergyConfID, Energy) 395 396 397 def GetEnergy(Mol, ConfID=None): 398 """Calculate energy.""" 399 400 Status = True 401 Energy = None 402 403 if ConfID is None: 404 ConfID = -1 405 406 if OptionsInfo["UseUFF"]: 407 UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId=ConfID) 408 if UFFMoleculeForcefield is None: 409 Status = False 410 else: 411 Energy = UFFMoleculeForcefield.CalcEnergy() 412 elif OptionsInfo["UseMMFF"]: 413 MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant=OptionsInfo["MMFFVariant"]) 414 MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId=ConfID) 415 if MMFFMoleculeForcefield is None: 416 Status = False 417 else: 418 Energy = MMFFMoleculeForcefield.CalcEnergy() 419 else: 420 MiscUtil.PrintError( 421 "Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"] 422 ) 423 424 return (Status, Energy) 425 426 427 def EmbedMolecule(Mol, MolNum=None): 428 """Embed conformations.""" 429 430 ConfIDs = [] 431 432 MaxConfs = OptionsInfo["MaxConfs"] 433 RandomSeed = OptionsInfo["RandomSeed"] 434 EnforceChirality = OptionsInfo["EnforceChirality"] 435 UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"] 436 ETVersion = OptionsInfo["ETVersion"] 437 UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"] 438 439 try: 440 ConfIDs = AllChem.EmbedMultipleConfs( 441 Mol, 442 numConfs=MaxConfs, 443 randomSeed=RandomSeed, 444 enforceChirality=EnforceChirality, 445 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs, 446 useBasicKnowledge=UseBasicKnowledge, 447 ETversion=ETVersion, 448 ) 449 except ValueError as ErrMsg: 450 if not OptionsInfo["QuietMode"]: 451 MolName = RDKitUtil.GetMolName(Mol, MolNum) 452 MiscUtil.PrintWarning("Embedding failed for molecule %s:\n%s\n" % (MolName, ErrMsg)) 453 ConfIDs = [] 454 455 return ConfIDs 456 457 458 def WriteMolecule(Writer, Mol, MolNum=None, ConfID=None, Energy=None): 459 """Write molecule.""" 460 461 if Energy is not None: 462 Mol.SetProp(OptionsInfo["EnergyLabel"], Energy) 463 464 try: 465 if ConfID is None: 466 Writer.write(Mol) 467 else: 468 Writer.write(Mol, confId=ConfID) 469 except (ValueError, RuntimeError) as ErrMsg: 470 if not OptionsInfo["QuietMode"]: 471 MolName = RDKitUtil.GetMolName(Mol, MolNum) 472 MiscUtil.PrintWarning("Failed to write molecule %s:\n%s\n" % (MolName, ErrMsg)) 473 return False 474 475 return True 476 477 478 def ProcesssConformerGeneratorOption(): 479 """Process comformer generator option.""" 480 481 ConfGenParams = MiscUtil.ProcessOptionConformerGenerator("--conformerGenerator", Options["--conformerGenerator"]) 482 483 OptionsInfo["ConformerGenerator"] = ConfGenParams["ConformerGenerator"] 484 OptionsInfo["SkipConformerGeneration"] = ConfGenParams["SkipConformerGeneration"] 485 OptionsInfo["UseBasicKnowledge"] = ConfGenParams["UseBasicKnowledge"] 486 OptionsInfo["UseExpTorsionAnglePrefs"] = ConfGenParams["UseExpTorsionAnglePrefs"] 487 OptionsInfo["ETVersion"] = ConfGenParams["ETVersion"] 488 489 490 def ProcessOptions(): 491 """Process and validate command line arguments and options.""" 492 493 MiscUtil.PrintInfo("Processing options...") 494 495 # Validate options... 496 ValidateOptions() 497 498 OptionsInfo["Infile"] = Options["--infile"] 499 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters( 500 "--infileParams", Options["--infileParams"], Options["--infile"] 501 ) 502 503 OptionsInfo["Outfile"] = Options["--outfile"] 504 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters( 505 "--outfileParams", Options["--outfileParams"] 506 ) 507 508 OptionsInfo["Overwrite"] = Options["--overwrite"] 509 510 OptionsInfo["AddHydrogens"] = True if re.match("^yes$", Options["--addHydrogens"], re.I) else False 511 512 ProcesssConformerGeneratorOption() 513 514 if re.match("^UFF$", Options["--forceField"], re.I): 515 ForceField = "UFF" 516 UseUFF = True 517 UseMMFF = False 518 elif re.match("^MMFF$", Options["--forceField"], re.I): 519 ForceField = "MMFF" 520 UseUFF = False 521 UseMMFF = True 522 else: 523 MiscUtil.PrintError( 524 'The value, %s, specified for "--forceField" is not supported.' % (Options["--forceField"],) 525 ) 526 527 MMFFVariant = "MMFF94" if re.match("^MMFF94$", Options["--forceFieldMMFFVariant"], re.I) else "MMFF94s" 528 529 OptionsInfo["ForceField"] = ForceField 530 OptionsInfo["MMFFVariant"] = MMFFVariant 531 OptionsInfo["UseMMFF"] = UseMMFF 532 OptionsInfo["UseUFF"] = UseUFF 533 534 OptionsInfo["EnergyOut"] = True if re.match("^yes$", Options["--energyOut"], re.I) else False 535 if UseMMFF: 536 OptionsInfo["EnergyLabel"] = "%s_Energy" % MMFFVariant 537 else: 538 OptionsInfo["EnergyLabel"] = "%s_Energy" % ForceField 539 540 OptionsInfo["EnforceChirality"] = True if re.match("^yes$", Options["--enforceChirality"], re.I) else False 541 542 OptionsInfo["MaxIters"] = int(Options["--maxIters"]) 543 OptionsInfo["MaxConfs"] = int(Options["--maxConfs"]) 544 545 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False 546 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"]) 547 548 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False 549 550 RandomSeed = -1 551 if not re.match("^auto$", Options["--randomSeed"], re.I): 552 RandomSeed = int(Options["--randomSeed"]) 553 OptionsInfo["RandomSeed"] = RandomSeed 554 555 OptionsInfo["RemoveHydrogens"] = True if re.match("^yes$", Options["--removeHydrogens"], re.I) else False 556 557 558 def RetrieveOptions(): 559 """Retrieve command line arguments and options.""" 560 561 # Get options... 562 global Options 563 Options = docopt(_docoptUsage_) 564 565 # Set current working directory to the specified directory... 566 WorkingDir = Options["--workingdir"] 567 if WorkingDir: 568 os.chdir(WorkingDir) 569 570 # Handle examples option... 571 if "--examples" in Options and Options["--examples"]: 572 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_)) 573 sys.exit(0) 574 575 576 def ValidateOptions(): 577 """Validate option values.""" 578 579 MiscUtil.ValidateOptionTextValue("-a, --addHydrogens", Options["--addHydrogens"], "yes no") 580 MiscUtil.ValidateOptionTextValue( 581 "-c, --conformerGenerator", Options["--conformerGenerator"], "SDG ETDG KDG ETKDG ETKDGv2 None" 582 ) 583 584 MiscUtil.ValidateOptionTextValue("-f, --forceField", Options["--forceField"], "UFF MMFF") 585 MiscUtil.ValidateOptionTextValue(" --forceFieldMMFFVariant", Options["--forceFieldMMFFVariant"], "MMFF94 MMFF94s") 586 587 MiscUtil.ValidateOptionTextValue("--energyOut", Options["--energyOut"], "yes no") 588 MiscUtil.ValidateOptionTextValue("--enforceChirality ", Options["--enforceChirality"], "yes no") 589 590 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"]) 591 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv") 592 593 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd") 594 MiscUtil.ValidateOptionsOutputFileOverwrite( 595 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"] 596 ) 597 MiscUtil.ValidateOptionsDistinctFileNames( 598 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"] 599 ) 600 601 MiscUtil.ValidateOptionIntegerValue("--maxConfs", Options["--maxConfs"], {">": 0}) 602 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0}) 603 604 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no") 605 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no") 606 607 if not re.match("^auto$", Options["--randomSeed"], re.I): 608 MiscUtil.ValidateOptionIntegerValue("--randomSeed", Options["--randomSeed"], {}) 609 610 MiscUtil.ValidateOptionTextValue("-r, --removeHydrogens", Options["--removeHydrogens"], "yes no") 611 612 613 # Setup a usage string for docopt... 614 _docoptUsage_ = """ 615 RDKitPerformMinimization.py - Perform structure minimization 616 617 Usage: 618 RDKitPerformMinimization.py [--addHydrogens <yes or no>] [--conformerGenerator <SDG, KDG, ETDG, ETKDG, ETKDGv2, None>] 619 [--forceField <UFF, or MMFF>] [--forceFieldMMFFVariant <MMFF94 or MMFF94s>] 620 [--energyOut <yes or no>] [--enforceChirality <yes or no>] [--infileParams <Name,Value,...>] 621 [--maxConfs <number>] [--maxIters <number>] [--mp <yes or no>] [--mpParams <Name,Value,...>] 622 [ --outfileParams <Name,Value,...> ] [--overwrite] [--quiet <yes or no>] [ --removeHydrogens <yes or no>] 623 [--randomSeed <number>] [-w <dir>] -i <infile> -o <outfile> 624 RDKitPerformMinimization.py -h | --help | -e | --examples 625 626 Description: 627 Generate 3D structures for molecules using a combination of distance geometry 628 and forcefield minimization or minimize existing 3D structures using a specified 629 forcefield. 630 631 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi 632 .csv, .tsv, .txt) 633 634 The supported output file formats are: SD (.sdf, .sd) 635 636 Options: 637 -a, --addHydrogens <yes or no> [default: yes] 638 Add hydrogens before minimization. 639 -c, --conformerGenerator <text or None> [default: ETKDGv2] 640 Conformation generation methodology for generating initial 3D coordinates. The 641 possible values along with a brief description are shown below: 642 643 SDG: Standard Distance Geometry 644 KDG: basic Knowledge-terms with Distance Geometry 645 ETDG: Experimental Torsion-angle preference with Distance Geometry 646 ETKDG: Experimental Torsion-angle preference along with basic 647 Knowledge-terms and Distance Geometry [Ref 129] 648 ETKDGv2: Experimental Torsion-angle preference along with basic 649 Knowledge-terms and Distance Geometry [Ref 167] 650 None: No conformation generation 651 652 The conformation generation step may be skipped by specifying 'None' value to 653 perform only forcefield minimization of molecules with 3D structures in input 654 file. This doesn't work for molecules in SMILES file or molecules in SD/MOL files 655 containing 2D structures. 656 -f, --forceField <UFF, MMFF> [default: MMFF] 657 Forcefield method to use for energy minimization. Possible values: Universal Force 658 Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force Field [ Ref 83-87 ] . 659 --forceFieldMMFFVariant <MMFF94 or MMFF94s> [default: MMFF94] 660 Variant of MMFF forcefield to use for energy minimization. 661 --energyOut <yes or no> [default: No] 662 Write out energy values. 663 --enforceChirality <yes or no> [default: Yes] 664 Enforce chirality for defined chiral centers. 665 -e, --examples 666 Print examples. 667 -h, --help 668 Print this help message. 669 -i, --infile <infile> 670 Input file name. 671 --infileParams <Name,Value,...> [default: auto] 672 A comma delimited list of parameter name and value pairs for reading 673 molecules from files. The supported parameter names for different file 674 formats, along with their default values, are shown below: 675 676 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes 677 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space, 678 smilesTitleLine,auto,sanitize,yes 679 680 Possible values for smilesDelimiter: space, comma or tab. 681 --maxConfs <number> [default: 250] 682 Maximum number of conformations to generate for each molecule by conformation 683 generation methodology for initial 3D coordinates. The conformations are minimized 684 using the specified forcefield and the lowest energy conformation is written to the 685 output file. This option is ignored during 'None' value of '-c --conformerGenerator' 686 option. 687 --maxIters <number> [default: 500] 688 Maximum number of iterations to perform for each molecule during forcefield 689 minimization. 690 --mp <yes or no> [default: no] 691 Use multiprocessing. 692 693 By default, input data is retrieved in a lazy manner via mp.Pool.imap() 694 function employing lazy RDKit data iterable. This allows processing of 695 arbitrary large data sets without any additional requirements memory. 696 697 All input data may be optionally loaded into memory by mp.Pool.map() 698 before starting worker processes in a process pool by setting the value 699 of 'inputDataMode' to 'InMemory' in '--mpParams' option. 700 701 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input 702 data mode may adversely impact the performance. The '--mpParams' section 703 provides additional information to tune the value of 'chunkSize'. 704 --mpParams <Name,Value,...> [default: auto] 705 A comma delimited list of parameter name and value pairs to configure 706 multiprocessing. 707 708 The supported parameter names along with their default and possible 709 values are shown below: 710 711 chunkSize, auto 712 inputDataMode, Lazy [ Possible values: InMemory or Lazy ] 713 numProcesses, auto [ Default: mp.cpu_count() ] 714 715 These parameters are used by the following functions to configure and 716 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and 717 mp.Pool.imap(). 718 719 The chunkSize determines chunks of input data passed to each worker 720 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions. 721 The default value of chunkSize is dependent on the value of 'inputDataMode'. 722 723 The mp.Pool.map() function, invoked during 'InMemory' input data mode, 724 automatically converts RDKit data iterable into a list, loads all data into 725 memory, and calculates the default chunkSize using the following method 726 as shown in its code: 727 728 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4) 729 if extra: chunkSize += 1 730 731 For example, the default chunkSize will be 7 for a pool of 4 worker processes 732 and 100 data items. 733 734 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs 735 'lazy' RDKit data iterable to retrieve data as needed, without loading all the 736 data into memory. Consequently, the size of input data is not known a priori. 737 It's not possible to estimate an optimal value for the chunkSize. The default 738 chunkSize is set to 1. 739 740 The default value for the chunkSize during 'Lazy' data mode may adversely 741 impact the performance due to the overhead associated with exchanging 742 small chunks of data. It is generally a good idea to explicitly set chunkSize to 743 a larger value during 'Lazy' input data mode, based on the size of your input 744 data and number of processes in the process pool. 745 746 The mp.Pool.map() function waits for all worker processes to process all 747 the data and return the results. The mp.Pool.imap() function, however, 748 returns the the results obtained from worker processes as soon as the 749 results become available for specified chunks of data. 750 751 The order of data in the results returned by both mp.Pool.map() and 752 mp.Pool.imap() functions always corresponds to the input data. 753 -o, --outfile <outfile> 754 Output file name. 755 --outfileParams <Name,Value,...> [default: auto] 756 A comma delimited list of parameter name and value pairs for writing 757 molecules to files. The supported parameter names for different file 758 formats, along with their default values, are shown below: 759 760 SD: kekulize,yes,forceV3000,no 761 762 --overwrite 763 Overwrite existing files. 764 -q, --quiet <yes or no> [default: no] 765 Use quiet mode. The warning and information messages will not be printed. 766 -r, --removeHydrogens <yes or no> [default: Yes] 767 Remove hydrogens after minimization. 768 --randomSeed <number> [default: auto] 769 Seed for the random number generator for reproducing 3D coordinates. 770 Default is to use a random seed. 771 -w, --workingdir <dir> 772 Location of working directory which defaults to the current directory. 773 774 Examples: 775 To generate up to 250 conformations using ETKDG methodology followed by MMFF 776 forcefield minimization for a maximum of 500 iterations for molecules in a SMILES file 777 and write out a SD file containing minimum energy structure corresponding to each 778 molecule, type: 779 780 % RDKitPerformMinimization.py -i Sample.smi -o SampleOut.sdf 781 782 To rerun the first example in a quiet mode and write out a SD file, type: 783 784 % RDKitPerformMinimization.py -q yes -i Sample.smi -o SampleOut.sdf 785 786 To run the first example in multiprocessing mode on all available CPUs 787 without loading all data into memory and write out a SD file, type: 788 789 % RDKitPerformMinimization.py --mp yes -i Sample.smi -o SampleOut.sdf 790 791 To run the first example in multiprocessing mode on all available CPUs 792 by loading all data into memory and write out a SD file, type: 793 794 % RDKitPerformMinimization.py --mp yes --mpParams "inputDataMode, 795 InMemory" -i Sample.smi -o SampleOut.sdf 796 797 To run the first example in multiprocessing mode on specific number of 798 CPUs and chunk size without loading all data into memory and write out a SD file, 799 type: 800 801 % RDKitPerformMinimization.py --mp yes --mpParams "inputDataMode,Lazy, 802 numProcesses,4,chunkSize,8" -i Sample.smi -o SampleOut.sdf 803 804 To generate up to 150 conformations using ETKDG methodology followed by MMFF 805 forcefield minimization for a maximum of 250 iterations along with a specified random 806 seed for molecules in a SMILES file and write out a SD file containing minimum energy 807 structures corresponding to each molecule, type 808 809 % RDKitPerformMinimization.py --maxConfs 150 --randomSeed 201780117 810 --maxIters 250 -i Sample.smi -o SampleOut.sdf 811 812 To minimize structures in a 3D SD file using UFF forcefield for a maximum of 150 813 iterations without generating any conformations and write out a SD file containing 814 minimum energy structures corresponding to each molecule, type 815 816 % RDKitPerformMinimization.py -c None -f UFF --maxIters 150 817 -i Sample3D.sdf -o SampleOut.sdf 818 819 To generate up to 50 conformations using SDG methodology followed 820 by UFF forcefield minimization for a maximum of 50 iterations for 821 molecules in a CSV SMILES file, SMILES strings in column 1, name in 822 column 2, and write out a SD file, type: 823 824 % RDKitPerformMinimization.py --maxConfs 50 --maxIters 50 -c SDG 825 -f UFF --infileParams "smilesDelimiter,comma,smilesTitleLine,yes, 826 smilesColumn,1,smilesNameColumn,2" -i SampleSMILES.csv 827 -o SampleOut.sdf 828 829 Author: 830 Manish Sud(msud@san.rr.com) 831 832 See also: 833 RDKitCalculateRMSD.py, RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py, 834 RDKitConvertFileFormat.py, RDKitGenerateConformers.py, RDKitPerformConstrainedMinimization.py 835 836 Copyright: 837 Copyright (C) 2026 Manish Sud. All rights reserved. 838 839 The functionality available in this script is implemented using RDKit, an 840 open source toolkit for cheminformatics developed by Greg Landrum. 841 842 This file is part of MayaChemTools. 843 844 MayaChemTools is free software; you can redistribute it and/or modify it under 845 the terms of the GNU Lesser General Public License as published by the Free 846 Software Foundation; either version 3 of the License, or (at your option) any 847 later version. 848 849 """ 850 851 if __name__ == "__main__": 852 main()
