diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f index 8c5f20a5..8e8d495a 100644 --- a/src/Determinants/H_apply.irp.f +++ b/src/Determinants/H_apply.irp.f @@ -189,11 +189,15 @@ subroutine copy_H_apply_buffer_to_wf !$OMP BARRIER H_apply_buffer(j)%N_det = 0 !$OMP END PARALLEL - call normalize(psi_coef,N_det) SOFT_TOUCH N_det psi_det psi_coef logical :: found_duplicates call remove_duplicates_in_psi_det(found_duplicates) + do k=1,N_states + call normalize(psi_coef(1,k),N_det) + enddo + SOFT_TOUCH N_det psi_det psi_coef + end subroutine remove_duplicates_in_psi_det(found_duplicates) @@ -230,7 +234,7 @@ subroutine remove_duplicates_in_psi_det(found_duplicates) j = i+1 do while (bit_tmp(j)==bit_tmp(i)) if (duplicate(j)) then - j += 1 + j = j+1 if (j > N_det) then exit else @@ -246,8 +250,10 @@ subroutine remove_duplicates_in_psi_det(found_duplicates) endif enddo if (dup) then + !$OMP CRITICAL duplicate(j) = .True. found_duplicates = .True. + !$OMP END CRITICAL endif j += 1 if (j > N_det) then @@ -265,17 +271,20 @@ subroutine remove_duplicates_in_psi_det(found_duplicates) k += 1 psi_det(:,:,k) = psi_det_sorted_bit (:,:,i) psi_coef(k,:) = psi_coef_sorted_bit(i,:) -! else -! call debug_det(psi_det_sorted_bit(1,1,i),N_int) -! stop 'duplicates in psi_det' + else + if (sum(abs(psi_coef_sorted_bit(i,:))) /= 0.d0 ) then + psi_coef(k,:) = psi_coef_sorted_bit(i,:) + endif endif enddo N_det = k - call write_bool(6,found_duplicates,'Found duplicate determinants') psi_det_sorted_bit(:,:,1:N_det) = psi_det(:,:,1:N_det) psi_coef_sorted_bit(1:N_det,:) = psi_coef(1:N_det,:) - SOFT_TOUCH N_det psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit + TOUCH N_det psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit endif + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + TOUCH psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit deallocate (duplicate,bit_tmp) end diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 9aa7f631..0797e1c8 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -7,11 +7,8 @@ integer*8 function det_search_key(det,Nint) integer, intent(in) :: Nint integer(bit_kind), intent(in) :: det(Nint,2) integer :: i - det_search_key = iand(det(1,1),det(1,2)) - do i=2,Nint - det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2))) - enddo - det_search_key = iand(huge(det(1,1)),det_search_key) + i = shiftr(elec_alpha_num, bit_kind_shift)+1 + det_search_key = shiftr(ior(det(i,1),det(i,2)),1) end @@ -24,11 +21,8 @@ integer*8 function occ_pattern_search_key(det,Nint) 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 - occ_pattern_search_key = iand(huge(det(1,1)),occ_pattern_search_key) + i = shiftr(elec_alpha_num, bit_kind_shift)+1 + occ_pattern_search_key = shiftr(ior(det(i,1),det(i,2)),1) end diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 668a4cef..3c7a355e 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -336,19 +336,19 @@ end subroutine ! function. END_DOC - call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, & + call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, size(psi_coef,1), & psi_det_sorted_bit, psi_coef_sorted_bit, N_states) END_PROVIDER -subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out, N_st) +subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, sze, det_out, coef_out, N_st) use bitmasks implicit none - integer, intent(in) :: Ndet, N_st - integer(bit_kind), intent(in) :: det_in (N_int,2,psi_det_size) - double precision , intent(in) :: coef_in(psi_det_size,N_st) - integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size) - double precision , intent(out) :: coef_out(psi_det_size,N_st) + integer, intent(in) :: Ndet, N_st, sze + integer(bit_kind), intent(in) :: det_in (N_int,2,sze) + double precision , intent(in) :: coef_in(sze,N_st) + integer(bit_kind), intent(out) :: det_out (N_int,2,sze) + double precision , intent(out) :: coef_out(sze,N_st) BEGIN_DOC ! Determinants are sorted are sorted according to their det_search_key. ! Useful to accelerate the search of a random determinant in the wave diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f index 2649aa70..46f455ca 100644 --- a/src/Determinants/occ_pattern.irp.f +++ b/src/Determinants/occ_pattern.irp.f @@ -3,7 +3,12 @@ subroutine det_to_occ_pattern(d,o,Nint) use bitmasks implicit none BEGIN_DOC - ! Transform a determinant to an occupation pattern + ! Transforms a determinant to an occupation pattern + ! + ! occ(:,1) : Single occupations + ! + ! occ(:,2) : Double occupations + ! END_DOC integer ,intent(in) :: Nint integer(bit_kind),intent(in) :: d(Nint,2) @@ -17,6 +22,7 @@ subroutine det_to_occ_pattern(d,o,Nint) enddo end + subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) use bitmasks implicit none @@ -29,120 +35,77 @@ subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) integer :: amax,bmax,k double precision, external :: binom_func - amax = n_alpha bmax = 0 + amax = n_alpha do k=1,Nint bmax += popcnt( o(k,1) ) amax -= popcnt( o(k,2) ) enddo - sze = int( min(binom_func(bmax, amax), 1.d8) ) - sze = 2*sze*sze + 64 - + sze = binom_int(bmax, amax) 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 ,intent(in) :: Nint + integer ,intent(in) :: n_alpha ! Number of alpha electrons + integer ,intent(inout) :: sze ! Dimension of the output dets + integer(bit_kind),intent(in) :: o(Nint,2) ! Occ patters + integer(bit_kind),intent(out) :: d(Nint,2,sze) ! Output determinants - integer :: i, l, k, nt, na, nd, amax - integer :: list_todo(2*n_alpha) - integer :: list_a(2*n_alpha) - integer :: ishift + integer :: i, k, n, ispin - amax = n_alpha - do k=1,Nint - amax -= popcnt( o(k,2) ) + ! Extract list of singly occupied MOs as (int,pos) pairs + ! ------------------------------------------------------ + + integer :: iint(2*n_alpha), ipos(2*n_alpha) + integer(bit_kind) :: v, t, tt + integer :: n_alpha_in_single + + n=0 + n_alpha_in_single = n_alpha + do i=1,Nint + v = o(i,1) + do while(v /= 0_bit_kind) + n = n+1 + iint(n) = i + ipos(n) = trailz(v) + v = iand(v,v-1) + enddo + n_alpha_in_single = n_alpha_in_single - popcnt( o(i,2) ) enddo - call bitstring_to_list(o(1,1), list_todo, nt, Nint) + v = shiftl(1,n_alpha_in_single) - 1 + sze = binom_int(n,n_alpha_in_single) + do i=1,sze + ! Initialize with doubly occupied MOs + d(:,1,i) = o(:,2) + d(:,2,i) = o(:,2) - 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 - - integer :: ne(2) - l=0 - do i=1,nd - ne(1) = 0 - ne(2) = 0 - l=l+1 - ! Doubly occupied orbitals - do k=1,Nint - d(k,1,l) = ior(d(k,1,i),o(k,2)) - d(k,2,l) = ior(d(k,2,i),o(k,2)) - ne(1) += popcnt(d(k,1,l)) - ne(2) += popcnt(d(k,2,l)) - enddo - if ( (ne(1) /= elec_alpha_num).or.(ne(2) /= elec_beta_num) ) then - l = l-1 - endif - enddo - sze = l - -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) - integer :: iint, ipos, i,j,k - - if (na == amax) then - nd += 1 - if (nd > sze) then - print *, irp_here, ': nd = ', nd - print *, irp_here, ': sze = ', sze - stop 'bug in rec_occ_pattern_to_dets' - endif - do i=1,na - iint = ishft(list_a(i)-1,-bit_kind_shift) + 1 - ipos = list_a(i)-ishft((iint-1),bit_kind_shift)-1 - d(iint,1,nd) = ibset( d(iint,1,nd), ipos ) - enddo - do i=1,nt - iint = ishft(list_todo(i)-1,-bit_kind_shift) + 1 - ipos = list_todo(i)-ishft((iint-1),bit_kind_shift)-1 - d(iint,2,nd) = ibset( d(iint,2,nd), ipos ) - enddo - else - integer, allocatable :: list_todo_tmp(:) - allocate (list_todo_tmp(nt)) - do i=1,nt - if (na > 0) then - if (list_todo(i) < list_a(na)) then - cycle - endif + do k=1,n + if (btest(v,k-1)) then + ispin = 1 + else + ispin = 2 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) + d(iint(k), ispin, i) = ibset( d(iint(k), ispin, i), ipos(k) ) enddo - deallocate(list_todo_tmp) - endif + + ! Generate next permutation with Anderson's algorithm + t = ior(v,v-1) + tt = t+1 + v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) ) + enddo + end + + BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ] &BEGIN_PROVIDER [ integer, N_occ_pattern ] implicit none @@ -348,13 +311,14 @@ subroutine make_s2_eigenfunction 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) ) then - N_det_new += 1 - det_buffer(:,:,N_det_new) = d(:,:,j) - if (N_det_new == bufsze) then - call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread) - N_det_new = 0 - endif + if ( is_in_wavefunction(d(1,1,j), N_int) ) then + cycle + endif + N_det_new += 1 + det_buffer(:,:,N_det_new) = d(:,:,j) + if (N_det_new == bufsze) then + call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread) + N_det_new = 0 endif enddo enddo diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f index 58a71fbc..cdaef347 100644 --- a/src/Determinants/psi_cas.irp.f +++ b/src/Determinants/psi_cas.irp.f @@ -53,7 +53,7 @@ END_PROVIDER ! CAS determinants sorted to accelerate the search of a random determinant in the wave ! function. END_DOC - call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, & + call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, size(psi_cas_coef,1), & psi_cas_sorted_bit, psi_cas_coef_sorted_bit, N_states) END_PROVIDER @@ -106,7 +106,7 @@ END_PROVIDER ! CAS determinants sorted to accelerate the search of a random determinant in the wave ! function. END_DOC - call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, & + call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, size(psi_non_cas_coef,1), & psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit, N_states) END_PROVIDER diff --git a/src/FCI/selection_buffer.irp.f b/src/FCI/selection_buffer.irp.f index 17410b7b..5e57403e 100644 --- a/src/FCI/selection_buffer.irp.f +++ b/src/FCI/selection_buffer.irp.f @@ -137,3 +137,42 @@ subroutine sort_selection_buffer(b) b%cur = nmwen end subroutine +subroutine make_selection_buffer_s2(b) + use selection_types + type(selection_buffer), intent(inout) :: b + + integer(bit_kind), pointer :: o(:,:,:) + double precision, pointer :: val(:) + + integer :: n_d + integer :: i,k,sze,n_alpha,j,n + + ! Create occupation patterns + n_d = 0 + allocate(o(N_int,2,b%cur)) + do i=1,b%cur + do k=1,N_int + o(k,1,i) = ieor(b%det(k,1,i), b%det(k,2,i)) + o(k,2,i) = iand(b%det(k,1,i), b%det(k,2,i)) + enddo + call occ_pattern_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) + n_d = n_d + sze + enddo + deallocate(b%det) + + allocate(b%det(N_int,2,n_d), val(n_d)) + k=1 + do i=1,b%cur + n=n_d + call occ_pattern_to_dets(o(1,1,i),b%det(1,1,k),n,elec_alpha_num,N_int) + val(k) = b%val(i) + do j=k+1,k+n-1 + val(j) = 0.d0 + enddo + k = k+n + enddo + deallocate(o,b%val) + b%val => val + b%N = k-1 + b%cur = k-1 +end diff --git a/src/FCI/zmq_selection.irp.f b/src/FCI/zmq_selection.irp.f index 2a06564d..7cdaee7a 100644 --- a/src/FCI/zmq_selection.irp.f +++ b/src/FCI/zmq_selection.irp.f @@ -120,11 +120,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm) norm(i) = 0.d0 enddo if (N_in > 0) then + if (s2_eig.or.(N_states > 1) ) then + call make_selection_buffer_s2(b) + endif call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) call copy_H_apply_buffer_to_wf() - if (s2_eig.or.(N_states > 1) ) then - call make_s2_eigenfunction - endif call save_wavefunction endif call delete_selection_buffer(b) diff --git a/src/Psiref_Utils/psi_ref_utils.irp.f b/src/Psiref_Utils/psi_ref_utils.irp.f index a63b0ded..20cd7154 100644 --- a/src/Psiref_Utils/psi_ref_utils.irp.f +++ b/src/Psiref_Utils/psi_ref_utils.irp.f @@ -151,7 +151,7 @@ END_PROVIDER ! Reference determinants sorted to accelerate the search of a random determinant in the wave ! function. END_DOC - call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, & + call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, size(psi_non_ref_coef,1), & psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit, N_states) END_PROVIDER diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index 9f4f2168..f11e03a1 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -28,6 +28,27 @@ double precision function binom_func(i,j) end +FUNCTION n_combinations(n,k) RESULT(r) + IMPLICIT NONE + + INTEGER(4), INTENT(IN) :: n,k + INTEGER(8) :: r + INTEGER(4) :: d, n0 + + IF (k > n) THEN + r = 0 + RETURN + ELSE + r = 1 + ENDIF + + n0 = n + DO d=1, k + r = (r*n0) / d + n0 = n0 - 1 + ENDDO +END FUNCTION n_combinations + BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ] &BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ] implicit none @@ -45,6 +66,22 @@ end END_PROVIDER + BEGIN_PROVIDER [ integer*8, binom_int, (0:40,0:40) ] +&BEGIN_PROVIDER [ integer*8, binom_int_transp, (0:40,0:40) ] + implicit none + BEGIN_DOC + ! Binomial coefficients, as integers*8 + END_DOC + integer :: k,l + double precision :: logfact + do l=0,40 + do k=0,40 + binom_int(k,l) = int(binom(k,l)+0.1d0,8) + enddo + enddo +END_PROVIDER + + double precision function fact(n) implicit none