Frequency integration grid scaling factor as input to RPA calculations (pyscf#1670)

* scaling factor of frequency integration grid as input parameter in RPA calculations

* default values for nw and x0 in kernel

Author: EgorTrushin

pyscf/gw/rpa.py

@@ -39,14 +39,14 @@
 # core routines, kernel, rpa_ecorr, rho_response
 # ****************************************************************************
 
-def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=None, verbose=logger.NOTE):
+def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=40, x0=0.5, verbose=logger.NOTE):
     """
     RPA correlation and total energy
 
     Args:
         Lpq : density fitting 3-center integral in MO basis.
         nw : number of frequency point on imaginary axis.
-        vhf_df : using density fitting integral to compute HF exchange.
+        x0: scaling factor for frequency grid.
 
     Returns:
         e_tot : RPA total energy
@@ -63,7 +63,7 @@ def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=None, verbose=logger.NOTE):
         Lpq = rpa.ao2mo(mo_coeff)
 
     # Grids for integration on imaginary axis
-    freqs, wts = _get_scaled_legendre_roots(nw)
+    freqs, wts = _get_scaled_legendre_roots(nw, x0)
 
     # Compute HF exchange energy (EXX)
     dm = mf.make_rdm1()
@@ -222,13 +222,14 @@ def nmo(self, n):
     get_nmo = get_nmo
     get_frozen_mask = get_frozen_mask
 
-    def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40):
+    def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40, x0=0.5):
         """
         Args:
             mo_energy : 1D array (nmo), mean-field mo energy
             mo_coeff : 2D array (nmo, nmo), mean-field mo coefficient
             Lpq : 3D array (naux, nmo, nmo), 3-index ERI
             nw: interger, grid number
+            x0: real, scaling factor for frequency grid
 
         Returns:
             self.e_tot : RPA total eenrgy
@@ -243,7 +244,7 @@ def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40):
         cput0 = (logger.process_clock(), logger.perf_counter())
         self.dump_flags()
         self.e_tot, self.e_hf, self.e_corr = \
-                        kernel(self, mo_energy, mo_coeff, Lpq=Lpq, nw=nw, verbose=self.verbose)
+                        kernel(self, mo_energy, mo_coeff, Lpq=Lpq, nw=nw, x0=x0, verbose=self.verbose)
 
         logger.timer(self, 'RPA', *cput0)
         return self.e_corr

pyscf/gw/urpa.py

@@ -40,14 +40,14 @@
 # core routines, kernel, rpa_ecorr, rho_response
 # ****************************************************************************
 
-def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=None, verbose=logger.NOTE):
+def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=40, x0=0.5, verbose=logger.NOTE):
     """
     RPA correlation and total energy
 
     Args:
         Lpq : density fitting 3-center integral in MO basis.
         nw : number of frequency point on imaginary axis.
-        vhf_df : using density fitting integral to compute HF exchange.
+        x0: scaling factor for frequency grid.
 
     Returns:
         e_tot : RPA total energy
@@ -64,7 +64,7 @@ def kernel(rpa, mo_energy, mo_coeff, Lpq=None, nw=None, verbose=logger.NOTE):
         Lpq = rpa.ao2mo(mo_coeff)
 
     # Grids for integration on imaginary axis
-    freqs, wts = _get_scaled_legendre_roots(nw)
+    freqs, wts = _get_scaled_legendre_roots(nw, x0)
 
     # Compute HF exchange energy (EXX)
     dm = mf.make_rdm1()
@@ -153,13 +153,14 @@ def dump_flags(self):
     get_nmo = get_nmo
     get_frozen_mask = get_frozen_mask
 
-    def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40):
+    def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40, x0=0.5):
         """
         Args:
             mo_energy : 2D array (2, nmo), mean-field mo energy
             mo_coeff : 3D array (2, nmo, nmo), mean-field mo coefficient
             Lpq : 4D array (2, naux, nmo, nmo), 3-index ERI
             nw: interger, grid number
+            x0: real, scaling factor for frequency grid
 
         Returns:
             self.e_tot : RPA total eenrgy
@@ -174,7 +175,7 @@ def kernel(self, mo_energy=None, mo_coeff=None, Lpq=None, nw=40):
         cput0 = (logger.process_clock(), logger.perf_counter())
         self.dump_flags()
         self.e_tot, self.e_hf, self.e_corr = \
-                        kernel(self, mo_energy, mo_coeff, Lpq=Lpq, nw=nw, verbose=self.verbose)
+                        kernel(self, mo_energy, mo_coeff, Lpq=Lpq, nw=nw, x0=x0, verbose=self.verbose)
 
         logger.timer(self, 'RPA', *cput0)
         return self.e_corr