From edd6c51eaac4b732c538995bf0a9d6aff773aebc Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Mon, 13 Jan 2020 14:44:18 +0100
Subject: [PATCH 01/13] initial changes

might need to change the indices of the C-matrix?
---
 coupled_cluster/oatdcc.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/coupled_cluster/oatdcc.py b/coupled_cluster/oatdcc.py
index f5b6be5..d4a98f0 100644
--- a/coupled_cluster/oatdcc.py
+++ b/coupled_cluster/oatdcc.py
@@ -152,13 +152,15 @@ def __call__(self, prev_amp, current_time):
         here that expm refers to the matrix exponential which I can not find in
         numpy only in scipy.
         """
-        # rho_pq_inv = self.np.linalg.inv(self.rho_qp)
+        rho_pq_inv = self.np.linalg.inv(self.rho_qp)
 
         # Solve Q-space for C and C_tilde
+
+        """
         C_new = np.dot(C, eta)
         C_tilde_new = -np.dot(eta, C_tilde)
-
         """
+
         C_new = -1j * compute_q_space_ket_equations(
             C,
             C_tilde,
@@ -183,7 +185,6 @@ def __call__(self, prev_amp, current_time):
             self.rho_qspr,
             np=np,
         )
-        """
 
         self.last_timestep = current_time
 
@@ -196,12 +197,13 @@ def __call__(self, prev_amp, current_time):
 def compute_q_space_ket_equations(
     C, C_tilde, eta, h, h_tilde, u, u_tilde, rho_inv_pq, rho_qspr, np
 ):
+    o = self.system.o
     rhs = 1j * np.dot(C, eta)
 
     rhs += np.dot(h, C)
     rhs -= np.dot(C, h_tilde)
 
-    u_quart = np.einsum("rb,gq,ds,abgd->arqs", C_tilde, C, C, u, optimize=True)
+    u_quart = np.einsum("rb,gq,ds,abgd->arqs", C_tilde[:], C, C, u, optimize=True)
     u_quart -= np.tensordot(C, u_tilde, axes=((1), (0)))
 
     temp_ap = np.tensordot(u_quart, rho_qspr, axes=((1, 2, 3), (3, 0, 1)))

From 3a2d2c3bbe2d0fe43c7888b0a5fcf9181e3c0a2c Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Wed, 12 Feb 2020 15:09:58 +0100
Subject: [PATCH 02/13] working ccd imaginary time propagation

---
 coupled_cluster/ccd/itdccd.py | 36 +++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)
 create mode 100644 coupled_cluster/ccd/itdccd.py

diff --git a/coupled_cluster/ccd/itdccd.py b/coupled_cluster/ccd/itdccd.py
new file mode 100644
index 0000000..3317125
--- /dev/null
+++ b/coupled_cluster/ccd/itdccd.py
@@ -0,0 +1,36 @@
+from coupled_cluster import TDCCD
+from coupled_cluster.cc_helper import ( AmplitudeContainer)
+from coupled_cluster.integrators import RungeKutta4
+
+class ITDCCD(TDCCD):
+    """Performs imaginary time propagation for ground state calculations."""
+    def __call__(self, prev_amp, current_time):
+        o, v = self.system.o, self.system.v
+
+        prev_amp = AmplitudeContainer.from_array(self._amplitudes, prev_amp)
+        t_old, l_old = prev_amp
+
+        self.update_hamiltonian(current_time, prev_amp)
+
+        # Remove phase from t-amplitude list
+        t_old = t_old[1:]
+
+        t_new = [
+            - rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            for rhs_t_func in self.rhs_t_amplitudes()
+        ]
+
+        # Compute derivative of phase
+        t_0_new = - self.rhs_t_0_amplitude(
+            self.f, self.u, *t_old, self.o, self.v, np=self.np
+        )
+        t_new = [t_0_new, *t_new]
+
+        l_new = [
+            - rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            for rhs_l_func in self.rhs_l_amplitudes()
+        ]
+
+        self.last_timestep = current_time
+
+        return AmplitudeContainer(t=t_new, l=l_new, np=self.np).asarray()

From aa3380516969193ffe3f7e080e8e6ca5f1cdbdc5 Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Wed, 12 Feb 2020 15:29:55 +0100
Subject: [PATCH 03/13] black, __init__, remove import

