1
#!/usr/bin/env python
2 60
"""
3
initialize_confs.py
4

5
This script takes a list of smiles (or multi-molecule file) and for each:
6
 - generates conformers with Omega
7
 - resolves electrostatic clashes with Cartesian optimization using Szybki
8
 - performs quick optimization with steepest descent opt with Szybki
9
 - writes out multi-mol multi-conformer file
10

11
To use this script as an imported module:
12
- import initialize_confs
13
- initialize_confs.initialize_confs('filename', resolve_clash=True, do_opt=True)
14

15
By: Christopher I. Bayly, Victoria T. Lim
16

17
"""
18

19 60
import os, sys
20 60
import openeye.oechem as oechem
21 60
import openeye.oeomega as oeomega
22 60
import openeye.oeszybki as oeszybki
23

24
### ------------------- Functions -------------------
25

26

27 60
def generate_confs(mol):
28
    """
29
    Generate conformers of molecule from its SMILES string.
30

31
    Parameters
32
    ----------
33
    mol : OEChem molecule
34

35
    Returns
36
    -------
37
    molWithConfs : OEChem molecule with omega-generated conformers
38

39
    """
40 60
    molWithConfs = oechem.OEMol(mol)
41 60
    omega = oeomega.OEOmega()
42 60
    maxConfs = 0
43 60
    omega.SetMaxConfs(maxConfs)
44 60
    omega.SetStrictStereo(False)
45 60
    omega.SetSampleHydrogens(True)
46 60
    omega.SetEnumNitrogen(oeomega.OENitrogenEnumeration_All)
47 60
    if not omega(molWithConfs):
48 0
        print("omega failed on %s" % mol.GetTitle())
49 0
        return None
50
    else:
51 60
        return molWithConfs
52

53

54 60
def resolve_clashes(mol, clashfile):
55
    """
56
    Minimize conformers with severe steric interaction.
57

58
    Parameters
59
    ----------
60
    mol : single OEChem molecule (single conformer)
61
    clashfile : string
62
        name of file to write output
63

64
    Returns
65
    -------
66
    boolean
67
        True if completed successfully, False otherwise.
68

69
    """
70

71
    # set general energy options along with the single-point specification
72 60
    spSzybki = oeszybki.OESzybkiOptions()
73 60
    spSzybki.SetForceFieldType(oeszybki.OEForceFieldType_MMFF94S)
74 60
    spSzybki.SetSolventModel(oeszybki.OESolventModel_Sheffield)
75 60
    spSzybki.SetRunType(oeszybki.OERunType_SinglePoint)
76
    # generate the szybki MMFF94 engine for single points
77 60
    szSP = oeszybki.OESzybki(spSzybki)
78
    # construct minimiz options from single-points options to get general optns
79 60
    optSzybki = oeszybki.OESzybkiOptions(spSzybki)
80
    # now reset the option for minimization
81 60
    optSzybki.SetRunType(oeszybki.OERunType_CartesiansOpt)
82
    # generate szybki MMFF94 engine for minimization
83 60
    szOpt = oeszybki.OESzybki(optSzybki)
84
    # add strong harmonic restraints to nonHs
85 60
    szOpt.SetHarmonicConstraints(10.0)
86
    # construct a results object to contain the results of a szybki calculation
87 60
    szResults = oeszybki.OESzybkiResults()
88
    # work on a copy of the molecule
89 60
    tmpmol = oechem.OEMol(mol)
90 60
    if not szSP(tmpmol, szResults):
91 0
        print('szybki run failed for %s' % tmpmol.GetTitle())
92 0
        return False
93
    #Etotsp = szResults.GetTotalEnergy()
94 60
    Evdwsp = szResults.GetEnergyTerm(oeszybki.OEPotentialTerms_MMFFVdW)
95 60
    if Evdwsp > 35:
96 60
        if not szOpt(tmpmol, szResults):
97 0
            print('szybki run failed for %s' % tmpmol.GetTitle())
98 0
            return False
99
        #Etot = szResults.GetTotalEnergy()
100 60
        Evdw = szResults.GetEnergyTerm(oeszybki.OEPotentialTerms_MMFFVdW)
101 60
        wfile = open(clashfile, 'a')
102 60
        wfile.write(
103
            '%s resolved bad clash: initial vdW: %.4f ; '
104
            'resolved EvdW: %.4f\n' % (tmpmol.GetTitle(), Evdwsp, Evdw))
105 60
        wfile.close()
106 60
        mol.SetCoords(tmpmol.GetCoords())
107 60
    oechem.OESetSDData(mol, oechem.OESDDataPair('MM Szybki Single Point Energy'\
108
, "%.12f" % szResults.GetTotalEnergy()))
109 60
    return True
110

111

112 60
def quick_opt(mol):
113
    """
114
    Fast MM optimization to whittle down number of conformers before QM.
115
    Default Szybki OEOptType type set to steepest descent (SD) based on
116
       preliminary comparisons.
117

118
    Parameters
119
    ----------
120
    mol : single OEChem molecule (single conformer)
121

122
    Returns
123
    -------
124
    boolean
125
        True if completed successfully, False otherwise.
126

127
    """
