From 7b0815354d768279ed4f6d5b95e81767f8eb3318 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sun, 27 Jun 2021 15:22:17 +0200
Subject: [PATCH] Fixed qmc_e_curve

---
 devel/qmcchem/qmc_e_curve2.irp.f | 30 ++++++++++++------------------
 1 file changed, 12 insertions(+), 18 deletions(-)

diff --git a/devel/qmcchem/qmc_e_curve2.irp.f b/devel/qmcchem/qmc_e_curve2.irp.f
index 6386bfb..5add73d 100644
--- a/devel/qmcchem/qmc_e_curve2.irp.f
+++ b/devel/qmcchem/qmc_e_curve2.irp.f
@@ -4,6 +4,7 @@ program e_curve
   integer :: i,j,k, kk, nab, m, l
   double precision :: norm, E, hij, num, ci, cj
   double precision :: e_0(N_states)
+  double precision, allocatable :: psi_save(:)
   PROVIDE mo_two_e_integrals_in_map mo_one_e_integrals
 
   if (.not.read_wf) then
@@ -23,38 +24,31 @@ program e_curve
   double precision :: thresh
   integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
   thresh = 1.d-10
+  integer :: n_det_prime
 
   nab = n_det_alpha_unique+n_det_beta_unique
 
-
+  allocate(psi_save(1:N_det))
+  psi_save(1:N_det) = psi_coef(1:N_det,1)
   do while (thresh  < 1.d0)
     norm = 0.d0
+    n_det_prime = n_det
     do k=1,n_det
+      psi_coef(k,1) = psi_save(k)
       if (dabs(psi_coef(k,1)) < thresh) then
         psi_coef(k,1) = 0.d0
+        n_det_prime -= 1
       endif
       norm = norm + psi_coef(k,1)**2
     enddo
-    TOUCH psi_coef
+    TOUCH psi_coef psi_det
     norm = norm/dsqrt(norm)
+    psi_coef(:,1) = psi_coef(:,1)/norm
 
-    psi_coef(1:N_det,1) = psi_coef_sorted(1:N_det,1)
-    psi_det(1:N_int,1:2,1:N_det) = psi_det_sorted(1:N_int,1:2,1:N_det)
-    do k=1,n_det
-      if (psi_coef(k,1) == 0.d0) then
-        exit
-      endif
-    enddo
-    n_det = k-1
-    TOUCH n_det psi_coef psi_det
+    E = psi_energy(1)
 
-    j = n_det_alpha_unique+n_det_beta_unique
-    call u_0_H_u_0(E,psi_coef,n_det,psi_det,N_int,1,size(psi_coef,1))
-
-    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)',  thresh, n_det, &
-      dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
-      dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-j)), norm, &
-       psi_energy(1)
+    print '(E9.1,2X,I8,2X,2X,F10.8,2X,F15.10)',  thresh, n_det_prime, &
+      norm, E
     thresh = thresh * dsqrt(10.d0)
   enddo
   print *,  '=========================================================='