From fb07e986c9b251eac32dda63cdc87260f7f14e0b Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Thu, 11 Feb 2021 01:03:24 +0100
Subject: [PATCH] Improving weights

---
 src/cipsi/pt2_stoch_routines.irp.f |  2 +-
 src/cipsi/selection.irp.f          |  4 +--
 src/cipsi/selection_weight.irp.f   | 58 +++++++++++++++++++-----------
 3 files changed, 41 insertions(+), 23 deletions(-)

diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 428285c2..489ffaaf 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -352,9 +352,9 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
 
     state_average_weight(:) = state_average_weight_save(:)
     TOUCH state_average_weight
+    call update_pt2_and_variance_weights(pt2_data, N_states)
   endif
 
-  call update_pt2_and_variance_weights(pt2_data, N_states)
 
 end subroutine
 
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index 5f331fef..5afb514e 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -648,7 +648,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         do_cycle = .True.
         do k=1,N_dominant_dets_of_cfgs
           call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
-          do_cycle = do_cycle .and. (degree > excitation_alpha_max) 
+          do_cycle = do_cycle .and. (degree > excitation_alpha_max)
         enddo
         if (do_cycle) cycle
       endif
@@ -658,7 +658,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         do_cycle = .True.
         do k=1,N_dominant_dets_of_cfgs
           call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int)
-          do_cycle = do_cycle .and. (degree > excitation_beta_max) 
+          do_cycle = do_cycle .and. (degree > excitation_beta_max)
         enddo
         if (do_cycle) cycle
       endif
diff --git a/src/cipsi/selection_weight.irp.f b/src/cipsi/selection_weight.irp.f
index 24e2ace0..fb317f9b 100644
--- a/src/cipsi/selection_weight.irp.f
+++ b/src/cipsi/selection_weight.irp.f
@@ -33,43 +33,59 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st)
 
   double precision :: avg, element, dt, x
   integer          :: k
-  integer, save    :: i_iter=0
-  integer, parameter :: i_itermax = 1
-  double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:)
+!  integer, save    :: i_iter=0
+!  integer, parameter :: i_itermax = 1
+!  double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:)
 
   pt2(:)      = pt2_data % pt2(:)
   variance(:) = pt2_data % variance(:)
 
-  if (i_iter == 0) then
-    allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax))
-    memo_pt2(:,:) = 1.d0
-    memo_variance(:,:) = 1.d0
-  endif
+!  if (i_iter == 0) then
+!    allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax))
+!    memo_pt2(:,:) = 1.d0
+!    memo_variance(:,:) = 1.d0
+!  endif
+!
+!  i_iter = i_iter+1
+!  if (i_iter > i_itermax) then
+!    i_iter = 1
+!  endif
+!
+!  dt = 2.0d0
 
-  i_iter = i_iter+1
-  if (i_iter > i_itermax) then
-    i_iter = 1
-  endif
+  avg = sum(pt2(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero
+!  do k=1,N_st
+!    element = exp(dt*(pt2(k)/avg -1.d0))
+!    element = min(2.0d0 , element)
+!    element = max(0.5d0 , element)
+!    memo_pt2(k,i_iter) = element
+!    pt2_match_weight(k) *= product(memo_pt2(k,:))
+  !enddo
 
-  dt = 2.0d0
-
-  avg = sum(pt2(1:N_st)) / dble(N_st) - 1.d-32 ! Avoid future division by zero
+  dt = 1.0d0 * selection_factor
   do k=1,N_st
-    element = exp(dt*(pt2(k)/avg -1.d0))
+    element = exp(dt*(pt2(k)/avg - 1.d0))
     element = min(2.0d0 , element)
     element = max(0.5d0 , element)
-    memo_pt2(k,i_iter) = element
-    pt2_match_weight(k) *= product(memo_pt2(k,:))
+  print *, k, element
+    pt2_match_weight(k) *= element
   enddo
 
 
   avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero
+!  do k=1,N_st
+!    element = exp(dt*(variance(k)/avg -1.d0))
+!    element = min(2.0d0 , element)
+!    element = max(0.5d0 , element)
+!    memo_variance(k,i_iter) = element
+!    variance_match_weight(k) *= product(memo_variance(k,:))
+!  enddo
+
   do k=1,N_st
     element = exp(dt*(variance(k)/avg -1.d0))
     element = min(2.0d0 , element)
     element = max(0.5d0 , element)
-    memo_variance(k,i_iter) = element
-    variance_match_weight(k) *= product(memo_variance(k,:))
+    variance_match_weight(k) *= element
   enddo
 
   if (N_det < 100) then
@@ -78,6 +94,8 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st)
     variance_match_weight(:) = 1.d0
   endif
 
+print *, 'XXX', n_det, pt2_match_weight(1), pt2_match_weight(2)
+
   threshold_davidson_pt2 = min(1.d-6, &
      max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) )