diff --git a/plugins/dress_zmq/alpha_factory.irp.f b/plugins/dress_zmq/alpha_factory.irp.f index 8cb0dcac..f1a3a8e9 100644 --- a/plugins/dress_zmq/alpha_factory.irp.f +++ b/plugins/dress_zmq/alpha_factory.irp.f @@ -20,6 +20,23 @@ subroutine alpha_callback(delta_ij_loc, i_generator, subset,iproc) end subroutine + BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ] +&BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ] + implicit none + integer :: i,inpsisor + + idx_non_ref_from_sorted = 0 + psi_from_sorted = 0 + + do i=1,N_det + psi_from_sorted(psi_det_sorted_order(i)) = i + inpsisor = psi_det_sorted_order(i) + if(inpsisor <= 0) stop "idx_non_ref_from_sorted" + idx_non_ref_from_sorted(inpsisor) = idx_non_ref_rev(i) + end do +END_PROVIDER + + subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index, subset, iproc) use bitmasks implicit none @@ -371,17 +388,17 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index end subroutine -subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexes, indexes_end, abuf, siz, iproc) +subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexes, indexes_end, rabuf, siz, iproc) use bitmasks implicit none double precision,intent(inout) :: delta_ij_loc(N_states,N_det,2) integer, intent(in) :: sp, indexes(0:mo_tot_num, 0:mo_tot_num), siz, iproc - integer, intent(in) :: indexes_end(0:mo_tot_num, 0:mo_tot_num), abuf(*) + integer, intent(in) :: indexes_end(0:mo_tot_num, 0:mo_tot_num), rabuf(*) logical, intent(in) :: bannedOrb(mo_tot_num,2), banned(mo_tot_num, mo_tot_num) integer(bit_kind), intent(in) :: mask(N_int, 2) integer(bit_kind) :: alpha(N_int, 2) - integer, allocatable :: labuf(:) + integer, allocatable :: labuf(:), abuf(:) logical :: ok integer :: i,j,k,s,st1,st2,st3,st4 integer :: lindex(mo_tot_num,2), lindex_end(mo_tot_num, 2) @@ -390,14 +407,19 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe integer(bit_kind), allocatable :: det_minilist(:,:,:) - allocate(labuf(N_det), putten(N_det), det_minilist(N_int, 2, N_det)) + allocate(abuf(siz), labuf(N_det), putten(N_det), det_minilist(N_int, 2, N_det)) + + do i=1,siz + abuf(i) = psi_from_sorted(rabuf(i)) + end do + putten = .false. st1 = indexes_end(0,0)-1 !! if(st1 > 0) then labuf(:st1) = abuf(:st1) do i=1,st1 - det_minilist(:,:,i) = psi_det_sorted(:,:,labuf(i)) + det_minilist(:,:,i) = psi_det(:,:,labuf(i)) end do end if st1 += 1 @@ -421,66 +443,66 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe lindex_end(:,1) = indexes_end(1:, 0)-1 end if - do i=1,mo_tot_num - if(bannedOrb(i,s1)) cycle - if(lindex(i,s1) /= 0) then - st2 = st1 + 1 + lindex_end(i,s1)-lindex(i,s1) - labuf(st1:st2-1) = abuf(lindex(i,s1):lindex_end(i,s1)) - do j=st1,st2-1 - putten(labuf(j)) = .true. - det_minilist(:,:,j) = psi_det_sorted(:,:,labuf(j)) - end do - else - st2 = st1 - end if - - if(sp == 3) then - stamo = 1 - else - stamo = i+1 - end if - - do j=stamo,mo_tot_num - if(bannedOrb(j,s2) .or. banned(i,j)) cycle - if(lindex(j,s2) /= 0) then - st3 = st2 - do k=lindex(j,s2), lindex_end(j,s2) - if(.not. putten(abuf(k))) then - labuf(st3) = abuf(k) - det_minilist(:,:,st3) = psi_det_sorted(:,:,abuf(k)) - st3 += 1 - end if - end do - else - st3 = st2 - end if - - if(indexes(i,j) /= 0) then - st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j) -1!! - labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j)-1) !! - do k=st3, st4-1 - det_minilist(:,:,k) = psi_det_sorted(:,:,labuf(k)) - end do - else - st4 = st3 - end if - !APPLY PART - if(st4 > 1) then - call apply_particles(mask, s1, i, s2, j, alpha, ok, N_int) - !if(.not. ok) stop "non existing alpha......" - !print *, "willcall", st4-1, size(labuf) - call dress_with_alpha_buffer(delta_ij_loc, labuf, det_minilist, st4-1, alpha, iproc) - !call dress_with_alpha_buffer(delta_ij_loc, abuf, siz, alpha, 1) - end if + do i=1,mo_tot_num + if(bannedOrb(i,s1)) cycle + if(lindex(i,s1) /= 0) then + st2 = st1 + 1 + lindex_end(i,s1)-lindex(i,s1) + labuf(st1:st2-1) = abuf(lindex(i,s1):lindex_end(i,s1)) + do j=st1,st2-1 + putten(labuf(j)) = .true. + det_minilist(:,:,j) = psi_det(:,:,labuf(j)) end do - - if(lindex(i,s1) /= 0) then - do j=st1,st2-1 - putten(labuf(j)) = .false. + else + st2 = st1 + end if + + if(sp == 3) then + stamo = 1 + else + stamo = i+1 + end if + + do j=stamo,mo_tot_num + if(bannedOrb(j,s2) .or. banned(i,j)) cycle + if(lindex(j,s2) /= 0) then + st3 = st2 + do k=lindex(j,s2), lindex_end(j,s2) + if(.not. putten(abuf(k))) then + labuf(st3) = abuf(k) + det_minilist(:,:,st3) = psi_det(:,:,abuf(k)) + st3 += 1 + end if end do + else + st3 = st2 end if + if(indexes(i,j) /= 0) then + st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j) -1!! + labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j)-1) !! + do k=st3, st4-1 + det_minilist(:,:,k) = psi_det(:,:,labuf(k)) + end do + else + st4 = st3 + end if + !APPLY PART + if(st4 > 1) then + call apply_particles(mask, s1, i, s2, j, alpha, ok, N_int) + !if(.not. ok) stop "non existing alpha......" + !print *, "willcall", st4-1, size(labuf) + call dress_with_alpha_buffer(delta_ij_loc, labuf, det_minilist, st4-1, alpha, iproc) + !call dress_with_alpha_buffer(delta_ij_loc, abuf, siz, alpha, 1) + end if end do + + if(lindex(i,s1) /= 0) then + do j=st1,st2-1 + putten(labuf(j)) = .false. + end do + end if + + end do end subroutine @@ -543,7 +565,6 @@ subroutine count_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, countedGlob, mobMask(1,1) = iand(negMask(1,1), det(1,1,i)) mobMask(1,2) = iand(negMask(1,2), det(1,2,i)) - if(interesting(i) < i_gen) cycle nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) @@ -656,11 +677,11 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, indexes, ab !if (interesting(i) < 0) then ! stop 'prefetch interesting(i)' !endif + if(interesting(i) < i_gen) cycle mobMask(1,1) = iand(negMask(1,1), det(1,1,i)) mobMask(1,2) = iand(negMask(1,2), det(1,2,i)) - if(interesting(i) < i_gen) cycle nt = popcnt(mobMask(1, 1)) + popcnt(mobMask(1, 2)) if(nt > 4) cycle diff --git a/plugins/mrcc_sto/mrcc_sto.irp.f b/plugins/mrcc_sto/mrcc_sto.irp.f index a1fc237a..09b5d626 100644 --- a/plugins/mrcc_sto/mrcc_sto.irp.f +++ b/plugins/mrcc_sto/mrcc_sto.irp.f @@ -7,28 +7,11 @@ program mrcc_sto call dress_zmq() end - BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ] -&BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ] - implicit none - integer :: i,inpsisor - - idx_non_ref_from_sorted = 0 - psi_from_sorted = 0 - - do i=1,N_det - psi_from_sorted(psi_det_sorted_order(i)) = i - inpsisor = psi_det_sorted_order(i) - if(inpsisor <= 0) stop "idx_non_ref_from_sorted" - idx_non_ref_from_sorted(inpsisor) = idx_non_ref_rev(i) - end do -END_PROVIDER - BEGIN_PROVIDER [ double precision, hij_cache_, (N_det,Nproc) ] &BEGIN_PROVIDER [ double precision, sij_cache_, (N_det,Nproc) ] &BEGIN_PROVIDER [ double precision, dIa_hla_, (N_states,N_det,Nproc) ] &BEGIN_PROVIDER [ double precision, dIa_sla_, (N_states,N_det,Nproc) ] -&BEGIN_PROVIDER [ integer, idx_alpha_, (0:N_det,Nproc) ] BEGIN_DOC ! temporay arrays for dress_with_alpha_buffer. Avoids realocation. END_DOC @@ -71,7 +54,7 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil if(n_minilist == 1) return do i=1,n_minilist - if(idx_non_ref_from_sorted(minilist(i)) == 0) return + if(idx_non_ref_rev(minilist(i)) == 0) return end do shdress = 0d0 @@ -81,40 +64,20 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat endif - ll_sd = 0 do l_sd=1,n_minilist - ok = .true. - k_sd = minilist(l_sd) - !do i_I=1,N_det_ref - ! call get_excitation_degree(psi_det_sorted(1,1,k_sd),psi_ref(1,1,i_I),degree1,N_int) - ! if(degree1 == 0) then - ! ok = .false. - ! exit - ! end if - !end do - - !if(ok) then - ! call get_excitation(psi_det_sorted(1,1,k_sd),alpha,exc,degree1,phase,N_int) - ! if(degree1 == 0 .or. degree1 > 2) stop "minilist error" - ! call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2) - ! - ! ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. & - ! (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') - ! if(ok .and. degree1 == 2) then - ! ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. & - ! (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') - ! end if + !call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int) + !if(degree1 == 0 .or. degree1 > 2) stop "minilist error" + !call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2) + ! + !ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. & + ! (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') + !if(ok .and. degree1 == 2) then + ! ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. & + ! (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') !end if - - if(ok) then - ll_sd += 1 - idx_alpha_(ll_sd,iproc) = k_sd - call i_h_j(alpha,psi_det_sorted(1,1,k_sd),N_int,hij_cache_(k_sd,iproc)) - call get_s2(alpha,psi_det_sorted(1,1,k_sd),N_int,sij_cache_(k_sd,iproc)) - end if + call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc)) + call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc)) enddo - if(ll_sd <= 1) return - idx_alpha_(0,iproc) = ll_sd do i_I=1,N_det_ref @@ -127,62 +90,68 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil dIa(i_state) = 0.d0 enddo - do k_sd=1,idx_alpha_(0,iproc) - call get_excitation_degree(psi_ref(1,1,i_I),psi_det_sorted(1,1,idx_alpha_(k_sd,iproc)),degree,N_int) -! print *, "diden" + do k_sd=1,n_minilist + call get_excitation_degree(psi_ref(1,1,i_I),det_minilist(1,1,k_sd),degree,N_int) if (degree > 2) then cycle endif - call get_excitation(psi_det_sorted(1,1,idx_alpha_(k_sd,iproc)),alpha,exc,degree2,phase,N_int) - !print *, "DEG", degree2 - call decode_exc(exc,degree2,h1,p1,h2,p2,s1,s2) - do k=1,N_int - tmp_det(k,1) = psi_ref(k,1,i_I) - tmp_det(k,2) = psi_ref(k,2,i_I) - enddo + call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int) call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int) - ok2 = .false. - do i_state=1,N_states - dIK(i_state) = dij(i_I, idx_non_ref_from_sorted(idx_alpha_(k_sd,iproc)), i_state) - if(dIK(i_state) /= 0d0) then - ok2 = .true. - endif - enddo - if(.not. ok2) cycle - - ! + if((.not. ok) .and. (.not. perturbative_triples)) cycle + do i_state=1,N_states dka(i_state) = 0.d0 enddo + + ok2 = .false. + !do i_state=1,N_states + ! !if(dka(i_state) == 0) cycle + ! dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state) + ! if(dIk(i_state) /= 0d0) then + ! ok2 = .true. + ! endif + !enddo + !if(.not. ok2) cycle if (ok) then - do l_sd=k_sd+1,idx_alpha_(0,iproc) - call get_excitation_degree(tmp_det,psi_det_sorted(1,1,idx_alpha_(l_sd,iproc)),degree,N_int) + phase2 = 0d0 + do l_sd=k_sd+1,n_minilist + call get_excitation_degree(tmp_det,det_minilist(1,1,l_sd),degree,N_int) if (degree == 0) then - call get_excitation(psi_ref(1,1,i_I),psi_det_sorted(1,1,idx_alpha_(l_sd,iproc)),exc,degree,phase2,N_int) do i_state=1,N_states - dka(i_state) = dij(i_I, idx_non_ref_from_sorted(idx_alpha_(l_sd,iproc)), i_state) * phase * phase2 - enddo + dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state) + if(dIk(i_state) /= 0d0) then + if(phase2 == 0d0) call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,l_sd),exc,degree,phase2,N_int) + dka(i_state) = dij(i_I, idx_non_ref_rev(minilist(l_sd)), i_state) * phase * phase2 + end if + end do + + !call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,l_sd),exc,degree,phase2,N_int) + !do i_state=1,N_states + ! if(dIk(i_state) /= 0d0) dka(i_state) = dij(i_I, idx_non_ref_rev(minilist(l_sd)), i_state) * phase * phase2 + !enddo exit + endif enddo else if (perturbative_triples) then - hka = hij_cache_(idx_alpha_(k_sd,iproc),iproc) - if (dabs(hka) > 1.d-12) then - call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv) + hka = hij_cache_(k_sd,iproc) + if (dabs(hka) > 1.d-12) then + call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv) - do i_state=1,N_states - ASSERT (Delta_E_inv(i_state) < 0.d0) - dka(i_state) = hka / Delta_E_inv(i_state) - enddo - endif + do i_state=1,N_states + ASSERT (Delta_E_inv(i_state) < 0.d0) + dka(i_state) = hka / Delta_E_inv(i_state) + enddo + endif endif - + + if (perturbative_triples.and. (degree2 == 1) ) then call i_h_j(psi_ref(1,1,i_I),tmp_det,N_int,hka) - hka = hij_cache_(idx_alpha_(k_sd,iproc),iproc) - hka + hka = hij_cache_(k_sd,iproc) - hka if (dabs(hka) > 1.d-12) then call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv) do i_state=1,N_states @@ -202,32 +171,20 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil enddo if(.not. ok2) cycle - do l_sd=1,idx_alpha_(0,iproc) - k_sd = idx_alpha_(l_sd,iproc) - hla = hij_cache_(k_sd,iproc) - sla = sij_cache_(k_sd,iproc) + do l_sd=1,n_minilist + k_sd = minilist(l_sd) + hla = hij_cache_(l_sd,iproc) + sla = sij_cache_(l_sd,iproc) do i_state=1,N_states - ! dIa_hla_(i_state,k_sd,iproc) = dIa(i_state) * hla - ! dIa_sla_(i_state,k_sd,iproc) = dIa(i_state) * sla - ! enddo - ! enddo - ! do l_sd=1,idx_alpha_(0,iproc) - ! do i_state=1,N_states - k_sd = idx_alpha_(l_sd,iproc) - m_sd = psi_from_sorted(k_sd) hdress = dIa(i_state) * hla * psi_ref_coef(i_I,i_state) sdress = dIa(i_state) * sla * psi_ref_coef(i_I,i_state) - ! hdress = dIa_hla_(i_state,k_sd,iproc) * psi_ref_coef(i_I,i_state) - ! sdress = dIa_sla_(i_state,k_sd,iproc) * psi_ref_coef(i_I,i_state) !$OMP ATOMIC - delta_ij_loc(i_state,m_sd,1) += hdress + delta_ij_loc(i_state,k_sd,1) += hdress !$OMP ATOMIC - delta_ij_loc(i_state,m_sd,2) += sdress - !print *, "ENDRES" + delta_ij_loc(i_state,k_sd,2) += sdress enddo enddo enddo - end subroutine