THC TDA

.gitignore

@@ -127,3 +127,6 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# pycharm
+.idea/

examples/07-thc_tda.py

@@ -0,0 +1,42 @@
+"""Example of using THC integrals instead of using Cholesky decomposition.
+In this example, the THC integrals are imported from thc_eri_8.h5.
+"""
+
+
+from os.path import abspath, dirname, join
+
+import numpy as np
+from pyscf.pbc import df, dft, gto
+
+from momentGW.gw import GW
+
+cell = gto.Cell()
+cell.atom = """He 0 0 0; He 1 1 1"""
+cell.a = np.eye(3) * 3
+cell.basis = "6-31g"
+cell.verbose = 3
+cell.max_memory = 1e10
+cell.precision = 1e-6
+cell.build()
+
+kpts = cell.make_kpts([1, 1, 1])
+
+# To compare to the Cholesky decomposition code, we need to use
+# Gaussian density fitting and not FFT. The THC integrals are
+# based of the FFT representation, and so there will be some
+# additional error between these calculations.
+
+mf = dft.RKS(cell, xc="pbe")
+mf = mf.density_fit()
+mf.exxdiv = None
+mf.kernel()
+
+print("THC-TDA:")
+gw = GW(mf)
+gw.thc_opts = dict(file_path=abspath(join(dirname(__file__), "..", "examples/thc.h5")))
+gw.polarizability = "thc-dtda"
+gw.kernel(nmom_max=7)
+
+print("TDA:")
+gw = GW(mf)
+gw.kernel(nmom_max=7)

momentGW/base.py

@@ -23,7 +23,7 @@ class BaseGW(lib.StreamObject):
         Default value is `False`.
     polarizability : str, optional
         Type of polarizability to use, can be one of `("drpa",
-        "drpa-exact", "dtda").  Default value is `"drpa"`.
+        "drpa-exact", "dtda", "thc-dtda").  Default value is `"drpa"`.
     npoints : int, optional
         Number of numerical integration points.  Default value is `48`.
     optimise_chempot : bool, optional
@@ -46,6 +46,9 @@ class BaseGW(lib.StreamObject):
         self-consistent scheme.  Default value is `"ia"`.
     compression_tol : float, optional
         Tolerance for the compression.  Default value is `1e-10`.
+    thc_opts : dict, optional
+        Dictionary of options to be used for THC calculations. Current
+        implementation requires a filepath to import the THC integrals.
     {extra_parameters}
     """
 
@@ -66,6 +69,9 @@ class BaseGW(lib.StreamObject):
     )
     compression = "ia"
     compression_tol = 1e-10
+    thc_opts = dict(
+        file_path=None,
+    )
 
     _opts = [
         "diagonal_se",
@@ -76,6 +82,7 @@ class BaseGW(lib.StreamObject):
         "fock_opts",
         "compression",
         "compression_tol",
+        "thc_opts",
     ]
 
     def __init__(self, mf, **kwargs):

momentGW/gw.py

@@ -15,7 +15,7 @@
 from pyscf.ao2mo import _ao2mo
 from pyscf.lib import logger
 
-from momentGW import energy, util
+from momentGW import energy, thc, util
 from momentGW.base import BaseGW
 from momentGW.fock import fock_loop
 from momentGW.ints import Integrals
@@ -182,19 +182,32 @@ def build_se_moments(self, nmom_max, integrals, **kwargs):
             tda = TDA(self, nmom_max, integrals, **kwargs)
             return tda.kernel()
 
+        elif self.polarizability == "thc-dtda":
+            tda = thc.TDA(self, nmom_max, integrals, **kwargs)
+            return tda.kernel()
+
         else:
             raise NotImplementedError
 
     def ao2mo(self, transform=True):
         """Get the integrals."""
 
-        integrals = Integrals(
+        if self.polarizability.startswith("thc"):
+            cls = thc.Integrals
+            kwargs = self.thc_opts
+        else:
+            cls = Integrals
+            kwargs = dict(
+                compression=self.compression,
+                compression_tol=self.compression_tol,
+                store_full=self.fock_loop,
+            )
+
+        integrals = cls(
             self.with_df,
             self.mo_coeff,
             self.mo_occ,
-            compression=self.compression,
-            compression_tol=self.compression_tol,
-            store_full=self.has_fock_loop,
+            **kwargs,
         )
         if transform:
             integrals.transform()