From 30f9c8bee93a050cdd0bbd39aa248e0f3bd1551b Mon Sep 17 00:00:00 2001
From: Ollie Backhouse <olbackhouse@gmail.com>
Date: Mon, 7 Aug 2023 17:54:05 +0100
Subject: [PATCH] Fock changes

---
 momentGW/pbc/fock.py |  5 +++++
 momentGW/pbc/ints.py | 11 ++++++-----
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/momentGW/pbc/fock.py b/momentGW/pbc/fock.py
index ed753380..fffad886 100644
--- a/momentGW/pbc/fock.py
+++ b/momentGW/pbc/fock.py
@@ -68,6 +68,11 @@ def fock_loop(
             fock = integrals.get_fock(rdm1, h1e)
             fock = diis.update(fock, xerr=None)
 
+            rdm1_ao = lib.einsum("kij,kpi,kqj->kpq", rdm1, gw.mo_coeff, np.conj(gw.mo_coeff))
+            fock_ao = gw._scf.get_fock(dm=rdm1_ao)
+            fock_test = lib.einsum("kpq,kpi,kqj->kij", fock_ao, np.conj(gw.mo_coeff), gw.mo_coeff)
+            assert np.allclose(integrals.get_fock(rdm1, h1e), fock_test)
+
             nerr = nerr[np.argmax(np.abs(nerr))]
             if niter2 > 1:
                 derr = np.max(np.absolute(rdm1 - rdm1_prev))
diff --git a/momentGW/pbc/ints.py b/momentGW/pbc/ints.py
index 7d2e5d5d..0fee1f03 100644
--- a/momentGW/pbc/ints.py
+++ b/momentGW/pbc/ints.py
@@ -118,7 +118,8 @@ def transform(self, do_Lpq=None, do_Lpx=True, do_Lia=True):
         if self._rot is None:
             self._rot = self.get_compression_metric()
         rot = self._rot
-        dtype = complex
+        if rot is None:
+            rot = np.eye(self.naux_full)
 
         do_Lpq = self.store_full if do_Lpq is None else do_Lpq
         if not any([do_Lpq, do_Lpx, do_Lia]):
@@ -140,10 +141,10 @@ def transform(self, do_Lpq=None, do_Lpx=True, do_Lia=True):
             o0, o1 = 0, self.nmo
             p0, p1 = 0, self.nmo_g[ki]
             q0, q1 = 0, self.nocc_w[kj] * self.nvir_w[kj]
-            Lpq_k = np.zeros((self.naux_full, self.nmo, o1 - o0), dtype=dtype) if do_Lpq else None
-            Lpx_k = np.zeros((self.naux, self.nmo, p1 - p0), dtype=dtype) if do_Lpx else None
-            Lia_k = np.zeros((self.naux, q1 - q0), dtype=dtype) if do_Lia else None
-            Lai_k = np.zeros((self.naux, q1 - q0), dtype=dtype) if do_Lia else None
+            Lpq_k = np.zeros((self.naux_full, self.nmo, o1 - o0), dtype=complex) if do_Lpq else None
+            Lpx_k = np.zeros((self.naux, self.nmo, p1 - p0), dtype=complex) if do_Lpx else None
+            Lia_k = np.zeros((self.naux, q1 - q0), dtype=complex) if do_Lia else None
+            Lai_k = np.zeros((self.naux, q1 - q0), dtype=complex) if do_Lia else None
 
             # Build the integrals blockwise
             b1 = 0