diff --git a/data/ezfio_defaults b/data/ezfio_defaults index be7e7605..940b5e26 100644 --- a/data/ezfio_defaults +++ b/data/ezfio_defaults @@ -19,9 +19,10 @@ determinants n_states 1 n_states_diag determinants_n_states n_det_max_jacobi 5000 - threshold_generators 0.995 + threshold_generators 0.99 threshold_selectors 0.999 read_wf False + s2_eig True full_ci n_det_max_fci 1000000 diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 1b7e816c..4c0401c9 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -237,9 +237,9 @@ class H_apply(object): !$ call omp_set_lock(lck) do k=1,N_st norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k) - delta_pt2(k) = 0.d0 - pt2_old(k) = 0.d0 - pt2(k) = select_max(i_generator) +! delta_pt2(k) = 0.d0 +! pt2_old(k) = 0.d0 +! pt2(k) = select_max(i_generator) enddo !$ call omp_unset_lock(lck) cycle diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 57ed7e01..3a95424c 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -506,11 +506,14 @@ subroutine $subroutine($params_main) !$OMP END PARALLEL !$ call omp_destroy_lock(lck) + abort_here = abort_all call stop_progress $copy_buffer + if (s2_eig) then + call make_s2_eigenfunction + endif $generate_psi_guess - abort_here = abort_all end diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index 2f01a8e0..83d7e297 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/src/Dets/connected_to_ref.irp.f @@ -13,6 +13,23 @@ integer*8 function det_search_key(det,Nint) enddo end + +integer*8 function occ_pattern_search_key(det,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Return an integer*8 corresponding to a determinant index for searching + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det(Nint,2) + integer :: i + occ_pattern_search_key = ieor(det(1,1),det(1,2)) + do i=2,Nint + occ_pattern_search_key = ieor(occ_pattern_search_key,iand(det(i,1),det(i,2))) + enddo +end + + logical function is_in_wavefunction(key,Nint,Ndet) implicit none diff --git a/src/Dets/determinants.ezfio_config b/src/Dets/determinants.ezfio_config index 4641db0f..38ff88c9 100644 --- a/src/Dets/determinants.ezfio_config +++ b/src/Dets/determinants.ezfio_config @@ -15,4 +15,5 @@ determinants det_coef double precision (determinants_det_num) read_wf logical expected_s2 double precision + s2_eig logical diff --git a/src/Dets/determinants.irp.f b/src/Dets/determinants.irp.f index 073654ba..9f2e7e67 100644 --- a/src/Dets/determinants.irp.f +++ b/src/Dets/determinants.irp.f @@ -101,6 +101,117 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ] END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_occ_pattern ] + implicit none + BEGIN_DOC + ! array of the occ_pattern present in the wf + ! psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupation + ! psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation + END_DOC + integer :: i,j,k + + ! create + do i = 1, N_det + do k = 1, N_int + psi_occ_pattern(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i)) + psi_occ_pattern(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i)) + enddo + enddo + + ! Sort + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: occ_pattern_search_key + integer(bit_kind), allocatable :: tmp_array(:,:,:) + logical,allocatable :: duplicate(:) + + + allocate ( iorder(N_det), duplicate(N_det), bit_tmp(N_det), tmp_array(N_int,2,psi_det_size) ) + + do i=1,N_det + iorder(i) = i + !$DIR FORCEINLINE + bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + enddo + call i8sort(bit_tmp,iorder,N_det) + !DIR$ IVDEP + do i=1,N_det + do k=1,N_int + tmp_array(k,1,i) = psi_occ_pattern(k,1,iorder(i)) + tmp_array(k,2,i) = psi_occ_pattern(k,2,iorder(i)) + enddo + duplicate(i) = .False. + enddo + + i=1 + integer (bit_kind) :: occ_pattern_tmp + do i=1,N_det + duplicate(i) = .False. + enddo + + do i=1,N_det + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j+=1 + cycle + endif + duplicate(j) = .True. + do k=1,N_int + if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & + .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then + duplicate(j) = .False. + exit + endif + enddo + j+=1 + if (j>N_det) then + exit + endif + enddo + enddo + + N_occ_pattern=0 + do i=1,N_det + if (duplicate(i)) then + cycle + endif + N_occ_pattern += 1 + do k=1,N_int + psi_occ_pattern(k,1,N_occ_pattern) = tmp_array(k,1,i) + psi_occ_pattern(k,2,N_occ_pattern) = tmp_array(k,2,i) + enddo + enddo + + deallocate(iorder,duplicate,bit_tmp,tmp_array) +! !TODO DEBUG +! integer :: s +! do i=1,N_occ_pattern +! do j=i+1,N_occ_pattern +! s = 0 +! do k=1,N_int +! if((psi_occ_pattern(k,1,j) /= psi_occ_pattern(k,1,i)).or. & +! (psi_occ_pattern(k,2,j) /= psi_occ_pattern(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error : occ ', j, 'already in wf' +! call debug_det(psi_occ_pattern(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG +END_PROVIDER + + subroutine read_dets(det,Nint,Ndet) use bitmasks implicit none @@ -118,7 +229,6 @@ subroutine read_dets(det,Nint,Ndet) integer :: i,k equivalence (det_8, det_bk) -! print*,'coucou' call ezfio_get_determinants_N_int(N_int2) ASSERT (N_int2 == Nint) call ezfio_get_determinants_bit_kind(k) @@ -145,7 +255,6 @@ subroutine read_dets(det,Nint,Ndet) enddo enddo deallocate(psi_det_read) -! print*,'ciao' end @@ -278,7 +387,7 @@ END_PROVIDER enddo enddo - deallocate(iorder) + deallocate(iorder, bit_tmp) END_PROVIDER diff --git a/src/Dets/occ_pattern.irp.f b/src/Dets/occ_pattern.irp.f new file mode 100644 index 00000000..33a9dc75 --- /dev/null +++ b/src/Dets/occ_pattern.irp.f @@ -0,0 +1,140 @@ +subroutine det_to_occ_pattern(d,o,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Transform a determinant to an occupation pattern + END_DOC + integer ,intent(in) :: Nint + integer(bit_kind),intent(in) :: d(Nint,2) + integer(bit_kind),intent(out) :: o(Nint,2) + + integer :: k + + do k=1,Nint + o(k,1) = ieor(d(k,1),d(k,2)) + o(k,2) = iand(d(k,1),d(k,2)) + enddo +end + +subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Number of possible determinants for a given occ_pattern + END_DOC + integer ,intent(in) :: Nint, n_alpha + integer(bit_kind),intent(in) :: o(Nint,2) + integer, intent(out) :: sze + integer :: amax,bmax,k + double precision, external :: binom_func + + amax = n_alpha + bmax = 0 + do k=1,Nint + bmax += popcnt( o(k,1) ) + amax -= popcnt( o(k,2) ) + enddo + sze = int( min(binom_func(bmax, amax), 1.d8) ) + +end + +subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Generate all possible determinants for a give occ_pattern + END_DOC + integer ,intent(in) :: Nint, n_alpha + integer ,intent(inout) :: sze + integer(bit_kind),intent(in) :: o(Nint,2) + integer(bit_kind),intent(out) :: d(Nint,2,sze) + + integer :: i, k, nt, na, nd, amax + integer :: list_todo(n_alpha) + integer :: list_a(n_alpha) + + amax = n_alpha + do k=1,Nint + amax -= popcnt( o(k,2) ) + enddo + + call bitstring_to_list(o(1,1), list_todo, nt, Nint) + + na = 0 + nd = 0 + d = 0 + call rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint) + + sze = nd + + do i=1,nd + ! Doubly occupied orbitals + do k=1,Nint + d(k,1,i) = ior(d(k,1,i),o(k,2)) + d(k,2,i) = ior(d(k,2,i),o(k,2)) + enddo + enddo + +! !TODO DEBUG +! integer :: j,s +! do i=1,nd +! do j=1,i-1 +! na=0 +! do k=1,Nint +! if((d(k,1,j) /= d(k,1,i)).or. & +! (d(k,2,j) /= d(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( j== 0 ) then +! print *, 'det ',i,' and ',j,' equal:' +! call debug_det(d(1,1,j),Nint) +! call debug_det(d(1,1,i),Nint) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG +end + +recursive subroutine rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint) + use bitmasks + implicit none + + integer, intent(in) :: nt, sze, amax, Nint,na + integer,intent(inout) :: list_todo(nt) + integer, intent(inout) :: list_a(na+1),nd + integer(bit_kind),intent(inout) :: d(Nint,2,sze) + + if (na == amax) then + nd += 1 + if (na > 0) then + call list_to_bitstring( d(1,1,nd), list_a, na, Nint) + endif + if (nt > 0) then + call list_to_bitstring( d(1,2,nd), list_todo, nt, Nint) + endif + else + integer :: i, j, k + integer :: list_todo_tmp(nt) + do i=1,nt + if (na > 0) then + if (list_todo(i) < list_a(na)) then + cycle + endif + endif + list_a(na+1) = list_todo(i) + k=1 + do j=1,nt + if (i/=j) then + list_todo_tmp(k) = list_todo(j) + k += 1 + endif + enddo + call rec_occ_pattern_to_dets(list_todo_tmp,nt-1,list_a,na+1,d,nd,sze,amax,Nint) + enddo + endif + +end + diff --git a/src/Dets/options.irp.f b/src/Dets/options.irp.f index 3f35ba4d..0ff351ae 100644 --- a/src/Dets/options.irp.f +++ b/src/Dets/options.irp.f @@ -22,6 +22,11 @@ T.set_ezfio_name( "read_wf" ) T.set_output ( "output_dets" ) print T +T.set_name ( "s2_eig" ) +T.set_doc ( "Force the wave function to be an eigenfunction of S^2" ) +T.set_ezfio_name( "s2_eig" ) +print T + END_SHELL BEGIN_PROVIDER [ integer, N_states_diag ] diff --git a/src/Dets/slater_rules.irp.f b/src/Dets/slater_rules.irp.f index 7076df2f..56fff5bd 100644 --- a/src/Dets/slater_rules.irp.f +++ b/src/Dets/slater_rules.irp.f @@ -960,7 +960,7 @@ BEGIN_PROVIDER [ integer*8, det_connections, (N_con_int,N_det) ] integer :: j_int, j_k, j_l integer, allocatable :: idx(:) integer :: thread_num - !$ integer :: omp_get_thread_num + integer :: omp_get_thread_num PROVIDE progress_bar call start_progress(N_det,'Det connections',0.d0) diff --git a/src/Perturbation/selection.irp.f b/src/Perturbation/selection.irp.f index 9420d482..96f0e876 100644 --- a/src/Perturbation/selection.irp.f +++ b/src/Perturbation/selection.irp.f @@ -134,4 +134,88 @@ subroutine remove_small_contributions endif end +subroutine make_s2_eigenfunction + implicit none + integer :: i,j,k + integer :: smax, s + integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:) + integer :: N_det_new + integer, parameter :: bufsze = 1000 + logical, external :: is_in_wavefunction + +! !TODO DEBUG +! do i=1,N_det +! do j=i+1,N_det +! s = 0 +! do k=1,N_int +! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & +! (psi_det(k,2,j) /= psi_det(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error0: det ', j, 'already in wf' +! call debug_det(psi_det(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG + + allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) ) + smax = 1 + N_det_new = 0 + + do i=1,N_occ_pattern + call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int) + if (s > smax) then + deallocate(d) + allocate ( d(N_int,2,s) ) + endif + call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) + do j=1,s + if (.not. is_in_wavefunction( d(1,1,j), N_int, N_det)) then + N_det_new += 1 + do k=1,N_int + det_buffer(k,1,N_det_new) = d(k,1,j) + det_buffer(k,2,N_det_new) = d(k,2,j) + enddo + if (N_det_new == bufsze) then + call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0) + N_det_new = 0 + endif + endif + enddo + enddo + + if (N_det_new > 0) then + call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0) + endif + + deallocate(d,det_buffer) + call copy_H_apply_buffer_to_wf + +! !TODO DEBUG +! do i=1,N_det +! do j=i+1,N_det +! s = 0 +! do k=1,N_int +! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & +! (psi_det(k,2,j) /= psi_det(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error : det ', j, 'already in wf at ', i +! call debug_det(psi_det(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG + call write_int(output_dets,N_det_new, 'Added deteminants for S^2') + +end