1
#!/usr/bin/env python
2 100
"""
3
confs_to_psi.py
4

5
Purpose:    Generate Psi4 inputs for many molecules/conformers.
6
By:         Victoria T. Lim, Christopher I. Bayly
7
Version:    Nov 30 2018
8

9
"""
10

11 100
import os, sys
12 100
import openeye.oechem as oechem
13 100
import json
14 100
import quanformer.reader as reader
15

16

17 100
def check_calc(calctype):
18 100
    if calctype not in {'opt', 'spe', 'hess'}:
19 0
        raise ValueError("Specify a valid calculation type.")
20

21 100
def make_psi_input(mol, label, method, basisset, calctype='opt', mem=None):
22
    """
23
    Get coordinates from input mol, and generate/format input text for
24
    Psi4 calculation.
25

26
    Parameters
27
    ----------
28
    mol : OpenEye OEMol
29
        OEMol with coordinates
30
    label : string
31
        Name of molecule with integer identifier (for conformers).
32
    method: string
33
        Name of the method as understood by Psi4. Example: "mp2"
34
    basis : string
35
        Name of the basis set as understood by Psi4. Example: "def2-sv(p)"
36
    calctype : string
37
        What kind of Psi4 calculation to run. Supported inputs are:
38
        'opt' for geometry optimization,
39
        'spe' for single point energy calculation, and
40
        'hess' for Hessian calculation.
41
    memory : string
42
        How much memory each Psi4 job should take. If not specified, the
43
        default in Psi4 is 500 Mb. Examples: "2000 MB" "1.5 GB"
44
        http://www.psicode.org/psi4manual/master/psithoninput.html
45

46
    Returns
47
    -------
48
    inputstring : string
49
        Contents of Psi4 input file to be written out
50

51
    """
52

53
    # check that specified calctype is valid
54 100
    check_calc(calctype)
55

56 100
    inputstring = ""
57

58
    # specify memory requirements, if defined
59 100
    if mem != None:
60 100
        inputstring += "memory %s\n" % mem
61 100
    inputstring += ('molecule %s {\n' % label)
62

63
    # charge and multiplicity; multiplicity hardwired to singlet (usually is)
64 100
    netCharge = oechem.OENetCharge(mol)
65 100
    inputstring += ('  %s 1' % netCharge)
66

67
    # get atomic symbol and coordinates of each atom
68 100
    xyz = oechem.OEFloatArray(3)
69 100
    for atom in mol.GetAtoms():
70 100
        mol.GetCoords(atom, xyz)
71 100
        inputstring+=( '\n  %s %10.4f %10.4f  %10.4f' \
72
                       %(oechem.OEGetAtomicSymbol(atom.GetAtomicNum()),
73
                       xyz[0], xyz[1], xyz[2]) )
74 100
    inputstring += ('\n  units angstrom\n}')
75

76
    # check if mol has a "freeze" tag
77 100
    for x in oechem.OEGetSDDataPairs(mol):
78 100
        if calctype == "opt" and "atoms to freeze" in x.GetTag():
79 100
            b = x.GetValue()
80 100
            y = b.replace("[", "")
81 100
            z = y.replace("]", "")
82 100
            a = z.replace(" ", "")
83 100
            freeze_list = a.split(",")
84 100
            inputstring += (
85
                "\n\nfreeze_list = \"\"\"\n  {} xyz\n  {} xyz\n  {} "
86
                "xyz\n  {} xyz\n\"\"\"".format(freeze_list[0], freeze_list[1],
87
                                               freeze_list[2], freeze_list[3]))
88 100
            inputstring += "\nset optking frozen_cartesian $freeze_list"
89 100
            inputstring += (
90
                "\nset optking dynamic_level = 1\nset optking "
91
                "consecutive_backsteps = 2\nset optking intrafrag_step_limit = "
92
                "0.1\nset optking interfrag_step_limit = 0.1\n")
93

94
    # best practices for scf calculations
95
    # http://www.psicode.org/psi4manual/master/scf.html#recommendations
96
    # http://www.psicode.org/psi4manual/master/dft.html#recommendations
97 100
    inputstring += '\n\nset scf_type df'
98 100
    inputstring += '\nset guess sad'
99

100
    # explicitly specify MP2 RI-auxiliary basis for [Ahlrichs] basis set
101
    # http://www.psicode.org/psi4manual/master/basissets_byfamily.html
102
    # DFMP2 *should* get MP2 aux sets fine for [Pople and Dunning] sets
103
    # http://www.psicode.org/psi4manual/master/dfmp2.html
104 100
    if method.lower() == 'mp2' and 'def2' in basisset:
105 100
        if basisset.lower() == 'def2-sv(p)':
106 100
            inputstring += ('\nset df_basis_mp2 def2-sv_p_-ri')
107 100
        elif basisset.lower() != 'def2-qzvpd':  # no aux set for qzvpd 10-6-18