---
 coupled_cluster/__init__.py     |  3 ++-
 coupled_cluster/ccd/__init__.py |  3 +++
 coupled_cluster/ccd/itdccd.py   | 11 ++++++-----
 3 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/coupled_cluster/__init__.py b/coupled_cluster/__init__.py
index 1af549e..01e07d0 100644
--- a/coupled_cluster/__init__.py
+++ b/coupled_cluster/__init__.py
@@ -1,2 +1,3 @@
-from .ccd import CCD, CoupledClusterDoubles, OACCD, TDCCD, OATDCCD
+from .ccd import (CCD, CoupledClusterDoubles, OACCD, TDCCD, 
+                  OATDCCD, ITDCCD, OAITDCCD)
 from .ccsd import CCSD, CoupledClusterSinglesDoubles, TDCCSD
diff --git a/coupled_cluster/ccd/__init__.py b/coupled_cluster/ccd/__init__.py
index 4e0652d..1004588 100644
--- a/coupled_cluster/ccd/__init__.py
+++ b/coupled_cluster/ccd/__init__.py
@@ -1,4 +1,7 @@
 from .ccd import CCD, CoupledClusterDoubles
 from .oaccd import OACCD
 from .tdccd import TDCCD
+from .itdccd import ITDCCD
 from .oatdccd import OATDCCD
+from .oaitdccd import OAITDCCD
+
diff --git a/coupled_cluster/ccd/itdccd.py b/coupled_cluster/ccd/itdccd.py
index 3317125..37adf30 100644
--- a/coupled_cluster/ccd/itdccd.py
+++ b/coupled_cluster/ccd/itdccd.py
@@ -1,9 +1,10 @@
 from coupled_cluster import TDCCD
-from coupled_cluster.cc_helper import ( AmplitudeContainer)
-from coupled_cluster.integrators import RungeKutta4
+from coupled_cluster.cc_helper import AmplitudeContainer
+
 
 class ITDCCD(TDCCD):
     """Performs imaginary time propagation for ground state calculations."""
+
     def __call__(self, prev_amp, current_time):
         o, v = self.system.o, self.system.v
 
@@ -16,18 +17,18 @@ def __call__(self, prev_amp, current_time):
         t_old = t_old[1:]
 
         t_new = [
-            - rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            -rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
             for rhs_t_func in self.rhs_t_amplitudes()
         ]
 
         # Compute derivative of phase
-        t_0_new = - self.rhs_t_0_amplitude(
+        t_0_new = -self.rhs_t_0_amplitude(
             self.f, self.u, *t_old, self.o, self.v, np=self.np
         )
         t_new = [t_0_new, *t_new]
 
         l_new = [
-            - rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            -rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
             for rhs_l_func in self.rhs_l_amplitudes()
         ]
 

From c1863b783f7207e0a13b3397447eb37bad4ffae1 Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Wed, 12 Feb 2020 15:31:02 +0100
Subject: [PATCH 04/13] remove oaitdccd.py, as it is not implemented yet

---
 coupled_cluster/__init__.py     | 2 +-
 coupled_cluster/ccd/__init__.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/coupled_cluster/__init__.py b/coupled_cluster/__init__.py
index 01e07d0..9de1e1d 100644
--- a/coupled_cluster/__init__.py
+++ b/coupled_cluster/__init__.py
@@ -1,3 +1,3 @@
 from .ccd import (CCD, CoupledClusterDoubles, OACCD, TDCCD, 
-                  OATDCCD, ITDCCD, OAITDCCD)
+                  OATDCCD, ITDCCD)
 from .ccsd import CCSD, CoupledClusterSinglesDoubles, TDCCSD
diff --git a/coupled_cluster/ccd/__init__.py b/coupled_cluster/ccd/__init__.py
index 1004588..96c1cfc 100644
--- a/coupled_cluster/ccd/__init__.py
+++ b/coupled_cluster/ccd/__init__.py
@@ -3,5 +3,4 @@
 from .tdccd import TDCCD
 from .itdccd import ITDCCD
 from .oatdccd import OATDCCD
-from .oaitdccd import OAITDCCD
 

From a0cc674c0e8653fa0ab0c74799c152e0e8820215 Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Wed, 12 Feb 2020 15:42:00 +0100
Subject: [PATCH 05/13] oaitdccd with Q=0 working

