From 1469822e944d6f44620b726874775c7b9fdd65ad Mon Sep 17 00:00:00 2001 From: Yann Garniron Date: Thu, 31 Mar 2016 10:40:39 +0200 Subject: [PATCH] corrected bug in MRCC --- plugins/MRCC_Utils/mrcc_dress.irp.f | 102 ++++++++++++---------------- plugins/MRCC_Utils/mrcc_utils.irp.f | 90 +++++++++++++++++------- src/Determinants/slater_rules.irp.f | 52 ++++++++++++++ 3 files changed, 162 insertions(+), 82 deletions(-) diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 78533348..b2304818 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -25,8 +25,8 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) integer :: i,j,k,l,m - integer :: degree_alpha(psi_det_size), degree_alpha_tmp(psi_det_size) - integer :: idx_alpha(0:psi_det_size), idx_alpha_tmp(0:psi_det_size) + integer :: degree_alpha(psi_det_size) + integer :: idx_alpha(0:psi_det_size) logical :: good, fullMatch integer(bit_kind) :: tq(Nint,2,n_selected) @@ -74,11 +74,6 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge call find_triples_and_quadruples_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) -! microlist(:N_miniList) = minilist(:N_miniList) -! idx_microlist(:N_minilist) = idx_minilist(:) -! N_microlist(0) = N_miniList -! ptr_microlist(1) = 1 -! ptr_microlist(2) = N_miniList + 1 end if @@ -139,41 +134,41 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l) end do - call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha_tmp,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha_tmp) - do j=1,idx_alpha_tmp(0) - idx_alpha_tmp(j) = idx_microlist_zero(idx_alpha_tmp(j)) + call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha) + 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 - - idx_alpha(l) = idx_alpha_tmp(k) - degree_alpha(l) = degree_alpha_tmp(k) - 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) @@ -240,32 +235,17 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge tmp_det(k,1) = psi_ref(k,1,i_I) tmp_det(k,2) = psi_ref(k,2,i_I) enddo - ! Hole (see list_to_bitstring) - iint = ishft(h1-1,-bit_kind_shift) + 1 - ipos = h1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos) - ! Particle - iint = ishft(p1-1,-bit_kind_shift) + 1 - ipos = p1-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos) - if (degree_alpha(k_sd) == 2) then - ! Hole (see list_to_bitstring) - iint = ishft(h2-1,-bit_kind_shift) + 1 - ipos = h2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos) - - ! Particle - iint = ishft(p2-1,-bit_kind_shift) + 1 - ipos = p2-ishft((iint-1),bit_kind_shift)-1 - tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos) - endif + logical :: ok + call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint) + if(.not. ok) cycle ! do i_state=1,Nstates dka(i_state) = 0.d0 enddo do l_sd=k_sd+1,idx_alpha(0) + call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint) if (degree == 0) then @@ -279,11 +259,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge if (.not.loop) then call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint) hIl = hij_mrcc(idx_alpha(l_sd),i_I) - ! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) +! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl) do i_state=1,Nstates dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2 enddo endif + exit endif enddo @@ -324,6 +305,9 @@ end + + + subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList) use bitmasks diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index c41a3431..712daca1 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,37 +1,81 @@ +! BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] +! implicit none +! BEGIN_DOC +! ! cm/ or perturbative 1/Delta_E(m) +! END_DOC +! integer :: i,k +! double precision :: ihpsi(N_states),ihpsi_current(N_states) +! integer :: i_pert_count +! +! i_pert_count = 0 +! lambda_mrcc = 0.d0 +! +! do i=1,N_det_non_ref +! call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) +! do k=1,N_states +! if (ihpsi_current(k) == 0.d0) then +! ihpsi_current(k) = 1.d-32 +! endif +! if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-6) then +! i_pert_count +=1 +! else +! lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) +! endif +! enddo +! enddo +! +! print*,'N_det_non_ref = ',N_det_non_ref +! print*,'Number of ignored determinants = ',i_pert_count +! print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) +! END_PROVIDER + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] implicit none BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) END_DOC integer :: i,k - double precision :: ihpsi(N_states),ihpsi_current(N_states) - integer :: i_pert_count - - i_pert_count = 0 - lambda_mrcc = 0.d0 - - do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current) - do k=1,N_states - if (ihpsi_current(k) == 0.d0) then - ihpsi_current(k) = 1.d-32 + double precision :: ihpsi_current(N_states) + integer :: i_pert_count + double precision :: hii, lambda_pert + + i_pert_count = 0 + lambda_mrcc = 0.d0 + + do i=1,N_det_non_ref + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, & + size(psi_ref_coef,1), N_states,ihpsi_current) + call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) + do k=1,N_states + if (ihpsi_current(k) == 0.d0) then + ihpsi_current(k) = 1.d-32 + endif + lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) + if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-6 ) then + i_pert_count += 1 + lambda_mrcc(k,i) = 0.d0 +! lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) +! if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then +! lambda_mrcc(k,i) = 0.d0 +! endif endif - if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then - i_pert_count +=1 - else - lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) - endif - enddo - enddo + double precision, parameter :: x = 2.d0 + if (lambda_mrcc(k,i) > x) then + lambda_mrcc(k,i) = x + else if (lambda_mrcc(k,i) < -x) then + lambda_mrcc(k,i) = -x + endif + enddo + enddo + + print*,'N_det_non_ref = ',N_det_non_ref + 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 min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc)) - print*,'N_det_non_ref = ',N_det_non_ref - print*,'Number of ignored determinants = ',i_pert_count - print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1) END_PROVIDER - - !BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ] !implicit none !BEGIN_DOC diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index b63ae69f..4ab9deda 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1728,3 +1728,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) deallocate (shortcut, sort_idx, sorted, version) end + +subroutine apply_excitation(det, exc, res, ok, Nint) + use bitmasks + implicit none + + integer, intent(in) :: Nint + integer, intent(in) :: exc(0:2,2,2) + integer(bit_kind),intent(in) :: det(Nint, 2) + integer(bit_kind),intent(out) :: res(Nint, 2) + logical, intent(out) :: ok + integer :: h1,p1,h2,p2,s1,s2,degree + integer :: ii, pos + + + ok = .false. + degree = exc(0,1,1) + exc(0,1,2) + + if(.not. (degree > 0 .and. degree <= 2)) then + print *, degree + print *, "apply ex" + STOP + endif + + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + res = det + + 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) + + 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) + + if(degree == 2) then + 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) + + 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) + endif + + ok = .true. +end subroutine + +