1 #!/bin/env python
2 #
3 # File: RDKitGenerateConformers.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 from rdkit.Chem import Descriptors
43 except ImportError as ErrMsg:
44 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
45 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
46 sys.exit(1)
47
48 # MayaChemTools imports...
49 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
50 try:
51 from docopt import docopt
52 import MiscUtil
53 import RDKitUtil
54 except ImportError as ErrMsg:
55 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
56 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
57 sys.exit(1)
58
59 ScriptName = os.path.basename(sys.argv[0])
60 Options = {}
61 OptionsInfo = {}
62
63
64 def main():
65 """Start execution of the script."""
66
67 MiscUtil.PrintInfo(
68 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
69 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
70 )
71
72 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
73
74 # Retrieve command line arguments and options...
75 RetrieveOptions()
76
77 # Process and validate command line arguments and options...
78 ProcessOptions()
79
80 # Perform actions required by the script...
81 GenerateConformers()
82
83 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
84 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
85
86
87 def GenerateConformers():
88 """Generate conformers."""
89
90 # Setup a molecule reader...
91 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
92 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
93
94 # Set up a molecule writer...
95 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
96 if Writer is None:
97 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
98 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
99
100 MolCount, ValidMolCount, ConfGenFailedCount = ProcessMolecules(Mols, Writer)
101
102 if Writer is not None:
103 Writer.close()
104
105 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
106 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
107 MiscUtil.PrintInfo(
108 "Number of molecules failed during conformation generation or minimization: %d" % ConfGenFailedCount
109 )
110 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + ConfGenFailedCount))
111
112
113 def ProcessMolecules(Mols, Writer):
114 """Process molecules to generate conformers."""
115
116 if OptionsInfo["MPMode"]:
117 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer)
118 else:
119 return ProcessMoleculesUsingSingleProcess(Mols, Writer)
120
121
122 def ProcessMoleculesUsingSingleProcess(Mols, Writer):
123 """Process molecules to generate conformers using a single process."""
124
125 if OptionsInfo["SkipForceFieldMinimization"]:
126 MiscUtil.PrintInfo("\nGenerating conformers without performing energy minimization...")
127 else:
128 MiscUtil.PrintInfo("\nGenerating conformers and performing energy minimization...")
129
130 (MolCount, ValidMolCount, ConfGenFailedCount) = [0] * 3
131 for Mol in Mols:
132 MolCount += 1
133
134 if Mol is None:
135 continue
136
137 if RDKitUtil.IsMolEmpty(Mol):
138 if not OptionsInfo["QuietMode"]:
139 MolName = RDKitUtil.GetMolName(Mol, MolCount)
140 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
141 continue
142 ValidMolCount += 1
143
144 ConformerMol, CalcStatus, ConfIDs, ConfEnergies = GenerateMolConformers(Mol, MolCount)
145
146 if not CalcStatus:
147 ConfGenFailedCount += 1
148 continue
149
150 WriteMolConformers(Writer, ConformerMol, MolCount, ConfIDs, ConfEnergies)
151
152 return (MolCount, ValidMolCount, ConfGenFailedCount)
153
154
155 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer):
156 """Process and minimize molecules using multiprocessing."""
157
158 if OptionsInfo["SkipForceFieldMinimization"]:
159 MiscUtil.PrintInfo("\nGenerating conformers without performing energy minimization using multiprocessing...")
160 else:
161 MiscUtil.PrintInfo("\nGenerating conformers and performing energy minimization using multiprocessing...")
162
163 MPParams = OptionsInfo["MPParams"]
164
165 # Setup data for initializing a worker process...
166 InitializeWorkerProcessArgs = (
167 MiscUtil.ObjectToBase64EncodedString(Options),
168 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
169 )
170
171 # Setup a encoded mols data iterable for a worker process...
172 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
173
174 # Setup process pool along with data initialization for each process...
175 MiscUtil.PrintInfo(
176 "\nConfiguring multiprocessing using %s method..."
177 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
178 )
179 MiscUtil.PrintInfo(
180 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
181 % (
182 MPParams["NumProcesses"],
183 MPParams["InputDataMode"],
184 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
185 )
186 )
187
188 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
189
190 # Start processing...
191 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
192 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
193 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
194 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
195 else:
196 MiscUtil.PrintError(
197 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
198 )
199
200 (MolCount, ValidMolCount, ConfGenFailedCount) = [0] * 3
201 for Result in Results:
202 MolCount += 1
203 MolIndex, EncodedMol, CalcStatus, ConfIDs, ConfEnergies = Result
204
205 if EncodedMol is None:
206 continue
207 ValidMolCount += 1
208
209 if not CalcStatus:
210 ConfGenFailedCount += 1
211 continue
212
213 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
214 WriteMolConformers(Writer, Mol, MolCount, ConfIDs, ConfEnergies)
215
216 return (MolCount, ValidMolCount, ConfGenFailedCount)
217
218
219 def InitializeWorkerProcess(*EncodedArgs):
220 """Initialize data for a worker process."""
221
222 global Options, OptionsInfo
223
224 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
225
226 # Decode Options and OptionInfo...
227 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
228 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
229
230
231 def WorkerProcess(EncodedMolInfo):
232 """Process data for a worker process."""
233
234 MolIndex, EncodedMol = EncodedMolInfo
235
236 CalcStatus = False
237 ConfIDs = None
238 ConfEnergies = None
239
240 if EncodedMol is None:
241 return [MolIndex, None, CalcStatus, ConfIDs, ConfEnergies]
242
243 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
244 if RDKitUtil.IsMolEmpty(Mol):
245 if not OptionsInfo["QuietMode"]:
246 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
247 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
248 return [MolIndex, None, CalcStatus, ConfIDs, ConfEnergies]
249
250 Mol, CalcStatus, ConfIDs, ConfEnergies = GenerateMolConformers(Mol, (MolIndex + 1))
251
252 return [
253 MolIndex,
254 RDKitUtil.MolToBase64EncodedMolString(
255 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
256 ),
257 CalcStatus,
258 ConfIDs,
259 ConfEnergies,
260 ]
261
262
263 def GenerateMolConformers(Mol, MolNum=None):
264 """Generate conformers for a molecule."""
265
266 if OptionsInfo["SkipForceFieldMinimization"]:
267 return GenerateMolConformersWithoutMinimization(Mol, MolNum)
268 else:
269 return GenerateMolConformersWithMinimization(Mol, MolNum)
270
271
272 def GenerateMolConformersWithoutMinimization(Mol, MolNum=None):
273 """Generate conformers for a molecule without performing minimization."""
274
275 ConfIDs = EmbedMolecule(Mol, MolNum)
276 if not len(ConfIDs):
277 if not OptionsInfo["QuietMode"]:
278 MolName = RDKitUtil.GetMolName(Mol, MolNum)
279 MiscUtil.PrintWarning(
280 "Conformation generation couldn't be performed for molecule %s: Embedding failed...\n" % MolName
281 )
282 return [Mol, False, None, None]
283
284 if OptionsInfo["AlignConformers"]:
285 AllChem.AlignMolConformers(Mol)
286
287 if not OptionsInfo["QuietMode"]:
288 MolName = RDKitUtil.GetMolName(Mol, MolNum)
289 MiscUtil.PrintInfo("\nNumber of conformations generated for %s: %d" % (MolName, len(ConfIDs)))
290
291 # Convert ConfIDs into a list...
292 ConfIDsList = [ConfID for ConfID in ConfIDs]
293
294 # Setup conformation energies...
295 ConfEnergies = None
296 if OptionsInfo["EnergyOut"]:
297 ConfEnergies = []
298 for ConfID in ConfIDsList:
299 EnergyStatus, Energy = GetConformerEnergy(Mol, ConfID)
300
301 Energy = "%.2f" % Energy if EnergyStatus else "NotAvailable"
302 ConfEnergies.append(Energy)
303
304 if not EnergyStatus:
305 if not OptionsInfo["QuietMode"]:
306 MolName = RDKitUtil.GetMolName(Mol, MolNum)
307 MiscUtil.PrintWarning(
308 "Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n"
309 % (ConfID, MolName)
310 )
311
312 return [Mol, True, ConfIDsList, ConfEnergies]
313
314
315 def GenerateMolConformersWithMinimization(Mol, MolNum):
316 """Generate and mininize conformers for a molecule."""
317
318 if OptionsInfo["AddHydrogens"]:
319 Mol = Chem.AddHs(Mol)
320
321 ConfIDs = EmbedMolecule(Mol, MolNum)
322 if not len(ConfIDs):
323 if not OptionsInfo["QuietMode"]:
324 MolName = RDKitUtil.GetMolName(Mol, MolNum)
325 MiscUtil.PrintWarning(
326 "Conformation generation couldn't be performed for molecule %s: Embedding failed...\n" % MolName
327 )
328 return [Mol, False, None, None]
329
330 CalcEnergyMap = {}
331 for ConfID in ConfIDs:
332 try:
333 if OptionsInfo["UseUFF"]:
334 Status = AllChem.UFFOptimizeMolecule(Mol, confId=ConfID, maxIters=OptionsInfo["MaxIters"])
335 elif OptionsInfo["UseMMFF"]:
336 Status = AllChem.MMFFOptimizeMolecule(
337 Mol, confId=ConfID, maxIters=OptionsInfo["MaxIters"], mmffVariant=OptionsInfo["MMFFVariant"]
338 )
339 else:
340 MiscUtil.PrintError(
341 "Minimization couldn't be performed: Specified forcefield, %s, is not supported"
342 % OptionsInfo["ForceField"]
343 )
344 except (RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
345 if not OptionsInfo["QuietMode"]:
346 MolName = RDKitUtil.GetMolName(Mol, MolNum)
347 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg))
348 return [Mol, False, None, None]
349
350 EnergyStatus, Energy = GetConformerEnergy(Mol, ConfID)
351 if not EnergyStatus:
352 if not OptionsInfo["QuietMode"]:
353 MolName = RDKitUtil.GetMolName(Mol, MolNum)
354 MiscUtil.PrintWarning(
355 "Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n"
356 % (ConfID, MolName)
357 )
358 return [Mol, False, None, None]
359
360 if Status != 0:
361 if not OptionsInfo["QuietMode"]:
362 MolName = RDKitUtil.GetMolName(Mol, MolNum)
363 MiscUtil.PrintWarning(
364 'Minimization failed to converge for conformation number %d of molecule %s in %d steps. Try using higher value for "--maxIters" option...\n'
365 % (ConfID, MolName, OptionsInfo["MaxIters"])
366 )
367
368 CalcEnergyMap[ConfID] = Energy
369
370 if OptionsInfo["RemoveHydrogens"]:
371 Mol = Chem.RemoveHs(Mol)
372
373 # Align molecules after minimization...
374 if OptionsInfo["AlignConformers"]:
375 AllChem.AlignMolConformers(Mol)
376
377 SortedConfIDs = sorted(ConfIDs, key=lambda ConfID: CalcEnergyMap[ConfID])
378
379 MinEnergyConfID = SortedConfIDs[0]
380 MinConfEnergy = CalcEnergyMap[MinEnergyConfID]
381 EnergyWindow = OptionsInfo["EnergyWindow"]
382
383 EnergyRMSDCutoff = OptionsInfo["EnergyRMSDCutoff"]
384 ApplyEnergyRMSDCutoff = False
385 if EnergyRMSDCutoff > 0:
386 ApplyEnergyRMSDCutoff = True
387
388 EnergyRMSDCutoffLowest = OptionsInfo["EnergyRMSDCutoffModeLowest"]
389 EnergyRMSDCalcModeBest = OptionsInfo["EnergyRMSDCalcModeBest"]
390
391 PreAligned = False
392 if OptionsInfo["AlignConformers"]:
393 PreAligned = True
394
395 RefMol, ProbeMol = [None] * 2
396 if EnergyRMSDCalcModeBest:
397 # Copy molecules for best RMSD calculations to avoid change in the coordinates
398 # of the conformations...
399 RefMol = AllChem.Mol(Mol)
400 ProbeMol = AllChem.Mol(Mol)
401
402 # Track conformers with in the specified energy window from the lowest
403 # energy conformation along with applying RMSD cutoff as needed...
404 #
405 SelectedConfIDs = []
406
407 ConfCount = 0
408 IgnoredByEnergyConfCount = 0
409 IgnoredByRMSDConfCount = 0
410
411 FirstConf = True
412 for ConfID in SortedConfIDs:
413 if FirstConf:
414 FirstConf = False
415 ConfCount += 1
416 SelectedConfIDs.append(ConfID)
417 continue
418
419 ConfEnergyDiff = abs(CalcEnergyMap[ConfID] - MinConfEnergy)
420 if ConfEnergyDiff > EnergyWindow:
421 IgnoredByEnergyConfCount += 1
422 continue
423
424 if ApplyEnergyRMSDCutoff:
425 IgnoreConf = False
426 if EnergyRMSDCutoffLowest:
427 # Compare RMSD with the lowest energy conformation...
428 if EnergyRMSDCalcModeBest:
429 CalcRMSD = AllChem.GetBestRMS(ProbeMol, RefMol, prbId=ConfID, refId=MinEnergyConfID)
430 else:
431 CalcRMSD = AllChem.GetConformerRMS(Mol, MinEnergyConfID, ConfID, prealigned=PreAligned)
432 if CalcRMSD < EnergyRMSDCutoff:
433 IgnoreConf = True
434 else:
435 for SelectedConfID in SelectedConfIDs:
436 if EnergyRMSDCalcModeBest:
437 CalcRMSD = AllChem.GetBestRMS(ProbeMol, RefMol, prbId=ConfID, refId=SelectedConfID)
438 else:
439 CalcRMSD = AllChem.GetConformerRMS(Mol, SelectedConfID, ConfID, prealigned=PreAligned)
440 if CalcRMSD < EnergyRMSDCutoff:
441 IgnoreConf = True
442 break
443 if IgnoreConf:
444 IgnoredByRMSDConfCount += 1
445 continue
446
447 ConfCount += 1
448 SelectedConfIDs.append(ConfID)
449
450 if not OptionsInfo["QuietMode"]:
451 MolName = RDKitUtil.GetMolName(Mol, MolNum)
452 MiscUtil.PrintInfo("\nTotal Number of conformations generated for %s: %d" % (MolName, ConfCount))
453 MiscUtil.PrintInfo(
454 "Number of conformations ignored due to energy window cutoff: %d" % (IgnoredByEnergyConfCount)
455 )
456 if ApplyEnergyRMSDCutoff:
457 MiscUtil.PrintInfo(
458 "Number of conformations ignored due to energy RMSD cutoff: %d" % (IgnoredByRMSDConfCount)
459 )
460
461 SelectedConfEnergies = None
462 if OptionsInfo["EnergyOut"]:
463 SelectedConfEnergies = ["%.2f" % CalcEnergyMap[ConfID] for ConfID in SelectedConfIDs]
464
465 return [Mol, True, SelectedConfIDs, SelectedConfEnergies]
466
467
468 def GetConformerEnergy(Mol, ConfID):
469 """Calculate conformer energy."""
470
471 Status = True
472 Energy = 0.0
473
474 if OptionsInfo["UseUFF"]:
475 UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId=ConfID)
476 if UFFMoleculeForcefield is None:
477 Status = False
478 else:
479 Energy = UFFMoleculeForcefield.CalcEnergy()
480 elif OptionsInfo["UseMMFF"]:
481 MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant=OptionsInfo["MMFFVariant"])
482 MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId=ConfID)
483 if MMFFMoleculeForcefield is None:
484 Status = False
485 else:
486 Energy = MMFFMoleculeForcefield.CalcEnergy()
487 else:
488 MiscUtil.PrintError(
489 "Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"]
490 )
491
492 return (Status, Energy)
493
494
495 def EmbedMolecule(Mol, MolNum=None):
496 """Embed conformations."""
497
498 ConfIDs = []
499
500 # Figure out the number of conformations to embded...
501 if re.match("^Auto$", OptionsInfo["MaxConfs"], re.I):
502 NumOfRotBonds = Descriptors.NumRotatableBonds(Mol)
503 if NumOfRotBonds <= 5:
504 MaxConfs = 100
505 elif NumOfRotBonds >= 6 and NumOfRotBonds <= 10:
506 MaxConfs = 200
507 else:
508 MaxConfs = 300
509 else:
510 MaxConfs = int(OptionsInfo["MaxConfs"])
511
512 RandomSeed = OptionsInfo["RandomSeed"]
513 EnforceChirality = OptionsInfo["EnforceChirality"]
514 UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
515 UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
516 ETVersion = OptionsInfo["ETVersion"]
517 EmbedRMSDCutoff = OptionsInfo["EmbedRMSDCutoff"]
518
519 try:
520 ConfIDs = AllChem.EmbedMultipleConfs(
521 Mol,
522 numConfs=MaxConfs,
523 randomSeed=RandomSeed,
524 pruneRmsThresh=EmbedRMSDCutoff,
525 enforceChirality=EnforceChirality,
526 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs,
527 useBasicKnowledge=UseBasicKnowledge,
528 ETversion=ETVersion,
529 )
530 except ValueError as ErrMsg:
531 if not OptionsInfo["QuietMode"]:
532 MolName = RDKitUtil.GetMolName(Mol, MolNum)
533 MiscUtil.PrintWarning("Embedding failed for molecule %s:\n%s\n" % (MolName, ErrMsg))
534 ConfIDs = []
535
536 if not OptionsInfo["QuietMode"]:
537 if EmbedRMSDCutoff > 0:
538 MiscUtil.PrintInfo(
539 "Generating initial conformation ensemble by distance geometry for %s - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s"
540 % (RDKitUtil.GetMolName(Mol, MolNum), EmbedRMSDCutoff, MaxConfs, len(ConfIDs))
541 )
542 else:
543 MiscUtil.PrintInfo(
544 "Generating initial conformation ensemble by distance geometry for %s - EmbedRMSDCutoff: None; Size: %s"
545 % (RDKitUtil.GetMolName(Mol, MolNum), len(ConfIDs))
546 )
547
548 return ConfIDs
549
550
551 def WriteMolConformers(Writer, Mol, MolNum, ConfIDs, ConfEnergies=None):
552 """Write molecule coformers."""
553
554 if ConfIDs is None:
555 return
556
557 MolName = RDKitUtil.GetMolName(Mol, MolNum)
558
559 for Index, ConfID in enumerate(ConfIDs):
560 SetConfMolName(Mol, MolName, ConfID)
561
562 if ConfEnergies is not None:
563 Mol.SetProp(OptionsInfo["EnergyLabel"], ConfEnergies[Index])
564
565 Writer.write(Mol, confId=ConfID)
566
567
568 def SetConfMolName(Mol, MolName, ConfCount):
569 """Set conf mol name."""
570
571 ConfName = "%s_Conf%d" % (MolName, ConfCount)
572 Mol.SetProp("_Name", ConfName)
573
574
575 def ProcesssConformerGeneratorOption():
576 """Process comformer generator option."""
577
578 ConfGenParams = MiscUtil.ProcessOptionConformerGenerator("--conformerGenerator", Options["--conformerGenerator"])
579
580 OptionsInfo["ConformerGenerator"] = ConfGenParams["ConformerGenerator"]
581 OptionsInfo["UseBasicKnowledge"] = ConfGenParams["UseBasicKnowledge"]
582 OptionsInfo["UseExpTorsionAnglePrefs"] = ConfGenParams["UseExpTorsionAnglePrefs"]
583 OptionsInfo["ETVersion"] = ConfGenParams["ETVersion"]
584
585
586 def ProcessOptions():
587 """Process and validate command line arguments and options."""
588
589 MiscUtil.PrintInfo("Processing options...")
590
591 # Validate options...
592 ValidateOptions()
593
594 OptionsInfo["Infile"] = Options["--infile"]
595 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
596 "--infileParams", Options["--infileParams"], Options["--infile"]
597 )
598
599 OptionsInfo["Outfile"] = Options["--outfile"]
600 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
601 "--outfileParams", Options["--outfileParams"]
602 )
603
604 OptionsInfo["Overwrite"] = Options["--overwrite"]
605
606 OptionsInfo["AddHydrogens"] = True
607 if re.match("^no$", Options["--addHydrogens"], re.I):
608 OptionsInfo["AddHydrogens"] = False
609
610 OptionsInfo["AlignConformers"] = True
611 if re.match("^no$", Options["--alignConformers"], re.I):
612 OptionsInfo["AlignConformers"] = False
613
614 ProcesssConformerGeneratorOption()
615
616 if re.match("^UFF$", Options["--forceField"], re.I):
617 ForceField = "UFF"
618 UseUFF = True
619 UseMMFF = False
620 SkipForceFieldMinimization = False
621 elif re.match("^MMFF$", Options["--forceField"], re.I):
622 ForceField = "MMFF"
623 UseUFF = False
624 UseMMFF = True
625 SkipForceFieldMinimization = False
626 else:
627 ForceField = "None"
628 UseUFF = False
629 UseMMFF = False
630 SkipForceFieldMinimization = True
631
632 MMFFVariant = "MMFF94" if re.match("^MMFF94$", Options["--forceFieldMMFFVariant"], re.I) else "MMFF94s"
633
634 OptionsInfo["SkipForceFieldMinimization"] = SkipForceFieldMinimization
635 OptionsInfo["ForceField"] = ForceField
636 OptionsInfo["MMFFVariant"] = MMFFVariant
637 OptionsInfo["UseMMFF"] = UseMMFF
638 OptionsInfo["UseUFF"] = UseUFF
639
640 OptionsInfo["EnergyOut"] = True if re.match("^yes$", Options["--energyOut"], re.I) else False
641 if UseUFF:
642 EnergyLabel = "UFF_Energy"
643 elif UseMMFF:
644 EnergyLabel = "%s_Energy" % MMFFVariant
645 else:
646 EnergyLabel = "Energy"
647 OptionsInfo["EnergyLabel"] = EnergyLabel
648
649 OptionsInfo["EnforceChirality"] = True
650 if re.match("^no$", Options["--enforceChirality"], re.I):
651 OptionsInfo["EnforceChirality"] = False
652
653 OptionsInfo["EnergyWindow"] = float(Options["--energyWindow"])
654
655 EmbedRMSDCutoff = -1.0
656 if not re.match("^none$", Options["--embedRMSDCutoff"], re.I):
657 EmbedRMSDCutoff = float(Options["--embedRMSDCutoff"])
658 OptionsInfo["EmbedRMSDCutoff"] = EmbedRMSDCutoff
659
660 OptionsInfo["EnergyRMSDCalcMode"] = Options["--energyRMSDCalcMode"]
661 OptionsInfo["EnergyRMSDCalcModeBest"] = (
662 True if re.match("^BestRMSD$", Options["--energyRMSDCalcMode"], re.I) else False
663 )
664
665 EnergyRMSDCutoff = -1.0
666 if not re.match("^none$", Options["--energyRMSDCutoff"], re.I):
667 EnergyRMSDCutoff = float(Options["--energyRMSDCutoff"])
668 OptionsInfo["EnergyRMSDCutoff"] = EnergyRMSDCutoff
669
670 OptionsInfo["EnergyRMSDCutoffMode"] = Options["--energyRMSDCutoffMode"]
671 OptionsInfo["EnergyRMSDCutoffModeLowest"] = (
672 True if re.match("^Lowest$", Options["--energyRMSDCutoffMode"], re.I) else False
673 )
674
675 OptionsInfo["MaxIters"] = int(Options["--maxIters"])
676 OptionsInfo["MaxConfs"] = Options["--maxConfs"]
677
678 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
679 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
680
681 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
682
683 RandomSeed = -1
684 if not re.match("^auto$", Options["--randomSeed"], re.I):
685 RandomSeed = int(Options["--randomSeed"])
686 OptionsInfo["RandomSeed"] = RandomSeed
687
688 OptionsInfo["RemoveHydrogens"] = True
689 if re.match("^no$", Options["--removeHydrogens"], re.I):
690 OptionsInfo["RemoveHydrogens"] = False
691
692
693 def RetrieveOptions():
694 """Retrieve command line arguments and options."""
695
696 # Get options...
697 global Options
698 Options = docopt(_docoptUsage_)
699
700 # Set current working directory to the specified directory...
701 WorkingDir = Options["--workingdir"]
702 if WorkingDir:
703 os.chdir(WorkingDir)
704
705 # Handle examples option...
706 if "--examples" in Options and Options["--examples"]:
707 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
708 sys.exit(0)
709
710
711 def ValidateOptions():
712 """Validate option values."""
713
714 MiscUtil.ValidateOptionTextValue("-a, --addHydrogens", Options["--addHydrogens"], "yes no")
715 MiscUtil.ValidateOptionTextValue("--alignConformers", Options["--alignConformers"], "yes no")
716 MiscUtil.ValidateOptionTextValue(
717 "-c, --conformerGenerator", Options["--conformerGenerator"], "SDG KDG ETDG ETKDG ETKDGv2"
718 )
719
720 MiscUtil.ValidateOptionTextValue("-f, --forceField", Options["--forceField"], "UFF MMFF None")
721 MiscUtil.ValidateOptionTextValue(" --forceFieldMMFFVariant", Options["--forceFieldMMFFVariant"], "MMFF94 MMFF94s")
722
723 MiscUtil.ValidateOptionTextValue("--energyOut", Options["--energyOut"], "yes no")
724 MiscUtil.ValidateOptionTextValue("--enforceChirality ", Options["--enforceChirality"], "yes no")
725 MiscUtil.ValidateOptionFloatValue("--energyWindow", Options["--energyWindow"], {">": 0.0})
726
727 MiscUtil.ValidateOptionTextValue(" --energyRMSDCalcMode", Options["--energyRMSDCalcMode"], "RMSD BestRMSD")
728
729 if not re.match("^none$", Options["--embedRMSDCutoff"], re.I):
730 MiscUtil.ValidateOptionFloatValue("--embedRMSDCutoff", Options["--embedRMSDCutoff"], {">": 0.0})
731 MiscUtil.ValidateOptionTextValue(" --energyRMSDCutoffMode", Options["--energyRMSDCutoffMode"], "All Lowest")
732
733 if not re.match("^none$", Options["--energyRMSDCutoff"], re.I):
734 MiscUtil.ValidateOptionFloatValue("--energyRMSDCutoff", Options["--energyRMSDCutoff"], {">": 0.0})
735 # Make sure that the alignConformers option is being used...
736 if not re.match("^yes$", Options["--alignConformers"], re.I):
737 MiscUtil.PrintError(
738 '"%s" value of "--alignConformers" is not allowed for %s value of "--energyRMSDCutoff" option. It must be set to "yes".'
739 % (Options["--alignConformers"], Options["--energyRMSDCutoff"])
740 )
741
742 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
743 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
744
745 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
746 MiscUtil.ValidateOptionsOutputFileOverwrite(
747 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
748 )
749 MiscUtil.ValidateOptionsDistinctFileNames(
750 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
751 )
752
753 if not re.match("^auto$", Options["--maxConfs"], re.I):
754 MiscUtil.ValidateOptionIntegerValue("--maxConfs", Options["--maxConfs"], {">": 0})
755
756 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
757
758 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
759 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
760
761 if not re.match("^auto$", Options["--randomSeed"], re.I):
762 MiscUtil.ValidateOptionIntegerValue("--randomSeed", Options["--randomSeed"], {})
763
764 MiscUtil.ValidateOptionTextValue("-r, --removeHydrogens", Options["--removeHydrogens"], "yes no")
765
766
767 # Setup a usage string for docopt...
768 _docoptUsage_ = """
769 RDKitGenerateConformers.py - Generate molecular conformations
770
771 Usage:
772 RDKitGenerateConformers.py [--alignConformers <yes or no>] [--addHydrogens <yes or no>]
773 [--conformerGenerator <SDG, ETDG, KDG, ETKDG, ETKDGv2>] [--embedRMSDCutoff <number>]
774 [--energyOut <yes or no>] [--enforceChirality <yes or no>] [--energyRMSDCalcMode <RMSD or BestRMSD>]
775 [--energyRMSDCutoff <number>] [--energyRMSDCutoffMode <All or Lowest>] [--energyWindow <number>]
776 [--forceField <UFF, MMFF, None>] [--forceFieldMMFFVariant <MMFF94 or MMFF94s>]
777 [--infileParams <Name,Value,...>] [--maxConfs <number>]
778 [--mp <yes or no>] [--mpParams <Name,Value,...>]
779 [--maxIters <number>] [ --outfileParams <Name,Value,...> ] [--overwrite]
780 [--quiet <yes or no>] [ --removeHydrogens <yes or no>] [--randomSeed <number>]
781 [-w <dir>] -i <infile> -o <outfile>
782 RDKitGenerateConformers.py -h | --help | -e | --examples
783
784 Description:
785 Generate 3D conformes of molecules using a combination of distance geometry and
786 forcefield minimization. The forcefield minimization may be skipped to only generate
787 conformations by available distance geometry based methodologies.
788
789 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
790 .csv, .tsv, .txt)
791
792 The supported output file format are: SD (.sdf, .sd)
793
794 Options:
795 -a, --addHydrogens <yes or no> [default: yes]
796 Add hydrogens before minimization.
797 --alignConformers <yes or no> [default: yes]
798 Align conformers for each molecule.
799 -c, --conformerGenerator <text> [default: ETKDGv2]
800 Conformation generation methodology for generating initial 3D coordinates. The
801 possible values along with a brief description are shown below:
802
803 SDG: Standard Distance Geometry
804 KDG: basic Knowledge-terms with Distance Geometry
805 ETDG: Experimental Torsion-angle preference with Distance Geometry
806 ETKDG: Experimental Torsion-angle preference along with basic
807 Knowledge-terms and Distance Geometry [Ref 129]
808 ETKDGv2: Experimental Torsion-angle preference along with basic
809 Knowledge-terms and Distance Geometry [Ref 167]
810 --embedRMSDCutoff <number> [default: none]
811 RMSD cutoff for retaining conformations after embedding and before energy minimization.
812 All embedded conformations are kept by default. Otherwise, only those conformations
813 which are different from each other by the specified RMSD cutoff are kept. The first
814 embedded conformation is always retained.
815 --energyOut <yes or no> [default: No]
816 Write out energy values.
817 --enforceChirality <yes or no> [default: Yes]
818 Enforce chirality for defined chiral centers.
819 --energyRMSDCalcMode <RMSD or BestRMSD> [default: RMSD]
820 Methodology for calculating RMSD values during the application of RMSD
821 cutoff for retaining conformations after energy minimization. Possible
822 values: RMSD or BestRMSD. This option is ignore during 'None' value of
823 '--energyRMSDCutoff' option.
824
825 During BestRMSMode mode, the RDKit 'function AllChem.GetBestRMS' is used to
826 align and calculate RMSD. This function calculates optimal RMSD for aligning two
827 molecules, taking symmetry into account. Otherwise, the RMSD value is calculated
828 using 'AllChem.GetConformerRMS' without changing the atom order. A word to the
829 wise from RDKit documentation: The AllChem.GetBestRMS function will attempt to
830 align all permutations of matching atom orders in both molecules, for some molecules
831 it will lead to 'combinatorial explosion'.
832 --energyRMSDCutoff <number> [default: none]
833 RMSD cutoff for retaining conformations after energy minimization. By default,
834 all minimized conformations with in the specified energy window from the lowest
835 energy conformation are kept. Otherwise, only those conformations which are
836 different from the lowest energy conformation or all selected conformations
837 by the specified RMSD cutoff and are with in the specified energy window are
838 kept. The lowest energy conformation is always retained.
839 --energyRMSDCutoffMode <All or Lowest> [default: All]
840 RMSD cutoff mode for retaining conformations after energy minimization.
841 Possible values: All or Lowest. The RMSD values are compared against all
842 the selected conformations or the lowest energy conformation during 'All'
843 and 'Lowest' value of '--energyRMSDCutoffMode'. This option is ignored
844 during 'None' value of '--energyRMSDCutoff' option.
845
846 By default, only those conformations which all different from all selected
847 conformations by the specified RMSD cutoff and are with in the specified
848 energy window are kept.
849 --energyWindow <number> [default: 20]
850 Energy window in kcal/mol for selecting conformers. This option is ignored during
851 'None' value of '-f, --forcefield' option.
852 -e, --examples
853 Print examples.
854 -f, --forceField <UFF, MMFF, None> [default: MMFF]
855 Forcefield method to use for energy minimization. Possible values: Universal Force
856 Field (UFF) [Ref 81], Merck Molecular Mechanics Force Field (MMFF) [Ref 83-87] or
857 None.
858 --forceFieldMMFFVariant <MMFF94 or MMFF94s> [default: MMFF94]
859 Variant of MMFF forcefield to use for energy minimization.
860 -h, --help
861 Print this help message.
862 -i, --infile <infile>
863 Input file name.
864 --infileParams <Name,Value,...> [default: auto]
865 A comma delimited list of parameter name and value pairs for reading
866 molecules from files. The supported parameter names for different file
867 formats, along with their default values, are shown below:
868
869 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
870 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
871 smilesTitleLine,auto,sanitize,yes
872
873 Possible values for smilesDelimiter: space, comma or tab.
874 --maxConfs <number> [default: auto]
875 Maximum number of conformations to generate for each molecule by conformation
876 generation methodology. The conformations are minimized using the specified
877 forcefield as needed and written to the output file. The default value for maximum
878 number of conformations is dependent on the number of rotatable bonds in molecules:
879 RotBonds <= 5, maxConfs = 100; RotBonds >=6 and <= 10, MaxConfs = 200;
880 RotBonds >= 11, maxConfs = 300
881 --maxIters <number> [default: 250]
882 Maximum number of iterations to perform for each molecule during forcefield
883 minimization.
884 --mp <yes or no> [default: no]
885 Use multiprocessing.
886
887 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
888 function employing lazy RDKit data iterable. This allows processing of
889 arbitrary large data sets without any additional requirements memory.
890
891 All input data may be optionally loaded into memory by mp.Pool.map()
892 before starting worker processes in a process pool by setting the value
893 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
894
895 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
896 data mode may adversely impact the performance. The '--mpParams' section
897 provides additional information to tune the value of 'chunkSize'.
898 --mpParams <Name,Value,...> [default: auto]
899 A comma delimited list of parameter name and value pairs to configure
900 multiprocessing.
901
902 The supported parameter names along with their default and possible
903 values are shown below:
904
905 chunkSize, auto
906 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
907 numProcesses, auto [ Default: mp.cpu_count() ]
908
909 These parameters are used by the following functions to configure and
910 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
911 mp.Pool.imap().
912
913 The chunkSize determines chunks of input data passed to each worker
914 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
915 The default value of chunkSize is dependent on the value of 'inputDataMode'.
916
917 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
918 automatically converts RDKit data iterable into a list, loads all data into
919 memory, and calculates the default chunkSize using the following method
920 as shown in its code:
921
922 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
923 if extra: chunkSize += 1
924
925 For example, the default chunkSize will be 7 for a pool of 4 worker processes
926 and 100 data items.
927
928 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
929 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
930 data into memory. Consequently, the size of input data is not known a priori.
931 It's not possible to estimate an optimal value for the chunkSize. The default
932 chunkSize is set to 1.
933
934 The default value for the chunkSize during 'Lazy' data mode may adversely
935 impact the performance due to the overhead associated with exchanging
936 small chunks of data. It is generally a good idea to explicitly set chunkSize to
937 a larger value during 'Lazy' input data mode, based on the size of your input
938 data and number of processes in the process pool.
939
940 The mp.Pool.map() function waits for all worker processes to process all
941 the data and return the results. The mp.Pool.imap() function, however,
942 returns the the results obtained from worker processes as soon as the
943 results become available for specified chunks of data.
944
945 The order of data in the results returned by both mp.Pool.map() and
946 mp.Pool.imap() functions always corresponds to the input data.
947 -o, --outfile <outfile>
948 Output file name.
949 --outfileParams <Name,Value,...> [default: auto]
950 A comma delimited list of parameter name and value pairs for writing
951 molecules to files. The supported parameter names for different file
952 formats, along with their default values, are shown below:
953
954 SD: kekulize,yes,forceV3000,no
955
956 --overwrite
957 Overwrite existing files.
958 -q, --quiet <yes or no> [default: no]
959 Use quiet mode. The warning and information messages will not be printed.
960 -r, --removeHydrogens <yes or no> [default: Yes]
961 Remove hydrogens after minimization.
962 --randomSeed <number> [default: auto]
963 Seed for the random number generator for reproducing 3D coordinates.
964 Default is to use a random seed.
965 -w, --workingdir <dir>
966 Location of working directory which defaults to the current directory.
967
968 Examples:
969 To generate conformers using Experimental Torsion-angle preference along
970 with basic Knowledge-terms and Distance Geometry (ETKDG) followed by
971 MMFF minimization with automatic determination of maximum number of
972 conformers for each molecule and write out a SD file, type:
973
974 % RDKitGenerateConformers.py -i Sample.smi -o SampleOut.sdf
975
976 To rerun the first example in a quiet mode and write out a SD file, type:
977
978 % RDKitGenerateConformers.py -q yes -i Sample.smi -o SampleOut.sdf
979
980 To rerun the first example in multiprocessing mode on all available CPUs
981 without loading all data into memory and write out a SD file, type:
982
983 % RDKitGenerateConformers.py --mp yes -i Sample.smi -o SampleOut.sdf
984
985 To run the first example in multiprocessing mode on all available CPUs
986 by loading all data into memory and write out a SD file, type:
987
988 % RDKitGenerateConformers.py --mp yes --mpParams "inputDataMode,
989 InMemory" -i Sample.smi -o SampleOut.sdf
990
991 To rerun the first example in multiprocessing mode on specific number of
992 CPUs and chunk size without loading all data into memory and write out a SD file,
993 type:
994
995 % RDKitGenerateConformers.py --mp yes --mpParams "inputDataMode,Lazy,
996 numProcesses,4,chunkSize,8" -i Sample.smi -o SampleOut.sdf
997
998 To generate up to 150 conformers for each molecule using ETKDG and UFF forcefield
999 minimization along with conformers within 25 kcal/mol energy window and write out a
1000 SD file, type:
1001
1002 % RDKitGenerateConformers.py --energyWindow 25 -f UFF --maxConfs 150
1003 -i Sample.smi -o SampleOut.sdf
1004
1005 To generate up to 50 conformers for each molecule using KDG without any forcefield
1006 minimization and alignment of conformers and write out a SD file, type:
1007
1008 % RDKitGenerateConformers.py -f none --maxConfs 50 --alignConformers no
1009 -i Sample.sdf -o SampleOut.sdf
1010
1011 To generate up to 50 conformers using SDG without any forcefield minimization
1012 and alignment of conformers for molecules in a CSV SMILES file, SMILES strings
1013 in column 1, name in column 2, and write out a SD file, type:
1014
1015 % RDKitGenerateConformers.py --maxConfs 50 --maxIters 50 -c SDG
1016 --alignConformers no -f none --infileParams "smilesDelimiter,comma,
1017 smilesTitleLine,yes, smilesColumn,1,smilesNameColumn,2"
1018 -i SampleSMILES.csv -o SampleOut.sdf
1019
1020 Author:
1021 Manish Sud(msud@san.rr.com)
1022
1023 See also:
1024 RDKitCalculateRMSD.py, RDKitCalculateMolecularDescriptors.py,
1025 RDKitCompareMoleculeShapes.py, RDKitConvertFileFormat.py,
1026 RDKitGenerateConstrainedConformers.py, RDKitPerformMinimization.py
1027
1028 Copyright:
1029 Copyright (C) 2026 Manish Sud. All rights reserved.
1030
1031 The functionality available in this script is implemented using RDKit, an
1032 open source toolkit for cheminformatics developed by Greg Landrum.
1033
1034 This file is part of MayaChemTools.
1035
1036 MayaChemTools is free software; you can redistribute it and/or modify it under
1037 the terms of the GNU Lesser General Public License as published by the Free
1038 Software Foundation; either version 3 of the License, or (at your option) any
1039 later version.
1040
1041 """
1042
1043 if __name__ == "__main__":
1044 main()
