1 4
__author__ = 'Chaya D. Stern'
2

3 4
import numpy as np
4 4
import itertools
5 4
from math import radians, degrees
6 4
import copy
7

8 4
from . import utils, chemi
9 4
from cmiles.utils import mol_to_smiles, has_atom_map, get_atom_map
10

11

12 4
def find_torsions(molecule, restricted=True, terminal=True):
13
    #ToDo: Get rid of equivalent torsions. Ex H-C-C-C and C-C-C-H.
14
    """
15
    This function takes an OEMol (atoms must be tagged with index map) and finds the map indices for torsion that need
16
    to be driven.
17

18
    Parameters
19
    ----------
20
    molecule : OEMol
21
        The atoms in the molecule need to be tagged with map indices
22
    restricted: bool, optional, default True
23
        If True, will find restricted torsions such as torsions in rings and double bonds.
24
    terminal: bool, optional, default True
25
        If True, will find terminal torsions
26

27
    Returns
28
    -------
29
    needed_torsion_scans: dict
30
        a dictionary that maps internal, terminal and restricted torsions to map indices of torsion atoms
31

32
    """
33
    # Check if molecule has map
34 4
    from openeye import oechem
35 4
    is_mapped = has_atom_map(molecule)
36 4
    if not is_mapped:
37 4
        utils.logger().warning('Molecule does not have atom map. A new map will be generated. You might need a new tagged SMARTS if the ordering was changed')
38 4
        tagged_smiles = mol_to_smiles(molecule, isomeric=True, mapped=True, explicit_hydrogen=True)
39
        # Generate new molecule with tags
40 4
        molecule = chemi.smiles_to_oemol(tagged_smiles)
41 4
        utils.logger().warning('If you already have a tagged SMARTS, compare it with the new one to ensure the ordering did not change')
42 4
        utils.logger().warning('The new tagged SMARTS is: {}'.format(tagged_smiles))
43
        # ToDo: save the new tagged SMILES somewhere. Maybe return it?
44

45 4
    needed_torsion_scans = {'internal': {}, 'terminal': {}, 'restricted': {}}
46 4
    mol = oechem.OEMol(molecule)
47 4
    if restricted:
48 4
        smarts = '[*]~[C,c]=,@[C,c]~[*]' # This should capture double bonds (not capturing rings because OpenEye does not
49
                                       # generate skewed conformations. ToDo: use scan in geometric or something else to get this done.
50 4
        restricted_tors = _find_torsions_from_smarts(molecule=mol, smarts=smarts)
51 4
        if len(restricted_tors) > 0:
52 0
            restricted_tors_min = one_torsion_per_rotatable_bond(restricted_tors)
53 0
            for i, tor in enumerate(restricted_tors_min):
54 0
                tor_name = ((tor[0].GetMapIdx() - 1), (tor[1].GetMapIdx() - 1), (tor[2].GetMapIdx() - 1), (tor[3].GetMapIdx() - 1))
55 0
                needed_torsion_scans['restricted']['torsion_{}'.format(str(i))] = tor_name
56

57 4
    if terminal:
58 4
        smarts = '[*]~[*]-[X2H1,X3H2,X4H3]-[#1]' # This smarts should match terminal torsions such as -CH3, -NH2, -NH3+, -OH, and -SH
59 4
        h_tors = _find_torsions_from_smarts(molecule=mol, smarts=smarts)
60 4
        if len(h_tors) > 0:
61 4
            h_tors_min = one_torsion_per_rotatable_bond(h_tors)
62 4
            for i, tor in enumerate(h_tors_min):
63 4
                tor_name = ((tor[0].GetMapIdx() -1 ), (tor[1].GetMapIdx() - 1), (tor[2].GetMapIdx() - 1), (tor[3].GetMapIdx() - 1))
64 4
                needed_torsion_scans['terminal']['torsion_{}'.format(str(i))] = tor_name
65

66 4
    mid_tors = [[tor.a, tor.b, tor.c, tor.d ] for tor in oechem.OEGetTorsions(mol)]
