Merge pull request #68 from TimotheMelin/add_pyscf_implicit_solvent

Use PYSCF for PCM calculations

Author: leucinw

PoltypeModules/optimization.py

@@ -715,11 +715,23 @@ def GeometryOptimization(poltype,mol,totcharge,suffix='1',loose=False,checkbonds
     comoptfname=poltype.comoptfname.replace('_1','_'+suffix)
     chkoptfname=poltype.chkoptfname.replace('_1','_'+suffix)
 
-    if (poltype.use_gaus==True or poltype.use_gausoptonly==True): # try to use gaussian for opt
+    # Define which software will be used for the QM optimization
+    if not poltype.use_gaus or not poltype.use_gausoptonly:
+        if poltype.optpcm==True or (poltype.optpcm==-1 and poltype.pcm):
+            Soft = 'PySCF'
+        else:
+            Soft = 'Psi4'
+    else:
+        Soft = 'Gaussian'
+
+
+    if Soft == 'Gaussian': # try to use gaussian for opt
         term,error=poltype.CheckNormalTermination(logoptfname,errormessages=None,skiperrors=True)
         if not term or overridecheckterm==True:
 
             poltype.WriteToLog("QM Geometry Optimization")
+            poltype.WriteToLog(f'The software {Soft} will be used to do the QM optimization')
+
             gen_optcomfile(poltype,comoptfname,poltype.numproc,poltype.maxmem,poltype.maxdisk,chkoptfname,mol,modred,torsionrestraints)
             cmdstr = 'GAUSS_SCRDIR=' + poltype.scrtmpdirgau + ' ' + poltype.gausexe + " " + comoptfname
             jobtooutputlog={cmdstr:os.getcwd()+r'/'+logoptfname}
@@ -749,15 +761,103 @@ def GeometryOptimization(poltype,mol,totcharge,suffix='1',loose=False,checkbonds
         optmol =  load_structfile(poltype,logoptfname)
         optmol=rebuild_bonds(poltype,optmol,mol)
 
-           
-    else:
+    if Soft == 'PySCF':
+
+        import importlib
+            
+        try:
+            importlib.import_module('pyscf')
+            poltype.WriteToLog("PySCF is installed.\n")
+        except ImportError:
+            poltype.WriteToLog("PySCF is not installed.\n Make sure that pyscf version 2.7.0 is installed along with dft3, dft4 and pyscf-dispersion.\n You can install using: \n        `pip install pyscf==2.7.0 dft3 dft4 pyscf-dispersion`")
+            raise ImportError("PySCF is not installed.\n Make sure that pyscf version 2.7.0 is installed along with dft3, dft4 and pyscf-dispersion.\n You can install using: \n        `pip install pyscf==2.7.0 dft3 dft4 pyscf-dispersion`")
+    
+
+        gen_optcomfile(poltype,comoptfname,poltype.numproc,poltype.maxmem,poltype.maxdisk,chkoptfname,mol,modred,torsionrestraints)
+        term,error=poltype.CheckNormalTermination(logoptfname,errormessages=None,skiperrors=True)
+        modred=False
+
+        # Import both PySCF setup class
+        from pyscf_setup import PySCF_init_setup
+        from pyscf_setup import PySCF_post_run
+
+        # Instantiate the initial pyscf setup object
+        Opt_prep = PySCF_init_setup(poltype,os.getcwd(),comoptfname,mol,\
+                               skipscferror,charge)
+    
+        # Read the gaussian input coordinate
+        Opt_prep.read_Gauss_inp_coord()
+
+        # Write the initial xyz file to be used by PySCF
+        Opt_prep.write_xyz('init')
+
+        # Write the torsion constraints to be used 
+        # during the QM optimization
+        Opt_prep.write_geometric_tor_const('init',torsionrestraints)
+
+        # Write the PySCF input file
+        Opt_prep.write_PySCF_input(True,False)
+
+        # Define the cmd to run PySCF
+        cmd = f'python {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_inp_file} &> {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}'
+
+        # Create the dictionnary with the commnd to be used by CallJobsSeriallyLocalHost
+        cmd_to_run = {cmd: f'{Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}'}
+
+        # If the user only want to create input file
+        if poltype.checkinputonly==True:
+            sys.exit()
+
+        # Make sure to delete the PySCF output file 
+        # before running it
+        if os.path.isfile(f'{Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}'):
+            os.remove(f'{Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}')
+
+        # Run PySCF with either external API 
+        # or SerialLocalHost
+        if poltype.externalapi==None:
+            finishedjobs,errorjobs=poltype.CallJobsSeriallyLocalHost(cmd_to_run,skiperrors)
+        else:
+            jobtoinputfilepaths = {cmd: [f'{Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_inp_file}']}
+            jobtooutputfiles = {cmd: [f'{Opt_prep.PySCF_out_file}']}
+            jobtoabsolutebinpath = {cmd: poltype.which('pyscf')}
+            scratchdir = poltype.scrtmpdirpsi4
+            jobtologlistfilepathprefix = f'{Opt_prep.init_data["topdir"]}/optimization_jobtolog_{poltype.molecprefix}'
+       
+            call.CallExternalAPI(poltype,jobtoinputfilepaths,jobtooutputfiles,jobtoabsolutebinpath,scratchdir,jobtologlistfilepathprefix)
+            finishedjobs,errorjobs=poltype.WaitForTermination(f'{Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}',False) 
+
+        # Check normal termination of PySCF job
+        term,error=poltype.CheckNormalTermination(f'{Opt_prep.PySCF_out_file}',None,skiperrors)
+
+        if error and term==False and skiperrors==False:
+            # This poltype method does not exists...
+            poltype.RaiseOutputFileError(f'{Opt_prep.PySCF_out_file}')
+
+        # Instantiate the post process PySCF class
+        Opt_post = PySCF_post_run(poltype,os.getcwd(),f'{Opt_prep.PySCF_out_file}',mol)
+
+        # Read the output file from PySCF
+        Opt_post.read_out()
+
+        # Use OpenBabel to build a molecule out of the optimized structure
+        optmol =  load_structfile(poltype,Opt_post.final_opt_xyz)
+        optmol=rebuild_bonds(poltype,optmol,mol)
+    
+        # Here we redefine logoptfname and poltype.logoptfname
+        # to make sure that hereafter, the pyscf output is used!
+        logoptfname = Opt_prep.PySCF_out_file
+        poltype.logoptfname = Opt_prep.PySCF_out_file 
+
+    if Soft == 'Psi4':
         gen_optcomfile(poltype,comoptfname,poltype.numproc,poltype.maxmem,poltype.maxdisk,chkoptfname,mol,modred,torsionrestraints)
         term,error=poltype.CheckNormalTermination(logoptfname,errormessages=None,skiperrors=True)
         modred=False
 
         inputname,outputname=CreatePsi4OPTInputFile(poltype,comoptfname,comoptfname,mol,modred,bondanglerestraints,skipscferror,charge,loose,torsionrestraints)
         if term==False or overridecheckterm==True:
             poltype.WriteToLog("QM Geometry Optimization")
+            poltype.WriteToLog(f'The software {Soft} will be used to do the QM optimization')
             cmdstr=' '.join(['psi4 ',poltype.psi4_args,inputname,logoptfname])
             jobtooutputlog={cmdstr:os.getcwd()+r'/'+logoptfname}
             jobtolog={cmdstr:os.getcwd()+r'/'+poltype.logfname}
@@ -767,10 +867,10 @@ def GeometryOptimization(poltype,mol,totcharge,suffix='1',loose=False,checkbonds
             jobtoinputfilepaths={cmdstr:[inputfilepath]}
             jobtooutputfiles={cmdstr:[logoptfname]}
             jobtoabsolutebinpath={cmdstr:poltype.which('psi4')}
+    
             if poltype.checkinputonly==True:
                 sys.exit()
 
-
             if os.path.isfile(logoptfname):
                 os.remove(logoptfname)
             if poltype.externalapi==None:
@@ -787,7 +887,10 @@ def GeometryOptimization(poltype,mol,totcharge,suffix='1',loose=False,checkbonds
         optmol =  load_structfile(poltype,logoptfname.replace('.log','.xyz'))
         optmol=rebuild_bonds(poltype,optmol,mol)
     optmol.SetTotalCharge(totcharge)
-    GrabFinalXYZStructure(poltype,logoptfname,logoptfname.replace('.log','.xyz'),mol)
+
+    if Soft != 'PySCF':
+        GrabFinalXYZStructure(poltype,logoptfname,logoptfname.replace('.log','.xyz'),mol)
+
     return optmol,error,torsionrestraints
 
 
@@ -801,6 +904,7 @@ def load_structfile(poltype,structfname):
     Referenced By: run_gaussian, tor_opt_sp, compute_qm_tor_energy, compute_mm_tor_energy 
     Description: -
     """
+
     strctext = os.path.splitext(structfname)[1]
     tmpconv = openbabel.OBConversion()
     if strctext in '.fchk':

PoltypeModules/poltype.py

@@ -456,6 +456,9 @@ class PolarizableTyper():
         espbasisset:str="aug-cc-pVTZ"
         torspbasisset:str="6-311+G*"
         optmethod:str='MP2'
+        pyscf_opt_meth:str = 'wb97x_d3'
+        pyscf_sol_imp:str = 'IEF-PCM' # C-PCM, SS(V)PE, COSMO
+        pyscf_sol_eps:float = 78.3553 # Water
         toroptmethod:str='xtb'
         torspmethod:str='wB97X-D'
         dmamethod:str='MP2'
@@ -1189,6 +1192,12 @@ def __post_init__(self):
                                     warnings.warn(f"Could not locate prmmod file '{fpath1}'")
                         elif newline.startswith("optmethod"):
                             self.optmethod = a
+                        elif newline.startswith("pyscf_opt_met"):
+                            self.pyscf_opt_met = a
+                        elif newline.startswith("pyscf_sol_imp"):
+                            self.pyscf_sol_imp = a
+                        elif newline.startswith("pyscf_sol_eps"):
+                            self.pyscf_sol_eps = a
                         elif newline.startswith("espmethod"):
                             self.espmethod = a
                         elif "torspmethod" in newline:
@@ -3373,6 +3382,7 @@ def CheckNormalTermination(self,logfname,errormessages=None,skiperrors=False):
             5. If error occured, write out line containing error to poltype log file.
             6. xtb only outputs the same filename for optimization jobs, so be sure to copy to desired filename after job finishes. 
             """
+            
             error=False
             term=False
             lastupdatetofilename={}
@@ -3412,7 +3422,7 @@ def CheckNormalTermination(self,logfname,errormessages=None,skiperrors=False):
                     else:
                         if "Final optimized geometry" in line or "Electrostatic potential computed" in line or 'Psi4 exiting successfully' in line or "LBFGS  --  Normal Termination due to SmallGrad" in line or "Normal termination" in line or 'Normal Termination' in line or 'Total Potential Energy' in line or 'Psi4 stopped on' in line or 'finished run' in line or ('Converged! =D' in line):
                             term=True
-                        if ('Tinker is Unable to Continue' in line or 'error' in line or ' Error ' in line or ' ERROR ' in line or 'impossible' in line or 'software termination' in line or 'segmentation violation, address not mapped to object' in line or 'galloc:  could not allocate memory' in line or 'Erroneous write.' in line or 'Optimization failed to converge!' in line) and 'DIIS' not in line and 'mpi' not in line and 'RMS Error' not in line:
+                        if ('Tinker is Unable to Continue' in line or 'error' in line or ' Error ' in line or ' ERROR ' in line or 'impossible' in line or 'software termination' in line or 'segmentation violation, address not mapped to object' in line or 'galloc:  could not allocate memory' in line or 'Erroneous write.' in line or 'Optimization failed to converge!' in line or 'Geometry optimization is not converged' in line or 'Error' in line) and 'DIIS' not in line and 'mpi' not in line and 'RMS Error' not in line:
                             error=True
                             errorline=line
                         if 'segmentation violation' in line and 'address not mapped to object' not in line or 'Waiting' in line or ('OptimizationConvergenceError' in line and 'except' in line) or "Error on total polarization charges" in line or 'Erroneous write' in line:
@@ -3514,7 +3524,13 @@ def call_subsystem(self,cmdstrs,wait=False,skiperrors=False):
                     # STEP 1
                     self.WriteToLog("Submitting: " + cmdstr+' '+'path'+' = '+os.getcwd())
                     # STEP 2
-                    p = subprocess.Popen(cmdstr,shell=True,stdout=self.logfh, stderr=self.logfh)
+                    if 'pyscf' in cmdstr:
+                        out = cmdstr.split()[-1]
+                        out_pyscf = open(out, 'w')
+                        p = subprocess.Popen(cmdstr,shell=True,stdout=out_pyscf, stderr=out_pyscf)
+                    else:
+                        p = subprocess.Popen(cmdstr,shell=True,stdout=self.logfh, stderr=self.logfh)
+
                     procs.append(p)
                     self.cmdtopid[cmdstr]=p
 
@@ -4737,6 +4753,9 @@ def GenerateParameters(self):
             self.localframe1 = [ 0 ] * mol.NumAtoms()
             self.localframe2 = [ 0 ] * mol.NumAtoms()
             self.WriteToLog("Atom Type Classification")
+            self.WriteToLog('Poltype image build time: ')
+            self.WriteToLog('Tinker commit: c9698d2101c5f66ce1d413f4aa2d5f62e4c22df2')
+            self.WriteToLog('Poltype commit: 51b235f8f0af0975c9afae789f671d698b44d7e5')
             self.idxtosymclass,self.symmetryclass=symm.gen_canonicallabels(self,mol,None,self.usesymtypes,True)
             # STEP 15
             torgen.FindPartialDoubleBonds(self,m,mol) # need to find something to hardcode transfer for partial double amide/acid, currently will derive torsion parameters if doesnt find "good" match in torsion database
@@ -5155,6 +5174,9 @@ def GenerateParameters(self):
                 shutil.copy(self.xyzfname,self.xyzoutfile)
             shutil.copy(self.xyzoutfile,self.tmpxyzfile)
             shutil.copy(self.key7fname,self.tmpkeyfile)
+
+            shutil.copy(self.key7fname,'TEST_tim.key')
+            
             # STEP 51
             if self.writeoutpolarize and self.writeoutmultipole==True:
                 opt.StructureMinimization(self,torsionrestraints)

PoltypeModules/pyscf_for_frag.py

@@ -0,0 +1,106 @@
+# Import both PySCF setup class
+from pyscf_setup import PySCF_init_setup
+from pyscf_setup import PySCF_post_run
+import os
+
+def setup_frag_opt(poltype,newtorset,toranglelist,\
+                    pre,inputmol,inputname,inputstruc):
+    
+    """
+    This function setup the PySCF file to do
+    a geometry optimization.
+
+    Inputs:
+        -   poltype: poltype class object (class obj)
+        -   newtorset: list of torsion to restrain (list of list)
+        -   toranglelist: list of torsion angle to restrain (list of float)
+        -   pre: prefix name (string)
+        -   inputmol: molecule object (babel or rdkit mol obj)
+        -   inputname: name of the input file (string)
+        -   inputstruc: name if the file with input geometry (string)
+
+    Outputs:
+        -   Opt_prep.PySCF_inp_file: Name of PySCF input file
+        -   Opt_prep.PySCF_out_file: Name of PySCF output file
+        -   cmd: command to run        
+
+    """
+
+    # Define the charge of the molecule 
+    charge = inputmol.GetTotalCharge()
+
+    # Instantiate the PySCF initial setup object
+    Opt_prep = PySCF_init_setup(poltype,os.getcwd(),pre,inputmol,\
+                               False,charge)
+
+    # Define the input geometry file
+    Opt_prep.PySCF_inp_xyz = inputstruc
+
+    # Write the geometric constrain file
+    Opt_prep.write_geometric_tor_const(pre,newtorset)
+
+    # Write the PySCF input file: is_opt = True, is_frag = True
+    Opt_prep.write_PySCF_input(True,True)
+
+    # Define the command to run
+    cmd = f'python {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_inp_file} &> {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}'
+
+    return Opt_prep.PySCF_inp_file,Opt_prep.PySCF_out_file,cmd
+
+
+def read_frag_opt(poltype,outputlog,mol):
+   
+    """
+    This function creates the post run PySCf object and
+    reads the optimized geometry
+    
+    Inputs:
+        -   poltype: poltype class object (class obj)
+        -   outputlog: Name of the output file (string)
+        -   mol: molecule object (babel or rdkit mol obj)
+
+    """
+
+    # Instantiate the PySCF post run class 
+    Opt_post = PySCF_post_run(poltype,os.getcwd(),outputlog,mol)
+
+    # Read and create the optimized geometry file
+    Opt_post.read_out()
+
+
+def setup_frag_sp(poltype,inputmol,inputstruc,pre):
+
+    """
+    This function set up the fragment SP calculation
+
+    Inputs:
+        -   poltype: poltype class object (class obj)
+        -   inputmol: molecule object (babel or rdkit mol obj)
+        -   inputstruc: name if the file with input geometry (string)
+        -   pre: prefix name (string)
+
+    Outputs:
+        -   Opt_prep.PySCF_inp_file: Name of PySCF input file
+        -   Opt_prep.PySCF_out_file: Name of PySCF output file
+        -   cmd: command to run        
+
+    """
+
+    # Define the charge of the molecule
+    charge = inputmol.GetTotalCharge()
+
+    # Instantiate the PySCF initial setup object 
+    Opt_prep = PySCF_init_setup(poltype,os.getcwd(),pre,inputmol,\
+                               False,charge)
+
+    # Define the input geometry file
+    Opt_prep.PySCF_inp_xyz = inputstruc
+
+    # Write the PySCF input file: is_opt = False, is_frag = True
+    Opt_prep.write_PySCF_input(False,True)
+
+    # Define the command to run
+    cmd = f'python {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_inp_file} &> {Opt_prep.init_data["topdir"]}/{Opt_prep.PySCF_out_file}'
+
+    return Opt_prep.PySCF_inp_file,Opt_prep.PySCF_out_file,cmd
+