Fixes for momentGW

dyson/expressions/gw.py

@@ -33,7 +33,7 @@ def __init__(self, *args, **kwargs):
             raise ImportError("momentGW is required for GW expressions.")
 
     def get_static_part(self):
-        static = self._gw.build_se_static()
+        static = self._gw.build_se_static(self._gw.ao2mo())
 
         return static
 
@@ -103,7 +103,9 @@ def diagonal(self, static=None):
         ija = slice(self.nmo, self.nmo + self.nocc * self.nocc * self.nvir)
         iab = slice(self.nmo + self.nocc * self.nocc * self.nvir, None)
 
-        Lpq, Lia = self._gw.ao2mo(self.mo_coeff)
+        integrals = self._gw.ao2mo()
+        Lpq = integrals.Lpx
+        Lia = integrals.Lia
         Lia = Lpq[:, i, a]
         Lai = Lpq[:, a, i]
         Lij = Lpq[:, i, i]
@@ -142,8 +144,8 @@ def get_wavefunction(self, orb):
         return r
 
     def build_se_moments(self, nmom):
-        Lpq, Lia = self._gw.ao2mo(self.mo_coeff)
-        moments = self._gw.build_se_moments(nmom, Lpq, Lia)
+        integrals = self._gw.ao2mo()
+        moments = self._gw.build_se_moments(nmom, integrals)
 
         return moments
 

dyson/lehmann.py

@@ -3,6 +3,7 @@
 """
 
 import numpy as np
+from pyscf import lib
 
 
 class Lehmann:
@@ -350,10 +351,8 @@ def weights(self, occupancy=1):
             Weights of the states.
         """
 
-        if self.hermitian:
-            wt = np.sum(self.couplings**2, axis=0) * occupancy
-        else:
-            wt = np.sum(self.couplings[0] * self.couplings[1].conj(), axis=0) * occupancy
+        couplings_l, couplings_r = self._unpack_couplings()
+        wt = np.sum(couplings_l * couplings_r.conj(), axis=0) * occupancy
 
         return wt
 
@@ -391,7 +390,7 @@ def as_orbitals(self, occupancy=1, mo_coeff=None):
             orb_coeff = self.couplings
 
         orb_occ = np.zeros_like(orb_energy)
-        orb_occ[orb_energy < self.chempot] = self.occupied().weights(occupancy=occupancy)
+        orb_occ[orb_energy < self.chempot] = np.abs(self.occupied().weights(occupancy=occupancy))
 
         return orb_energy, orb_coeff, orb_occ
 
@@ -448,6 +447,44 @@ def as_perturbed_mo_energy(self):
 
         return mo_energy
 
+    def on_grid(self, grid, eta=1e-1, ordering="time-ordered"):
+        """
+        Return the Lehmann representation realised on a real frequency
+        grid.
+
+        Parameters
+        ----------
+        grid : numpy.ndarray
+            Array of real frequency points.
+        eta : float, optional
+            Broadening parameter.  Default value is `1e-1`.
+        ordering : str, optional
+            Time ordering.  Can be one of `{"time-ordered",
+            "advanced", "retarded"}`.  Default value is
+            `"time-ordered"`.
+
+        Returns
+        -------
+        f : numpy.ndarray
+            Lehmann representation realised at each frequency point.
+        """
+
+        if ordering == "time-ordered":
+            signs = np.sign(self.energies - self.chempot)
+        elif ordering == "advanced":
+            signs = -np.ones_like(self.energies)
+        elif ordering == "retarded":
+            signs = np.ones_like(self.energies)
+        else:
+            raise ValueError("ordering = {}".format(ordering))
+
+        couplings_l, couplings_r = self._unpack_couplings()
+
+        denom = 1.0 / lib.direct_sum("w+k-k->wk", grid, signs * 1.0j * eta, self.energies)
+        f = lib.einsum("pk,qk,wk->wpq", couplings_l, couplings_r.conj(), denom)
+
+        return f
+
     @property
     def hermitian(self):
         """Boolean flag for the Hermiticity."""

dyson/solvers/cpgf.py

@@ -128,6 +128,9 @@ def _kernel(self, iteration=None, trace=True):
                 integral = scipy.integrate.simps(as_trace(gf.imag), self.grid)
                 self.log.info("%4d %16.8g", niter, integral)
 
+        # Not sure why we need to do this...
+        gf = -gf.conj()
+
         self.log.info("-" * 21)
 
         return filter_type(gf)

tests/expressions/test_gw.py

@@ -51,8 +51,8 @@ def test_moment_gw(self):
 
         import momentGW
         gw_ref = momentGW.GW(self.mf)
-        _, gf_ref, se_ref = gw_ref.kernel(9)
-        gf_ref.remove_uncoupled(tol=0.1)
+        _, gf_ref, se_ref, _ = gw_ref.kernel(9)
+        gf_ref = gf_ref.physical(weight=0.1)
 
         np.testing.assert_allclose(gf_ref.moment(0), gf.moment(0), rtol=1e10, atol=1e-10)
         np.testing.assert_allclose(gf_ref.moment(1), gf.moment(1), rtol=1e10, atol=1e-10)