67 4
    if mid_tors:
68 4
        mid_tors_min = one_torsion_per_rotatable_bond(mid_tors)
69 4
        for i, tor in enumerate(mid_tors_min):
70 4
            tor_name = ((tor[0].GetMapIdx() - 1), (tor[1].GetMapIdx() - 1), (tor[2].GetMapIdx() - 1), (tor[3].GetMapIdx() - 1))
71 4
            needed_torsion_scans['internal']['torsion_{}'.format(str(i))] = tor_name
72

73
    # Check that there are no duplicate torsions in mid and h_torsions
74 4
    list_tor = list(needed_torsion_scans['internal'].values()) + list(needed_torsion_scans['terminal'].values())
75 4
    set_tor = set(list_tor)
76

77 4
    if not len(set_tor) == len(list_tor):
78 0
        raise Warning("There is a torsion defined in both mid and terminal torsions. This should not happen. Check "
79
                      "your molecule and the atom mapping")
80 4
    return needed_torsion_scans
81

82

83 4
def _find_torsions_from_smarts(molecule, smarts):
84
    """
85
    Do a substrcutre search on provided SMARTS to find torsions that match the SAMRTS
86

87
    Parameters
88
    ----------
89
    molecule: OEMol
90
        molecule to search on
91
    smarts: str
92
        SMARTS pattern to search for
93

94
    Returns
95
    -------
96
    tors: list
97
        list of torsions that match the SMARTS string
98

99
    """
100 4
    from openeye import oechem
101

102
    #ToDO use MDL aromaticity model
103 4
    qmol=oechem.OEQMol()
104 4
    if not oechem.OEParseSmarts(qmol, smarts):
105 0
        utils.logger().warning('OEParseSmarts failed')
106 4
    ss = oechem.OESubSearch(qmol)
107 4
    tors = []
108 4
    oechem.OEPrepareSearch(molecule, ss)
109 4
    unique = True
110 4
    for match in ss.Match(molecule, unique):
111 4
        tor = []
112 4
        for ma in match.GetAtoms():
113 4
            tor.append(ma.target)
114 4
        tors.append(tor)
115

116 4
    return tors
117

118

119 4
def one_torsion_per_rotatable_bond(torsion_list):
120
    """
121
    Keep only one torsion per rotatable bond
122
    Parameters
123
    ----------
124
    torsion_list: list
125
        list of torsion in molecule
126

127
    Returns
128
    -------
129
    list of only one torsion per rotatable bonds
130

131
    """
132

133 4
    central_bonds = np.zeros((len(torsion_list), 3), dtype=int)
134 4
    for i, tor in enumerate(torsion_list):
135 4
        central_bonds[i][0] = i
136 4
        central_bonds[i][1] = tor[1].GetIdx()
137 4
        central_bonds[i][2] = tor[2].GetIdx()
138

139 4
    grouped = central_bonds[central_bonds[:, 2].argsort()]
140 4
    sorted_tors = [torsion_list[i] for i in grouped[:, 0]]
141

142
    # Keep only one torsion per rotatable bond
143 4
    tors = []
144 4
    best_tor = [sorted_tors[0][0], sorted_tors[0][0], sorted_tors[0][0], sorted_tors[0][0]]
145 4
    best_tor_order = best_tor[0].GetAtomicNum() + best_tor[3].GetAtomicNum()
146 4
    first_pass = True
147 4
    for tor in sorted_tors:
148 4
        utils.logger().debug("Map Idxs: {} {} {} {}".format(tor[0].GetMapIdx(), tor[1].GetMapIdx(), tor[2].GetMapIdx(), tor[3].GetMapIdx()))
149 4
        utils.logger().debug("Atom Numbers: {} {} {} {}".format(tor[0].GetAtomicNum(), tor[1].GetAtomicNum(), tor[2].GetAtomicNum(), tor[3].GetAtomicNum()))
150 4
        if tor[1].GetMapIdx() != best_tor[1].GetMapIdx() or tor[2].GetMapIdx() != best_tor[2].GetMapIdx():
151
            #new_tor = True
152 4
            if not first_pass:
153 4
                utils.logger().debug("Adding to list: {} {} {} {}".format(best_tor[0].GetMapIdx(), best_tor[1].GetMapIdx(), best_tor[2].GetMapIdx(), best_tor[3].GetMapIdx()))
154 4
                tors.append(best_tor)
155 4
            first_pass = False
156 4
            best_tor = tor
157 4
            best_tor_order = tor[0].GetAtomicNum() + tor[3].GetAtomicNum()
158 4
            utils.logger().debug("new_tor with central bond across atoms: {} {}".format(tor[1].GetMapIdx(), tor[2].GetMapIdx()))
159
        else:
160 4
            utils.logger().debug("Not a new_tor but now with end atoms: {} {}".format(tor[0].GetMapIdx(), tor[3].GetMapIdx()))
161 4
            tor_order = tor[0].GetAtomicNum() + tor[3].GetAtomicNum()
162 4
            if tor_order > best_tor_order:
163 4
                best_tor = tor
164 4
                best_tor_order = tor_order
165 4
    utils.logger().debug("Adding to list: {} {} {} {}".format(best_tor[0].GetMapIdx(), best_tor[1].GetMapIdx(), best_tor[2].GetMapIdx(), best_tor[3].GetMapIdx()))
166 4
    tors.append(best_tor)
167

168 4
    utils.logger().debug("List of torsion to drive:")
169 4
    for tor in tors:
170 4
        utils.logger().debug("Idx: {} {} {} {}".format(tor[0].GetMapIdx(), tor[1].GetMapIdx(), tor[2].GetMapIdx(), tor[3].GetMapIdx()))
171 4
        utils.logger().debug("Atom numbers: {} {} {} {}".format(tor[0].GetAtomicNum(), tor[1].GetAtomicNum(), tor[2].GetAtomicNum(), tor[3].GetAtomicNum()))
172

173 4
    return tors
174

175 4
def find_torsion_around_bond(molecule, bond):
176
    """
177
    Find the torsion around a given bond
178
    Parameters
179
    ----------
180
    molecule : molecule with atom maps
181
    bond : tuple of map idx of bond atoms
182

183
    Returns
184
    -------
185
    list of 4 atom map idx (-1)
186

187
    Note:
188
    This returns the map indices of the torsion -1, not the atom indices.
189

190
    """
191 4
    from openeye import oechem
192
    #if not has_atom_map(molecule):
193
    #    raise ValueError("Molecule must have atom maps")
194

195 4
    terminal_smarts = '[*]~[*]-[X2H1,X3H2,X4H3]-[#1]'
196 4
    terminal_torsions = _find_torsions_from_smarts(molecule, terminal_smarts)
197 4
    mid_torsions = [[tor.a, tor.b, tor.c, tor.d] for tor in oechem.OEGetTorsions(molecule)]
198 4
    all_torsions = terminal_torsions + mid_torsions
199

200 4
    tors = one_torsion_per_rotatable_bond(all_torsions)
201

202 4
    tor_idx = [tuple(i.GetMapIdx() for i in tor) for tor in tors]
203 4
    central_bonds = [(tor[1], tor[2]) for tor in tor_idx]
204 4
    try:
205 4
        idx = central_bonds.index(bond)
206 0
    except ValueError:
207 0
        idx = central_bonds.index(tuple(reversed(bond)))
208

209 4
    torsion = [i-1 for i in tor_idx[idx]]
210 4
    return torsion
211

212

213 4
def find_equivelant_torsions(mapped_mol, restricted=False, central_bonds=None):
214
    """
215
    Final all torsions around a given central bond
216
    Parameters
217
    ----------
218
    mapped_mol: oemol. Must contaion map indices
219
    restricted: bool, optional, default False
220
        If True, will also find restricted torsions
221
    central_bonds: list of tuple of ints, optional, defualt None
222
        If provides, only torsions around those central bonds will be given. If None, all torsions in molecule will be found
223

224
    Returns
225
    -------
226
    eq_torsions: dict
227
        maps central bond to all equivelant torisons
228
    """