128
    # set general energy options along with the run type specification
129 60
    optSzybki = oeszybki.OESzybkiOptions()
130 60
    optSzybki.SetForceFieldType(oeszybki.OEForceFieldType_MMFF94S)
131 60
    optSzybki.SetSolventModel(oeszybki.OESolventModel_Sheffield)
132 60
    optSzybki.SetOptimizerType(oeszybki.OEOptType_SD)
133 60
    taglabel = 'MM Szybki SD Energy'
134

135
    # generate szybki MMFF94 engine for minimization
136 60
    szOpt = oeszybki.OESzybki(optSzybki)
137
    # construct a results object to contain the results of a szybki calculation
138 60
    szResults = oeszybki.OESzybkiResults()
139
    # work on a copy of the molecule
140 60
    tmpmol = oechem.OEMol(mol)
141 60
    if not szOpt(tmpmol, szResults):
142 0
        print('szybki run failed for %s' % tmpmol.GetTitle())
143 0
        return False
144 60
    mol.SetCoords(tmpmol.GetCoords())
145 60
    oechem.OESetSDData(mol, oechem.OESDDataPair(taglabel, "%.12f" \
146
        % szResults.GetTotalEnergy()))
147 60
    return True
148

149

150
### ------------------- Script -------------------
151

152

153 60
def initialize_confs(smiles, resolve_clash=True, do_opt=True):
154
    """
155
    From a file containing smiles strings, generate omega conformers,
156
       resolve steric clashes, do a quick MM opt, and write SDF output.
157

158
    Parameters
159
    ----------
160
    smiles : string
161
        Name of the input molecule file. Note: this isn't strictly limited
162
        to SMILES files, can also work with SDF, MOL2, etc. The file can be in
163
        a different location than the directory of which this function is
164
        called. The output SDF and txt files will be placed in the directory
165
        of which the function is called.
166
    resolve_clash : Boolean
167
        True to resolve steric clashes by geometry optimization using OESzybki
168
    do_opt : Boolean
169
        True to run quick geometry optimization using OESzybki with MMFF94S
170
        force field and Sheffield solvent model
171

172
    """
173 60
    base, extension = os.path.splitext(os.path.basename(smiles))
174 60
    if extension == '.sdf':
175 60
        base = base + '_quanformer'
176 60
    sdfout = base + '.sdf'
177

178
    ### Read in smiles file.
179 60
    ifs = oechem.oemolistream()
180 60
    if not ifs.open(smiles):
181 0
        oechem.OEThrow.Warning("Unable to open %s for reading" % smiles)
182

183
    ### Open output file to write molecules.
184 60
    ofs = oechem.oemolostream()
185 60
    if os.path.exists(sdfout):
186 0
        raise FileExistsError(
187
            "Output .sdf file already exists. Exiting initialize_confs.\n{}\n".format(
188
                os.path.abspath(sdfout)))
189 60
    if not ofs.open(sdfout):
190 0
        oechem.OEThrow.Fatal("Unable to open %s for writing" % sdfout)
191

192
    ### Output files detailing number of resolved clashes
193
    ###   and original number of conformers before MM opt.
194 60
    conffile = open('numConfs.txt', 'a')
195 60
    conffile.write("Number of original conformers\n")
196

197
    ### For each molecule: label atoms, generate confs, resolve clashes, optimize.
198 60
    for smimol in ifs.GetOEMols():
199 60
        oechem.OETriposAtomNames(smimol)
200 60
        oechem.OEAddExplicitHydrogens(smimol)
201 60
        mol = generate_confs(smimol)
202 60
        if mol is None:
203 0
            continue
204 60
        conffile.write("%s\t%s\n" % (mol.GetTitle(), mol.NumConfs()))
205

206 60
        for i, conf in enumerate(mol.GetConfs()):
207 60
            print(mol.GetTitle(), i + 1)
208
            ### Resolve bad clashes.
209 60
            if resolve_clash:
210 60
                print("  Resolving bad clashes...")
211 60
                if not resolve_clashes(conf, "numClashes.txt"):
212 0
                    print('  Resolving bad clashes failed for molecule '
213
                        f'{mol.GetTitle()} conformer {i+1}')
214 0
                    continue
215
            ### MM optimization.
216 60
            if do_opt:
217 60
                print("  Doing a quick MM (SD) optimization...")
218 60
                if not quick_opt(conf):
219 0
                    print('  Quick optimization failed for molecule '
220
                        f'{mol.GetTitle()} conformer {i+1}')
221 0
                    continue
222 60
        oechem.OEWriteConstMolecule(ofs, mol)
223

224
    ### Close files.
225 60
    ifs.close()
226 60
    ofs.close()
227 60
    conffile.close()
228

229

230 60
if __name__ == "__main__":
231 0
    initialize_confs(sys.argv[1], sys.argv[2], sys.argv[3])

Read our documentation on viewing source code .

Loading