1 #
2 # File: RDKitUtil.py
3 # Author: Manish Sud <msud@san.rr.com>
4 #
5 # Copyright (C) 2026 Manish Sud. All rights reserved.
6 #
7 # The functionality available in this file is implemented using RDKit, an
8 # open source toolkit for cheminformatics developed by Greg Landrum.
9 #
10 # This file is part of MayaChemTools.
11 #
12 # MayaChemTools is free software; you can redistribute it and/or modify it under
13 # the terms of the GNU Lesser General Public License as published by the Free
14 # Software Foundation; either version 3 of the License, or (at your option) any
15 # later version.
16 #
17 # MayaChemTools is distributed in the hope that it will be useful, but without
18 # any warranty; without even the implied warranty of merchantability of fitness
19 # for a particular purpose. See the GNU Lesser General Public License for more
20 # details.
21 #
22 # You should have received a copy of the GNU Lesser General Public License
23 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
24 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
25 # Boston, MA, 02111-1307, USA.
26 #
27
28 from __future__ import print_function
29
30 import re
31 import base64
32 import json
33
34 from rdkit import Chem
35 from rdkit.Chem import AllChem
36 from rdkit.Chem import Draw
37
38 import MiscUtil
39
40 __all__ = [
41 "AreAtomIndicesSequentiallyConnected",
42 "AreAtomMapNumbersPresentInMol",
43 "AreHydrogensMissingInMolecule",
44 "CalculateFormalCharge",
45 "CalculateSpinMultiplicity",
46 "ClearAtomMapNumbers",
47 "ConstrainAndEmbed",
48 "FilterSubstructureMatchByAtomMapNumbers",
49 "FilterSubstructureMatchesByAtomMapNumbers",
50 "GetAtomIndices",
51 "GetAtomMapIndices",
52 "GetAtomMapIndicesAndMapNumbers",
53 "GetAtomSymbols",
54 "GetAtomPositions",
55 "GetFormalCharge",
56 "GetHeavyAtomNeighbors",
57 "GetInlineSVGForMolecule",
58 "GetInlineSVGForMolecules",
59 "GetMolName",
60 "GetNumFragments",
61 "GetNumHeavyAtomNeighbors",
62 "GetSpinMultiplicity",
63 "GetSVGForMolecule",
64 "GetSVGForMolecules",
65 "GetPsi4XYZFormatString",
66 "GetTorsionsAroundRotatableBonds",
67 "GenerateBase64EncodedMolStrings",
68 "GenerateBase64EncodedMolStringWithConfIDs",
69 "IsAtomSymbolPresentInMol",
70 "IsMolEmpty",
71 "IsValidElementSymbol",
72 "IsValidAtomIndex",
73 "MolFromBase64EncodedMolString",
74 "GenerateBase64EncodedMolStringsWithIDs",
75 "MolToBase64EncodedMolString",
76 "MolFromSubstructureMatch",
77 "MolsFromSubstructureMatches",
78 "ReadMolecules",
79 "ReadAndValidateMolecules",
80 "ReadMoleculesFromSDFile",
81 "ReadMoleculesFromMolFile",
82 "ReadMoleculesFromMol2File",
83 "ReadMoleculesFromPDBFile",
84 "ReadMoleculesFromSMILESFile",
85 "ReorderAtomIndicesInSequentiallyConnectedManner",
86 "SetAtomPositions",
87 "SetupHTMLForTorsionScanViewer",
88 "SetWriterMolProps",
89 "ValidateElementSymbols",
90 "WriteMolecules",
91 ]
92
93
94 def GetMolName(Mol, MolNum=None):
95 """Get molecule name.
96
97 Arguments:
98 Mol (object): RDKit molecule object.
99 MolNum (int or None): Molecule number in input file.
100
101 Returns:
102 str : Molname corresponding to _Name property of a molecule, generated
103 from specieid MolNum using the format "Mol%d" % MolNum, or an
104 empty string.
105
106 """
107
108 MolName = ""
109 if Mol.HasProp("_Name"):
110 MolName = Mol.GetProp("_Name")
111
112 if not len(MolName):
113 if MolNum is not None:
114 MolName = "Mol%d" % MolNum
115
116 return MolName
117
118
119 def GetInlineSVGForMolecule(
120 Mol,
121 Width,
122 Height,
123 Legend=None,
124 AtomListToHighlight=None,
125 BondListToHighlight=None,
126 BoldText=True,
127 Base64Encoded=True,
128 ):
129 """Get SVG image text for a molecule suitable for inline embedding into a HTML page.
130
131 Arguments:
132 Mol (object): RDKit molecule object.
133 Width (int): Width of a molecule image in pixels.
134 Height (int): Height of a molecule image in pixels.
135 Legend (str): Text to display under the image.
136 AtomListToHighlight (list): List of atoms to highlight.
137 BondListToHighlight (list): List of bonds to highlight.
138 BoldText (bool): Flag to make text bold in the image of molecule.
139 Base64Encoded (bool): Flag to return base64 encoded string.
140
141 Returns:
142 str : SVG image text for inline embedding into a HTML page using "img"
143 tag: <img src="data:image/svg+xml;charset=UTF-8,SVGImageText> or
144 tag: <img src="data:image/svg+xml;base64,SVGImageText>
145
146 """
147
148 SVGText = GetSVGForMolecule(Mol, Width, Height, Legend, AtomListToHighlight, BondListToHighlight, BoldText)
149 return _ModifySVGForInlineEmbedding(SVGText, Base64Encoded)
150
151
152 def GetInlineSVGForMolecules(
153 Mols,
154 MolsPerRow,
155 MolWidth,
156 MolHeight,
157 Legends=None,
158 AtomListsToHighlight=None,
159 BondListsToHighLight=None,
160 BoldText=True,
161 Base64Encoded=True,
162 ):
163 """Get SVG image text for molecules suitable for inline embedding into a HTML page.
164
165 Arguments:
166 Mols (list): List of RDKit molecule objects.
167 MolsPerRow (int): Number of molecules per row.
168 Width (int): Width of a molecule image in pixels.
169 Height (int): Height of a molecule image in pixels.
170 Legends (list): List containing strings to display under images.
171 AtomListsToHighlight (list): List of lists containing atoms to highlight
172 for molecules.
173 BondListsToHighlight (list): List of lists containing bonds to highlight
174 for molecules
175 BoldText (bool): Flag to make text bold in the image of molecules.
176 Base64Encoded (bool): Flag to return base64 encoded string.
177
178 Returns:
179 str : SVG image text for inline embedding into a HTML page using "img"
180 tag: <img src="data:image/svg+xml;charset=UTF-8,SVGImageText> or
181 tag: <img src="data:image/svg+xml;base64,SVGImageText>
182
183 """
184
185 SVGText = GetSVGForMolecules(
186 Mols, MolsPerRow, MolWidth, MolHeight, Legends, AtomListsToHighlight, BondListsToHighLight, BoldText
187 )
188 return _ModifySVGForInlineEmbedding(SVGText, Base64Encoded)
189
190
191 def _ModifySVGForInlineEmbedding(SVGText, Base64Encoded):
192 """Modify SVG for inline embedding into a HTML page using "img" tag
193 along with performing base64 encoding.
194 """
195
196 # Take out all tags till the start of '<svg' tag...
197 Pattern = re.compile("^.*<svg", re.I | re.S)
198 SVGText = Pattern.sub("<svg", SVGText)
199
200 # Add an extra space before the "width=..." tag. Otherwise, inline embedding may
201 # cause the following XML error on some browsers due to start of the "width=..."
202 # at the begining of the line in <svg ...> tag:
203 #
204 # XML5607: Whitespace expected.
205 #
206 SVGText = re.sub("width='", " width='", SVGText, flags=re.I)
207
208 # Take out trailing new line...
209 SVGText = SVGText.strip()
210
211 # Perform base64 encoding by turning text into byte stream using string
212 # encode and transform byte stream returned by b64encode into a string
213 # by string decode...
214 #
215 if Base64Encoded:
216 SVGText = base64.b64encode(SVGText.encode()).decode()
217
218 return SVGText
219
220
221 def GetSVGForMolecule(
222 Mol, Width, Height, Legend=None, AtomListToHighlight=None, BondListToHighlight=None, BoldText=True
223 ):
224 """Get SVG image text for a molecule suitable for viewing in a browser.
225
226 Arguments:
227 Mol (object): RDKit molecule object.
228 Width (int): Width of a molecule image in pixels.
229 Height (int): Height of a molecule image in pixels.
230 Legend (str): Text to display under the image.
231 AtomListToHighlight (list): List of atoms to highlight.
232 BondListToHighlight (list): List of bonds to highlight.
233 BoldText (bool): Flag to make text bold in the image of molecule.
234
235 Returns:
236 str : SVG image text for writing to a SVG file for viewing in a browser.
237
238 """
239
240 Mols = [Mol]
241
242 MolsPerRow = 1
243 MolWidth = Width
244 MolHeight = Height
245
246 Legends = [Legend] if Legend is not None else None
247 AtomListsToHighlight = [AtomListToHighlight] if AtomListToHighlight is not None else None
248 BondListsToHighLight = [BondListToHighlight] if BondListToHighlight is not None else None
249
250 return GetSVGForMolecules(
251 Mols, MolsPerRow, MolWidth, MolHeight, Legends, AtomListsToHighlight, BondListsToHighLight, BoldText
252 )
253
254
255 def GetSVGForMolecules(
256 Mols,
257 MolsPerRow,
258 MolWidth,
259 MolHeight,
260 Legends=None,
261 AtomListsToHighlight=None,
262 BondListsToHighlight=None,
263 BoldText=True,
264 ):
265 """Get SVG image text for molecules suitable for viewing in a browser.
266
267 Arguments:
268 Mols (list): List of RDKit molecule objects.
269 MolsPerRow (int): Number of molecules per row.
270 Width (int): Width of a molecule image in pixels.
271 Height (int): Height of a molecule image in pixels.
272 Legends (list): List containing strings to display under images.
273 AtomListsToHighlight (list): List of lists containing atoms to highlight
274 for molecules.
275 BondListsToHighlight (list): List of lists containing bonds to highlight
276 for molecules
277 BoldText (bool): Flag to make text bold in the image of molecules.
278
279 Returns:
280 str : SVG image text for writing to a SVG file for viewing in a browser.
281
282 """
283
284 SVGText = Draw.MolsToGridImage(
285 Mols,
286 molsPerRow=MolsPerRow,
287 subImgSize=(MolWidth, MolHeight),
288 legends=Legends,
289 highlightAtomLists=AtomListsToHighlight,
290 highlightBondLists=BondListsToHighlight,
291 useSVG=True,
292 )
293
294 return _ModifySVGForBrowserViewing(SVGText, BoldText)
295
296
297 def _ModifySVGForBrowserViewing(SVGText, BoldText=True):
298 """Modify SVG for loading into a browser."""
299
300 # It appears that the string 'xmlns:svg' needs to be replaced with 'xmlns' in the
301 # SVG image string generated by older versions of RDKit. Otherwise, the image
302 # doesn't load in web browsers.
303 #
304 if re.search("xmlns:svg", SVGText, re.I):
305 SVGText = re.sub("xmlns:svg", "xmlns", SVGText, flags=re.I)
306
307 # Make text bold...
308 if BoldText:
309 SVGText = re.sub("font-weight:normal;", "font-weight:bold;", SVGText, flags=re.I)
310
311 return SVGText
312
313
314 def IsMolEmpty(Mol):
315 """Check for the presence of atoms in a molecule.
316
317 Arguments:
318 Mol (object): RDKit molecule object.
319
320 Returns:
321 bool : True - No atoms in molecule; Otherwise, false.
322
323 """
324
325 Status = False if Mol.GetNumAtoms() else True
326
327 return Status
328
329
330 def IsAtomSymbolPresentInMol(Mol, AtomSymbol, IgnoreCase=True):
331 """Check for the presence of an atom symbol in a molecule.
332
333 Arguments:
334 Mol (object): RDKit molecule object.
335 AtomSymbol (str): Atom symbol.
336
337 Returns:
338 bool : True - Atom symbol in molecule; Otherwise, false.
339
340 """
341
342 for Atom in Mol.GetAtoms():
343 Symbol = Atom.GetSymbol()
344 if IgnoreCase:
345 if re.match("^%s$" % AtomSymbol, Symbol, re.I):
346 return True
347 else:
348 if re.match("^%s$" % AtomSymbol, Symbol):
349 return True
350
351 return False
352
353
354 def ValidateElementSymbols(ElementSymbols):
355 """Validate element symbols.
356
357 Arguments:
358 ElementSymbols (list): List of element symbols to validate.
359
360 Returns:
361 bool : True - All element symbols are valid; Otherwise, false.
362
363 """
364 for ElementSymbol in ElementSymbols:
365 if not IsValidElementSymbol(ElementSymbol):
366 return False
367
368 return True
369
370
371 def GetAtomPositions(Mol, ConfID=-1):
372 """Retrieve a list of lists containing coordinates of all atoms in a
373 molecule.
374
375 Arguments:
376 Mol (object): RDKit molecule object.
377 ConfID (int): Conformer number.
378
379 Returns:
380 list : List of lists containing atom positions.
381
382 Examples:
383
384 for AtomPosition in RDKitUtil.GetAtomPositions(Mol):
385 print("X: %s; Y: %s; Z: %s" % (AtomPosition[0], AtomPosition[1], AtomPosition[2]))
386
387 """
388
389 return Mol.GetConformer(id=ConfID).GetPositions().tolist()
390
391
392 def SetAtomPositions(Mol, AtomPositions, ConfID=-1):
393 """Set atom positions of all atoms in a molecule.
394
395 Arguments:
396 Mol (object): RDKit molecule object.
397 AtomPositions (object): List of lists containing atom positions.
398 ConfID (int): Conformer number.
399
400 Returns:
401 object : RDKit molecule object.
402
403 """
404
405 MolConf = Mol.GetConformer(ConfID)
406
407 for Index in range(len(AtomPositions)):
408 MolConf.SetAtomPosition(Index, tuple(AtomPositions[Index]))
409
410 return Mol
411
412
413 def GetAtomSymbols(Mol):
414 """Retrieve a list containing atom symbols of all atoms a molecule.
415
416 Arguments:
417 Mol (object): RDKit molecule object.
418
419 Returns:
420 list : List of atom symbols.
421
422 """
423
424 return [Atom.GetSymbol() for Atom in Mol.GetAtoms()]
425
426
427 def GetAtomIndices(Mol):
428 """Retrieve a list containing atom indices of all atoms a molecule.
429
430 Arguments:
431 Mol (object): RDKit molecule object.
432
433 Returns:
434 list : List of atom indices.
435
436 """
437
438 return [Atom.GetIdx() for Atom in Mol.GetAtoms()]
439
440
441 def GetFormalCharge(Mol, CheckMolProp=True):
442 """Get formal charge of a molecule. The formal charge is either retrieved
443 from 'FormalCharge' molecule property or calculated using RDKit function
444 Chem.GetFormalCharge(Mol).
445
446 The 'FormalCharge' molecule property may contain multiple space delimited
447 values. The total formal charge corresponds to the sum of the specified formal
448 charge values.
449
450 Arguments:
451 Mol (object): RDKit molecule object.
452 CheckMolProp (bool): Check 'FormalCharge' molecule property to
453 retrieve formal charge.
454
455 Returns:
456 int : Formal charge.
457
458 """
459
460 Name = "FormalCharge"
461 if CheckMolProp and Mol.HasProp(Name):
462 FormalCharge = Mol.GetProp(Name)
463 Values = FormalCharge.split()
464 if len(Values) > 1:
465 MiscUtil.PrintWarning(
466 "RDKitUtil.GetFormalCharge: Molecule property, %s, contains multiple values, %s. Formal charge corresponds to sum of the specified values..."
467 % (Name, FormalCharge)
468 )
469 FormalCharge = 0.0
470 for Value in Values:
471 FormalCharge += float(Value)
472 FormalCharge = int(FormalCharge)
473 else:
474 FormalCharge = int(float(FormalCharge))
475 else:
476 FormalCharge = CalculateFormalCharge(Mol)
477
478 return int(FormalCharge)
479
480
481 def CalculateFormalCharge(Mol):
482 """Calculate formal charge of a molecule. The formal charge is calculated
483 using RDKit function Chem.GetFormalCharge(Mol).
484
485 Arguments:
486 Mol (object): RDKit molecule object.
487 retrieve formal charge.
488
489 Returns:
490 int : Formal charge.
491
492 """
493
494 return int(Chem.GetFormalCharge(Mol))
495
496
497 def GetSpinMultiplicity(Mol, CheckMolProp=True):
498 """Get spin multiplicity of a molecule. The spin multiplicity is either
499 retrieved from 'SpinMultiplicity' molecule property or calculated
500 from the number of free radical electrons using Hund's rule of maximum
501 multiplicity defined as 2S + 1 where S is the total electron spin. The
502 total spin is 1/2 the number of free radical electrons in a molecule.
503
504 The 'SpinMultiplicity' molecule property may contain multiple space delimited
505 values. The total spin multiplicity corresponds to the total number of free radical
506 electrons which are calculated for each specified value.
507
508 Arguments:
509 Mol (object): RDKit molecule object.
510 CheckMolProp (bool): Check 'SpinMultiplicity' molecule property to
511 retrieve spin multiplicity.
512
513 Returns:
514 int : Spin multiplicity.
515
516 """
517
518 Name = "SpinMultiplicity"
519 if CheckMolProp and Mol.HasProp(Name):
520 SpinMultiplicity = Mol.GetProp(Name)
521 Values = SpinMultiplicity.split()
522 if len(Values) > 1:
523 MiscUtil.PrintWarning(
524 "RDKitUtil.GetSpinMultiplicity: Molecule property, %s, contains multiple values, %s. Calculating spin multiplicity corresponding to total number of free radical electrons for each specified value..."
525 % (Name, SpinMultiplicity)
526 )
527 NumRadicalElectrons = 0
528 for Value in Values:
529 NumRadicalElectrons += int(float(Value)) - 1
530
531 TotalElectronicSpin = NumRadicalElectrons / 2
532 SpinMultiplicity = 2 * TotalElectronicSpin + 1
533 else:
534 SpinMultiplicity = int(float(SpinMultiplicity))
535 else:
536 SpinMultiplicity = CalculateSpinMultiplicity(Mol)
537
538 return int(SpinMultiplicity)
539
540
541 def CalculateSpinMultiplicity(Mol):
542 """Calculate spin multiplicity of a molecule. The spin multiplicity is calculated
543 from the number of free radical electrons using Hund's rule of maximum
544 multiplicity defined as 2S + 1 where S is the total electron spin. The
545 total spin is 1/2 the number of free radical electrons in a molecule.
546
547 Arguments:
548 Mol (object): RDKit molecule object.
549
550 Returns:
551 int : Spin multiplicity.
552
553 """
554
555 # Calculate spin multiplicity using Hund's rule of maximum multiplicity...
556 NumRadicalElectrons = 0
557 for Atom in Mol.GetAtoms():
558 NumRadicalElectrons += Atom.GetNumRadicalElectrons()
559
560 TotalElectronicSpin = NumRadicalElectrons / 2
561 SpinMultiplicity = 2 * TotalElectronicSpin + 1
562
563 return int(SpinMultiplicity)
564
565
566 def GetPsi4XYZFormatString(
567 Mol,
568 ConfID=-1,
569 FormalCharge="auto",
570 SpinMultiplicity="auto",
571 Symmetry="auto",
572 NoCom=False,
573 NoReorient=False,
574 CheckFragments=False,
575 ):
576 """Retrieve geometry string of a molecule in Psi4ish XYZ format to perform
577 Psi4 quantum chemistry calculations.
578
579 You may explicit specify multiple space delimited values for formal charge
580 and spin multiplicity. Otherwise, these values are either automatically
581 retrieved from 'FormalCharge' and 'SpinMultiplicity' molecule properties or
582 calculated using RDKit. The number of specified values for these properties
583 must match the number of fragments in the molecule during the processing
584 of the fragments.
585
586 Arguments:
587 Mol (object): RDKit molecule object.
588 ConfID (int): Conformer number.
589 FormalCharge (str): Specified formal charge or 'auto' to calculate
590 its value by RDKit.
591 SpinMultiplicity (str): Specified spin multiplicity or 'auto' to calculate
592 its value by RDKit.
593 Symmetry (str): Specified symmetry or 'auto' to calculate its value by
594 Psi4.
595 NoCom (bool): Flag to disable recentering of a molecule by Psi4.
596 NoReorient (bool): Flag to disable reorientation of a molecule by Psi4.
597 CheckFragments (bool): Check for fragments and setup geometry string
598 using -- separator between fragments.
599
600 Returns:
601 str : Geometry string of a molecule in Psi4ish XYZ format.
602
603 """
604
605 # Check for fragments...
606 Mols = [Mol]
607 if CheckFragments:
608 Fragments = list(Chem.rdmolops.GetMolFrags(Mol, asMols=True))
609 if len(Fragments) > 1:
610 Mols = Fragments
611
612 FragMolFormalCharges = _SetFormalChargesForPsi4XYZFormatString(Mol, Mols, FormalCharge, CheckFragments)
613 FragMolSpinMultiplicities = _SetSpinMultiplicitiesForPsi4XYZFormatString(
614 Mol, Mols, SpinMultiplicity, CheckFragments
615 )
616
617 # Setup geometry string for Ps4...
618 GeometryList = []
619 FragMolCount = 0
620
621 for FragMolIndex, FragMol in enumerate(Mols):
622 FragMolCount += 1
623 if FragMolCount > 1:
624 GeometryList.append("--")
625
626 FragMolFormalCharge = FragMolFormalCharges[FragMolIndex]
627 FragMolSpinMultiplicity = FragMolSpinMultiplicities[FragMolIndex]
628 if FragMolFormalCharge is None or FragMolSpinMultiplicity is None:
629 MiscUtil.PrintInfo("")
630 MiscUtil.PrintWarning(
631 "RDKitUtil.GetPsi4XYZFormatString: Failed to set formal charge and spin multiplicity values. Both formal charge, %s, and spin multiplicity, %s, must be valid values. These values are either specified explicitly or automatically calculated..."
632 % (FragMolFormalCharge, FragMolSpinMultiplicity)
633 )
634 else:
635 GeometryList.append("%s %s" % (FragMolFormalCharge, FragMolSpinMultiplicity))
636
637 AtomSymbols = GetAtomSymbols(FragMol)
638 AtomPositions = GetAtomPositions(FragMol, ConfID)
639
640 for AtomSymbol, AtomPosition in zip(AtomSymbols, AtomPositions):
641 GeometryList.append("%s %s %s %s" % (AtomSymbol, AtomPosition[0], AtomPosition[1], AtomPosition[2]))
642
643 GeometryList.append("units angstrom")
644
645 if not re.match("^auto$", Symmetry, re.I):
646 Name = "Symmetry"
647 if Mol.HasProp(Name):
648 Symmetry = Mol.GetProp(Name)
649 GeometryList.append("symmetry %s" % Symmetry)
650
651 if NoCom:
652 GeometryList.append("no_com")
653
654 if NoReorient:
655 GeometryList.append("no_reorient")
656
657 Geometry = "\n".join(GeometryList)
658
659 return Geometry
660
661
662 def _SetFormalChargesForPsi4XYZFormatString(Mol, FragMols, FormalCharge, CheckFragments):
663 """Setup formal charges for Psi4 XYZ format string."""
664
665 if not CheckFragments:
666 if re.match("^auto$", FormalCharge, re.I):
667 MolFormalCharge = GetFormalCharge(Mol)
668 else:
669 MolFormalCharge = int(FormalCharge)
670 return [MolFormalCharge]
671
672 FragMolsCount = len(FragMols)
673 FormalCharges = [None] * FragMolsCount
674
675 if re.match("^auto$", FormalCharge, re.I):
676 PropName = "FormalCharge"
677 if Mol.HasProp(PropName):
678 FormalCharge = Mol.GetProp(PropName)
679 FormalChargeWords = FormalCharge.split()
680 if len(FormalChargeWords) == FragMolsCount:
681 FormalCharges = [int(float(FormalCharge)) for FormalCharge in FormalChargeWords]
682 else:
683 MiscUtil.PrintWarning(
684 "RDKitUtil.GetPsi4XYZFormatString: Ignoring specified value, %s, for FormalCharge molecule property. The number of space delimted specified values, %s, must match number of fragments, %s, in the molecule..."
685 % (FormalCharge, len(FormalChargeWords), FragMolsCount)
686 )
687 else:
688 FormalCharges = [CalculateFormalCharge(FragMol) for FragMol in FragMols]
689 else:
690 FormalChargeWords = FormalCharge.split()
691 if len(FormalChargeWords) != FragMolsCount:
692 MiscUtil.PrintWarning(
693 "RDKitUtil.GetPsi4XYZFormatString: Ignoring specified value, %s, for FormalCharge paramater. The number of space delimted specified values, %s, must match number of fragments, %s, in the molecule..."
694 % (FormalCharge, len(FormalChargeWords), FragMolsCount)
695 )
696 else:
697 FormalCharges = [int(FormalCharge) for FormalCharge in FormalChargeWords]
698
699 return FormalCharges
700
701
702 def _SetSpinMultiplicitiesForPsi4XYZFormatString(Mol, FragMols, SpinMultiplicity, CheckFragments):
703 """Setup spin multiplicites for Psi4 XYZ format string."""
704
705 if not CheckFragments:
706 if re.match("^auto$", SpinMultiplicity, re.I):
707 MolSpinMultiplicity = GetSpinMultiplicity(Mol)
708 else:
709 MolSpinMultiplicity = int(SpinMultiplicity)
710 return [MolSpinMultiplicity]
711
712 FragMolsCount = len(FragMols)
713 SpinMultiplicities = [None] * FragMolsCount
714
715 if re.match("^auto$", SpinMultiplicity, re.I):
716 PropName = "SpinMultiplicity"
717 if Mol.HasProp(PropName):
718 SpinMultiplicity = Mol.GetProp(PropName)
719 SpinMultiplicityWords = SpinMultiplicity.split()
720 if len(SpinMultiplicityWords) == FragMolsCount:
721 SpinMultiplicities = [int(float(SpinMultiplicity)) for SpinMultiplicity in SpinMultiplicityWords]
722 else:
723 MiscUtil.PrintWarning(
724 "RDKitUtil.GetPsi4XYZFormatString: Ignoring specified value, %s, for SpinMultiplicity molecule property. The number of space delimted specified values, %s, must match number of fragments, %s, in the molecule..."
725 % (SpinMultiplicity, len(SpinMultiplicityWords), FragMolsCount)
726 )
727 else:
728 SpinMultiplicities = [CalculateSpinMultiplicity(FragMol) for FragMol in FragMols]
729 else:
730 SpinMultiplicityWords = SpinMultiplicity.split()
731 if len(SpinMultiplicityWords) != FragMolsCount:
732 MiscUtil.PrintWarning(
733 "RDKitUtil.GetPsi4XYZFormatString: Ignoring specified value, %s, for SpinMultiplicity paramater. The number of space delimted specified values, %s, must match number of fragments, %s, in the molecule..."
734 % (SpinMultiplicity, len(SpinMultiplicityWords), FragMolsCount)
735 )
736 else:
737 SpinMultiplicities = [int(SpinMultiplicity) for SpinMultiplicity in SpinMultiplicityWords]
738
739 return SpinMultiplicities
740
741
742 def GetTorsionsAroundRotatableBonds(Mol, RotBondsPatternMol, IgnoreHydrogens=True):
743 """Identify torsions around rotatable bonds and return a list of lists
744 containing atom indices of torsions.
745
746 Arguments:
747 Mol (object): RDKit molecule object.
748 PatternMol (object): RDKit molecule object corresponding to SMARTS
749 pattern to identify rotatable bonds.
750 IgnoreHydrogens (bool): Flag to include torsions around rotatable bonds
751 containing hydrogens.
752
753 Returns:
754 list : List of lists containing atom indices of torsions around
755 rotatable bonds.
756
757 """
758
759 # Match rotatable bonds...
760 RotBondsMatches = FilterSubstructureMatchesByAtomMapNumbers(
761 Mol, RotBondsPatternMol, Mol.GetSubstructMatches(RotBondsPatternMol, useChirality=False)
762 )
763 if not len(RotBondsMatches):
764 return None
765
766 # Identify torsions...
767 TorsioAtomIndicesList = []
768 for RotBondMatch in RotBondsMatches:
769 if len(RotBondMatch) != 2:
770 continue
771
772 Atom1Index, Atom2Index = RotBondMatch
773 Atom1NbrIndices = _GetRotBondAtomNeighbors(Mol, Atom1Index, Atom2Index, IgnoreHydrogens)
774 Atom2NbrIndices = _GetRotBondAtomNeighbors(Mol, Atom2Index, Atom1Index, IgnoreHydrogens)
775
776 TorsionAtomIndices = []
777 for Atom1NbrIndex in Atom1NbrIndices:
778 if Atom1NbrIndex in Atom2NbrIndices:
779 continue
780 for Atom2NbrIndex in Atom2NbrIndices:
781 if Atom2NbrIndex in Atom1NbrIndices:
782 continue
783 TorsionAtomIndices.append([Atom1NbrIndex, Atom1Index, Atom2Index, Atom2NbrIndex])
784
785 TorsioAtomIndicesList.extend(TorsionAtomIndices)
786
787 return TorsioAtomIndicesList
788
789
790 def _GetRotBondAtomNeighbors(Mol, AtomIndex, BondedAtomIndex, IgnoreHydrogens):
791 """Get atom neigbors around a rotatable bond."""
792
793 Atom = Mol.GetAtomWithIdx(AtomIndex)
794 AtomNeighbors = []
795
796 for AtomNbr in Atom.GetNeighbors():
797 AtomNbrIndex = AtomNbr.GetIdx()
798 if AtomNbrIndex == BondedAtomIndex:
799 continue
800 if IgnoreHydrogens:
801 if AtomNbr.GetAtomicNum() == 1:
802 continue
803
804 AtomNeighbors.append(AtomNbrIndex)
805
806 return AtomNeighbors
807
808
809 def GetNumFragments(Mol):
810 """Get number of fragment in a molecule.
811
812 Arguments:
813 Atom (object): RDKit molecule object.
814
815 Returns:
816 int : Number of fragments.
817
818 """
819
820 Fragments = Chem.rdmolops.GetMolFrags(Mol, asMols=False)
821
822 return len(Fragments) if Fragments is not None else 0
823
824
825 def GetNumHeavyAtomNeighbors(Atom):
826 """Get number of heavy atom neighbors.
827
828 Arguments:
829 Atom (object): RDKit atom object.
830
831 Returns:
832 int : Number of neighbors.
833
834 """
835
836 NbrCount = 0
837 for AtomNbr in Atom.GetNeighbors():
838 if AtomNbr.GetAtomicNum() > 1:
839 NbrCount += 1
840
841 return NbrCount
842
843
844 def GetHeavyAtomNeighbors(Atom):
845 """Get a list of heavy atom neighbors.
846
847 Arguments:
848 Atom (object): RDKit atom object.
849
850 Returns:
851 list : List of heavy atom neighbors.
852
853 """
854
855 AtomNeighbors = []
856 for AtomNbr in Atom.GetNeighbors():
857 if AtomNbr.GetAtomicNum() > 1:
858 AtomNeighbors.append(AtomNbr)
859
860 return AtomNeighbors
861
862
863 def IsValidElementSymbol(ElementSymbol):
864 """Validate element symbol.
865
866 Arguments:
867 ElementSymbol (str): Element symbol
868
869 Returns:
870 bool : True - Valid element symbol; Otherwise, false.
871
872 """
873
874 try:
875 AtomicNumber = Chem.GetPeriodicTable().GetAtomicNumber(ElementSymbol)
876 Status = True if AtomicNumber > 0 else False
877 except Exception:
878 Status = False
879
880 return Status
881
882
883 def IsValidAtomIndex(Mol, AtomIndex):
884 """Validate presence atom index in a molecule.
885
886 Arguments:
887 Mol (object): RDKit molecule object.
888 AtomIndex (int): Atom index.
889
890 Returns:
891 bool : True - Valid atom index; Otherwise, false.
892
893 """
894 for Atom in Mol.GetAtoms():
895 if AtomIndex == Atom.GetIdx():
896 return True
897
898 return False
899
900
901 def AreHydrogensMissingInMolecule(Mol):
902 """Check for any missing hydrogens in in a molecue.
903
904 Arguments:
905 Mol (object): RDKit molecule object.
906
907 Returns:
908 bool : True - Missing hydrogens; Otherwise, false.
909
910 """
911
912 for Atom in Mol.GetAtoms():
913 NumExplicitAndImplicitHs = Atom.GetNumExplicitHs() + Atom.GetNumImplicitHs()
914 if NumExplicitAndImplicitHs > 0:
915 return True
916
917 return False
918
919
920 def AreAtomIndicesSequentiallyConnected(Mol, AtomIndices):
921 """Check for the presence bonds between sequential pairs of atoms in a
922 molecule.
923
924 Arguments:
925 Mol (object): RDKit molecule object.
926 AtomIndices (list): List of atom indices.
927
928 Returns:
929 bool : True - Sequentially connected; Otherwise, false.
930
931 """
932
933 for Index in range(0, (len(AtomIndices) - 1)):
934 Bond = Mol.GetBondBetweenAtoms(AtomIndices[Index], AtomIndices[Index + 1])
935 if Bond is None:
936 return False
937
938 if Bond.GetIdx() is None:
939 return False
940
941 return True
942
943
944 def ReorderAtomIndicesInSequentiallyConnectedManner(Mol, AtomIndices):
945 """Check for the presence of sequentially connected list of atoms in an
946 arbitray list of atoms in molecule.
947
948 Arguments:
949 Mol (object): RDKit molecule object.
950 AtomIndices (list): List of atom indices.
951
952 Returns:
953 bool : True - Sequentially connected list found; Otherwise, false.
954 list : List of seqeuntially connected atoms or None.
955
956 """
957
958 # Count the number of neighbors for specified atom indices ensuring
959 # that the neighbors are also part of atom indices...
960 AtomNbrsCount = {}
961 for AtomIndex in AtomIndices:
962 Atom = Mol.GetAtomWithIdx(AtomIndex)
963
964 AtomNbrsCount[AtomIndex] = 0
965 for AtomNbr in Atom.GetNeighbors():
966 AtomNbrIndex = AtomNbr.GetIdx()
967 if AtomNbrIndex not in AtomIndices:
968 continue
969 AtomNbrsCount[AtomIndex] += 1
970
971 # Number of neighbors for each specified atom indices must be 1 or 2
972 # for sequentially connected list of atom indices...
973 AtomsWithOneNbr = []
974 for AtomIndex, NbrsCount in AtomNbrsCount.items():
975 if not (NbrsCount == 1 or NbrsCount == 2):
976 return (False, None)
977
978 if NbrsCount == 1:
979 AtomsWithOneNbr.append(AtomIndex)
980
981 # A sequentially connected list of indices must have two atom indices with
982 # exactly # one neighbor...
983 if len(AtomsWithOneNbr) != 2:
984 return (False, None)
985
986 # Setup a reordered list of sequentially connected atoms...
987 ReorderedAtomIndices = []
988
989 AtomIndex1, AtomIndex2 = AtomsWithOneNbr
990 AtomIndex = AtomIndex1 if AtomIndex1 < AtomIndex2 else AtomIndex2
991 ReorderedAtomIndices.append(AtomIndex)
992
993 while len(ReorderedAtomIndices) < len(AtomIndices):
994 Atom = Mol.GetAtomWithIdx(AtomIndex)
995
996 for AtomNbr in Atom.GetNeighbors():
997 AtomNbrIndex = AtomNbr.GetIdx()
998 if AtomNbrIndex not in AtomIndices:
999 continue
1000
1001 if AtomNbrIndex in ReorderedAtomIndices:
1002 continue
1003
1004 # Treat neighbor as next connected atom...
1005 AtomIndex = AtomNbrIndex
1006 ReorderedAtomIndices.append(AtomIndex)
1007 break
1008
1009 # Check reorderd list size...
1010 if len(ReorderedAtomIndices) != len(AtomIndices):
1011 return (False, None)
1012
1013 # A final check to validate reorderd list...
1014 if not AreAtomIndicesSequentiallyConnected(Mol, ReorderedAtomIndices):
1015 return (False, None)
1016
1017 return (True, ReorderedAtomIndices)
1018
1019
1020 def MolToBase64EncodedMolString(Mol, PropertyPickleFlags=Chem.PropertyPickleOptions.AllProps):
1021 """Encode RDkit molecule object into a base64 encoded string. The properties
1022 can be optionally excluded.
1023
1024 The molecule is pickled using RDKit Mol.ToBinary() function before
1025 their encoding.
1026
1027 Arguments:
1028 Mol (object): RDKit molecule object.
1029 PropertyPickleFlags: RDKit property pickle options.
1030
1031 Returns:
1032 str : Base64 encode molecule string or None.
1033
1034 Notes:
1035 The following property pickle flags are currently available in RDKit:
1036
1037 Chem.PropertyPickleOptions.NoProps
1038 Chem.PropertyPickleOptions.MolProps
1039 Chem.PropertyPickleOptions.AtomProps
1040 Chem.PropertyPickleOptions.BondProps
1041 Chem.PropertyPickleOptions.PrivateProps
1042 Chem.PropertyPickleOptions.AllProps
1043
1044 """
1045
1046 return None if Mol is None else base64.b64encode(Mol.ToBinary(PropertyPickleFlags)).decode()
1047
1048
1049 def MolFromBase64EncodedMolString(EncodedMol):
1050 """Generate a RDKit molecule object from a base64 encoded string.
1051
1052 Arguments:
1053 str: Base64 encoded molecule string.
1054
1055 Returns:
1056 object : RDKit molecule object or None.
1057
1058 """
1059
1060 return None if EncodedMol is None else Chem.Mol(base64.b64decode(EncodedMol))
1061
1062
1063 def GenerateBase64EncodedMolStrings(Mols, PropertyPickleFlags=Chem.PropertyPickleOptions.AllProps):
1064 """Setup an iterator for generating base64 encoded molecule string
1065 from a RDKit molecule iterator. The iterator returns a list containing
1066 a molecule index and encoded molecule string or None.
1067
1068 The molecules are pickled using RDKit Mol.ToBinary() function
1069 before their encoding.
1070
1071 Arguments:
1072 iterator: RDKit molecules iterator.
1073 PropertyFlags: RDKit property pickle options.
1074
1075 Returns:
1076 object : Base64 endcoded molecules iterator. The iterator returns a
1077 list containing a molecule index and an encoded molecule string
1078 or None.
1079
1080 Notes:
1081 The following property pickle flags are currently available in RDKit:
1082
1083 Chem.PropertyPickleOptions.NoProps
1084 Chem.PropertyPickleOptions.MolProps
1085 Chem.PropertyPickleOptions.AtomProps
1086 Chem.PropertyPickleOptions.BondProps
1087 Chem.PropertyPickleOptions.PrivateProps
1088 Chem.PropertyPickleOptions.AllProps
1089
1090 Examples:
1091
1092 EncodedMolsInfo = GenerateBase64EncodedMolStrings(Mols)
1093 for MolIndex, EncodedMol in EncodedMolsInfo:
1094 if EncodeMol is not None:
1095 Mol = MolFromBase64EncodedMolString(EncodedMol)
1096
1097 """
1098 for MolIndex, Mol in enumerate(Mols):
1099 yield [MolIndex, None] if Mol is None else [MolIndex, MolToBase64EncodedMolString(Mol, PropertyPickleFlags)]
1100
1101
1102 def GenerateBase64EncodedMolStringsWithIDs(Mols, MolIDs, PropertyPickleFlags=Chem.PropertyPickleOptions.AllProps):
1103 """Setup an iterator for generating base64 encoded molecule string
1104 from a RDKit molecule iterator. The iterator returns a list containing
1105 a molecule ID and encoded molecule string or None.
1106
1107 The molecules are pickled using RDKit Mol.ToBinary() function
1108 before their encoding.
1109
1110 Arguments:
1111 iterator: RDKit molecules iterator.
1112 MolIDs (list): Molecule IDs.
1113 PropertyFlags: RDKit property pickle options.
1114
1115 Returns:
1116 object : Base64 endcoded molecules iterator. The iterator returns a
1117 list containing a molecule ID and an encoded molecule string
1118 or None.
1119
1120 Notes:
1121 The following property pickle flags are currently available in RDKit:
1122
1123 Chem.PropertyPickleOptions.NoProps
1124 Chem.PropertyPickleOptions.MolProps
1125 Chem.PropertyPickleOptions.AtomProps
1126 Chem.PropertyPickleOptions.BondProps
1127 Chem.PropertyPickleOptions.PrivateProps
1128 Chem.PropertyPickleOptions.AllProps
1129
1130 Examples:
1131
1132 EncodedMolsInfo = GenerateBase64EncodedMolStringsWithIDs(Mols)
1133 for MolID, EncodedMol in EncodedMolsInfo:
1134 if EncodeMol is not None:
1135 Mol = MolFromBase64EncodedMolString(EncodedMol)
1136
1137 """
1138 for MolIndex, Mol in enumerate(Mols):
1139 yield (
1140 [MolIDs[MolIndex], None]
1141 if Mol is None
1142 else [MolIDs[MolIndex], MolToBase64EncodedMolString(Mol, PropertyPickleFlags)]
1143 )
1144
1145
1146 def GenerateBase64EncodedMolStringWithConfIDs(
1147 Mol, MolIndex, ConfIDs, PropertyPickleFlags=Chem.PropertyPickleOptions.AllProps
1148 ):
1149 """Setup an iterator generating base64 encoded molecule string for a
1150 molecule. The iterator returns a list containing a molecule index, an encoded
1151 molecule string, and conf ID.
1152
1153 The molecules are pickled using RDKit Mol.ToBinary() function
1154 before their encoding.
1155
1156 Arguments:
1157 Mol (object): RDKit molecule object.
1158 MolIndex (int): Molecule index.
1159 ConfIDs (list): Conformer IDs.
1160 PropertyFlags: RDKit property pickle options.
1161
1162 Returns:
1163 object : Base64 endcoded molecules iterator. The iterator returns a
1164 list containing a molecule index, an encoded molecule string, and
1165 conf ID.
1166
1167 Notes:
1168 The following property pickle flags are currently available in RDKit:
1169
1170 Chem.PropertyPickleOptions.NoProps
1171 Chem.PropertyPickleOptions.MolProps
1172 Chem.PropertyPickleOptions.AtomProps
1173 Chem.PropertyPickleOptions.BondProps
1174 Chem.PropertyPickleOptions.PrivateProps
1175 Chem.PropertyPickleOptions.AllProps
1176
1177 Examples:
1178
1179 EncodedMolsInfo = GenerateBase64EncodedMolStringWithConfIDs(Mol, MolIndex, ConfIDs)
1180 for MolIndex, EncodedMol, ConfID in EncodedMolsInfo:
1181 if EncodeMol is not None:
1182 Mol = MolFromBase64EncodedMolString(EncodedMol)
1183
1184 """
1185 for ConfID in ConfIDs:
1186 yield (
1187 [MolIndex, None, ConfID]
1188 if Mol is None
1189 else [MolIndex, MolToBase64EncodedMolString(Mol, PropertyPickleFlags), ConfID]
1190 )
1191
1192
1193 def AreAtomMapNumbersPresentInMol(Mol):
1194 """Check for the presence of atom map numbers in a molecue.
1195
1196 Arguments:
1197 Mol (object): RDKit molecule object.
1198
1199 Returns:
1200 bool : True - Atom map numbers present; Otherwise, false.
1201
1202 """
1203
1204 return False if _GetAtomMapIndices(Mol) is None else True
1205
1206
1207 def ClearAtomMapNumbers(Mol, AllowImplicitValence=True, ClearRadicalElectrons=True):
1208 """Check and clear atom map numbers in a molecule. In addition, allow implicit
1209 valence and clear radical electrons for atoms with associated map numbers.
1210
1211 For example, the following atomic properties are assigned by RDKit to atom
1212 map number 1 in a molecule corresponding to SMILES C[C:1](C)C:
1213
1214 NoImplicit: True; ImplicitValence: 0; ExplicitValence: 3; NumExplicitHs: 0;
1215 NumImplicitHs: 0; NumRadicalElectrons: 1
1216
1217 This function clears atoms map numbers in the molecule leading to SMILES
1218 CC(C)C, along with optionally updating atomic properties as shown below:
1219
1220 NoImplicit: False; ImplicitValence: 1; ExplicitValence: 3; NumExplicitHs: 0;
1221 NumImplicitHs: 1; NumRadicalElectrons: 0
1222
1223 Arguments:
1224 Mol (object): RDKit molecule object.
1225
1226 Returns:
1227 Mol (object): RDKit molecule object.
1228
1229 """
1230
1231 AtomMapIndices = GetAtomMapIndices(Mol)
1232
1233 if AtomMapIndices is None:
1234 return Mol
1235
1236 for AtomMapIndex in AtomMapIndices:
1237 Atom = Mol.GetAtomWithIdx(AtomMapIndex)
1238
1239 # Clear map number property 'molAtomMapNumber'...
1240 Atom.SetAtomMapNum(0)
1241
1242 # Allow implit valence...
1243 if AllowImplicitValence:
1244 Atom.SetNoImplicit(False)
1245
1246 # Set number of electrons to 0...
1247 if ClearRadicalElectrons:
1248 Atom.SetNumRadicalElectrons(0)
1249
1250 Atom.UpdatePropertyCache()
1251
1252 Mol.UpdatePropertyCache()
1253
1254 return Mol
1255
1256
1257 def GetAtomMapIndices(Mol):
1258 """Get a list of available atom indices corresponding to atom map numbers
1259 present in a SMILES/SMARTS pattern used for creating a molecule. The list of
1260 atom indices is sorted in ascending order by atom map numbers.
1261
1262 Arguments:
1263 Mol (object): RDKit molecule object.
1264
1265 Returns:
1266 list : List of atom indices sorted in the ascending order of atom map
1267 numbers or None.
1268
1269 """
1270
1271 return _GetAtomMapIndices(Mol)
1272
1273
1274 def GetAtomMapIndicesAndMapNumbers(Mol):
1275 """Get lists of available atom indices and atom map numbers present in a
1276 SMILES/SMARTS pattern used for creating a molecule. Both lists are sorted
1277 in ascending order by atom map numbers.
1278
1279 Arguments:
1280 Mol (object): RDKit molecule object.
1281
1282 Returns:
1283 list : List of atom indices sorted in the ascending order of atom map
1284 numbers or None.
1285 list : List of atom map numbers sorted in the ascending order or None.
1286
1287 """
1288
1289 return _GetAtomMapIndicesAndMapNumbers(Mol)
1290
1291
1292 def MolFromSubstructureMatch(Mol, PatternMol, AtomIndices, FilterByAtomMapNums=False):
1293 """Generate a RDKit molecule object for a list of matched atom indices
1294 present in a pattern molecule. The list of atom indices correspond to a
1295 list retrieved by RDKit function GetSubstructureMatches using SMILES/SMARTS
1296 pattern. The atom indices are optionally filtered by mapping atom numbers
1297 to appropriate atom indices during the generation of the molecule.
1298 For example: [O:1]=[S:2](=[O])[C:3][C:4].
1299
1300 Arguments:
1301 Mol (object): RDKit molecule object.
1302 PatternMol (object): RDKit molecule object for a SMILES/SMARTS pattern.
1303 AtomIndices (list): Atom indices.
1304 FilterByAtomMapNums (bool): Filter matches by atom map numbers.
1305
1306 Returns:
1307 object : RDKit molecule object or None.
1308
1309 """
1310
1311 AtomMapIndices = _GetAtomMapIndices(PatternMol) if FilterByAtomMapNums else None
1312
1313 return _MolFromSubstructureMatch(Mol, PatternMol, AtomIndices, AtomMapIndices)
1314
1315
1316 def MolsFromSubstructureMatches(Mol, PatternMol, AtomIndicesList, FilterByAtomMapNums=False):
1317 """Generate a list of RDKit molecule objects for a list containing lists of
1318 matched atom indices present in a pattern molecule. The list of atom indices
1319 correspond to a list retrieved by RDKit function GetSubstructureMatches using
1320 SMILES/SMARTS pattern. The atom indices are optionally filtered by mapping
1321 atom numbers to appropriate atom indices during the generation of the
1322 molecule. For example: [O:1]=[S:2](=[O])[C:3][C:4].
1323
1324 Arguments:
1325 Mol (object): RDKit molecule object.
1326 PatternMol (object): RDKit molecule object for a SMILES/SMARTS pattern.
1327 AtomIndicesList (list): A list of lists containing atom indices.
1328 FilterByAtomMapNums (bool): Filter matches by atom map numbers.
1329
1330 Returns:
1331 list : A list of lists containg RDKit molecule objects or None.
1332
1333 """
1334
1335 AtomMapIndices = _GetAtomMapIndices(PatternMol) if FilterByAtomMapNums else None
1336
1337 Mols = []
1338 for AtomIndices in AtomIndicesList:
1339 Mols.append(_MolFromSubstructureMatch(Mol, PatternMol, AtomIndices, AtomMapIndices))
1340
1341 return Mols if len(Mols) else None
1342
1343
1344 def FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices):
1345 """Filter a list of matched atom indices by map atom numbers present in a
1346 pattern molecule. The list of atom indices correspond to a list retrieved by
1347 RDKit function GetSubstructureMatches using SMILES/SMARTS pattern. The
1348 atom map numbers are mapped to appropriate atom indices during the generation
1349 of molecules. For example: [O:1]=[S:2](=[O])[C:3][C:4].
1350
1351 Arguments:
1352 Mol (object): RDKit molecule object.
1353 PatternMol (object): RDKit molecule object for a SMILES/SMARTS pattern.
1354 AtomIndices (list): Atom indices.
1355
1356 Returns:
1357 list : A list of filtered atom indices.
1358
1359 """
1360 AtomMapIndices = _GetAtomMapIndices(PatternMol)
1361
1362 return _FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices, AtomMapIndices)
1363
1364
1365 def FilterSubstructureMatchesByAtomMapNumbers(Mol, PatternMol, AtomIndicesList):
1366 """Filter a list of lists containing matched atom indices by map atom numbers
1367 present in a pattern molecule. The list of atom indices correspond to a list retrieved by
1368 RDKit function GetSubstructureMatches using SMILES/SMARTS pattern. The
1369 atom map numbers are mapped to appropriate atom indices during the generation
1370 of molecules. For example: [O:1]=[S:2](=[O])[C:3][C:4].
1371
1372 Arguments:
1373 Mol (object): RDKit molecule object.
1374 PatternMol (object): RDKit molecule object for a SMILES/SMARTS pattern.
1375 AtomIndicesList (list): A list of lists containing atom indices.
1376
1377 Returns:
1378 list : A list of lists containing filtered atom indices.
1379
1380 """
1381 AtomMapIndices = _GetAtomMapIndices(PatternMol)
1382
1383 MatchedAtomIndicesList = []
1384 for AtomIndices in AtomIndicesList:
1385 MatchedAtomIndicesList.append(
1386 _FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices, AtomMapIndices)
1387 )
1388
1389 return MatchedAtomIndicesList
1390
1391
1392 def _MolFromSubstructureMatch(Mol, PatternMol, AtomIndices, AtomMapIndices):
1393 """Generate a RDKit molecule object for a list of matched atom indices and available
1394 atom map indices.
1395 """
1396
1397 if AtomMapIndices is not None:
1398 MatchedAtomIndices = [AtomIndices[Index] for Index in AtomMapIndices]
1399 else:
1400 MatchedAtomIndices = list(AtomIndices)
1401
1402 return _GetMolFromAtomIndices(Mol, MatchedAtomIndices)
1403
1404
1405 def _GetAtomMapIndices(Mol):
1406 """Get a list of available atom indices corresponding to sorted atom map
1407 numbers present in a SMILES/SMARTS pattern used for creating a molecule.
1408 """
1409
1410 AtomMapIndices, AtomMapNumbers = _GetAtomMapIndicesAndMapNumbers(Mol)
1411
1412 return AtomMapIndices
1413
1414
1415 def _GetAtomMapIndicesAndMapNumbers(Mol):
1416 """Get a list of available atom indices and atom map numbers present
1417 in a SMILES/SMARTS pattern used for creating a molecule. Both lists
1418 are sorted in ascending order by atom map numbers.
1419 """
1420
1421 # Setup a atom map number to atom indices map..
1422 AtomMapNumToIndices = {}
1423 for Atom in Mol.GetAtoms():
1424 AtomMapNum = Atom.GetAtomMapNum()
1425
1426 if AtomMapNum:
1427 AtomMapNumToIndices[AtomMapNum] = Atom.GetIdx()
1428
1429 # Setup atom indices corresponding to sorted atom map numbers...
1430 AtomMapIndices = None
1431 AtomMapNumbers = None
1432 if len(AtomMapNumToIndices):
1433 AtomMapNumbers = sorted(AtomMapNumToIndices)
1434 AtomMapIndices = [AtomMapNumToIndices[AtomMapNum] for AtomMapNum in AtomMapNumbers]
1435
1436 return (AtomMapIndices, AtomMapNumbers)
1437
1438
1439 def _FilterSubstructureMatchByAtomMapNumbers(Mol, PatternMol, AtomIndices, AtomMapIndices):
1440 """Filter substructure match atom indices by atom map indices corresponding to
1441 atom map numbers.
1442 """
1443
1444 if AtomMapIndices is None:
1445 return list(AtomIndices)
1446
1447 return [AtomIndices[Index] for Index in AtomMapIndices]
1448
1449
1450 def _GetMolFromAtomIndices(Mol, AtomIndices):
1451 """Generate a RDKit molecule object from atom indices returned by
1452 substructure search.
1453 """
1454
1455 BondIndices = []
1456 for AtomIndex in AtomIndices:
1457 Atom = Mol.GetAtomWithIdx(AtomIndex)
1458
1459 for AtomNbr in Atom.GetNeighbors():
1460 AtomNbrIndex = AtomNbr.GetIdx()
1461 if AtomNbrIndex not in AtomIndices:
1462 continue
1463
1464 BondIndex = Mol.GetBondBetweenAtoms(AtomIndex, AtomNbrIndex).GetIdx()
1465 if BondIndex in BondIndices:
1466 continue
1467
1468 BondIndices.append(BondIndex)
1469
1470 MatchedMol = Chem.PathToSubmol(Mol, BondIndices) if len(BondIndices) else None
1471
1472 return MatchedMol
1473
1474
1475 def ConstrainAndEmbed(
1476 mol,
1477 core,
1478 coreMatchesMol=None,
1479 useTethers=True,
1480 coreConfId=-1,
1481 randomseed=2342,
1482 getForceField=AllChem.UFFGetMoleculeForceField,
1483 **kwargs,
1484 ):
1485 """
1486 The function is a local copy of RDKit fucntion AllChem.ConstrainedEmbed().
1487 It has been enhanced to support an explicit list of core matches corresponding
1488 to the matched atom indices in the molecule. The number of matched atom indices
1489 must be equal to the number of atoms in core molecule.
1490
1491 Arguments:
1492 mol (object): RDKit molecule object to embed.
1493 core (object): RDKit molecule to use as a source of constraints.
1494 coreMatchesMol (list): A list matches atom indices in mol.
1495 useTethers: (bool) if True, the final conformation will be optimized
1496 subject to a series of extra forces that pull the matching atoms to
1497 the positions of the core atoms. Otherwise simple distance
1498 constraints based on the core atoms will be used in the
1499 optimization.
1500 coreConfId (int): ID of the core conformation to use.
1501 randomSeed (int): Seed for the random number generator
1502
1503 Returns:
1504 mol (object): RDKit molecule object.
1505
1506 """
1507 if coreMatchesMol is None:
1508 match = mol.GetSubstructMatch(core)
1509 if not match:
1510 raise ValueError("Molecule doesn't match the core.")
1511 else:
1512 if core.GetNumAtoms() != len(coreMatchesMol):
1513 raise ValueError(
1514 "Number of atoms, %s, in core molecule must match number of atom indices, %s, specified in the list coreMatchesMol."
1515 % (core.GetNumAtoms(), len(coreMatchesMol))
1516 )
1517 # Check specified matched atom indices in coreMatchesMol and use the match atom
1518 # indices returned by GetSubstructMatches() for embedding...
1519 coreMatch = None
1520 matches = mol.GetSubstructMatches(core)
1521 for match in matches:
1522 if len(match) != len(coreMatchesMol):
1523 continue
1524 matchFound = True
1525 for atomIndex in match:
1526 if atomIndex not in coreMatchesMol:
1527 matchFound = False
1528 break
1529 if matchFound:
1530 coreMatch = match
1531 break
1532 if coreMatch is None:
1533 raise ValueError("Molecule doesn't match the atom indices specified in the list coreMatchesMol.")
1534 match = coreMatch
1535
1536 coordMap = {}
1537 coreConf = core.GetConformer(coreConfId)
1538 for i, idxI in enumerate(match):
1539 corePtI = coreConf.GetAtomPosition(i)
1540 coordMap[idxI] = corePtI
1541
1542 ci = AllChem.EmbedMolecule(mol, coordMap=coordMap, randomSeed=randomseed, **kwargs)
1543 if ci < 0:
1544 raise ValueError("Could not embed molecule.")
1545
1546 algMap = [(j, i) for i, j in enumerate(match)]
1547
1548 if not useTethers:
1549 # clean up the conformation
1550 ff = getForceField(mol, confId=0)
1551 for i, idxI in enumerate(match):
1552 for j in range(i + 1, len(match)):
1553 idxJ = match[j]
1554 d = coordMap[idxI].Distance(coordMap[idxJ])
1555 ff.AddDistanceConstraint(idxI, idxJ, d, d, 100.0)
1556 ff.Initialize()
1557 n = 4
1558 more = ff.Minimize()
1559 while more and n:
1560 more = ff.Minimize()
1561 n -= 1
1562 # rotate the embedded conformation onto the core:
1563 rms = AllChem.AlignMol(mol, core, atomMap=algMap)
1564 else:
1565 # rotate the embedded conformation onto the core:
1566 rms = AllChem.AlignMol(mol, core, atomMap=algMap)
1567 ff = getForceField(mol, confId=0)
1568 conf = core.GetConformer()
1569 for i in range(core.GetNumAtoms()):
1570 p = conf.GetAtomPosition(i)
1571 pIdx = ff.AddExtraPoint(p.x, p.y, p.z, fixed=True) - 1
1572 ff.AddDistanceConstraint(pIdx, match[i], 0, 0, 100.0)
1573 ff.Initialize()
1574 n = 4
1575 more = ff.Minimize(energyTol=1e-4, forceTol=1e-3)
1576 while more and n:
1577 more = ff.Minimize(energyTol=1e-4, forceTol=1e-3)
1578 n -= 1
1579 # realign
1580 rms = AllChem.AlignMol(mol, core, atomMap=algMap)
1581
1582 mol.SetProp("EmbedRMS", str(rms))
1583 return mol
1584
1585
1586 def ReadAndValidateMolecules(FileName, **KeyWordArgs):
1587 """Read molecules from an input file, validate all molecule objects, and return
1588 a list of valid molecule objects along with the count of valid and non-valid molecule
1589 objects.
1590
1591 Arguments:
1592 FileName (str): Name of a file with complete path.
1593 **KeyWordArgs (dictionary) : Parameter name and value pairs for reading and
1594 processing molecules.
1595
1596 Returns:
1597 list : List of valid RDKit molecule objects.
1598 int : Number of total molecules in input file.
1599 int : Number of valid molecules in input file.
1600
1601 Notes:
1602 The file extension is used to determine type of the file and set up an appropriate
1603 file reader.
1604
1605 """
1606
1607 AllowEmptyMols = True
1608 if "AllowEmptyMols" in KeyWordArgs:
1609 AllowEmptyMols = KeyWordArgs["AllowEmptyMols"]
1610
1611 Mols = ReadMolecules(FileName, **KeyWordArgs)
1612
1613 if AllowEmptyMols:
1614 ValidMols = [Mol for Mol in Mols if Mol is not None]
1615 else:
1616 ValidMols = []
1617 MolCount = 0
1618 for Mol in Mols:
1619 MolCount += 1
1620 if Mol is None:
1621 continue
1622
1623 if IsMolEmpty(Mol):
1624 MolName = GetMolName(Mol, MolCount)
1625 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
1626 continue
1627
1628 ValidMols.append(Mol)
1629
1630 MolCount = len(Mols)
1631 ValidMolCount = len(ValidMols)
1632
1633 return (ValidMols, MolCount, ValidMolCount)
1634
1635
1636 def ReadMolecules(FileName, **KeyWordArgs):
1637 """Read molecules from an input file without performing any validation
1638 and creation of molecule objects.
1639
1640 Arguments:
1641 FileName (str): Name of a file with complete path.
1642 **KeyWordArgs (dictionary) : Parameter name and value pairs for reading and
1643 processing molecules.
1644
1645 Returns:
1646 list : List of RDKit molecule objects.
1647
1648 Notes:
1649 The file extension is used to determine type of the file and set up an appropriate
1650 file reader.
1651
1652 """
1653
1654 # Set default values for possible arguments...
1655 ReaderArgs = {
1656 "Sanitize": True,
1657 "RemoveHydrogens": True,
1658 "StrictParsing": True,
1659 "SMILESDelimiter": " ",
1660 "SMILESColumn": 1,
1661 "SMILESNameColumn": 2,
1662 "SMILESTitleLine": True,
1663 }
1664
1665 # Set specified values for possible arguments...
1666 for Arg in ReaderArgs:
1667 if Arg in KeyWordArgs:
1668 ReaderArgs[Arg] = KeyWordArgs[Arg]
1669
1670 # Modify specific valeus for SMILES...
1671 if MiscUtil.CheckFileExt(FileName, "smi csv tsv txt"):
1672 Args = ["Sanitize", "SMILESTitleLine"]
1673 for Arg in Args:
1674 if ReaderArgs[Arg] is True:
1675 ReaderArgs[Arg] = 1
1676 else:
1677 ReaderArgs[Arg] = 0
1678
1679 Mols = []
1680 if MiscUtil.CheckFileExt(FileName, "sdf sd"):
1681 return ReadMoleculesFromSDFile(
1682 FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"], ReaderArgs["StrictParsing"]
1683 )
1684 elif MiscUtil.CheckFileExt(FileName, "mol"):
1685 return ReadMoleculesFromMolFile(
1686 FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"], ReaderArgs["StrictParsing"]
1687 )
1688 elif MiscUtil.CheckFileExt(FileName, "mol2"):
1689 return ReadMoleculesFromMol2File(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"])
1690 elif MiscUtil.CheckFileExt(FileName, "pdb"):
1691 return ReadMoleculesFromPDBFile(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"])
1692 elif MiscUtil.CheckFileExt(FileName, "smi txt csv tsv"):
1693 SMILESColumnIndex = ReaderArgs["SMILESColumn"] - 1
1694 SMILESNameColumnIndex = ReaderArgs["SMILESNameColumn"] - 1
1695 return ReadMoleculesFromSMILESFile(
1696 FileName,
1697 ReaderArgs["SMILESDelimiter"],
1698 SMILESColumnIndex,
1699 SMILESNameColumnIndex,
1700 ReaderArgs["SMILESTitleLine"],
1701 ReaderArgs["Sanitize"],
1702 )
1703 else:
1704 MiscUtil.PrintWarning("RDKitUtil.ReadMolecules: Non supported file type: %s" % FileName)
1705
1706 return Mols
1707
1708
1709 def ReadMoleculesFromSDFile(FileName, Sanitize=True, RemoveHydrogens=True, StrictParsing=True):
1710 """Read molecules from a SD file.
1711
1712 Arguments:
1713 FileName (str): Name of a file with complete path.
1714 Sanitize (bool): Sanitize molecules.
1715 RemoveHydrogens (bool): Remove hydrogens from molecules.
1716 StrictParsing (bool): Perform strict parsing.
1717
1718 Returns:
1719 list : List of RDKit molecule objects.
1720
1721 """
1722 return Chem.SDMolSupplier(FileName, sanitize=Sanitize, removeHs=RemoveHydrogens, strictParsing=StrictParsing)
1723
1724
1725 def ReadMoleculesFromMolFile(FileName, Sanitize=True, RemoveHydrogens=True, StrictParsing=True):
1726 """Read molecule from a MDL Mol file.
1727
1728 Arguments:
1729 FileName (str): Name of a file with complete path.
1730 Sanitize (bool): Sanitize molecules.
1731 RemoveHydrogens (bool): Remove hydrogens from molecules.
1732 StrictParsing (bool): Perform strict parsing.
1733
1734 Returns:
1735 list : List of RDKit molecule objects.
1736
1737 """
1738
1739 Mols = []
1740 Mols.append(Chem.MolFromMolFile(FileName, sanitize=Sanitize, removeHs=RemoveHydrogens, strictParsing=StrictParsing))
1741 return Mols
1742
1743
1744 def ReadMoleculesFromMol2File(FileName, Sanitize=True, RemoveHydrogens=True):
1745 """Read molecules from a Tripos Mol2 file. The first call to the function
1746 creates and returns a generator object using Python yield statement. The
1747 molecules are created during the subsequent iteration by the generator object.
1748
1749 Arguments:
1750 FileName (str): Name of a file with complete path.
1751 Sanitize (bool): Sanitize molecules.
1752 RemoveHydrogens (bool): Remove hydrogens from molecules.
1753
1754 Returns:
1755 list : A Python generator object for iterating over the molecules.
1756
1757 """
1758
1759 return _Mol2MolSupplier(FileName, Sanitize, RemoveHydrogens)
1760
1761
1762 def _Mol2MolSupplier(FileName, Sanitize=True, RemoveHydrogens=True):
1763 """Read molecules from a Tripos Mol2 file."""
1764
1765 fh = open(FileName, "r")
1766
1767 FirstMol = True
1768 ProcessingMol = False
1769
1770 for Line in fh:
1771 if re.match("^#", Line, re.I):
1772 continue
1773
1774 if re.match("^@<TRIPOS>MOLECULE", Line, re.I):
1775 ProcessingMol = True
1776
1777 if FirstMol:
1778 FirstMol = False
1779
1780 MolLines = []
1781 MolLines.append(Line)
1782 continue
1783
1784 # Process lines for existing molecule...
1785 MolBlock = "".join(MolLines)
1786
1787 Mol = Chem.MolFromMol2Block(MolBlock, sanitize=Sanitize, removeHs=RemoveHydrogens)
1788 yield Mol
1789
1790 # Track lines for next molecule...
1791 MolLines = []
1792 MolLines.append(Line)
1793 continue
1794
1795 if not ProcessingMol:
1796 continue
1797
1798 MolLines.append(Line)
1799
1800 fh.close
1801
1802 # Process last molecule...
1803 if len(MolLines):
1804 MolBlock = "".join(MolLines)
1805 Mol = Chem.MolFromMol2Block(MolBlock, sanitize=Sanitize, removeHs=RemoveHydrogens)
1806 yield Mol
1807
1808
1809 def ReadMoleculesFromPDBFile(FileName, Sanitize=True, RemoveHydrogens=True):
1810 """Read molecule from a PDB file.
1811
1812 Arguments:
1813 FileName (str): Name of a file with complete path.
1814 Sanitize (bool): Sanitize molecules.
1815 RemoveHydrogens (bool): Remove hydrogens from molecules.
1816
1817 Returns:
1818 list : List of RDKit molecule objects.
1819
1820 """
1821
1822 Mols = []
1823 Mols.append(Chem.MolFromPDBFile(FileName, sanitize=Sanitize, removeHs=RemoveHydrogens))
1824 return Mols
1825
1826
1827 def ReadMoleculesFromSMILESFile(
1828 FileName, SMILESDelimiter=" ", SMILESColIndex=0, SMILESNameColIndex=1, SMILESTitleLine=1, Sanitize=1
1829 ):
1830 """Read molecules from a SMILES file.
1831
1832 Arguments:
1833 SMILESDelimiter (str): Delimiter for parsing SMILES line
1834 SMILESColIndex (int): Column index containing SMILES string.
1835 SMILESNameColIndex (int): Column index containing molecule name.
1836 SMILESTitleLine (int): Flag to indicate presence of title line.
1837 Sanitize (int): Sanitize molecules.
1838
1839 Returns:
1840 list : List of RDKit molecule objects.
1841
1842 """
1843
1844 return Chem.SmilesMolSupplier(
1845 FileName,
1846 delimiter=SMILESDelimiter,
1847 smilesColumn=SMILESColIndex,
1848 nameColumn=SMILESNameColIndex,
1849 titleLine=SMILESTitleLine,
1850 sanitize=Sanitize,
1851 )
1852
1853
1854 def MoleculesWriter(FileName, **KeyWordArgs):
1855 """Set up a molecule writer.
1856
1857 Arguments:
1858 FileName (str): Name of a file with complete path.
1859 **KeyWordArgs (dictionary) : Parameter name and value pairs for writing and
1860 processing molecules.
1861
1862 Returns:
1863 RDKit object : Molecule writer.
1864
1865 Notes:
1866 The file extension is used to determine type of the file and set up an appropriate
1867 file writer.
1868
1869 """
1870
1871 # Set default values for possible arguments...
1872 WriterArgs = {
1873 "Compute2DCoords": False,
1874 "Kekulize": True,
1875 "ForceV3000": False,
1876 "SMILESKekulize": False,
1877 "SMILESDelimiter": " ",
1878 "SMILESIsomeric": True,
1879 "SMILESTitleLine": True,
1880 "SMILESMolName": True,
1881 }
1882
1883 # Set specified values for possible arguments...
1884 for Arg in WriterArgs:
1885 if Arg in KeyWordArgs:
1886 WriterArgs[Arg] = KeyWordArgs[Arg]
1887
1888 Writer = None
1889 if MiscUtil.CheckFileExt(FileName, "sdf sd"):
1890 Writer = Chem.SDWriter(FileName)
1891 Writer.SetKekulize(WriterArgs["Kekulize"])
1892 Writer.SetForceV3000(WriterArgs["ForceV3000"])
1893 elif MiscUtil.CheckFileExt(FileName, "pdb"):
1894 Writer = Chem.PDBWriter(FileName)
1895 elif MiscUtil.CheckFileExt(FileName, "smi"):
1896 # Text for the name column in the title line. Blank indicates not to include name column
1897 # in the output file...
1898 NameHeader = "Name" if WriterArgs["SMILESMolName"] else ""
1899 Writer = Chem.SmilesWriter(
1900 FileName,
1901 delimiter=WriterArgs["SMILESDelimiter"],
1902 nameHeader=NameHeader,
1903 includeHeader=WriterArgs["SMILESTitleLine"],
1904 isomericSmiles=WriterArgs["SMILESIsomeric"],
1905 kekuleSmiles=WriterArgs["SMILESKekulize"],
1906 )
1907 else:
1908 MiscUtil.PrintWarning("RDKitUtil.WriteMolecules: Non supported file type: %s" % FileName)
1909
1910 return Writer
1911
1912
1913 def WriteMolecules(FileName, Mols, **KeyWordArgs):
1914 """Write molecules to an output file.
1915
1916 Arguments:
1917 FileName (str): Name of a file with complete path.
1918 Mols (list): List of RDKit molecule objects.
1919 **KeyWordArgs (dictionary) : Parameter name and value pairs for writing and
1920 processing molecules.
1921
1922 Returns:
1923 int : Number of total molecules.
1924 int : Number of processed molecules written to output file.
1925
1926 Notes:
1927 The file extension is used to determine type of the file and set up an appropriate
1928 file writer.
1929
1930 """
1931
1932 Compute2DCoords = False
1933 if "Compute2DCoords" in KeyWordArgs:
1934 Compute2DCoords = KeyWordArgs["Compute2DCoords"]
1935
1936 SetSMILESMolProps = KeyWordArgs["SetSMILESMolProps"] if "SetSMILESMolProps" in KeyWordArgs else False
1937
1938 MolCount = 0
1939 ProcessedMolCount = 0
1940
1941 Writer = MoleculesWriter(FileName, **KeyWordArgs)
1942
1943 if Writer is None:
1944 return (MolCount, ProcessedMolCount)
1945
1946 FirstMol = True
1947 for Mol in Mols:
1948 MolCount += 1
1949 if Mol is None:
1950 continue
1951
1952 if FirstMol:
1953 FirstMol = False
1954 if SetSMILESMolProps:
1955 SetWriterMolProps(Writer, Mol)
1956
1957 ProcessedMolCount += 1
1958 if Compute2DCoords:
1959 AllChem.Compute2DCoords(Mol)
1960
1961 Writer.write(Mol)
1962
1963 Writer.close()
1964
1965 return (MolCount, ProcessedMolCount)
1966
1967
1968 def SetWriterMolProps(Writer, Mol):
1969 """Setup molecule properties for a writer to output.
1970
1971 Arguments:
1972 Writer (object): RDKit writer object.
1973 Mol (object): RDKit molecule object.
1974
1975 Returns:
1976 object : Writer object.
1977
1978 """
1979 PropNames = list(Mol.GetPropNames())
1980 if len(PropNames):
1981 Writer.SetProps(PropNames)
1982
1983 return Writer
1984
1985
1986 def SetupHTMLForTorsionScanViewer(
1987 MolName, TorsionMols, TorsionEnergies, RelativeTorsionEnergies, TorsionAngles, Units, PlotHeight, TitleLine
1988 ):
1989 """Setup HTML for torsion scan viewer for writing to a text file. The viewer
1990 consists of plot and structure viewer panels. The plot panel contains an
1991 interactive relative energy rendered using Plotly. The structure of the
1992 molecules, corresponding to the scanned torsion angles, are shown in the
1993 structure panel employing 3DMol.js. In addition, an animation button is
1994 available to automatically scan the torsion angles and view the results.
1995
1996 Arguments:
1997 MolName (str): Molecule name.
1998 TorsionEnergies (list): List of energies for scanned torsion angles.
1999 RelativeTorsionEnergies (list): List of relative energies for scanned
2000 torsion angles.
2001 TorsionAngles (list): List of scanned torsion angles.
2002 Units (str): Energy units.
2003 PlotHeight (int): Energy plot height in pixels.
2004 TitleLine (str): Title line torsion scan viewer.
2005
2006 Returns:
2007 str: HTML text for torsion scan viewer.
2008
2009 """
2010 TorsionDataFrames = _SetupTorsionScanViewerDataFrames(
2011 MolName, TorsionMols, TorsionEnergies, RelativeTorsionEnergies, TorsionAngles
2012 )
2013
2014 # Setup PlotHeight and strucure viewer minimum height...
2015 StrViewerMinHeight = PlotHeight - 90
2016 if StrViewerMinHeight <= 0:
2017 StrViewerMinHeight = PlotHeight
2018
2019 PlotHeight = "%spx" % PlotHeight
2020 StrViewerMinHeight = "%spx" % StrViewerMinHeight
2021
2022 # Setup footer line...
2023 FooterLine = "Visualization: Plotly · 3Dmol.js; Calculations: RDKit · MayaChemTools"
2024
2025 # Setup a JSON data string...
2026 TorsionDataFramesJSONString = json.dumps(TorsionDataFrames, separators=(",", ":"))
2027
2028 # Setup viewer HTML template string...
2029 ViewerTemplateHTML = _SetupHTMLTorsionScanViewerTemplate()
2030
2031 # Update viewer HTML template string...
2032 ViewerHTML = ViewerTemplateHTML.replace("__DATA__", TorsionDataFramesJSONString)
2033
2034 ViewerHTML = ViewerHTML.replace("__TITLE__", TitleLine).replace("__FOOTER__", FooterLine)
2035 ViewerHTML = ViewerHTML.replace("__UNITS__", Units)
2036 ViewerHTML = ViewerHTML.replace("__PLOT_HEIGHT__", PlotHeight).replace("__VIEWER_MIN_HEIGHT__", StrViewerMinHeight)
2037
2038 return ViewerHTML
2039
2040
2041 def _SetupTorsionScanViewerDataFrames(MolName, TorsionMols, TorsionEnergies, RelativeTorsionEnergies, TorsionAngles):
2042 """Setup torsion scan viewer data frames."""
2043
2044 TorsionDataFrames = []
2045
2046 for Index, TorsionMol in enumerate(TorsionMols):
2047 TorsionAngle = "%s" % TorsionAngles[Index]
2048
2049 TorsionMolName = "%s_Deg%s" % (MolName, TorsionAngle)
2050 TorsionMol.SetProp("_Name", TorsionMolName)
2051
2052 # Setup molblock with standardized new lines...
2053 TorsionMolBlock = Chem.MolToMolBlock(TorsionMol)
2054 TorsionMolBlock = re.sub(r"(\r\n|\r)", r"\n", TorsionMolBlock)
2055
2056 TorsionEnergy = "%.2f" % TorsionEnergies[Index]
2057 RelativeTorsionEnergy = "%.2f" % RelativeTorsionEnergies[Index]
2058
2059 TorsionDataFrames.append(
2060 {
2061 "title": TorsionMolName,
2062 "angle": float(TorsionAngle),
2063 "energy": float(TorsionEnergy),
2064 "relEnergy": float(RelativeTorsionEnergy),
2065 "mol": TorsionMolBlock,
2066 }
2067 )
2068
2069 return TorsionDataFrames
2070
2071
2072 def _SetupHTMLTorsionScanViewerTemplate():
2073 """Setup torsion scan viewer template."""
2074
2075 ViewerTemplate = r"""<!DOCTYPE html>
2076 <html lang="en">
2077 <head>
2078 <meta charset="utf-8">
2079 <meta name="viewport" content="width=device-width, initial-scale=1">
2080 <title>Torsion Profile Viewer</title>
2081 <script src="https://cdn.plot.ly/plotly-2.35.2.min.js" charset="utf-8"></script>
2082 <script src="https://3Dmol.org/build/3Dmol-min.js"></script>
2083 <style>
2084 :root { --bg:#f7f7fa; --panel:#fff; --line:#4c72b0; --accent:#dd6b3f; --ink:#262626; }
2085 * { box-sizing: border-box; }
2086 body { margin:0; font-family:sans-serif,Helvetica,Arial;
2087 background:var(--bg); color:var(--ink); }
2088 header { padding:14px 22px 6px; }
2089 header h1 { margin:0; font-size:18px; font-weight:700; }
2090 header p { margin:4px 0 0; font-size:13px; color:#666; }
2091 .wrap { display:flex; flex-wrap:wrap; gap:16px; padding:14px 22px 22px; align-items:stretch; }
2092
2093 .panel { background:var(--panel); border:1px solid #e3e3ea; border-radius:10px;
2094 box-shadow:0 1px 3px rgba(0,0,0,.05); }
2095 #plotPanel { flex:1 1 460px; min-width:380px; padding:8px; }
2096 #plot { width:100%; height:__PLOT_HEIGHT__; }
2097 #structPanel { flex:1 1 420px; min-width:340px; display:flex; flex-direction:column; padding:12px; }
2098 #viewerHolder { position:relative; flex:1 1 auto; min-height:__VIEWER_MIN_HEIGHT__; border-radius:8px;
2099 overflow:hidden; background:#fdfdfd; border:1px solid #eee; }
2100 #viewer { position:absolute; inset:0; }
2101
2102 .readout { display:flex; gap:18px; flex-wrap:wrap; margin-top:10px; font-size:13px; }
2103 .readout b { font-size:15px; font-weight:700; display:block; }
2104 .readout .lab { color:#888; font-size:12px; letter-spacing:.04em; }
2105
2106 .controls { display:flex; align-items:center; gap:12px; margin-top:12px; }
2107 .controls button { font-size:15px; border:1px solid #ccd; background:#fff; border-radius:7px;
2108 padding:6px 14px; cursor:pointer; }
2109 .controls button:hover { background:#eef; }
2110
2111 #slider { flex:1 1 auto; }
2112
2113 .opts { display:flex; gap:16px; align-items:center; margin-top:6px; font-size:12.5px; color:#555; }
2114 .opts label { cursor:pointer; }
2115
2116 footer { padding:0 22px 18px; font-size:11.5px; color:#999; }
2117 </style>
2118 </head>
2119 <body>
2120 <header>
2121 <h1>Torsional Energy Profile</h1>
2122 <p>__TITLE__ </p>
2123 </header>
2124
2125 <div class="wrap">
2126 <div class="panel" id="plotPanel"><div id="plot"></div></div>
2127
2128 <div class="panel" id="structPanel">
2129 <div id="viewerHolder"><div id="viewer"></div></div>
2130 <div class="readout">
2131 <div><span class="lab">Torsion Angle</span><b><span id="rAngle">–</span>°</b></div>
2132 <div><span class="lab">Relative Energy</span><b><span id="rRelEnergy">–</span></b><span class="lab">__UNITS__</span></div>
2133 <div><span class="lab">Energy</span><b><span id="rAbs">–</span></b><span class="lab">__UNITS__</span></div>
2134 </div>
2135 <div class="controls">
2136 <button id="play">▶ Play</button>
2137 <input type="range" id="slider" min="0" max="0" value="0" step="1">
2138 </div>
2139 <div class="opts">
2140 <label><input type="checkbox" id="optSticks" checked> Sticks</label>
2141 <label><input type="checkbox" id="optH"> Show H</label>
2142 <label><input type="checkbox" id="optSpin"> Spin</label>
2143 <label><input type="checkbox" id="optFreeze"> Lock view while scanning</label>
2144 </div>
2145 </div>
2146 </div>
2147 <footer>__FOOTER__</footer>
2148
2149 <script>
2150 const FRAMES = __DATA__;
2151 const angles = FRAMES.map(f => f.angle);
2152 const relEnergy = FRAMES.map(f => f.relEnergy);
2153
2154 let current = 0;
2155 let playing = false;
2156 let timer = null;
2157
2158 const viewer = $3Dmol.createViewer("viewer", { backgroundColor: "white" });
2159 let model = null;
2160
2161 function renderStructure(idx, keepView) {
2162 const f = FRAMES[idx];
2163 const showH = document.getElementById("optH").checked;
2164
2165 viewer.removeAllModels();
2166
2167 model = viewer.addModel(f.mol, "sdf");
2168 const sticks = document.getElementById("optSticks").checked;
2169 const base = sticks
2170 ? { stick: { radius: 0.20 }, sphere: { scale: 0.0 } }
2171 : { stick: { radius: 0.20 }, sphere: { scale: 0.30 } };
2172 viewer.setStyle({}, base);
2173
2174 if (!showH) viewer.setStyle({ elem: "H" }, { stick: { hidden: true }, sphere: { hidden: true } });
2175
2176 if (!keepView) viewer.zoomTo();
2177
2178 viewer.render();
2179 }
2180
2181 function buildPlot() {
2182 const profile = {
2183 x: angles, y: relEnergy, mode: "lines+markers", type: "scatter",
2184 line: { color: "#4c72b0", width: 2 },
2185 marker: { color: "#4c72b0", size: 6, line: { color: "#fff", width: 1 } },
2186 hovertemplate: "%{x:.0f}°<br>%{y:.2f} __UNITS__<extra></extra>",
2187 name: "RelativeEnergy"
2188 };
2189
2190 const cursor = {
2191 x: [angles[current]], y: [relEnergy[current]], mode: "markers", type: "scatter",
2192 marker: { color: "#dd6b3f", size: 10, line: { color: "#fff", width: 2 }, symbol: "circle" },
2193 hoverinfo: "skip", showlegend: false, name: "selected"
2194 };
2195
2196 const layout = {
2197 margin: { l: 56, r: 16, t: 10, b: 48 },
2198 paper_bgcolor: "#fff", plot_bgcolor: "#eaeaf2",
2199 xaxis: { title: "Torsion Angle (degrees)", gridcolor: "#fff", zeroline: false },
2200 yaxis: { title: "Relative Energy (__UNITS__)", gridcolor: "#fff", zeroline: false },
2201 showlegend: false, hovermode: "closest"
2202 };
2203
2204 Plotly.newPlot("plot", [profile, cursor], layout, { displayModeBar: false, responsive: true });
2205
2206 const plotDiv = document.getElementById("plot");
2207
2208 function nearestByAngle(xVal) {
2209 let bestIdx = 0, bestDiff = Infinity;
2210 for (let i = 0; i < angles.length; i++) {
2211 const a = angles[i];
2212 if (a == null) continue;
2213 const d = Math.abs(a - xVal);
2214 if (d < bestDiff) { bestDiff = d; bestIdx = i; }
2215 }
2216 return bestIdx;
2217 }
2218
2219 plotDiv.on("plotly_click", e => { if (e.points.length) select(e.points[0].pointNumber); });
2220 plotDiv.on("plotly_hover", e => {
2221 if (!playing && e.points.length && e.points[0].curveNumber === 0)
2222 select(e.points[0].pointNumber);
2223 });
2224
2225 // Fallback: clicks landing in the plot area but not on a data point
2226 // still snap to the nearest rotamer by torsion angle.
2227 plotDiv.addEventListener("click", evt => {
2228 const xaxis = plotDiv._fullLayout && plotDiv._fullLayout.xaxis;
2229 if (!xaxis || typeof xaxis.p2d !== "function") return;
2230
2231 const bbox = plotDiv.getBoundingClientRect();
2232 const xPixel = evt.clientX - bbox.left - xaxis._offset;
2233 if (xPixel < 0 || xPixel > xaxis._length) return;
2234
2235 select(nearestByAngle(xaxis.p2d(xPixel)));
2236 });
2237 }
2238
2239 function moveCursor() {
2240 Plotly.restyle("plot", { x: [[angles[current]]], y: [[relEnergy[current]]] }, [1]);
2241 }
2242
2243 function select(idx) {
2244 if (idx == null || idx < 0 || idx >= FRAMES.length) return;
2245
2246 current = idx;
2247 const f = FRAMES[idx];
2248
2249 document.getElementById("rAngle").textContent = f.angle != null ? f.angle.toFixed(0) : "NA";
2250 document.getElementById("rRelEnergy").textContent = f.relEnergy != null ? f.relEnergy.toFixed(2) : "NA";
2251 document.getElementById("rAbs").textContent = f.energy!= null ? f.energy.toFixed(2): "NA";
2252 document.getElementById("slider").value = idx;
2253
2254 moveCursor();
2255
2256 renderStructure(idx, document.getElementById("optFreeze").checked);
2257 }
2258
2259 function play() {
2260 playing = true;
2261 document.getElementById("play").innerHTML = "❙❙ Pause";
2262 document.getElementById("optFreeze").checked = true; // steadier during animation
2263
2264 timer = setInterval(() => select((current + 1) % FRAMES.length), 120);
2265 }
2266
2267 function pause() {
2268 playing = false;
2269 document.getElementById("play").innerHTML = "▶ Play";
2270
2271 clearInterval(timer);
2272 }
2273
2274 window.addEventListener("DOMContentLoaded", () => {
2275 const slider = document.getElementById("slider");
2276 slider.max = FRAMES.length - 1;
2277
2278 slider.addEventListener("input", e => { if (playing) pause(); select(+e.target.value); });
2279
2280 document.getElementById("play").addEventListener("click", () => playing ? pause() : play());
2281 ["optSticks", "optH"].forEach(id =>
2282 document.getElementById(id).addEventListener("change", () => renderStructure(current, true)));
2283
2284 document.getElementById("optSpin").addEventListener("change", e => {
2285 viewer.spin(e.target.checked ? "y" : false);
2286 });
2287
2288 buildPlot();
2289 select(0);
2290 });
2291
2292 </script>
2293 </body>
2294 </html>
2295 """
2296 return ViewerTemplate
