From 33bd506328796e9008bf25cbf0202570bfe61291 Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Fri, 20 May 2016 09:44:22 +0200 Subject: [PATCH] working mrcc --- config/gfortran_debug.cfg | 2 +- plugins/MRCC_Utils/davidson.irp.f | 2 +- plugins/MRCC_Utils/mrcc_dress.irp.f | 78 +-- plugins/MRCC_Utils/mrcc_utils.irp.f | 455 ++++++++++++++- plugins/mrcepa0/dressing.irp.f | 521 ++++++------------ .../mrcepa0/{mrsc2sub.irp.f => mrcc.irp.f} | 2 +- plugins/mrcepa0/mrcepa0.irp.f | 2 +- plugins/mrcepa0/mrcepa0_general.irp.f | 1 + plugins/mrcepa0/mrsc2.irp.f | 3 +- src/Determinants/H_apply.template.f | 4 +- src/Determinants/connected_to_ref.irp.f | 25 +- src/Determinants/determinants.irp.f | 41 ++ src/Determinants/slater_rules.irp.f | 10 +- src/Determinants/spindeterminants.irp.f | 20 +- 14 files changed, 733 insertions(+), 433 deletions(-) rename plugins/mrcepa0/{mrsc2sub.irp.f => mrcc.irp.f} (88%) diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 72084241..03663eea 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -51,7 +51,7 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [DEBUG] -FCFLAGS : -g -pedantic -msse4.2 +FCFLAGS : -g -msse4.2 # OpenMP flags ################# diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index 15077481..66f4975a 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -269,7 +269,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin to_print(2,k) = residual_norm(k) enddo - write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st) + write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))') iter, to_print(:,1:N_st) call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) if (converged) then exit diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 1c2e8b74..412c52e2 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -51,9 +51,9 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:) integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) integer :: mobiles(2), smallerlist + logical, external :: is_generable - - + print *, i_generator leng = max(N_det_generators, N_det_non_ref) allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref)) @@ -69,7 +69,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge allocate( microlist(Nint,2,N_minilist*4), & idx_microlist(N_minilist*4)) - if(key_mask(1,1) /= 0) then + if(key_mask(1,1) /= 0_8) then call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) else @@ -87,6 +87,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge ! |alpha> if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint) if(N_minilist == 0) return @@ -117,8 +118,18 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge - do i_alpha=1,N_tq +! ok = .false. +! do i=N_det_generators, 1, -1 +! if(is_generable(psi_det_generators(1,1,i), tq(1,1,i_alpha), Nint)) then +! ok = .true. +! exit +! end if +! end do +! if(.not. ok) then +! cycle +! end if + if(key_mask(1,1) /= 0) then call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) @@ -138,37 +149,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge do j=1,idx_alpha(0) idx_alpha(j) = idx_microlist_zero(idx_alpha(j)) end do - - -! i = 1 -! j = 2 -! do j = 2, idx_alpha_tmp(0) -! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit -! end do -! -! m = j -! -! idx_alpha(0) = idx_alpha_tmp(0) -! -! do l = 1, idx_alpha(0) -! if(j > idx_alpha_tmp(0)) then -! k = i -! i += 1 -! else if(i >= m) then -! k = j -! j += 1 -! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then -! k = i -! i += 1 -! else -! k = j -! j += 1 -! end if -! ! k=l -! idx_alpha(l) = idx_alpha_tmp(k) -! degree_alpha(l) = degree_alpha_tmp(k) -! end do -! else call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) do j=1,idx_alpha(0) @@ -177,12 +157,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge end if -! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) -! do j=1,idx_alpha(0) -! idx_alpha(j) = idx_miniList(idx_alpha(j)) -! end do - !print *, idx_alpha(:idx_alpha(0)) - do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hij_cache(k_sd)) @@ -285,33 +259,31 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge enddo enddo call omp_set_lock( psi_ref_lock(i_I) ) + + do i_state=1,Nstates if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) - delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) enddo else - delta_ii_(i_state,i_I) = 0.d0 + !delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) enddo endif enddo call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo - !deallocate (dIa_hla,hij_cache) - !deallocate(miniList, idx_miniList) + deallocate (dIa_hla,hij_cache) + deallocate(miniList, idx_miniList) end - - - - subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) use bitmasks @@ -360,7 +332,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq endif enddo if (good) then - if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then N_tq += 1 do k=1,N_int tq(k,1,N_tq) = det_buffer(k,1,i) @@ -437,7 +409,7 @@ subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,N endif enddo if (good) then - if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then N_tq += 1 do k=1,N_int tq(k,1,N_tq) = det_buffer(k,1,i) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index cfaf7a67..8873a940 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -112,20 +112,12 @@ END_PROVIDER lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i endif endif -! j = int(lambda_mrcc(k,i) * 100) -! if(j < -200) j = -200 -! if(j > 200) j = 200 -! histo(j) += 1 enddo enddo lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 end if - -! do i=-200,200 -! print *, i, histo(i) -! end do print*,'N_det_non_ref = ',N_det_non_ref - !print*,'Number of ignored determinants = ',i_pert_count + print*,'Number of ignored determinants = ',i_pert_count print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) print*,'lambda max = ',maxval(dabs(lambda_mrcc)) @@ -163,6 +155,7 @@ END_PROVIDER call H_apply_mrcc(delta_ij,delta_ii,N_states,N_det_non_ref,N_det_ref) END_PROVIDER + BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ] implicit none @@ -288,3 +281,447 @@ subroutine diagonalize_CI_dressed(lambda) SOFT_TOUCH psi_coef end + + +logical function is_generable(det1, det2, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) + integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t + integer, external :: searchExc + logical, external :: excEq + double precision :: phase + + is_generable = .false. + call get_excitation(det1, det2, exc, degree, phase, Nint) + if(degree == -1) return + if(degree == 0) then + is_generable = .true. + return + end if + if(degree > 2) stop "?22??" + !!!!! +! call dec_exc(exc, h1, h2, p1, p2) +! f = searchExc(toutmoun(1,1), (/h1, h2, p1, p2/), hh_shortcut(hh_shortcut(0)+1)-1) +! !print *, toutmoun(:,1), hh_shortcut(hh_shortcut(0)+1)-1, (/h1, h2, p1, p2/) +! if(f /= -1) then +! is_generable = .true. +! if(.not. excEq(toutmoun(1,f), (/h1, h2, p1, p2/))) stop "????" +! end if +! ! print *, f +! return + + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + + if(degree == 1) then + h2 = h1 + p2 = p1 + s2 = s1 + h1 = 0 + p1 = 0 + s1 = 0 + end if + + if(h1 + s1*mo_tot_num < h2 + s2*mo_tot_num) then + f = searchExc(hh_exists(1,1), (/s1, h1, s2, h2/), hh_shortcut(0)) + else + f = searchExc(hh_exists(1,1), (/s2, h2, s1, h1/), hh_shortcut(0)) + end if + if(f == -1) return + + if(p1 + s1*mo_tot_num < p2 + s2*mo_tot_num) then + f = searchExc(pp_exists(1,hh_shortcut(f)), (/s1, p1, s2, p2/), hh_shortcut(f+1)-hh_shortcut(f)) + else + f = searchExc(pp_exists(1,hh_shortcut(f)), (/s2, p2, s1, p1/), hh_shortcut(f+1)-hh_shortcut(f)) + end if + + if(f /= -1) is_generable = .true. +end function + + + +integer function searchDet(dets, det, n, Nint) + implicit none + use bitmasks + + integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2) + integer, intent(in) :: nint, n + integer :: l, h, c + integer, external :: detCmp + logical, external :: detEq + + l = 1 + h = n + do while(.true.) + searchDet = (l+h)/2 + c = detCmp(dets(1,1,searchDet), det(:,:), Nint) + if(c == 0) return + if(c == 1) then + h = searchDet-1 + else + l = searchDet+1 + end if + if(l > h) then + searchDet = -1 + return + end if + + end do +end function + + +integer function searchExc(excs, exc, n) + implicit none + use bitmasks + + integer, intent(in) :: n + integer,intent(in) :: excs(4,n), exc(4) + integer :: l, h, c + integer, external :: excCmp + logical, external :: excEq + + l = 1 + h = n + do + searchExc = (l+h)/2 + c = excCmp(excs(1,searchExc), exc(:)) + if(c == 0) return + if(c == 1) then + h = searchExc-1 + else + l = searchExc+1 + end if + if(l > h) then + searchExc = -1 + return + end if + end do +end function + + +subroutine sort_det(key, idx, N_key, Nint) + implicit none + + + integer, intent(in) :: Nint, N_key + integer(8),intent(inout) :: key(Nint,2,N_key) + integer,intent(out) :: idx(N_key) + integer(8) :: tmp(Nint, 2) + integer :: tmpidx,i,ni + + do i=1,N_key + idx(i) = i + end do + + do i=N_key/2,1,-1 + call tamiser(key, idx, i, N_key, Nint, N_key) + end do + + do i=N_key,2,-1 + do ni=1,Nint + tmp(ni,1) = key(ni,1,i) + tmp(ni,2) = key(ni,2,i) + key(ni,1,i) = key(ni,1,1) + key(ni,2,i) = key(ni,2,1) + key(ni,1,1) = tmp(ni,1) + key(ni,2,1) = tmp(ni,2) + enddo + + tmpidx = idx(i) + idx(i) = idx(1) + idx(1) = tmpidx + call tamiser(key, idx, 1, i-1, Nint, N_key) + end do +end subroutine + + +subroutine sort_exc(key, N_key) + implicit none + + + integer, intent(in) :: N_key + integer,intent(inout) :: key(4,N_key) + integer :: tmp(4) + integer :: i,ni + + + do i=N_key/2,1,-1 + call tamise_exc(key, i, N_key, N_key) + end do + + do i=N_key,2,-1 + do ni=1,4 + tmp(ni) = key(ni,i) + key(ni,i) = key(ni,1) + key(ni,1) = tmp(ni) + enddo + + call tamise_exc(key, 1, i-1, N_key) + end do +end subroutine + + +logical function exc_inf(exc1, exc2) + implicit none + integer,intent(in) :: exc1(4), exc2(4) + integer :: i + exc_inf = .false. + do i=1,4 + if(exc1(i) < exc2(i)) then + exc_inf = .true. + return + else if(exc1(i) > exc2(i)) then + return + end if + end do +end function + + +subroutine tamise_exc(key, no, n, N_key) + use bitmasks + implicit none + + BEGIN_DOC +! Uncodumented : TODO + END_DOC + integer,intent(in) :: no, n, N_key + integer,intent(inout) :: key(4, N_key) + integer :: k,j + integer :: tmp(4) + logical :: exc_inf + integer :: ni + + k = no + j = 2*k + do while(j <= n) + if(j < n) then + if (exc_inf(key(1,j), key(1,j+1))) then + j = j+1 + endif + endif + if(exc_inf(key(1,k), key(1,j))) then + do ni=1,4 + tmp(ni) = key(ni,k) + key(ni,k) = key(ni,j) + key(ni,j) = tmp(ni) + enddo + k = j + j = k+k + else + return + endif + enddo +end subroutine + + +subroutine dec_exc(exc, h1, h2, p1, p2) + implicit none + integer :: exc(0:2,2,2), s1, s2, degree + integer, intent(out) :: h1, h2, p1, p2 + + degree = exc(0,1,1) + exc(0,1,2) + + h1 = 0 + h2 = 0 + p1 = 0 + p2 = 0 + + if(degree == 0) return + + call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2) + + h1 += mo_tot_num * (s1-1) + p1 += mo_tot_num * (s1-1) + + if(degree == 2) then + h2 += mo_tot_num * (s2-1) + p2 += mo_tot_num * (s2-1) + if(h1 > h2) then + s1 = h1 + h1 = h2 + h2 = s1 + end if + if(p1 > p2) then + s1 = p1 + p1 = p2 + p2 = s1 + end if + else + h2 = h1 + p2 = p1 + p1 = 0 + h1 = 0 + end if +end subroutine + + + BEGIN_PROVIDER [ integer, hh_exists, (4, N_det_ref * N_det_non_ref) ] +&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_det_ref * N_det_non_ref + 1) ] +&BEGIN_PROVIDER [ integer, pp_exists, (4, N_det_ref * N_det_non_ref) ] + implicit none + integer,allocatable :: num(:,:) + integer :: exc(0:2, 2, 2), degree, n, on, s, h1, h2, p1, p2, l, i + double precision :: phase + logical, external :: excEq + + allocate(num(4, N_det_ref * N_det_non_ref)) + + hh_shortcut = 0 + hh_exists = 0 + pp_exists = 0 + num = 0 + + n = 0 + do i=1, N_det_ref + do l=1, N_det_non_ref + call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int) + if(degree == -1) cycle + call dec_exc(exc, h1, h2, p1, p2) + n += 1 + num(:, n) = (/h1, h2, p1, p2/) + end do + end do + + call sort_exc(num, n) + + hh_shortcut(0) = 1 + hh_shortcut(1) = 1 + hh_exists(:,1) = (/1, num(1,1), 1, num(2,1)/) + pp_exists(:,1) = (/1, num(3,1), 1, num(4,1)/) + s = 1 + do i=2,n + if(.not. excEq(num(1,i), num(1,s))) then + s += 1 + num(:, s) = num(:, i) + pp_exists(:,s) = (/1, num(3,s), 1, num(4,s)/) + if(hh_exists(2, hh_shortcut(0)) /= num(1,s) .or. & + hh_exists(4, hh_shortcut(0)) /= num(2,s)) then + hh_shortcut(0) += 1 + hh_shortcut(hh_shortcut(0)) = s + hh_exists(:,hh_shortcut(0)) = (/1, num(1,s), 1, num(2,s)/) + end if + end if + end do + hh_shortcut(hh_shortcut(0)+1) = s+1 + + do s=2,4,2 + do i=1,hh_shortcut(0) + if(hh_exists(s, i) == 0) then + hh_exists(s-1, i) = 0 + else if(hh_exists(s, i) > mo_tot_num) then + hh_exists(s, i) -= mo_tot_num + hh_exists(s-1, i) = 2 + end if + end do + + do i=1,hh_shortcut(hh_shortcut(0)+1)-1 + if(pp_exists(s, i) == 0) then + pp_exists(s-1, i) = 0 + else if(pp_exists(s, i) > mo_tot_num) then + pp_exists(s, i) -= mo_tot_num + pp_exists(s-1, i) = 2 + end if + end do + end do +END_PROVIDER + + +logical function excEq(exc1, exc2) + implicit none + integer, intent(in) :: exc1(4), exc2(4) + integer :: i + excEq = .false. + do i=1, 4 + if(exc1(i) /= exc2(i)) return + end do + excEq = .true. +end function + + +integer function excCmp(exc1, exc2) + implicit none + integer, intent(in) :: exc1(4), exc2(4) + integer :: i + excCmp = 0 + do i=1, 4 + if(exc1(i) > exc2(i)) then + excCmp = 1 + return + else if(exc1(i) < exc2(i)) then + excCmp = -1 + return + end if + end do +end function + + +subroutine apply_hole(det, exc, res, ok, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: exc(4) + integer :: s1, s2, h1, h2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos + + ok = .false. + s1 = exc(1) + h1 = exc(2) + s2 = exc(3) + h2 = exc(4) + res = det + + if(h1 /= 0) then + ii = (h1-1)/bit_kind_size + 1 + pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s1) = ibclr(res(ii, s1), pos) + end if + + ii = (h2-1)/bit_kind_size + 1 + pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return + res(ii, s2) = ibclr(res(ii, s2), pos) + + + ok = .true. +end subroutine + + +subroutine apply_particle(det, exc, res, ok, Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: exc(4) + integer :: s1, s2, p1, p2 + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: ii, pos + + ok = .false. + s1 = exc(1) + p1 = exc(2) + s2 = exc(3) + p2 = exc(4) + res = det + + if(p1 /= 0) then + ii = (p1-1)/bit_kind_size + 1 + pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) + if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s1) = ibset(res(ii, s1), pos) + end if + + ii = (p2-1)/bit_kind_size + 1 + pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) + if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return + res(ii, s2) = ibset(res(ii, s2), pos) + + + ok = .true. +end subroutine + diff --git a/plugins/mrcepa0/dressing.irp.f b/plugins/mrcepa0/dressing.irp.f index e0689642..5d8d4850 100644 --- a/plugins/mrcepa0/dressing.irp.f +++ b/plugins/mrcepa0/dressing.irp.f @@ -1,227 +1,32 @@ use bitmasks -subroutine dec_exc(exc, h1, h2, p1, p2) - implicit none - integer :: exc(0:2,2,2), s1, s2, degree - integer, intent(out) :: h1, h2, p1, p2 - - degree = exc(0,1,1) + exc(0,1,2) - - h1 = 0 - h2 = 0 - p1 = 0 - p2 = 0 - - if(degree == 0) return - - call decode_exc(exc, degree, h1, p1, h2, p2, s1, s2) - - h1 += mo_tot_num * (s1-1) - p1 += mo_tot_num * (s1-1) - - if(degree == 2) then - h2 += mo_tot_num * (s2-1) - p2 += mo_tot_num * (s2-1) - if(h1 > h2) then - s1 = h1 - h1 = h2 - h2 = s1 - end if - if(p1 > p2) then - s1 = p1 - p1 = p2 - p2 = s1 - end if - else - h2 = h1 - p2 = p1 - p1 = 0 - h1 = 0 - end if -end subroutine - - - - BEGIN_PROVIDER [ integer, hh_exists, (4, N_det_ref * N_det_non_ref) ] -&BEGIN_PROVIDER [ integer, hh_shortcut, (0:N_det_ref * N_det_non_ref + 1) ] -&BEGIN_PROVIDER [ integer, pp_exists, (4, N_det_ref * N_det_non_ref) ] - implicit none - integer :: num(0:mo_tot_num*2, 0:mo_tot_num*2) - integer :: exc(0:2, 2, 2), degree, n, on, s, h1, h2, p1, p2, l, i - double precision :: phase - - hh_shortcut = 0 - hh_exists = 0 - pp_exists = 0 - num = 0 - - do i=1, N_det_ref - do l=1, N_det_non_ref - call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int) - if(degree == -1) cycle - call dec_exc(exc, h1, h2, p1, p2) - num(h1, h2) += 1 - end do - end do - - n = 1 - do l=0,mo_tot_num*2 - do i=0,l - on = num(i,l) - if(on /= 0) then - hh_shortcut(0) += 1 - hh_shortcut(hh_shortcut(0)) = n - hh_exists(:, hh_shortcut(0)) = (/1, i, 1, l/) - end if - - num(i,l) = n - n += on - end do - end do - - hh_shortcut(hh_shortcut(0)+1) = n - - do i=1, N_det_ref - do l=1, N_det_non_ref - call get_excitation(psi_ref(1,1,i), psi_non_ref(1,1,l), exc, degree, phase, N_int) - if(degree == -1) cycle - call dec_exc(exc, h1, h2, p1, p2) - pp_exists(:, num(h1, h2)) = (/1,p1,1,p2/) - num(h1, h2) += 1 - end do - end do - - do s=2,4,2 - do i=1,hh_shortcut(0) - if(hh_exists(s, i) == 0) then - hh_exists(s-1, i) = 0 - else if(hh_exists(s, i) > mo_tot_num) then - hh_exists(s, i) -= mo_tot_num - hh_exists(s-1, i) = 2 - end if - end do - - do i=1,hh_shortcut(hh_shortcut(0)+1)-1 - if(pp_exists(s, i) == 0) then - pp_exists(s-1, i) = 0 - else if(pp_exists(s, i) > mo_tot_num) then - pp_exists(s, i) -= mo_tot_num - pp_exists(s-1, i) = 2 - end if - end do - end do - -END_PROVIDER - - - - -subroutine apply_hole(det, exc, res, ok, Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer, intent(in) :: exc(4) - integer :: s1, s2, h1, h2 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos - - ok = .false. - s1 = exc(1) - h1 = exc(2) - s2 = exc(3) - h2 = exc(4) - res = det - - if(h1 /= 0) then - ii = (h1-1)/bit_kind_size + 1 - pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64 - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s1) = ibclr(res(ii, s1), pos) - end if - - ii = (h2-1)/bit_kind_size + 1 - pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return - res(ii, s2) = ibclr(res(ii, s2), pos) - - - ok = .true. -end subroutine - - -subroutine apply_particle(det, exc, res, ok, Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer, intent(in) :: exc(4) - integer :: s1, s2, p1, p2 - integer(bit_kind),intent(in) :: det(Nint, 2) - integer(bit_kind),intent(out) :: res(Nint, 2) - logical, intent(out) :: ok - integer :: ii, pos - - ok = .false. - s1 = exc(1) - p1 = exc(2) - s2 = exc(3) - p2 = exc(4) - res = det - - if(p1 /= 0) then - ii = (p1-1)/bit_kind_size + 1 - pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1) - if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s1) = ibset(res(ii, s1), pos) - end if - - ii = (p2-1)/bit_kind_size + 1 - pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1) - if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return - res(ii, s2) = ibset(res(ii, s2), pos) - - - ok = .true. -end subroutine - BEGIN_PROVIDER [ double precision, delta_ij_mrcc, (N_states,N_det_non_ref,N_det_ref) ] &BEGIN_PROVIDER [ double precision, delta_ii_mrcc, (N_states, N_det_ref) ] use bitmasks implicit none - integer :: gen, h, p, i_state, n, t + integer :: gen, h, p, i_state, n, t, i, h1, h2, p1, p2, s1, s2 integer(bit_kind) :: mask(N_int, 2), omask(N_int, 2), buf(N_int, 2, N_det_non_ref) logical :: ok + logical, external :: detEq delta_ij_mrcc = 0d0 delta_ii_mrcc = 0d0 i_state = 1 - do gen=1, N_det_generators - !print *, gen, "/", N_det_generators + do gen= 1, N_det_generators + print *, gen, "/", N_det_generators do h=1, hh_shortcut(0) call apply_hole(psi_det_generators(1,1,gen), hh_exists(1, h), mask, ok, N_int) if(.not. ok) cycle - omask = 0 + omask = 0_bit_kind if(hh_exists(1, h) /= 0) omask = mask - !-459.6378590456251 - !-199.0659502581943 n = 1 - ploop : do p=hh_shortcut(h), hh_shortcut(h+1)-1 - - do t=hh_shortcut(h), p-1 - if(pp_exists(1, p) == pp_exists(1,t) .and. & - pp_exists(2, p) == pp_exists(2,t) .and. & - pp_exists(3, p) == pp_exists(3,t) .and. & - pp_exists(4, p) == pp_exists(4,t)) cycle ploop - end do + do p=hh_shortcut(h), hh_shortcut(h+1)-1 call apply_particle(mask, pp_exists(1, p), buf(1,1,n), ok, N_int) - !-459.6379081607463 - !-199.0659982685706 if(ok) n = n + 1 - end do ploop + end do n = n - 1 if(n /= 0) call mrcc_part_dress(delta_ij_mrcc, delta_ii_mrcc,gen,n,buf,N_int,omask) end do @@ -229,7 +34,6 @@ end subroutine END_PROVIDER - subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffer,Nint,key_mask) use bitmasks implicit none @@ -258,7 +62,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe integer :: i_state, k_sd, l_sd, i_I, i_alpha integer(bit_kind),allocatable :: miniList(:,:,:) - integer(bit_kind),intent(in) :: key_mask(Nint, 2) + integer(bit_kind) :: key_mask(Nint, 2) integer,allocatable :: idx_miniList(:) integer :: N_miniList, ni, leng double precision, allocatable :: hij_cache(:) @@ -266,18 +70,18 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:) integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:) integer :: mobiles(2), smallerlist - - - + logical, external :: detEq, is_generable + + leng = max(N_det_generators, N_det_non_ref) allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref)) !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint) call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint) - if(fullMatch) then - return - end if +! if(fullMatch) then +! return +! end if allocate(ptr_microlist(0:mo_tot_num*2+1), & N_microlist(0:mo_tot_num*2) ) @@ -286,9 +90,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe if(key_mask(1,1) /= 0) then call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint) - call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + call filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) else - call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) + call filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist) end if @@ -332,7 +136,7 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe do i_alpha=1,N_tq - if(key_mask(1,1) /= 0) then + if(key_mask(1,1) /= 0) then call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then @@ -463,26 +267,27 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,i_generator,n_selected,det_buffe enddo enddo call omp_set_lock( psi_ref_lock(i_I) ) + do i_state=1,N_states if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) - delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd) enddo else delta_ii_(i_state,i_I) = 0.d0 do l_sd=1,idx_alpha(0) k_sd = idx_alpha(l_sd) - delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) + delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd) enddo endif enddo call omp_unset_lock( psi_ref_lock(i_I) ) enddo enddo - !deallocate (dIa_hla,hij_cache) - !deallocate(miniList, idx_miniList) + deallocate (dIa_hla,hij_cache) + deallocate(miniList, idx_miniList) end @@ -494,29 +299,30 @@ end implicit none integer :: i, j, i_state - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc do i_state = 1, N_states if(mrmode == 3) then - do i = 1, N_det_ref - delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) - do j = 1, N_det_non_ref - delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) - end do - end do - else if(mrmode == 2) then - do i = 1, N_det_ref - delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) - do j = 1, N_det_non_ref - delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) - end do + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_mrcc(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_mrcc(i_state,j,i) end do + end do +! ! do i = 1, N_det_ref -! delta_ii(i_state,i)= delta_ii_old(i_state,i) +! delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) - delta_sub_ii(i,i_state) ! do j = 1, N_det_non_ref -! delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) +! delta_ij(i_state,j,i) = delta_mrcepa0_ij(i,j,i_state) - delta_sub_ij(i,j,i_state) ! end do ! end do + else if(mrmode == 2) then + do i = 1, N_det_ref + delta_ii(i_state,i)= delta_ii_old(i_state,i) + do j = 1, N_det_non_ref + delta_ij(i_state,j,i) = delta_ij_old(i_state,j,i) + end do + end do else if(mrmode == 1) then do i = 1, N_det_ref delta_ii(i_state,i)= delta_mrcepa0_ii(i,i_state) @@ -646,7 +452,6 @@ END_PROVIDER end do end do end do - print *, "pre done" END_PROVIDER @@ -708,21 +513,6 @@ END_PROVIDER END_PROVIDER -logical function detEq(a,b,Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) - integer :: ni, i - - detEq = .false. - do i=1,2 - do ni=1,Nint - if(a(ni,i) /= b(ni,i)) return - end do - end do - detEq = .true. -end function logical function isInCassd(a,Nint) @@ -793,106 +583,6 @@ subroutine getHP(a,h,p,Nint) !isInCassd = .true. end function -integer function detCmp(a,b,Nint) - use bitmasks - implicit none - integer, intent(in) :: Nint - integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) - integer :: ni, i - - detCmp = 0 - do i=1,2 - do ni=Nint,1,-1 - - if(a(ni,i) < b(ni,i)) then - detCmp = -1 - return - else if(a(ni,i) > b(ni,i)) then - detCmp = 1 - return - end if - - end do - end do -end function - - -integer function searchDet(dets, det, n, Nint) - implicit none - use bitmasks - - integer(bit_kind),intent(in) :: dets(Nint,2,n), det(Nint,2) - integer, intent(in) :: nint, n - integer :: l, h, c - integer, external :: detCmp - logical, external :: detEq - - !do l=1,n - ! if(detEq(det(1,1), dets(1,1,l),Nint)) then - ! searchDet = l - ! return - ! end if - !end do - !searchDet = -1 - !return - - - l = 1 - h = n - do while(.true.) - searchDet = (l+h)/2 - c = detCmp(dets(1,1,searchDet), det(:,:), Nint) - if(c == 0) return - if(c == 1) then - h = searchDet-1 - else - l = searchDet+1 - end if - if(l > h) then - searchDet = -1 - return - end if - - end do -end function - - -subroutine sort_det(key, idx, N_key, Nint) - implicit none - - - integer, intent(in) :: Nint, N_key - integer(8),intent(inout) :: key(Nint,2,N_key) - integer,intent(out) :: idx(N_key) - integer(8) :: tmp(Nint, 2) - integer :: tmpidx,i,ni - - do i=1,N_key - idx(i) = i - end do - - do i=N_key/2,1,-1 - call tamiser(key, idx, i, N_key, Nint, N_key) - end do - - do i=N_key,2,-1 - do ni=1,Nint - tmp(ni,1) = key(ni,1,i) - tmp(ni,2) = key(ni,2,i) - key(ni,1,i) = key(ni,1,1) - key(ni,2,i) = key(ni,2,1) - key(ni,1,1) = tmp(ni,1) - key(ni,2,1) = tmp(ni,2) - enddo - - tmpidx = idx(i) - idx(i) = idx(1) - idx(1) = tmpidx - call tamiser(key, idx, 1, i-1, Nint, N_key) - end do -end subroutine - - BEGIN_PROVIDER [ double precision, delta_mrcepa0_ij, (N_det_ref,N_det_non_ref,N_states) ] &BEGIN_PROVIDER [ double precision, delta_mrcepa0_ii, (N_det_ref,N_states) ] @@ -1116,6 +806,7 @@ end subroutine BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ] + implicit none integer :: i,j do i=1,N_det_ref do j=1,N_det_non_ref @@ -1125,3 +816,135 @@ BEGIN_PROVIDER [ double precision, h_, (N_det_ref,N_det_non_ref) ] END_PROVIDER + +subroutine filter_tq(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + integer(bit_kind),intent(in) :: miniList(Nint,2,N_det_generators) + integer,intent(in) :: N_miniList + + + N_tq = 0 + + i_loop : do i=1,N_selected + do k=1, N_minilist + if(is_generable(miniList(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + + ! Select determinants that are triple or quadruple excitations + ! from the ref + good = .True. + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo i_loop +end + + +subroutine filter_tq_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask) + + use bitmasks + implicit none + + integer, intent(in) :: i_generator,n_selected, Nint + + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k,m + logical :: is_in_wavefunction + integer :: degree(psi_det_size) + integer :: idx(0:psi_det_size) + logical :: good + + integer(bit_kind), intent(out) :: tq(Nint,2,n_selected) + integer, intent(out) :: N_tq + + integer :: nt,ni + logical, external :: is_connected_to, is_generable + + integer(bit_kind),intent(in) :: microlist(Nint,2,*) + integer,intent(in) :: ptr_microlist(0:*) + integer,intent(in) :: N_microlist(0:*) + integer(bit_kind),intent(in) :: key_mask(Nint, 2) + + integer :: mobiles(2), smallerlist + + + N_tq = 0 + + i_loop : do i=1,N_selected + call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint) + if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then + smallerlist = mobiles(1) + else + smallerlist = mobiles(2) + end if + + if(N_microlist(smallerlist) > 0) then + do k=ptr_microlist(smallerlist), ptr_microlist(smallerlist)+N_microlist(smallerlist)-1 + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + end if + + if(N_microlist(0) > 0) then + do k=1, N_microlist(0) + if(is_generable(microlist(1,1,k), det_buffer(1,1,i), Nint)) cycle i_loop + end do + end if + + ! Select determinants that are triple or quadruple excitations + ! from the ref + good = .True. + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector + do k=1,idx(0) + if (degree(k) < 3) then + good = .False. + exit + endif + enddo + if (good) then + if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint)) then + N_tq += 1 + do k=1,N_int + tq(k,1,N_tq) = det_buffer(k,1,i) + tq(k,2,N_tq) = det_buffer(k,2,i) + enddo + endif + endif + enddo i_loop +end + + + + diff --git a/plugins/mrcepa0/mrsc2sub.irp.f b/plugins/mrcepa0/mrcc.irp.f similarity index 88% rename from plugins/mrcepa0/mrsc2sub.irp.f rename to plugins/mrcepa0/mrcc.irp.f index 07a07c83..91592e62 100644 --- a/plugins/mrcepa0/mrsc2sub.irp.f +++ b/plugins/mrcepa0/mrcc.irp.f @@ -3,7 +3,7 @@ program mrsc2sub double precision, allocatable :: energy(:) allocate (energy(N_states)) - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc mrmode = 3 read_wf = .True. diff --git a/plugins/mrcepa0/mrcepa0.irp.f b/plugins/mrcepa0/mrcepa0.irp.f index 9473361b..34d3dec5 100644 --- a/plugins/mrcepa0/mrcepa0.irp.f +++ b/plugins/mrcepa0/mrcepa0.irp.f @@ -3,7 +3,7 @@ program mrcepa0 double precision, allocatable :: energy(:) allocate (energy(N_states)) - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc mrmode = 1 read_wf = .True. diff --git a/plugins/mrcepa0/mrcepa0_general.irp.f b/plugins/mrcepa0/mrcepa0_general.irp.f index df10de34..82b1fc9b 100644 --- a/plugins/mrcepa0/mrcepa0_general.irp.f +++ b/plugins/mrcepa0/mrcepa0_general.irp.f @@ -15,6 +15,7 @@ subroutine run(N_st,energy) integer :: n_it_mrcc_max double precision :: thresh_mrcc + thresh_mrcc = 1d-7 n_it_mrcc_max = 10 diff --git a/plugins/mrcepa0/mrsc2.irp.f b/plugins/mrcepa0/mrsc2.irp.f index d4e1b1d4..d0f44a33 100644 --- a/plugins/mrcepa0/mrsc2.irp.f +++ b/plugins/mrcepa0/mrsc2.irp.f @@ -3,9 +3,8 @@ program mrsc2 double precision, allocatable :: energy(:) allocate (energy(N_states)) - !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrsc2 sub + !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc mrmode = 2 - read_wf = .True. SOFT_TOUCH read_wf call print_cas_coefs diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 4e419af5..c46a5bb0 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -11,7 +11,7 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl integer(bit_kind), intent(in) :: key_prev(N_int, 2, *) PROVIDE N_int PROVIDE N_det - + $declarations @@ -184,7 +184,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl $initialization - + $omp_parallel !$ iproc = omp_get_thread_num() allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 7a54bdbc..2f53c799 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -165,7 +165,7 @@ logical function is_connected_to(key,keys,Nint,Ndet) integer :: i, l integer :: degree_x2 - + logical, external :: is_generable_cassd ASSERT (Nint > 0) ASSERT (Nint == N_int) @@ -183,12 +183,35 @@ logical function is_connected_to(key,keys,Nint,Ndet) if (degree_x2 > 4) then cycle else +! if(.not. is_generable_cassd(keys(1,1,i), key(1,1), Nint)) cycle !!!Nint==1 !!!!! is_connected_to = .true. return endif enddo end + +logical function is_generable_cassd(det1, det2, Nint) !!! TEST Cl HARD !!!!! + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind) :: det1(Nint, 2), det2(Nint, 2) + integer :: degree, f, exc(0:2, 2, 2), h1, h2, p1, p2, s1, s2, t + double precision :: phase + + is_generable_cassd = .false. + call get_excitation(det1, det2, exc, degree, phase, Nint) + if(degree == -1) return + if(degree == 0) then + is_generable_cassd = .true. + return + end if + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + if(degree == 1 .and. h1 <= 11) is_generable_cassd = .true. + if(degree == 2 .and. h1 <= 11 .and. h2 <= 11) is_generable_cassd = .true. +end function + + logical function is_connected_to_by_mono(key,keys,Nint,Ndet) use bitmasks implicit none diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 4476ed45..400345c1 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -664,3 +664,44 @@ subroutine save_wavefunction_specified(ndet,nstates,psidet,psicoef,ndetsave,inde end +logical function detEq(a,b,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detEq = .false. + do i=1,2 + do ni=1,Nint + if(a(ni,i) /= b(ni,i)) return + end do + end do + detEq = .true. +end function + + +integer function detCmp(a,b,Nint) + use bitmasks + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2) + integer :: ni, i + + detCmp = 0 + do i=1,2 + do ni=Nint,1,-1 + + if(a(ni,i) < b(ni,i)) then + detCmp = -1 + return + else if(a(ni,i) > b(ni,i)) then + detCmp = 1 + return + end if + + end do + end do +end function + + diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 0b456751..2b5ae4f1 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -914,7 +914,6 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis fullMatch = .false. N_miniList = 0 N_subList = 0 - l = popcnt(key_mask(1,1)) + popcnt(key_mask(1,2)) do ni = 2,Nint l = l + popcnt(key_mask(ni,1)) + popcnt(key_mask(ni,2)) @@ -947,8 +946,13 @@ subroutine create_minilist_find_previous(key_mask, fullList, miniList, N_fullLis miniList(ni,2,N_minilist) = fullList(ni,2,i) enddo else if(k == 0) then - fullMatch = .true. - return + N_minilist += 1 + do ni=1,Nint + miniList(ni,1,N_minilist) = fullList(ni,1,i) + miniList(ni,2,N_minilist) = fullList(ni,2,i) + enddo +! fullMatch = .true. +! return end if end do end if diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f index 8d5726f5..2eec0dea 100644 --- a/src/Determinants/spindeterminants.irp.f +++ b/src/Determinants/spindeterminants.irp.f @@ -10,7 +10,7 @@ integer*8 function spin_det_search_key(det,Nint) use bitmasks implicit none BEGIN_DOC -! Return an integer*8 corresponding to a determinant index for searching +! Return an integer(8) corresponding to a determinant index for searching END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: det(Nint) @@ -64,9 +64,9 @@ BEGIN_TEMPLATE integer :: i,j,k integer, allocatable :: iorder(:) - integer*8, allocatable :: bit_tmp(:) - integer*8 :: last_key - integer*8, external :: spin_det_search_key + integer(8), allocatable :: bit_tmp(:) + integer(8) :: last_key + integer(8), external :: spin_det_search_key logical,allocatable :: duplicate(:) allocate ( iorder(N_det), bit_tmp(N_det), duplicate(N_det) ) @@ -149,8 +149,8 @@ integer function get_index_in_psi_det_alpha_unique(key,Nint) integer(bit_kind), intent(in) :: key(Nint) integer :: i, ibegin, iend, istep, l - integer*8 :: det_ref, det_search - integer*8, external :: spin_det_search_key + integer(8) :: det_ref, det_search + integer(8), external :: spin_det_search_key logical :: in_wavefunction in_wavefunction = .False. @@ -231,8 +231,8 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) integer(bit_kind), intent(in) :: key(Nint) integer :: i, ibegin, iend, istep, l - integer*8 :: det_ref, det_search - integer*8, external :: spin_det_search_key + integer(8) :: det_ref, det_search + integer(8), external :: spin_det_search_key logical :: in_wavefunction in_wavefunction = .False. @@ -305,10 +305,10 @@ end subroutine write_spindeterminants use bitmasks implicit none - integer*8, allocatable :: tmpdet(:,:) + integer(8), allocatable :: tmpdet(:,:) integer :: N_int2 integer :: i,j,k - integer*8 :: det_8(100) + integer(8) :: det_8(100) integer(bit_kind) :: det_bk((100*8)/bit_kind) equivalence (det_8, det_bk)