forked from luancarvalhomartins/PyAutoFEP
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmol_util.py
1683 lines (1467 loc) · 87.7 KB
/
mol_util.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#! /usr/bin/env python3
#
# mol_util.py
#
# Copyright 2019 Luan Carvalho Martins <[email protected]>
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA.
#
#
import os
import shutil
import itertools
import rdkit
import rdkit.Chem
import rdkit.Chem.AllChem
import numpy
import os_util
def rwmol_to_obmol(rdkit_rwmol, verbosity=0):
""" Converts a rdkit.RWMol to openbabel.OBMol
:param rdkit.Chem.rdchem.Mol rdkit_rwmol: the ROMol to be converted
:param int verbosity: be verbosity
:rtype: pybel.ob.OBMol
"""
try:
from openbabel import pybel
except ImportError:
import pybel
if verbosity < os_util.verbosity_level.extra_debug:
pybel.ob.obErrorLog.SetOutputLevel(pybel.ob.obError)
else:
os_util.local_print('OpenBabel warning messages are on, expect a lot of output.',
msg_verbosity=os_util.verbosity_level.extra_debug, current_verbosity=verbosity)
if isinstance(rdkit_rwmol, pybel.ob.OBMol):
os_util.local_print('Molecule {} (SMILES={}) is already a pybel.ob.OBMol'
''.format(rdkit_rwmol, pybel.Molecule(rdkit_rwmol).write('smi')),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.warning)
return rdkit_rwmol
if isinstance(rdkit_rwmol, pybel.Molecule):
os_util.local_print('Molecule {} (SMILES={}) is already a a pybel.Molecule, converting to pybel.ob.OBMol only'
''.format(rdkit_rwmol, rdkit.Chem.MolToSmiles(rdkit_rwmol)),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.warning)
return rdkit_rwmol.OBMol
# Set some lookups
_bondorders = {rdkit.Chem.BondType.SINGLE: 1,
rdkit.Chem.rdchem.BondType.UNSPECIFIED: 1,
rdkit.Chem.BondType.DOUBLE: 2,
rdkit.Chem.BondType.TRIPLE: 3,
rdkit.Chem.BondType.AROMATIC: 5}
_bondstereo = {rdkit.Chem.rdchem.BondStereo.STEREONONE: 0,
rdkit.Chem.rdchem.BondStereo.STEREOE: 1,
rdkit.Chem.rdchem.BondStereo.STEREOZ: 2}
new_obmol = pybel.ob.OBMol()
new_obmol.BeginModify()
# Assigning atoms
for index, each_atom in enumerate(rdkit_rwmol.GetAtoms()):
new_atom = new_obmol.NewAtom()
new_atom.SetAtomicNum(each_atom.GetAtomicNum())
new_atom.SetFormalCharge(each_atom.GetFormalCharge())
try:
# Update to OpenBabel 3.X, SetImplicitValence was removed
new_atom.SetImplicitHCount(each_atom.GetNumImplicitHs())
except AttributeError:
# This is OpenBabel 2.4
new_atom.SetImplicitValence(each_atom.GetImplicitValence())
if each_atom.GetIsAromatic():
new_atom.SetAromatic()
new_atom.SetVector(rdkit_rwmol.GetConformer().GetAtomPosition(index).x,
rdkit_rwmol.GetConformer().GetAtomPosition(index).y,
rdkit_rwmol.GetConformer().GetAtomPosition(index).z)
# Assigning bonds
for each_bond in rdkit_rwmol.GetBonds():
new_obmol.AddBond(each_bond.GetBeginAtomIdx() + 1, each_bond.GetEndAtomIdx() + 1,
_bondorders[each_bond.GetBondType()])
if each_bond.GetIsAromatic():
new_obmol.GetBond(each_bond.GetBeginAtomIdx() + 1, each_bond.GetEndAtomIdx() + 1).SetAromatic()
# Copy molecule data
for k, v in rdkit_rwmol.GetPropsAsDict().items():
new_data = pybel.ob.OBPairData()
try:
new_data.SetAttribute(k)
new_data.SetValue(v)
except TypeError:
# Key or value maybe a boost-derived type which cannot be used in SetValue/SetAttribute, ignore and move on
os_util.local_print('Molecule data {}: {} cannot be converted to a RDKit compatible type'.format(k, v),
msg_verbosity=os_util.verbosity_level.info, current_verbosity=verbosity)
else:
new_obmol.CloneData(new_data)
# FIXME: assign stereochemistry
new_obmol.EndModify()
os_util.local_print('Converted rdkit molecule SMILES {} to an openbabel molecule SMILES: {}'
''.format(rdkit.Chem.MolToSmiles(rdkit_rwmol), pybel.Molecule(new_obmol).write('smi')),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.info)
return new_obmol
def obmol_to_rwmol(openbabel_obmol, verbosity=0):
"""Converts an openbabel.OBMol to rdkit.RWMol
Parameters
----------
openbabel_obmol : pybel.ob.OBMol
The OBMol to be converted
verbosity : int
Sets verbosity level
Returns
-------
rdkit.Chem.Mol
Converted molecule
"""
try:
from openbabel import pybel
except ImportError:
import pybel
if verbosity < os_util.verbosity_level.extra_debug:
pybel.ob.obErrorLog.SetOutputLevel(pybel.ob.obError)
else:
os_util.local_print('OpenBabel warning messages are on, expect a lot of output.',
msg_verbosity=os_util.verbosity_level.extra_debug, current_verbosity=verbosity)
if isinstance(openbabel_obmol, rdkit.Chem.Mol):
os_util.local_print('Entering obmol_to_rwmol. Molecule {} (Props: {}) is already a rdkit.Chem.Mol object!'
''.format(openbabel_obmol, openbabel_obmol.GetPropsAsDict()),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.warning)
return openbabel_obmol
elif isinstance(openbabel_obmol, pybel.Molecule):
openbabel_obmol = openbabel_obmol.OBMol
elif not isinstance(openbabel_obmol, pybel.ob.OBMol):
os_util.local_print('Entering obmol_to_rwmol. Molecule {} is a {}, but pybel.Molecule or pybel.ob.OBMol '
'required.'
''.format(openbabel_obmol, type(openbabel_obmol)),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.error)
raise ValueError('pybel.Molecule or pybel.ob.OBMol expected, got {} instead'.format(type(openbabel_obmol)))
# Set some lookups
_bondtypes = {0: rdkit.Chem.BondType.UNSPECIFIED,
1: rdkit.Chem.BondType.SINGLE,
2: rdkit.Chem.BondType.DOUBLE,
3: rdkit.Chem.BondType.TRIPLE,
5: rdkit.Chem.BondType.AROMATIC}
_bondstereo = {0: rdkit.Chem.rdchem.BondStereo.STEREONONE,
1: rdkit.Chem.rdchem.BondStereo.STEREOE,
2: rdkit.Chem.rdchem.BondStereo.STEREOZ}
_atomstereo = {1: rdkit.Chem.rdchem.ChiralType.CHI_TETRAHEDRAL_CW,
2: rdkit.Chem.rdchem.ChiralType.CHI_TETRAHEDRAL_CCW,
3: rdkit.Chem.rdchem.ChiralType.CHI_UNSPECIFIED}
_bondtypes_names = {0: "UNSPECIFIED", 1: "SINGLE", 2: "DOUBLE", 3: "TRIPLE", 5: "AROMATIC"}
_implicit_h_ref = 4294967294
rdmol = rdkit.Chem.Mol()
rdedmol = rdkit.Chem.RWMol(rdmol)
# Use pybel write to trigger residue data evaluation, otherwise we get and StopIteration error
pybel.Molecule(openbabel_obmol).write('pdb')
try:
residue_iter = pybel.ob.OBResidueIter(openbabel_obmol).__next__()
except StopIteration:
os_util.local_print('Could not read atom names from molecule "{}" (Smiles: {})'
''.format(openbabel_obmol.GetTitle(), pybel.Molecule(openbabel_obmol).write('smi')),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.warning)
residue_iter = None
ob_to_rw_atom_map = {}
# Assigning atoms
dummy_atoms = set()
for index, each_atom in enumerate(pybel.ob.OBMolAtomIter(openbabel_obmol)):
is_dummy_atom = False
if residue_iter is not None and residue_iter.GetAtomID(each_atom)[0:2].upper() in ['LP', 'XX'] \
and each_atom.GetAtomicMass() == 0:
is_dummy_atom = True
rdatom = rdkit.Chem.MolFromSmarts('*').GetAtomWithIdx(0)
os_util.local_print('Atom {} was detected as a lone pair because of its name {} and its mass {}'
''.format(index, residue_iter.GetAtomID(each_atom), each_atom.GetAtomicMass()),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.info)
elif residue_iter is None and each_atom.GetAtomicMass() == 0:
is_dummy_atom = True
rdatom = rdkit.Chem.MolFromSmarts('*').GetAtomWithIdx(0)
os_util.local_print('Atom {} was detected as a lone pair because of its mass {} (Note: it was not possible '
'to read atom name)'
''.format(index, each_atom.GetAtomicMass()),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.info)
else:
rdatom = rdkit.Chem.Atom(each_atom.GetAtomicNum())
new_atom = rdedmol.AddAtom(rdatom)
rdedmol.GetAtomWithIdx(new_atom).SetFormalCharge(each_atom.GetFormalCharge())
if residue_iter is not None:
rdedmol.SetProp('_TriposAtomName', residue_iter.GetAtomID(each_atom))
if each_atom.IsAromatic():
rdedmol.GetAtomWithIdx(new_atom).SetIsAromatic(True)
ob_to_rw_atom_map[each_atom.GetId()] = new_atom
if is_dummy_atom:
dummy_atoms.add(new_atom)
rw_to_ob_atom_map = {v: k for k, v in ob_to_rw_atom_map.items()}
if dummy_atoms:
os_util.local_print('These are the dummy atoms detected: dummy_atoms={}'.format(dummy_atoms),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.debug)
else:
os_util.local_print('No dummy atoms detected.'.format(dummy_atoms),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.debug)
# Assigning bonds
for each_bond in pybel.ob.OBMolBondIter(openbabel_obmol):
rdedmol.AddBond(ob_to_rw_atom_map[each_bond.GetBeginAtom().GetId()],
ob_to_rw_atom_map[each_bond.GetEndAtom().GetId()],
_bondtypes[each_bond.GetBondOrder()])
this_bond = rdedmol.GetBondBetweenAtoms(ob_to_rw_atom_map[each_bond.GetBeginAtom().GetId()],
ob_to_rw_atom_map[each_bond.GetEndAtom().GetId()])
if each_bond.IsAromatic():
this_bond.SetIsAromatic(True)
this_bond.SetBondType(_bondtypes[5])
# This bond contains a dummy atom, converting bond to a UNSPECIFIED
if dummy_atoms.intersection({ob_to_rw_atom_map[each_bond.GetBeginAtom().GetId()],
ob_to_rw_atom_map[each_bond.GetEndAtom().GetId()]}):
this_bond.SetBondType(_bondtypes[0])
os_util.local_print('Bond between atoms {} and {} converted to {} type'
''.format(ob_to_rw_atom_map[each_bond.GetBeginAtom().GetId()],
ob_to_rw_atom_map[each_bond.GetEndAtom().GetId()],
_bondtypes_names[each_bond.GetBondOrder()]),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.debug)
stereofacade = pybel.ob.OBStereoFacade(openbabel_obmol)
for each_bond in pybel.ob.OBMolBondIter(openbabel_obmol):
this_bond = rdedmol.GetBondBetweenAtoms(ob_to_rw_atom_map[each_bond.GetBeginAtom().GetId()],
ob_to_rw_atom_map[each_bond.GetEndAtom().GetId()])
if stereofacade.HasCisTransStereo(each_bond.GetId()):
bond_info = stereofacade.GetCisTransStereo(each_bond.GetId())
bond_config = bond_info.GetConfig()
if not bond_info.IsSpecified():
os_util.local_print('Bond {} has cis/trans stereochemistry, but configuration is undefined. This '
'should not happen. I will try to guess the configuration from the 3D structure.',
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
continue
for (i, j) in itertools.combinations(bond_config.refs, r=2):
if not bond_info.IsOnSameAtom(i, j):
rdkit.RDLogger.DisableLog('rdApp.*')
try:
this_bond.SetStereoAtoms(ob_to_rw_atom_map[openbabel_obmol.GetAtomById(i).GetId()],
ob_to_rw_atom_map[openbabel_obmol.GetAtomById(j).GetId()])
except (RuntimeError, AttributeError):
rdkit.RDLogger.EnableLog('rdApp.*')
continue
else:
os_util.local_print('Stereochemistry of the bond between atoms {} and {} is defined by atoms '
'{} and {}.'
''.format(this_bond.GetBeginAtom().GetIdx(),
this_bond.GetEndAtom().GetIdx(),
ob_to_rw_atom_map[openbabel_obmol.GetAtomById(i).GetId()],
ob_to_rw_atom_map[openbabel_obmol.GetAtomById(j).GetId()]),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
rdkit.RDLogger.EnableLog('rdApp.*')
this_bond.SetStereo(rdkit.Chem.rdchem.BondStereo.STEREOTRANS if bond_info.IsTrans(i, j)
else rdkit.Chem.rdchem.BondStereo.STEREOCIS)
os_util.local_print('Setting bond between atoms {} and {} to be a {}.'
''.format(this_bond.GetBeginAtom().GetIdx(), this_bond.GetEndAtom().GetIdx(),
this_bond.GetStereo()),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
break
else:
os_util.local_print('Bond between atoms {} and {} is a double bond, but configuration is undefined. '
'This should not happen. I will try to guess the configuration from the 3D '
'structure, if possible.'
''.format(bond_config.begin, bond_config.end),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
rdkit.Chem.AssignStereochemistry(rdedmol)
for each_atom in pybel.ob.OBMolAtomIter(openbabel_obmol):
if stereofacade.HasTetrahedralStereo(each_atom.GetId()):
if openbabel_obmol.Has3D():
os_util.local_print('Atom {} {} is a stereocenter. Because there are 3D coordinates, I will guess the '
'stereochemistry from it when converting OBMol to RDKit Mol.'
''.format(pybel.Atom(each_atom).type, pybel.Atom(each_atom).idx),
msg_verbosity=os_util.verbosity_level.info, current_verbosity=verbosity)
continue
rdkit_atom = rdedmol.GetAtomWithIdx(ob_to_rw_atom_map[each_atom.GetId()])
tetstereo = stereofacade.GetTetrahedralStereo(each_atom.GetId())
if not tetstereo.GetConfig().specified:
os_util.local_print('Atom {} {} is a stereocenter, but configuration is undefined. Make sure that this '
'is what you want.'
''.format(pybel.Atom(each_atom).type, pybel.Atom(each_atom).idx),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
rdkit_atom.SetChiralTag(_atomstereo[3])
continue
# This code tries to get the CW/ACW info using the same view RDKit uses to do so [1, 2]. This is done via
# the GetConfig [3, 4], then the sequence is tested to match the bond sequence in RDKit.
# References:
# [1] https://sourceforge.net/p/rdkit/mailman/message/36796561/
# [2] https://www.rdkit.org/docs/cppapi/classRDKit_1_1Atom.html
# [3] https://open-babel.readthedocs.io/en/latest/Stereochemistry/stereo.html
# [4] https://openbabel.org/dev-api/structOpenBabel_1_1OBTetrahedralStereo_1_1Config.shtml
bonded_atoms_data = {b.GetBeginAtomIdx() if b.GetBeginAtomIdx() != rdkit_atom.GetIdx()
else b.GetEndAtomIdx(): b.GetIdx() for b in rdkit_atom.GetBonds()}
bonded_atoms_by_bond = sorted(bonded_atoms_data, key=bonded_atoms_data.get)
first_bonded_atom = bonded_atoms_by_bond[0]
input_config = tetstereo.GetConfig(rw_to_ob_atom_map[first_bonded_atom], pybel.ob.OBStereo.Clockwise,
pybel.ob.OBStereo.ViewFrom)
# Now test whether the 2nd and 3rd bonds are in the CW sequence. First, get all the three possible sequences
# using RDKit ids.
atom_sequences = []
for i, j in [[0, 1], [1, 2], [2, 0]]:
try:
atom_sequences.append([ob_to_rw_atom_map[input_config.refs[i]],
ob_to_rw_atom_map[input_config.refs[j]]])
except KeyError as error:
if error.args[0] != _implicit_h_ref:
raise error
# Now get the 2nd and 3rd atoms. Note that implicit Hs are assumed to be in the end of the bond list (see
# [2] above.
if [bonded_atoms_by_bond[1], bonded_atoms_by_bond[2]] in atom_sequences:
# This is means that RDKit would assign CW to this center
this_stereo = _atomstereo[pybel.ob.OBStereo.Clockwise]
else:
this_stereo = _atomstereo[pybel.ob.OBStereo.AntiClockwise]
rdkit_atom.SetChiralTag(this_stereo)
os_util.local_print('Setting stereocenter in atom {} {} to {}'
''.format(pybel.Atom(each_atom).type, pybel.Atom(each_atom).idx, this_stereo),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
rdmol = rdedmol.GetMol()
# Fix the formal charge of N. From Greg Landrum's Gist at
# https://gist.github.com/greglandrum/7f546d1e35c2df537c68a64d887793b8
for each_atom in rdmol.GetAtoms():
if each_atom.GetAtomicNum() == 7 and sum([b.GetBondTypeAsDouble() for b in each_atom.GetBonds()]) == 4 \
and each_atom.GetFormalCharge() == 0:
each_atom.SetFormalCharge(1)
os_util.local_print('Setting formal charge of N{} to 1.'.format(each_atom.GetIdx()),
msg_verbosity=os_util.verbosity_level.info, current_verbosity=verbosity)
try:
rdmol.UpdatePropertyCache()
except ValueError as error:
os_util.local_print('Failed to convert molecule {} to RWMol. Retrying with strict=False. Check your input. '
'Error was {}.'.format(openbabel_obmol.GetTitle(), error),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
try:
rdmol.UpdatePropertyCache(strict=False)
except ValueError as error:
os_util.local_print('Failed to convert molecule {} to RWMol. Cannot go on. Please, check your input. Error '
'was {}.'.format(openbabel_obmol.GetTitle(), error),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
raise error
rdkit.RDLogger.DisableLog('rdApp.*')
# Copy coordinates, first generate at least one conformer
rdkit.Chem.AllChem.EmbedMolecule(rdmol, useRandomCoords=True, maxAttempts=1000, enforceChirality=True,
ignoreSmoothingFailures=True)
rdkit.RDLogger.EnableLog('rdApp.*')
if rdmol.GetNumConformers() != 1:
os_util.local_print('Failed to generate coordinates to molecule',
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.error)
raise ValueError
for atom_obmol in pybel.ob.OBMolAtomIter(openbabel_obmol):
atom_rdkit = ob_to_rw_atom_map[atom_obmol.GetId()]
this_position = rdkit.Geometry.rdGeometry.Point3D()
this_position.x = atom_obmol.x()
this_position.y = atom_obmol.y()
this_position.z = atom_obmol.z()
rdmol.GetConformer().SetAtomPosition(atom_rdkit, this_position)
# Copy data
[rdmol.SetProp(k, v) for k, v in pybel.MoleculeData(openbabel_obmol).items()]
rdmol.SetProp('_Name', openbabel_obmol.GetTitle())
for each_atom in rdmol.GetAtoms():
if each_atom.GetBonds() != ():
continue
import numpy
dist_list = numpy.argsort(numpy.array(rdkit.Chem.AllChem.Get3DDistanceMatrix(rdmol)[each_atom.GetIdx()]))
try:
closer_atom = int(dist_list[1])
except IndexError:
# Is this a lone atom?
continue
rdedmol = rdkit.Chem.RWMol(rdmol)
rdedmol.AddBond(each_atom.GetIdx(), closer_atom)
rdmol = rdedmol.GetMol()
rdkit.Chem.SanitizeMol(rdmol)
os_util.local_print('Atom id: {} is not explicitly bonded to any atom in molecule, connecting it to the closer '
'atom id: {}'.format(each_atom.GetIdx(), closer_atom),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.warning)
rdkit.Chem.SanitizeMol(rdmol)
rdkit.Chem.AssignStereochemistry(rdmol)
rdkit.Chem.AssignStereochemistryFrom3D(rdmol)
os_util.local_print("obmol_to_rwmol converted molecule {} (name: {}). Pybel SMILES: {} to rdkit SMILES: {}"
"".format(openbabel_obmol,
openbabel_obmol.GetTitle() if openbabel_obmol.GetTitle() else '<<UNNAMED>>',
pybel.Molecule(openbabel_obmol).write('smi').replace('\n', ''),
rdkit.Chem.MolToSmiles(rdmol)),
current_verbosity=verbosity, msg_verbosity=os_util.verbosity_level.debug)
return rdmol
def rdmol_com(input_mol, conformer=-1):
""" Calculates COM of a rdkit.Chem.Mol
:param rdkit.Chem.Mol input_mol: molecule
:param conformer: COM of this conformer (-1: auto)
:rtype: numpy.array
"""
# TODO: rewrite this for speed using Get
return numpy.reshape(numpy.array([each_atom_pos * each_atom.GetMass()
for each_atom in input_mol.GetAtoms()
for each_atom_pos in
input_mol.GetConformer(conformer).GetAtomPosition(each_atom.GetIdx())]),
[-1, 3]).mean(0)
def initialise_neutralization_reactions():
""" Prepare SMARTS patterns to neutralize_molecule. Borrowed from rdkit cookbook
(http://www.rdkit.org/docs/Cookbook.html#neutralizing-charged-molecules)
:rtype: list
"""
patts = (
# Imidazoles
('[n+;H]', 'n'),
# Amines
('[N+;!H0]', 'N'),
# Carboxylic acids and alcohols
('[$([O-]);!$([O-][#7])]', 'O'),
# Thiols
('[S-;X1]', 'S'),
# Sulfonamides
('[$([N-;X2]S(=O)=O)]', 'N'),
# Enamines
('[$([N-;X2][C,N]=C)]', 'N'),
# Tetrazoles
('[n-]', '[nH]'),
# Sulfoxides
('[$([S-]=O)]', 'S'),
# Amides
('[$([N-]C=O)]', 'N'),
)
return [(rdkit.Chem.MolFromSmarts(x), rdkit.Chem.MolFromSmiles(y, False)) for x, y in patts]
def neutralize_molecule(mol, reactions=None):
""" Neutralizes a molecule. Borrowed from rdkit cookbook
(http://www.rdkit.org/docs/Cookbook.html#neutralizing-charged-molecules)
:param rdkit.Chem.Mol mol: molecule to be neutralized
:param list reactions: use these reactions instead of default ones
:rtype: set
"""
global _reactions
if reactions is None:
if _reactions is None:
_reactions = initialise_neutralization_reactions()
reactions = _reactions
replaced = False
for i, (reactant, product) in enumerate(reactions):
while mol.HasSubstructMatch(reactant):
replaced = True
rms = rdkit.Chem.AllChem.ReplaceSubstructs(mol, reactant, product)
mol = rms[0]
return mol, replaced
_reactions = None
def verify_molecule_name(molecule, moldict, new_default_name=None, verbosity=0):
""" Verify the a molecule name exists and is unique and return a valid name the molecule
:param [rdkit.Chem.Mol, str] molecule: molecule to be verified
:param dict moldict: dict of read molecules
:param str new_default_name: if molecule lacks a name, use this name instead (Default: generate a random name)
:param int verbosity: controls the verbosity level
:rtype: str
"""
if isinstance(molecule, rdkit.Chem.Mol):
try:
this_mol_name = molecule.GetProp('_Name')
except KeyError:
this_mol_name = None
else:
if not molecule:
this_mol_name = None
else:
this_mol_name = molecule
if this_mol_name is None:
if new_default_name:
this_mol_name = new_default_name
else:
this_mol_name = '(mol_{})'.format(numpy.random.randint(1, 999999999))
while this_mol_name in moldict:
this_mol_name = '(mol_{})'.format(numpy.random.randint(1, 999999999))
if isinstance(molecule, rdkit.Chem.Mol):
os_util.local_print('Molecule {} have no name. Molecule name is used to save molecule data '
'and serves as an index. I will generate a random name for it, namely: {}'
''.format(rdkit.Chem.MolToSmiles(molecule), this_mol_name),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
else:
os_util.local_print('A molecule have no name. Molecule name is used to save molecule data '
'and serves as an index. I will generate a random name for it, namely: {}'
''.format(this_mol_name),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
this_mol_name = os.path.splitext(os.path.basename(this_mol_name))[0]
if this_mol_name in moldict:
colliding_name = this_mol_name
this_mol_name = '{}_1'.format(this_mol_name)
while this_mol_name in moldict:
this_mol_name = this_mol_name[:-1] + str(int(this_mol_name[-1]) + 1)
if isinstance(molecule, rdkit.Chem.Mol):
os_util.local_print('Two molecules (Smiles: {} and {}) have the same name {}. Molecule name is used to '
'save molecule data and serves as an index. I will rename molecule {} to {}'
''.format(rdkit.Chem.MolToSmiles(moldict[colliding_name]),
rdkit.Chem.MolToSmiles(molecule), colliding_name,
rdkit.Chem.MolToSmiles(molecule), this_mol_name),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
else:
os_util.local_print('Two molecules have the same name {}. Molecule name is used to '
'save molecule data and serves as an index. I will rename the last molecule {}'
''.format(colliding_name, this_mol_name),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
if isinstance(molecule, rdkit.Chem.Mol):
molecule.SetProp('_Name', this_mol_name)
return this_mol_name
def process_dummy_atoms(molecule, verbosity=0):
""" Sanitizes dummy atoms in a rdkit.Chem.Mol
:param rdkit.Chem.Mol molecule: molecule to be verified
:param int verbosity: controls the verbosity level
:rtype: rdkit.Chem.Mol
"""
os_util.local_print('Entering process_dummy_atoms(molecule=({}; SMILES={}), verbosity={})'
''.format(molecule, rdkit.Chem.MolToSmiles(molecule), verbosity),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
# Iterates over a copy of molecule ahd convert query atoms to dummy atoms, adding bonds if necessary
temp_mol = rdkit.Chem.Mol(molecule)
for atom_idx, each_atom in enumerate(temp_mol.GetAtoms()):
if isinstance(each_atom, rdkit.Chem.rdchem.QueryAtom):
newdummy = rdkit.Chem.Atom(0)
rdedmol = rdkit.Chem.RWMol(molecule)
rdedmol.ReplaceAtom(atom_idx, newdummy, preserveProps=True)
molecule = rdedmol.GetMol()
if each_atom.GetProp('_TriposAtomName')[:2] == 'LP':
os_util.local_print('Lone pair found. Atom with id {} was assumed a lone pair by its name ({}) and '
'its type ({}). If this is wrong, please change the atom name.'
''.format(atom_idx, each_atom.GetProp('_TriposAtomName'),
each_atom.GetProp('_TriposAtomType')),
msg_verbosity=os_util.verbosity_level.info, current_verbosity=verbosity)
if each_atom.GetBonds() == ():
if temp_mol.GetNumConformers() == 0:
os_util.local_print('Disconnected lone pair atom found in a molecule with no 3D coordinates. '
'3D coordinates are used to guess the LP host, but are absent in molecule '
'{}. I cannot continue.'.format(temp_mol),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
raise SystemExit(1)
rdkit.Chem.rdMolTransforms.CanonicalizeConformer(temp_mol.GetConformer(0))
# LP is not bonded to any other atom. Connect it to the closer one
import numpy
temp_mol.GetConformer(0)
dist_list = numpy.argsort(numpy.array(rdkit.Chem.Get3DDistanceMatrix(temp_mol)[atom_idx]))
closer_atom = int(dist_list[1])
rdedmol = rdkit.Chem.RWMol(molecule)
rdedmol.AddBond(atom_idx, closer_atom)
molecule = rdedmol.GetMol()
rdkit.Chem.SanitizeMol(molecule)
os_util.local_print('Lonepair {} (id: {}) is not explicitly bonded to any atom in molecule, '
'connecting it to the closer atom {} (id: {}). Please, check the output'
''.format(molecule.GetAtomWithIdx(atom_idx).GetProp('_TriposAtomName'),
atom_idx,
molecule.GetAtomWithIdx(closer_atom).GetProp('_TriposAtomName'),
closer_atom),
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
else:
# FIXME: support other dummy atoms (eg: in linear molecules)
os_util.local_print('The molecule {} contains dummy atoms which are not lonepairs. This is not '
'supported.'.format(molecule),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
raise SystemExit(1)
return molecule
@os_util.trace
def adjust_query_properties(query_molecule, generic_atoms=False, ignore_charge=True, ignore_isotope=True, verbosity=0):
""" Adjust query settings removing all charges, isotope, aromaticity and valence info from core_structure SMARTS
:param rdkit.Chem.Mol query_molecule: query molecule
:param bool generic_atoms: make atoms generic
:param bool ignore_charge: set all atomic charges to 0
:param bool ignore_isotope: ignore atomic isotopes
:param int verbosity: controls the verbosity level
:rtype: rdkit.Chem.Mol
"""
new_query_molecule = rdkit.Chem.Mol(query_molecule)
# Parameters to GetSubstructMatch
query_m = rdkit.Chem.rdmolops.AdjustQueryParameters()
query_m.makeBondsGeneric = True
query_m.makeDummiesQueries = True
query_m.adjustDegree = False
if generic_atoms:
query_m.makeAtomsGeneric = True
else:
if ignore_isotope:
[a0.SetQuery(rdkit.Chem.MolFromSmarts('[#{}]'.format(a0.GetAtomicNum())).GetAtomWithIdx(0))
for a0 in new_query_molecule.GetAtoms() if isinstance(a0, rdkit.Chem.QueryAtom)]
if ignore_charge:
[a0.SetFormalCharge(0) for a0 in new_query_molecule.GetAtoms()]
new_query_molecule = rdkit.Chem.AdjustQueryProperties(new_query_molecule, query_m)
os_util.local_print('The molecule {} (SMARTS={}) was altered by adjust_query_properties to {} (SMARTS={})'
''.format(query_molecule, rdkit.Chem.MolToSmarts(query_molecule),
new_query_molecule, rdkit.Chem.MolToSmarts(new_query_molecule)),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
return new_query_molecule
def num_explict_hydrogens(mol):
""" Return the number of explicit hydrogens in molecular graph
:param rdkit.Chem.Mol mol: the input molecule
:rtype: int
"""
return sum([1 for i in mol.GetAtoms() if i.GetAtomicNum() == 1])
def num_implicit_hydrogens(mol):
"""Return the number of implicit hydrogens in the molecular graph
Parameters
----------
mol : rdkit.Chem.Mol
Input molecule
Returns
-------
int
"""
return sum([each_atom.GetNumImplicitHs() for each_atom in mol.GetAtoms()])
def loose_replace_side_chains(mol, core_query, use_chirality=False, verbosity=True):
""" Reconstruct a molecule based on common core. First, try to use the regular query. If fails, fallback to
generalized bonds then generalized atoms.
:param rdkit.Chem.Mol mol: the molecule to be modified
:param rdkit.Chem.Mol core_query: the molecule to be used as a substructure query for recognizing the core
:param bool use_chirality: match the substructure query using chirality
:param int verbosity: set verbosity level
:rtype: rdkit.Chem.Mol
"""
temp_core_structure = rdkit.Chem.Mol(core_query)
if num_explict_hydrogens(core_query) > 0 and num_explict_hydrogens(mol) == 0:
os_util.local_print('loose_replace_side_chains was called with a mol without explict hydrogens and a '
'core_query with {} explict hydrogens. Removing core_query explict Hs.'
''.format(num_explict_hydrogens(core_query)),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
editable_core = rdkit.Chem.EditableMol(core_query)
hydrogen_atoms = [each_atom.GetIdx() for each_atom in core_query.GetAtoms() if each_atom.GetAtomicNum() == 1]
for idx in sorted(hydrogen_atoms, reverse=True):
editable_core.RemoveAtom(idx)
temp_core_structure = editable_core.GetMol()
rdkit.Chem.SanitizeMol(temp_core_structure, catchErrors=True)
result_core_structure = rdkit.Chem.ReplaceSidechains(mol, temp_core_structure, useChirality=use_chirality)
if result_core_structure is None:
os_util.local_print('rdkit.Chem.ReplaceSidechains failed with mol={} (SMILES="{}") and coreQuery={} '
'(SMARTS="{}"). Retrying with adjust_query_properties.'
''.format(mol, rdkit.Chem.MolToSmiles(mol), temp_core_structure,
rdkit.Chem.MolToSmarts(temp_core_structure)),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
temp_core_mol = adjust_query_properties(temp_core_structure, verbosity=verbosity)
result_core_structure = rdkit.Chem.ReplaceSidechains(mol, temp_core_mol, useChirality=use_chirality)
if result_core_structure is None:
os_util.local_print('rdkit.Chem.ReplaceSidechains failed with mol={} (SMILES="{}") and coreQuery={} '
'(SMARTS="{}"). Retrying with adjust_query_properties setting generic_atoms=True.'
''.format(mol, rdkit.Chem.MolToSmiles(mol), temp_core_structure,
rdkit.Chem.MolToSmarts(temp_core_structure)),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
temp_core_mol = adjust_query_properties(temp_core_structure, generic_atoms=True, verbosity=verbosity)
result_core_structure = rdkit.Chem.ReplaceSidechains(mol, temp_core_mol, useChirality=use_chirality)
return result_core_structure
def has_3d(temp_mol, conf_id=-1, tolerance=1e-5, verbosity=0):
"""Check if molecules has 3D coordinates. Based upon OpenBabel OBMol::Has3D()
Parameters
----------
temp_mol : rdkit.Chem.Mol
Molecule to be checked.
conf_id : int
Test 3D for this conformation.
tolerance : float
Tolerance, in angstroms, for a value to be taken as not null.
verbosity : int
Sets the verbosity level.
Returns
-------
bool
True if molecule has 3D coordinates, False otherwise.
"""
positions_array = temp_mol.GetConformer(conf_id).GetPositions()
return not numpy.allclose(positions_array, 0.0, atol=tolerance, rtol=0.0)
def translation_to_4by4_mat(translation, verbosity=0):
""" Convert a translation to a 4-by-4 matrix format
Parameters
----------
translation : numpy.array
Operation to be converted
verbosity : int
Sets the verbosity level
Returns
-------
numpy.array
4-by-4 translation matrix
"""
if translation.shape == (3,):
ret_mat = numpy.identity(4)
ret_mat[:3, 3] = translation[:3]
elif translation.shape == (4, 4):
if numpy.identity(4)[: 4, :3] != translation[: 4, :3]:
os_util.local_print('{} is not a valid translation matrix!',
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
ret_mat = translation
else:
os_util.local_print('Cannot understand translation {}. It must be either a 3-element vector or a 4-by-4 '
'matrix'.format(translation), msg_verbosity=os_util.verbosity_level.error,
current_verbosity=verbosity)
raise TypeError
return ret_mat
def rotation_to_4by4_mat(rotation, verbosity=0):
""" Convert a rotation to a 4-by-4 matrix format
Parameters
----------
rotation : numpy.array
Operation to be converted
verbosity : int
Sets the verbosity level
Returns
-------
numpy.array
4-by-4 rotation matrix
"""
if rotation.shape == (3, 3):
ret_mat = numpy.identity(4)
ret_mat[:3, :3] = rotation
elif rotation.shape == (4, 4):
if numpy.allclose(numpy.identity(4)[3, :], rotation[:, 3]) \
and numpy.allclose(numpy.identity(4)[3, :], rotation[3, :]):
os_util.local_print('{} is not a valid rotation matrix!',
msg_verbosity=os_util.verbosity_level.warning, current_verbosity=verbosity)
ret_mat = rotation
else:
os_util.local_print('Cannot understand rotation {}. It must be either a 3-by-3 or a 4-by-4 matrix'
''.format(rotation),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
raise TypeError
return ret_mat
def parameterize_small_molecule(input_molecule, param_type='acpype', executable=None, charge_method=None,
atom_type=None, output_dir=None, verbosity=0, **kwargs):
""" Parameterize small molecules using AcPYPE.
Parameters
----------
input_molecule : rdkit.Chem.Mol
Molecule to be parameterized
param_type : str
Which parameterization software to use. Currently, only "acpype" is supported
executable : str
Use this executable to parameterize
charge_method : str
Selects the charge method. Choices are "bcc", "user", and "gas" for AcPYPE
atom_type : str
Which atom type to use. Only relevant for AcPYPE. Choices are 'gaff', 'amber', 'gaff2', 'amber2'
output_dir : str
Save generated topology files to this directory
verbosity : int
Sets the verbosity level
Returns
-------
dict
'topology': list containing the generated topology files for this ligand; 'keepfiles': list of files extra files
copied.
"""
from distutils.dir_util import copy_tree
if not output_dir:
output_dir = os.getcwd()
return_data = {'topology': []}
top_files = return_data['topology']
timeout = kwargs.get('timeout', None)
die_on_error = kwargs.get('die_on_error', True)
# Extra files to keep from the acpype result
keepfiles = kwargs.get('keepfiles', [])
if param_type == 'acpype':
# TODO: use AcPYPE API instead of subprocess. But it will require acpype to be installed in the same env as
# PyAutoFEP, which may or may not be the best option for the user. Maybe using the API, then falling back to
# subprocess would be ideal.
from subprocess import run, CalledProcessError, TimeoutExpired
from tempfile import TemporaryDirectory
ligand_file = '{}.mdl'.format(input_molecule.GetProp('_Name'))
# Assemble AcPYPE command line using:
# acpype -i _file_ -c _string_ -n _int_ -a _string_
if not executable:
executable = 'acpype'
# Total charge is pre-calculated using rdkit. In case acpype guesses it wrong, topology would be wrong.
cmd_line = [executable, '-i', ligand_file, '-n', str(rdkit.Chem.GetFormalCharge(input_molecule)), '-o', 'gmx',
'-b', input_molecule.GetProp('_Name')]
# Parse charge method
if charge_method is None:
charge_method = 'bcc'
if charge_method in ['gas', 'bcc', 'user']:
cmd_line += ['-c', charge_method]
elif charge_method:
# User supplied a charge_method, but acpype won't support it
os_util.local_print('AcPYPE allowed charge methods are "gas","bcc", and "user". Charge method {} cannot be'
'used.'.format(charge_method),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
if die_on_error:
raise ValueError('Incompatible charge method {}.'.format(charge_method))
else:
return False
# Parse atom types
if atom_type is None:
atom_type = 'gaff2'
if atom_type in ['gaff', 'amber', 'gaff2', 'amber2']:
cmd_line += ['-a', atom_type]
elif atom_type:
os_util.local_print('AcPYPE allowed atom types are "gaff", "amber", "gaff2", "amber2". Atom typing {} '
'cannot be used.'.format(atom_type),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
if die_on_error:
raise ValueError('Incompatible atom type {}.'.format(atom_type))
else:
return False
for each_option in ['qprog', 'max_time']:
if each_option in kwargs:
cmd_line.extend(['--{}'.format(each_option), str(kwargs[each_option])])
# Check whether ligand parameters are already there. It is done here to make sure we have already set
# charge_method and atom_type vars
new_files = [os.path.join(output_dir, input_molecule.GetProp('_Name') + '_GMX' + each_ext)
for each_ext in ['.itp', '.top']]
new_files.append(os.path.join(output_dir, 'posre_' + input_molecule.GetProp('_Name') + '.itp'))
if all([os.path.isfile(f) for f in new_files]):
# Parameter files exists in output_dir. Check for extra files.
if keepfiles:
new_keepfiles = []
for each_extra in keepfiles:
if each_extra == 'mol2':
mol2_file = '{}_{}_{}.mol2'.format(input_molecule.GetProp('_Name'), charge_method, atom_type)
if os.path.isfile(os.path.join(output_dir, mol2_file)):
new_keepfiles.append(mol2_file)
else:
os_util.local_print('Output topologies for ligand {} are aleady present in {}, but the '
'acpype mol2 file {} was no found. Parameterizing molecule again.'
''.format(input_molecule.GetProp('_Name'), output_dir, mol2_file),
msg_verbosity=os_util.verbosity_level.warning,
current_verbosity=verbosity)
break
if each_extra == 'all':
alldir = '{}.acpype'.format(input_molecule.GetProp('_Name'))
if os.path.isdir(os.path.join(output_dir, alldir)):
new_keepfiles.append(alldir)
else:
os_util.local_print('Output topologies for ligand {} are aleady present in {}, but the '
'acpype output directory {} was not found. Parameterizing molecule '
'again.'
''.format(input_molecule.GetProp('_Name'), output_dir, alldir),
msg_verbosity=os_util.verbosity_level.warning,
current_verbosity=verbosity)
break
else:
return_data['topology'] = new_files
return_data['keepfiles'] = new_keepfiles
return return_data
else:
return_data['topology'] = new_files
return return_data
# Make sure sqm use a single thread. Parallelization of sqm is less efficient than running several instances in
# parallel
this_env = os.environ.copy()
this_env['OMP_NUM_THREADS'] = '1'
os_util.local_print("Parameterizing small molecule {} using {} and the following options: executable={}, "
"charge_method={}, atom_type={}, output_dir={}, verbosity={}. Command line is: \"{}\""
"".format(input_molecule.GetProp('_Name'), param_type, executable, charge_method,
atom_type, output_dir, verbosity, ' '.join(cmd_line)),
msg_verbosity=os_util.verbosity_level.debug, current_verbosity=verbosity)
# Run acpype in a temporary directory
with TemporaryDirectory() as this_tmp_dir:
rdkit.Chem.MolToMolFile(input_molecule, os.path.join(this_tmp_dir, ligand_file))
try:
acpype_run = run(cmd_line, capture_output=True, text=True, cwd=this_tmp_dir, timeout=timeout,
env=this_env)
except TimeoutExpired as error:
os_util.local_print('AcPYPE run failed due to timeout (timeout={}). The complete run command was {}. '
'Error was: {}.'
''.format(timeout, ' '.join(cmd_line), error),
msg_verbosity=os_util.verbosity_level.error, current_verbosity=verbosity)
if die_on_error:
raise error
else: