diff --git a/ForceBalance_12mol_plots/allData/process_QMMM.py b/ForceBalance_12mol_plots/allData/process_QMMM.py new file mode 100644 index 0000000..9332fef --- /dev/null +++ b/ForceBalance_12mol_plots/allData/process_QMMM.py @@ -0,0 +1,480 @@ +### This script will take an .xyz (scan-final.xyz) from the finished QM torsion scan and generate/process data to create a plot of QM and MM energies. +#Imports +from openforcefield.typing.engines.smirnoff import * +from openforcefield.utils import get_data_filename, extractPositionsFromOEMol, generateTopologyFromOEMol +from openeye.oechem import * +from openeye.oeomega import * # conformer generation +from openeye.oequacpac import * #for partial charge assignment +from calc_improper import find_improper_angles, calc_improper_angles, angle_betwen +import numpy as np +import matplotlib.pyplot as plt + +#Functions +def SDF2oemol(sdffile): + """ + Description: + Takes in the sdf file and converts it to an oemol + + Input: + sdffile: .sdf file + + Returns: + Oemol: Oemol object of the molecule stored in the input .sdf file + """ + #open file + ifs = oechem.oemolistream() + oemol = oechem.OECreateOEGraphMol() + ifs.open(sdffile) + #empty oemol list + molecules = list() + oechem.OEReadMolecule(ifs, oemol) + ifs.close() + return oemol + +def smiles2oemol(smiles): + """ + Description: + Takes in an SMILES string and returns an oemol. + + Input: + smiles: SMILES string for molecule + + Returns: + mol1: OEMol object + """ + mol1 = OEMol() + OESmilesToMol(mol1, smiles) + + #assign charges, necessary to keep? + chargeEngine = OEAM1BCCCharges() + OEAssignCharges(mol1, chargeEngine) + + return mol1 + + +def QM2Oemol(xyzfile, oemol): + """ + Description: + Takes in an xyz file and creates oemol objects with the coordinates in the original .xyz file + with QM energies. The .xyz files are formatted as an output geomeTRIC xyz file, scan-final.xyz. The original energies are stored in Hartree but are converted to kcal/mol in the QM energy tag. + Input: + xyzfile:Take in scan-final.xyz file which contains the geometry and energy outputs from a + geomeTRIC torsion scan + oemol:The oemol of the molecule involved in the torsion scan of the .xyz file + + Returns: + oemolList: A list of OEmols that contain the QM data tagged as "QMEng" + """ + + #open input xyz file + ifs = oechem.oemolistream() + ifs.open(xyzfile) + + #generate empty list of oemols + oemolList=list() + + #iterate through oemols and set coordinates to original oemol + for mol in ifs.GetOEGraphMols(): + #get coordinates + coords = mol.GetCoords() + #set coordinates of original oemol with correct bond connectivity + newmol = copy.deepcopy(oemol) + newmol.SetCoords(coords) + #get the energies from the .xyz file title and save in oemol data + title = mol.GetTitle() + energy = float(str.split(title)[-1]) + energy_kcal= energy * 627.509 + #set the QM energy in a tag called qm + newmol.SetData("QM", energy_kcal) + + #append the list with the new mol that has stored QM energies + oemolList.append(newmol) + + return oemolList + + +def GetMM(oemol, FF, tag): + """ + Description: + Takes in an oemol and calculates the MM energies with two forcefields, one which has removed + nitrogen improper parameters. The data from these calculations is stored in the oemol object + as MM and MMRmNit. The units of energy are KCal/mol. + + Input: + oemol: A single oemol object + FF: .offxml file of smirnoff99Frosst.offxml + tag: The name to tag the energy as, ex: "MM" + + Return: + oemol: An oemol with the MM energies stored in tag in units kcal/mol + """ + + #prep both force fields, create system and topology, get MM energies and store in data tag + ff = ForceField(FF) + topology = generateTopologyFromOEMol(oemol) + system = ff.createSystem(topology, [oemol]) + positions = extractPositionsFromOEMol(oemol) + integrator = openmm.VerletIntegrator(2.0*unit.femtoseconds) + simulation = app.Simulation(topology, system, integrator) + simulation.context.setPositions(positions) + state = simulation.context.getState(getEnergy = True, getPositions=True) + energy = state.getPotentialEnergy() / unit.kilocalories_per_mole + oemol.SetData(tag, energy) + + return oemol + +def findAngles(oemol, constraint): + """ + Description: + Calculates the improper (and potentially valence) angles of an oemol and stores the value in + tag "improper" and "valence" for an oemol. + + Input: + oemol: An oemol object + constraint: An constraints.txt file used in an input for an geomeTRIC scan + + Return: + oemol: An oemol with the calculated improper or valence angles with corresponding tags + "improper" or "valence" + """ + + #determine the constraints + #open the constraint file + constraintfile = open(constraint, "r") + f = open(constraint) + lines = f.readlines() + f.close() + coords = oemol.GetCoords() + + for l in lines: + if "dihedral" in l: + split = l.split() + atoms_imp = (int(split[1])-1, int(split[2])-1, int(split[3])-1, int(split[4])-1) + crd1 = np.asarray(coords[atoms_imp[0]]) + crd2 = np.asarray(coords[atoms_imp[1]]) + crd3 = np.asarray(coords[atoms_imp[2]]) + crd4 = np.asarray(coords[atoms_imp[3]]) + angle_imp = calc_improper_angle(crd1, crd2, crd3, crd4, True) + if angle_imp < 90 and angle_imp > 0: + oemol.SetData("improper", angle_imp) + if angle_imp < 0: + if angle_imp < -90: + #angle_imp = 180 + angle_imp + oemol.SetData("improper", angle_imp) + if angle_imp > 90: + angle_imp = angle_imp - 180 + oemol.SetData("improper", angle_imp) + + #calculate the valence angle and store in tag (2-d scan) + if "angle" in l: + split = l.split() + atoms_val = (int(split[1])-1, int(split[2])-1, int(split[3])-1) + crd1 = np.asarray(coords[atoms_val[0]]) + crd2 = np.asarray(coords[atoms_val[1]]) + crd3 = np.asarray(coords[atoms_val[2]]) + angle_val= calc_valence_angle(crd1, crd2, crd3) + oemol.SetData("valence", angle_val) + + return oemol + + + +#TODO +def makeOEB(oemolList, tag): + """ + Description: + Takes in an oemol list and creates an output OEB file. + + Input: + oemolList: A list of oemols + tag: The title of the OEB file. + + Return: + + """ + ofile = oemolostream(tag+'.oeb') + for mol in oemolList: + OEWriteConstMolecule(ofile, mol) + ofile.close() + return + + + + +def adjust_energy(oemolList, qm_tag, mm_tag): + """ + Description: + Iterates through oemols in a list and normalizes the energies to the lowest + QM energy in the list. The corresponding MM energy geometry is subtracted from the + MM energies tagged in the oemol. + + Input: + oemolList: A list of oemols + qm_tag: The tag for the qm energy + mm_tag: The tag for the mm energy being normalized + + Return: + oemolList with the updated energies. + """ + + #sort the oemols based on the lowest QM energy, this oemol is now "low_mol" + low_mol = sorted(oemolList, key=lambda x: x.GetData(qm_tag))[0] + #get corresponding lowest qm energy + low_qm = low_mol.GetData(qm_tag) + + #get correspoding lowest mm energy + #low_mol_mm = sorted(oemolList, key=lambda x: x.GetData(mm_tag))[0] + low_mm = low_mol.GetData(mm_tag) + + #iterate through the list and subtract lowest mm and qm energies + for m in oemolList: + qm = m.GetData(qm_tag) - low_qm + m.SetData(qm_tag, qm) + mm = m.GetData(mm_tag) - low_mm + m.SetData(mm_tag, mm) + + return oemolList + + + + +def plotResultsMulti(oemolList, names, tagList, colorList, markers): + """ + Description: This function plots data for multiple molecules in a tag list. + + Input: + oemolList: List of oemols with tags of energy data + names: List of molecule names to title plot + tagList: List of tagged data of interest + colorList: List of colors for plot corresponding with the tags + markers: Specified marker types for the plot + """ + for mols, title in zip(oemolList, names): + sort_oemol = sorted(mols, key=lambda x: x.GetData('improper')) + xs = [m.GetData('improper') for m in sort_oemol] + qm = [m.GetData('QM') for m in sort_oemol] + plt.plot(xs, qm, color="royalblue", label='QM', marker= "^", markersize="8", linewidth="3") + for tag, col, mar in zip(tagList, colorList, markers): + trend = str(tag) + name = [m.GetData(tag) for m in sort_oemol] + plt.plot(xs, name, color=col, label=trend, marker=mar,markersize="8", linewidth="3") + plt.title(title) + plt.legend() + plt.show() + + + +def oeb2mollist(oeb): + """ + Description: + Takes in oeb file and creates oemolList + + Input: + oeb: oeb file + + Return: + oemolList: a list of eomols contained in the .oeb file + """ + + #open input xyz file + ifs = oechem.oemolistream() + ifs.open(oeb) + + #generate empty list of oemols + oemolList=list() + + #iterate through oemols and set coordinates to original oemol + for mol in ifs.GetOEGraphMols(): + oemolList.append(oechem.OEGraphMol(mol)) + + return oemolList + + + +def compareGromacs(oemolList, FF, topfile, grofile): + """ + Description: + Compares energies of openMM to gromacs + + input: + oemolList: List of oemol objects to compare + FF: openmm force field, .offxml file + topfile: GROMACS topology file name to write + grofile: GROMACS coordinate file name (.gro format) to write + + + """ + + for idx, mol in enumerate(oemolList): + name = str(idx) + "_" + topfile + name_gro = str(idx) + "_" + grofile + ff = ForceField(FF) + topology = generateTopologyFromOEMol(mol) + system = ff.createSystem(topology, [mol]) + positions = extractPositionsFromOEMol(mol) + save_system_to_gromacs(topology, system, positions, name, name_gro) + + + top = parmed.load_file(name) + gromacssys = top.createSystem(nonbondedMethod= app.NoCutoff, constraints = None, implicitSolvent = None) + gro = parmed.load_file(name_gro) + + #smirnoff + smirfftop, smirffsys, smirffpos = create_system_from_molecule(ff, mol, verbose = False) + + print(oechem.OEMolToSmiles(mol)) + #compare + try: + groups0, groups1, energy0, energy1 = compare_system_energies( smirfftop, top.topology, smirffsys, gromacssys, positions, verbose = False) + print(groups0, groups1, energy0, energy1) + except: + print(oechem.OEMolToSmiles(mol)) + print("Failed") + + + + + +def sdf2mol2(sdf, output): + """ + Takes in .sdf file and writes to a .mol2 file + input: + sdf: .sdf file + output: name of output .mol2 file + """ + + ifs = oechem.oemolistream(sdf) + ofs = oechem.oemolostream(output) + + ifs.SetFormat(oechem.OEFormat_SDF) + ofs.SetFormat(oechem.OEFormat_MOL2) + + for mol in ifs.GetOEGraphMols(): + oechem.OEWriteMolecule(ofs, mol) + + + +def oeb2oemol(oebfile): + """ + Takes in oebfile and generates oemolList + input: + oebfile: Oebfile (.oeb) + + return: + mollist: List of oemols + + """ + + ifs = oechem.oemolistream(oebfile) + mollist = [] + + for mol in ifs.GetOEGraphMols(): + mollist.append(oechem.OEGraphMol(mol)) + + return mollist + + +def processData(smiles, xyzfile, constraintFile, moltitle, FF, tags, color): + """ + Description: + Takes in innitial .sdf file, xyzfile from geomeTRIC output, constraint file and a title to + create an .oeb file with oemols with QM, and MM data. + + + input: + #UPDATE to .sdf or .mol2 file + smiles: Smiles string for molecule + xyzfile: output .xyz file from geomeTRIC scan, scan-final.xyz that contains QM energies in title + constraintFile: constraint.txt for geomeTRIC input that contains the indicies constrained in scan + molTitle: title for the output .oeb file + FF: Forcefield list + tags: List of tags corresponding to force field + color: List of colors used for plot + + Return: + none, generates .oeb file and plot of data + """ + + #generate lsit of oemols with stored QM energies + oemolList = QM2Oemol(xyzfile, SDF2oemol(smiles)) + + #get MM energies from force fields + for f, tag in zip(FF,tags): + for mol in oemolList: + print("THIS IS THE FORCE FIELD" + str(f)) + GetMM(mol, f, tag) + findAngles(mol, constraintFile) + + #normalize eneriges: + for tag in tags: + adjust_energy(oemolList, 'QM', tag) + + #write out oeb file + makeOEB(oemolList, moltitle) + + return oemolList + + +def makeInputs(molName): + ''' + This function makes the input files in a format to process MM and QM energies of groups of molecules + + Input: + molName: List of molecules name in a set of molecules + + Return + molName: Original molName list + SDFList: List of sdf files for molecule in list + contList: List of constraint files for molecule + XYZList: List of .xyz scan files + ''' + SDFList = [] + contList = [] + XYZList = [] + + for name in molName: + sdfName = str(name) + '.sdf' + SDFList.append(sdfName) + + contName = 'c-' + str(name) + '.txt' + contList.append(contName) + + XYZName = 'scan-final-' + str(name) + '.xyz' + XYZList.append(XYZName) + + return molName, SDFList, contList, XYZList + + +def processSet(filesList, molNames, FFList, tagNames, colorList, markerList): + ''' + This function takes in information for a set of molecules and generates + subplots for all of the molecules from MM energies from the specified FFs. + + input: + filesList: A list of molNames, SDF file of starting molecule, constraints file list, xyz files with QM scans + molNames: List of oemol names + FFList : List of FFs of interest + tagNames: name of tags associated with FFs + colorList: List of colors for the plot + markerList: Type of markers for pots + ''' + + names = filesList[0] + sdf = filesList[1] + cont = filesList[2] + xyz = filesList[3] + + mols = [] + + i = 0 + while i < len(filesList[0]): + doneList = processData(sdf[i], xyz[i], cont[i], names[i], FFList, tagNames, colorList) + mols.append(doneList) + i += 1 + + #plot the list of list of oemols + plotResultsMulti(mols, molNames, tagNames, colorList, markerList) + + + + diff --git a/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/calc_improper.py b/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/calc_improper.py new file mode 100644 index 0000000..04ffee2 --- /dev/null +++ b/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/calc_improper.py @@ -0,0 +1,173 @@ +#!/usr/bin/env python + +#============================================================================================= +# MODULE DOCSTRING +#============================================================================================= + +""" +calc_improper.py + +Find and calculate improper dihedral angles in a given molecule. +Code loosely follows OpenMM: https://tinyurl.com/y8mhwxlv + +Let's say we have j as central atom and call addTorsion(j, i, k, l). +Then we compute the following vectors: + + v0 = j-i + v1 = k-i + v2 = k-l + w0 = v0 x v1 + w1 = v1 x v2 + +The final improper angle is computed as the angle between w0 and w1. + +By: Victoria Lim + +""" + +#============================================================================================= +# GLOBAL IMPORTS +#============================================================================================= + +import numpy as np +import openeye.oechem as oechem + +#============================================================================================= +# PRIVATE SUBROUTINES +#============================================================================================= + +def angle_between(v1, v2): + """ + Calculate the angle in degrees between vectors 'v1' and 'v2'. + Modified from: https://tinyurl.com/yb89sstz + + Parameters + ---------- + v1 : tuple, list, or numpy array + v2 : tuple, list, or numpy array + + Returns + ------- + float + Angle in degrees. + + """ + v1_u = v1/np.linalg.norm(v1) + v2_u = v2/np.linalg.norm(v2) + return np.degrees(np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))) + +def calc_valence_angle(atom0, atom1, atom2): + """ + Calculate the valence angle of three atoms. + + Parameters + ---------- + atom0 : numpy array + CENTRAL atom coordinates + atom1 : numpy array + outer atom coordinates + atom2 : numpy array + outer atom coordinates + + Returns + ------- + float + Angle in degrees. + """ + v1 = atom1-atom0 + v2 = atom2-atom0 + return(angle_between(v1, v2)) + +def calc_improper_angle(atom0, atom1, atom2, atom3, translate=False): + """ + Calculate the improper dihedral angle of a set of given four atoms. + + Parameters + ---------- + atom0 : numpy array + CENTRAL atom coordinates + atom1 : numpy array + outer atom coordinates + atom2 : numpy array + outer atom coordinates + atom3 : numpy array + outer atom coordinates + translate : bool + True to translate central atom to origin, False to keep as is. + This should not affect the results of the calculation. + + Returns + ------- + float + Angle in degrees. + """ + if translate: + atom1 = atom1 - atom0 + atom2 = atom2 - atom0 + atom3 = atom3 - atom0 + atom0 = atom0 - atom0 # central must be moved last + + # calculate vectors + v0 = atom0-atom1 + v1 = atom2-atom1 + v2 = atom2-atom3 + w1 = np.cross(v0, v1) + w2 = np.cross(v1, v2) + angle = angle_between(w1,w2) # this angle should be in range [0,90] + + # compute distance from plane to central atom + # eq 6 from http://mathworld.wolfram.com/Point-PlaneDistance.html + # here I'm using atom1 for (x,y,z), but could also use atom2 or atom3 + numer = w2[0]*(atom0[0]-atom1[0]) + w2[1]*(atom0[1]-atom1[1]) + w2[2]*(atom0[2]-atom1[2]) + denom = np.sqrt(w2[0]**2 + w2[1]**2 + w2[2]**2) + dist = numer/denom + # set reference so that if central atom is above plane, angle -> [90,180] + if dist > 0: + angle = 180-angle + + return angle + +def find_improper_angles(mol): + """ + Find the improper dihedral angles in some molecule. Currently supports + those with a central trivalent nitrogen atom. + + Parameters + ---------- + mol : OpenEye oemol + oemol in which to look for improper angles + + Returns + ------- + list + Each element in the list is a 4-tuple of the coordinates for the + atoms involved in the improper. The central atom is listed first + in the tuple. Each member of the tuple is a numpy array. + list + List of strings for atoms in the improper, central atom is first. + """ + + mol_coords = mol.GetCoords() + crdlist = [] + Idxlist = [] + for atom in mol.GetAtoms(oechem.OEIsInvertibleNitrogen()): + # central atom + aidx = atom.GetIdx() + crd0 = np.asarray(mol_coords[aidx]) + # sort the neighbors + nbors = sorted(list(atom.GetAtoms())) + #check if there are 3 atoms connected to central atom in improper + if len(nbors) != 3: + return crdlist, namelist + crd1 = np.asarray(mol_coords[nbors[0].GetIdx()]) + crd2 = np.asarray(mol_coords[nbors[1].GetIdx()]) + crd3 = np.asarray(mol_coords[nbors[2].GetIdx()]) + # store coordinates + crdlist.append([crd0, crd1, crd2, crd3]) + Idxlist.append([atom.GetIdx(), nbors[0].GetIdx(), nbors[1].GetIdx(),nbors[2].GetIdx()]) + + + return crdlist, Idxlist + + + diff --git a/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/process_QMMM.py b/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/process_QMMM.py new file mode 100644 index 0000000..0f50168 --- /dev/null +++ b/mmCalc/mol2_geometric/mol6_scan_40_40/Data_Processing/process_QMMM.py @@ -0,0 +1,392 @@ +###WIP: +### This script will take an .xyz (scan-final.xyz) from the finished QM torsion scan and generate/process data to create a plot of QM and MM energies. + + +#Imports +from openforcefield.typing.engines.smirnoff import * +from openforcefield.utils import get_data_filename, extractPositionsFromOEMol, generateTopologyFromOEMol +from openeye.oechem import * +#import oenotebook as oenb +from openeye.oeomega import * # conformer generation +from openeye.oequacpac import * #for partial charge assignment +from calc_improper import * +import numpy as np +import matplotlib.pyplot as plt + + +#Functions +def SDF2oemol(sdffile): + """ + Description: + Takes in the sdf file and converts it to an oemol + + Input: + sdffile: .sdf file + + Returns: + Oemol: Oemol object of the molecule stored in the input .sdf file + """ + #open file + ifs = oechem.oemolistream() + oemol = oechem.OECreateOEGraphMol() + ifs.open(sdffile) + #empty oemol list + molecules = list() + oechem.OEReadMolecule(ifs, oemol) + ifs.close() + print(oemol.GetCoords()) + return oemol + +def smiles2oemol(smiles): + """ + Description: + Takes in an SMILES string and returns an oemol. + + Input: + smiles: SMILES string for molecule + + Returns: + mol1: OEMol object + """ + mol1 = OEMol() + OESmilesToMol(mol1, smiles) + + #assign charges, necessary to keep? + chargeEngine = OEAM1BCCCharges() + OEAssignCharges(mol1, chargeEngine) + + return mol1 + + +def QM2Oemol(xyzfile, oemol): + """ + Description: + Takes in an xyz file and creates oemol objects with the coordinates in the original .xyz file + with QM energies. The .xyz files are formatted as an output geomeTRIC xyz file, scan-final.xyz. The original energies are stored in Hartree but are converted to kcal/mol in the QM energy tag. + Input: + xyzfile:Take in scan-final.xyz file which contains the geometry and energy outputs from a + geomeTRIC torsion scan + oemol:The oemol of the molecule involved in the torsion scan of the .xyz file + + Returns: + oemolList: A list of OEmols that contain the QM data tagged as "QMEng" + """ + + #open input xyz file + ifs = oechem.oemolistream() + ifs.open(xyzfile) + + #generate empty list of oemols + oemolList=list() + + #iterate through oemols and set coordinates to original oemol + for mol in ifs.GetOEGraphMols(): + #get coordinates + coords = mol.GetCoords() + #set coordinates of original oemol with correct bond connectivity + newmol = copy.deepcopy(oemol) + newmol.SetCoords(coords) + #get the energies from the .xyz file title and save in oemol data + title = mol.GetTitle() + energy = float(str.split(title)[-1]) + energy_kcal= energy * 627.509 + #set the QM energy in a tag called qm + newmol.SetData("QM", energy_kcal) + + #append the list with the new mol that has stored QM energies + oemolList.append(newmol) + + return oemolList + + +def GetMM(oemol, FF, FFRmNit): + """ + Description: + Takes in an oemol and calculates the MM energies with two forcefields, one which has removed + nitrogen improper parameters. The data from these calculations is stored in the oemol object + as MM and MMRmNit. The units of energy are KCal/mol. + + Input: + oemol: A single oemol object + FF: .offxml file of smirnoff99Frosst.offxml + FFRmNit: .offxml file of smirnoff99Frosst.offxml with the removed nitrogen improper parameter + + Return: + oemol: An oemol with the MM energies stored in tags MM and MMRmNit in kcal/mol + """ + + #prep both force fields, create system and topology, get MM energies and store in data tag + ff = ForceField(FF) + topology = generateTopologyFromOEMol(oemol) + system = ff.createSystem(topology, [oemol]) + positions = extractPositionsFromOEMol(oemol) + integrator = openmm.VerletIntegrator(2.0*unit.femtoseconds) + simulation = app.Simulation(topology, system, integrator) + simulation.context.setPositions(positions) + state = simulation.context.getState(getEnergy = True, getPositions=True) + energy = state.getPotentialEnergy() / unit.kilocalories_per_mole + print("this is the original mm energy: " + str(energy)) + oemol.SetData("MM", energy) + + #repeat with removed nitrogen + ffn = ForceField(FFRmNit) + topology = generateTopologyFromOEMol(oemol) + system_n = ffn.createSystem(topology, [oemol]) + positions = extractPositionsFromOEMol(oemol) + integrator = openmm.VerletIntegrator(2.0*unit.femtoseconds) + simulation = app.Simulation(topology, system_n, integrator) + simulation.context.setPositions(positions) + state = simulation.context.getState(getEnergy = True, getPositions=True) + energy_rm = state.getPotentialEnergy() / unit.kilocalories_per_mole + print("this is the original mm energy: " + str(energy)) + oemol.SetData("MMRmNit", energy_rm) + + return oemol + +def findAngles(oemol, constraint): + """ + Description: + Calculates the improper (and potentially valence) angles of an oemol and stores the value in + tag "improper" and "valence" for an oemol. + + Input: + oemol: An oemol object + constraint: An constraints.txt file used in an input for an geomeTRIC scan + + Return: + oemol: An oemol with the calculated improper or valence angles with corresponding tags + "improper" or "valence" + """ + + #determine the constraints + #open the constraint file + constraintfile = open(constraint, "r") + f = open(constraint) + lines = f.readlines() + f.close() + coords = oemol.GetCoords() + + for l in lines: + if "dihedral" in l: + split = l.split() + atoms_imp = (int(split[1])-1, int(split[2])-1, int(split[3])-1, int(split[4])-1) + print(atoms_imp) + crd1 = np.asarray(coords[atoms_imp[0]]) + crd2 = np.asarray(coords[atoms_imp[1]]) + crd3 = np.asarray(coords[atoms_imp[2]]) + crd4 = np.asarray(coords[atoms_imp[3]]) + angle_imp = calc_improper_angle(crd1, crd2, crd3, crd4, True) + print(angle_imp) + if angle_imp < 90 and angle_imp > 0: + oemol.SetData("improper", angle_imp) + if angle_imp < 0: + if angle_imp < -90: + #angle_imp = 180 + angle_imp + oemol.SetData("improper", angle_imp) + if angle_imp > 90: + angle_imp = angle_imp - 180 + oemol.SetData("improper", angle_imp) + print(angle_imp) + + #calculate the valence angle and store in tag (2-d scan) + if "angle" in l: + split = l.split() + atoms_val = (int(split[1])-1, int(split[2])-1, int(split[3])-1) + crd1 = np.asarray(coords[atoms_val[0]]) + crd2 = np.asarray(coords[atoms_val[1]]) + crd3 = np.asarray(coords[atoms_val[2]]) + angle_val= calc_valence_angle(crd1, crd2, crd3) + oemol.SetData("valence", angle_val) + + return oemol + + + +#TODO +def makeOEB(oemolList, tag): + """ + Description: + Takes in an oemol list and creates an output OEB file. + + Input: + oemolList: A list of oemols + tag: The title of the OEB file. + + Return: + + """ + ofile = oemolostream(tag+'.oeb') + for mol in oemolList: + OEWriteConstMolecule(ofile, mol) + ofile.close() + return + + + + + +def adjust_energy(oemolList, qm_tag, mm_tag): + """ + Description: + Iterates through oemols in a list and normalizes the energies to the lowest + QM energy in the list. The corresponding MM energy geometry is subtracted from the + MM energies tagged in the oemol. + + Input: + oemolList: A list of oemols + qm_tag: The tag for the qm energy + mm_tag: The tag for the mm energy being normalized + + Return: + oemolList with the updated energies. + """ + + #sort the oemols based on the lowest QM energy, this oemol is now "low_mol" + low_mol = sorted(oemolList, key=lambda x: x.GetData(qm_tag))[0] + #get corresponding lowest qm energy + low_qm = low_mol.GetData(qm_tag) + + #get correspoding lowest mm energy + #low_mol_mm = sorted(oemolList, key=lambda x: x.GetData(mm_tag))[0] + low_mm = low_mol.GetData(mm_tag) + print("this is the low mm value:" + str(low_mm)) + + #iterate through the list and subtract lowest mm and qm energies + for m in oemolList: + qm = m.GetData(qm_tag) - low_qm + m.SetData(qm_tag, qm) + mm = m.GetData(mm_tag) - low_mm + m.SetData(mm_tag, mm) + + return oemolList + + +def plotResults(oemolList): + """ + Description: This function plots the QM and MM data for 1d torsion scans + + Input: + oemolList: List of eomols with stored angles, MM, and QM energies + """ + #sort oemolList based on the improper angles + sort_oemol = sorted(oemolList, key=lambda x: x.GetData('improper')) + + xs = [m.GetData('improper') for m in sort_oemol] + qm = [m.GetData('QM') for m in sort_oemol] + mm = [m.GetData('MM') for m in sort_oemol] + mm_nit = [m.GetData('MMRmNit') for m in sort_oemol] + + plt.plot(xs, qm, color="royalblue", label='QM', marker='^',markersize="8", linewidth="3") + plt.plot(xs, mm, color="red", label='MM', marker='^',markersize="8", linewidth="3") + plt.plot(xs, mm_nit, color="orange", label='MM removed Imp.', marker='^',markersize="8", linewidth="3") + #plt.plot(xs, qm, color="royalblue", label='Carbon 1', marker='^',markersize="12", linewidth="5") + plt.legend() + plt.show() + + +def oeb2mollist(oeb): + """ + Description: + Takes in oeb file and creates oemolList + + Input: + oeb: oeb file + + Return: + oemolList: a list of eomols contained in the .oeb file + """ + + #open input xyz file + ifs = oechem.oemolistream() + ifs.open(oeb) + + #generate empty list of oemols + oemolList=list() + + #iterate through oemols and set coordinates to original oemol + for mol in ifs.GetOEGraphMols(): + oemolList.append(oechem.OEGraphMol(mol)) + + return oemolList + + + +def compareGromacs(oemolList, FF, top, gro): + """ + Description: + Compares energies of openMM to gromacs + + input: + oemolList: List of oemol objects to compare + FF: openmm force field, .offxml file + top: GROMACS topology file name to write + gro: GROMACS coordinate file name (.gro format) to write + + + """ + for mol in oemolList: + ff = ForceField(FF) + topology = generateTopologyFromOEMol(mol) + system = ff.createSystem(topology, [mol]) + positions = extractPositionsFromOEMol(mol) + save_system_to_gromacs( topology, system, positions, top, gro ) + + + + +#Main +def processData(smiles, xyzfile, constraintFile, moltitle, FF, FFRmNit): + """ + Description: + Takes in innitial .sdf file, xyzfile from geomeTRIC output, constraint file and a title to + create an .oeb file with oemols with QM, and MM data. + Also generates plot comparing QM and MM energies. + + + input: + #UPDATE to .sdf or .mol2 file + smiles: Smiles string for molecule + xyzfile: output .xyz file from geomeTRIC scan, scan-final.xyz that contains QM energies in title + constraintFile: constraint.txt for geomeTRIC input that contains the indicies constrained in scan + molTitle: title for the output .oeb file + FF: Forcefield + FFRmNit: Force field with removed nitrogens + + Return: + none, generates .oeb file and plot of data + """ + + #generate lsit of oemols with stored QM energies + oemolList = QM2Oemol(xyzfile, SDF2oemol(smiles)) + + #get MM energies and improper angles from geometries (and potentially valence depending no the scan) + for mol in oemolList: + GetMM(mol, FF, FFRmNit) + findAngles(mol, constraintFile) + + #normalize eneriges: + adjust_energy(oemolList, 'QM', 'MM') + adjust_energy(oemolList, 'QM', 'MMRmNit') + + #write out oeb file + makeOEB(oemolList, moltitle) + #plot results + plotResults(oemolList) + + return + + + + + + + +#plotResults(oeb2mollist('results.oeb')) + +compareGromacs(QM2Oemol('scan-final.xyz', SDF2oemol('mol6.sdf')), 'smirnoff99Frosst.offxml', 'topology', 'coordfile') + +#test +#processData('mol6.sdf', 'scan-final.xyz', 'constraints.txt', 'fixresults', 'smirnoff99Frosst.offxml', 'nitRem.offxml') + + + diff --git a/off_nitrogens/README.md b/off_nitrogens/README.md new file mode 100644 index 0000000..c7c2b98 --- /dev/null +++ b/off_nitrogens/README.md @@ -0,0 +1,10 @@ +# OFF: Nitrogens Project +README last updated: 2018-07-06 + +This repository contains the codes for generating different geometries of trivalently bonded nitrogens, calculating the improper and valence angle in preparation for QM caluclations to map the potential energy function of an improper torsion. The purpose of these codes is to parameterize the force fields improper torsions. + +## I. Repository contents + +| Script | Stage | Brief description | +| ---------------------|---------------|----------------------------------------------------------------------------| +| `perturb_angle.py` | | Create the various geometry moelcules. | diff --git a/off_nitrogens/Update/README.md b/off_nitrogens/Update/README.md new file mode 100644 index 0000000..5a2d77e --- /dev/null +++ b/off_nitrogens/Update/README.md @@ -0,0 +1,18 @@ +# OFF: Nitrogens Project +README last updated: 2018-10-18 + +This repository contains the codes for generating different geometries of trivalently bonded nitrogens, calculating the improper and valence angle in preparation for QM caluclations to map the potential energy function of an improper torsion. The purpose of these codes is to parameterize the force fields improper torsions. + +## I. The Pipeline +![alt text](http://url/to/img.png) + + +## I. Repository contents + +| Script | Stage | Brief description | +| ---------------------|---------------|----------------------------------------------------------------------------| +| `perturb_angle.py` | | Create the various geometry moelcules. | +| | | + + + diff --git a/off_nitrogens/Update/mmEval.py b/off_nitrogens/Update/mmEval.py new file mode 100755 index 0000000..9ab1c8d --- /dev/null +++ b/off_nitrogens/Update/mmEval.py @@ -0,0 +1,119 @@ +from openforcefield.typing.engines.smirnoff import * +from openforcefield.utils import get_data_filename, extractPositionsFromOEMol, generateTopologyFromOEMol +from openeye.oechem import * +#import oenotebook as oenb +from openeye.oeomega import * # conformer generation +from openeye.oequacpac import * #for partial charge assignment + + +bondList = [] +def energyminimization(oemol, nsteps, pert, outprefix='molecule'): + """ Energy minimization calculation on oemol + Arguments: + _________ + oemol: (OpenEye OEMol object) + OpenEye OEMol of the molecule to simulate. Must have all hydrogens and have 3D conformation. + nsteps: integer + Number of 2 femtosecond timesteps to take + outprefix (optional, default 'molecule'): string + Prefix for output files (trajectory/Topology/etc.). + Returns: + ________ + status: (bool) + True/False depending on whether task succeeded + system: (OpenMM System) + system as simulated + topology: (OpenMM Topology) + Topology for system + positions: (simtk.unit position array) + final position array + + """ + + # Prep forcefield, create system and Topology + ff = ForceField('smirnoff99Frosst.offxml') + + topology = generateTopologyFromOEMol(oemol) + system = ff.createSystem(topology, [oemol]) + + positions = extractPositionsFromOEMol(oemol) + + # Energy minimize + # Even though we're just going to minimize, we still have to set up an integrator, since a Simulation needs one + integrator = openmm.VerletIntegrator(2.0*unit.femtoseconds) + # Prep the Simulation using the parameterized system, the integrator, and the topology + simulation = app.Simulation(topology, system, integrator) + # Copy in the positions + simulation.context.setPositions(positions) + + # Get initial state and energy; print + state = simulation.context.getState(getEnergy = True, getPositions=True) + energy = state.getPotentialEnergy() / unit.kilocalories_per_mole + print(str(pert) + " %.7g" % energy) + + # Write out a PDB + from oeommtools.utils import openmmTop_to_oemol + outmol = openmmTop_to_oemol( topology, state.getPositions()) + ofile = oemolostream(outprefix+'_initial_newfield.pdb') + OEWriteMolecule(ofile, outmol) + ofile.close() + + + # Return + return True, system, topology, state.getPositions() + + + + + +def input_energy_minimization(smiles, name): + """Reads in SMILE strings of molecules, creates OEmols and runs energy minimization + """ + mol = OEMol() + OESmilesToMol(mol, smiles) + omega = OEOmega() + omega.SetMaxConfs(100) #Generate up to 100 conformers since we'll use for docking + omega.SetIncludeInput(False) + omega.SetStrictStereo(False) #Pick random stereoisomer if stereochemistry not provided + + #Initialize charge generation + chargeEngine = OEAM1BCCCharges() + + + status = omega(mol) + if not status: + print("error generating conformers.") + OEAssignCharges(mol, chargeEngine) + + energyminimization(OEMol(mol), 100000, outprefix=name) + +#list of my molecules +#molecules = {1:'CNC', 2:'CNC(=O)C', 3:'CNC(=O)OC', 4:'CNC(=O)NC', 5:'CNS(=O)(=O)C', 6:'CS(=O)(=O)Nc1ncncc1', 7:'CS(=O)(=O)Nc1ccccc1', 8:'CNc1ccc([O-])cc1', 9:'CNc1ccc(N)cc1', 10:'CNc1ccccc1', 11:'CNc1ncncc1', 12:'CNc1ccncc1'} + + +#for molecules, use input energy minimization and run energy minization +#for k in molecules: + #input_energy_minimization(molecules[k] , name = str(k)+ '_' + molecules[k] ) + + + +def xyz2mol(xyz): + ifs = oechem.oemolistream() + ifs.open(xyz) + for mol in ifs.GetOEGraphMols(): + molobj = oechem.OEGraphMol(mol) + return molobj + + +k = 1 +m = -10 +while k <= 21: + filename = 'pyrnit_1_constituent_0_improper_' + str(k) +'.xyz' + name = 'pyrnit_1_constituent_0_improper_' + str(k) + energyminimization(xyz2mol(filename), 5000, m, name) + m += 1 + k += 1 + +#running script + +#energyminimization(xyz2mol('pyrnit_1_constituent_9_improper_21.xyz'), 5000, 'molecule') diff --git a/off_nitrogens/Update/one_psi.slurm b/off_nitrogens/Update/one_psi.slurm new file mode 100644 index 0000000..0636a45 --- /dev/null +++ b/off_nitrogens/Update/one_psi.slurm @@ -0,0 +1,41 @@ +#!/bin/bash + +#SBATCH --job-name="perturb_valence_60" +#SBATCH --partition=mf_nes2.8 +#SBATCH --nodes=1 +#SBATCH --tasks-per-node=1 +#SBATCH --cpus-per-task=1 +#SBATCH --mem-per-cpu=10gb +# one hour +#SBATCH --time=2:00:00 +#SBATCH --distribution=block:cyclic +#-------------- + +# Informational output +echo "=================================== SLURM JOB ===================================" +echo +echo "The job will be started on the following node(s):" +echo $SLURM_JOB_NODELIST +echo +echo "Slurm User: $SLURM_JOB_USER" +echo "Run Directory: $(pwd)" +echo "Job ID: $SLURM_JOB_ID" +echo "Job Name: $SLURM_JOB_NAME" +echo "Partition: $SLURM_JOB_PARTITION" +echo "Number of nodes: $SLURM_JOB_NUM_NODES" +echo "Number of tasks: $SLURM_NTASKS" +echo "Submitted From: $SLURM_SUBMIT_HOST" +echo "Submit directory: $SLURM_SUBMIT_DIR" +echo "=================================== SLURM JOB ===================================" +echo + + +cd $SLURM_SUBMIT_DIR +echo 'Working Directory:' +pwd + +date + +/beegfs/DATA/mobley/limvt/local/miniconda3/envs/oepython3/bin/psi4 -i input.dat -o output.dat + +date diff --git a/off_nitrogens/Update/optimizer.ipynb b/off_nitrogens/Update/optimizer.ipynb new file mode 100644 index 0000000..4f3a072 --- /dev/null +++ b/off_nitrogens/Update/optimizer.ipynb @@ -0,0 +1,135 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import scipy.optimize as opt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[0.067717792501689, 0, 0.353868942785597, 0.575564030264505, 0.457620451913129, 0.0696958207232113, 0.118606344568143, 0.705373266009391, 0.400032195845465, 0.00745054048863659, 0.198262812667296, 0.125492746306463, 0.26167665898275, 0.678608755286266, 0.0304297412223832, 0.183790772163396, 0.0304386242402993, 0.289900123135208, 0.529374195299206, 0.00779152636820928, 0.848325689701502]\n" + ] + } + ], + "source": [ + "pert_Deg= [-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10]\n", + "QM_Data=[0.705373266009391, 0.575564030264505, 0.457620451913129, 0.353868942785597, 0.26167665898275, 0.183790772163396, 0.118606344568143, 0.067717792501689, 0.0304297412223832, 0.00779152636820928, 0, 0.00745054048863659, 0.0304386242402993, 0.0696958207232113, 0.125492746306463, 0.198262812667296, 0.289900123135208, 0.400032195845465, 0.529374195299206, 0.678608755286266, 0.848325689701502]\n", + "MM_Data=[2.288293, 2.223676, 2.179397, 2.152777, 2.140672, 2.156729, 2.185541, 2.236188, 2.310113, 2.397269, 2.509073, 2.63514, 2.796373, 2.971429, 3.178807, 3.401509, 3.649064, 3.922269, 4.215558, 4.551246, 4.891297]\n", + "\n", + "#improper removed from smirnoffxml.offxml file for the updated MM data\n", + "update=[2.288293, 2.223676, 2.179397, 2.152777, 2.140672, 2.156729, 2.185541, 2.236188, 2.310113, 2.397269, 2.509073, 2.63514, 2.796373, 2.971429, 3.178807, 3.401509, 3.649064, 3.922269, 4.215558, 4.551246, 4.891297]\n", + "\n", + "#subtracted QM and MM data\n", + "def Diff(li1, li2): \n", + " return (list(set(li1) - set(li2))) \n", + "subtract = Diff(QM_Data, update)\n", + "print subtract" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8XNWZ8PHfmSJpRr1LLrIt926DCLCEZogxkACbBoQk7JIsb5JlUzZ503ZDsmHzbiAhsNlsNusQUgiBEBK6AROaaTZY7rZc5SarWl2aPnPeP+6MZkYaSaMyGo30fD8ffSTduXfmzJ07zz33OeUqrTVCCCFShynZBRBCCDEyEriFECLFSOAWQogUI4FbCCFSjARuIYRIMRK4hRAixUjgFkKIFCOBWwghUowEbiGESDGWRDxpUVGRnjt3biKeWgghpqTq6uqzWuvieNZNSOCeO3cu27dvT8RTCyHElKSUOhnvunEFbqXUCaAb8AM+rXXV6IomhBBirEZS475ca302YSURQggRF2mcFEKIFBNv4NbAZqVUtVLq9kQWSAghxNDiTZVcpLWuV0qVAC8ppQ5qrbdErhAM6LcDVFRUjHMxhRBChMQVuLXW9cHfzUqpJ4D3AVv6rbMR2AhQVVUld2cQQkwbm+uq2XhoE83Odkps+dy++BrWzzo3Ya83bKpEKZWplMoO/Q2sB/YlrERCCJFCNtdVc8/exyju3s/nXVtocrRxz97H2FxXnbDXjCfHXQq8qZTaDbwLPKe1fiFhJRJCiBSy8dAmznUd5L8cf+RmTzUbe/+A2+9l46FNCXvNYVMlWutaYHXCSiCEECns3M43+L+uv/YF0yWBZnIDTpqdiXvNhIycFEKIqU4HAux64zN8y/VSeBmw2bKETpONUlt+wl5bArcQQoyQ3+9lz1+vZ23r833LDptK+Jr9BtpNmaSbrdy++JqEvb4EbiGEGAGPp5fDL17O2u73+pbtS1vAD3M+RofHRekE9CqRwC2EEHHqcTRS/8KlrHAd7lu2N/t8Fl/1Mr9Py5ywckjgFkKIOLS2H6b3r5exyNvQt2xn0bWsvvJJTKaJDaUyV4kQQgyjoXEb/hfPpyIyaM/6LGuufHrCgzZI4BZCiCEdP/EM9lfXURLoAMCPYvei77D2kl+iTMkJoZIqEUKIQRyseYA5O7+ADS8ALiwcX/1TVi//fFLLJYFbCCFi2FP97yw79F0sBADoUjbOnv8wSyv/Nsklk8AthBAD7HzzC6w99T99/7eYcvFc8gyVMy5OYqnCJHALIQTGZFG/PPgMH+54lps94QmiTlvKsF/5MjMLliWxdNEkcAshpr3NddXcu/sR/tnxHFd5D/Ytr0mby6wNb5KdNTOJpRtIepUIIaa93x74C3f1PB4VtLdY5vP93E9OuqANUuMWQkxzbR1H+G7bAywKNPctCwDfsX0Q7e5NXsGGIDVuIcS0Vd/4Dr4X3jcgaD9pXUVAmShJ4Ax/YyE1biHEtFRb+wSF224hVxsTZ/tQ/DjjSp5LWwGQ8Bn+xkICtxBi2jmw97+o3PvPZOADwImVv877JtsdVtQE3TdyLCRwCyGmld3bvsGKYz/CjHFP805lp/WCP/ChedfzoSSXLV4SuIUQ04IOBNj5+qc4p+EPfcsazfnoS5+nsuz8JJZs5CRwCyGmPJ/Pxf7NV3NOx2t9y05YZ5G7/jXyc+cnrVyjJYFbCDGlOV1tnHjhUlY79vUtO2hbRsWGLdhthUks2ehJ4BZCTFkdnbW0v3QZSz2n+5btzr2E5Ve9iMWSkcSSjY0EbiHElNTQtA312jXM87f1LdtZdiNrLvtD0ubRHi+pXXohhIjh+IlnsL2yjrJg0Paj2FX5VdauezTlgzZIjVsIMYVsrqvm7V338o2uP2EL9tF2Y+Hoyh+zZuWXkly68SOBWwgxJWw+/R773vs2/+p6GUuwj3YX6by++Ad8aAoFbZDALYSYAlzuDth2G//sCfcc6cXKFzJvxNXhSJmBNfGSwC2ESGmNze/hef161kfcgR0gAx8nzYUoZ3uSSpY4qZ+lF0JMWwf2/zeZL19CRUTQPqny8aN4yroKYNLO8DcWUuMWQqQcv9/L7i23ck7DI33LPJj5WcY6nrAuB6WAyT3D31hI4BZCpJSu7lPU//VqznEe6FvWYsql+/zfssI6i7cPbaI5BWb4GwsJ3EKIlHHixLNkbv0kSwKdfcsO2ZZQduXzVGbPpRKmZKDuL+4ct1LKrJTaqZR6NpEFEkKIWHa/+21mvH0DxRFBe2fZjSy4bg+52XOTV7AkGEmN+0tADZCToLIIIcQAbncnB1++gdURM/v1qHROrfgRa1f+U/IKlkRx1biVUrOAa4EHElscIYQIa2rZQePTy6OC9ilLOd3rXmfZNA3aEH+N+37g60D2YCsopW4HbgeoqKgYe8mEENPS5rpqNh7axNyuHdzpfJ5S3H2P7cm5iEUfeJaM9LwkljD5hq1xK6U+CDRrrauHWk9rvVFrXaW1riouLh63Agohpo/NddX8aM8fubpjE/c4nyQnGLS9mNg1759Zec2WaR+0Ib4a90XAdUqpa4AMIEcp9Xut9ScTWzQhxHTzyIFH+ffuP3G+/2TfsiaVxY9zP86PLrw3iSWbXIYN3FrrbwHfAlBKXQZ8TYK2EGK8Hax5gHtbf06BdvQt08BnM2+hM2BPXsEmIenHLYRIKq/Xwb7XbmFty5NRywPAU9ZVdJjslE7BYetjMaLArbV+DXgtISURQkw7DY3bcL/xUdZ66/qWtSo7d9k2UG2ZA0zdYetjIZNMCSEmnA4E2LP9u+S9cjFzI4J2jX0FW6v+SF32GhRQasvn6ys/Pi1GQ46EpEqEEBOq19FM7Ss3sKrrnb5lHszsn/MF1lx4P0tNJq5d+MEklnDyk8AthJgwtcefJPPd21jpD8+RXW8uxnvh71hbsSGJJUstEriFEAkXCPjY/ebnWFn3aywE+pbvzrucResex5ZRkMTSpR4J3EKIhGprP8jZV29gretQ37JelU7t0u+xes03k1iy1CWBWwgxrkJD1pud7VxJPV/qfppF2tn3+LH0uWRf8gQri9cksZSpTQK3EGLcbK6r5p69j6F9Tr7ofoOPenb1PRYAdpfdxMpLfoPFkp68Qk4BEriFEONm46FNzPHU8ZPev5CLq295s8qi89xfsHbRLUks3dQhgVsIMS7cnm6ub3+Gmz3bsaD7lr9hmc/dGR/gWQna40YCtxBizGprnyB9++f4lK+5b5kGdptn8m3bhyi1S6+R8SSBWwgxai53BzVbbmV1y9NRw7B3mmdxt+0DnDHlyZD1BJDALYQYlSNHHiZ7xxdZ62/rW+YgjS1lt/CrQAVNri5Kp/Cd1pNJArcQYkQczhYOv/5p1rS9ELX8oG0ZBRf/gQ1Fq5ExkIklgVsIEbdDBx8kf9dXWRPo6FvWo9I5Nv9rrKr6Psok89ZNBAncQohh9TgaOfbaLazueCVq+QH7KkoveYTVBcuSVLLpSQK3EGKAyNGPV9DAP3Y/x2rd3fd4l7JxYtG3Wbn221LLTgIJ3EKIKKHRj2m+bv7F9SpXeQ9GPb4vq4qZlz7Kqtz5SSqhkMAthIjyy4PPcLlzF191vUIGvr7lbcpO/dJ/Y8WaryWxdAIkcAshIpw49Tx3tm5kpb8havlm6xL+K/0ynpGgPSlI4BZC0Nl9ghNvfpaV7S9HDaTRwFbzXO6yXS037J1EJHALMY35/V72vvtN5p/4Oat1eFIoLyYeTTuXh9Lfh1OlyejHSUYCtxDT1NFjj2Hd8RXWeOujltfYV3B04Td4qvE4Lme7jH6chCRwCzHNtLYf5sxbt7Gq662o5Y3mfNqW38XSZZ9nqcnEh5YnqYBiWBK4hZgmvF4H+7Z+hcWnf8MqPH3LXViomflJll34U8rSspNYQhEvCdxCTFGRg2guo4nbe19krb81ap192edR+je/Ym3hyiSVUoyGBG4hpqDQIJp871nucm3hUt/RqMfrLCU4Vv+IFYs/naQSirGQwC3EFPS7A3/m1t5XudnzXtSXvIc0jsz5B1ae/2MsloyklU+MjQRuIaYQn8/Fvne/yc9aN5IXcWd1gE3WZWxMfz9PXvSzJJVOjBcJ3EJMAToQ4MC++ymo+X+s6ZfH1sDrlgX8h+0qGUQzRUjgFiLFHav9C+z8Gsvdx6OWN6gcNmZcxMuWxWilZBDNFCKBW4gU1dC0jdZtd7CiZ3vU8h6VzpHZf0fTrE+w99jrIINoppxhA7dSKgPYAqQH139ca/3dRBdMCBFbZ/cJat/+R1a2vkA5gb7lXkzsK7qW+Rf+jLXZFQBsmHtJsoopEiieGrcbWKe17lFKWYE3lVLPa623JrhsQogIbncnB7Z9lUVnfs9a7Y56bG/2+ZSc/9+sLZEa9XQwbODWWmugJ/ivNfijE1koIYRhc101vzz4DKu7tvIP7rdZG3EXGoCj6ZWYzrmXlfNuSFIJRTLEleNWSpmBamAB8N9a620JLZUQghdPb+ONHT/kB863WBRoiXqs3lxE+7J/Zdnyf5Jbh01DcQVurbUfWKOUygOeUEqt0Frvi1xHKXU7cDtARUXFuBdUiOnC63VwYPcPWXT051wViO7a165s/MF+Cf/n2ieZIQNopq0R9SrRWncopV4DNgD7+j22EdgIUFVVJakUIUbI4Wzh0I7vMvv0w6wOdA14PADclPX3OFU6/yhBe1qLp1dJMeANBm0bcCVwd8JLJsQ00dFZy/Hqb7Og6SnWRtzMAMCBlTpTHvMDZ3nKugqHSpdBNCKuGnc58NtgntsEPKa1fjaxxRJi6mtqrqZhx7dZ2vYKayNuygvQoTLZXrCBn3ln0oq5b7kMohEQX6+SPcDaCSiLENPCyVMv0LX7eyzrfpfSfh20Gs0FNM65jaVr/5Ur03MJREzNWiKDaESQjJwUIoH65sR2tHGJOsuN7m2sdB8ZsN4J6yx6Fv4TS1d+hTKztW/5+lnnSqAWA0jgFiJBNtdVc+/uR7jAU8Nd7mqWBpoGrHPItgSWfoNFiz4t3fpE3CRwC5EAHZ21NLz3z/zeVU2x7o16zI+iJvs8slbdyeI51yaphCKVSeAWYhydOPkcnfvuYWnnW9yKP+oxDTxlXcUf08/lkQ89mJwCiilBArcQY+TzuajZez/22o3M7ze1KoATC+n4eN6yjHttV0h3PjFmEriFGKWOzlqO7/o+FQ1/YWWge8DjR6yzeMyyir9aFuBTRpc+6c4nxoMEbiGGsblfl7xbCwqYU/9HlnS+zdp+6RAvJmpyLiBr+f9l4bwbOK+ump3SnU+MMwncQgwhdLd0n8/FOt9RPtK7k5VNDQPWazNlcbLsBuasuZNVeQv7lkt3PpEIEriFGITP52bLrvu5w7mbq70HSO9XuwaoTZtDb+U/sGTlV1hrtSehlGI6ksAtRAS/38uxo4/gOP575ra/xb9rx4B1vJh4xbqIhVU/pHLe9UkopZjuJHCLac/v93Ls2KM4ao1gvahfv+sQDdSYSvmW/XqsmbN4XIK2SBIJ3GJa8vu91B57jN7ah5jT/uagwbpVZfGqdSGvWBayzzyj727pX5eeISKJJHCLaWFzXTUPHHyG0p4DXOk/xsWeQywcJFi3mbI4mX8x2fNvpbLyI+TU76b50Ca5W7qYNCRwiynN4WzllR33oOo38T++Wgpj5KwB2lUmJ/MvJnP+p6ms/ChrZaInMYlJ4BZTTlPLDuqP/A5b00tUOg/xwRi9QQDalJ1T+ReTOf9TVFZ+nDURwVqIyUwCt0h5Pp+bEyeepPvEYxS3vcksXzOlQ6wfAL5s/yh7zDN5bcP9E1VMIcaNBG6RMiJHMM5NS+NGm5vZne8yt3snC7Rz0O2OmYp5yzKPikAbF/uO8ZR1FTsts2XOEJGyJHCLlPDCybd4Zvf/sM57igt9x1nZdQYLse9J7cZCrX0p7rL1zFx4K8ecHn639zHcfm/fOjJniEhlErjFpNTdc4a6k8/ganyFrI6dXO4+zoZBctUAraZsTueeT1rFDVQuuIWl6Xl9j60P/pZbgImpQgK3SDodCNB8didNp59DN2+hsHsfs3xNLB1qG+CAuYy3LZWsP++7VMy6ksIh7iAjPUPEVDJpAnf/GdikRjQ1ba6r5lc1T5PjOMJ5uoVzdTPz3LWUBnqGbFAEOKNy8aGYrTt53rKUH9qvotSWzz9UrB9mSyGmlkkRuDfXVfPrXb/EHPCAyqXJ2c49ex8DkOCd4nw+F/WNb9Le8DpdTW9S3nuY3/qbycA39HaYOJ02k46cNVhLL+Vkxnzuq31H8tRCMEkC98ZDm7jRtY2PeXahgXZsbLXO4+R71Rzsvoni0r+hIG+J3Ex1kvN4eqlveI3Oxi3oth3k9h5ipqeeCvxUDLNtj0qnLmMBzoIqsso/wKyKq5mXUdD3+ArAmjlbrsqEYJIE7mZnO5X+VgAUUICTa7wHwHsAdr4MQLfKoCltBj2Z8yF3BZlF51JadjE52RWSZplAoX3d6WhmjdnBFRl+yl0nyes9zExvI3MJxP1cAeC+jHXsMc/kwQ/+iiXDDICRPLUQhkkRuEts+XQ4bLiwDHoJna1dZLtrwV0LbS9B8NZ+bSqTAlMhHzcXcsqUT6M3h0d3Hsfv+QRXV145ge9i6ul1NHP27A662/fj7TqEo72GAucp/jPQSZnuwjxId7z+mk15tNgr2e7PYh8FXOY9zAZfDU9ZV/Fk2mpKbfmYZdSiEHGbFIH79sXX8B97e3D7veQFHMwLtDI/0EaVxUOpt54yzxmytDvmtgW6lwJ/L1X+U9EPbP013dsyaDfn05NWjCdjBjpzNpasSuw5C8jNX0p+7iJebtgz7WrroVpzi6ONhWmK6/NLqKCHQPcRLI5TZLnPUOA9S652kjnC5240F3DWPh9v7mrsJRdSNnMdJdlzKQHO1FXz0N7H2Gqdxw+5CpA8tRCjobSOr9Y0ElVVVXr79u0j2maodIcOBGjtqOFs41s423Zg6qwh11FLmbdx2EauofhQtKhsmkzZtKgsupSNblMG5bkLKM9fiDWjmDRbKTZ7GXZ7OVn2GVgs6XGVeazveTTbXlqyiJ6e0/Q6zuDqbcDjbMTnaibgakG5WzF728F9Fqu3k1ztpEj3kjZE3+jhaOAlyxIOm0tYv/KzlM+8kuzM8oS9ZyGmMqVUtda6Kq51J0vgHg2/38sdL32VfOcJbnG/y9JAEy1k4jFZKQl0x7zV1Fj1qHR6TXZ6lI02baVLpdOtMnArCx5lpTxrBoX2YjDbUOYMTBY7JrMdk9WG2ZyJ2WrHYslkX2cTj558B4cGK35s2ku2CrC+dBmV9mwC3m783m60z4H29aB8veDvxeRz4vZ04PF0kqE9ZGs3OdpFrnZixzv8GxgBD2bOmgvoTC/FZZvNu04vx7WdK7wHWec7wlPWVdxnW0epLZ/Hr/jOuL62ENPNtAncEL6Za/9uYl9f8VHel1NIR/sBersO4+muhd7TWF31ZLmbyPO1kTfIfMzTSWggS73KpahoDabshdjzlpBfuIaCvMWYTOFs2qD7euXHpdYsxBiNJHBPihz3WIQCxmCX33m5lYNue/NL3wFHHV9wb+Ei33H2mMrZYamgQHmZn56G1dtFur+bDH8vmQEHWdrJZO6Q6EPRrez0mjNxmnPwWHPxWfMIpBVAehHmjGKebDzEGZ/iw56dXOk7HF1rXjd0rXm4fS2EmBgpX+Mei5HWIAMBH72OJnod9fx4+/+iPW3c6K6myn+KfaZytlrnkaVgaXYJyu9CBdyogBtzwI0p4MYc8GLWHswBL/gdpGkf+dpBBj66SOeMKR+3KY30tDz8ZhuB4A9mO9qSibJmoixZ/LXpEG1+zTXefVzoO8FfLYu5z7aOLFspj19557i+ZyHExBjXGrdSajbwO6AMo+vtRq31f46tiJPDSGuQJpOF7KyZZGfNZP0qE/fsfYx3LXP7Hg8FwLVxBMCxBNDm4LZvWBdEbXv7kmuHfV2pNQuR+oatcSulyoFyrfUOpVQ2UA3coLU+MNg2qVLjHqvJ1qtEgq8QqSuhjZNKqaeAn2mtXxpsnekSuIUQYryMJHCPqK1NKTUXWAtsG3mxhBBCjIe4A7dSKgv4M/BlrXVXjMdvV0ptV0ptb2lpGc8yCiGEiBBX4FZKWTGC9sNa67/EWkdrvVFrXaW1riouLh7PMgohhIgwbOBWSingV0CN1voniS+SEEKIocRT474I+BSwTim1K/gjswIJIUSSDNuPW2v9JsY02UIIISaByTyCWwghRAwSuIUQIsVI4BZCiBQjgVsIIVKMBG4hhEgxEriFECLFSOAWQogUI4FbCCFSjARuIYRIMRK4hRAixUjgFkKIFCOBWwghUowEbiGESDESuIUQIsVI4BZCiBQjgVsIIVKMBG4hhEgxEriFECLFSOAWQogUI4FbCCFSjARuIYRIMRK4hRAixUjgFkKIFCOBWwghUowEbiGESDESuIUQIsVI4BZCiBQjgVsIIVKMBG4hhEgxEriFECLFSOAWQogUI4FbCCFSzLCBWyn1oFKqWSm1byIKJIQQYmjx1Lh/A2xIcDmEEELEadjArbXeArRNQFmEEELEQXLcQgiRYsYtcCulbldKbVdKbW9paRmvpxVCiJTQ2unnyz9poq3Tn/DXGrfArbXeqLWu0lpXFRcXj9fTCiFESnhoUyd7j7l56PnOhL+WpEqEEGKMWjv9vLC1B63h+Xd6E17rjqc74CPAO8BipVSdUuozCS2REEKkmIc2dTIz04dCEwjohNe6LcOtoLW+OaElEEKIFHa2w8cr27p4+Np2Otwm/nzIxuat8KmrcynINSfkNSdVqmQik/tCCDEWXp/mxa09fOYHDXxgjovsNM3sbD+fXdWD35/YWvekCtwTmdwXQoj+4qk8djsC/OHFTj7xnXru/l0bvY4AH1nk6HvcbgWPXyU01z1pAvdEJ/eFEKK/oSqPDWd9/OyxNm78lzM88FQnrcEYddFMNzOyAgBoDc8dywBIaK570gTuhzZ1EjDe+4Qk94UQItJglcea427+7YGzfOq79fzltR5cbh213ccXh2vbD9fYua86BwCfP3GV0GEbJydCa6ef59/pwRd8f6E3nMjkvhijQABMk+a8L8SY9a883v3QWVxu2HvMPWDdeTOs5GaZ8Lf3sLzIB4DHD08esUWtF6qEfummgnEt66T45j20qRN/v5OS36/53aaO5BRIDC4QgLd3wRs74L399B3pQqSwUG07svL43gH3gKB97pIM7r6jmLvvKObAcTcfXhiubb9yKoM2V3RFM1G17qQH7tAOC0RffeAPwNNv9PKTh1upOeFGax37CSKeR3qkTICT9eA1ahg4nLD3SPh/Mf7qGmHLdnh9Oxw+mezSTFk/f7wd3yCHsdkE68/P5JffLuNHXyzhvGU2fv98FyUZft4/MxzY/3TIHnP7RKR+k54qibw8ieXZt3p59q1eygvNXF6VyboqO/NmWFFKDXieUKPCeF+WiKAeB5xuil7W0Q07DsDyBZAV+8AVo1TfDMfqwv83tMCiOckrzxTj9gR4fYeDJ1/v4eBJT8x1zCb4+TdKWTg7PWr5/uNurl/gwBys+m5vtHK8M3Y49flhX+3AdMtYqOFqsqNRVVWlt2/fPux6rZ1+brnzDB7vyJ5/TrmVdefaubzKzqwSa9TzpFkVf/j+DMmNjzetYWcNdAcvDdPTwB1xsJtMsGQuFMtJc1y0tMOBY9HLTAouXAMWObaH0trp565fneXOzxTFjAO1Zzw8+2YPf323lx7n0PHPYoZrL8oaWBn0+eCdPeFU4cqFUJA7pnIrpaq11lXxrJvUGvdwtW0ApYyzni8iA3Kywcuvn+3k1892sqgiDYuZAT1SpNY9zuqawkFbKVi1yEiVHDxu5LUCAThQCxVOmDvDWEeMTkc31NQOXB7QcKoBKmdNfJlSSKyrb6c7wGvVDp59s4eaE7Fr17EM2lGiviUcdOwZkJ8znm9hWEnLcfdvDBiM1mAyKb7+qQIur7KTkRYdEA6f8nDguCeqUWHT2z0jynWPOj/u98PxM1B9wMg/tnVOzcY6pxtO1If/nzvDOFiL8mHtUrBFXEaeaoD9xxj2gxWx9Thg31HjwAdj386PCNR1TeB0JadsKaB/l773Dji575E2PvatM/zo920DgrY9Q2Eapo4xIEcdCMCZ5vD/s8omvKKStMAdT207JBDQHDrp4Tu3FfHnu2fyr7cVctEqG9ZBrhe8Pvj7u+p58JkODp/yDNuwOeIRm1pDUyu8u88IVD0OI/+49wi8vduoLTW3TY3gpTUcPhE+IWXaYFZp+PFMmxG8I2scrR1GWsUhAWZEXG7jGAp1sUqzwspFMLMUsjONZVpH571FlMi44vFqvvGzFp55oweHKxwDLGa4/Fw7d36mEJ9fD+gY0d+AniEt7fTld60WKJ34q/ukpUr2H3fHHdcik/u2dBPrqjJZV5XJqUYvn/1BQ8zn6XZofv98l9H6m2/motU2LlptZ9WCdCzm8Nmx/xl62L7jPQ44cgq6emI/7vcbQbu5zTgL5+dAUR4U5hlfxFTT1Gpcuocsnjuw/7bVYuT4auuMGiEYQXtnDSytHHPub1rwemHPkXBAMJuNfRq6mlkwG3YeNP5u7TCu7mS/9vF4Na+910NOaxMPXuXiQKuVe9/LxhsIf9dnl1q49qIs1p+fSV62mfsfaRtR5fGh5zv50o354WMcYGZJUsYzJC1w//Lb5WN+jr+82j38SkBzu58nXuvhidd6yLabuGBFBu9fY6dqaUbMEZsx8+NeH5w4Y+S2IpmUkXu02yDgB1fEpZjWxhesrRM4CTmZRnqhKA9sGaN70xPJ44Vjp8P/z4qo+fWnFMyfbfQsOXzC2Cc+v1GDrJxlbCt579j8fth7NJwCUQpWzI/upZOTBaWFxokUjM8lL3vKDoIaroERjLEeOw65eLXaQfPJDr6wsov1K41a3KxsP5W5Xm7fXMjMEjNfu6WQVQvSo3qjjary2NFtVN7A2PczknPTmKR3BxyteHLkJgUZ6dFX7N2OAC+96+Cldx1YLcYHEsqkxGyI0NpIgxw/E536UMqOlcyKAAAdZ0lEQVQ4286ZEW7l1xp6nXC2A1rboccZXaCuXuOnti6cIy7KM76gkzGoHT0dfs8ZaUZuezilhcZ7238U3MHaY22dcbAvmktf/ylhCASM3iPdveFlS+dBXozGrspZcLbdaAx2uIxKRGTaagoZrHtvIKDZe8zNq9sdbNnpwO3y8dlVvXzt/c4BueoF+X7OK3OzuzWD2SUDuxCPqvK490j479JCsCbnKjplv0Xx5MhNJriiKpMff7GEGy7Nojgv+szt9YWDdojPp7nv0VZ8fg2dwT7KR05FB+38HKhaBvNn09pLuGFTKSMIz50B5y6H81catdDc7IGFc7iM/PiOGthSDYeOj3JPJEhrB7S0hf9fOMe4fI9Hdiacs8yoJYY0t8Gug0YeVxi0DjZqd4WXLagYvEtlmtWoKIScqGfEfWlTQP/0ZWuHj5oTbn7+eDs3/Ws9X7mvmaff6GFZjpNfX93G3y4MB22nV9HhCgfob5zfRY7VPz4DYBzO4NVzUBJPmkntxz1aI+n/HdmvW2vNkdNe3tjl4PUdDuqaYw+VKrL5+fzaHi6f3S/IZKQZgbgwr6+GfP8jbTzzZg/XXRyjr2ckrxdaO43aeHsnMVtELlo7Ofro+vywfV+4xlxaCEvmjfx5AgE4egoazoaXWS2wbL5xmT/d1dbB6cbw/xXlMG/m0NsEArB9v9HTB6C8yLiSmULuf6SNTW8bV9MmZXztHBFfxYIMP/90Tg+X9vt+bqtP4/7qbFx+xQNXtVFoM2p22xrS+N47eTz8/ZljG99x+ET4WC7MgxULRv9cMYykH3dK1rhH2iMldLZVSrGoIo3PXJfHOYszBsRIq0lz05Jefnt1W1TQdvng8dpsfnZ8Bm/Vp+MIzg42oqlorVYoKzI+7L9ZA8vnQ0m/QB/ZoyCZjteFg7bVYpysRsNkMoLKwopwKsjrgz2HjfTTdFbXFB20y4riS0WZTNGfR8PZcM51EhlNF1u/X/PWbgfPvRVOgQZ0OGgrNNdWOvntNW1RQdvhN/H/tuXwrTdyaXKY6XSb+OG2cKrp/HIP1893jK3W7fGG2xcg6SmqlMxxj7ZHSkis/Pj55W6+sLaH2dnRT/zaqXR+sTuLZocZcPCXLcYw1+WV6Xi8Gv9oBv6YzcH8dj49aXay6oLdu7p6jD68KxYmLxfc2R3dALuggkH7XcZrRonReHvgWDg/dfgkNLXBmsVje+5U1NwW3ehbmGcMZY+3naMg10jXtQdTLEdPwerFk6qdJN4pKDp7/Lx3wMXWfU6217jo6o1dI5ub6+PfLnMwOyO6i6mroIBbf2eitTf6+1LdlMajNXZuWmqc1D6zoocvv5ZGW+coZxytbwlfJWfZITdr6PUTLCUD91h7pETW2DMsAb75vm4u6XfZdbzTzLON+dR77XT53EA4teEPwJ6j0ev7/PDsWz1cvMbGmkUZmIbr1R/0wLtppDVn8YU1we6FHd1Gw96KBRPfYyAQiJ7IqCAXivPH57nzsuGcpcZ7CzXadnZDcyuUFI7Pa6SC9i5jtGlITqbRGDmSoKuU0T1w+wHjJNjZY/Qt7n8FlyRDdbHVWlN7xsvWfU627nNSc9wzZD9qi0lz0xIHn1zWS1pkvM1Ih0Vz+MULPjpdsbvmPrgvk7WlHhYX+LCa4ZvndfLIC+38441FI3tDgYAxb0zI7OT3kErJwD0WkbXtLGuA/7iko28+XYAej+LX+zJ5+qgNs8XEH75fSJbdxL5aN9trXGw/4ORoXezkut8PX/tpC1k2xcoFGaxemM7qheksmJWG2Tzwgw6VxeO1k2GB21YED8D2LqN2umz+xAbvUw3hLjhmk9EgOZ4HaEa6UTPcto++qdhqjhvvsWicThCTWXevceIKtSvZM4JXV6OoAdptxpXMmWCf4to6KMwd3XONs/5dbH/9bAcXrLCxdZ+TbftdnO2I73J5WaGXr57XxbzcfuvPLjMaac0m9h+PPY4DwBdQ/Ps7OWxc347Nqpmd42eVsxkYYeBuag3PgJlunRTHako2To5FqOEj2xLgnss6mJ8XDtoBDR99qogOtxEsB5tg5lidh8/f3Rh3usaeoVhRmc6qhemsWpDB4jlpWC0qqhHGYobvX+Xlguz28IZF+bCsctDgGU9f17j1Oo2h+6HjYUGF0d0xEbxe2HUofJJQyrjCmMoDSpwuYwBNZABYs9RoeRstn88YvRt6zjnlMHeYxs0RGM3xNdKJ45SCpXPTWLUgnT+/2o3XB3ZLgM+u6uW6BdFd/A61WSi7cD65ZSNs2G48C4dOhP9fVhn/ZGhaG43BoWO1cpZx4kiAlJlkKhn2H3dTmO7nnss6ovLZAQ1PH7X1BW0YfDrGZ94YZNTkIBwuzbsHXLx7wAV0km5VLKqwcuC4py9H7vPD916y8sTnSrE1B2tRZ9uNy+olsS+lx20qW62NAzsUtHMyEzuwwGo1JqnafcjoHaG1URNduTB2/+VU5wmOigwFWIvZGMo+lqANYLEYvVBC6a3TjUYjZ0b60NvFKd7jKxDQHK/3svuImz+/0jVs0M6yKc5bbuOC5TbetzyD3CxjFKPWsKbEwzfe10VpZjjX7fQqHtyXyTPHbFzt9vKlm0b4RkoLjS6Xoe6th08aXVbj2U9tXdFXoeUjrK0nyLSrceNwGb0aIqckXTzXOODjEE+NwmqB2z6Ux7E6D7uPuGmJ89JQKVi1II07L3eR3xnRha600ChjRPAe16ls65rCjWVKwbnLjDlIEs3lMfp2hz4Lk8kI6Elu+BkXWhu9PZpajYat0PfMpGDV4vF7j1obYwFCPUuK840U2xgNdXz5/ZrDpz3sOeJmz1E3e4+6hp0eVQHXX5rF5efaWTYvPSp12Nrp5++/V8enlvbwscXRg9a0hk88W0iTw3jtUR/rPp9xRRka2ZyTZTSMD5cK3H0oPOXDzBLjSjRBpMY9mB6HEbRDNR+ljMumEeSs4umKqDU0tvr49t8XobWmsdXP7iMudh9xs+eIi4bW2IFca9h9xMNHjii+UmXjQ/ODB3FTK90uyFw5B1Owt0ncQ/UHEboM/t6ncsg7fib8wJzyiQnaYNQ4Vy8y0iYer9EItPeIsWywofWTncttBOvmttiTbC2dP74nplBD5a5Dxv8t7UagCfaTH206rf/xdd+jrSyuSGfPUTf7at0Dbpg7HLPZOL5XLhg41cOLLzXz08vbmBuRy3b7wGqGp49m9AXtUFlGdYVpscCSSqOiAEYPrpMNQ3fB7HFEz9MziUapTp/A3dkT3U/aZDL6Uo8grxrvVLT9h86XF1koL8piw4XGF/Y/fnOWV7Y7+tIkAynu356FSWmurTS+/NmdrTz9WwevdRZSUWZl09u9Y7q5snEZ7KJjx1nyMiLmFU5Q/m5Qtoxw8Pb6jM9nz2FYs2R8TyChfrg+P2TZjC5dGenj0/jq9RlpraZW4zgbTGjCsfGWm230KGkOpgKOnjKumpQacTpNa82BE56+thcwdtlbu128tXvw2R5zsxTdvYPPtDfodBKnGvh4YT2WiDb4bfVp/Pi9bFpdA4/lMd1NJjfLCNShKYpP1kN+duyRzRA9mVRx/riloMbD9AjcbZ3GHNGhKoTFbLTmj7DmM5qBP/2/LK2dfl7fOVTQNmgU923PxmqC9XONL8x18504atrY+EYmxsVnmNen+f6vWrjtujzmzUgj2z54b5TQCWhdhZu5GRGXprFm/htk+3FrFAWjh0Qo5+3zGz+7DxnB2z7Gybi8PiP3e6Z54FzpZnM4iGfZjVq+PSO+YB4IGMdVU6sxIjZWytEc7C1TWmjUgEd5kohrf1fOMkblBgJGQ3NDC62ZhUPOfOkPaE43+ThyysPROk/f7wx8fGKJi48s6iXLCkc7zNz1di51PeFwUZJv7mtsX70wncdf6eL5t3sJDFGpifpOOFxG+013bzhoBwcXnX9JEX+6OUHd7SrKjV5boRNszXFj+gpLv1Do9oRPhDCpatswHQJ3S7sxP3boi2W1GEFiFPdHHOvAH4gv+FvMcOV5mVx6rp3DJ3PZ01nPqlwjh3nTUgeeAPxmX/RJR2vYc9TDl39i9DctyjMzb4Y14ieNOWUW0tNMPLSpk2xLgDvWRlwGziyJnltkCAm5v2eW3Wiw23PI6Cjv9Rl/r1kyupqO12fUmM40MehZ0u83vsCRtWSTgkx7OJhn2Y3gbjIZO7mr1wjWLUPMt16QYwTrwrxx6Z4X1/5OT4OKsnBt8ng9j51WUemO/3qsjapltr4AfazOi9trfC8UmnPLPHz9HCcXzvBEjf9amO/nd9e2saPJiqe4iLkrSygrCk+u1Nrp58WtvXFdib7wTg+frfKQ2VAffSLNyTQa4RM9a6ZSRsqker9RILfHaKxc2q/31pnmiMb6rLi/GxNlagfu/t2A0tOMoD3KWtxYB/6MJNXySrWDz16fx/nLbRDIQR+oRbV2APDp5Q68fsXDNYPngc92+DnbYYxKC1EK5hUrTB4Pd6x1kJtuHJhNDhPWgjLiCcEjnr98kOeIWYPMyTR6luw5Ynyp3d5wzTs9zh4YvmDArmseOH2AxWzs3DSr8aWMdXf6gDb6W0fO1gdG2sYfGHySrGy7MZCopGDAvOtjuUIZ0f6eVYZuOItye8Dno6ynGZ/fSAP4/PD6Tiev74xu/MtND7BhnpMPzXcyI2voGsU5pV6gAY6che4io4dFRnrcV6JFNj/fOL+LzDMRLftKGemL2RN4F5mMNGMqhtA9PVvaIb813GPE74+ekmH25Kptw1QO3GeajGlJQ2zpRtBOYp5q1KkWkwm1rBLP7qOkdRnDnD+zqhdvQPHYoegrB6tJU57tpyTDT3lWgLJMP+WZfsqyjN+hYB3pvvey2bGpkcqZaZQXWZhRZDF+Fxt/F+ebMQc71I61UTT0HIPWIHOzjT7de48YwdXlgd2HjR4AadbBg6DPb3zmdU0Da8L2DJg7g1ZrDnc92Gpsm2My8t7dDqMRKvTjHuR+hL3OgcvS04yadUnBkPn4sVyhxNrfd3wsn8ZWH6ebfdQ1e6lrCv5u9rE4M4N/u8h4D9fNd/LMURsnuvp/zTUri7xct8DJJbPdWGNkx3Y2WXn6mA2XT/HB+U4uKI+ohXu8xmCtUw1QkIupw0IgYKJ/+i7y9a6Y4+ZL53STlRZx/GXajFr2KK5+x6w43wjUoUmjjp4yUqf2DGhsjZjOON24cppkJk/g7ugyZpPJshkf6GgvMYMNHlH3SMwM5lCTfAeaMaVaTCb+92ABF5mdwZoPfG5ND7OzfUb30kw/5Vl+imyBYe+hFymg4d1G42R2+JSHw6cGBi6LGcoKLRTmmtl7zN0XSHx+eP7tHm6+KoeS/PgOpbhqkPk5RsPx/mPG5+kMduFcvZiHNnVFB0Gf3xiOfLpxYMC2GQGb4nyjoe6Rtuht09OMn8gGQ2+MYO4Mfw5ag5pRZNSuc7OGrSWO5golENB09gY4VjewkfDpLT0880bPoI2Aze3p7GyysrbUi9kEX1jbzddfzwMUmdYAV811cctqN/nmgf1ZXQETm46l89QRG6e7w5/ntoZ0im1+rp3v5NpKV9+sewC0dfLFZfDFtWlQXmwEw8jvmdcHR05CS8TUtWDUsOfOSO6NIObPNtJkDpdxhVdTa1zdRTZKTtIbgMTVj1sptQH4T8AMPKC1/uFQ64+qH/ehE0ZqIyQjPRzEM4N5xuF6AWgdffssMC6/Vywc+0RJSRbqV6sCmh9e0sHqktHNw+zxQ2OvmTSTpiQzwHO1Gdy3feyDXvKyTBTkmikM/eSYo//PNVOQY+bnj7dHjRaNNTK1T0ubcef4IJ/Nzsf/aKfDaSI7Ax75nAV7S3N4+HyILd0YEl1S0He8jKXfe1ubhwMvHOTCGR6eP27jbz68JO5t+4+OveqCTD5+ZQ6tnf6+dNbZDh9nQ/93+mnr9I/6dqUmEywp9PGfl7X11ZD/d1cms3L8rKtwYYv1NcjOpCevkFt+5qTbNXSQMivNJRUevvmBANbuGHegUso4EZYXh+e+iRz0kJEGi+dNnml9exxGP/hQHMzODKfJLGa4YNWETSMwkn7cwwZupZQZOAx8AKgD3gNu1lofGGybUQXuHQeMms5QTCYjkGcFg3lmMLBbLcaOP3Iyeu7nvOBl9ySYv2GsIgOAzRLgnkuj51gJCWg46zThMacxqyI4OsyWTqfPwud/2kZTtwk9yCWt1QL/58P5dPX4aTjroz74096VuDvXmxR84Hw7RbkWMm0mMm2KLJuJTJuJLLuJEk8nhY2n+kq876yVN+rSuXlJL3kZ/Y7djHSjH3pp4YATfP8AOuQJo5/+2264IJO/+1AePc4AvY6A8dsZoMepw8tcAdo6/Ly5xxmzw8lY5WebmFNmZVaplVklFuOn1EqaRfF336/ncyu7uWFhjPROkDaZUKWFRoDNtke9x+H07b/r7EYuuLF14MkzlvIiqJw9OeacjxQ5AC3S7DKjt84EGe8BOO8Djmqta4NP/ihwPTBo4B6VkuBlr3OIPpqBQOyGo/Q042CIzEMW5RktxVPknnyRaRanz8TXX8/jR5d1sKTAx55mKw/XZNLQa6LZYcYXUMyfZeWX68ONqb9+pI2zDjNDxRCt4XSjd0BAc7oD/OihVt7Y5Ry2G+NIBTS8uHXoE/a1ldl89TyjdreiyMuKouirjbMuMy/U51Ddnol6x4/V0oLZBBazwmpR+Pyat3Y7+9ILPr8xbUFLhx+zyfjf79f4/DrYG1Hj94PXr3F7ApxpCUczYxbIXp59q98xOM5C7aiDPXbJWnvME0/oBri/2ZfJugoXOf3aNI51mHmu1oZ1RhGfvzg8WnhUaTx7gZFumDfTaOCrb4l9E22rxehqOglzxYDRo6q9K/ruNqFbE05S8QTumUDk6agOOH/cSzKrzPiBYF9UF/Q6jGDcE/wdqxcAGA1KkfE+xhDxVDdUj5Y1wZ/BjHbgUIjDpXl779BBO80Cd32uGJ/feL3In7ZOP83tPtq7Rxf1n6u1kW7W3HFOdFBo6jXx+wOZvHgiA19AAYM0LMYQ0PD2nsFrpIm2eE4aZYUWivLMFOWao34DfOYHDTDYrHeDfE6Rn3OX38T91dl858IuVPB+1l98OZ8DrRZAkXbCyY0b/H3bj6nHlMlkfOdKC43vakNL9Jzu5y1P2r0Z46KUES+27w/HmIy0+HsyJUE8gTtW9BtQcVNK3Q7cDlBRMcbx/CaT0b0qu19rs8cbDuK9zmBgdw0c/DDFgvZYjXXgUDzbhwLhYOmHoS7FTSZYONuYIc5IOQToDaUdnAG6e/385Ygdk4LPr+npC0Sf3lSIN5C8z7m80ExOppks+8AUz3v7ndSc8MQ82VnMsGRO2pD7atj9Hcfn9NrpDM4t87Bhrotnjtk40GodcvtxkWU3pgMOaKPNakbx5A7aIWlWY46X3cGpA4a68p8E4gncdUDkvatmAfX9V9JabwQ2gpHjHpfS9ZdmNYaoRw5TDwSMnXz0tNEzpbxYgnY/Y+nNMtbaejzPEQjA8XovP/hcccxGv1DQf/ywndnZPq6pNAKRN6Awm+H9q2x89IqcqBRH6O9n3uhm12F3zABqNsHaxelcf0k2ZrPCYgarWfX9/ceXunhzjzPm3eQsZnjfclvMwNfa6efhFzoHvUIZy74a6jlifc73vpfDve8NbHwe09DxeCyea/ykkrzscBfBRM6OOQ7iaZy0YDROXgGcwWic/ITWev9g20zq2QHFiIyq0apfMIvnOQbbNp7ZGAfrJZKsbcfyfsdjf4vUNK43C9Za+4A7gBeBGuCxoYK2mFoScX/PwbaNdcPluNI0ETeETva2Y32/4zGtgpj64urcrLXeBGxKcFnEJDSe9/ccTv+861jSBsnadqztCWPd32J6mBp95cSkNZYa5GiCYLK3lRqzmAipPZxQTHpjqUGOJQgma1upMYuJMP1uXSaEEJPQuDZOCiGEmFwkcAshRIpJSKpEKdUCnBzl5kXA2WHXmnhSrpGRco2MlGtkpmK55mit4xr5k5DAPRZKqe3x5nkmkpRrZKRcIyPlGpnpXi5JlQghRIqRwC2EEClmMgbujckuwCCkXCMj5RoZKdfITOtyTboctxBCiKFNxhq3EEKIISQlcCulPqaU2q+UCiilqvo99i2l1FGl1CGl1FWDbD9PKbVNKXVEKfVHpdS436oi+Ly7gj8nlFK7BlnvhFJqb3C9hA8XVUp9Tyl1JqJs1wyy3obgPjyqlPrmBJTrR0qpg0qpPUqpJ5RSMe9TNVH7a7j3r5RKD37GR4PH0txElSXiNWcrpV5VStUEj/8vxVjnMqVUZ8Tne2eiyxV83SE/F2X4aXB/7VFKnTMBZVocsR92KaW6lFJf7rfOhOwvpdSDSqlmpdS+iGUFSqmXgnHoJaVU/iDb3hpc54hS6tZxKZDWesJ/gKXAYuA1oCpi+TJgN5AOzAOOAeYY2z8G3BT8+xfA5xNc3nuBOwd57ARQNIH77nvA14ZZxxzcd5VAWnCfLktwudYDluDfdwN3J2t/xfP+gS8Avwj+fRPwxwn47MqBc4J/Z2PMc9+/XJcBz07U8RTv5wJcAzyPcUesC4BtE1w+M9CI0dd5wvcXcAlwDrAvYtk9wDeDf38z1jEPFAC1wd/5wb/zx1qepNS4tdY1WutDMR66HnhUa+3WWh8HjmLcrLiPUkoB64DHg4t+C9yQqLIGX+/jwCOJeo0E6LvBs9baA4Ru8JwwWuvN2pi7HWArxp2SkiWe9389xrEDxrF0RfCzThitdYPWekfw726M+e1nJvI1x9H1wO+0YSuQp5SayBm1rgCOaa1HO7BvTLTWW4C2fosjj6HB4tBVwEta6zatdTvwErBhrOWZbDnuWDcm7n9gFwIdEUEi1jrj6WKgSWt9ZJDHNbBZKVUdvO/mRLgjeLn64CCXZ/Hsx0S6DaN2FstE7K943n/fOsFjqRPj2JoQwdTMWmBbjIcvVErtVko9r5RaPkFFGu5zSfYxdRODV56Ssb8ASrXWDWCclIFYt4VPyH5L2LSuSqm/AmUxHvoXrfVTg20WY1n/bi9x3bw4HnGW8WaGrm1fpLWuV0qVAC8ppQ4Gz86jNlS5gP8B7sJ4z3dhpHFu6/8UMbYdc/ehePaXUupfAB/w8CBPM+77K1ZRYyxL2HE0UkqpLODPwJe11l39Ht6BkQ7oCbZfPAksnIBiDfe5JHN/pQHXAd+K8XCy9le8ErLfEha4tdZXjmKzeG5MfBbjMs0SrCnFvHnxeJRRGffb/DBw7hDPUR/83ayUegLjMn1MgSjefaeU+iXwbIyH4rrB83iXK9jw8kHgCh1M8MV4jnHfXzHE8/5D69QFP+dcBl4KjzullBUjaD+stf5L/8cjA7nWepNS6udKqSKtdULn5Yjjc0nIMRWnq4EdWuum/g8ka38FNSmlyrXWDcG0UXOMdeow8vAhszDa9sZksqVKngZuCrb4z8M4c74buUIwILwKfDS46FZgsBr8WF0JHNRa18V6UCmVqZTKDv2N0UC3L9a646VfXvFvB3m994CFyuh9k4Zxmfl0gsu1AfgGcJ3W2jHIOhO1v+J5/09jHDtgHEuvDHayGS/BHPqvgBqt9U8GWacslGtXSr0P4zvamuByxfO5PA18Oti75AKgM5QmmACDXvUmY39FiDyGBotDLwLrlVL5wbTm+uCysUl0a+wgLbR/i3EmcgNNwIsRj/0LRo+AQ8DVEcs3ATOCf1diBPSjwJ+A9ASV8zfA5/otmwFsiijH7uDPfoyUQaL33UPAXmBP8MAp71+u4P/XYPRaODZB5TqKkcvbFfz5Rf9yTeT+ivX+ge9jnFgAMoLHztHgsVQ5Afvo/RiXyXsi9tM1wOdCxxnGjbn3B/fRVuBvJqBcMT+XfuVSwH8H9+deInqDJbhsdoxAnBuxbML3F8aJowHwBmPXZzDaRF4GjgR/FwTXrQIeiNj2tuBxdhT4+/Eoj4ycFEKIFDPZUiVCCCGGIYFbCCFSjARuIYRIMRK4hRAixUjgFkKIFCOBWwghUowEbiGESDESuIUQIsX8f/kCrggLCS/NAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(pert_Deg, QM_Data, color=\"royalblue\", label='QM', marker='^',markersize=\"12\", linewidth=\"3\")\n", + "plt.plot(pert_Deg, MM_Data, color = \"mediumseagreen\", label='MM', marker ='o', markersize=\"6\", linewidth=\"3\")\n", + "plt.plot(pert_Deg, update, color = \"orange\", label='MM missing Imp', marker ='o', markersize=\"2\", linewidth=\"3\")\n", + "plt.plot(pert_Deg, subtract, color = \"pink\", label='QM - MM missing Imp', marker ='o', markersize=\"2\", linewidth=\"3\")\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([1.56469718]), array([[0.00019625]]))" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "def improper(k, theta):\n", + " return k*np.cos(theta)\n", + "opt.curve_fit(improper, pert_Deg, subtract, bounds=[-10,10])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.15" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/off_nitrogens/Update/restrainQM.py b/off_nitrogens/Update/restrainQM.py new file mode 100644 index 0000000..99d22b9 --- /dev/null +++ b/off_nitrogens/Update/restrainQM.py @@ -0,0 +1,5 @@ +import sys +sys.path.insert(0, '/export/home/jmaat/psi4_clone/01_scripts/') + +import confs2psi +confs2psi.confs2psi('/export/home/jmaat/off_nitrogens/QM_inputs/pyrnit.sdf','mp2','def2-SV(P)',False,"1.5 Gb") diff --git a/off_nitrogens/Update/run_multiple_directories.ssh b/off_nitrogens/Update/run_multiple_directories.ssh new file mode 100755 index 0000000..bf9643e --- /dev/null +++ b/off_nitrogens/Update/run_multiple_directories.ssh @@ -0,0 +1,6 @@ +#!/bin/bash + + + +for i in */; do (cp one_psi.slurm ./"$i"/1); done +for i in */ ; do (cd "$i"/1 && sbatch one_psi.slurm); done diff --git a/off_nitrogens/Update/run_perturb.py b/off_nitrogens/Update/run_perturb.py new file mode 100644 index 0000000..e173fd4 --- /dev/null +++ b/off_nitrogens/Update/run_perturb.py @@ -0,0 +1,14 @@ +#!/usr/bin/env python +# I m p o r t s +from perturb_angle import * + + + +# run the ~script!!~ +mol2_files = {3:"molClass_pyrnit_molecule_3.xyz", 4:"molClass_pyrnit_molecule_4.xyz", 5:"molClass_pyrnit_molecule_5.xyz", 6:"molClass_pyrnit_molecule_6.xyz", 7:"molClass_pyrnit_molecule_7.xyz", 8:"molClass_pyrnit_molecule_8.xyz", 9:"molClass_pyrnit_molecule_9.xyz"} + +input = input2mol(mol2_files) +print(input) +oemol_perturb(input, True, "pyrnit", 20, 1) + + diff --git a/off_nitrogens/calc_improper.py b/off_nitrogens/calc_improper.py index 7f34e6e..2482aec 100644 --- a/off_nitrogens/calc_improper.py +++ b/off_nitrogens/calc_improper.py @@ -36,6 +36,55 @@ # PRIVATE SUBROUTINES #============================================================================================= + +def translate(atom0, atom1, atom2, atom3, to_origin=True, old_central_atom=[]): + """ + Translate the central atom to the origin. Or, move it back + to its original position based on the old central atom's coordinates. + + Important for perturb_valence formula to shift valence angle of + outer atom around central atom, else shifts around incorrect origin + (evidenced by changed bond length). + + Parameters + ---------- + atom0 : numpy array + CENTRAL atom coordinates + atom1 : numpy array + outer atom coordinates + atom2 : numpy array + outer atom coordinates + atom3 : numpy array + outer atom coordinates + to_origin : Boolean + True to move central atom to origin and other atoms by same translation + False to reverse move from origin back to original central atom coordinates + old_central_atom : numpy array + CENTRAL atom coordinates of old system + + Returns + ------- + the four numpy arrays, translated: atom0, atom1, atom2, atom3 + + """ + if to_origin: + atom1 = atom1 - atom0 + atom2 = atom2 - atom0 + atom3 = atom3 - atom0 + atom0 = atom0 - atom0 # central must be moved last + else: + old_central_atom = np.asarray(old_central_atom) + if old_central_atom.shape[0] == 0: + # throw exception? (TODO) + return + atom1 = atom1 + old_central_atom + atom2 = atom2 + old_central_atom + atom3 = atom3 + old_central_atom + atom0 = atom0 + old_central_atom + + return atom0, atom1, atom2, atom3 + + def angle_between(v1, v2): """ Calculate the angle in degrees between vectors 'v1' and 'v2'. @@ -78,7 +127,8 @@ def calc_valence_angle(atom0, atom1, atom2): v2 = atom2-atom0 return(angle_between(v1, v2)) -def calc_improper_angle(atom0, atom1, atom2, atom3, translate=False): + +def calc_improper_angle(atom0, atom1, atom2, atom3): """ Calculate the improper dihedral angle of a set of given four atoms. @@ -92,20 +142,12 @@ def calc_improper_angle(atom0, atom1, atom2, atom3, translate=False): outer atom coordinates atom3 : numpy array outer atom coordinates - translate : bool - True to translate central atom to origin, False to keep as is. - This should not affect the results of the calculation. Returns ------- float Angle in degrees. """ - if translate: - atom1 = atom1 - atom0 - atom2 = atom2 - atom0 - atom3 = atom3 - atom0 - atom0 = atom0 - atom0 # central must be moved last # calculate vectors v0 = atom0-atom1 @@ -164,7 +206,7 @@ def find_improper_angles(mol): crd3 = np.asarray(mol_coords[nbors[2].GetIdx()]) # store coordinates crdlist.append([crd0, crd1, crd2, crd3]) - namelist.append([atom.GetName(), nbors[0].GetName(), nbors[1].GetName(),nbors[2].GetName()]) - #print(atom.GetName(),nbors[0].GetName(),nbors[1].GetName(),nbors[2].GetName()) + Idxlist.append([atom.GetIdx(), nbors[0].GetIdx(), nbors[1].GetIdx(),nbors[2].GetIdx()]) + - return crdlist, namelist + return crdlist, Idxlist diff --git a/off_nitrogens/perturb_angle.py b/off_nitrogens/perturb_angle.py index 1ed1dfd..cfd6d14 100644 --- a/off_nitrogens/perturb_angle.py +++ b/off_nitrogens/perturb_angle.py @@ -11,7 +11,7 @@ 1. Change valence angle without changing improper 2. Change improper angle without changing valence angles -For use in generating geometries to parameterize valence and improper angles +For use in generating geometries to parameterize valence and improper angles for an oemol By: Victoria Lim and Jessica Maat @@ -24,14 +24,18 @@ import numpy as np import math import sys +from openeye import oeomega +from openeye import oechem +from calc_improper import * +#from off_nitrogens.calc_improper import * -#from calc_improper import * -from off_nitrogens.calc_improper import * #============================================================================================= # PRIVATE SUBROUTINES #============================================================================================= + + def rotation_matrix(axis, theta): """ Return the rotation matrix associated with counterclockwise rotation about @@ -55,7 +59,6 @@ def rotation_matrix(axis, theta): """ axis = np.asarray(axis) theta = np.radians(theta) - axis = axis/math.sqrt(np.dot(axis, axis)) a = math.cos(theta/2.0) b, c, d = -axis*math.sin(theta/2.0) @@ -65,6 +68,10 @@ def rotation_matrix(axis, theta): [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)], [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]]) + + + + def perturb_valence(atom0, atom1, atom2, atom3, theta, verbose=False): """ Rotate atom3 in the plane of atom1, atom2, and atom3 by theta degrees. @@ -93,107 +100,226 @@ def perturb_valence(atom0, atom1, atom2, atom3, theta, verbose=False): rot_mat : numpy array three-dimension rotation matrix, returned so that the same matrix can be applied to any atoms attached to atom3. - """ # outer atoms are atom1 atom2 atom3. get normal vector to that plane. - v1 = atom2-atom1 - v2 = atom2 - atom3 + v1 = np.asarray(atom2)-np.asarray(atom1) + v2 = np.asarray(atom2) - np.asarray(atom3) w2 = np.cross(v1, v2) # calculate rotation matrix rot_mat = rotation_matrix(w2, theta) + length = np.linalg.norm(atom0-atom3) atom3_rot = np.dot(rot_mat, atom3) - # print details of geometry - if verbose: - print("\n>>> Perturbing valence while maintaining improper angle...") - # atom being moved is atom3. - print("\natom0 original coords:\t",atom0) - print("atom1 original coords:\t",atom1) - print("atom2 original coords:\t",atom2) - print("atom3 original coords:\t",atom3) - print("atom3 {:.2f} deg rotated: {}".format(60.,atom3_rot)) - - # check improper but make sure the central is first and moved is last - print("\nImproper angle, before: ", calc_improper_angle(atom0, atom1, atom2, atom3)) - print("Improper angle, before: ", calc_improper_angle(atom0, atom2, atom3, atom1)) - print("Improper angle, before: ", calc_improper_angle(atom0, atom3, atom1, atom2)) - print("Improper angle, after: ", calc_improper_angle(atom0, atom1, atom2, atom3_rot)) - print("Improper angle, after: ", calc_improper_angle(atom0, atom2, atom3_rot, atom1)) - print("Improper angle, after: ", calc_improper_angle(atom0, atom3_rot, atom1, atom2)) - - ## check valences - print("\nvalence angle 1, before: ", calc_valence_angle(atom0, atom1, atom2)) - print("valence angle 2, before: ", calc_valence_angle(atom0, atom1, atom3)) - print("valence angle 3, before: ", calc_valence_angle(atom0, atom2, atom3)) - print("valence angle 1, after: ", calc_valence_angle(atom0, atom1, atom2)) - print("valence angle 2, after: ", calc_valence_angle(atom0, atom1, atom3_rot)) - print("valence angle 3, after: ", calc_valence_angle(atom0, atom2, atom3_rot)) - print() + new_length = np.linalg.norm(atom0-atom3_rot) return atom0, atom1, atom2, atom3_rot -def oemol_perturb_valence(mol, central_atom, outer_atom, theta): + + + + + +def perturb_improper(atom0, atom1, atom2, atom3, theta, verbose=False): """ - From an OpenEye OEMol, specify the improper angle and the specific atom of - that improper that should be perturbed. The improper angles are obtained - from the find_improper_angles function in the calc_improper script. + Rotate atom3 out of the plane of atom1, atom2 by theta degrees while keeping valence angle the same. Parameters - --------- - mol : OpenEye OEMol - molecule from which to generate perturbed geometry - central_atom : string - atom name in the mol which is central to the improper of interest - Ex., "N1" - outer_atom : string - atom name in the mol which is to be rotated - Ex., "N1" + ---------- + atom0 : numpy array + CENTRAL atom coordinates + atom1 : numpy array + outer atom coordinates + atom2 : numpy array + outer atom coordinates + atom3 : numpy array + outer atom coordinates TO BE MOVED theta : float - how many degrees by which to rotate + how many degrees by which to move atom upwards + + Returns + ------- + atom* : numpy arrays + the coordinates in the same order as given in parameters + rot_mat : numpy array + three-dimension rotation matrix, returned so that the same + matrix can be applied to any atoms attached to atom3. + + """ + + # get the vector between atom2 and atom1 and then rotate atom3 about that plane by angle theta. + v1 = atom2-atom1 + # calculate rotation matrix + rot_mat = rotation_matrix(v1, theta) + atom3_rot = np.dot(rot_mat, atom3) + + return atom0, atom1, atom2, atom3_rot - [TODO] +def input2mol(xyz): """ - # todo [1] - # calc_improper.py: change find_improper_angles to return the name of the central atom as 5th element in crdlist - # - under aidx, add something like: aname = atom.GetName() - # - then update this line: crdlist.append((crd0, crd1, crd2, crd3, aname)) - # - update returns section in docstring - # - make sure tests are still passing - - # call find_improper_angle function to get coordinates for some specified improper angle - # todo [2] - - # call perturb_valence on those coordinates - # todo [3] - - # from the MATRIX returned by perturb_valence, update coordinates using your new function of todo [4] - # note: this means you probably won't need the COORDINATES returned by the perturb_valence function - # todo [5] - - return # placeholder - -# todo [4] -# write a new function to set the new coordinates back to the OEMol -# def update_oemol_coordinates(mol, moved_atom, move_matrix) -# 1. use the name moved_atom to get the OEAtom (atom_moved) -# 2. loop over all atoms connected to moved_atom -# 3. apply the move_matrix on those coordinates: something like -# -# for connected_atom in ___: -# old_coords = connected_atom.GetCoords() -# connected_atom.SetCoords(np.dot(move_matrix, old_coords)) -# -# check the syntax though, I'm just writing pseudo-code -# + d e s c r i p t i o n : + Takes an xyz file dictionary and converts the coordinates to an eomol. + Function returns a dictionary of oemols. + + p a r a m e t e r s : + xyz: dictionary of xyz files + -def perturb_improper(atom0, atom1, atom2, atom3, theta, verbose=False): """ - [TODO] [7] + ifs = oechem.oemolistream() + for key, value in xyz.items(): + ifs.open(value) + for mol in ifs.GetOEGraphMols(): + xyz[key] = oechem.OEGraphMol(mol) + return xyz + +""" +ifs = oechem.oemolistream() +ofs = oechem.oemolostream() + +if ifs.open("molClass_pyrnit_molecule_1.xyz"): + if ofs.open("output.mol2"): + for mol in ifs.GetOEGraphMols(): + oechem.OEWriteMolecule(ofs, mol) +""" +def smileslist2mol(smilesList): + """ + This function creates a dictionary of oemol with keys that are """ - # todo [6] - return # placeholder + mol = oechem.OEMol() + oechem.OESmilesToMol(mol, smiles) + oechem.OEAddExplicitHydrogens(mol) + omega = oeomega.OEOmega() + omega.SetMaxConfs(1) + omega(mol) + + + +def oemol_perturb(molList, angle_type, molClass, pertRange, pertIncr): + """ + The function takes in a list of smiles strings, converts the smiles strings into OpenEye + OEMols and then performs either a valence or improper perturbation to the molecules geometry + by a specified angle "theta". The function creates an output .sdf file which contains all the + molecules from the molList and a SD tag that contains the the indices of the atoms that + are in the improper molecule center. The first index is the center of the improper and the last + index indicates the atom that was perturbed in the geometry change. The function utilizes + find_improper_angles from the calc_improper.py script to identify the improper locations in the molecule. + + The code will iterate through each trivalent center, and perturb the centers individually in the output + .sdf file. The code will also generate 3 perturbed geometries for each nitrogen center where each of the + individual constiutuents will be perturbed individually. + + + Parameters + --------- + molList : Dictionary of oemol objects + angle_tyle: + True: Boolean + True = Improper perturbation + False = Valence perturbation + Ex., True + molClass : String + Type of molecules + Ex., pyrnitrogens + PertRang : Int + The range of degrees the attached atom will be perturbed by + pertIncr: Int + The increment of degrees that the molecule with be perturbed by + + Returns + -------- + .sdf file for each molecule in smiles string that perturbs the improper or valence by increments specified in specified + range with tags that include the frozen indices of the molecule + + """ + + #Set perturbation type + if angle_type == True: + angle_type = perturb_improper + perturbation = "improper" + else: + angle_type = perturb_valence + perturbation = "valence" + #convert molList into OEMols + for key, mol in molList.items(): + impCent = find_improper_angles(mol) + centcount = 0 + constcount = 0 + for center in impCent[1]: + for constituent in center: + constcount +=1 + #determine which improper angle on the atom is of interest and store the coordinates of the improper in various variables + #cmol is the parent mol which we will add conformers to + cmol = oechem.OEMol(mol) + mol_pert = str(constituent) + center_atom = cmol.GetAtom(oechem.OEHasAtomIdx(center[0])) + if constituent == center[0]: + continue + move_atom = cmol.GetAtom(oechem.OEHasAtomIdx(constituent)) + center_coord = np.array(cmol.GetCoords(center_atom)) + move_coord = np.array(cmol.GetCoords(move_atom)) + + + #center_coord = center[0] + #move_atom = constituent + other_coords = list() + + #if constituent != center_coord: + # if constituent != move_atom: + # other_coords.append(constituent) + + + for neighbor in center_atom.GetAtoms(): + if neighbor.GetIdx() != constituent: + new_coord = np.array(cmol.GetCoords(neighbor)) + other_coords.append(new_coord) + + # rotate atom by desired increment, and write out each perturbation to the .sdf file + # DEBUGGING + theta = 0 + print(pertRange) + + oemol_list = [cmol] + + #in the list adding tags for the indices around the nitrogen center, improver or valence move, and angle of perturbation + ofile = oechem.oemolostream(str(molClass) + "_" + str(key) + "_constituent_" + mol_pert +"_" + perturbation +'.sdf') + count = 1 + theta = (360 - (pertRange/2)) + print("starting theta:" + str(theta)) + maxRange = ((pertRange/2)+360) + print("Max range:" + str(maxRange)) + pertTheta = -(pertRange/2) + while theta < maxRange: + #move in direction of theta + atom0, atom1, atom2, atom3_rot = angle_type(center_coord, other_coords[0], other_coords[1], move_coord, theta) + move_mol = oechem.OEMol(oemol_list[0]) + move_mol.SetCoords(move_atom, oechem.OEFloatArray(atom3_rot)) + move_mol.SetTitle(str(molClass) + "_" + str(key) + "_constituent_" + mol_pert +"_" + perturbation) + oechem.OEAddSDData(move_mol, "Index list of atoms to freeze", str(center)) + oechem.OEAddSDData(move_mol, "Angle OEMol is perturbed by", str(pertTheta)) + oechem.OEAddSDData(move_mol, "Type of perturbation (improper or valence)", perturbation) + oemol_list.append(move_mol) + oechem.OEWriteConstMolecule(ofile, move_mol) + + #seperate .mol2 for each perturbation + mfile = oechem.oemolostream(str(molClass) + "_" + str(key) + "_constituent_" + mol_pert + "_" + perturbation + '_angle_'+ str(pertTheta) +'.mol2') + oechem.OEWriteConstMolecule(mfile, move_mol) + mfile.close() + + + + #count for while loop + theta += pertIncr + count += 1 + pertTheta += pertIncr + + + ofile.close() + print("end of loop") + + return + + -# todo [8] write a test for your function diff --git a/off_nitrogens/tests/test_improper.py b/off_nitrogens/tests/test_improper.py index 07fb0d5..e2afc5d 100644 --- a/off_nitrogens/tests/test_improper.py +++ b/off_nitrogens/tests/test_improper.py @@ -1,5 +1,6 @@ from off_nitrogens.calc_improper import * #from calc_improper import * + import numpy as np def test_two_vectors():