108 0
            inputstring += ('\nset df_basis_mp2 %s-ri' % (basisset))
109

110
    # check for None or empty string such as in case of psi4 cbs energy
111 100
    if basisset:
112 100
        inputstring += ('\n\nset basis %s' % (basisset))
113

114 100
    inputstring += ('\nset freeze_core True')
115
    # specify command for type of calculation
116 100
    if calctype == 'opt':
117 100
        inputstring += ('\noptimize(\'%s\')\n\n' % (method))
118 100
    elif calctype == 'spe':
119 100
        if method[:3] == 'cbs':
120 0
            inputstring += ('\nenergy(%s)\n\n' % (method))
121
        else:
122 100
            inputstring += ('\nenergy(\'%s\')\n\n' % (method))
123 100
    elif calctype == 'hess':
124 100
        inputstring += (
125
            '\nH, wfn = hessian(\'%s\', return_wfn=True)\nwfn.hessian().print_out()\n\n'
126
            % (method))
127

128 100
    return inputstring
129

130

131 100
def make_psi_json(mol, label, method, basisset, calctype='opt', mem=None):
132
    """
133
    THIS FUNCTION IS A WORK IN PROGRESS.
134

135
    Get coordinates from input mol, and generate/format input text for
136
    Psi4 calculation via JSON wrapper.
137

138
    Parameters
139
    ----------
140
    mol : OpenEye OEMol
141
        OEMol with coordinates
142
    label : string
143
        Name of molecule with integer identifier (for conformers).
144
    method: string
145
        Name of the method as understood by Psi4. Example: "mp2"
146
    basis : string
147
        Name of the basis set as understood by Psi4. Example: "def2-sv(p)"
148
    calctype : string
149
        What kind of Psi4 calculation to run. Supported inputs are:
150
        'opt' for geometry optimization,
151
        'spe' for single point energy calculation, and
152
        'hess' for Hessian calculation.
153
    memory : string
154
        How much memory each Psi4 job should take. If not specified, the
155
        default in Psi4 is 500 Mb. Examples: "2000 MB" "1.5 GB"
156
        http://www.psicode.org/psi4manual/master/psithoninput.html
157

158
    Returns
159
    -------
160
    inputstring : string
161
        Contents of Psi4 input file to be written out
162

163
    """
164
    # check that specified calctype is valid
165 100
    check_calc(calctype)
166

167 100
    inputdict = {}
168 100
    moldict = {}
169 100
    modeldict = {}
170 100
    keydict = {}
171 100
    inputdict["schema_name"] = "qc_schema_input"
172 100
    inputdict["schema_version"] = 1
173

174
    # specify memory requirements, if defined
175 100
    if mem != None:
176 0
        inputdict["memory"] = mem
177

178
    #TODO -- json version
179
    # charge and multiplicity; multiplicity hardwired to singlet (usually is)
180
    #inputdict["charge"] = oechem.OENetCharge( mol)
181
    #inputdict["multiplicity"] = 1
182

183
    # get atomic symbol and coordinates of each atom
184 100
    geom_list = []
185 100
    elem_list = []
186 100
    xyz = oechem.OEFloatArray(3)
187 100
    for atom in mol.GetAtoms():
188 100
        mol.GetCoords(atom, xyz)
189 100
        geom_list.append(xyz[0])
190 100
        geom_list.append(xyz[1])
191 100
        geom_list.append(xyz[2])
192 100
        elem_list.append(oechem.OEGetAtomicSymbol(atom.GetAtomicNum()))
193 100
    moldict["geometry"] = geom_list
194 100
    moldict["symbols"] = elem_list
195 100
    inputdict["molecule"] = moldict
196

197
    #TODO -- json version
198
    ## check if mol has a "freeze" tag
199
    #for x in oechem.OEGetSDDataPairs(mol):
200
    #    if calctype=="opt" and "atoms to freeze" in x.GetTag():
201
    #        b = x.GetValue()
202
    #        y = b.replace("[", "")
203
    #        z = y.replace("]", "")
204
    #        a = z.replace(" ", "")
205
    #        freeze_list = a.split(",")
206
    #        inputstring += ("\n\nfreeze_list = \"\"\"\n  {} xyz\n  {} xyz\n  {} "
207
    #                       "xyz\n  {} xyz\n\"\"\"".format(freeze_list[0],
208
    #                       freeze_list[1], freeze_list[2], freeze_list[3]))
209
    #        inputstring += "\nset optking frozen_cartesian $freeze_list"
210
    #        inputstring += ("\nset optking dynamic_level = 1\nset optking "
211
    #            "consecutive_backsteps = 2\nset optking intrafrag_step_limit = "
212
    #            "0.1\nset optking interfrag_step_limit = 0.1\n")
213

214
    #TODO -- json version
215
    ## explicitly specify MP2 RI-auxiliary basis for Ahlrichs basis set
216
    ## http://www.psicode.org/psi4manual/master/basissets_byfamily.html