229
    #ToDo check that mol has mapping
230 4
    from openeye import oechem
231 4
    mol = oechem.OEMol(mapped_mol)
232 4
    if not has_atom_map(mol):
233 0
        raise ValueError("OEMol must have map indices")
234

235 4
    terminal_smarts = '[*]~[*]-[X2H1,X3H2,X4H3]-[#1]'
236 4
    terminal_torsions = _find_torsions_from_smarts(mol, terminal_smarts)
237 4
    mid_torsions = [[tor.a, tor.b, tor.c, tor.d] for tor in oechem.OEGetTorsions(mapped_mol)]
238 4
    all_torsions = terminal_torsions + mid_torsions
239 4
    if restricted:
240 0
        restricted_smarts = '[*]~[C,c]=,@[C,c]~[*]'
241 0
        restricted_torsions = _find_torsions_from_smarts(mol, restricted_smarts)
242 0
        all_torsions = all_torsions + restricted_torsions
243

244 4
    tor_idx = []
245 4
    for tor in all_torsions:
246 4
        tor_name = (tor[0].GetMapIdx()-1, tor[1].GetMapIdx()-1, tor[2].GetMapIdx()-1, tor[3].GetMapIdx()-1)
247 4
        tor_idx.append(tor_name)
248

249 4
    if central_bonds:
250 4
        if not isinstance(central_bonds, list):
251 4
            central_bonds = [central_bonds]
252 4
    if not central_bonds:
253 4
        central_bonds = set((tor[1], tor[2]) for tor in tor_idx)
254

255 4
    eq_torsions = {cb : [tor for tor in tor_idx if cb == (tor[1], tor[2]) or  cb ==(tor[2], tor[1])] for cb in
256
              central_bonds}
257

258 4
    return eq_torsions
259

260

261 4
def define_torsiondrive_jobs(needed_torsion_drives, internal_torsion_resolution=30, terminal_torsion_resolution=0,
262
                     scan_internal_terminal_combination=0, scan_dimension=2):
263
    """
264
    define crank jobs with torsions to drive and resolution to drive them at.
265

266
    Parameters
267
    ----------
268
    fragment_data: dict
269
        dictionary that maps fragment to needed torsions
270
    internal_torsion_resolution: int, optional. Default 15
271
        interval in degrees for torsion scan. If 0, internal torsions will not be driven
272
    terminal_torsion_resolution: int, optional. Default 0
273
        interval in degrees for torsion scans. If 0, terminal torsions will not be driven
274
    scan_internal_terminal_combination: int, optional. Default 0
275
        flag if internal and terminal torsions should be combined for higher dimension. If 0, only internal torsions will
276
        be driven. If 1, terminal and internal torsions will be scanned together.
277
    scan_dimension: int, optional. Default 2
278
        dimension of torsion scan. Combinations of torsions at the specified dimension will be generated as separate crank jobs
279
    qc_program: str, optional. Default Psi4
280
    method: str, optional. Default B3LYP
281
    basis: str, optional. Default aug-cc-pVDZ
282
    kwargs: optional keywords for psi4
283

284
    Returns
285
    -------
286
    fragment_data: dict
287
        dictionary that maps fragment to crank torsion jobs specifications.
288

289
    """
290

291 4
    if not internal_torsion_resolution and not terminal_torsion_resolution:
292 0
        utils.logger().warning("Resolution for internal and terminal torsions are 0. No torsions will be driven")
293

294 4
    if scan_internal_terminal_combination and (not internal_torsion_resolution or not terminal_torsion_resolution):
295 0
        raise Warning("If you are not scanning internal or terminal torsions, you must set scan_internal_terminal_"
296
                      "combinations to 0")
297

298 4
    internal_torsions = needed_torsion_drives['internal']
299 4
    terminal_torsions = needed_torsion_drives['terminal']
