From b7787f5e6dce198bee06eb92f69b9904a7448bea Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 7 May 2024 19:43:05 +0200 Subject: [PATCH] trying to speed up the PT2 in TC by transposing the array of tc integrals --- .../local/bi_ort_ints/total_twoe_pot.irp.f | 8 +- .../cipsi_tc_bi_ortho/get_d0_transp.irp.f | 140 +++++++++++ .../local/cipsi_tc_bi_ortho/get_d2_good.irp.f | 3 - .../cipsi_tc_bi_ortho/get_d2_transp.irp.f | 235 ++++++++++++++++++ plugins/local/cipsi_tc_bi_ortho/pt2.irp.f | 1 + .../local/cipsi_tc_bi_ortho/selection.irp.f | 94 +------ .../cipsi_tc_bi_ortho/stochastic_cipsi.irp.f | 3 + plugins/local/fci_tc_bi/pt2_tc.irp.f | 2 + .../local/tc_bi_ortho/e_corr_bi_ortho.irp.f | 1 - plugins/local/tc_keywords/EZFIO.cfg | 11 +- 10 files changed, 404 insertions(+), 94 deletions(-) create mode 100644 plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f create mode 100644 plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f diff --git a/plugins/local/bi_ort_ints/total_twoe_pot.irp.f b/plugins/local/bi_ort_ints/total_twoe_pot.irp.f index 1e127fac..71269fdc 100644 --- a/plugins/local/bi_ort_ints/total_twoe_pot.irp.f +++ b/plugins/local/bi_ort_ints/total_twoe_pot.irp.f @@ -259,15 +259,21 @@ BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_transp, (mo_num, mo_num, END_DOC integer :: i,j,k,l + print*,'Providing mo_bi_ortho_tc_two_e_transp' + double precision :: t0,t1 + call wall_time(t0) do i = 1, mo_num do j = 1, mo_num do k = 1, mo_num do l = 1, mo_num - mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e_transp(k,l,i,j) + mo_bi_ortho_tc_two_e_transp(i,j,k,l) = mo_bi_ortho_tc_two_e(k,l,i,j) enddo enddo enddo enddo + call wall_time(t1) + + print *, ' WALL TIME for PROVIDING mo_bi_ortho_tc_two_e_transp (min)', (t1-t0)/60.d0 END_PROVIDER ! --- diff --git a/plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f b/plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f new file mode 100644 index 00000000..56238e13 --- /dev/null +++ b/plugins/local/cipsi_tc_bi_ortho/get_d0_transp.irp.f @@ -0,0 +1,140 @@ +subroutine get_d0_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs) + !todo: indices/conjg should be okay for complex + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2) + integer(bit_kind), intent(in) :: phasemask(N_int,2) + logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) + integer(bit_kind) :: det(N_int, 2) + double precision, intent(in) :: coefs(N_states,2) + double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num) + double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num) + integer, intent(in) :: h(0:2,2), p(0:4,2), sp + + integer :: i, j, k, s, h1, h2, p1, p2, puti, putj, mm + double precision :: phase + double precision :: hij,hji + double precision, external :: get_phase_bi + logical :: ok + + integer, parameter :: bant=1 + double precision, allocatable :: hij_cache1(:), hij_cache2(:) + allocate (hij_cache1(mo_num),hij_cache2(mo_num)) + double precision, allocatable :: hji_cache1(:), hji_cache2(:) + allocate (hji_cache1(mo_num),hji_cache2(mo_num)) +! print*,'in get_d0_new' +! call debug_det(gen,N_int) +! print*,'coefs',coefs(1,:) + + 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_complex(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) + do mm = 1, mo_num + hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) + hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,p1,h2,h1) + enddo + !!!!!!!!!! + do p2=1, mo_num + if(bannedOrb(p2,2)) cycle + if(banned(p1, p2, bant)) cycle ! rentable? + if(p1 == h1 .or. p2 == h2) then + call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) + ! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this +! call i_h_j_complex(det, gen, N_int, hij) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) + else + phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) + hij = hij_cache1(p2) * phase + end if + if (hij == (0.d0,0.d0)) cycle + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,2) * hij ! HOTSPOT + enddo + end do + !!!!!!!!!! + do p2=1, mo_num + if(bannedOrb(p2,2)) cycle + if(banned(p1, p2, bant)) cycle ! rentable? + if(p1 == h1 .or. p2 == h2) then + call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) + ! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this +! call i_h_j_complex(det, gen, N_int, hij) + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) + else + phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) + hji = hji_cache1(p2) * phase + end if + if (hji == (0.d0,0.d0)) cycle + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,1) * hji ! HOTSPOT + enddo + end do + end do + + else ! AA BB + p1 = p(1,sp) + p2 = p(2,sp) + do puti=1, mo_num + if(bannedOrb(puti, sp)) cycle +! call get_mo_two_e_integrals_complex(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map,mo_integrals_map_2) +! call get_mo_two_e_integrals_complex(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map,mo_integrals_map_2) + do mm = 1, mo_num + hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,puti,p2,p1) + hij_cache2(mm) = mo_bi_ortho_tc_two_e(mm,puti,p1,p2) + hji_cache1(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p2,p1) + hji_cache2(mm) = mo_bi_ortho_tc_two_e_transp(mm,puti,p1,p2) + enddo + !!!!!!!!!! + do putj=puti+1, mo_num + if(bannedOrb(putj, sp)) cycle + if(banned(puti, putj, bant)) cycle ! rentable? + if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then + call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) + !call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this +! call i_h_j_complex(det, gen, N_int, hij) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) + if (hij == 0.d0) cycle + else +! hij = (mo_two_e_integral_complex(p1, p2, puti, putj) - mo_two_e_integral_complex(p2, p1, puti, putj)) +! hij = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) - mo_bi_ortho_tc_two_e(p2, p1, puti, putj)) + hij = (mo_bi_ortho_tc_two_e(puti, putj, p1, p2) - mo_bi_ortho_tc_two_e(puti, putj, p2, p1)) + if (hij == 0.d0) cycle + hij = (hij) * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) + end if + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij + enddo + end do + + !!!!!!!!!! + do putj=puti+1, mo_num + if(bannedOrb(putj, sp)) cycle + if(banned(puti, putj, bant)) cycle ! rentable? + if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then + call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int) + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) + if (hji == 0.d0) cycle + else +! hji = (mo_bi_ortho_tc_two_e( p1, p2, puti, putj) - mo_bi_ortho_tc_two_e( p2, p1, puti, putj)) + hji = (mo_bi_ortho_tc_two_e_transp(puti, putj, p1, p2 ) - mo_bi_ortho_tc_two_e_transp( puti, putj, p2, p1)) + if (hji == 0.d0) cycle + hji = (hji) * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) + end if + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji + enddo + end do + end do + end if + + deallocate(hij_cache1,hij_cache2) +end + diff --git a/plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f b/plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f index d01ed433..86922ae9 100644 --- a/plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f @@ -25,9 +25,6 @@ subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, integer :: bant bant = 1 -! print*, 'in get_d2_new' -! call debug_det(gen,N_int) -! print*,'coefs',coefs(1,:) tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles diff --git a/plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f b/plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f new file mode 100644 index 00000000..b2a7ea31 --- /dev/null +++ b/plugins/local/cipsi_tc_bi_ortho/get_d2_transp.irp.f @@ -0,0 +1,235 @@ + +subroutine get_d2_new_transp(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs) + !todo: indices/conjg should be correct for complex + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2) + integer(bit_kind), intent(in) :: phasemask(N_int,2) + logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) + double precision, intent(in) :: coefs(N_states,2) + double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num) + double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num) + integer, intent(in) :: h(0:2,2), p(0:4,2), sp + + double precision, external :: get_phase_bi + + integer :: i, j, k, tip, ma, mi, puti, putj + integer :: h1, h2, p1, p2, i1, i2 + double precision :: phase + double precision :: hij,hji + + integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/)) + integer, parameter :: turn2(2) = (/2, 1/) + integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/)) + + integer :: bant + bant = 1 + + tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles + + ma = sp !1:(alpha,alpha); 2:(b,b); 3:(a,b) + if(p(0,1) > p(0,2)) ma = 1 ! more alpha particles than beta particles + if(p(0,1) < p(0,2)) ma = 2 ! fewer alpha particles than beta particles + mi = mod(ma, 2) + 1 + + if(sp == 3) then ! if one alpha and one beta xhole + !(where xholes refer to the ionizations from the generator, not the holes occupied in the ionized generator) + if(ma == 2) bant = 2 ! if more beta particles than alpha particles + + if(tip == 3) then ! if 3 of one particle spin and 1 of the other particle spin + puti = p(1, mi) + if(bannedOrb(puti, mi)) return + h1 = h(1, ma) + h2 = h(2, ma) + + !! + do i = 1, 3 ! loop over all 3 combinations of 2 particles with spin ma + putj = p(i, ma) + if(banned(putj,puti,bant)) cycle + i1 = turn3(1,i) + i2 = turn3(2,i) + p1 = p(i1, ma) + p2 = p(i2, ma) + + ! |G> = |psi_{gen,i}> + ! |G'> = a_{x1} a_{x2} |G> + ! |alpha> = a_{puti}^{\dagger} a_{putj}^{\dagger} |G'> + ! |alpha> = t_{x1,x2}^{puti,putj} |G> + ! hij = + ! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}> + !todo: = ( - ) * phase + ! += dconjg(c_i) * + ! = ( - ) * phase + ! += * c_i + +!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! + ! take the transpose of what's written above because later use the complex conjugate + +! hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e( h1, h2, p2, p1) +! hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e_transp( h1, h2, p2, p1) + hij = mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp( p1, p2, h2, h1) + hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e( p1, p2, h2, h1) + if (hij == 0.d0.or.hji==0.d0) cycle + + ! take conjugate to get contribution to instead of +! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + hij = hij * phase + hji = hji * phase + + if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij + mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji + enddo + else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji + enddo + end if + end do + else ! if 2 alpha and 2 beta particles + h1 = h(1,1) + h2 = h(1,2) + !! + do j = 1,2 ! loop over all 4 combinations of one alpha and one beta particle + putj = p(j, 2) + if(bannedOrb(putj, 2)) cycle + p2 = p(turn2(j), 2) + do i = 1,2 + puti = p(i, 1) + if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle + p1 = p(turn2(i), 1) + ! hij = +! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) +!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! + ! take the transpose of what's written above because later use the complex conjugate +! hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2 ) +! hji = mo_bi_ortho_tc_two_e_transp(h1, h2, p1, p2 ) + hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2 ) + hji = mo_bi_ortho_tc_two_e( p1, p2, h1, h2) + if (hij /= 0.d0.or.hji==0.d0) then + ! take conjugate to get contribution to instead of +! hij = dconjg(hij) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) + phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) + hij = hij * phase + hji = hji * phase + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji + enddo + endif + end do + end do + end if + + else ! if holes are (a,a) or (b,b) + if(tip == 0) then ! if particles are (a,a,a,a) or (b,b,b,b) + h1 = h(1, ma) + h2 = h(2, ma) + !! + do i=1,3 + puti = p(i, ma) + if(bannedOrb(puti,ma)) cycle + do j=i+1,4 + putj = p(j, ma) + if(bannedOrb(putj,ma)) cycle + if(banned(puti,putj,1)) cycle + + i1 = turn2d(1, i, j) + i2 = turn2d(2, i, j) + p1 = p(i1, ma) + p2 = p(i2, ma) +! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2) +!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! + ! take the transpose of what's written above because later use the complex conjugate + hij = mo_bi_ortho_tc_two_e_transp(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e_transp(p1, p2, h2,h1 ) + hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p1, p2, h2,h1 ) + if (hij == 0.d0.or.hji == 0.d0) cycle + + ! take conjugate to get contribution to instead of +! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + phase = get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) + hij = hij * phase + hji = hji * phase + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hij + mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hji + enddo + end do + end do + else if(tip == 3) then ! if particles are (a,a,a,b) (ma=1,mi=2) or (a,b,b,b) (ma=2,mi=1) + h1 = h(1, mi) + h2 = h(1, ma) + p1 = p(1, mi) + !! + do i=1,3 + puti = p(turn3(1,i), ma) + if(bannedOrb(puti,ma)) cycle + putj = p(turn3(2,i), ma) + if(bannedOrb(putj,ma)) cycle + if(banned(puti,putj,1)) cycle + p2 = p(i, ma) + +! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) +!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! + ! take the transpose of what's written above because later use the complex conjugate + hij = mo_bi_ortho_tc_two_e_transp(p1, p2 ,h1, h2) + hji = mo_bi_ortho_tc_two_e(p1, p2,h1, h2 ) + if (hij == 0.d0) cycle + + ! take conjugate to get contribution to instead of +! hij = dconjg(hij) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int) + phase = get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int) + hij = hij * phase + hji = hji * phase + if (puti < putj) then + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji + enddo + else + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hij + mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hji + enddo + endif + end do + else ! tip == 4 (a,a,b,b) + puti = p(1, sp) + putj = p(2, sp) + if(.not. banned(puti,putj,1)) then + p1 = p(1, mi) + p2 = p(2, mi) + h1 = h(1, mi) + h2 = h(2, mi) + !! +! hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)) +!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!! + ! take the transpose of what's written above because later use the complex conjugate + hij = (mo_bi_ortho_tc_two_e_transp(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e_transp(p2,p1,h1, h2)) + hji = (mo_bi_ortho_tc_two_e(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e(p2,p1,h1, h2)) + if (hij /= 0.d0.or.hji==0.d0) then + ! take conjugate to get contribution to instead of +! hij = dconjg(hij) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) + phase = get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) + hij = hij * phase + hji = hji* phase + !DIR$ LOOP COUNT AVG(4) + do k=1,N_states + mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hij + mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hji + enddo + end if + end if + end if + end if +end diff --git a/plugins/local/cipsi_tc_bi_ortho/pt2.irp.f b/plugins/local/cipsi_tc_bi_ortho/pt2.irp.f index 833cc0ea..ada19c6b 100644 --- a/plugins/local/cipsi_tc_bi_ortho/pt2.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/pt2.irp.f @@ -67,6 +67,7 @@ subroutine tc_pt2 call pt2_alloc(pt2_data_err, N_states) call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) + call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2) end diff --git a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f index 0b4345d5..0f785ba2 100644 --- a/plugins/local/cipsi_tc_bi_ortho/selection.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/selection.irp.f @@ -636,10 +636,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere negMask(i,2) = not(mask(i,2)) end do -! print*,'in selection ' do i = 1, N_sel -! call debug_det(det(1,1,i),N_int) -! print*,i,dabs(psi_selectors_coef_transp_tc(1,2,i) * psi_selectors_coef_transp_tc(1,1,i)) if(interesting(i) < 0) then stop 'prefetch interesting(i) and det(i)' endif @@ -691,11 +688,19 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) if(nt == 4) then - call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + if(transpose_two_e_int)then + call get_d2_new_transp(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + else + call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + endif elseif(nt == 3) then call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) else + if(transpose_two_e_int)then + call get_d0_transp (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + else call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) + endif endif elseif(nt == 4) then call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) @@ -887,79 +892,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii) do istate = 1,N_states delta_E = E0(istate) - Hii + E_shift - double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp, error - if(debug_tc_pt2 == 1)then !! Using the old version - psi_h_alpha = 0.d0 - alpha_h_psi = 0.d0 - do iii = 1, N_det_selectors - call htilde_mu_mat_opt_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_opt_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - call get_excitation_degree(psi_selectors(1,1,iii), det,degree,N_int) - if(degree == 0)then - print*,'problem !!!' - print*,'a determinant is already in the wave function !!' - print*,'it corresponds to the selector number ',iii - call debug_det(det,N_int) - stop - endif -! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) -! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function - alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function - enddo - else if(debug_tc_pt2 == 2)then !! debugging the new version -! psi_h_alpha_tmp = 0.d0 -! alpha_h_psi_tmp = 0.d0 -! do iii = 1, N_det_selectors ! old version -! call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) -! call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) -! psi_h_alpha_tmp += i_h_alpha * psi_selectors_coef_tc(iii,1,1) ! left function -! alpha_h_psi_tmp += alpha_h_i * psi_selectors_coef_tc(iii,2,1) ! right function -! enddo - psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version - alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version - psi_h_alpha = 0.d0 - alpha_h_psi = 0.d0 - do iii = 1, N_det ! old version - call htilde_mu_mat_opt_bi_ortho_no_3e(psi_det(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_det(1,1,iii), N_int, alpha_h_i) - psi_h_alpha += i_h_alpha * psi_l_coef_bi_ortho(iii,1) ! left function - alpha_h_psi += alpha_h_i * psi_r_coef_bi_ortho(iii,1) ! right function - enddo - if(dabs(psi_h_alpha*alpha_h_psi/delta_E).gt.1.d-10)then - error = dabs(psi_h_alpha * alpha_h_psi - psi_h_alpha_tmp * alpha_h_psi_tmp)/dabs(psi_h_alpha * alpha_h_psi) - if(error.gt.1.d-2)then - call debug_det(det, N_int) - print*,'error =',error,psi_h_alpha * alpha_h_psi/delta_E,psi_h_alpha_tmp * alpha_h_psi_tmp/delta_E - print*,psi_h_alpha , alpha_h_psi - print*,psi_h_alpha_tmp , alpha_h_psi_tmp - print*,'selectors ' - do iii = 1, N_det_selectors ! old version - print*,'iii',iii,psi_selectors_coef_tc(iii,1,1),psi_selectors_coef_tc(iii,2,1) - call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha) - call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i) - print*,i_h_alpha,alpha_h_i - call debug_det(psi_selectors(1,1,iii),N_int) - enddo -! print*,'psi_det ' -! do iii = 1, N_det! old version -! print*,'iii',iii,psi_l_coef_bi_ortho(iii,1),psi_r_coef_bi_ortho(iii,1) -! call debug_det(psi_det(1,1,iii),N_int) -! enddo - stop - endif - endif - else - psi_h_alpha = mat_l(istate, p1, p2) - alpha_h_psi = mat_r(istate, p1, p2) - endif + psi_h_alpha = mat_l(istate, p1, p2) + alpha_h_psi = mat_r(istate, p1, p2) val = 4.d0 * psi_h_alpha * alpha_h_psi tmp = dsqrt(delta_E * delta_E + val) -! if (delta_E < 0.d0) then -! tmp = -tmp -! endif e_pert(istate) = 0.25 * val / delta_E -! e_pert(istate) = 0.5d0 * (tmp - delta_E) if(dsqrt(tmp).gt.1.d-4.and.dabs(psi_h_alpha).gt.1.d-4)then coef(istate) = e_pert(istate) / psi_h_alpha else @@ -976,15 +913,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if(e_pert(istate).gt.0.d0)e_pert(istate)=0.d0 endif -! if(selection_tc == 1 )then -! if(e_pert(istate).lt.0.d0)then -! e_pert(istate) = 0.d0 -! endif -! else if(selection_tc == -1)then -! if(e_pert(istate).gt.0.d0)then -! e_pert(istate) = 0.d0 -! endif -! endif enddo diff --git a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 99a8de7e..bb5a89a1 100644 --- a/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -88,6 +88,9 @@ subroutine run_stochastic_cipsi call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) + if(transpose_two_e_int)then + provide mo_bi_ortho_tc_two_e_transp + endif call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection ! stop diff --git a/plugins/local/fci_tc_bi/pt2_tc.irp.f b/plugins/local/fci_tc_bi/pt2_tc.irp.f index 390042bf..3c07e367 100644 --- a/plugins/local/fci_tc_bi/pt2_tc.irp.f +++ b/plugins/local/fci_tc_bi/pt2_tc.irp.f @@ -13,6 +13,8 @@ program tc_pt2_prog pruning = -1.d0 touch pruning + read_wf = .True. + touch read_wf ! pt2_relative_error = 0.01d0 ! touch pt2_relative_error diff --git a/plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f b/plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f index 4abdc25b..5a3971c5 100644 --- a/plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f +++ b/plugins/local/tc_bi_ortho/e_corr_bi_ortho.irp.f @@ -27,7 +27,6 @@ if(degree == 1)then e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij else -! print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij endif endif diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg index 1e89eaa4..39968ec8 100644 --- a/plugins/local/tc_keywords/EZFIO.cfg +++ b/plugins/local/tc_keywords/EZFIO.cfg @@ -184,12 +184,6 @@ doc: Read/Write normal_two_body_bi_orth from/to disk [ Write | Read | None ] interface: ezfio,provider,ocaml default: None -[debug_tc_pt2] -type: integer -doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body -interface: ezfio,provider,ocaml -default: -1 - [only_spin_tc_right] type: logical doc: If |true|, only the right part of WF is used to compute spin dens @@ -268,3 +262,8 @@ doc: Thresholds on the Imag part of TC energy interface: ezfio,provider,ocaml default: 1.e-7 +[transpose_two_e_int] +type: logical +doc: If |true|, you duplicate the two-electron TC integrals with the transpose matrix. Acceleates the PT2. +interface: ezfio,provider,ocaml +default: False