From b430721eb42c2a36bc16f33e5fc5812752057990 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Sun, 27 Jun 2021 23:45:30 +0200
Subject: [PATCH] Estimate cost

---
 devel/qmcchem/qmc_e_curve.irp.f  | 16 ++++++++++++++--
 devel/qmcchem/qmc_e_curve2.irp.f | 24 +++++++++++++++++++++---
 2 files changed, 35 insertions(+), 5 deletions(-)

diff --git a/devel/qmcchem/qmc_e_curve.irp.f b/devel/qmcchem/qmc_e_curve.irp.f
index e6135b0..f9bc210 100644
--- a/devel/qmcchem/qmc_e_curve.irp.f
+++ b/devel/qmcchem/qmc_e_curve.irp.f
@@ -45,16 +45,21 @@ program e_curve
   print *,  '=========================================================='
   double precision :: thresh
   integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  integer :: na, nb
   thresh = 1.d-10
+  na = n_det_alpha_unique
+  nb = n_det_beta_unique
   do j=0,nab
     i = iorder(j)
     if (i<0) then
+      nb -= 1
       do k=1,n_det
         if (psi_bilinear_matrix_columns(k) == -i) then
           psi_bilinear_matrix_values(k,1) = 0.d0
         endif
       enddo
     else
+      na -= 1
       do k=1,n_det
         if (psi_bilinear_matrix_rows(k) ==  i) then
           psi_bilinear_matrix_values(k,1) = 0.d0
@@ -82,11 +87,18 @@ program e_curve
        exit
     endif
     E = E_0(1) + nuclear_repulsion
+
+    double precision :: cost0, cost
+    cost0 = elec_alpha_num**3 + elec_beta_num**3
+    cost = (na-1) * elec_alpha_num**2 + &
+           (nb-1) * elec_beta_num**2 + &
+            elec_alpha_num**3 + elec_beta_num**3
+    cost = cost/cost0
+
     double precision :: u_dot_u
     norm = dsqrt(u_dot_u(psi_coef(1,1),N_det))
     print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)',  thresh, m, &
-      dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
-      dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, psi_energy(1)
+      cost,  norm, psi_energy(1)
     thresh = thresh * dsqrt(10.d0)
   enddo
   print *,  '=========================================================='
diff --git a/devel/qmcchem/qmc_e_curve2.irp.f b/devel/qmcchem/qmc_e_curve2.irp.f
index 5add73d..5ae388d 100644
--- a/devel/qmcchem/qmc_e_curve2.irp.f
+++ b/devel/qmcchem/qmc_e_curve2.irp.f
@@ -45,10 +45,28 @@ program e_curve
     norm = norm/dsqrt(norm)
     psi_coef(:,1) = psi_coef(:,1)/norm
 
+    integer :: na, nb
+    na = 0
+    do k=1,n_det_alpha_unique
+      if (det_alpha_norm(k) > 0.d0) then
+        na += 1
+      endif
+    enddo
+    nb = 0
+    do k=1,n_det_beta_unique
+      if (det_beta_norm(k) > 0.d0) then
+        nb += 1
+      endif
+    enddo
     E = psi_energy(1)
-
-    print '(E9.1,2X,I8,2X,2X,F10.8,2X,F15.10)',  thresh, n_det_prime, &
-      norm, E
+    double precision :: cost0, cost
+    cost0 = elec_alpha_num**3 + elec_beta_num**3
+    cost = (na-1) * elec_alpha_num**2 + &
+           (nb-1) * elec_beta_num**2 + &
+            elec_alpha_num**3 + elec_beta_num**3
+    cost = cost/cost0
+    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F15.10)',  thresh, n_det_prime, &
+      cost, norm, E
     thresh = thresh * dsqrt(10.d0)
   enddo
   print *,  '=========================================================='