217
    ## DFMP2 *should* get MP2 aux sets fine for Pople/Dunning
218
    ## http://www.psicode.org/psi4manual/master/dfmp2.html
219
    #if method.lower()=='mp2' and 'def2' in basisset:
220
    #    if basisset.lower()=='def2-sv(p)':
221
    #        inputstring+=('\nset df_basis_mp2 def2-sv_p_-ri')
222
    #    elif basisset.lower()!='def2-qzvpd':  # no aux set for qzvpd 10-6-18
223
    #        inputstring+=('\nset df_basis_mp2 %s-ri' % (basisset))
224

225 100
    modeldict["basis"] = basisset
226 100
    modeldict["method"] = method
227 100
    inputdict["model"] = modeldict
228

229
    #inputstring+=('\nset freeze_core True')
230
    # specify command for type of calculation
231 100
    if calctype == 'opt':
232
        # TODO
233 0
        pass
234 100
    elif calctype == 'spe':
235 100
        inputdict["driver"] = 'energy'
236 0
    elif calctype == 'hess':
237 0
        inputdict["driver"] = 'hessian'
238 0
        keydict["return_wfn"] = True
239 100
    inputdict["keywords"] = keydict
240

241 100
    return inputdict
242

243

244 100
def confs_to_psi(insdf,
245
                 method,
246
                 basis,
247
                 calctype='opt',
248
                 memory=None,
249
                 via_json=False):
250
    """
251
    Read in molecule(s) (and conformers, if present) in insdf file. Create
252
    Psi4 input calculations for each structure.
253

254
    Parameters
255
    ----------
256
    insdf: string
257
        Name of the molecule file for which to create Psi4 input file.
258
        SDF format can contain multiple molecules and multiple conformers per
259
        molecule in a single file.
260
    method: string
261
        Name of the method as understood by Psi4. Example: "mp2"
262
    basis : string
263
        Name of the basis set as understood by Psi4. Example: "def2-sv(p)"
264
    calctype : string
265
        What kind of Psi4 calculation to run. Supported inputs are:
266
        'opt' for geometry optimization,
267
        'spe' for single point energy calculation, and
268
        'hess' for Hessian calculation.
269
    memory : string
270
        How much memory each Psi4 job should take. If not specified, the
271
        default in Psi4 is 500 Mb. Examples: "2000 MB" "1.5 GB"
272
        http://www.psicode.org/psi4manual/master/psithoninput.html
273
    via_json : Boolean
274
        If True, use JSON wrapper for Psi4 input and output.
275
        - Psi4 input would be in "input.py", called with python
276
        - Psi4 output would be in "output.json"
277
        If False, use normal text files for Psi4 input and output.
278
        - Psi4 input would be in "input.dat"
279
        - Psi4 output would be in "output.dat"
280
    """
281 100
    wdir = os.getcwd()
282

283
    # open molecules
284 100
    molecules = reader.read_mols(insdf)
285

286
    ### For each molecule: for each conf, generate input
287 100
    for mol in molecules:
288 100
        print(mol.GetTitle(), mol.NumConfs())
289 100
        if not mol.GetTitle():
290 0
            sys.exit("ERROR: OEMol must have title assigned! Exiting.")
291 100
        for i, conf in enumerate(mol.GetConfs()):
292
            # change into subdirectory ./mol/conf/
293 100
            subdir = os.path.join(wdir, "%s/%s" % (mol.GetTitle(), i + 1))
294 100
            if not os.path.isdir(subdir):
295 100
                os.makedirs(subdir)
296 100
            if os.path.exists(os.path.join(subdir, 'input.dat')):
297 0
                print("Input file already exists. Skipping.\n{}\n".format(
298
                    os.path.join(subdir, 'input.dat')))
299 0
                continue
300 100
            label = mol.GetTitle() + '_' + str(i + 1)
301 100
            if via_json:
302 100
                ofile = open(os.path.join(subdir, 'input.py'), 'w')
303 100
                ofile.write("# molecule {}\n\nimport numpy as np\nimport psi4"
304
                            "\nimport json\n\njson_data = ".format(label))
305 100
                json.dump(
306
                    make_psi_json(conf, label, method, basis, calctype,
307
                                  memory),
308
                    ofile,
309
                    indent=4,
310
                    separators=(',', ': '))
311 100
                ofile.write(
312
                    "\njson_ret = psi4.json_wrapper.run_json(json_data)\n\n")
313 100
                ofile.write("with open(\"output.json\", \"w\") as ofile:\n\t"
314
                            "json.dump(json_ret, ofile, indent=2)\n\n")
315
            else:
316 100
                ofile = open(os.path.join(subdir, 'input.dat'), 'w')
317 100
                ofile.write(
318
                    make_psi_input(conf, label, method, basis, calctype,
319
                                   memory))
320 100
            ofile.close()

Read our documentation on viewing source code .

Loading