From ef5e9a0859c8fc32d029b9284113692c2729e322 Mon Sep 17 00:00:00 2001 From: SangniXun Date: Wed, 18 Sep 2024 10:57:51 -0400 Subject: [PATCH 1/4] update conda env for orca esp --- devtools/conda-envs/doc_env.yaml | 1 + devtools/conda-envs/test_env.yaml | 3 ++- devtools/conda-envs/test_env_WIN.yaml | 1 + 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/devtools/conda-envs/doc_env.yaml b/devtools/conda-envs/doc_env.yaml index 81fd688..449036b 100644 --- a/devtools/conda-envs/doc_env.yaml +++ b/devtools/conda-envs/doc_env.yaml @@ -30,5 +30,6 @@ dependencies: - imolecule - pillow - pubchempy + - git+https://github.com/Liu-group/pyvdwsurface_moreelements.git # - codecov diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 3e14168..0e42d96 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -1,4 +1,4 @@ -name: autosolvate +name: autosolvate_add_orca channels: - conda-forge dependencies: @@ -36,5 +36,6 @@ dependencies: - imolecule - pillow - pubchempy + - git+https://github.com/Liu-group/pyvdwsurface_moreelements.git # - codecov diff --git a/devtools/conda-envs/test_env_WIN.yaml b/devtools/conda-envs/test_env_WIN.yaml index 4c312d5..8183599 100644 --- a/devtools/conda-envs/test_env_WIN.yaml +++ b/devtools/conda-envs/test_env_WIN.yaml @@ -33,5 +33,6 @@ dependencies: - imolecule - pillow - pubchempy + - git+https://github.com/Liu-group/pyvdwsurface_moreelements.git # - codecov From f54fdfa4bb584f0d6c222f39668b550d296789a9 Mon Sep 17 00:00:00 2001 From: SangniXun Date: Wed, 18 Sep 2024 11:02:13 -0400 Subject: [PATCH 2/4] update autosolvate orca resp charge --- autosolvate/autosolvate.py | 19 +- autosolvate/globs.py | 4 +- autosolvate/resp_classes/gen_esp.py | 68 +++++++ autosolvate/resp_classes/resp_abstract.py | 2 + autosolvate/resp_classes/resp_factory.py | 9 +- autosolvate/resp_classes/resp_orca.py | 217 ++++++++++++++++++++++ 6 files changed, 308 insertions(+), 11 deletions(-) create mode 100644 autosolvate/resp_classes/gen_esp.py create mode 100644 autosolvate/resp_classes/resp_orca.py diff --git a/autosolvate/autosolvate.py b/autosolvate/autosolvate.py index 88f8fed..3a337ca 100644 --- a/autosolvate/autosolvate.py +++ b/autosolvate/autosolvate.py @@ -46,7 +46,7 @@ class solventBoxBuilder(): def __init__(self, xyzfile, solvent='water', slu_netcharge=0, cube_size=54, - charge_method="resp", slu_spinmult=1, outputFile="", + charge_method="resp", slu_spinmult=1, outputFile="", nprocs='8', srun_use=False, qm_program='gaussian', qm_exe=None, qm_dir=None, amberhome=None, closeness=0.8, solvent_off="", solvent_frcmod="",rundir=""): @@ -54,6 +54,7 @@ def __init__(self, xyzfile, solvent='water', slu_netcharge=0, cube_size=54, self.solute = pybel.readfile('xyz', xyzfile).__next__() self.slu_netcharge = slu_netcharge self.slu_spinmult = slu_spinmult + self.nprocs = nprocs # currently hard coded. Can be changed later to be determined automatically based on the density of the solute self.solvent = solvent if closeness=='automated': @@ -245,8 +246,8 @@ def getFrcmod(self): self.removeConectFromPDB() if self.charge_method == "resp": - myresp = resp_factory(pdbfile="solute.xyz.pdb", charge=self.slu_netcharge, - spinmult=self.slu_spinmult, qm_program=self.qm_program, + myresp = resp_factory(pdbfile="solute.xyz.pdb", charge=self.slu_netcharge,xyzfile = self.xyz, + spinmult=self.slu_spinmult, qm_program=self.qm_program, nprocs = self.nprocs, qm_exe=self.qm_exe, qm_dir=self.qm_dir,srun_use=self.srun_use,rundir=self.rundir) myresp.run() @@ -653,8 +654,8 @@ def startboxgen(argumentList): Generates the structure files and save as ```.pdb```. Generates the MD parameter-topology and coordinates files and saves as ```.prmtop``` and ```.inpcrd``` """ #print(argumentList) - options = "hm:n:s:o:c:b:g:u:rq:e:d:a:t:l:p:vD:" - long_options = ["help", "main","solutename", "solvent", "output", "charge", "cubesize", "chargemethod", "spinmultiplicity", "srunuse","qmprogram","qmexe", "qmdir", "amberhome", "closeness","solventoff","solventfrcmod", "validation", "runningdirectory"] + options = "hm:n:s:o:c:b:g:u:rq:e:d:a:t:l:p:vD:y:" + long_options = ["help", "main","solutename", "solvent", "output", "charge", "cubesize", "chargemethod", "spinmultiplicity", "srunuse","qmprogram","qmexe", "qmdir", "amberhome", "closeness","solventoff","solventfrcmod", "validation", "runningdirectory","nprocs"] arguments, values = getopt.getopt(argumentList, options, long_options) solutename = "" solutexyz="" @@ -675,6 +676,7 @@ def startboxgen(argumentList): solvent_off="" solvent_frcmod="" rundir = "" + nprocs = 8 for currentArgument, currentValue in arguments: if currentArgument in ("-h", "--help"): print('Usage: autosolvate boxgen [OPTIONS]') @@ -696,7 +698,9 @@ def startboxgen(argumentList): print(' -p, --solventfrcmod path to the custom solvent .frcmod file') print(' -n, --solutename initialize suggested parameter for given solute') print(' -v, --validation option to run validation step for input parameters') + print(' -y, --nprocs procs to run orca resp charge calculation, if -q orca') print(' -h, --help short usage description') + print exit() elif currentArgument in ("-m", "--main"): print ("Main/solutexyz", currentValue) @@ -722,6 +726,9 @@ def startboxgen(argumentList): elif currentArgument in ("-s", "--solvent"): print ("Solvent:", currentValue) solvent=str(currentValue) + elif currentArgument in ("-y","--nprocs"): + nprocs = str(currentValue) + print('Use',nprocs,'proces to run orca respfitting if -g orca') elif currentArgument in ("-o", "--output"): print ("Output:", currentValue) outputFile=str(currentValue) @@ -794,7 +801,7 @@ def startboxgen(argumentList): exit() builder = solventBoxBuilder(solutexyz, solvent=solvent, slu_netcharge=slu_netcharge, - cube_size=cube_size, charge_method=charge_method, + cube_size=cube_size, charge_method=charge_method, nprocs=nprocs, slu_spinmult=slu_spinmult, outputFile=outputFile, srun_use=srun_use, qm_program=qmprogram, qm_exe=qmexe, qm_dir=qmdir, amberhome=amberhome, closeness=closeness, diff --git a/autosolvate/globs.py b/autosolvate/globs.py index 516afc9..fbc4490 100644 --- a/autosolvate/globs.py +++ b/autosolvate/globs.py @@ -2,8 +2,8 @@ Global variables of AutoSolvate ''' -available_qm_programs = ["gaussian", "gamess", "terachem"] -keywords_avail = ["molecule", "pdbfile", "molname", +available_qm_programs = ["gaussian", "gamess", "terachem","orca"] +keywords_avail = ["molecule", "pdbfile", "molname", "nprocs","xyzfile", "qm_program", "qm_exe", "qm_dir", "charge", "spinmult", "rundir","srun_use","nnodes","ncpus","gamessversion","rundir"] available_charge_methods = ['resp','bcc'] diff --git a/autosolvate/resp_classes/gen_esp.py b/autosolvate/resp_classes/gen_esp.py new file mode 100644 index 0000000..49af153 --- /dev/null +++ b/autosolvate/resp_classes/gen_esp.py @@ -0,0 +1,68 @@ +import numpy as np +import subprocess +from pyvdwsurface import vdwsurface + +AtoB = 1/0.52917721092 + +def read_xyz(xyz): + crds = [] + elements = [] + ifile = open(xyz,'r') + data = ifile.readlines() + ifile.close() + for line in data[2:]: + ele = str(line.split()[0]) + x = float(line.split()[1]) + y = float(line.split()[2]) + z = float(line.split()[3]) + crds.append([x,y,z]) + elements.append(ele.capitalize()) + # print(crds, elements) + return crds, elements + +def gen_grids(crds,elements,orcapath,gbw,denisty,out): + gridall = [] + atoms = np.array(crds,dtype=float) + elements_bytes = [element.encode('utf-8') for element in elements] + points_1 = vdwsurface(atoms, elements=elements_bytes, density=1,scale_factor = 2.0) + new_point_1 = [[x *AtoB for x in sublist] for sublist in points_1] + gridall.extend(new_point_1) + + points_2 = vdwsurface(atoms, elements=elements_bytes, density=1,scale_factor = 1.4) + new_point_2= [[x *AtoB for x in sublist] for sublist in points_2] + gridall.extend(new_point_2) + + points_3 = vdwsurface(atoms, elements=elements_bytes, density=1,scale_factor = 1.6) + new_point_3= [[x *AtoB for x in sublist] for sublist in points_3] + gridall.extend(new_point_3) + + points_4 = vdwsurface(atoms, elements=elements_bytes, density=1,scale_factor = 1.8) + new_point_4= [[x *AtoB for x in sublist] for sublist in points_4] + gridall.extend(new_point_4) + + # print(gridall) + + with open('esp.xyz','w') as f: + f.write(str(len(gridall))+'\n') + for line in gridall: + for crd in line: + f.write(str(crd) + ' ') + f.write('\n') + cmd = orcapath + ' ' + gbw + ' ' + denisty + ' esp.xyz ' + out + print(cmd) + with open('esp_gen.log','w') as f: + subprocess.call(cmd,shell=True,stdout=f, stderr=subprocess.STDOUT) + +def convertoesp(espin,espout,crds): + ifile = open(espin,'r') + data = ifile.readlines() + ifile.close() + length_crds = len(crds) + length_esp = int(data[0]) + with open(espout,'w') as f: + print("%5d%5d%5d" %(length_crds, length_esp, 0), file=f) + for i in crds: + print("%16s %15.7E %15.7E %15.7E" %(' ', float(i[0])*AtoB, float(i[1])*AtoB, float(i[2])*AtoB), file=f) + for i in data[1:]: + line = i.split() + print("%16.7E %15.7E %15.7E %15.7E" %(float(line[3]), float(line[0]), float(line[1]), float(line[2])), file=f) diff --git a/autosolvate/resp_classes/resp_abstract.py b/autosolvate/resp_classes/resp_abstract.py index 1543781..cb88b24 100644 --- a/autosolvate/resp_classes/resp_abstract.py +++ b/autosolvate/resp_classes/resp_abstract.py @@ -28,6 +28,8 @@ def __init__(self, **kwargs): self.rundir = os.path.abspath(self.rundir) self.resp_scr_dir = os.path.join(os.path.abspath(self.rundir), 'resp_scr') self.initialization_check(**kwargs) + self.xyzfile = kwargs["xyzfile"] if "xyzfile" in kwargs.keys() else "undef" + self.nprocs = kwargs["nprocs"] if "nprocs" in kwargs.keys() else "undef" def initialization_check(self, **kwargs): for key, value in kwargs.items(): diff --git a/autosolvate/resp_classes/resp_factory.py b/autosolvate/resp_classes/resp_factory.py index 0d020b0..de88e66 100644 --- a/autosolvate/resp_classes/resp_factory.py +++ b/autosolvate/resp_classes/resp_factory.py @@ -6,6 +6,8 @@ """ from autosolvate.resp_classes.resp_gaussian import RespGaussian from autosolvate.resp_classes.resptools.resp_gamess import RespGAMESS +from autosolvate.resp_classes.resp_orca import RespORCA +from autosolvate.resp_classes.resp_orca import get_crds_from_xyz def resp_factory(**kwargs): """ @@ -19,21 +21,22 @@ def resp_factory(**kwargs): qm_program: str Quantum chemistry pacakge used to calculat RESP - Currently support: ['gamess','gaussian'] + Currently support: ['gamess','gaussian','orca'] Returns ------- cls_instance: instance of the class that performs RESP fitting with the request quantum chemistry program """ - cls_dict = dict(gaussian=RespGaussian, gamess=RespGAMESS) + cls_dict = dict(gaussian=RespGaussian, gamess=RespGAMESS, orca=RespORCA) qm_program = kwargs["qm_program"] if "qm_program" in kwargs.keys() else "gamess" if qm_program not in cls_dict.keys(): raise Exception("The request quantum chemistry program, {}, is not supported.".format(qm_program)) - + if "pdbfile" not in kwargs.keys(): raise KeyError("Error: funciton resp_factory requeests a pdbfile file in the keyword argument for resp charge fitting.") + cls = cls_dict[qm_program] cls_instance = cls(**kwargs) diff --git a/autosolvate/resp_classes/resp_orca.py b/autosolvate/resp_classes/resp_orca.py new file mode 100644 index 0000000..3b6fe9e --- /dev/null +++ b/autosolvate/resp_classes/resp_orca.py @@ -0,0 +1,217 @@ +import autosolvate.resp_classes.gen_esp as gen_esp +from autosolvate.resp_classes.resp_abstract import RespABC +import os +import subprocess +#import os, re, subprocess, shutil, glob, sys +def get_crds_from_xyz(inp): + crds = [] + with open(inp,'r') as f: + data = f.readlines() + for line in data[2:]: + if len(line) > 3: + line = line.strip() + atom = line.split()[0].capitalize() + x = line.split()[1] + y = line.split()[2] + z = line.split()[3] + crds.append([atom,x,y,z]) + return crds + + +class RespORCA(RespABC): + def __init__(self, **kwargs): + print("*"*40) + print("Run Orca to generate RESP charge and mol2 file".center(40," ") ) + print("*"*40) + self.srun_use = kwargs["srun_use"] if "srun_use" in kwargs else False + # TODO: read in srun_use option from kwargs + RespABC.__init__(self, **kwargs) + + if self.qm_exe == None: + print("WARNING: Orca executable name is not specified for RESP charge fitting!") + print("WARNING: Using orca by default. If failed later,\ + please rerun with the option -e specified!") + self.qm_exe = 'orca' + if self.qm_dir == None: + print("WARNING: orca executable directory is not specified for RESP charge fitting!") + print("WARNING: Setting to default path: /opt/orca/5.0.2/") + print("WARNING: If failed later, please rerun with the option -d specified!") + self.qm_dir = '/opt/orca/5.0.2/' + + self.orcapath = os.path.join(self.qm_dir, self.qm_exe) + + def writeOrcaInput(self): + """ + Set up Orca calculation to compute electrostatic potential. + + optional arguments: + calculation - one of 'optimize' (hours), or 'singlepoint' (minutes) (defaults to 'singlepoint') + + """ + print("-"*40) + print("Preparing Orca input files...".center(40," ")) + print("-"*40) + + charge = self.molecule.GetTotalCharge() + multiplicity = self.molecule.GetTotalSpinMultiplicity() + + orca_inp = self.molname + "_orca.in" + orca_out = self.molname + '_orca.out' + xyzfile = self.xyzfile + orca_crd = get_crds_from_xyz(xyzfile) + + with open(orca_inp,'w') as ofile: + ofile.write('! b3lyp 6-31g* keepdens\n\n') + ofile.write("%pal\nnprocs " + str(self.nprocs) + "\nend\n\n") + ofile.write('* xyz ' + str(charge) + ' ' + str(multiplicity) + '\n') + for line in orca_crd: + atomname = line[0] + x = line[1] + y = line[2] + z = line[3] + ofile.write(atomname + ' ' + x + ' ' + y + ' ' + z + '\n') + ofile.write('*') + + cmd1 = self.orcapath + ' ' + orca_inp + ' > ' + orca_out + + ## check convergence of orcaout ## + convergence = False + if os.path.exists(orca_out): + with open(orca_out) as f: + data = f.read() + if 'ORCA TERMINATED NORMALLY' in data: + print('Orca calculation is finished for esp...') + convergence = True + + if convergence == False: + if self.srun_use: + cmd1='srun -n ' + self.nprocs + ' ' +cmd1 + print(cmd1) + subprocess.call(cmd1, shell=True) + + def generateESP(self): + vpotpath = self.orcapath + '_vpot' + orca_inp = self.molname + "_orca.in" + orca_out = self.molname + '_orca.out' + gbw = self.molname + '_orca.gbw' + density = self.molname + '_orca.scfp' + espf = self.molname + '_orca.esp' + out = self.molname + '_orca.espout' + crds = get_crds_from_xyz(self.xyzfile) + elements = [] + newcrds = [] + for line in crds: + newcrds.append([float(line[1]),float(line[2]),float(line[3])]) + elements.append(line[0].capitalize()) + + gen_esp.gen_grids(crds=newcrds,elements=elements,orcapath=vpotpath,gbw=gbw,denisty=density,out=out) + gen_esp.convertoesp(espin=out,espout=espf,crds=newcrds) + + # os.system("resp -O -i resp1.in -o resp1.out -p resp1.pch -t resp1.chg \ + # -e %s -s resp1_calc.esp" %espf) + # os.system("resp -O -i resp2.in -o resp2.out -p resp2.pch -q resp1.chg \ + # -t resp2.chg -e %s -s resp2_calc.esp" %espf) + + + + def respFit(self): + Method1 = True + print("-"*40) + print("Start RESP fitting with ESP generaged by orca".center(40," ") ) + print("-"*40) + + ac_filename = self.molname + ".ac" + print("Generating Antechamber file: " + ac_filename) + + cmd = "antechamber -i {pdbfile} -fi pdb -o {ac} -fo ac -c dc -nc {charge} -pl 30".format(pdbfile=self.pdbfile, ac=ac_filename, charge=self.charge) + print(cmd) + with open('antechamber_acfile.log', 'w') as f: + subprocess.call(cmd, shell=True, stdout=f, stderr=subprocess.STDOUT) + + resp1_filename = self.molname + '.resp1' + print("Using respgen to generate stage 1 RESP fitting input file: " + resp1_filename) + + cmd = "respgen -i {ac} -o {resp1} -f resp1 -l 10".format(ac=ac_filename, resp1=resp1_filename) + print(cmd) + subprocess.call(cmd,shell=True) + + # Stage 2 + resp2_filename = self.molname + '.resp2' + print("Using respgen to generate stage 2 RESP fitting input file: " + resp2_filename) + + cmd = "respgen -i {ac} -o {resp2} -f resp2 -l 10".format(ac=ac_filename, resp2=resp2_filename) + print(cmd) + subprocess.call(cmd,shell=True) + + esp_filename = self.molname + '_orca.esp' + + # Perform RESP fit. + #Stage 1 + qout1_filename = self.molname + '.qout_stage1' + respout1_filename = self.molname + '.respout1' + print("Using resp to perform stage 1 RESP fitting to generate: " + qout1_filename) + + cmd = "resp -O -i {resp1} -o {respout1} -e {esp} -t {qout1}".format( + resp1=resp1_filename, respout1=respout1_filename, + esp=esp_filename, qout1=qout1_filename) + print(cmd) + subprocess.call(cmd,shell=True) + + #Stage 2 + qout2_filename = self.molname + '.qout_stage2' + respout2_filename = self.molname + '.respout2' + print("Using resp to perform stage 2 RESP fitting to generate: " + qout2_filename) + + cmd = "resp -O -i {resp2} -o {respout2} -e {esp} -q {qout1} -t {qout2}".format( + resp2=resp2_filename, respout2=respout2_filename, esp=esp_filename, + qout1=qout1_filename, qout2=qout2_filename) + print(cmd) + subprocess.call(cmd,shell=True) + + # Write molecule with updated charges to mol2 file. + mol2_filename = self.molname + ".mol2" + print("Writing out the mol2 file with resp charge: " + mol2_filename) + cmd = "antechamber -i {ac} -fi ac -o {mol2} -fo mol2 -c rc -cf {qout2} -pl 30".format( + ac=ac_filename, mol2=mol2_filename, qout2=qout2_filename) + print(cmd) + subprocess.call(cmd,shell=True) + + + # with open('resp1.in','w') as fresp1: + # fresp1.write() + + def run(self): + # if not os.path.isdir(self.resp_scr_dir): + # print("Creating the scratch folder for RESP fitting: ", self.resp_scr_dir) + # os.mkdir(self.resp_scr_dir) + + # print("Copying over the pdb file", self.pdbfile, " to ", self.resp_scr_dir) + # shutil.copy(os.path.join(self.rundir,self.pdbfile), self.resp_scr_dir) + + # print("Navigating to the scratch folder for RESP fitting: ", self.resp_scr_dir) + # os.chdir(self.resp_scr_dir) + + self.writeOrcaInput() + self.generateESP() + self.respFit() + + # print("Navigating back to the folder for AutoSolvate run: ", self.rundir) + # os.chdir(self.rundir) + + # mol2_filename = self.molname + ".mol2" + # print("Copying the generated mol2 file ", mol2_filename, " to rundir", self.rundir) + # shutil.copy(os.path.join(self.resp_scr_dir, mol2_filename), self.rundir) + + + + + + + + + + + + + + \ No newline at end of file From ca9d410ae4a4a4936a715d3b0302f890b21d35ee Mon Sep 17 00:00:00 2001 From: SangniXun <111337156+SangniXun@users.noreply.github.com> Date: Thu, 19 Sep 2024 13:56:33 -0400 Subject: [PATCH 3/4] Update test_env.yaml --- devtools/conda-envs/test_env.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 0e42d96..62d7c00 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -6,6 +6,7 @@ dependencies: - python=3.8 - pip - setuptools + - cython # Testing - pytest From ad124c0372fe6313212319094d038a5c19492d7d Mon Sep 17 00:00:00 2001 From: SangniXun <111337156+SangniXun@users.noreply.github.com> Date: Thu, 19 Sep 2024 16:04:58 -0400 Subject: [PATCH 4/4] Update test_env.yaml --- devtools/conda-envs/test_env.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/devtools/conda-envs/test_env.yaml b/devtools/conda-envs/test_env.yaml index 62d7c00..c296541 100644 --- a/devtools/conda-envs/test_env.yaml +++ b/devtools/conda-envs/test_env.yaml @@ -1,4 +1,4 @@ -name: autosolvate_add_orca +name: autosolvate channels: - conda-forge dependencies: