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/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 diff --git a/coupled_cluster/ccd/__init__.py b/coupled_cluster/ccd/__init__.py index 4e0652d..9077af1 100644 --- a/coupled_cluster/ccd/__init__.py +++ b/coupled_cluster/ccd/__init__.py @@ -1,4 +1,6 @@ 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 new file mode 100644 index 0000000..9c5313f --- /dev/null +++ b/coupled_cluster/ccd/itdccd.py @@ -0,0 +1,37 @@ +from . import TDCCD +from ..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/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( diff --git a/coupled_cluster/ccd/oaitdccd.py b/coupled_cluster/ccd/oaitdccd.py new file mode 100644 index 0000000..b6fa48f --- /dev/null +++ b/coupled_cluster/ccd/oaitdccd.py @@ -0,0 +1,144 @@ +from . import OATDCCD +from ..oatdcc import ( + compute_q_space_ket_equations, + compute_q_space_bra_equations, +) +from ..cc_helper import OACCVector + + +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 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 + + 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() 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))) 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):