---
 coupled_cluster/__init__.py     |  2 +-
 coupled_cluster/ccd/__init__.py |  1 +
 coupled_cluster/ccd/oaitdccd.py | 95 +++++++++++++++++++++++++++++++++
 3 files changed, 97 insertions(+), 1 deletion(-)
 create mode 100644 coupled_cluster/ccd/oaitdccd.py

diff --git a/coupled_cluster/__init__.py b/coupled_cluster/__init__.py
index 9de1e1d..01e07d0 100644
--- a/coupled_cluster/__init__.py
+++ b/coupled_cluster/__init__.py
@@ -1,3 +1,3 @@
 from .ccd import (CCD, CoupledClusterDoubles, OACCD, TDCCD, 
-                  OATDCCD, ITDCCD)
+                  OATDCCD, ITDCCD, OAITDCCD)
 from .ccsd import CCSD, CoupledClusterSinglesDoubles, TDCCSD
diff --git a/coupled_cluster/ccd/__init__.py b/coupled_cluster/ccd/__init__.py
index 96c1cfc..1004588 100644
--- a/coupled_cluster/ccd/__init__.py
+++ b/coupled_cluster/ccd/__init__.py
@@ -3,4 +3,5 @@
 from .tdccd import TDCCD
 from .itdccd import ITDCCD
 from .oatdccd import OATDCCD
+from .oaitdccd import OAITDCCD
 
diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
new file mode 100644
index 0000000..9446c71
--- /dev/null
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -0,0 +1,95 @@
+from coupled_cluster import OATDCCD
+from coupled_cluster.cc_helper import OACCVector
+
+class OAITDCCD(OATDCCD):
+    def __call__(self, prev_amp, current_time):
+        np = self.np
+        o, v = self.o, self.v
+
+        prev_amp = OACCVector.from_array(self._amplitudes, prev_amp)
+        t_old, l_old, C, C_tilde = prev_amp
+
+        self.update_hamiltonian(current_time, prev_amp)
+
+        # Remove t_0 phase as this is not used in any of the equations
+        t_old = t_old[1:]
+
+        # OATDCC procedure:
+        # Do amplitude step
+        t_new = [
+            - rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            for rhs_t_func in self.rhs_t_amplitudes()
+        ]
+
+        # Compute derivative of phase
+        t_0_new = - self.rhs_t_0_amplitude(
+            self.f, self.u, *t_old, self.o, self.v, np=self.np
+        )
+
+        t_new = [t_0_new, *t_new]
+
+        l_new = [
+            - rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            for rhs_l_func in self.rhs_l_amplitudes()
+        ]
+
+        # Compute density matrices
+        self.rho_qp = self.one_body_density_matrix(t_old, l_old)
+        self.rho_qspr = self.two_body_density_matrix(t_old, l_old)
+
+        # Solve P-space equations for eta
+        # divide by imaginary number
+        eta = -1j*self.compute_p_space_equations() 
+        # TODO: move 1j out of compute_p_space_equations 
+
+        # Compute the inverse of rho_qp needed in Q-space eqs.
+        """
+        If rho_qp is singular we can regularize it as,
+
+        rho_qp_reg = rho_qp + eps*expm( -(1.0/eps) * rho_qp) Eq [3.14]
+        Multidimensional Quantum Dynamics, Meyer
+
+        with eps = 1e-8 (or some small number). It seems like it is standard in
+        the MCTDHF literature to always work with the regularized rho_qp. Note
+        here that expm refers to the matrix exponential which I can not find in
+        numpy only in scipy.
+        """
+        # rho_pq_inv = self.np.linalg.inv(self.rho_qp)
+
+        # Solve Q-space for C and C_tilde
+        C_new = np.dot(C, eta)
+        C_tilde_new = -np.dot(eta, C_tilde)
+
+        """
+        C_new = -1j * compute_q_space_ket_equations(
+            C,
+            C_tilde,
+            eta,
+            self.h_orig,
+            self.h,
+            self.u_orig,
+            self.u,
+            rho_pq_inv,
+            self.rho_qspr,
+            np=np,
+        )
+        C_tilde_new = 1j * compute_q_space_bra_equations(
+            C,
+            C_tilde,
+            eta,
+            self.h_orig,
+            self.h,
+            self.u_orig,
+            self.u,
+            rho_pq_inv,
+            self.rho_qspr,
+            np=np,
+        )
+        """
+
+        self.last_timestep = current_time
+
+        # Return amplitudes and C and C_tilde
+        return OACCVector(
+            t=t_new, l=l_new, C=C_new, C_tilde=C_tilde_new, np=self.np
+        ).asarray()