300 4
    internal_dimension = len(internal_torsions)
301 4
    terminal_dimension = len(terminal_torsions)
302 4
    torsion_dimension = internal_dimension + terminal_dimension
303

304 4
    crank_job = 0
305 4
    crank_jobs = dict()
306

307 4
    if not scan_internal_terminal_combination:
308 4
        if internal_torsion_resolution:
309 4
            for comb in itertools.combinations(internal_torsions, scan_dimension):
310 0
                dihedrals = [internal_torsions[torsion] for torsion in comb]
311 0
                grid = [internal_torsion_resolution]*len(dihedrals)
312 0
                crank_jobs['crank_job_{}'.format(crank_job)] = {'dihedrals': dihedrals, 'grid_spacing': grid}
313 0
                crank_job +=1
314 4
            if internal_dimension < scan_dimension and internal_dimension > 0:
315 4
                dihedrals = [internal_torsions[torsion] for torsion in internal_torsions]
316 4
                grid = [internal_torsion_resolution]*len(dihedrals)
317 4
                crank_jobs['crank_job_{}'.format(crank_job)] = {'dihedrals': dihedrals, 'grid_spacing': grid}
318 4
                crank_job +=1
319

320 4
        if terminal_torsion_resolution:
321
            # If scanning terminal torsions separately, only do 1D torsion scans
322 4
            for torsion in terminal_torsions:
323 4
                dihedrals = [terminal_torsions[torsion]]
324 4
                grid = [terminal_torsion_resolution]
325 4
                crank_jobs['crank_job_{}'.format(crank_job)] = {'dihedrals': dihedrals, 'grid_spacing': grid}
326 4
                crank_job +=1
327
            # for comb in itertools.combinations(terminal_torsions, scan_dimension):
328
            #     dihedrals = [terminal_torsions[torsion] for torsion in comb]
329
            #     grid = [terminal_torsion_resolution]*scan_dimension
330
            #     crank_jobs['crank_job_{}'.format(crank_job)] = {'dihedrals': dihedrals, 'grid_spacing': grid}
331
            #     crank_job +=1
332
            # if terminal_dimension < scan_dimension and terminal_dimension > 0:
333
            #     dihedrals = [terminal_torsions[torsion] for torsion in terminal_torsions]
334
            #     grid = [terminal_torsion_resolution]*len(dihedrals)
335
            #     crank_jobs['crank_job_{}'.format(crank_job)] = {'dihedrals': dihedrals, 'grid_spacing': grid}
336
            #     crank_job +=1
337
    else:
338
        # combine both internal and terminal torsions
339 0
        all_torsion_idx = np.arange(0, torsion_dimension)
340 0
        for comb in itertools.combinations(all_torsion_idx, scan_dimension):
341 0
            internal_torsions = [internal_torsions['torsion_{}'.format(i)] for i in comb if i < internal_dimension]
342 0
            terminal_torsions = [terminal_torsions['torsion_{}'.format(i-internal_dimension)] for i in comb if i >= internal_dimension]
343 0
            grid = [internal_torsion_resolution]*len(internal_torsions)
344 0
            grid.extend([terminal_torsion_resolution]*len(terminal_torsions))
345 0
            dihedrals = internal_torsions + terminal_torsions
346 0
            crank_jobs['crank_job_{}'.format(crank_job)] = {'diherals': dihedrals, 'grid_spacing': grid}
347 0
            crank_job += 1
348 0
        if torsion_dimension < scan_dimension:
349 0
            internal_torsions = [internal_torsions['torsion_{}'.format(i)] for i in all_torsion_idx if i < internal_dimension]
350 0
            terminal_torsions = [terminal_torsions['torsion_{}'.format(i-internal_dimension)] for i in all_torsion_idx if i >= internal_dimension]
351 0
            grid = [internal_torsion_resolution]*len(internal_torsions)
352 0
            grid.extend([terminal_torsion_resolution]*len(terminal_torsions))
353 0
            dihedrals = internal_torsions + terminal_torsions
