Added switch for multiple selection weights, including variance

src/cipsi/pt2_stoch_routines.irp.f

@@ -333,13 +333,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
     pt2(k) = 0.d0
   enddo
 
-  ! Adjust PT2 weights for next selection
+  call update_pt2_and_variance_weights(pt2, variance, norm)
-  double precision :: pt2_avg
-  pt2_avg = sum(pt2) / dble(N_states)
-  do k=1,N_states
-    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
-  enddo
-  SOFT_TOUCH pt2_match_weight
 
 end subroutine
 

src/cipsi/selection.irp.f

@@ -9,12 +9,73 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
  pt2_match_weight = 1.d0
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weights adjusted along the selection to make the variances 
+ ! of each state coincide.
+ END_DOC
+ variance_match_weight = 1.d0
+END_PROVIDER
+
+subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_states)
+  implicit none
+  BEGIN_DOC
+! Updates the rPT2- and Variance- matching weights.
+  END_DOC
+  integer, intent(in)          :: N_st
+  double precision, intent(in) :: pt2(N_st)
+  double precision, intent(in) :: variance(N_st)
+  double precision, intent(in) :: norm(N_st)
+
+  double precision :: avg, rpt2(N_st)
+  integer          :: k
+
+  do k=1,N_st
+    rpt2(k) = pt2(k)/(1.d0 + norm(k))                                                     
+  enddo
+
+  avg = sum(rpt2(1:N_st)) / dble(N_st)
+  do k=1,N_states
+    pt2_match_weight(k) *= (rpt2(k)/avg)**2
+  enddo
+
+  avg = sum(variance(1:N_st)) / dble(N_st)
+  do k=1,N_states
+    variance_match_weight(k) *= (variance(k)/avg)**2
+  enddo
+
+  SOFT_TOUCH pt2_match_weight variance_match_weight
+end
+
+
 BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
    implicit none
    BEGIN_DOC
    ! Weights used in the selection criterion
    END_DOC
-   selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+   select (weight_selection)
+
+     case (0)
+      selection_weight(1:N_states) = weight_one_e_dm(1:N_states)
+
+     case (1)
+      selection_weight(1:N_states) = c0_weight(1:N_states) 
+
+     case (2)
+      selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (3)
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+     case (4)
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (5)
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+    end select
+
 END_PROVIDER
 
 

src/determinants/EZFIO.cfg

@@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2|
 interface: ezfio,provider,ocaml
 default: True
 
-[used_weight]
+[weight_one_e_dm]
 type: integer
 doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
 interface: ezfio,provider,ocaml
 default: 1
 
+[weight_selection]
+type: integer
+doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching                                  
+interface: ezfio,provider,ocaml
+default: 2
+
 [threshold_generators]
 type: Threshold
 doc: Thresholds on generators (fraction of the square of the norm) 

src/determinants/density_matrix.irp.f

@@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
    logical                        :: exists
 
    state_average_weight(:) = 1.d0
-   if (used_weight == 0) then
+   if (weight_one_e_dm == 0) then
      state_average_weight(:) = c0_weight(:)
-   else if (used_weight == 1) then
+   else if (weight_one_e_dm == 1) then
      state_average_weight(:) = 1./N_states
    else
      call ezfio_has_determinants_state_average_weight(exists)