From 5aed62450e000fa8d1840f4287559c2b314241a7 Mon Sep 17 00:00:00 2001
From: AbdAmmar <abd.ammar.phys@gmail.com>
Date: Wed, 5 Apr 2023 15:59:38 +0200
Subject: [PATCH] print tc energy: OK

---
 src/tc_bi_ortho/print_tc_energy.irp.f | 32 --------------------------
 src/tc_bi_ortho/tc_bi_ortho.irp.f     | 33 ---------------------------
 src/tc_bi_ortho/tc_utils.irp.f        | 32 ++++++++++++++++++++++++++
 3 files changed, 32 insertions(+), 65 deletions(-)
 create mode 100644 src/tc_bi_ortho/tc_utils.irp.f

diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f
index c9f8cecb..e5f123a7 100644
--- a/src/tc_bi_ortho/print_tc_energy.irp.f
+++ b/src/tc_bi_ortho/print_tc_energy.irp.f
@@ -13,35 +13,3 @@ program print_tc_energy
  call write_tc_energy
 end
 
-subroutine write_tc_energy()
-
-  implicit none
-  integer          :: i, j, k
-  double precision :: hmono, htwoe, htot
-  double precision :: E_TC, O_TC
-
-  do k = 1, n_states
-
-    E_TC = 0.d0
-    do i = 1, N_det
-      do j = 1, N_det
-        call hmat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, htot)
-        E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
-        !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
-      enddo
-    enddo
-
-    O_TC = 0.d0
-    do i = 1, N_det
-      O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
-      !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
-    enddo
-
-    print *, ' state :', k
-    print *, " E_TC = ", E_TC
-    print *, " O_TC = ", O_TC
-
-  enddo
-
-end
-
diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f
index ef2e5659..98b83329 100644
--- a/src/tc_bi_ortho/tc_bi_ortho.irp.f
+++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f
@@ -89,36 +89,3 @@ end
 
 
 
-subroutine write_tc_energy()
-
-  implicit none
-  integer          :: i, j, k
-  double precision :: hmono, htwoe, htot
-  double precision :: E_TC, O_TC
-
-  do k = 1, n_states
-
-    E_TC = 0.d0
-    do i = 1, N_det
-      do j = 1, N_det
-        htot = htilde_matrix_elmt_bi_ortho(j,i)
-        !call hmat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, htot)
-        E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
-        !E_TC = E_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(j,k) * htot
-      enddo
-    enddo
-
-    O_TC = 0.d0
-    do i = 1, N_det
-      O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
-      !O_TC = O_TC + leigvec_tc_bi_orth(i,k) * reigvec_tc_bi_orth(i,k)
-    enddo
-
-    print *, ' state :', k
-    print *, " E_TC = ", E_TC
-    print *, " O_TC = ", O_TC
-
-  enddo
-
-end
-
diff --git a/src/tc_bi_ortho/tc_utils.irp.f b/src/tc_bi_ortho/tc_utils.irp.f
new file mode 100644
index 00000000..92e8639d
--- /dev/null
+++ b/src/tc_bi_ortho/tc_utils.irp.f
@@ -0,0 +1,32 @@
+
+subroutine write_tc_energy()
+
+  implicit none
+  integer          :: i, j, k
+  double precision :: hmono, htwoe, hthree, htot
+  double precision :: E_TC, O_TC
+
+  do k = 1, n_states
+
+    E_TC = 0.d0
+    do i = 1, N_det
+      do j = 1, N_det
+        !htot = htilde_matrix_elmt_bi_ortho(i,j)
+        call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,j), N_int, hmono, htwoe, hthree, htot)
+        E_TC = E_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(j,k) * htot
+      enddo
+    enddo
+
+    O_TC = 0.d0
+    do i = 1, N_det
+      O_TC = O_TC + psi_l_coef_bi_ortho(i,k) * psi_r_coef_bi_ortho(i,k)
+    enddo
+
+    print *, ' state :', k
+    print *, " E_TC = ", E_TC / O_TC
+    print *, " O_TC = ", O_TC
+
+  enddo
+
+end
+