1
#!/usr/bin/env python
2 100
"""
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 100
import os, sys
20 100
import openeye.oechem as oechem
21 100
import openeye.oeomega as oeomega
22 100
import openeye.oeszybki as oeszybki
23

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

26

27 100
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 100
    molWithConfs = oechem.OEMol(mol)
41 100
    omega = oeomega.OEOmega()
42 100
    maxConfs = 0
43 100
    omega.SetMaxConfs(maxConfs)
44 100
    omega.SetStrictStereo(False)
45 100
    omega.SetSampleHydrogens(True)
46 100
    omega.SetEnumNitrogen(oeomega.OENitrogenEnumeration_All)
47 100
    if not omega(molWithConfs):
48 0
        print("omega failed on %s" % mol.GetTitle())
49 0
        return None
50
    else:
51 100
        return molWithConfs
52

53

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

111

112 100
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 100
    optSzybki = oeszybki.OESzybkiOptions()
130 100
    optSzybki.SetForceFieldType(oeszybki.OEForceFieldType_MMFF94S)
131 100
    optSzybki.SetSolventModel(oeszybki.OESolventModel_Sheffield)
132 100
    optSzybki.SetOptimizerType(oeszybki.OEOptType_SD)
133 100
    taglabel = 'MM Szybki SD Energy'
134

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

149

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

152

153 100
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 100
    base, extension = os.path.splitext(os.path.basename(smiles))
174 100
    if extension == '.sdf':
175 100
        base = base + '_quanformer'
176 100
    sdfout = base + '.sdf'
177

178
    ### Read in smiles file.
179 100
    ifs = oechem.oemolistream()
180 100
    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 100
    ofs = oechem.oemolostream()
185 100
    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 100
    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 100
    conffile = open('numConfs.txt', 'a')
195 100
    conffile.write("Number of original conformers\n")
196

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

206 100
        for i, conf in enumerate(mol.GetConfs()):
207 100
            print(mol.GetTitle(), i + 1)
208
            ### Resolve bad clashes.
209 100
            if resolve_clash:
210 100
                print("  Resolving bad clashes...")
211 100
                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 100
            if do_opt:
217 100
                print("  Doing a quick MM (SD) optimization...")
218 100
                if not quick_opt(conf):
219 0
                    print('  Quick optimization failed for molecule '
220
                        f'{mol.GetTitle()} conformer {i+1}')
221 0
                    continue
222 100
        oechem.OEWriteConstMolecule(ofs, mol)
223

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

229

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

Read our documentation on viewing source code .

Loading