Autosolvate add orca

autosolvate/autosolvate.py

@@ -46,14 +46,15 @@ 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=""):
         self.xyz = xyzfile
         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, 

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']

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)

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():

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)