Truncate eigenfnuction of S^2

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -414,6 +414,7 @@ end function
       pt2_J(i) = i
   end do
 
+  integer :: m
   integer, allocatable :: seed(:)
   call random_seed(size=m)
   allocate(seed(m))

plugins/Full_CI_ZMQ/selection.irp.f

@@ -412,10 +412,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
   integer :: nb_count, maskInd_save, monoBdo_save
   logical :: found
-
   do s1=1,2
     do i1=N_holes(s1),1,-1   ! Generate low excitations first
 
+      found = .False.
       monoBdo_save = monoBdo
       maskInd_save = maskInd
       do s2=s1,2

plugins/Perturbation/EZFIO.cfg

@@ -15,7 +15,7 @@ default: 0.0001
 type: Normalized_float
 doc: Stop stochastic PT2 when the relative error is smaller than PT2_relative_error
 interface: ezfio,provider,ocaml
-default: 0.001
+default: 0.01
 
 [correlation_energy_ratio_max]
 type: Normalized_float

plugins/shiftedbk/EZFIO.cfg

@@ -30,5 +30,5 @@ default: EN
 type: Normalized_float
 doc: Stop stochastic dressing when the relative error is smaller than PT2_relative_error
 interface: ezfio,provider,ocaml
-default: 0.001
+default: 0.01
 

src/Determinants/occ_pattern.irp.f

@@ -248,6 +248,52 @@ end
 
 END_PROVIDER 
 
+BEGIN_PROVIDER [ integer, det_to_occ_pattern, (N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Returns the index of the occupation pattern for each determinant
+ END_DOC
+ integer :: i,j,k
+ integer(bit_kind) :: occ(N_int,2)
+ logical :: found
+ do i=1,N_det
+    do k = 1, N_int
+      occ(k,1) = ieor(psi_det(k,1,i),psi_det(k,2,i))
+      occ(k,2) = iand(psi_det(k,1,i),psi_det(k,2,i))
+    enddo
+    do j=1,N_occ_pattern
+      found = .True.
+      do k=1,N_int
+        if ( (occ(k,1) /= psi_occ_pattern(k,1,j)) &
+        .or. (occ(k,2) /= psi_occ_pattern(k,2,j)) ) then
+          found = .False.
+          exit
+        endif
+      enddo
+      if (found) then
+        det_to_occ_pattern(i) = j
+        exit
+      endif
+    enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, weight_occ_pattern, (N_occ_pattern,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weight of the occupation patterns in the wave function
+ END_DOC
+ integer :: i,j,k
+ weight_occ_pattern = 0.d0
+ do i=1,N_det
+  j = det_to_occ_pattern(i)
+  do k=1,N_states
+    weight_occ_pattern(j,k) += psi_coef(i,k) * psi_coef(i,k)
+  enddo
+ enddo
+END_PROVIDER
+
+
 subroutine make_s2_eigenfunction
   implicit none
   integer                        :: i,j,k

src/Determinants/truncate_wf.irp.f

@@ -1,12 +1,17 @@
 program s2_eig_restart
  implicit none
  read_wf = .True.
+ if (s2_eig) then
+  call routine_s2
+ else
   call routine
+ endif
 end
+
 subroutine routine
  implicit none
  integer :: ndet_max
- print*, 'How many determinants would you like ?'
+ print*, 'Max number of determinants ?'
  read(5,*)ndet_max
  integer(bit_kind), allocatable :: psi_det_tmp(:,:,:)
  double precision, allocatable :: psi_coef_tmp(:,:)
@@ -37,3 +42,50 @@ subroutine routine
  call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,N_det_max,psi_coef_tmp)
  
 end
+
+subroutine routine_s2
+ implicit none
+ integer :: ndet_max
+ double precision :: wmin
+ integer(bit_kind), allocatable :: psi_det_tmp(:,:,:)
+ double precision, allocatable :: psi_coef_tmp(:,:)
+ integer :: i,j,k
+ double precision :: accu(N_states)
+
+ print*, 'Min weight of the occupation pattern ?'
+ read(5,*) wmin
+
+ ndet_max = 0
+ do i=1,N_det
+   if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle
+   ndet_max = ndet_max+1
+ enddo
+
+ allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states))
+ 
+ accu = 0.d0
+ k=0
+ do i = 1, N_det
+  if (maxval(weight_occ_pattern( det_to_occ_pattern(i),:)) < wmin) cycle
+  k = k+1
+  do j = 1, N_int
+   psi_det_tmp(j,1,k) = psi_det(j,1,i)
+   psi_det_tmp(j,2,k) = psi_det(j,2,i)
+  enddo
+  do j = 1, N_states
+   psi_coef_tmp(k,j) = psi_coef(i,j)
+   accu(j) += psi_coef_tmp(k,j) **2
+  enddo
+ enddo
+ do j = 1, N_states
+  accu(j) = 1.d0/dsqrt(accu(j))
+ enddo
+ do j = 1, N_states
+  do i = 1, ndet_max
+   psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j)
+  enddo
+ enddo
+
+ call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,N_det_max,psi_coef_tmp)
+ 
+end