Merge pull request #6 from BoothGroup/ccsd_lmatvec

Support imaginary frequency grids

dyson/expressions/expression.py

@@ -152,7 +152,7 @@ def build_gf_moments(self, nmom, store_vectors=True, left=False):
                     u = apply_hamiltonian(u)
 
         if left:
-            t = t.transpose(1, 2).conj()
+            t = t.transpose(0, 2, 1).conj()
 
         return t
 

dyson/lehmann.py

@@ -447,21 +447,25 @@ def as_perturbed_mo_energy(self):
 
         return mo_energy
 
-    def on_grid(self, grid, eta=1e-1, ordering="time-ordered"):
+    def on_grid(self, grid, eta=1e-1, ordering="time-ordered", axis="real"):
         """
-        Return the Lehmann representation realised on a real frequency
+        Return the Lehmann representation realised on a frequency
         grid.
 
         Parameters
         ----------
         grid : numpy.ndarray
-            Array of real frequency points.
+            Array of frequency points.
         eta : float, optional
             Broadening parameter.  Default value is `1e-1`.
+            Only relevant for real axis.
         ordering : str, optional
             Time ordering.  Can be one of `{"time-ordered",
             "advanced", "retarded"}`.  Default value is
             `"time-ordered"`.
+        axis : str, optional
+            Frequency axis.  Can be one of `{"real", "imag"}`.  Default
+            value is `"real"`.
 
         Returns
         -------
@@ -480,7 +484,12 @@ def on_grid(self, grid, eta=1e-1, ordering="time-ordered"):
 
         couplings_l, couplings_r = self._unpack_couplings()
 
-        denom = 1.0 / lib.direct_sum("w+k-k->wk", grid, signs * 1.0j * eta, self.energies)
+        if axis == "real":
+            denom = 1.0 / lib.direct_sum("w+k-k->wk", grid, signs * 1.0j * eta, self.energies)
+        elif axis == "imag":
+            denom = 1.0 / lib.direct_sum("w-k->wk", 1j * grid, self.energies)
+        else:
+            raise ValueError("axis = {}".format(axis))
         f = lib.einsum("pk,qk,wk->wpq", couplings_l, couplings_r.conj(), denom)
 
         return f

dyson/util/moments.py

@@ -129,13 +129,15 @@ def build_block_tridiagonal(on_diagonal, off_diagonal_upper, off_diagonal_lower=
     m = np.block(
         [
             [
-                on_diagonal[i]
-                if i == j
-                else off_diagonal_upper[j]
-                if j == i - 1
-                else off_diagonal_lower[i]
-                if i == j - 1
-                else zero
+                (
+                    on_diagonal[i]
+                    if i == j
+                    else (
+                        off_diagonal_upper[j]
+                        if j == i - 1
+                        else off_diagonal_lower[i] if i == j - 1 else zero
+                    )
+                )
                 for j in range(len(on_diagonal))
             ]
             for i in range(len(on_diagonal))

pyproject.toml

@@ -44,7 +44,7 @@ build-backend = "setuptools.build_meta"
 
 [project.optional-dependencies]
 dev = [
-    "black==23.1.0",
+    "black>=22.6.0",
     "isort>=5.10.1",
     "coverage[toml]",
     "pytest",