354 0
            crank_jobs['crank_job_{}'.format(crank_job)] = {'diherals': dihedrals, 'grid_spacing': grid}
355 0
            crank_job += 1
356

357 4
    return crank_jobs
358

359

360 4
def define_restricted_drive(qc_molecule, restricted_dihedrals, steps=6, maximum_rotation=30, scan_dimension=1):
361
    """
362

363
    Parameters
364
    ----------
365
    qc_molecule: molecule in QC_JSON format. This comes with a geometry, connectivity table, identifiers that also has
366
    a mapped SMILES so it can be checked.
367
    needed_torsion_drives
368
    grid_resolution
369
    maximum_rotation
370

371
    Returns
372
    -------
373

374
    """
375
    #ToDo extend to multi-dimensional scans
376
    #natoms = len(qc_molecule['symbols'])
377
    # Convert to 3D shape for
378
    #coords = np.array(qc_molecule['geometry'], dtype=float).reshape(natoms, 3) * utils.BOHR_2_ANGSTROM
379 0
    connectivity = np.asarray(qc_molecule['connectivity'])
380

381
    # Check dihedral indices are connected
382 0
    bond_tuples = list(zip(connectivity[:, :2].T[0], connectivity[:, :2].T[1]))
383 0
    optimization_jobs = {}
384 0
    i = 0
385 0
    for torsion in restricted_dihedrals:
386 0
        for a1, a2 in zip(restricted_dihedrals[torsion], restricted_dihedrals[torsion][1:]):
387 0
            if (a1, a2) not in bond_tuples and (a2, a1) not in bond_tuples:
388 0
                utils.logger().warning("torsion {} is not bonded. Skipping this torsion")
389 0
                continue
390
        # measure dihedral angle
391 0
        dihedral_angle = measure_dihedral_angle(restricted_dihedrals[torsion], qc_molecule['geometry'])
392 0
        t_tuple = restricted_dihedrals[torsion]
393 0
        angle = round(dihedral_angle)
394 0
        optimization_jobs['{}_{}'.format(t_tuple, i)] = {
395
            'type': 'optimization_input',
396
            'initial_molecule': qc_molecule,
397
            'dihedrals': [restricted_dihedrals[torsion]],
398
            'constraints': {'scan': [('dihedral', str(t_tuple[0]), str(t_tuple[1]), str(t_tuple[2]),
399
                                               str(t_tuple[3]), str(angle), str(angle + maximum_rotation),
400
                                               str(steps))]}}
401 0
        optimization_jobs['{}_{}'.format(t_tuple, i+1)] = {
402
            'type': 'optimization_input',
403
            'initial_molecule': qc_molecule,
404
            'dihedrals': [restricted_dihedrals[torsion]],
405
            'constraints': {'scan': [('dihedral', str(t_tuple[0]), str(t_tuple[1]), str(t_tuple[2]),
406
                                               str(t_tuple[3]), str(angle), str(angle - maximum_rotation),
407
                                               str(steps))]}}
408

409 0
    return optimization_jobs
410

411

412 4
def generate_constraint_opt_input(qc_molecule, dihedrals, maximum_rotation=30, interval=5, filename=None):
413
    """
414

415
    Parameters
416
    ----------
417
    qc_molecule
418
    dihedrals
419

420
    Returns
421
    -------
422
    QCFractal optimization jobs input
423

424
    """
425 0
    from openeye import oechem
426

427 0
    optimization_jobs = {}
428 0
    tagged_smiles = qc_molecule['identifiers']['canonical_isomeric_explicit_hydrogen_mapped_smiles']
429 0
    mol = oechem.OEMol()
430 0
    oechem.OESmilesToMol(mol, tagged_smiles)
431 0
    atom_map = get_atom_map(mol, tagged_smiles)
432

433 0
    coords = chemi.from_mapped_xyz_to_mol_idx_order(qc_molecule['geometry'], atom_map)
434

435
    # convert coord to Angstrom
436 0
    coords = coords * utils.BOHR_2_ANGSTROM
