From 9eba8d692d549e24d30e3236d432e800b7551249 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 20 Jan 2023 16:33:37 +0100 Subject: [PATCH] optimized all matrix elements with three body terms --- src/tc_bi_ortho/slater_tc_opt.irp.f | 12 +- src/tc_bi_ortho/slater_tc_opt_double.irp.f | 197 +++++++++++++++++++++ src/tc_bi_ortho/test_normal_order.irp.f | 4 +- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 3 +- 4 files changed, 205 insertions(+), 11 deletions(-) diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index c334b274..8ab3388c 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -28,19 +28,13 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree.gt.2) return if(degree == 0)then - call diag_htilde_mu_mat_fock_bi_ortho(Nint, key_i, hmono, htwoe, hthree, htot) + call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1)then - call single_htilde_mu_mat_fock_bi_ortho (Nint,key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint,key_j, key_i , hmono, htwoe, hthree, htot) else if(degree == 2)then - call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) - if(three_body_h_tc) then - if(.not.double_normal_ord) then - call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) - endif - endif + call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) endif - htot = hmono + htwoe + hthree if(degree==0) then htot += nuclear_repulsion endif diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index ef319c47..ca1d0eea 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -1,3 +1,200 @@ + +subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) + + BEGIN_DOC + ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int,phase + + + call get_excitation_degree(key_i, key_j, degree, Nint) + + hmono = 0.d0 + htwoe = 0.d0 + hthree = 0.d0 + htot = 0.d0 + + if(degree.ne.2)then + return + endif + + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + call get_double_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + + if(s1.ne.s2)then + ! opposite spin two-body +! key_j, key_i + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + if(three_body_h_tc)then + if(.not.double_normal_ord)then + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + elseif(double_normal_ord.and.+Ne(1).gt.2)then + htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? + endif + endif + else + ! same spin two-body + ! direct terms + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + ! exchange terms + htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) + if(three_body_h_tc)then + if(.not.double_normal_ord)then + call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + elseif(double_normal_ord.and.+Ne(1).gt.2)then + htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? + htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? + endif + endif + endif + hthree *= phase + htwoe *= phase + htot = htwoe + hthree + +end + + + +subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) + implicit none + integer(bit_kind), intent(in) :: key_i(N_int,2) + integer, intent(in) :: h1,h2,p1,p2,s1,s2 + double precision, intent(out) :: hthree + integer :: nexc(2),i,ispin,na,nb + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: particle(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_particle(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) + integer(bit_kind) :: det_tmp(N_int,2) + integer :: ipart, ihole + double precision :: direct_int, exchange_int, three_e_double_parrallel_spin + + nexc(1) = 0 + nexc(2) = 0 + !! Get all the holes and particles of key_i with respect to the ROHF determinant + do i=1,N_int + hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),key_i(i,1)) + particle(i,2) = iand(hole(i,2),key_i(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) = nexc(1) + popcnt(hole(i,1)) + nexc(2) = nexc(2) + popcnt(hole(i,2)) + enddo + integer :: tmp(2) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha + ASSERT (tmp(2) == nexc(2)) ! Number of particle beta + !DIR$ FORCEINLINE + call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) + ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha + ASSERT (tmp(2) == nexc(2)) ! Number of holes beta + if(s1==s2.and.s1==1)then + !!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc + hthree = eff_2_e_from_3_e_aa(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 1 ! i==alpha ==> pure same spin terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ipart=occ_particle(i,ispin) + hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + ihole=occ_hole(i,ispin) + hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + enddo + ispin = 2 ! i==beta ==> alpha/alpha/beta terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and (h2,p2) + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + elseif(s1==s2.and.s1==2)then + !!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc + hthree = eff_2_e_from_3_e_bb(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 2 ! i==beta ==> pure same spin terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ipart=occ_particle(i,ispin) + hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1) + ihole=occ_hole(i,ispin) + hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1) + enddo + ispin = 1 ! i==alpha==> beta/beta/alpha terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and (h2,p2) + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + else ! (h1,p1) == alpha/(h2,p2) == beta + hthree = eff_2_e_from_3_e_ab(p2,p1,h2,h1) + if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant + !!!!!!!! the matrix element is already exact + !!!!!!!! else you need to take care of holes and particles + !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!! + ispin = 1 ! i==alpha ==> alpha/beta/alpha terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h1,p1) and i + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch13_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch13_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + ispin = 2 ! i==beta ==> alpha/beta/beta terms + do i = 1, nexc(ispin) ! number of couple of holes/particles + ! exchange between (h2,p2) and i + ipart=occ_particle(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch23_bi_ort(ipart,p2,h2,p1,h1) + hthree += direct_int - exchange_int + ihole=occ_hole(i,ispin) + direct_int = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch23_bi_ort(ihole,p2,h2,p1,h1) + hthree -= direct_int - exchange_int + enddo + endif +end + + BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)] implicit none BEGIN_DOC diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f index 46705f5f..118e481a 100644 --- a/src/tc_bi_ortho/test_normal_order.irp.f +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -40,13 +40,15 @@ subroutine test call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) hthree *= phase ! !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) - normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) + call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal) +! normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1) accu += dabs(hthree-normal) enddo enddo enddo enddo print*,'accu opposite spin = ',accu +stop ! p2=6 ! p1=5 diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index 8afb3b25..551eced2 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -50,13 +50,14 @@ subroutine test_slater_tc_opt accu += dabs(htot-hnewtot) ! if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then - print*,j,i + print*,j,i,degree call debug_det(psi_det(1,1,i),N_int) call debug_det(psi_det(1,1,j),N_int) print*,htot,hnewtot,dabs(htot-hnewtot) ! print*,hthree,hnewthree,dabs(hthree-hnewthree) stop endif + print*,htot,hnewtot,dabs(htot-hnewtot) endif enddo enddo