Skip to content

Commit

Permalink
Removed hij_cache giving wrong PT2
Browse files Browse the repository at this point in the history
  • Loading branch information
@scemama
scemama committed Oct 18, 2024
1 parent fda1d77 commit d77f48b
Showing 1 changed file with 23 additions and 21 deletions.
44 changes: 23 additions & 21 deletions src/cipsi/selection.irp.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
logical, allocatable :: banned(:,:,:), bannedOrb(:,:)
double precision, allocatable :: coef_fullminilist_rev(:,:)
double precision, allocatable :: mat(:,:,:), hij_cache(:,:,:)
double precision, allocatable :: mat(:,:,:)


PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
Expand All @@ -205,7 +205,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
! Removed to avoid introducing determinants already presents in the wf
!double precision, parameter :: norm_thr = 1.d-16

allocate (indices(N_det), hij_cache(mo_num,mo_num,2), &
allocate (indices(N_det), &
exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))

! Pre-compute excitation degrees wrt alpha determinants
Expand Down Expand Up @@ -511,15 +511,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc

maskInd = maskInd + 1
if(mod(maskInd, csubset) == (subset-1)) then
call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache(1,1,1),mo_integrals_map)
if (sp /= 3) then ! AA or BB
call get_mo_two_e_integrals_ij(h1,h2,mo_num,hij_cache(1,1,2),mo_integrals_map)
endif

call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
if(fullMatch) cycle

call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, hij_cache)
call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)

call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
end if
Expand All @@ -535,7 +530,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
enddo
enddo
deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
deallocate(banned, bannedOrb, mat, hij_cache)
deallocate(banned, bannedOrb, mat)
end subroutine

BEGIN_TEMPLATE
Expand Down Expand Up @@ -765,6 +760,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
enddo
do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1
!do_diag = .False.
double precision :: eigvalues(N_states+1)
double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2)
Expand Down Expand Up @@ -918,7 +914,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
END_TEMPLATE
subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, hij_cache)
subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
use bitmasks
implicit none
BEGIN_DOC
Expand All @@ -930,7 +926,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
integer, intent(in) :: sp, i_gen, N_sel
integer, intent(in) :: interesting(0:N_sel)
integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
Expand Down Expand Up @@ -1000,9 +995,9 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
if(nt == 4) then
call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
else if(nt == 3) then
call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
else
call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
end if
else if(nt == 4) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
Expand Down Expand Up @@ -1203,7 +1198,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
end
subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
Expand All @@ -1214,7 +1209,6 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok
Expand Down Expand Up @@ -1256,11 +1250,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache1(1,1),mo_integrals_map)
call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache1(1,2),mo_integrals_map)
tmp_row = 0d0
do putj=1, hfix-1
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(hfix,putj,1) - hij_cache(putj,hfix,1)
hij = hij_cache1(putj,1) - hij_cache1(putj,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
Expand All @@ -1272,7 +1268,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
do putj=hfix+1, mo_num
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,hfix,1) - hij_cache(hfix,putj,1)
hij = hij_cache1(putj,2) - hij_cache1(putj,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
Expand Down Expand Up @@ -1463,7 +1459,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
Expand All @@ -1474,7 +1470,6 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, intent(in) :: hij_cache(mo_num, mo_num, 2)
integer :: i, j, k, s, h1, h2, p1, p2, puti, putj
double precision :: hij, phase
Expand All @@ -1483,13 +1478,17 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
integer, parameter :: bant=1
double precision, allocatable :: hij_cache1(:), hij_cache2(:)
allocate (hij_cache1(mo_num),hij_cache2(mo_num))
PROVIDE mo_integrals_threshold
if(sp == 3) then ! AB
h1 = p(1,1)
h2 = p(1,2)
do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle
call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
Expand All @@ -1498,7 +1497,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
call i_h_j(gen, det, N_int, hij)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij_cache(p2,p1,1) * phase
hij = hij_cache1(p2) * phase
end if
if (dabs(hij) < mo_integrals_threshold) cycle
!DIR$ LOOP COUNT AVG(4)
Expand All @@ -1513,6 +1512,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
p2 = p(2,sp)
do puti=1, mo_num
if (bannedOrb(puti, sp)) cycle
call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map)
call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map)
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
Expand All @@ -1521,7 +1522,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
call i_h_j(gen, det, N_int, hij)
if (dabs(hij) < mo_integrals_threshold) cycle
else
hij = hij_cache(putj,puti,1) - hij_cache(putj,puti,2)
hij = hij_cache1(putj) - hij_cache2(putj)
if (dabs(hij) < mo_integrals_threshold) cycle
hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if
Expand All @@ -1533,6 +1534,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
end do
end if
deallocate(hij_cache1,hij_cache2)
end
Expand Down

0 comments on commit d77f48b

Please sign in to comment.