1 #!/bin/env python
2 #
3 # File: Psi4GenerateConformers.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 Psi4, an
9 # open source quantum chemistry software package, and RDKit, an open
10 # source toolkit for cheminformatics developed by Greg Landrum.
11 #
12 # This file is part of MayaChemTools.
13 #
14 # MayaChemTools is free software; you can redistribute it and/or modify it under
15 # the terms of the GNU Lesser General Public License as published by the Free
16 # Software Foundation; either version 3 of the License, or (at your option) any
17 # later version.
18 #
19 # MayaChemTools is distributed in the hope that it will be useful, but without
20 # any warranty; without even the implied warranty of merchantability of fitness
21 # for a particular purpose. See the GNU Lesser General Public License for more
22 # details.
23 #
24 # You should have received a copy of the GNU Lesser General Public License
25 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
26 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
27 # Boston, MA, 02111-1307, USA.
28 #
29
30 from __future__ import print_function
31
32 import os
33 import sys
34 import time
35 import re
36 import shutil
37 import multiprocessing as mp
38
39 # Psi4 imports...
40 if hasattr(shutil, "which") and shutil.which("psi4") is None:
41 sys.stderr.write("\nWarning: Failed to find 'psi4' in your PATH indicating potential issues with your\n")
42 sys.stderr.write("Psi4 environment. The 'import psi4' directive in the global scope of the script\n")
43 sys.stderr.write("interferes with the multiprocessing functionality. It is imported later in the\n")
44 sys.stderr.write("local scope during the execution of the script and may fail. Check/update your\n")
45 sys.stderr.write("Psi4 environment and try again.\n\n")
46
47 # RDKit imports...
48 try:
49 from rdkit import rdBase
50 from rdkit import Chem
51 from rdkit.Chem import AllChem
52 except ImportError as ErrMsg:
53 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
54 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
55 sys.exit(1)
56
57 # MayaChemTools imports...
58 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
59 try:
60 from docopt import docopt
61 import MiscUtil
62 import Psi4Util
63 import RDKitUtil
64 except ImportError as ErrMsg:
65 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
66 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
67 sys.exit(1)
68
69 ScriptName = os.path.basename(sys.argv[0])
70 Options = {}
71 OptionsInfo = {}
72
73
74 def main():
75 """Start execution of the script."""
76
77 MiscUtil.PrintInfo(
78 "\n%s (Psi4: Imported later; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
79 % (ScriptName, rdBase.rdkitVersion, MiscUtil.GetMayaChemToolsVersion(), time.asctime())
80 )
81
82 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
83
84 # Retrieve command line arguments and options...
85 RetrieveOptions()
86
87 # Process and validate command line arguments and options...
88 ProcessOptions()
89
90 # Perform actions required by the script...
91 GenerateConformers()
92
93 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
94 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
95
96
97 def GenerateConformers():
98 """Generate conformers."""
99
100 # Setup a molecule reader...
101 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
102 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
103
104 # Set up a molecule writer...
105 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
106 if Writer is None:
107 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
108 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
109
110 MolCount, ValidMolCount, ConfGenFailedCount = ProcessMolecules(Mols, Writer)
111
112 if Writer is not None:
113 Writer.close()
114
115 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
116 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
117 MiscUtil.PrintInfo(
118 "Number of molecules failed during conformation generation or minimization: %d" % ConfGenFailedCount
119 )
120 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + ConfGenFailedCount))
121
122
123 def ProcessMolecules(Mols, Writer):
124 """Process molecules to generate conformers."""
125
126 if OptionsInfo["MPMode"]:
127 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer)
128 else:
129 return ProcessMoleculesUsingSingleProcess(Mols, Writer)
130
131
132 def ProcessMoleculesUsingSingleProcess(Mols, Writer):
133 """Process molecules to generate conformers using a single process."""
134
135 # Intialize Psi4...
136 MiscUtil.PrintInfo("\nInitializing Psi4...")
137 Psi4Handle = Psi4Util.InitializePsi4(
138 Psi4RunParams=OptionsInfo["Psi4RunParams"],
139 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
140 PrintVersion=True,
141 PrintHeader=True,
142 )
143 OptionsInfo["psi4"] = Psi4Handle
144
145 # Setup max iterations global variable...
146 Psi4Util.UpdatePsi4OptionsParameters(Psi4Handle, {"GEOM_MAXITER": OptionsInfo["MaxIters"]})
147
148 # Setup conversion factor for energy units...
149 SetupEnergyConversionFactor(Psi4Handle)
150
151 MiscUtil.PrintInfo("\nGenerating conformers and performing energy minimization...")
152
153 (MolCount, ValidMolCount, ConfGenFailedCount) = [0] * 3
154 for Mol in Mols:
155 MolCount += 1
156
157 if not CheckAndValidateMolecule(Mol, MolCount):
158 continue
159
160 # Setup 2D coordinates for SMILES input file...
161 if OptionsInfo["SMILESInfileStatus"]:
162 AllChem.Compute2DCoords(Mol)
163
164 ValidMolCount += 1
165
166 ConformerMol, CalcStatus, ConfIDs, ConfEnergies = GenerateMolConformers(Psi4Handle, Mol, MolCount)
167
168 if not CalcStatus:
169 if not OptionsInfo["QuietMode"]:
170 MiscUtil.PrintWarning(
171 "Failed to calculate energy for molecule %s" % RDKitUtil.GetMolName(Mol, MolCount)
172 )
173
174 ConfGenFailedCount += 1
175 continue
176
177 WriteMolConformers(Writer, ConformerMol, MolCount, ConfIDs, ConfEnergies)
178
179 return (MolCount, ValidMolCount, ConfGenFailedCount)
180
181
182 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer):
183 """Process and minimize molecules using multiprocessing."""
184
185 if OptionsInfo["MPLevelConformersMode"]:
186 return ProcessMoleculesUsingMultipleProcessesAtConformersLevel(Mols, Writer)
187 elif OptionsInfo["MPLevelMoleculesMode"]:
188 return ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols, Writer)
189 else:
190 MiscUtil.PrintError('The value, %s, option "--mpLevel" is not supported.' % (OptionsInfo["MPLevel"]))
191
192
193 def ProcessMoleculesUsingMultipleProcessesAtMoleculesLevel(Mols, Writer):
194 """Process molecules to generate conformers using multiprocessing at molecules level."""
195
196 MiscUtil.PrintInfo(
197 "\nGenerating conformers and performing energy minimization using multiprocessing at molecules level..."
198 )
199
200 MPParams = OptionsInfo["MPParams"]
201
202 # Setup data for initializing a worker process...
203 InitializeWorkerProcessArgs = (
204 MiscUtil.ObjectToBase64EncodedString(Options),
205 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
206 )
207
208 # Setup a encoded mols data iterable for a worker process...
209 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
210
211 # Setup process pool along with data initialization for each process...
212 if not OptionsInfo["QuietMode"]:
213 MiscUtil.PrintInfo(
214 "\nConfiguring multiprocessing using %s method..."
215 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
216 )
217 MiscUtil.PrintInfo(
218 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
219 % (
220 MPParams["NumProcesses"],
221 MPParams["InputDataMode"],
222 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
223 )
224 )
225
226 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
227
228 # Start processing...
229 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
230 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
231 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
232 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
233 else:
234 MiscUtil.PrintError(
235 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
236 )
237
238 # Print out Psi4 version in the main process...
239 MiscUtil.PrintInfo("\nInitializing Psi4...\n")
240 Psi4Handle = Psi4Util.InitializePsi4(PrintVersion=True, PrintHeader=False)
241 OptionsInfo["psi4"] = Psi4Handle
242
243 (MolCount, ValidMolCount, ConfGenFailedCount) = [0] * 3
244 for Result in Results:
245 MolCount += 1
246 MolIndex, EncodedMol, CalcStatus, ConfIDs, ConfEnergies = Result
247
248 if EncodedMol is None:
249 continue
250 ValidMolCount += 1
251
252 if not CalcStatus:
253 ConfGenFailedCount += 1
254 continue
255
256 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
257 WriteMolConformers(Writer, Mol, MolCount, ConfIDs, ConfEnergies)
258
259 return (MolCount, ValidMolCount, ConfGenFailedCount)
260
261
262 def InitializeWorkerProcess(*EncodedArgs):
263 """Initialize data for a worker process."""
264
265 global Options, OptionsInfo
266
267 if not OptionsInfo["QuietMode"]:
268 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
269
270 # Decode Options and OptionInfo...
271 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
272 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
273
274 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
275 # output files for each process...
276 OptionsInfo["Psi4Initialized"] = False
277
278
279 def WorkerProcess(EncodedMolInfo):
280 """Process data for a worker process."""
281
282 if not OptionsInfo["Psi4Initialized"]:
283 InitializePsi4ForWorkerProcess()
284
285 MolIndex, EncodedMol = EncodedMolInfo
286 MolNum = MolIndex + 1
287
288 CalcStatus = False
289 ConfIDs = None
290 ConfEnergies = None
291
292 if EncodedMol is None:
293 return [MolIndex, None, CalcStatus, ConfIDs, ConfEnergies]
294
295 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
296
297 if not CheckAndValidateMolecule(Mol, MolNum):
298 return [MolIndex, None, CalcStatus, ConfIDs, ConfEnergies]
299
300 # Setup 2D coordinates for SMILES input file...
301 if OptionsInfo["SMILESInfileStatus"]:
302 AllChem.Compute2DCoords(Mol)
303
304 Mol, CalcStatus, ConfIDs, ConfEnergies = GenerateMolConformers(OptionsInfo["psi4"], Mol, MolNum)
305
306 return [
307 MolIndex,
308 RDKitUtil.MolToBase64EncodedMolString(
309 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
310 ),
311 CalcStatus,
312 ConfIDs,
313 ConfEnergies,
314 ]
315
316
317 def ProcessMoleculesUsingMultipleProcessesAtConformersLevel(Mols, Writer):
318 """Process molecules to generate conformers using multiprocessing at conformers level."""
319
320 MiscUtil.PrintInfo(
321 "\nPerforming minimization with generation of conformers using multiprocessing at conformers level..."
322 )
323
324 (MolCount, ValidMolCount, ConfGenFailedCount) = [0] * 3
325 for Mol in Mols:
326 MolCount += 1
327
328 if not CheckAndValidateMolecule(Mol, MolCount):
329 continue
330
331 # Setup 2D coordinates for SMILES input file...
332 if OptionsInfo["SMILESInfileStatus"]:
333 AllChem.Compute2DCoords(Mol)
334
335 ValidMolCount += 1
336
337 Mol, CalcStatus, ConfIDs, ConfEnergies = ProcessConformersUsingMultipleProcesses(Mol, MolCount)
338
339 if not CalcStatus:
340 ConfGenFailedCount += 1
341 continue
342
343 WriteMolConformers(Writer, Mol, MolCount, ConfIDs, ConfEnergies)
344
345 return (MolCount, ValidMolCount, ConfGenFailedCount)
346
347
348 def ProcessConformersUsingMultipleProcesses(Mol, MolNum):
349 """Generate coformers and minimize them using multiple processes."""
350
351 # Add hydrogens...
352 Mol = Chem.AddHs(Mol)
353
354 # Setup conformers...
355 ConfIDs = EmbedMolecule(Mol, MolNum)
356 if not len(ConfIDs):
357 if not OptionsInfo["QuietMode"]:
358 MolName = RDKitUtil.GetMolName(Mol, MolNum)
359 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s: Embedding failed...\n" % MolName)
360 return (Mol, False, None, None)
361
362 MPParams = OptionsInfo["MPParams"]
363
364 # Setup data for initializing a worker process...
365 InitializeWorkerProcessArgs = (
366 MiscUtil.ObjectToBase64EncodedString(Options),
367 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
368 )
369
370 # Setup a encoded mols data iterable for a worker process...
371 MolIndex = MolNum - 1
372 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStringWithConfIDs(Mol, MolIndex, ConfIDs)
373
374 # Setup process pool along with data initialization for each process...
375 if not OptionsInfo["QuietMode"]:
376 MiscUtil.PrintInfo(
377 "\nConfiguring multiprocessing using %s method..."
378 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
379 )
380 MiscUtil.PrintInfo(
381 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
382 % (
383 MPParams["NumProcesses"],
384 MPParams["InputDataMode"],
385 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
386 )
387 )
388
389 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeConformerWorkerProcess, InitializeWorkerProcessArgs)
390
391 # Start processing...
392 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
393 Results = ProcessPool.imap(ConformerWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
394 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
395 Results = ProcessPool.map(ConformerWorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
396 else:
397 MiscUtil.PrintError(
398 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
399 )
400
401 CalcEnergyMap = {}
402 CalcFailedCount = 0
403 for Result in Results:
404 MolIndex, EncodedMol, CalcStatus, ConfID, Energy = Result
405
406 if EncodedMol is None:
407 CalcFailedCount += 1
408 continue
409
410 if not CalcStatus:
411 CalcFailedCount += 1
412 continue
413
414 # Retrieve minimized atom positions...
415 MinimizedMol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
416 AtomPositions = RDKitUtil.GetAtomPositions(MinimizedMol, ConfID=ConfID)
417
418 # Update atom positions...
419 RDKitUtil.SetAtomPositions(Mol, AtomPositions, ConfID=ConfID)
420
421 CalcEnergyMap[ConfID] = Energy
422
423 if CalcFailedCount:
424 return (Mol, False, None, None)
425
426 # Align molecules after minimization...
427 if OptionsInfo["ConfGenerationParams"]["AlignConformers"]:
428 AllChem.AlignMolConformers(Mol)
429
430 # Filter conformers...
431 SelectedConfIDs, SelectedConfEnergies = FilterMolConformers(Mol, MolNum, ConfIDs, CalcEnergyMap)
432
433 return [Mol, True, SelectedConfIDs, SelectedConfEnergies]
434
435
436 def InitializeConformerWorkerProcess(*EncodedArgs):
437 """Initialize data for a conformer worker process."""
438
439 global Options, OptionsInfo
440
441 if not OptionsInfo["QuietMode"]:
442 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
443
444 # Decode Options and OptionInfo...
445 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
446 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
447
448 # Psi4 is initialized in the worker process to avoid creation of redundant Psi4
449 # output files for each process...
450 OptionsInfo["Psi4Initialized"] = False
451
452
453 def ConformerWorkerProcess(EncodedMolInfo):
454 """Process data for a conformer worker process."""
455
456 if not OptionsInfo["Psi4Initialized"]:
457 InitializePsi4ForWorkerProcess()
458
459 MolIndex, EncodedMol, ConfID = EncodedMolInfo
460
461 MolNum = MolIndex + 1
462
463 CalcStatus = False
464 Energy = None
465
466 if EncodedMol is None:
467 return [MolIndex, None, CalcStatus, ConfID, Energy]
468
469 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
470 MolName = RDKitUtil.GetMolName(Mol, MolNum)
471
472 if not OptionsInfo["QuietMode"]:
473 MiscUtil.PrintInfo("Processing conformer ID %s for molecule %s..." % (ConfID, MolName))
474
475 Status, ConvergeStatus = MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID)
476 if not Status:
477 return [MolIndex, EncodedMol, CalcStatus, ConfID, Energy]
478
479 if ConvergeStatus != 0:
480 if not OptionsInfo["QuietMode"]:
481 MolName = RDKitUtil.GetMolName(Mol, (MolIndex + 1))
482 MiscUtil.PrintWarning(
483 'Minimization using forcefield failed to converge for molecule %s in %d steps. Try using higher value for "maxIters" in "--confParams" option...\n'
484 % (MolName, OptionsInfo["ConfGenerationParams"]["MaxIters"])
485 )
486
487 # Perform Psi4 minimization...
488 CalcStatus, Energy = MinimizeMoleculeUsingPsi4(OptionsInfo["psi4"], Mol, MolNum, ConfID)
489 if not CalcStatus:
490 if not OptionsInfo["QuietMode"]:
491 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s\n" % (MolName))
492 return [MolIndex, EncodedMol, False, ConfID, Energy]
493
494 return [
495 MolIndex,
496 RDKitUtil.MolToBase64EncodedMolString(
497 Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
498 ),
499 CalcStatus,
500 ConfID,
501 Energy,
502 ]
503
504
505 def InitializePsi4ForWorkerProcess():
506 """Initialize Psi4 for a worker process."""
507
508 if OptionsInfo["Psi4Initialized"]:
509 return
510
511 OptionsInfo["Psi4Initialized"] = True
512
513 if OptionsInfo["MPLevelConformersMode"] and re.match(
514 "auto", OptionsInfo["Psi4RunParams"]["OutputFileSpecified"], re.I
515 ):
516 # Run Psi4 in quiet mode during multiprocessing at Conformers level for 'auto' OutputFile...
517 OptionsInfo["Psi4RunParams"]["OutputFile"] = "quiet"
518 else:
519 # Update output file...
520 OptionsInfo["Psi4RunParams"]["OutputFile"] = Psi4Util.UpdatePsi4OutputFileUsingPID(
521 OptionsInfo["Psi4RunParams"]["OutputFile"], os.getpid()
522 )
523
524 # Intialize Psi4...
525 OptionsInfo["psi4"] = Psi4Util.InitializePsi4(
526 Psi4RunParams=OptionsInfo["Psi4RunParams"],
527 Psi4OptionsParams=OptionsInfo["Psi4OptionsParams"],
528 PrintVersion=False,
529 PrintHeader=True,
530 )
531
532 # Setup max iterations global variable...
533 Psi4Util.UpdatePsi4OptionsParameters(OptionsInfo["psi4"], {"GEOM_MAXITER": OptionsInfo["MaxIters"]})
534
535 # Setup conversion factor for energy units...
536 SetupEnergyConversionFactor(OptionsInfo["psi4"])
537
538
539 def GenerateMolConformers(Psi4Handle, Mol, MolNum):
540 """Generate and mininize conformers for a molecule."""
541
542 # Add hydrogens..
543 Mol = Chem.AddHs(Mol)
544
545 MolName = RDKitUtil.GetMolName(Mol, MolNum)
546
547 # Setup conformers...
548 ConfIDs = EmbedMolecule(Mol, MolNum)
549 if not len(ConfIDs):
550 if not OptionsInfo["QuietMode"]:
551 MiscUtil.PrintWarning(
552 "Conformation generation couldn't be performed for molecule %s: Embedding failed...\n" % MolName
553 )
554 return (Mol, False, None, None)
555
556 # Minimize conformers...
557 CalcEnergyMap = {}
558 for ConfID in ConfIDs:
559 if not OptionsInfo["QuietMode"]:
560 MiscUtil.PrintInfo("Processing conformer ID %s for molecule %s..." % (ConfID, MolName))
561
562 # Perform forcefield minimization...
563 Status, ConvergeStatus = MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID)
564 if not Status:
565 return (Mol, False, None, None)
566
567 if ConvergeStatus != 0:
568 if not OptionsInfo["QuietMode"]:
569 MiscUtil.PrintWarning(
570 'Minimization using forcefield failed to converge for molecule %s in %d steps. Try using higher value for "maxIters" in "--confParams" option...\n'
571 % (MolName, OptionsInfo["ConfGenerationParams"]["MaxIters"])
572 )
573
574 # Perform Psi4 minimization...
575 CalcStatus, Energy = MinimizeMoleculeUsingPsi4(Psi4Handle, Mol, MolNum, ConfID)
576 if not CalcStatus:
577 if not OptionsInfo["QuietMode"]:
578 MiscUtil.PrintWarning("Minimization couldn't be performed for molecule %s\n" % (MolName))
579 return (Mol, False, None, None)
580
581 CalcEnergyMap[ConfID] = Energy
582
583 # Align molecules after minimization...
584 if OptionsInfo["ConfGenerationParams"]["AlignConformers"]:
585 AllChem.AlignMolConformers(Mol)
586
587 # Filter conformers...
588 SelectedConfIDs, SelectedConfEnergies = FilterMolConformers(Mol, MolNum, ConfIDs, CalcEnergyMap)
589
590 return [Mol, True, SelectedConfIDs, SelectedConfEnergies]
591
592
593 def FilterMolConformers(Mol, MolNum, ConfIDs, CalcEnergyMap):
594 """Filter conformers for a molecule."""
595
596 SortedConfIDs = sorted(ConfIDs, key=lambda ConfID: CalcEnergyMap[ConfID])
597
598 MinEnergyConfID = SortedConfIDs[0]
599 MinConfEnergy = CalcEnergyMap[MinEnergyConfID]
600 EnergyWindow = OptionsInfo["EnergyWindow"]
601
602 EnergyRMSDCutoff = OptionsInfo["EnergyRMSDCutoff"]
603 ApplyEnergyRMSDCutoff = False
604 if EnergyRMSDCutoff > 0:
605 ApplyEnergyRMSDCutoff = True
606
607 EnergyRMSDCutoffLowest = OptionsInfo["EnergyRMSDCutoffModeLowest"]
608 EnergyRMSDCalcModeBest = OptionsInfo["EnergyRMSDCalcModeBest"]
609
610 PreAligned = False
611 if OptionsInfo["ConfGenerationParams"]["AlignConformers"]:
612 PreAligned = True
613
614 RefMol, ProbeMol = [None] * 2
615 if EnergyRMSDCalcModeBest:
616 # Copy molecules for best RMSD calculations to avoid change in the coordinates
617 # of the conformations...
618 RefMol = AllChem.Mol(Mol)
619 ProbeMol = AllChem.Mol(Mol)
620
621 # Track conformers with in the specified energy window from the lowest
622 # energy conformation along with applying RMSD cutoff as needed...
623 #
624 SelectedConfIDs = []
625
626 ConfCount = 0
627 IgnoredByEnergyConfCount = 0
628 IgnoredByRMSDConfCount = 0
629
630 FirstConf = True
631
632 for ConfID in SortedConfIDs:
633 if FirstConf:
634 FirstConf = False
635 ConfCount += 1
636 SelectedConfIDs.append(ConfID)
637 continue
638
639 ConfEnergyDiff = abs(CalcEnergyMap[ConfID] - MinConfEnergy)
640 if ConfEnergyDiff > EnergyWindow:
641 IgnoredByEnergyConfCount += 1
642 continue
643
644 if ApplyEnergyRMSDCutoff:
645 IgnoreConf = False
646 if EnergyRMSDCutoffLowest:
647 # Compare RMSD with the lowest energy conformation...
648 if EnergyRMSDCalcModeBest:
649 CalcRMSD = AllChem.GetBestRMS(ProbeMol, RefMol, prbId=ConfID, refId=MinEnergyConfID)
650 else:
651 CalcRMSD = AllChem.GetConformerRMS(Mol, MinEnergyConfID, ConfID, prealigned=PreAligned)
652 if CalcRMSD < EnergyRMSDCutoff:
653 IgnoreConf = True
654 else:
655 for SelectedConfID in SelectedConfIDs:
656 if EnergyRMSDCalcModeBest:
657 CalcRMSD = AllChem.GetBestRMS(ProbeMol, RefMol, prbId=ConfID, refId=SelectedConfID)
658 else:
659 CalcRMSD = AllChem.GetConformerRMS(Mol, SelectedConfID, ConfID, prealigned=PreAligned)
660 if CalcRMSD < EnergyRMSDCutoff:
661 IgnoreConf = True
662 break
663 if IgnoreConf:
664 IgnoredByRMSDConfCount += 1
665 continue
666
667 ConfCount += 1
668 SelectedConfIDs.append(ConfID)
669
670 if not OptionsInfo["QuietMode"]:
671 MiscUtil.PrintInfo(
672 "\nTotal Number of conformations generated for %s: %d" % (RDKitUtil.GetMolName(Mol, MolNum), ConfCount)
673 )
674 MiscUtil.PrintInfo(
675 "Number of conformations ignored due to energy window cutoff: %d" % (IgnoredByEnergyConfCount)
676 )
677 if ApplyEnergyRMSDCutoff:
678 MiscUtil.PrintInfo(
679 "Number of conformations ignored due to energy RMSD cutoff: %d" % (IgnoredByRMSDConfCount)
680 )
681
682 SelectedConfEnergies = None
683 if OptionsInfo["EnergyOut"]:
684 SelectedConfEnergies = []
685 for ConfID in SelectedConfIDs:
686 Energy = "%.*f" % (OptionsInfo["Precision"], CalcEnergyMap[ConfID])
687 SelectedConfEnergies.append(Energy)
688
689 return [SelectedConfIDs, SelectedConfEnergies]
690
691
692 def MinimizeMoleculeUsingPsi4(Psi4Handle, Mol, MolNum, ConfID=-1):
693 """Minimize molecule using Psi4."""
694
695 # Setup a Psi4Mol...
696 Psi4Mol = SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID)
697 if Psi4Mol is None:
698 return (False, None)
699
700 # Setup reference wave function...
701 Reference = SetupReferenceWavefunction(Mol)
702 Psi4Handle.set_options({"Reference": Reference})
703
704 # Setup method name and basis set...
705 MethodName, BasisSet = SetupMethodNameAndBasisSet(Mol)
706
707 # Optimize geometry...
708 Status, Energy, WaveFunction = Psi4Util.PerformGeometryOptimization(
709 Psi4Handle, Psi4Mol, MethodName, BasisSet, ReturnWaveFunction=True, Quiet=OptionsInfo["QuietMode"]
710 )
711
712 if not Status:
713 PerformPsi4Cleanup(Psi4Handle)
714 return (False, None)
715
716 # Update atom positions...
717 AtomPositions = Psi4Util.GetAtomPositions(Psi4Handle, WaveFunction, InAngstroms=True)
718 RDKitUtil.SetAtomPositions(Mol, AtomPositions, ConfID=ConfID)
719
720 # Convert energy units...
721 if OptionsInfo["ApplyEnergyConversionFactor"]:
722 Energy = Energy * OptionsInfo["EnergyConversionFactor"]
723
724 # Clean up
725 PerformPsi4Cleanup(Psi4Handle)
726
727 return (True, Energy)
728
729
730 def MinimizeMoleculeUsingForceField(Mol, MolNum, ConfID=-1):
731 """Minimize molecule using forcefield available in RDKit."""
732
733 try:
734 if OptionsInfo["ConfGenerationParams"]["UseUFF"]:
735 ConvergeStatus = AllChem.UFFOptimizeMolecule(
736 Mol, confId=ConfID, maxIters=OptionsInfo["ConfGenerationParams"]["MaxIters"]
737 )
738 elif OptionsInfo["ConfGenerationParams"]["UseMMFF"]:
739 ConvergeStatus = AllChem.MMFFOptimizeMolecule(
740 Mol,
741 confId=ConfID,
742 maxIters=OptionsInfo["ConfGenerationParams"]["MaxIters"],
743 mmffVariant=OptionsInfo["ConfGenerationParams"]["ForceFieldMMFFVariant"],
744 )
745 else:
746 MiscUtil.PrintError(
747 "Minimization couldn't be performed: Specified forcefield, %s, is not supported"
748 % OptionsInfo["ConfGenerationParams"]["ForceField"]
749 )
750 except (ValueError, RuntimeError, Chem.rdchem.KekulizeException) as ErrMsg:
751 if not OptionsInfo["QuietMode"]:
752 MolName = RDKitUtil.GetMolName(Mol, MolNum)
753 MiscUtil.PrintWarning(
754 "Minimization using forcefield couldn't be performed for molecule %s:\n%s\n" % (MolName, ErrMsg)
755 )
756 return (False, None)
757
758 return (True, ConvergeStatus)
759
760
761 def EmbedMolecule(Mol, MolNum=None):
762 """Embed conformations"""
763
764 ConfIDs = []
765
766 MaxConfs = OptionsInfo["ConfGenerationParams"]["MaxConfs"]
767 RandomSeed = OptionsInfo["ConfGenerationParams"]["RandomSeed"]
768 EnforceChirality = OptionsInfo["ConfGenerationParams"]["EnforceChirality"]
769 UseExpTorsionAnglePrefs = OptionsInfo["ConfGenerationParams"]["UseExpTorsionAnglePrefs"]
770 ETVersion = OptionsInfo["ConfGenerationParams"]["ETVersion"]
771 UseBasicKnowledge = OptionsInfo["ConfGenerationParams"]["UseBasicKnowledge"]
772 EmbedRMSDCutoff = OptionsInfo["ConfGenerationParams"]["EmbedRMSDCutoff"]
773
774 try:
775 ConfIDs = AllChem.EmbedMultipleConfs(
776 Mol,
777 numConfs=MaxConfs,
778 randomSeed=RandomSeed,
779 pruneRmsThresh=EmbedRMSDCutoff,
780 enforceChirality=EnforceChirality,
781 useExpTorsionAnglePrefs=UseExpTorsionAnglePrefs,
782 useBasicKnowledge=UseBasicKnowledge,
783 ETversion=ETVersion,
784 )
785 except ValueError as ErrMsg:
786 if not OptionsInfo["QuietMode"]:
787 MolName = RDKitUtil.GetMolName(Mol, MolNum)
788 MiscUtil.PrintWarning("Embedding failed for molecule %s:\n%s\n" % (MolName, ErrMsg))
789 ConfIDs = []
790
791 if not OptionsInfo["QuietMode"]:
792 if EmbedRMSDCutoff > 0:
793 MiscUtil.PrintInfo(
794 "Generating initial conformation ensemble by distance geometry for %s - EmbedRMSDCutoff: %s; Size: %s; Size after RMSD filtering: %s"
795 % (RDKitUtil.GetMolName(Mol, MolNum), EmbedRMSDCutoff, MaxConfs, len(ConfIDs))
796 )
797 else:
798 MiscUtil.PrintInfo(
799 "Generating initial conformation ensemble by distance geometry for %s - EmbedRMSDCutoff: None; Size: %s"
800 % (RDKitUtil.GetMolName(Mol, MolNum), len(ConfIDs))
801 )
802
803 return ConfIDs
804
805
806 def SetupPsi4Mol(Psi4Handle, Mol, MolNum, ConfID=-1):
807 """Setup a Psi4 molecule object."""
808
809 if OptionsInfo["RecenterAndReorient"]:
810 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, ConfID=ConfID, NoCom=False, NoReorient=False)
811 else:
812 MolGeometry = RDKitUtil.GetPsi4XYZFormatString(Mol, ConfID=ConfID, NoCom=True, NoReorient=True)
813
814 try:
815 Psi4Mol = Psi4Handle.geometry(MolGeometry)
816 except Exception as ErrMsg:
817 Psi4Mol = None
818 if not OptionsInfo["QuietMode"]:
819 MiscUtil.PrintWarning("Failed to create Psi4 molecule from geometry string: %s\n" % ErrMsg)
820 MolName = RDKitUtil.GetMolName(Mol, MolNum)
821 MiscUtil.PrintWarning("Ignoring molecule: %s" % MolName)
822
823 if OptionsInfo["Symmetrize"]:
824 Psi4Mol.symmetrize(OptionsInfo["SymmetrizeTolerance"])
825
826 return Psi4Mol
827
828
829 def PerformPsi4Cleanup(Psi4Handle):
830 """Perform clean up."""
831
832 # Clean up after Psi4 run...
833 Psi4Handle.core.clean()
834
835 # Clean up any leftover scratch files...
836 if OptionsInfo["MPMode"]:
837 Psi4Util.RemoveScratchFiles(Psi4Handle, OptionsInfo["Psi4RunParams"]["OutputFile"])
838
839
840 def CheckAndValidateMolecule(Mol, MolCount=None):
841 """Validate molecule for Psi4 calculations."""
842
843 if Mol is None:
844 if not OptionsInfo["QuietMode"]:
845 MiscUtil.PrintInfo("\nProcessing molecule number %s..." % MolCount)
846 return False
847
848 MolName = RDKitUtil.GetMolName(Mol, MolCount)
849 if not OptionsInfo["QuietMode"]:
850 MiscUtil.PrintInfo("\nProcessing molecule %s..." % MolName)
851
852 if RDKitUtil.IsMolEmpty(Mol):
853 if not OptionsInfo["QuietMode"]:
854 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
855 return False
856
857 if not RDKitUtil.ValidateElementSymbols(RDKitUtil.GetAtomSymbols(Mol)):
858 if not OptionsInfo["QuietMode"]:
859 MiscUtil.PrintWarning("Ignoring molecule containing invalid element symbols: %s\n" % MolName)
860 return False
861
862 return True
863
864
865 def SetupMethodNameAndBasisSet(Mol):
866 """Setup method name and basis set."""
867
868 MethodName = OptionsInfo["MethodName"]
869 if OptionsInfo["MethodNameAuto"]:
870 MethodName = "B3LYP"
871
872 BasisSet = OptionsInfo["BasisSet"]
873 if OptionsInfo["BasisSetAuto"]:
874 BasisSet = "6-31+G**" if RDKitUtil.IsAtomSymbolPresentInMol(Mol, "S") else "6-31G**"
875
876 return (MethodName, BasisSet)
877
878
879 def SetupReferenceWavefunction(Mol):
880 """Setup reference wavefunction."""
881
882 Reference = OptionsInfo["Reference"]
883 if OptionsInfo["ReferenceAuto"]:
884 Reference = "UHF" if (RDKitUtil.GetSpinMultiplicity(Mol) > 1) else "RHF"
885
886 return Reference
887
888
889 def SetupEnergyConversionFactor(Psi4Handle):
890 """Setup converstion factor for energt units. The Psi4 energy units are Hartrees."""
891
892 EnergyUnits = OptionsInfo["EnergyUnits"]
893
894 ApplyConversionFactor = True
895 if re.match(r"^kcal\/mol$", EnergyUnits, re.I):
896 ConversionFactor = Psi4Handle.constants.hartree2kcalmol
897 elif re.match(r"^kJ\/mol$", EnergyUnits, re.I):
898 ConversionFactor = Psi4Handle.constants.hartree2kJmol
899 elif re.match("^eV$", EnergyUnits, re.I):
900 ConversionFactor = Psi4Handle.constants.hartree2ev
901 else:
902 ApplyConversionFactor = False
903 ConversionFactor = 1.0
904
905 OptionsInfo["ApplyEnergyConversionFactor"] = ApplyConversionFactor
906 OptionsInfo["EnergyConversionFactor"] = ConversionFactor
907
908
909 def WriteMolConformers(Writer, Mol, MolNum, ConfIDs, ConfEnergies=None):
910 """Write molecule coformers."""
911
912 if ConfIDs is None:
913 return True
914
915 MolName = RDKitUtil.GetMolName(Mol, MolNum)
916
917 for Index, ConfID in enumerate(ConfIDs):
918 SetConfMolName(Mol, MolName, ConfID)
919
920 if OptionsInfo["EnergyOut"] and ConfEnergies is not None:
921 Mol.SetProp(OptionsInfo["EnergyDataFieldLabel"], ConfEnergies[Index])
922
923 Writer.write(Mol, confId=ConfID)
924
925
926 def SetConfMolName(Mol, MolName, ConfCount):
927 """Set conf mol name."""
928
929 ConfName = "%s_Conf%d" % (MolName, ConfCount)
930 Mol.SetProp("_Name", ConfName)
931
932
933 def ProcessOptions():
934 """Process and validate command line arguments and options."""
935
936 MiscUtil.PrintInfo("Processing options...")
937
938 # Validate options...
939 ValidateOptions()
940
941 OptionsInfo["Infile"] = Options["--infile"]
942 OptionsInfo["SMILESInfileStatus"] = True if MiscUtil.CheckFileExt(Options["--infile"], "smi csv tsv txt") else False
943 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
944 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
945 "--infileParams",
946 Options["--infileParams"],
947 InfileName=Options["--infile"],
948 ParamsDefaultInfo=ParamsDefaultInfoOverride,
949 )
950
951 OptionsInfo["Outfile"] = Options["--outfile"]
952 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
953 "--outfileParams", Options["--outfileParams"]
954 )
955
956 OptionsInfo["Overwrite"] = Options["--overwrite"]
957
958 # Method, basis set, and reference wavefunction...
959 OptionsInfo["BasisSet"] = Options["--basisSet"]
960 OptionsInfo["BasisSetAuto"] = True if re.match("^auto$", Options["--basisSet"], re.I) else False
961
962 OptionsInfo["MethodName"] = Options["--methodName"]
963 OptionsInfo["MethodNameAuto"] = True if re.match("^auto$", Options["--methodName"], re.I) else False
964
965 OptionsInfo["Reference"] = Options["--reference"]
966 OptionsInfo["ReferenceAuto"] = True if re.match("^auto$", Options["--reference"], re.I) else False
967
968 # Run and options parameters...
969 OptionsInfo["Psi4OptionsParams"] = Psi4Util.ProcessPsi4OptionsParameters(
970 "--psi4OptionsParams", Options["--psi4OptionsParams"]
971 )
972 OptionsInfo["Psi4RunParams"] = Psi4Util.ProcessPsi4RunParameters(
973 "--psi4RunParams", Options["--psi4RunParams"], InfileName=OptionsInfo["Infile"]
974 )
975
976 # Conformer generation paramaters...
977 ParamsDefaultInfoOverride = {"MaxConfs": 50, "MaxIters": 250}
978 OptionsInfo["ConfGenerationParams"] = MiscUtil.ProcessOptionConformerParameters(
979 "--confParams", Options["--confParams"], ParamsDefaultInfoOverride
980 )
981
982 # Energy parameters...
983 OptionsInfo["EnergyOut"] = True if re.match("^yes$", Options["--energyOut"], re.I) else False
984 OptionsInfo["EnergyUnits"] = Options["--energyUnits"]
985
986 EnergyDataFieldLabel = Options["--energyDataFieldLabel"]
987 if re.match("^auto$", EnergyDataFieldLabel, re.I):
988 EnergyDataFieldLabel = "Psi4_Energy (%s)" % Options["--energyUnits"]
989 OptionsInfo["EnergyDataFieldLabel"] = EnergyDataFieldLabel
990
991 OptionsInfo["EnergyRMSDCalcMode"] = Options["--energyRMSDCalcMode"]
992 OptionsInfo["EnergyRMSDCalcModeBest"] = (
993 True if re.match("^BestRMSD$", Options["--energyRMSDCalcMode"], re.I) else False
994 )
995
996 OptionsInfo["EnergyRMSDCutoff"] = float(Options["--energyRMSDCutoff"])
997
998 OptionsInfo["EnergyRMSDCutoffMode"] = Options["--energyRMSDCutoffMode"]
999 OptionsInfo["EnergyRMSDCutoffModeLowest"] = (
1000 True if re.match("^Lowest$", Options["--energyRMSDCutoffMode"], re.I) else False
1001 )
1002
1003 if OptionsInfo["EnergyRMSDCutoff"] > 0:
1004 # Make sure that the alignConformers option is being used...
1005 if not OptionsInfo["ConfGenerationParams"]["AlignConformers"]:
1006 MiscUtil.PrintError(
1007 'The value for "alignConformers" specified using "--confParams" must be set to " yes" for non-zero values of "--energyRMSDCutoff" '
1008 )
1009
1010 # Process energy window...
1011 EnergyWindow = Options["--energyWindow"]
1012 if re.match("^auto$", EnergyWindow, re.I):
1013 # Set default energy window based on units...
1014 EnergyUnits = Options["--energyUnits"]
1015 if re.match(r"^kcal\/mol$", EnergyUnits, re.I):
1016 EnergyWindow = 20
1017 elif re.match(r"^kJ\/mol$", EnergyUnits, re.I):
1018 EnergyWindow = 83.68
1019 elif re.match("^eV$", EnergyUnits, re.I):
1020 EnergyWindow = 0.8673
1021 elif re.match("^Hartrees$", EnergyUnits, re.I):
1022 EnergyWindow = 0.03188
1023 else:
1024 MiscUtil.PrintError(
1025 'Failed to set default value for "--energyWindow". The value, %s, specified for option "--energyUnits" is not valid. '
1026 % EnergyUnits
1027 )
1028 else:
1029 EnergyWindow = float(EnergyWindow)
1030 OptionsInfo["EnergyWindow"] = EnergyWindow
1031
1032 OptionsInfo["MaxIters"] = int(Options["--maxIters"])
1033
1034 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1035 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1036
1037 # Multiprocessing level...
1038 MPLevelMoleculesMode = False
1039 MPLevelConformersMode = False
1040 MPLevel = Options["--mpLevel"]
1041 if re.match("^Molecules$", MPLevel, re.I):
1042 MPLevelMoleculesMode = True
1043 elif re.match("^Conformers$", MPLevel, re.I):
1044 MPLevelConformersMode = True
1045 else:
1046 MiscUtil.PrintError('The value, %s, specified for option "--mpLevel" is not valid. ' % MPLevel)
1047 OptionsInfo["MPLevel"] = MPLevel
1048 OptionsInfo["MPLevelMoleculesMode"] = MPLevelMoleculesMode
1049 OptionsInfo["MPLevelConformersMode"] = MPLevelConformersMode
1050
1051 OptionsInfo["Precision"] = int(Options["--precision"])
1052 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1053
1054 OptionsInfo["RecenterAndReorient"] = True if re.match("^yes$", Options["--recenterAndReorient"], re.I) else False
1055
1056 Symmetrize = Options["--symmetrize"]
1057 if re.match("^auto$", Symmetrize, re.I):
1058 Symmetrize = True if OptionsInfo["RecenterAndReorient"] else False
1059 else:
1060 Symmetrize = True if re.match("^yes$", Symmetrize, re.I) else False
1061 OptionsInfo["Symmetrize"] = Symmetrize
1062
1063 OptionsInfo["SymmetrizeTolerance"] = float(Options["--symmetrizeTolerance"])
1064
1065
1066 def RetrieveOptions():
1067 """Retrieve command line arguments and options."""
1068
1069 # Get options...
1070 global Options
1071 Options = docopt(_docoptUsage_)
1072
1073 # Set current working directory to the specified directory...
1074 WorkingDir = Options["--workingdir"]
1075 if WorkingDir:
1076 os.chdir(WorkingDir)
1077
1078 # Handle examples option...
1079 if "--examples" in Options and Options["--examples"]:
1080 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1081 sys.exit(0)
1082
1083
1084 def ValidateOptions():
1085 """Validate option values."""
1086
1087 MiscUtil.ValidateOptionTextValue("--energyOut", Options["--energyOut"], "yes no")
1088 MiscUtil.ValidateOptionTextValue("--energyUnits", Options["--energyUnits"], "Hartrees kcal/mol kJ/mol eV")
1089
1090 MiscUtil.ValidateOptionTextValue(" --energyRMSDCalcMode", Options["--energyRMSDCalcMode"], "RMSD BestRMSD")
1091
1092 MiscUtil.ValidateOptionFloatValue("--energyRMSDCutoff", Options["--energyRMSDCutoff"], {">=": 0})
1093 MiscUtil.ValidateOptionTextValue(" --energyRMSDCutoffMode", Options["--energyRMSDCutoffMode"], "All Lowest")
1094
1095 if not re.match("^auto$", Options["--energyWindow"], re.I):
1096 MiscUtil.ValidateOptionFloatValue("--energyWindow", Options["--energyWindow"], {">": 0})
1097
1098 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1099 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi txt csv tsv")
1100
1101 MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
1102
1103 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
1104 MiscUtil.ValidateOptionsOutputFileOverwrite(
1105 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
1106 )
1107 MiscUtil.ValidateOptionsDistinctFileNames(
1108 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
1109 )
1110
1111 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1112 MiscUtil.ValidateOptionTextValue("--mpLevel", Options["--mpLevel"], "Molecules Conformers")
1113
1114 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
1115 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1116
1117 MiscUtil.ValidateOptionTextValue("--recenterAndReorient", Options["--recenterAndReorient"], "yes no")
1118 MiscUtil.ValidateOptionTextValue("--symmetrize", Options["--symmetrize"], "yes no auto")
1119 MiscUtil.ValidateOptionFloatValue("--symmetrizeTolerance", Options["--symmetrizeTolerance"], {">": 0})
1120
1121
1122 # Setup a usage string for docopt...
1123 _docoptUsage_ = """
1124 Psi4GenerateConformers.py - Generate molecular conformations
1125
1126 Usage:
1127 Psi4GenerateConformers.py [--basisSet <text>] [--confParams <Name,Value,...>] [--energyOut <yes or no>]
1128 [--energyDataFieldLabel <text>] [--energyUnits <text>] [--energyRMSDCalcMode <RMSD or BestRMSD>]
1129 [--energyRMSDCutoff <number>] [--energyRMSDCutoffMode <All or Lowest>] [--energyWindow <number>]
1130 [--infileParams <Name,Value,...>] [--maxIters <number>]
1131 [--methodName <text>] [--mp <yes or no>] [--mpLevel <Molecules or Conformers>]
1132 [--mpParams <Name, Value,...>] [ --outfileParams <Name,Value,...> ] [--overwrite] [--precision <number>]
1133 [--psi4OptionsParams <Name,Value,...>] [--psi4RunParams <Name,Value,...>]
1134 [--quiet <yes or no>] [--reference <text>] [--recenterAndReorient <yes or no>]
1135 [--symmetrize <yes or no>] [--symmetrizeTolerance <number>] [-w <dir>] -i <infile> -o <outfile>
1136 Psi4GenerateConformers.py -h | --help | -e | --examples
1137
1138 Description:
1139 Generate 3D conformers of molecules using a combination of distance geometry
1140 and forcefield minimization followed by geometry optimization using a quantum
1141 chemistry method. A set of initial 3D structures are generated for a molecule
1142 employing distance geometry. The 3D structures in the conformation ensemble
1143 are sequentially minimized using forcefield and a quantum chemistry method.
1144
1145 A Psi4 XYZ format geometry string is automatically generated for each molecule
1146 in input file. It contains atom symbols and 3D coordinates for each atom in a
1147 molecule. In addition, the formal charge and spin multiplicity are present in the
1148 the geometry string. These values are either retrieved from molecule properties
1149 named 'FormalCharge' and 'SpinMultiplicty' or dynamically calculated for a
1150 molecule.
1151
1152 The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
1153 .csv, .tsv .txt)
1154
1155 The supported output file formats are: SD (.sdf, .sd)
1156
1157 Options:
1158 -b, --basisSet <text> [default: auto]
1159 Basis set to use for energy minimization. Default: 6-31+G** for sulfur
1160 containing molecules; Otherwise, 6-31G** [ Ref 150 ]. The specified
1161 value must be a valid Psi4 basis set. No validation is performed.
1162
1163 The following list shows a representative sample of basis sets available
1164 in Psi4:
1165
1166 STO-3G, 6-31G, 6-31+G, 6-31++G, 6-31G*, 6-31+G*, 6-31++G*,
1167 6-31G**, 6-31+G**, 6-31++G**, 6-311G, 6-311+G, 6-311++G,
1168 6-311G*, 6-311+G*, 6-311++G*, 6-311G**, 6-311+G**, 6-311++G**,
1169 cc-pVDZ, cc-pCVDZ, aug-cc-pVDZ, cc-pVDZ-DK, cc-pCVDZ-DK, def2-SVP,
1170 def2-SVPD, def2-TZVP, def2-TZVPD, def2-TZVPP, def2-TZVPPD
1171
1172 --confParams <Name,Value,...> [default: auto]
1173 Generate an initial 3D conformation ensemble using distance geometry and
1174 forcefield minimization before final geometry optimization by a specified
1175 method name and basis set. Possible values: yes or no.
1176
1177 A comma delimited list of parameter name and value pairs for generating
1178 initial sets of 3D conformations for molecules. The 3D conformation ensemble
1179 is generated using distance geometry and forcefield functionality available
1180 in RDKit. The 3D structures in the conformation ensemble are subsequently
1181 minimized by a quantum chemistry method available in Psi4.
1182
1183 The supported parameter names along with their default values are shown
1184 below:
1185
1186 confMethod,ETKDGv2,
1187 forceField,MMFF, forceFieldMMFFVariant,MMFF94,
1188 enforceChirality,yes,alignConformers,yes, embedRMSDCutoff,0.5,
1189 maxConfs,50,maxIters,250,randomSeed,auto
1190
1191 confMethod,ETKDGv2 [ Possible values: SDG, KDG, ETDG,
1192 ETKDG , or ETKDGv2]
1193 forceField, MMFF [ Possible values: UFF or MMFF ]
1194 forceFieldMMFFVariant,MMFF94 [ Possible values: MMFF94 or MMFF94s ]
1195 enforceChirality,yes [ Possible values: yes or no ]
1196 alignConformers,yes [ Possible values: yes or no ]
1197 embedRMSDCutoff,0.5 [ Possible values: number or None]
1198
1199 confMethod: Conformation generation methodology for generating initial 3D
1200 coordinates. Possible values: Standard Distance Geometry (SDG), Experimental
1201 Torsion-angle preference with Distance Geometry (ETDG), basic Knowledge-terms
1202 with Distance Geometry (KDG) and Experimental Torsion-angle preference
1203 along with basic Knowledge-terms with Distance Geometry (ETKDG or
1204 ETKDGv2) [Ref 129, 167] .
1205
1206 forceField: Forcefield method to use for energy minimization. Possible
1207 values: Universal Force Field (UFF) [ Ref 81 ] or Merck Molecular Mechanics
1208 Force Field [ Ref 83-87 ] .
1209
1210 enforceChirality: Enforce chirality for defined chiral centers during
1211 forcefield minimization.
1212
1213 alignConformers: Align conformers for each molecule.
1214
1215 maxConfs: Maximum number of conformations to generate for each molecule
1216 during the generation of an initial 3D conformation ensemble using
1217 conformation generation methodology. The conformations are minimized
1218 using the specified forcefield and a quantum chemistry method. The lowest
1219 energy conformation is written to the output file.
1220
1221 embedRMSDCutoff: RMSD cutoff for retaining initial set of conformers embedded
1222 using distance geometry and before forcefield minimization. All embedded
1223 conformers are kept for 'None' value. Otherwise, only those conformers which
1224 are different from each other by the specified RMSD cutoff, 0.5 by default,
1225 are kept. The first embedded conformer is always retained.
1226
1227 maxIters: Maximum number of iterations to perform for each conformation
1228 during forcefield minimization.
1229
1230 randomSeed: Seed for the random number generator for reproducing initial
1231 3D coordinates in a conformation ensemble. Default is to use a random seed.
1232 --energyOut <yes or no> [default: yes]
1233 Write out energy values.
1234 --energyDataFieldLabel <text> [default: auto]
1235 Energy data field label for writing energy values. Default: Psi4_Energy (<Units>).
1236 --energyUnits <text> [default: kcal/mol]
1237 Energy units. Possible values: Hartrees, kcal/mol, kJ/mol, or eV.
1238 --energyRMSDCalcMode <RMSD or BestRMSD> [default: RMSD]
1239 Methodology for calculating RMSD values during the application of RMSD
1240 cutoff for retaining conformations after the final energy minimization. Possible
1241 values: RMSD or BestRMSD. This option is ignore during 'None' value of
1242 '--energyRMSDCutoff' option.
1243
1244 During BestRMSMode mode, the RDKit 'function AllChem.GetBestRMS' is used to
1245 align and calculate RMSD. This function calculates optimal RMSD for aligning two
1246 molecules, taking symmetry into account. Otherwise, the RMSD value is calculated
1247 using 'AllChem.GetConformerRMS' without changing the atom order. A word to the
1248 wise from RDKit documentation: The AllChem.GetBestRMS function will attempt to
1249 align all permutations of matching atom orders in both molecules, for some molecules
1250 it will lead to 'combinatorial explosion'.
1251 --energyRMSDCutoff <number> [default: 0.5]
1252 RMSD cutoff for retaining conformations after the final energy minimization.
1253 By default, only those conformations which are different from the lowest
1254 energy conformation by the specified RMSD cutoff and are with in the
1255 specified energy window are kept. The lowest energy conformation is always
1256 retained. A value of zero keeps all minimized conformations with in the
1257 specified energy window from the lowest energy.
1258 --energyRMSDCutoffMode <All or Lowest> [default: All]
1259 RMSD cutoff mode for retaining conformations after the final energy
1260 minimization. Possible values: All or Lowest. The RMSD values are compared
1261 against all the selected conformations or the lowest energy conformation during
1262 'All' and 'Lowest' value of '--energyRMSDCutoffMode'. This option is ignored
1263 during 'None' value of --energyRMSDCutoff.
1264
1265 By default, only those conformations which all different from all selected
1266 conformations by the specified RMSD cutoff and are with in the specified
1267 energy window are kept.
1268 --energyWindow <number> [default: auto]
1269 Psi4 Energy window for selecting conformers after the final energy minimization.
1270 The default value is dependent on '--energyUnits': 20 kcal/mol, 83.68 kJ/mol,
1271 0.8673 ev, or 0.03188 Hartrees. The specified value must be in '--energyUnits'.
1272 -e, --examples
1273 Print examples.
1274 -h, --help
1275 Print this help message.
1276 -i, --infile <infile>
1277 Input file name.
1278 --infileParams <Name,Value,...> [default: auto]
1279 A comma delimited list of parameter name and value pairs for reading
1280 molecules from files. The supported parameter names for different file
1281 formats, along with their default values, are shown below:
1282
1283 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
1284 SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
1285 smilesTitleLine,auto,sanitize,yes
1286
1287 Possible values for smilesDelimiter: space, comma or tab.
1288 --maxIters <number> [default: 50]
1289 Maximum number of iterations to perform for each molecule or conformer
1290 during energy minimization by a quantum chemistry method.
1291 -m, --methodName <text> [default: auto]
1292 Method to use for energy minimization. Default: B3LYP [ Ref 150 ]. The
1293 specified value must be a valid Psi4 method name. No validation is
1294 performed.
1295
1296 The following list shows a representative sample of methods available
1297 in Psi4:
1298
1299 B1LYP, B2PLYP, B2PLYP-D3BJ, B2PLYP-D3MBJ, B3LYP, B3LYP-D3BJ,
1300 B3LYP-D3MBJ, CAM-B3LYP, CAM-B3LYP-D3BJ, HF, HF-D3BJ, HF3c, M05,
1301 M06, M06-2x, M06-HF, M06-L, MN12-L, MN15, MN15-D3BJ,PBE, PBE0,
1302 PBEH3c, PW6B95, PW6B95-D3BJ, WB97, WB97X, WB97X-D, WB97X-D3BJ
1303
1304 --mp <yes or no> [default: no]
1305 Use multiprocessing.
1306
1307 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1308 function employing lazy RDKit data iterable. This allows processing of
1309 arbitrary large data sets without any additional requirements memory.
1310
1311 All input data may be optionally loaded into memory by mp.Pool.map()
1312 before starting worker processes in a process pool by setting the value
1313 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1314
1315 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1316 data mode may adversely impact the performance. The '--mpParams' section
1317 provides additional information to tune the value of 'chunkSize'.
1318 --mpLevel <Molecules or Conformers> [default: Molecules]
1319 Perform multiprocessing at molecules or conformers level. Possible values:
1320 Molecules or Conformers. The 'Molecules' value starts a process pool at the
1321 molecules level. All conformers of a molecule are processed in a single
1322 process. The 'Conformers' value, however, starts a process pool at the
1323 conformers level. Each conformer of a molecule is processed in an individual
1324 process in the process pool. The default Psi4 'OutputFile' is set to 'quiet'
1325 using '--psi4RunParams' for 'Conformers' level. Otherwise, it may generate
1326 a large number of Psi4 output files.
1327 --mpParams <Name,Value,...> [default: auto]
1328 A comma delimited list of parameter name and value pairs to configure
1329 multiprocessing.
1330
1331 The supported parameter names along with their default and possible
1332 values are shown below:
1333
1334 chunkSize, auto
1335 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1336 numProcesses, auto [ Default: mp.cpu_count() ]
1337
1338 These parameters are used by the following functions to configure and
1339 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1340 mp.Pool.imap().
1341
1342 The chunkSize determines chunks of input data passed to each worker
1343 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1344 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1345
1346 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1347 automatically converts RDKit data iterable into a list, loads all data into
1348 memory, and calculates the default chunkSize using the following method
1349 as shown in its code:
1350
1351 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1352 if extra: chunkSize += 1
1353
1354 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1355 and 100 data items.
1356
1357 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1358 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1359 data into memory. Consequently, the size of input data is not known a priori.
1360 It's not possible to estimate an optimal value for the chunkSize. The default
1361 chunkSize is set to 1.
1362
1363 The default value for the chunkSize during 'Lazy' data mode may adversely
1364 impact the performance due to the overhead associated with exchanging
1365 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1366 a larger value during 'Lazy' input data mode, based on the size of your input
1367 data and number of processes in the process pool.
1368
1369 The mp.Pool.map() function waits for all worker processes to process all
1370 the data and return the results. The mp.Pool.imap() function, however,
1371 returns the the results obtained from worker processes as soon as the
1372 results become available for specified chunks of data.
1373
1374 The order of data in the results returned by both mp.Pool.map() and
1375 mp.Pool.imap() functions always corresponds to the input data.
1376 -o, --outfile <outfile>
1377 Output file name.
1378 --outfileParams <Name,Value,...> [default: auto]
1379 A comma delimited list of parameter name and value pairs for writing
1380 molecules to files. The supported parameter names for different file
1381 formats, along with their default values, are shown below:
1382
1383 SD: kekulize,yes,forceV3000,no
1384
1385 --overwrite
1386 Overwrite existing files.
1387 --precision <number> [default: 6]
1388 Floating point precision for writing energy values.
1389 --psi4OptionsParams <Name,Value,...> [default: none]
1390 A comma delimited list of Psi4 option name and value pairs for setting
1391 global and module options. The names are 'option_name' for global options
1392 and 'module_name__option_name' for options local to a module. The
1393 specified option names must be valid Psi4 names. No validation is
1394 performed.
1395
1396 The specified option name and value pairs are processed and passed to
1397 psi4.set_options() as a dictionary. The supported value types are float,
1398 integer, boolean, or string. The float value string is converted into a float.
1399 The valid values for a boolean string are yes, no, true, false, on, or off.
1400 --psi4RunParams <Name,Value,...> [default: auto]
1401 A comma delimited list of parameter name and value pairs for configuring
1402 Psi4 jobs.
1403
1404 The supported parameter names along with their default and possible
1405 values are shown below:
1406
1407 MemoryInGB, 1
1408 NumThreads, 1
1409 OutputFile, auto [ Possible values: stdout, quiet, or FileName ]
1410 ScratchDir, auto [ Possivle values: DirName]
1411 RemoveOutputFile, yes [ Possible values: yes, no, true, or false]
1412
1413 These parameters control the runtime behavior of Psi4.
1414
1415 The default file name for 'OutputFile' is <InFileRoot>_Psi4.out. The PID
1416 is appended to output file name during multiprocessing as shown:
1417 <InFileRoot>_Psi4_<PIDNum>.out. The 'stdout' value for 'OutputType'
1418 sends Psi4 output to stdout. The 'quiet' or 'devnull' value suppresses
1419 all Psi4 output. The 'OutputFile' is set to 'quiet' for 'auto' value during
1420 'Conformers' of '--mpLevel' option.
1421
1422 The default 'Yes' value of 'RemoveOutputFile' option forces the removal
1423 of any existing Psi4 before creating new files to append output from
1424 multiple Psi4 runs.
1425
1426 The option 'ScratchDir' is a directory path to the location of scratch
1427 files. The default value corresponds to Psi4 default. It may be used to
1428 override the deafult path.
1429 -q, --quiet <yes or no> [default: no]
1430 Use quiet mode. The warning and information messages will not be printed.
1431 -r, --reference <text> [default: auto]
1432 Reference wave function to use for energy calculation. Default: RHF or UHF.
1433 The default values are Restricted Hartree-Fock (RHF) for closed-shell molecules
1434 with all electrons paired and Unrestricted Hartree-Fock (UHF) for open-shell
1435 molecules with unpaired electrons.
1436
1437 The specified value must be a valid Psi4 reference wave function. No validation
1438 is performed. For example: ROHF, CUHF, RKS, etc.
1439
1440 The spin multiplicity determines the default value of reference wave function
1441 for input molecules. It is calculated from number of free radical electrons using
1442 Hund's rule of maximum multiplicity defined as 2S + 1 where S is the total
1443 electron spin. The total spin is 1/2 the number of free radical electrons in a
1444 molecule. The value of 'SpinMultiplicity' molecule property takes precedence
1445 over the calculated value of spin multiplicity.
1446 --recenterAndReorient <yes or no> [default: yes]
1447 Recenter and reorient a molecule during creation of a Psi4 molecule from
1448 a geometry string.
1449
1450 The 'No' values allows the minimization of a molecule in its initial 3D
1451 coordinate space generated by RDKit.
1452 --symmetrize <yes or no> [default: auto]
1453 Symmetrize molecules before energy minimization. Default: 'Yes' during
1454 'Yes' value of '--recenterAndReorient'; Otherwise, 'No'. The psi4 function,
1455 psi4mol.symmetrize( SymmetrizeTolerance), is called to symmetrize
1456 the molecule before calling psi4.optimize().
1457
1458 The 'No' value of '--symmetrize' during 'Yes' value of '--recenterAndReorient'
1459 may cause psi4.optimize() to fail with a 'Point group changed...' error
1460 message.
1461 --symmetrizeTolerance <number> [default: 0.01]
1462 Symmetry tolerance for '--symmetrize'.
1463 -w, --workingdir <dir>
1464 Location of working directory which defaults to the current directory.
1465
1466 Examples:
1467 To generate an initial conformer ensemble of up to 50 conformations using a
1468 combination of ETKDGv2 distance geometry methodology, applying embed RMSD
1469 cutoff of 0.5 and MMFF forcefield minimization, followed by energy minimization
1470 using B3LYP/6-31G** or B3LYP/6-31+G** (sulfur containing), selecting a final set
1471 of minimized conformers for molecules in a SMILES file, applying energy RMSD
1472 cutoff of 0.5 and energy window value value of 20 kcal/mol, and write out a SD
1473 file containing minimized conformers, type:
1474
1475 % Psi4GenerateConformers.py -i Psi4Sample.smi -o Psi4SampleOut.sdf
1476
1477 To run the first example in a quiet mode and write out a SD file, type:
1478
1479 % Psi4GenerateConformers.py -q yes -i Psi4Sample.smi -o
1480 Psi4SampleOut.sdf
1481
1482 To run the first example in multiprocessing mode at molecules level on all
1483 available CPUs without loading all data into memory and write out a SD file,
1484 type:
1485
1486 % Psi4GenerateConformers.py --mp yes -i Psi4Sample.smi -o
1487 Psi4SampleOut.sdf
1488
1489 To run the first example in multiprocessing mode at conformers level on all
1490 available CPUs without loading all data into memory and write out a SD file,
1491 type:
1492
1493 % Psi4GenerateConformers.py --mp yes --mpLevel Conformers
1494 -i Psi4Sample.smi -o Psi4SampleOut.sdf
1495
1496 To run the first example in multiprocessing mode at molecules level on specific
1497 number of CPUs and chunk size without loading all data into memory and write
1498 out a SD file, type:
1499
1500 % Psi4GenerateConformers.py --mp yes --mpParams "inputDataMode,Lazy,
1501 numProcesses,4,chunkSize,8" -i Psi4Sample.smi -o Psi4SampleOut.sdf
1502
1503 To run the first example by using an explicit set of specific parameters, and
1504 write out a SD file, type
1505
1506 % Psi4GenerateConformers.py --confParams "confMethod,ETKDGv2,
1507 forceField,MMFF, forceFieldMMFFVariant,MMFF94s, maxConfs,20,
1508 embedRMSDCutoff,0.25" --energyUnits "kJ/mol" -m B3LYP
1509 -b "6-31+G**" --maxIters 20 -i Psi4Sample.smi -o Psi4SampleOut.sdf
1510
1511 To run the first example for molecules in a CSV SMILES file, SMILES strings
1512 in column 1, name column 2, and write out a SD file, type:
1513
1514 % Psi4GenerateConformers.py --infileParams "smilesDelimiter,comma,
1515 smilesTitleLine,yes,smilesColumn,1,smilesNameColumn,2"
1516 -i Psi4Sample.csv -o Psi4SampleOut.sdf
1517
1518 Author:
1519
1520 Manish Sud(msud@san.rr.com)
1521
1522 See also:
1523 Psi4CalculateEnergy.py, Psi4CalculatePartialCharges.py, Psi4PerformMinimization.py
1524
1525 Copyright:
1526 Copyright (C) 2026 Manish Sud. All rights reserved.
1527
1528 The functionality available in this script is implemented using Psi4, an
1529 open source quantum chemistry software package, and RDKit, an open
1530 source toolkit for cheminformatics developed by Greg Landrum.
1531
1532 This file is part of MayaChemTools.
1533
1534 MayaChemTools is free software; you can redistribute it and/or modify it under
1535 the terms of the GNU Lesser General Public License as published by the Free
1536 Software Foundation; either version 3 of the License, or (at your option) any
1537 later version.
1538
1539 """
1540
1541 if __name__ == "__main__":
1542 main()
