1 #!/bin/env python
2 #
3 # File: RDKitGenerateConstrainedConformers.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 rdFMCS
43 from rdkit.Chem import rdMolAlign
44 except ImportError as ErrMsg:
45 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
46 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
47 sys.exit(1)
48
49 # MayaChemTools imports...
50 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
51 try:
52 from docopt import docopt
53 import MiscUtil
54 import RDKitUtil
55 except ImportError as ErrMsg:
56 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
57 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
58 sys.exit(1)
59
60 ScriptName = os.path.basename(sys.argv[0])
61 Options = {}
62 OptionsInfo = {}
63
64
65 def main():
66 """Start execution of the script."""
67
68 MiscUtil.PrintInfo(
69 "\n%s (RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
70 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
71 )
72
73 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
74
75 # Retrieve command line arguments and options...
76 RetrieveOptions()
77
78 # Process and validate command line arguments and options...
79 ProcessOptions()
80
81 # Perform actions required by the script...
82 GenerateConstrainedConformers()
83
84 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
85 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
86
87
88 def GenerateConstrainedConformers():
89 """Generate constrained conformers."""
90
91 # Read and validate reference molecule...
92 RefMol = RetrieveReferenceMolecule()
93
94 # Setup a molecule reader for input file...
95 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
96 OptionsInfo["InfileParams"]["AllowEmptyMols"] = True
97 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
98
99 # Set up a molecule writer...
100 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
101 if Writer is None:
102 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
103 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
104
105 MolCount, ValidMolCount, CoreScaffoldMissingCount, ConfGenFailedCount = ProcessMolecules(RefMol, Mols, Writer)
106
107 if Writer is not None:
108 Writer.close()
109
110 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
111 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
112 MiscUtil.PrintInfo("Number of molecules with missing core scaffold: %d" % CoreScaffoldMissingCount)
113 MiscUtil.PrintInfo(
114 "Number of molecules failed during conformation generation or minimization: %d" % ConfGenFailedCount
115 )
116 MiscUtil.PrintInfo(
117 "Number of ignored molecules: %d" % (MolCount - ValidMolCount + CoreScaffoldMissingCount + ConfGenFailedCount)
118 )
119
120
121 def ProcessMolecules(RefMol, Mols, Writer):
122 """Process molecules to generate constrained conformers."""
123
124 if OptionsInfo["MPMode"]:
125 return ProcessMoleculesUsingMultipleProcesses(RefMol, Mols, Writer)
126 else:
127 return ProcessMoleculesUsingSingleProcess(RefMol, Mols, Writer)
128
129
130 def ProcessMoleculesUsingSingleProcess(RefMol, Mols, Writer):
131 """Process molecules to generate constrained conformers using a single process."""
132
133 (MolCount, ValidMolCount, CoreScaffoldMissingCount, ConfGenFailedCount) = [0] * 4
134
135 for Mol in Mols:
136 MolCount += 1
137
138 if Mol is None:
139 continue
140
141 if RDKitUtil.IsMolEmpty(Mol):
142 if not OptionsInfo["QuietMode"]:
143 MolName = RDKitUtil.GetMolName(Mol, MolCount)
144 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
145 continue
146 ValidMolCount += 1
147
148 # Setup a reference molecule core containing common scaffold atoms...
149 RefMolCore = SetupCoreScaffold(RefMol, Mol, MolCount)
150 if RefMolCore is None:
151 CoreScaffoldMissingCount += 1
152 continue
153
154 ConfMols, CalcStatus, ConfIDs, ConfEnergyValues, ConfScaffoldEmbedRMSDValues = GenerateMolConformers(
155 Mol, RefMolCore, MolCount
156 )
157
158 if not CalcStatus:
159 ConfGenFailedCount += 1
160 continue
161
162 WriteMolConformers(Writer, Mol, MolCount, ConfMols, ConfIDs, ConfEnergyValues, ConfScaffoldEmbedRMSDValues)
163
164 return (MolCount, ValidMolCount, CoreScaffoldMissingCount, ConfGenFailedCount)
165
166
167 def ProcessMoleculesUsingMultipleProcesses(RefMol, Mols, Writer):
168 """Process molecules to generate constrained conformers using multiprocessing."""
169
170 MPParams = OptionsInfo["MPParams"]
171
172 # Setup data for initializing a worker process...
173 MiscUtil.PrintInfo("Encoding options info and reference molecule...")
174
175 OptionsInfo["EncodedRefMol"] = RDKitUtil.MolToBase64EncodedMolString(RefMol)
176 InitializeWorkerProcessArgs = (
177 MiscUtil.ObjectToBase64EncodedString(Options),
178 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
179 )
180
181 # Setup a encoded mols data iterable for a worker process...
182 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
183
184 # Setup process pool along with data initialization for each process...
185 MiscUtil.PrintInfo(
186 "\nConfiguring multiprocessing using %s method..."
187 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
188 )
189 MiscUtil.PrintInfo(
190 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
191 % (
192 MPParams["NumProcesses"],
193 MPParams["InputDataMode"],
194 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
195 )
196 )
197
198 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
199
200 # Start processing...
201 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
202 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
203 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
204 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
205 else:
206 MiscUtil.PrintError(
207 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
208 )
209
210 (MolCount, ValidMolCount, CoreScaffoldMissingCount, ConfGenFailedCount) = [0] * 4
211 for Result in Results:
212 MolCount += 1
213 (
214 MolIndex,
215 EncodedMol,
216 EncodedConfMols,
217 CoreScaffoldMissingStatus,
218 CalcStatus,
219 ConfIDs,
220 ConfEnergyValues,
221 ConfScaffoldEmbedRMSDValues,
222 ) = Result
223
224 if EncodedMol is None:
225 continue
226 ValidMolCount += 1
227
228 if CoreScaffoldMissingStatus:
229 CoreScaffoldMissingCount += 1
230 continue
231
232 if not CalcStatus:
233 ConfGenFailedCount += 1
234 continue
235
236 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
237 ConfMols = [RDKitUtil.MolFromBase64EncodedMolString(EncodedConfMol) for EncodedConfMol in EncodedConfMols]
238
239 WriteMolConformers(Writer, Mol, MolCount, ConfMols, ConfIDs, ConfEnergyValues, ConfScaffoldEmbedRMSDValues)
240
241 return (MolCount, ValidMolCount, CoreScaffoldMissingCount, ConfGenFailedCount)
242
243
244 def InitializeWorkerProcess(*EncodedArgs):
245 """Initialize data for a worker process."""
246
247 global Options, OptionsInfo
248
249 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
250
251 # Decode Options and OptionInfo...
252 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
253 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
254
255 # Decode RefMol...
256 OptionsInfo["RefMol"] = RDKitUtil.MolFromBase64EncodedMolString(OptionsInfo["EncodedRefMol"])
257
258
259 def WorkerProcess(EncodedMolInfo):
260 """Process data for a worker process."""
261
262 MolIndex, EncodedMol = EncodedMolInfo
263
264 ConfMols = None
265 CoreScaffoldMissingStatus = False
266 CalcStatus = False
267 ConfIDs = None
268 ConfEnergyValues = None
269 ConfScaffoldEmbedRMSDValues = None
270
271 if EncodedMol is None:
272 return [
273 MolIndex,
274 None,
275 ConfMols,
276 CoreScaffoldMissingStatus,
277 CalcStatus,
278 ConfIDs,
279 ConfEnergyValues,
280 ConfScaffoldEmbedRMSDValues,
281 ]
282
283 RefMol = OptionsInfo["RefMol"]
284
285 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
286 if RDKitUtil.IsMolEmpty(Mol):
287 if not OptionsInfo["QuietMode"]:
288 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
289 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
290 return [
291 MolIndex,
292 None,
293 ConfMols,
294 CoreScaffoldMissingStatus,
295 CalcStatus,
296 ConfIDs,
297 ConfEnergyValues,
298 ConfScaffoldEmbedRMSDValues,
299 ]
300
301 # Setup a reference molecule core containing common scaffold atoms...
302 RefMolCore = SetupCoreScaffold(RefMol, Mol, (MolIndex + 1))
303 if RefMolCore is None:
304 CoreScaffoldMissingStatus = True
305 return [
306 MolIndex,
307 EncodedMol,
308 ConfMols,
309 CoreScaffoldMissingStatus,
310 CalcStatus,
311 ConfIDs,
312 ConfEnergyValues,
313 ConfScaffoldEmbedRMSDValues,
314 ]
315
316 ConfMols, CalcStatus, ConfIDs, ConfEnergyValues, ConfScaffoldEmbedRMSDValues = GenerateMolConformers(
317 Mol, RefMolCore, (MolIndex + 1)
318 )
319
320 EncodedConfMols = None
321 if ConfMols is not None:
322 EncodedConfMols = [
323 RDKitUtil.MolToBase64EncodedMolString(
324 ConfMol,
325 PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps,
326 )
327 for ConfMol in ConfMols
328 ]
329
330 return [
331 MolIndex,
332 EncodedMol,
333 EncodedConfMols,
334 CoreScaffoldMissingStatus,
335 CalcStatus,
336 ConfIDs,
337 ConfEnergyValues,
338 ConfScaffoldEmbedRMSDValues,
339 ]
340
341
342 def RetrieveReferenceMolecule():
343 """Retrieve and validate reference molecule."""
344
345 RefFile = OptionsInfo["RefFile"]
346
347 MiscUtil.PrintInfo("\nProcessing file %s..." % (RefFile))
348 OptionsInfo["InfileParams"]["AllowEmptyMols"] = False
349 ValidRefMols, RefMolCount, ValidRefMolCount = RDKitUtil.ReadAndValidateMolecules(
350 RefFile, **OptionsInfo["InfileParams"]
351 )
352
353 if ValidRefMolCount == 0:
354 MiscUtil.PrintError("The reference file, %s, contains no valid molecules." % RefFile)
355 elif ValidRefMolCount > 1:
356 MiscUtil.PrintWarning(
357 "The reference file, %s, contains, %d, valid molecules. Using first molecule as the reference molecule..."
358 % (RefFile, ValidRefMolCount)
359 )
360
361 RefMol = ValidRefMols[0]
362
363 if OptionsInfo["UseScaffoldSMARTS"]:
364 ScaffoldPatternMol = Chem.MolFromSmarts(OptionsInfo["ScaffoldSMARTS"])
365 if ScaffoldPatternMol is None:
366 MiscUtil.PrintError(
367 'Failed to create scaffold pattern molecule. The scaffold SMARTS pattern, %s, specified using "-s, --scaffold" option is not valid.'
368 % (OptionsInfo["ScaffoldSMARTS"])
369 )
370
371 if not RefMol.HasSubstructMatch(ScaffoldPatternMol):
372 MiscUtil.PrintError(
373 'The scaffold SMARTS pattern, %s, specified using "-s, --scaffold" option, is missing in the first valid reference molecule.'
374 % (OptionsInfo["ScaffoldSMARTS"])
375 )
376
377 return RefMol
378
379
380 def SetupCoreScaffold(RefMol, Mol, MolCount):
381 """Setup a reference molecule core containing common scaffold atoms between
382 a pair of molecules."""
383
384 if OptionsInfo["UseScaffoldMCS"]:
385 return SetupCoreScaffoldByMCS(RefMol, Mol, MolCount)
386 elif OptionsInfo["UseScaffoldSMARTS"]:
387 return SetupCoreScaffoldBySMARTS(RefMol, Mol, MolCount)
388 else:
389 MiscUtil.PrintError(
390 'The value, %s, specified for "-s, --scaffold" option is not supported.' % (OptionsInfo["Scaffold"])
391 )
392
393
394 def SetupCoreScaffoldByMCS(RefMol, Mol, MolCount):
395 """Setup a reference molecule core containing common scaffold atoms between
396 a pair of molecules using MCS."""
397
398 MCSParams = OptionsInfo["MCSParams"]
399 Mols = [RefMol, Mol]
400
401 MCSResultObject = rdFMCS.FindMCS(
402 Mols,
403 maximizeBonds=MCSParams["MaximizeBonds"],
404 threshold=MCSParams["Threshold"],
405 timeout=MCSParams["TimeOut"],
406 verbose=MCSParams["Verbose"],
407 matchValences=MCSParams["MatchValences"],
408 ringMatchesRingOnly=MCSParams["RingMatchesRingOnly"],
409 completeRingsOnly=MCSParams["CompleteRingsOnly"],
410 matchChiralTag=MCSParams["MatchChiralTag"],
411 atomCompare=MCSParams["AtomCompare"],
412 bondCompare=MCSParams["BondCompare"],
413 seedSmarts=MCSParams["SeedSMARTS"],
414 )
415
416 if MCSResultObject.canceled:
417 if not OptionsInfo["QuietMode"]:
418 MiscUtil.PrintWarning(
419 'MCS failed to identify a common core scaffold between reference moecule and input molecule %s. Specify a different set of parameters using "-m, --mcsParams" option and try again.'
420 % (RDKitUtil.GetMolName(Mol, MolCount))
421 )
422 return None
423
424 CoreNumAtoms = MCSResultObject.numAtoms
425 CoreNumBonds = MCSResultObject.numBonds
426
427 SMARTSCore = MCSResultObject.smartsString
428
429 if not len(SMARTSCore):
430 if not OptionsInfo["QuietMode"]:
431 MiscUtil.PrintWarning(
432 'MCS failed to identify a common core scaffold between reference moecule and input molecule %s. Specify a different set of parameters using "-m, --mcsParams" option and try again.'
433 % (RDKitUtil.GetMolName(Mol, MolCount))
434 )
435 return None
436
437 if CoreNumAtoms < MCSParams["MinNumAtoms"]:
438 if not OptionsInfo["QuietMode"]:
439 MiscUtil.PrintWarning(
440 'Number of atoms, %d, in core scaffold identified by MCS is less than, %d, as specified by "minNumAtoms" parameter in "-m, --mcsParams" option.'
441 % (CoreNumAtoms, MCSParams["MinNumAtoms"])
442 )
443 return None
444
445 if CoreNumBonds < MCSParams["MinNumBonds"]:
446 if not OptionsInfo["QuietMode"]:
447 MiscUtil.PrintWarning(
448 'Number of bonds, %d, in core scaffold identified by MCS is less than, %d, as specified by "minNumBonds" parameter in "-m, --mcsParams" option.'
449 % (CoreNumBonds, MCSParams["MinNumBonds"])
450 )
451 return None
452
453 return GenerateCoreMol(RefMol, SMARTSCore)
454
455
456 def SetupCoreScaffoldBySMARTS(RefMol, Mol, MolCount):
457 """Setup a reference molecule core containing common scaffold atoms between
458 a pair of molecules using specified SMARTS."""
459
460 if OptionsInfo["ScaffoldPatternMol"] is None:
461 OptionsInfo["ScaffoldPatternMol"] = Chem.MolFromSmarts(OptionsInfo["ScaffoldSMARTS"])
462
463 if not Mol.HasSubstructMatch(OptionsInfo["ScaffoldPatternMol"]):
464 if not OptionsInfo["QuietMode"]:
465 MiscUtil.PrintWarning(
466 'The scaffold SMARTS pattern, %s, specified using "-s, --scaffold" option is missing in input molecule, %s.'
467 % (OptionsInfo["ScaffoldSMARTS"], RDKitUtil.GetMolName(Mol, MolCount))
468 )
469 return None
470
471 return GenerateCoreMol(RefMol, OptionsInfo["ScaffoldSMARTS"])
472
473
474 def GenerateCoreMol(RefMol, SMARTSCore):
475 """Generate core molecule for embedding."""
476
477 # Create a molecule corresponding to core atoms...
478 SMARTSCoreMol = Chem.MolFromSmarts(SMARTSCore)
479
480 # Setup a ref molecule containing core atoms with dummy atoms as
481 # attachment points for atoms around the core atoms...
482 Core = AllChem.ReplaceSidechains(Chem.RemoveHs(RefMol), SMARTSCoreMol)
483
484 # Delete any substructures containing dummy atoms..
485 RefMolCore = AllChem.DeleteSubstructs(Core, Chem.MolFromSmiles("*"))
486 RefMolCore.UpdatePropertyCache()
487
488 return RefMolCore
489
490
491 def GenerateMolConformers(Mol, RefMolCore, MolNum=None):
492 """Generate constrained conformers."""
493
494 if OptionsInfo["AddHydrogens"]:
495 Mol = Chem.AddHs(Mol, addCoords=True)
496
497 # Setup forcefield function to use for constrained minimization...
498 ForceFieldFunction = None
499 ForceFieldName = None
500 if OptionsInfo["UseUFF"]:
501 ForceFieldFunction = lambda mol, confId=-1: AllChem.UFFGetMoleculeForceField(mol, confId=confId)
502 ForceFieldName = "UFF"
503 else:
504 ForceFieldFunction = lambda mol, confId=-1: AllChem.MMFFGetMoleculeForceField(
505 mol, AllChem.MMFFGetMoleculeProperties(mol, mmffVariant=OptionsInfo["MMFFVariant"]), confId=confId
506 )
507 ForceFieldName = "MMFF"
508
509 if ForceFieldFunction is None:
510 if not OptionsInfo["QuietMode"]:
511 MiscUtil.PrintWarning(
512 "Failed to setup forcefield %s for molecule: %s\n" % (ForceFieldName, RDKitUtil.GetMolName(Mol, MolNum))
513 )
514 return (None, False, None, None, None)
515
516 MaxConfs = OptionsInfo["MaxConfs"]
517 EnforceChirality = OptionsInfo["EnforceChirality"]
518 UseExpTorsionAnglePrefs = OptionsInfo["UseExpTorsionAnglePrefs"]
519 ETVersion = OptionsInfo["ETVersion"]
520 UseBasicKnowledge = OptionsInfo["UseBasicKnowledge"]
521 UseTethers = OptionsInfo["UseTethers"]
522
523 CalcEnergyMap = {}
524 MolConfsMap = {}
525 CalcRMSDMap = {}
526
527 ScaffoldEmbedRMSDMap = {}
528
529 ConfIDs = [ConfID for ConfID in range(0, MaxConfs)]
530 for ConfID in ConfIDs:
531 try:
532 MolConf = Chem.Mol(Mol)
533 AllChem.ConstrainedEmbed(
534 MolConf,
535 RefMolCore,
536 useTethers=UseTethers,
537 coreConfId=-1,
538 randomseed=ConfID,
539 getForceField=ForceFieldFunction,
540 enforceChirality=EnforceChirality,
541 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs,
542 useBasicKnowledge=UseBasicKnowledge,
543 ETversion=ETVersion,
544 )
545 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
546 if not OptionsInfo["QuietMode"]:
547 MolName = RDKitUtil.GetMolName(Mol, MolNum)
548 MiscUtil.PrintWarning(
549 "Constrained embedding coupldn't be performed for molecule %s:\n%s\n"
550 % (RDKitUtil.GetMolName(Mol, MolNum), ErrMsg)
551 )
552 return (None, False, None, None, None)
553
554 EnergyStatus, Energy = GetEnergy(MolConf)
555
556 if not EnergyStatus:
557 if not OptionsInfo["QuietMode"]:
558 MolName = RDKitUtil.GetMolName(Mol, MolNum)
559 MiscUtil.PrintWarning(
560 "Failed to retrieve calculated energy for conformation number %d of molecule %s. Try again after removing any salts or cleaing up the molecule...\n"
561 % (ConfID, MolName)
562 )
563 return (None, False, None, None, None)
564
565 CalcEnergyMap[ConfID] = Energy
566
567 if OptionsInfo["ScaffoldRMSDOut"]:
568 ScaffoldEmbedRMSDMap[ConfID] = "%.4f" % float(MolConf.GetProp("EmbedRMS"))
569 MolConf.ClearProp("EmbedRMS")
570
571 if OptionsInfo["RemoveHydrogens"]:
572 MolConf = Chem.RemoveHs(MolConf)
573 MolConfsMap[ConfID] = MolConf
574
575 # Sort conformers by energy...
576 SortedConfIDs = sorted(ConfIDs, key=lambda ConfID: CalcEnergyMap[ConfID])
577
578 MinEnergyConfID = SortedConfIDs[0]
579 EnergyWindow = OptionsInfo["EnergyWindow"]
580
581 MinConfEnergy = CalcEnergyMap[MinEnergyConfID]
582 MinEnergyMolConf = MolConfsMap[MinEnergyConfID]
583
584 # Calculate RMSD values for conformers...
585 CalcRMSDMap = {}
586
587 EnergyRMSDCutoff = OptionsInfo["EnergyRMSDCutoff"]
588 ApplyEnergyRMSDCutoff = True if EnergyRMSDCutoff > 0 else False
589
590 if ApplyEnergyRMSDCutoff:
591 FirstConf = True
592 for ConfID in SortedConfIDs:
593 if FirstConf:
594 FirstConf = False
595 CalcRMSDMap[ConfID] = 0.0
596 continue
597
598 # Make a copy for probe molecule. It gets updated during the calculation...
599 ProbeMolConf = Chem.Mol(MolConfsMap[ConfID])
600 RMSD = rdMolAlign.AlignMol(ProbeMolConf, MinEnergyMolConf)
601 CalcRMSDMap[ConfID] = RMSD
602
603 # Track conformers with in the specified energy window from the lowest
604 # energy conformation along with applying RMSD cutoff as needed...
605 #
606 SelectedConfIDs = []
607
608 ConfCount = 0
609 IgnoredByEnergyConfCount = 0
610 IgnoredByRMSDConfCount = 0
611
612 FirstConf = True
613 for ConfID in SortedConfIDs:
614 if FirstConf:
615 ConfCount += 1
616 FirstConf = False
617 SelectedConfIDs.append(ConfID)
618 continue
619
620 ConfEnergyDiff = abs(CalcEnergyMap[ConfID] - MinConfEnergy)
621 if ConfEnergyDiff > EnergyWindow:
622 IgnoredByEnergyConfCount += 1
623 continue
624
625 if ApplyEnergyRMSDCutoff:
626 if CalcRMSDMap[ConfID] < EnergyRMSDCutoff:
627 IgnoredByRMSDConfCount += 1
628 continue
629
630 ConfCount += 1
631 SelectedConfIDs.append(ConfID)
632
633 if not OptionsInfo["QuietMode"]:
634 MolName = RDKitUtil.GetMolName(Mol, MolNum)
635 MiscUtil.PrintInfo("\nTotal Number of conformations generated for %s: %d" % (MolName, ConfCount))
636 MiscUtil.PrintInfo(
637 "Number of conformations ignored due to energy window cutoff: %d" % (IgnoredByEnergyConfCount)
638 )
639 if ApplyEnergyRMSDCutoff:
640 MiscUtil.PrintInfo(
641 "Number of conformations ignored due to energy RMSD cutoff: %d" % (IgnoredByRMSDConfCount)
642 )
643
644 # Setup selected conformer molecules...
645 SelectedConfMols = [MolConfsMap[ConfID] for ConfID in SelectedConfIDs]
646
647 # Setup selected conformer energy values...
648 SelectedConfEnergyValues = None
649 if OptionsInfo["EnergyOut"]:
650 SelectedConfEnergyValues = ["%.2f" % CalcEnergyMap[ConfID] for ConfID in SelectedConfIDs]
651
652 # Setup selected conformer scaffold RMSD values...
653 SelectedConfScaffoldEmbedRMSDValues = None
654 if OptionsInfo["ScaffoldRMSDOut"]:
655 SelectedConfScaffoldEmbedRMSDValues = []
656 for ConfID in SelectedConfIDs:
657 SelectedConfScaffoldEmbedRMSDValues.append(ScaffoldEmbedRMSDMap[ConfID])
658
659 return (SelectedConfMols, True, SelectedConfIDs, SelectedConfEnergyValues, SelectedConfScaffoldEmbedRMSDValues)
660
661
662 def GetEnergy(Mol, ConfID=None):
663 """Calculate energy."""
664
665 Status = True
666 Energy = None
667
668 if ConfID is None:
669 ConfID = -1
670
671 if OptionsInfo["UseUFF"]:
672 UFFMoleculeForcefield = AllChem.UFFGetMoleculeForceField(Mol, confId=ConfID)
673 if UFFMoleculeForcefield is None:
674 Status = False
675 else:
676 Energy = UFFMoleculeForcefield.CalcEnergy()
677 elif OptionsInfo["UseMMFF"]:
678 MMFFMoleculeProperties = AllChem.MMFFGetMoleculeProperties(Mol, mmffVariant=OptionsInfo["MMFFVariant"])
679 MMFFMoleculeForcefield = AllChem.MMFFGetMoleculeForceField(Mol, MMFFMoleculeProperties, confId=ConfID)
680 if MMFFMoleculeForcefield is None:
681 Status = False
682 else:
683 Energy = MMFFMoleculeForcefield.CalcEnergy()
684 else:
685 MiscUtil.PrintError(
686 "Couldn't retrieve conformer energy: Specified forcefield, %s, is not supported" % OptionsInfo["ForceField"]
687 )
688
689 return (Status, Energy)
690
691
692 def WriteMolConformers(Writer, Mol, MolNum, ConfMols, ConfIDs, ConfEnergyValues=None, ConfScaffoldEmbedRMSDValues=None):
693 """Write molecule coformers."""
694
695 if ConfMols is None:
696 return
697
698 for Index, ConfMol in enumerate(ConfMols):
699 ConfMolName = RDKitUtil.GetMolName(Mol, MolNum)
700 SetConfMolName(ConfMol, ConfMolName, ConfIDs[Index])
701
702 if ConfScaffoldEmbedRMSDValues is not None:
703 ConfMol.SetProp("CoreScaffoldEmbedRMSD", ConfScaffoldEmbedRMSDValues[Index])
704
705 if ConfEnergyValues is not None:
706 ConfMol.SetProp(OptionsInfo["EnergyLabel"], ConfEnergyValues[Index])
707
708 Writer.write(ConfMol)
709
710
711 def SetConfMolName(Mol, MolName, ConfCount):
712 """Set conf mol name."""
713
714 ConfName = "%s_Conf%d" % (MolName, ConfCount)
715 Mol.SetProp("_Name", ConfName)
716
717
718 def ProcessMCSParameters():
719 """Set up and process MCS parameters."""
720
721 SetupMCSParameters()
722 ProcessSpecifiedMCSParameters()
723
724
725 def SetupMCSParameters():
726 """Set up default MCS parameters."""
727
728 OptionsInfo["MCSParams"] = {
729 "MaximizeBonds": True,
730 "Threshold": 0.9,
731 "TimeOut": 3600,
732 "Verbose": False,
733 "MatchValences": True,
734 "MatchChiralTag": False,
735 "RingMatchesRingOnly": True,
736 "CompleteRingsOnly": True,
737 "AtomCompare": rdFMCS.AtomCompare.CompareElements,
738 "BondCompare": rdFMCS.BondCompare.CompareOrder,
739 "SeedSMARTS": "",
740 "MinNumAtoms": 1,
741 "MinNumBonds": 0,
742 }
743
744
745 def ProcessSpecifiedMCSParameters():
746 """Process specified MCS parameters."""
747
748 if re.match("^auto$", OptionsInfo["SpecifiedMCSParams"], re.I):
749 # Nothing to process...
750 return
751
752 # Parse specified parameters...
753 MCSParams = re.sub(" ", "", OptionsInfo["SpecifiedMCSParams"])
754 if not MCSParams:
755 MiscUtil.PrintError('No valid parameter name and value pairs specified using "-m, --mcsParams" option.')
756
757 MCSParamsWords = MCSParams.split(",")
758 if len(MCSParamsWords) % 2:
759 MiscUtil.PrintError(
760 'The number of comma delimited paramater names and values, %d, specified using "-m, --mcsParams" option must be an even number.'
761 % (len(MCSParamsWords))
762 )
763
764 # Setup canonical parameter names...
765 ValidParamNames = []
766 CanonicalParamNamesMap = {}
767 for ParamName in sorted(OptionsInfo["MCSParams"]):
768 ValidParamNames.append(ParamName)
769 CanonicalParamNamesMap[ParamName.lower()] = ParamName
770
771 # Validate and set paramater names and value...
772 for Index in range(0, len(MCSParamsWords), 2):
773 Name = MCSParamsWords[Index]
774 Value = MCSParamsWords[Index + 1]
775
776 CanonicalName = Name.lower()
777 if CanonicalName not in CanonicalParamNamesMap:
778 MiscUtil.PrintError(
779 'The parameter name, %s, specified using "-m, --mcsParams" option is not a valid name. Supported parameter names: %s'
780 % (Name, " ".join(ValidParamNames))
781 )
782
783 ParamName = CanonicalParamNamesMap[CanonicalName]
784 if re.match("^Threshold$", ParamName, re.I):
785 Value = float(Value)
786 if Value <= 0.0 or Value > 1.0:
787 MiscUtil.PrintError(
788 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: > 0 and <= 1.0'
789 % (Value, Name)
790 )
791 ParamValue = Value
792 elif re.match("^Timeout$", ParamName, re.I):
793 Value = int(Value)
794 if Value <= 0:
795 MiscUtil.PrintError(
796 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: > 0'
797 % (Value, Name)
798 )
799 ParamValue = Value
800 elif re.match("^MinNumAtoms$", ParamName, re.I):
801 Value = int(Value)
802 if Value < 1:
803 MiscUtil.PrintError(
804 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: >= 1'
805 % (Value, Name)
806 )
807 ParamValue = Value
808 elif re.match("^MinNumBonds$", ParamName, re.I):
809 Value = int(Value)
810 if Value < 0:
811 MiscUtil.PrintError(
812 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: >=0 '
813 % (Value, Name)
814 )
815 ParamValue = Value
816 elif re.match("^AtomCompare$", ParamName, re.I):
817 if re.match("^CompareAny$", Value, re.I):
818 ParamValue = rdFMCS.AtomCompare.CompareAny
819 elif re.match("^CompareElements$", Value, re.I):
820 ParamValue = Chem.rdFMCS.AtomCompare.CompareElements
821 elif re.match("^CompareIsotopes$", Value, re.I):
822 ParamValue = Chem.rdFMCS.AtomCompare.CompareIsotopes
823 else:
824 MiscUtil.PrintError(
825 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: CompareAny CompareElements CompareIsotopes'
826 % (Value, Name)
827 )
828 elif re.match("^BondCompare$", ParamName, re.I):
829 if re.match("^CompareAny$", Value, re.I):
830 ParamValue = Chem.rdFMCS.BondCompare.CompareAny
831 elif re.match("^CompareOrder$", Value, re.I):
832 ParamValue = rdFMCS.BondCompare.CompareOrder
833 elif re.match("^CompareOrderExact$", Value, re.I):
834 ParamValue = rdFMCS.BondCompare.CompareOrderExact
835 else:
836 MiscUtil.PrintError(
837 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: CompareAny CompareOrder CompareOrderExact'
838 % (Value, Name)
839 )
840 elif re.match("^SeedSMARTS$", ParamName, re.I):
841 if not len(Value):
842 MiscUtil.PrintError(
843 'The parameter value specified using "-m, --mcsParams" option for parameter, %s, is empty. '
844 % (Name)
845 )
846 ParamValue = Value
847 else:
848 if not re.match("^(Yes|No|True|False)$", Value, re.I):
849 MiscUtil.PrintError(
850 'The parameter value, %s, specified using "-m, --mcsParams" option for parameter, %s, is not a valid value. Supported values: Yes No True False'
851 % (Value, Name)
852 )
853 ParamValue = False
854 if re.match("^(Yes|True)$", Value, re.I):
855 ParamValue = True
856
857 # Set value...
858 OptionsInfo["MCSParams"][ParamName] = ParamValue
859
860
861 def ProcesssConformerGeneratorOption():
862 """Process comformer generator option."""
863
864 ConfGenParams = MiscUtil.ProcessOptionConformerGenerator("--conformerGenerator", Options["--conformerGenerator"])
865
866 OptionsInfo["ConformerGenerator"] = ConfGenParams["ConformerGenerator"]
867 OptionsInfo["UseBasicKnowledge"] = ConfGenParams["UseBasicKnowledge"]
868 OptionsInfo["UseExpTorsionAnglePrefs"] = ConfGenParams["UseExpTorsionAnglePrefs"]
869 OptionsInfo["ETVersion"] = ConfGenParams["ETVersion"]
870
871
872 def ProcessOptions():
873 """Process and validate command line arguments and options."""
874
875 MiscUtil.PrintInfo("Processing options...")
876
877 # Validate options...
878 ValidateOptions()
879
880 OptionsInfo["Infile"] = Options["--infile"]
881 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
882 "--infileParams", Options["--infileParams"], Options["--infile"]
883 )
884
885 OptionsInfo["RefFile"] = Options["--reffile"]
886
887 OptionsInfo["Scaffold"] = Options["--scaffold"]
888 if re.match("^auto$", Options["--scaffold"], re.I):
889 UseScaffoldMCS = True
890 UseScaffoldSMARTS = False
891 ScaffoldSMARTS = None
892 else:
893 UseScaffoldMCS = False
894 UseScaffoldSMARTS = True
895 ScaffoldSMARTS = OptionsInfo["Scaffold"]
896
897 OptionsInfo["UseScaffoldMCS"] = UseScaffoldMCS
898 OptionsInfo["UseScaffoldSMARTS"] = UseScaffoldSMARTS
899 OptionsInfo["ScaffoldSMARTS"] = ScaffoldSMARTS
900 OptionsInfo["ScaffoldPatternMol"] = None
901
902 OptionsInfo["SpecifiedMCSParams"] = Options["--mcsParams"]
903 ProcessMCSParameters()
904
905 OptionsInfo["Outfile"] = Options["--outfile"]
906 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
907 "--outfileParams", Options["--outfileParams"]
908 )
909
910 OptionsInfo["Overwrite"] = Options["--overwrite"]
911
912 OptionsInfo["AddHydrogens"] = True if re.match("^yes$", Options["--addHydrogens"], re.I) else False
913
914 ProcesssConformerGeneratorOption()
915
916 if re.match("^UFF$", Options["--forceField"], re.I):
917 ForceField = "UFF"
918 UseUFF = True
919 UseMMFF = False
920 elif re.match("^MMFF$", Options["--forceField"], re.I):
921 ForceField = "MMFF"
922 UseUFF = False
923 UseMMFF = True
924 else:
925 MiscUtil.PrintError(
926 'The value, %s, specified for "--forceField" is not supported.' % (Options["--forceField"],)
927 )
928
929 MMFFVariant = "MMFF94" if re.match("^MMFF94$", Options["--forceFieldMMFFVariant"], re.I) else "MMFF94s"
930
931 OptionsInfo["ForceField"] = ForceField
932 OptionsInfo["MMFFVariant"] = MMFFVariant
933 OptionsInfo["UseMMFF"] = UseMMFF
934 OptionsInfo["UseUFF"] = UseUFF
935
936 OptionsInfo["ScaffoldRMSDOut"] = True if re.match("^yes$", Options["--scaffoldRMSDOut"], re.I) else False
937
938 OptionsInfo["EnergyOut"] = True if re.match("^yes$", Options["--energyOut"], re.I) else False
939 if UseMMFF:
940 OptionsInfo["EnergyLabel"] = "%s_Energy" % MMFFVariant
941 else:
942 OptionsInfo["EnergyLabel"] = "%s_Energy" % ForceField
943
944 OptionsInfo["EnforceChirality"] = True if re.match("^yes$", Options["--enforceChirality"], re.I) else False
945 EnergyRMSDCutoff = -1.0
946 if not re.match("^none$", Options["--energyRMSDCutoff"], re.I):
947 EnergyRMSDCutoff = float(Options["--energyRMSDCutoff"])
948 OptionsInfo["EnergyRMSDCutoff"] = EnergyRMSDCutoff
949
950 OptionsInfo["EnergyWindow"] = float(Options["--energyWindow"])
951
952 OptionsInfo["MaxConfs"] = int(Options["--maxConfs"])
953
954 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
955 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
956
957 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
958
959 OptionsInfo["RemoveHydrogens"] = True if re.match("^yes$", Options["--removeHydrogens"], re.I) else False
960 OptionsInfo["UseTethers"] = True if re.match("^yes$", Options["--useTethers"], re.I) else False
961
962
963 def RetrieveOptions():
964 """Retrieve command line arguments and options."""
965
966 # Get options...
967 global Options
968 Options = docopt(_docoptUsage_)
969
970 # Set current working directory to the specified directory...
971 WorkingDir = Options["--workingdir"]
972 if WorkingDir:
973 os.chdir(WorkingDir)
974
975 # Handle examples option...
976 if "--examples" in Options and Options["--examples"]:
977 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
978 sys.exit(0)
979
980
981 def ValidateOptions():
982 """Validate option values."""
983
984 MiscUtil.ValidateOptionTextValue("-a, --addHydrogens", Options["--addHydrogens"], "yes no")
985 MiscUtil.ValidateOptionTextValue(
986 "-c, --conformerGenerator", Options["--conformerGenerator"], "SDG KDG ETDG ETKDG ETKDGv2"
987 )
988
989 MiscUtil.ValidateOptionTextValue("-f, --forceField", Options["--forceField"], "UFF MMFF")
990 MiscUtil.ValidateOptionTextValue(" --forceFieldMMFFVariant", Options["--forceFieldMMFFVariant"], "MMFF94 MMFF94s")
991
992 MiscUtil.ValidateOptionFloatValue("--energyWindow", Options["--energyWindow"], {">": 0.0})
993 if not re.match("^none$", Options["--energyRMSDCutoff"], re.I):
994 MiscUtil.ValidateOptionFloatValue("--energyRMSDCutoff", Options["--energyRMSDCutoff"], {">": 0.0})
995
996 MiscUtil.ValidateOptionTextValue("--scaffoldRMSDOut", Options["--scaffoldRMSDOut"], "yes no")
997
998 MiscUtil.ValidateOptionTextValue("--energyOut", Options["--energyOut"], "yes no")
999 MiscUtil.ValidateOptionTextValue("--enforceChirality ", Options["--enforceChirality"], "yes no")
1000
1001 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1002 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
1003
1004 MiscUtil.ValidateOptionFilePath("-r, --reffile", Options["--reffile"])
1005 MiscUtil.ValidateOptionFileExt("-r, --reffile", Options["--reffile"], "sdf sd mol")
1006
1007 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
1008 MiscUtil.ValidateOptionsOutputFileOverwrite(
1009 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
1010 )
1011 MiscUtil.ValidateOptionsDistinctFileNames(
1012 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
1013 )
1014
1015 MiscUtil.ValidateOptionIntegerValue("--maxConfs", Options["--maxConfs"], {">": 0})
1016
1017 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1018 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1019
1020 MiscUtil.ValidateOptionTextValue("-r, --removeHydrogens", Options["--removeHydrogens"], "yes no")
1021
1022 MiscUtil.ValidateOptionTextValue("-u, --useTethers", Options["--useTethers"], "yes no")
1023
1024
1025 # Setup a usage string for docopt...
1026 _docoptUsage_ = """
1027 RDKitGenerateConstrainedConformers.py - Generate constrained molecular conformations
1028
1029 Usage:
1030 RDKitGenerateConstrainedConformers.py [--addHydrogens <yes or no>] [--conformerGenerator <SDG, KDG, ETDG, ETKDG, ETKDGv2>]
1031 [--forceField <UFF, or MMFF>] [--forceFieldMMFFVariant <MMFF94 or MMFF94s>]
1032 [--energyOut <yes or no>] [--enforceChirality <yes or no>] [--energyRMSDCutoff <number>]
1033 [--energyWindow <number> ] [--infileParams <Name,Value,...>] [--maxConfs <number>]
1034 [--mcsParams <Name,Value,...>] [--mp <yes or no>] [--mpParams <Name,Value,...>]
1035 [ --outfileParams <Name,Value,...> ] [--overwrite] [--quiet <yes or no>] [ --removeHydrogens <yes or no>]
1036 [--scaffold <auto or SMARTS>] [--scaffoldRMSDOut <yes or no>] [--useTethers <yes or no>]
1037 [-w <dir>] -i <infile> -r <reffile> -o <outfile>
1038 RDKitGenerateConstrainedConformers.py -h | --help | -e | --examples
1039
1040 Description:
1041 Generate molecular conformations by performing a constrained energy minimization
1042 against a reference molecule. An initial set of 3D conformers are generated for the input
1043 molecules using distance geometry. A common core scaffold, corresponding to
1044 a Maximum Common Substructure (MCS) or an explicit SMARTS pattern, is identified
1045 between a pair of input and reference molecules. The core scaffold atoms in input
1046 molecules are aligned against the same atoms in the reference molecule. The energy
1047 of aligned structures are minimized using the forcefield to generate the final 3D structures.
1048
1049 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
1050 .csv, .tsv, .txt)
1051
1052 The supported output file formats are: SD (.sdf, .sd)
1053
1054 Options:
1055 -a, --addHydrogens <yes or no> [default: yes]
1056 Add hydrogens before minimization.
1057 -c, --conformerGenerator <text> [default: ETKDGv2]
1058 Conformation generation methodology for generating initial 3D coordinates
1059 for molecules in input file. A common core scaffold is identified between a
1060 a pair of input and reference molecules. The atoms in common core scaffold
1061 of input molecules are aligned against the reference molecule followed by
1062 energy minimization to generate final 3D structure.
1063
1064 The possible values along with a brief description are shown below:
1065
1066 SDG: Standard Distance Geometry
1067 KDG: basic Knowledge-terms with Distance Geometry
1068 ETDG: Experimental Torsion-angle preference with Distance Geometry
1069 ETKDG: Experimental Torsion-angle preference along with basic
1070 Knowledge-terms and Distance Geometry [Ref 129]
1071 ETKDGv2: Experimental Torsion-angle preference along with basic
1072 Knowledge-terms and Distance Geometry [Ref 167]
1073 -f, --forceField <UFF, MMFF> [default: MMFF]
1074 Forcefield method to use for constrained energy minimization. Possible values:
1075 Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics Force
1076 Field [ Ref 83-87 ] .
1077 --forceFieldMMFFVariant <MMFF94 or MMFF94s> [default: MMFF94]
1078 Variant of MMFF forcefield to use for energy minimization.
1079 --energyOut <yes or no> [default: No]
1080 Write out energy values.
1081 --enforceChirality <yes or no> [default: Yes]
1082 Enforce chirality for defined chiral centers.
1083 --energyRMSDCutoff <number> [default: 0.5]
1084 RMSD cutoff for retaining conformations after embedding and energy minimization.
1085 Possible values: A number or None
1086
1087 The default is to keep only those conformations which are different from the
1088 lowest energy conformation by the specified RMSD cutoff. The None value may
1089 be used to keep all minimized conformations with in the specified energy window
1090 from the lowest energy conformation. The lowest energy conformation is always
1091 retained.
1092 --energyWindow <number> [default: 20]
1093 Energy window in kcal/mol for selecting conformers.
1094 -e, --examples
1095 Print examples.
1096 -h, --help
1097 Print this help message.
1098 -i, --infile <infile>
1099 Input file name.
1100 --infileParams <Name,Value,...> [default: auto]
1101 A comma delimited list of parameter name and value pairs for reading
1102 molecules from files. The supported parameter names for different file
1103 formats, along with their default values, are shown below:
1104
1105 SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
1106
1107 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
1108 smilesTitleLine,auto,sanitize,yes
1109
1110 Possible values for smilesDelimiter: space, comma or tab.
1111 --maxConfs <number> [default: 50]
1112 Maximum number of conformations to generate for each molecule by conformation
1113 generation methodology for initial 3D coordinates. A constrained minimization is
1114 performed using the specified forcefield and the lowest energy conformation is written
1115 to the output file.
1116 --mcsParams <Name,Value,...> [default: auto]
1117 Parameter values to use for identifying a maximum common substructure
1118 (MCS) in between a pair of reference and input molecules.In general, it is a
1119 comma delimited list of parameter name and value pairs. The supported
1120 parameter names along with their default values are shown below:
1121
1122 atomCompare,CompareElements,bondCompare,CompareOrder,
1123 maximizeBonds,yes,matchValences,yes,matchChiralTag,no,
1124 minNumAtoms,1,minNumBonds,0,ringMatchesRingOnly,yes,
1125 completeRingsOnly,yes,threshold,1.0,timeOut,3600,seedSMARTS,none
1126
1127 Possible values for atomCompare: CompareAny, CompareElements,
1128 CompareIsotopes. Possible values for bondCompare: CompareAny,
1129 CompareOrder, CompareOrderExact.
1130
1131 A brief description of MCS parameters taken from RDKit documentation is
1132 as follows:
1133
1134 atomCompare - Controls match between two atoms
1135 bondCompare - Controls match between two bonds
1136 maximizeBonds - Maximize number of bonds instead of atoms
1137 matchValences - Include atom valences in the MCS match
1138 matchChiralTag - Include atom chirality in the MCS match
1139 minNumAtoms - Minimum number of atoms in the MCS match
1140 minNumBonds - Minimum number of bonds in the MCS match
1141 ringMatchesRingOnly - Ring bonds only match other ring bonds
1142 completeRingsOnly - Partial rings not allowed during the match
1143 threshold - Fraction of the dataset that must contain the MCS
1144 seedSMARTS - SMARTS string as the seed of the MCS
1145 timeout - Timeout for the MCS calculation in seconds
1146
1147 --mp <yes or no> [default: no]
1148 Use multiprocessing.
1149
1150 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1151 function employing lazy RDKit data iterable. This allows processing of
1152 arbitrary large data sets without any additional requirements memory.
1153
1154 All input data may be optionally loaded into memory by mp.Pool.map()
1155 before starting worker processes in a process pool by setting the value
1156 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1157
1158 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1159 data mode may adversely impact the performance. The '--mpParams' section
1160 provides additional information to tune the value of 'chunkSize'.
1161 --mpParams <Name,Value,...> [default: auto]
1162 A comma delimited list of parameter name and value pairs to configure
1163 multiprocessing.
1164
1165 The supported parameter names along with their default and possible
1166 values are shown below:
1167
1168 chunkSize, auto
1169 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1170 numProcesses, auto [ Default: mp.cpu_count() ]
1171
1172 These parameters are used by the following functions to configure and
1173 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1174 mp.Pool.imap().
1175
1176 The chunkSize determines chunks of input data passed to each worker
1177 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1178 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1179
1180 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1181 automatically converts RDKit data iterable into a list, loads all data into
1182 memory, and calculates the default chunkSize using the following method
1183 as shown in its code:
1184
1185 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1186 if extra: chunkSize += 1
1187
1188 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1189 and 100 data items.
1190
1191 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1192 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1193 data into memory. Consequently, the size of input data is not known a priori.
1194 It's not possible to estimate an optimal value for the chunkSize. The default
1195 chunkSize is set to 1.
1196
1197 The default value for the chunkSize during 'Lazy' data mode may adversely
1198 impact the performance due to the overhead associated with exchanging
1199 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1200 a larger value during 'Lazy' input data mode, based on the size of your input
1201 data and number of processes in the process pool.
1202
1203 The mp.Pool.map() function waits for all worker processes to process all
1204 the data and return the results. The mp.Pool.imap() function, however,
1205 returns the the results obtained from worker processes as soon as the
1206 results become available for specified chunks of data.
1207
1208 The order of data in the results returned by both mp.Pool.map() and
1209 mp.Pool.imap() functions always corresponds to the input data.
1210 -o, --outfile <outfile>
1211 Output file name.
1212 --outfileParams <Name,Value,...> [default: auto]
1213 A comma delimited list of parameter name and value pairs for writing
1214 molecules to files. The supported parameter names for different file
1215 formats, along with their default values, are shown below:
1216
1217 SD: kekulize,yes,forceV3000,no
1218
1219 --overwrite
1220 Overwrite existing files.
1221 -q, --quiet <yes or no> [default: no]
1222 Use quiet mode. The warning and information messages will not be printed.
1223 -r, --reffile <reffile>
1224 Reference input file name containing a 3D reference molecule. A common
1225 core scaffold must be present in a pair of an input and reference molecules.
1226 Otherwise, no constrained minimization is performed on the input molecule.
1227 --removeHydrogens <yes or no> [default: Yes]
1228 Remove hydrogens after minimization.
1229 -s, --scaffold <auto or SMARTS> [default: auto]
1230 Common core scaffold between a pair of input and reference molecules used for
1231 constrained minimization of molecules in input file. Possible values: Auto or a
1232 valid SMARTS pattern. The common core scaffold is automatically detected
1233 corresponding to the Maximum Common Substructure (MCS) between a pair of
1234 reference and input molecules. A valid SMARTS pattern may be optionally specified
1235 for the common core scaffold.
1236 --scaffoldRMSDOut <yes or no> [default: No]
1237 Write out RMSD value for common core alignment between a pair of input and
1238 reference molecules.
1239 -u, --useTethers <yes or no> [default: yes]
1240 Use tethers to optimize the final conformation by applying a series of extra forces
1241 to align matching atoms to the positions of the core atoms. Otherwise, use simple
1242 distance constraints during the optimization.
1243 -w, --workingdir <dir>
1244 Location of working directory which defaults to the current directory.
1245
1246 Examples:
1247 To generate conformers by performing constrained energy minimization for molecules
1248 in a SMILES file against a reference 3D molecule in a SD file using a common core
1249 scaffold between pairs of input and reference molecules identified using MCS,
1250 generating up to 50 conformations using ETKDG methodology followed by MMFF
1251 forcefield minimization within energy window of 20 kcal/mol and RMSD of greater
1252 than 0.5 from the lowest energy conformation, and write out a SD file:
1253
1254 % RDKitGenerateConstrainedConformers.py -i SampleSeriesD3R.smi
1255 -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1256
1257 To rerun the first example in a quiet mode and write out a SD file, type:
1258
1259 % RDKitGenerateConstrainedConformers.py -q yes -i SampleSeriesD3R.smi
1260 -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1261
1262 To rerun the first example in multiprocessing mode on all available CPUs
1263 without loading all data into memory and write out a SD file, type:
1264
1265 % RDKitGenerateConstrainedConformers.py --mp yes -i SampleSeriesD3R.smi
1266 -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1267
1268 To run the first example in multiprocessing mode on all available CPUs
1269 by loading all data into memory and write out a SD file, type:
1270
1271 % RDKitGenerateConstrainedConformers.py --mp yes --mpParams
1272 "inputDataMode,InMemory" -i SampleSeriesD3R.smi
1273 -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1274
1275 To rerun the first example in multiprocessing mode on specific number of
1276 CPUs and chunk size without loading all data into memory and write out a SD file,
1277 type:
1278
1279 % RDKitGenerateConstrainedConformers.py --mp yes --mpParams
1280 "inputDataMode,Lazy,numProcesses,4,chunkSize,8"
1281 -i SampleSeriesD3R.smi -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1282
1283 To rerun the first example using an explicit SMARTS string for a common core
1284 scaffold and write out a SD file, type:
1285
1286 % RDKitGenerateConstrainedConformers.py --scaffold
1287 "c2cc(-c3nc(N)ncc3)cn2" -i SampleSeriesD3R.smi
1288 -r SampleSeriesRef3D.sdf -o SampleOut.sdf
1289
1290 To rerun the first example using molecules in a CSV SMILES file, SMILES
1291 strings in column 1, name in column2, and write out a SD file, type:
1292
1293 % RDKitGenerateConstrainedConformers.py --infileParams
1294 "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1,
1295 smilesNameColumn,2" -i SampleSeriesD3R.csv -r SampleSeriesRef3D.sdf
1296 -o SampleOut.sdf
1297
1298 To generate constrained conformers for molecules in a SD file against a reference
1299 3D molecule in a SD file using a common core scaffold between pairs of input and
1300 reference molecules identified using MCS, generating up to 10 conformations
1301 using SDG methodology followed by UFF forcefield minimization, conformations
1302 with in an energy window of 10 kcal/mol and RMSD of greater that 1, and write out
1303 a SD file containing minimum energy structure along with energy and embed RMS
1304 values corresponding to each constrained molecule, type:
1305
1306 % RDKitGenerateConstrainedConformers.py --maxConfs 10 -c SDG -f UFF
1307 --scaffoldRMSDOut yes --energyOut yes --energyRMSDCutoff 1.0
1308 --energyWindow 10 -i SampleSeriesD3R.sdf -r SampleSeriesRef3D.sdf
1309 -o SampleOut.sdf
1310
1311 Author:
1312 Manish Sud(msud@san.rr.com)
1313
1314 See also:
1315 RDKitCalculateRMSD.py, RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py,
1316 RDKitConvertFileFormat.py, RDKitGenerateConformers.py, RDKitPerformConstrainedMinimization.py
1317
1318 Copyright:
1319 Copyright (C) 2026 Manish Sud. All rights reserved.
1320
1321 The functionality available in this script is implemented using RDKit, an
1322 open source toolkit for cheminformatics developed by Greg Landrum.
1323
1324 This file is part of MayaChemTools.
1325
1326 MayaChemTools is free software; you can redistribute it and/or modify it under
1327 the terms of the GNU Lesser General Public License as published by the Free
1328 Software Foundation; either version 3 of the License, or (at your option) any
1329 later version.
1330
1331 """
1332
1333 if __name__ == "__main__":
1334 main()
