From 50b12bcacbaef29de0f3dc9b6d3a128cf803fc3a Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Fri, 25 May 2018 15:03:12 +0200
Subject: [PATCH] L3 weight

---
 plugins/Full_CI_ZMQ/selection.irp.f     | 24 +---------------
 src/Determinants/EZFIO.cfg              |  7 +++++
 src/Determinants/density_matrix.irp.f   | 38 +++++++++++++++++++++----
 src/Determinants/spindeterminants.irp.f |  6 ++--
 4 files changed, 43 insertions(+), 32 deletions(-)

diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f
index 095f7104..516c4567 100644
--- a/plugins/Full_CI_ZMQ/selection.irp.f
+++ b/plugins/Full_CI_ZMQ/selection.irp.f
@@ -1,26 +1,4 @@
 use bitmasks
-BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
- implicit none
- BEGIN_DOC
- ! Weight of the states in the selection : 1/(sum_i c_i^4)
- END_DOC
- integer :: i,k
- double precision :: c
- print *,  'Selection weights'
- print *,  '-----------------'
- do i=1,N_states
-   selection_weight(i) = 0.d0
-   do k=1,N_det
-    c = psi_coef(k,i)*psi_coef(k,i)
-    selection_weight(i) = selection_weight(i) + c*abs(psi_coef(k,i))
-   enddo
-   selection_weight(i) = 1.d0/selection_weight(i)
-   print *,  i, selection_weight(i)
- enddo
- print *,  '-----------------'
-
-END_PROVIDER
-
 
 BEGIN_PROVIDER [ integer, fragment_count ]
   implicit none
@@ -654,7 +632,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         endif
         e_pert = 0.5d0 * (tmp - delta_E)
         pt2(istate) = pt2(istate) + e_pert
-        sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
+        sum_e_pert = sum_e_pert + e_pert * state_average_weight(istate)
       end do
       
       if(sum_e_pert <= buf%mini) then
diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg
index c86f831d..d55dd612 100644
--- a/src/Determinants/EZFIO.cfg
+++ b/src/Determinants/EZFIO.cfg
@@ -40,6 +40,13 @@ doc: Force the wave function to be an eigenfunction of S^2
 interface: ezfio,provider,ocaml
 default: True
 
+[use_l3_weight]
+type: logical
+doc: If true, set the state-averaged weight to 1/(Norm_L3(Psi))
+interface: ezfio,provider,ocaml
+default: False
+
+
 [threshold_generators]
 type: Threshold
 doc: Thresholds on generators (fraction of the norm) 
diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f
index bd5f0741..44d26346 100644
--- a/src/Determinants/density_matrix.irp.f
+++ b/src/Determinants/density_matrix.irp.f
@@ -361,6 +361,28 @@ subroutine save_natural_mos
 end
 
 
+BEGIN_PROVIDER [ double precision, l3_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weight of the states in the selection : 1/(sum_i c_i^4)
+ END_DOC
+ integer :: i,k
+ double precision :: c
+ do i=1,N_states
+   l3_weight(i) = 1.d-31
+   do k=1,N_det
+    c = psi_coef(k,i)*psi_coef(k,i)
+    l3_weight(i) = l3_weight(i) + c*abs(psi_coef(k,i))
+   enddo
+   l3_weight(i) = min(1.d0/l3_weight(i), 100.d0)
+ enddo
+ print *,  'L3 weights'
+ print *,  '----------'
+ print *,  l3_weight(1:N_states)
+
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
  implicit none
  BEGIN_DOC
@@ -368,13 +390,17 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
  END_DOC
  logical :: exists
 
- state_average_weight(:) = 1.d0
- call ezfio_has_determinants_state_average_weight(exists)
- if (exists) then
-  call ezfio_get_determinants_state_average_weight(state_average_weight)
+ if (use_l3_weight) then
+  state_average_weight(:) = l3_weight(:)
+ else
+  state_average_weight(:) = 1.d0
+  call ezfio_has_determinants_state_average_weight(exists)
+  if (exists) then
+    call ezfio_get_determinants_state_average_weight(state_average_weight)
+  endif
+  state_average_weight(:) = state_average_weight(:)+1.d-31
+  state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:)))
  endif
- state_average_weight(:) = state_average_weight(:)+1.d-31
- state_average_weight(:) = state_average_weight(:)/(sum(state_average_weight(:)))
 END_PROVIDER
 
 
diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f
index 136872b4..1f076057 100644
--- a/src/Determinants/spindeterminants.irp.f
+++ b/src/Determinants/spindeterminants.irp.f
@@ -367,13 +367,13 @@ end
    j = psi_bilinear_matrix_columns(k)
    f = 0.d0
    do l=1,N_states
-    f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l)
+    f += psi_bilinear_matrix_values(k,l)*psi_bilinear_matrix_values(k,l) * state_average_weight(l)
    enddo
    det_alpha_norm(i) += f
    det_beta_norm(j)  += f
  enddo
- det_alpha_norm = det_alpha_norm / dble(N_states)
- det_beta_norm = det_beta_norm / dble(N_states)
+ det_alpha_norm = det_alpha_norm 
+ det_beta_norm = det_beta_norm 
 
 END_PROVIDER