From 85a3168153171722a741753ce89aa31eebbf026c Mon Sep 17 00:00:00 2001 From: Pavlo Pelikh <51415689+Pavlo3P@users.noreply.github.com> Date: Wed, 15 May 2024 00:50:56 +0000 Subject: [PATCH 1/4] Added T2' condition --- pyci/rdm/constraints.py | 75 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 73 insertions(+), 2 deletions(-) diff --git a/pyci/rdm/constraints.py b/pyci/rdm/constraints.py index f26fe08..7d87e87 100644 --- a/pyci/rdm/constraints.py +++ b/pyci/rdm/constraints.py @@ -309,6 +309,77 @@ def calc_T2(gamma, N, conjugate=False): a_bar - (term_6 - term_7 - term_8 + term_9) -def calc_T2_prime(): - pass +def calc_T2_prime(gamma, N, conjugate=False): + """ + Calculating T2' tensor + + Parameters + ---------- + gamma: np.ndarray + 1DM tensor + N: int + number of electrons in the system + conjugate: bool + conjugate or regular condition + + Returns + ------- + np.ndarray + Notes + ----- + T2' is defined as: + + .. math:: + \begin{aligned} + \mathcal{T}'_{2}(\Gamma) + = \left( \begin{matrix} + \mathcal{T}_{2}(\Gamma)_{\alpha \beta \gamma; \delta \epsilon \zeta} & (\Gamma_\omega)_{\alpha \beta \gamma; \nu} \\ + (\Gamma_\omega^{\dagger})_{\mu; \delta \epsilon \zeta} & (\Gamma_\rho)_{\mu \nu} + \end{matrix} \right) + \end{aligned} + + \begin{aligned} + \mathcal{T}'^{\dagger}_{2}(\Gamma)_{\alpha \beta; \gamma \delta} + = & + \mathcal{T}^{\dagger}_{2}(A_{\mathcal{T}}) + + (\Gamma_{\omega})_{\alpha \beta \delta; \gamma} + + (\Gamma_{\omega})_{\gamma \delta \beta; \alpha} + \\ + & + - (\Gamma_{\omega})_{\alpha \beta \gamma; \delta} + - (\Gamma_{\omega})_{\gamma \delta \alpha; \beta} + \\ + & + + \frac{1}{N-1} + \left( + \delta_{\beta \delta} (\Gamma_{\rho})_{\gamma \alpha} + - \delta_{\alpha \delta} (\Gamma_{\rho})_{\gamma \beta} + - \delta_{\beta \gamma} (\Gamma_{\rho})_{\delta \alpha} + + \delta_{\alpha \gamma} (\Gamma_{\rho})_{\delta \beta} + \right) + \end{aligned} + """ + omega = np.einsum('abgd -> abdg', gamma) + rho = 1/(N-1) * np.einsum('abgb -> ag', gamma) + + if not conjugate: + t2 = calc_T2(gamma, N, False) + omega_d = omega.conj().T + n = t2.shape[0] + return np.block([ + [t2.reshape(n**3, n**3), omega.reshape((n**3, n))], + [omega_d.reshape(n, n**3), rho] + ]) + else: + t2 = calc_T2(gamma, N, False) + t2_d = calc_T2(t2, N, True) + eye = np.eye(N) + + term1 = np.einsum('abgd, abdg, gdba, abgd, gdab -> abgd', t2_d, omega, omega, -omega, -omega) + term2 = np.einsum('bd, ga -> abgd', eye, rho) + term3 = np.einsum('ad, gb -> abgd', eye, rho) + term4 = np.einsum('bg, da -> abgd', eye, rho) + term5 = np.einsum('ag, db -> abgd', eye, rho) + + return term1 + 1 / (N - 1) * (term2 - term3 - term4 + term5) From 25c900c9433d313c867191f3c3ad49c9d69187d5 Mon Sep 17 00:00:00 2001 From: Pavlo Pelikh <51415689+Pavlo3P@users.noreply.github.com> Date: Wed, 15 May 2024 01:29:27 +0000 Subject: [PATCH 2/4] Added P condition --- pyci/rdm/constraints.py | 26 +++++++++++++++++++++++--- 1 file changed, 23 insertions(+), 3 deletions(-) diff --git a/pyci/rdm/constraints.py b/pyci/rdm/constraints.py index 7d87e87..908d9f5 100644 --- a/pyci/rdm/constraints.py +++ b/pyci/rdm/constraints.py @@ -78,8 +78,28 @@ def calculate_shift(eigenvalues, alpha): res = root(constraint, 0) return res.x -def calc_P(): - pass +def calc_P(gamma, N, conjugate=False): + """ + Calculating P tensor + + Parameters + ---------- + gamma: np.ndarray + 1DM tensor + N: int + number of electrons in the system + conjugate: bool + conjugate or regular condition + + Returns + ------- + np.ndarray + + Notes + ----- + + """ + return gamma def calc_Q(): pass @@ -381,5 +401,5 @@ def calc_T2_prime(gamma, N, conjugate=False): term3 = np.einsum('ad, gb -> abgd', eye, rho) term4 = np.einsum('bg, da -> abgd', eye, rho) term5 = np.einsum('ag, db -> abgd', eye, rho) - + return term1 + 1 / (N - 1) * (term2 - term3 - term4 + term5) From c3c82181aa78a1a26426c583adacc7f915b89162 Mon Sep 17 00:00:00 2001 From: Pavlo Pelikh <51415689+Pavlo3P@users.noreply.github.com> Date: Thu, 16 May 2024 21:18:21 +0000 Subject: [PATCH 3/4] Fixed T2' --- pyci/rdm/constraints.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/pyci/rdm/constraints.py b/pyci/rdm/constraints.py index 908d9f5..c2c71b3 100644 --- a/pyci/rdm/constraints.py +++ b/pyci/rdm/constraints.py @@ -396,10 +396,14 @@ def calc_T2_prime(gamma, N, conjugate=False): t2_d = calc_T2(t2, N, True) eye = np.eye(N) - term1 = np.einsum('abgd, abdg, gdba, abgd, gdab -> abgd', t2_d, omega, omega, -omega, -omega) - term2 = np.einsum('bd, ga -> abgd', eye, rho) - term3 = np.einsum('ad, gb -> abgd', eye, rho) - term4 = np.einsum('bg, da -> abgd', eye, rho) - term5 = np.einsum('ag, db -> abgd', eye, rho) - - return term1 + 1 / (N - 1) * (term2 - term3 - term4 + term5) + term1 = t2_d + term2 = np.einsum('abdg -> abgd', omega) + term3 = np.einsum('gdba -> abgd', omega) + term4 = omega + term5 = np.einsum('gdab -> abgd', omega) + term6 = np.einsum('bd, ga -> abgd', eye, rho) + term7 = np.einsum('ad, gb -> abgd', eye, rho) + term8 = np.einsum('bg, da -> abgd', eye, rho) + term9 = np.einsum('ag, db -> abgd', eye, rho) + + return term1 + term2 + term3 - term4 - term5 + 1 / (N - 1) * (term6 - term7 - term8 + term9) From 39ff1c36553ec5962e6977a57a4ab0707e590aef Mon Sep 17 00:00:00 2001 From: Pavlo Pelikh <51415689+Pavlo3P@users.noreply.github.com> Date: Sun, 19 May 2024 16:31:15 +0000 Subject: [PATCH 4/4] Fixed T2' Used formula (2.67) from Ward Poelmans' thesis for computing T2', flat_tensor function replaced reshaping in the formula. --- pyci/rdm/constraints.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/pyci/rdm/constraints.py b/pyci/rdm/constraints.py index c2c71b3..b636699 100644 --- a/pyci/rdm/constraints.py +++ b/pyci/rdm/constraints.py @@ -18,6 +18,7 @@ import numpy as np from scipy.optimize import root +from pyci.rdm.tools import flat_tensor __all__ = [ @@ -355,7 +356,7 @@ def calc_T2_prime(gamma, N, conjugate=False): \mathcal{T}'_{2}(\Gamma) = \left( \begin{matrix} \mathcal{T}_{2}(\Gamma)_{\alpha \beta \gamma; \delta \epsilon \zeta} & (\Gamma_\omega)_{\alpha \beta \gamma; \nu} \\ - (\Gamma_\omega^{\dagger})_{\mu; \delta \epsilon \zeta} & (\Gamma_\rho)_{\mu \nu} + (\Gamma_\omega)_{\mu; \delta \epsilon \zeta} & (\Gamma_\rho)_{\mu \nu} \end{matrix} \right) \end{aligned} @@ -385,11 +386,10 @@ def calc_T2_prime(gamma, N, conjugate=False): if not conjugate: t2 = calc_T2(gamma, N, False) - omega_d = omega.conj().T n = t2.shape[0] - return np.block([ - [t2.reshape(n**3, n**3), omega.reshape((n**3, n))], - [omega_d.reshape(n, n**3), rho] + return np.block([ + [flat_tensor(t2, (n**3, n**3)), flat_tensor(omega, (n**3, n))], + [flat_tensor(omega, (n, n**3)), rho] ]) else: t2 = calc_T2(gamma, N, False)