diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f index 05ab7b07..7c48c57b 100644 --- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f +++ b/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)) diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 35898a46..3febff69 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/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 diff --git a/plugins/Perturbation/EZFIO.cfg b/plugins/Perturbation/EZFIO.cfg index 485e15cd..89125b37 100644 --- a/plugins/Perturbation/EZFIO.cfg +++ b/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 diff --git a/plugins/shiftedbk/EZFIO.cfg b/plugins/shiftedbk/EZFIO.cfg index 6069c855..77e97cf7 100644 --- a/plugins/shiftedbk/EZFIO.cfg +++ b/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 diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index 8250823a..016b546b 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/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 diff --git a/src/Determinants/truncate_wf.irp.f b/src/Determinants/truncate_wf.irp.f index 320e07c2..619489aa 100644 --- a/src/Determinants/truncate_wf.irp.f +++ b/src/Determinants/truncate_wf.irp.f @@ -1,12 +1,17 @@ program s2_eig_restart implicit none read_wf = .True. - call routine + 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