From d3edb70b76471f13882b67fd3ced43a194293ddf Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Wed, 12 Feb 2020 15:51:50 +0100
Subject: [PATCH 06/13] q-space imaginary time propagation workinf for Q=0

blows up if not stopped sensibly
---
 coupled_cluster/ccd/oaitdccd.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
index 9446c71..c2b8faa 100644
--- a/coupled_cluster/ccd/oaitdccd.py
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -1,8 +1,10 @@
 from coupled_cluster import OATDCCD
 from coupled_cluster.cc_helper import OACCVector
 
+
 class OAITDCCD(OATDCCD):
     def __call__(self, prev_amp, current_time):
+        print(1)
         np = self.np
         o, v = self.o, self.v
 
@@ -17,19 +19,19 @@ def __call__(self, prev_amp, current_time):
         # OATDCC procedure:
         # Do amplitude step
         t_new = [
-            - rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            -rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
             for rhs_t_func in self.rhs_t_amplitudes()
         ]
 
         # Compute derivative of phase
-        t_0_new = - self.rhs_t_0_amplitude(
+        t_0_new = -self.rhs_t_0_amplitude(
             self.f, self.u, *t_old, self.o, self.v, np=self.np
         )
 
         t_new = [t_0_new, *t_new]
 
         l_new = [
-            - rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            -rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
             for rhs_l_func in self.rhs_l_amplitudes()
         ]
 
@@ -39,8 +41,8 @@ def __call__(self, prev_amp, current_time):
 
         # Solve P-space equations for eta
         # divide by imaginary number
-        eta = -1j*self.compute_p_space_equations() 
-        # TODO: move 1j out of compute_p_space_equations 
+        eta = -1j * self.compute_p_space_equations()
+        # TODO: move 1j out of compute_p_space_equations
 
         # Compute the inverse of rho_qp needed in Q-space eqs.
         """
@@ -57,11 +59,10 @@ def __call__(self, prev_amp, current_time):
         # rho_pq_inv = self.np.linalg.inv(self.rho_qp)
 
         # Solve Q-space for C and C_tilde
-        C_new = np.dot(C, eta)
-        C_tilde_new = -np.dot(eta, C_tilde)
+        # C_new = np.dot(C, eta)
+        # C_tilde_new = -np.dot(eta, C_tilde)
 
-        """
-        C_new = -1j * compute_q_space_ket_equations(
+        C_new = -compute_q_space_ket_equations(
             C,
             C_tilde,
             eta,
@@ -73,7 +74,7 @@ def __call__(self, prev_amp, current_time):
             self.rho_qspr,
             np=np,
         )
-        C_tilde_new = 1j * compute_q_space_bra_equations(
+        C_tilde_new = -compute_q_space_bra_equations(
             C,
             C_tilde,
             eta,
@@ -85,7 +86,6 @@ def __call__(self, prev_amp, current_time):
             self.rho_qspr,
             np=np,
         )
-        """
 
         self.last_timestep = current_time
 

From 25ed1769256cd68411d4af5a1cb29ee49e57b6ba Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Thu, 13 Feb 2020 13:44:28 +0100
Subject: [PATCH 07/13] black and add q-space equations

---
 coupled_cluster/ccd/oaitdccd.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
index c2b8faa..95c7837 100644
--- a/coupled_cluster/ccd/oaitdccd.py
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -4,7 +4,6 @@
 
 class OAITDCCD(OATDCCD):
     def __call__(self, prev_amp, current_time):
-        print(1)
         np = self.np
         o, v = self.o, self.v
 
@@ -40,7 +39,7 @@ def __call__(self, prev_amp, current_time):
         self.rho_qspr = self.two_body_density_matrix(t_old, l_old)
 
         # Solve P-space equations for eta
-        # divide by imaginary number
+        # as in regular td, but divided by imaginary number
         eta = -1j * self.compute_p_space_equations()
         # TODO: move 1j out of compute_p_space_equations
 
@@ -59,8 +58,10 @@ def __call__(self, prev_amp, current_time):
         # rho_pq_inv = self.np.linalg.inv(self.rho_qp)
 
         # Solve Q-space for C and C_tilde
-        # C_new = np.dot(C, eta)
-        # C_tilde_new = -np.dot(eta, C_tilde)
+        """
+        C_new = -1j*(np.dot(C, eta))
+        C_tilde_new = 1j* (-np.dot(eta, C_tilde))
+        """
 
         C_new = -compute_q_space_ket_equations(
             C,

From 66d32627898b3334705e18adb7ccdee18d62921f Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Fri, 14 Feb 2020 13:41:22 +0100
Subject: [PATCH 08/13] add residuals

---
 coupled_cluster/cc_helper.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/coupled_cluster/cc_helper.py b/coupled_cluster/cc_helper.py
index 4dea8fe..5c35e96 100644
--- a/coupled_cluster/cc_helper.py
+++ b/coupled_cluster/cc_helper.py
@@ -157,6 +157,12 @@ def from_array(u, arr):
 
         return type(u)(*args, np=np)
 
+    def residuals(self):
+        return [
+            [np.linalg.norm(t) for t in self.t],
+            [np.linalg.norm(l) for l in self.l],
+        ]
+
 
 class OACCVector(AmplitudeContainer):
     """Container for OA amplitudes

From 6a2274524f7cf2048edb813d733ae0d671b5dc13 Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Fri, 14 Feb 2020 14:04:40 +0100
Subject: [PATCH 09/13] recover untracked, stashed files operation

---
 coupled_cluster/ccd/itdccd.py   |  37 ++++++++++++
 coupled_cluster/ccd/oaitdccd.py | 100 ++++++++++++++++++++++++++++++++
 2 files changed, 137 insertions(+)
 create mode 100644 coupled_cluster/ccd/itdccd.py
 create mode 100644 coupled_cluster/ccd/oaitdccd.py

diff --git a/coupled_cluster/ccd/itdccd.py b/coupled_cluster/ccd/itdccd.py
new file mode 100644
index 0000000..d991cc2
--- /dev/null
+++ b/coupled_cluster/ccd/itdccd.py
@@ -0,0 +1,37 @@
+from coupled_cluster.ccd.tdccd import TDCCD
+from coupled_cluster.cc_helper import AmplitudeContainer
+
+
+class ITDCCD(TDCCD):
+    """Performs imaginary time propagation for ground state calculations."""
+
+    def __call__(self, prev_amp, current_time):
+        o, v = self.system.o, self.system.v
+
+        prev_amp = AmplitudeContainer.from_array(self._amplitudes, prev_amp)
+        t_old, l_old = prev_amp
+
+        self.update_hamiltonian(current_time, prev_amp)
+
+        # Remove phase from t-amplitude list
+        t_old = t_old[1:]
+
+        t_new = [
+            -rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            for rhs_t_func in self.rhs_t_amplitudes()
+        ]
+
+        # Compute derivative of phase
+        t_0_new = -self.rhs_t_0_amplitude(
+            self.f, self.u, *t_old, self.o, self.v, np=self.np
+        )
+        t_new = [t_0_new, *t_new]
+
+        l_new = [
+            -rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            for rhs_l_func in self.rhs_l_amplitudes()
+        ]
+
+        self.last_timestep = current_time
+
+        return AmplitudeContainer(t=t_new, l=l_new, np=self.np).asarray()
diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
new file mode 100644
index 0000000..d309add
--- /dev/null
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -0,0 +1,100 @@
+from coupled_cluster.ccd.oatdccd import OATDCCD
+from coupled_cluster.oatdcc import (
+    compute_q_space_ket_equations,
+    compute_q_space_bra_equations,
+)
+from coupled_cluster.cc_helper import OACCVector
+
+
+class OAITDCCD(OATDCCD):
+    def __call__(self, prev_amp, current_time):
+        np = self.np
+        o, v = self.o, self.v
+
+        prev_amp = OACCVector.from_array(self._amplitudes, prev_amp)
+        t_old, l_old, C, C_tilde = prev_amp
+
+        self.update_hamiltonian(current_time, prev_amp)
+
+        # Remove t_0 phase as this is not used in any of the equations
+        t_old = t_old[1:]
+
+        # OATDCC procedure:
+        # Do amplitude step
+        t_new = [
+            -rhs_t_func(self.f, self.u, *t_old, o, v, np=self.np)
+            for rhs_t_func in self.rhs_t_amplitudes()
+        ]
+
+        # Compute derivative of phase
+        t_0_new = -self.rhs_t_0_amplitude(
+            self.f, self.u, *t_old, self.o, self.v, np=self.np
+        )
+
+        t_new = [t_0_new, *t_new]
+
+        l_new = [
+            -rhs_l_func(self.f, self.u, *t_old, *l_old, o, v, np=self.np)
+            for rhs_l_func in self.rhs_l_amplitudes()
+        ]
+
+        # Compute density matrices
+        self.rho_qp = self.one_body_density_matrix(t_old, l_old)
+        self.rho_qspr = self.two_body_density_matrix(t_old, l_old)
+
+        # Solve P-space equations for eta
+        # as in regular td, but divided by imaginary number
+        eta = -1j * self.compute_p_space_equations()
+        # TODO: move 1j out of compute_p_space_equations
+
+        # Compute the inverse of rho_qp needed in Q-space eqs.
+        """
+        If rho_qp is singular we can regularize it as,
+
+        rho_qp_reg = rho_qp + eps*expm( -(1.0/eps) * rho_qp) Eq [3.14]
+        Multidimensional Quantum Dynamics, Meyer
+
+        with eps = 1e-8 (or some small number). It seems like it is standard in
+        the MCTDHF literature to always work with the regularized rho_qp. Note
+        here that expm refers to the matrix exponential which I can not find in
+        numpy only in scipy.
+        """
+        # rho_pq_inv = self.np.linalg.inv(self.rho_qp)
+
+        # Solve Q-space for C and C_tilde
+        C_new = -1j * (np.dot(C, eta))
+        C_tilde_new = 1j * (-np.dot(eta, C_tilde))
+        """
+
+        C_new = -compute_q_space_ket_equations(
+            C,
+            C_tilde,
+            eta,
+            self.h_orig,
+            self.h,
+            self.u_orig,
+            self.u,
+            rho_pq_inv,
+            self.rho_qspr,
+            np=np,
+        )
+        C_tilde_new = -compute_q_space_bra_equations(
+            C,
+            C_tilde,
+            eta,
+            self.h_orig,
+            self.h,
+            self.u_orig,
+            self.u,
+            rho_pq_inv,
+            self.rho_qspr,
+            np=np,
+        )
+        """
+
+        self.last_timestep = current_time
+
+        # Return amplitudes and C and C_tilde
+        return OACCVector(
+            t=t_new, l=l_new, C=C_new, C_tilde=C_tilde_new, np=self.np
+        ).asarray()

From 0357e9193f8a41152dadad646febb9835d840a2e Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Tue, 25 Feb 2020 13:00:40 +0100
Subject: [PATCH 10/13] Relative imports

---
 coupled_cluster/ccd/itdccd.py   | 4 ++--
 coupled_cluster/ccd/oaitdccd.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/coupled_cluster/ccd/itdccd.py b/coupled_cluster/ccd/itdccd.py
index 37adf30..9c5313f 100644
--- a/coupled_cluster/ccd/itdccd.py
+++ b/coupled_cluster/ccd/itdccd.py
@@ -1,5 +1,5 @@
-from coupled_cluster import TDCCD
-from coupled_cluster.cc_helper import AmplitudeContainer
+from . import TDCCD
+from ..cc_helper import AmplitudeContainer
 
 
 class ITDCCD(TDCCD):
diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
index c2b8faa..9b8b83d 100644
--- a/coupled_cluster/ccd/oaitdccd.py
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -1,5 +1,5 @@
-from coupled_cluster import OATDCCD
-from coupled_cluster.cc_helper import OACCVector
+from . import OATDCCD
+from ..cc_helper import OACCVector
 
 
 class OAITDCCD(OATDCCD):

From 861a8ba71b961b0ee99bfac8a9374391aeb7f85c Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Tue, 25 Feb 2020 13:13:23 +0100
Subject: [PATCH 11/13] New version of quantum-system

---
 coupled_cluster/ccd/oaccd.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/coupled_cluster/ccd/oaccd.py b/coupled_cluster/ccd/oaccd.py
index 892d9f7..910017d 100644
--- a/coupled_cluster/ccd/oaccd.py
+++ b/coupled_cluster/ccd/oaccd.py
@@ -153,7 +153,7 @@ def compute_ground_state(
             if self.verbose:
                 print("Changing system basis...")
 
-            self.system.change_basis(c=S, c_tilde=S_inv)
+            self.system.change_basis(C=S, C_tilde=S_inv)
 
         if self.verbose:
             print(

From 8dc6febe8a5994ca310a6a564e06c6d6dbee7987 Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Tue, 25 Feb 2020 13:16:35 +0100
Subject: [PATCH 12/13] mostly black

---
 coupled_cluster/ccd/__init__.py |  1 -
 coupled_cluster/ccd/oaitdccd.py | 14 +++++++++++++-
 coupled_cluster/tdcc.py         |  1 +
 3 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/coupled_cluster/ccd/__init__.py b/coupled_cluster/ccd/__init__.py
index 1004588..9077af1 100644
--- a/coupled_cluster/ccd/__init__.py
+++ b/coupled_cluster/ccd/__init__.py
@@ -4,4 +4,3 @@
 from .itdccd import ITDCCD
 from .oatdccd import OATDCCD
 from .oaitdccd import OAITDCCD
-
diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
index d309add..eaa9c05 100644
--- a/coupled_cluster/ccd/oaitdccd.py
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -7,6 +7,18 @@
 
 
 class OAITDCCD(OATDCCD):
+    def compute_ground_state(self, *args, **kwargs):
+        if "change_system_basis" not in kwargs:
+            kwargs["change_system_basis"] = True
+
+        # self.cc.compute_ground_state(*args, **kwargs)
+        self.h_orig = self.system.h
+        self.u_orig = self.system.u
+
+        self.h = self.system.h
+        self.u = self.system.u
+        self.f = self.system.construct_fock_matrix(self.h, self.u)
+
     def __call__(self, prev_amp, current_time):
         np = self.np
         o, v = self.o, self.v
@@ -64,8 +76,8 @@ def __call__(self, prev_amp, current_time):
         # Solve Q-space for C and C_tilde
         C_new = -1j * (np.dot(C, eta))
         C_tilde_new = 1j * (-np.dot(eta, C_tilde))
-        """
 
+        """
         C_new = -compute_q_space_ket_equations(
             C,
             C_tilde,
diff --git a/coupled_cluster/tdcc.py b/coupled_cluster/tdcc.py
index 1d172db..4a0bf58 100644
--- a/coupled_cluster/tdcc.py
+++ b/coupled_cluster/tdcc.py
@@ -95,6 +95,7 @@ def amplitudes(self):
         return self._amplitudes
 
     def solve(self, time_points, timestep_tol=1e-8):
+
         n = len(time_points)
 
         for i in range(n - 1):

From 6dff85cb047df1a64f5916410ac090d9e73eedea Mon Sep 17 00:00:00 2001
From: halvarsu <halvard.sutterud@gmail.com>
Date: Thu, 5 Mar 2020 13:19:47 +0100
Subject: [PATCH 13/13] add solve method

---
 coupled_cluster/ccd/oaitdccd.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py
index d7dd2da..b6fa48f 100644
--- a/coupled_cluster/ccd/oaitdccd.py
+++ b/coupled_cluster/ccd/oaitdccd.py
@@ -19,6 +19,38 @@ def compute_ground_state(self, *args, **kwargs):
         self.u = self.system.u
         self.f = self.system.construct_fock_matrix(self.h, self.u)
 
+    def solve(self, time_points, timestep_tol=1e-8):
+
+        n = len(time_points)
+
+        for i in range(n - 1):
+            dt = time_points[i + 1] - time_points[i]
+            amp_vec = self.integrator.step(
+                self._amplitudes.asarray(), time_points[i], dt
+            )
+
+            self._amplitudes = type(self._amplitudes).from_array(
+                self._amplitudes, amp_vec
+            )
+
+            # from coupledcluster.tools import biortonormalize
+
+            # t,l,C,Ctilde = self._amplitudes
+
+            # print("AAA", self.np.linalg.norm(Ctilde @ C - self.np.eye(C.shape[0])))
+            # C, Ctilde = biortonormalize(C, Ctilde)
+            # print("BBB", self.np.linalg.norm(Ctilde @ C - self.np.eye(C.shape[0])))
+
+            # self._amplitudes = OACCVector(t, l, C ,Ctilde, np=self.np)
+
+
+            if abs(self.last_timestep - (time_points[i] + dt)) > timestep_tol:
+                self.update_hamiltonian(time_points[i] + dt, self._amplitudes)
+                self.last_timestep = time_points[i] + dt
+
+            yield self._amplitudes
+
+
     def __call__(self, prev_amp, current_time):
         np = self.np
         o, v = self.o, self.v