437 0
    conf = mol.GetConfs().next()
438 0
    conf.SetCoords(oechem.OEFloatArray(coords))
439

440
    # new molecule for setting dihedral angles
441 0
    mol_2 = oechem.OEMol(mol)
442 0
    conf_2 = mol_2.GetConfs().next()
443 0
    coords_2 = oechem.OEFloatArray(conf_2.GetMaxAtomIdx()*3)
444 0
    conf.GetCoords(coords_2)
445 0
    mol_2.DeleteConfs()
446

447 0
    interval = radians(interval)
448 0
    max_rot = radians(maximum_rotation)
449 0
    for dihedral in dihedrals:
450
        #j = 0
451 0
        dih_idx = dihedrals[dihedral]
452 0
        tor = []
453 0
        for i in dih_idx:
454 0
            a = mol.GetAtom(oechem.OEHasMapIdx(i+1))
455 0
            tor.append(a)
456 0
        dih_angle = oechem.OEGetTorsion(conf, tor[0], tor[1], tor[2], tor[3])
457 0
        for i, angle in enumerate(np.arange(dih_angle-max_rot, dih_angle+max_rot, interval)):
458 0
            newconf = mol.NewConf(coords_2)
459 0
            oechem.OESetTorsion(newconf, tor[0], tor[1], tor[2], tor[3], angle)
460
            #new_angle = oechem.OEGetTorsion(newconf, tor[0], tor[1], tor[2], tor[3])
461
            # if new_angle == dih_angle:
462
            #     j += 1
463
            #     if j > 1:
464
            #         # One equivalent angle should be generated.
465
            #         logger().warning("Openeye did not generate a new conformer for torsion and angle {} {}. Will not generate"
466
            #                      "qcfractal optimizaiton input".format(dih_idx, angle))
467
            #         break
468 0
            if filename:
469 0
                pdb = oechem.oemolostream("{}_{}.pdb".format(filename, i))
470 0
                oechem.OEWritePDBFile(pdb, newconf)
471 0
            symbols, geometry = chemi.to_mapped_geometry(newconf, atom_map)
472 0
            qc_molecule = copy.deepcopy(qc_molecule)
473 0
            qc_molecule['geometry'] = geometry
474 0
            qc_molecule['symbols'] = symbols
475 0
            degree = degrees(angle)
476 0
            optimization_jobs['{}_{}'.format(dih_idx, int(round(degree)))] = {
477
                'type': 'optimization_input',
478
                'initial_molecule': qc_molecule,
479
                'dihedral': dih_idx,
480
                'constraints': {
481
                    "set": [{
482
                        "type": "dihedral",
483
                        "indices": dih_idx,
484
                        "value": degree
485
                    }]
486
                }
487
            }
488 0
    return optimization_jobs
489

490

491 4
def measure_dihedral_angle(dihedral, coords):
492
    """
493
    calculate the dihedral angle in degrees
494

495
    Parameters
496
    ----------
497
    dihedral
498
    coords
499

500
    Returns
501
    -------
502

503
    """
504 0
    reshape = True
505 0
    if isinstance(coords, np.ndarray):
506 0
        if coords.shape[-1] == 3:
507 0
            reshape = False
508 0
    if reshape:
509 0
        coords = np.array(coords, dtype=float).reshape((len(coords)/3), 3) * utils.BOHR_2_ANGSTROM
510 0
    a = coords[dihedral[0]]
511 0
    b = coords[dihedral[1]]
512 0
    c = coords[dihedral[2]]
513 0
    d = coords[dihedral[3]]
514 0
    v1 = b-a
515 0
    v2 = c-b
516 0
    v3 = d-c
517 0
    t1 = np.linalg.norm(v2)*np.dot(v1, np.cross(v2, v3))
518 0
    t2 = np.dot(np.cross(v1, v2), np.cross(v2, v3))
519 0
    phi = np.arctan2(t1, t2)
520 0
    degree = phi * 180 / np.pi
521 0
    return degree

Read our documentation on viewing source code .

Loading