From f0178d09a2572990111211fbe4b894bde9fa05ee Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 19 Jan 2023 22:34:11 +0100 Subject: [PATCH] diagonal matrix elements work with 3-e a la fock --- src/bi_ort_ints/three_body_ijm.irp.f | 2 +- src/tc_bi_ortho/slater_tc_3e.irp.f | 13 ----- src/tc_bi_ortho/slater_tc_opt.irp.f | 78 ++++++++++++++++++++++++-- src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 18 ++++-- 4 files changed, 88 insertions(+), 23 deletions(-) diff --git a/src/bi_ort_ints/three_body_ijm.irp.f b/src/bi_ort_ints/three_body_ijm.irp.f index 0e42264b..4d21cb93 100644 --- a/src/bi_ort_ints/three_body_ijm.irp.f +++ b/src/bi_ort_ints/three_body_ijm.irp.f @@ -60,7 +60,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num ! ! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation ! - ! three_e_3_idx_direct_bi_ort(m,j,i) = + ! three_e_3_idx_cycle_1_bi_ort(m,j,i) = ! ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign ! diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index a56a432f..0d5f8542 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -49,8 +49,6 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) if(Ne(1)+Ne(2).ge.3)then !! ! alpha/alpha/beta three-body - double precision :: accu - accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) do j = i+1, Ne(1) @@ -62,14 +60,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR hthree += direct_int - exchange_int - accu += direct_int - exchange_int enddo enddo enddo - !print*,'aab = ',accu ! beta/beta/alpha three-body - accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) do j = i+1, Ne(2) @@ -79,14 +74,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) hthree += direct_int - exchange_int - accu += direct_int - exchange_int enddo enddo enddo - !print*,'abb = ',accu ! alpha/alpha/alpha three-body - accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) ! 1 do j = i+1, Ne(1) @@ -95,14 +87,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,1) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS - accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo - !print*,'aaa = ',accu ! beta/beta/beta three-body - accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) ! 1 do j = i+1, Ne(2) @@ -111,11 +100,9 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,2) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS - accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo - !print*,'bbb = ',accu endif end diff --git a/src/tc_bi_ortho/slater_tc_opt.irp.f b/src/tc_bi_ortho/slater_tc_opt.irp.f index 50886fe3..4048f481 100644 --- a/src/tc_bi_ortho/slater_tc_opt.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt.irp.f @@ -55,6 +55,9 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, enddo if (nexc(1)+nexc(2) == 0) then + hmono = ref_tc_energy_1e + htwoe = ref_tc_energy_2e + hthree= ref_tc_energy_3e htot = ref_tc_energy_tot return endif @@ -75,7 +78,6 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, hmono = ref_tc_energy_1e htwoe = ref_tc_energy_2e hthree= ref_tc_energy_3e - do ispin=1,2 na = elec_num_tab(ispin) nb = elec_num_tab(iand(ispin,1)+1) @@ -110,7 +112,9 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer :: occ(Nint*bit_kind_size,2) integer :: other_spin - integer :: k,l,i + integer :: k,l,i,jj,mm,j,m + double precision :: three_e_diag_parrallel_spin, direct_int, exchange_int + if (iorb < 1) then print *, irp_here, ': iorb < 1' @@ -151,6 +155,39 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) do i=1,nb htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo + + if(three_body_h_tc)then + !!!!! 3-e part + !! same-spin/same-spin + do j = 1, na + jj = occ(j,ispin) + do m = j+1, na + mm = occ(m,ispin) + hthree += three_e_diag_parrallel_spin(mm,jj,iorb) + enddo + enddo + !! same-spin/oposite-spin + do j = 1, na + jj = occ(j,ispin) + do m = 1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo + enddo + !! oposite-spin/opposite-spin + do j = 1, nb + jj = occ(j,other_spin) + do m = j+1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo + enddo + endif + na = na+1 end @@ -172,10 +209,11 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) integer, intent(inout) :: na, nb integer(bit_kind), intent(inout) :: key(Nint,2) double precision, intent(inout) :: hmono,htwoe,hthree - + + double precision :: direct_int, exchange_int, three_e_diag_parrallel_spin integer :: occ(Nint*bit_kind_size,2) integer :: other_spin - integer :: k,l,i + integer :: k,l,i,jj,mm,j,m integer :: tmp(2) ASSERT (iorb > 0) @@ -205,5 +243,37 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo + if(three_body_h_tc)then + !!!!! 3-e part + !! same-spin/same-spin + do j = 1, na + jj = occ(j,ispin) + do m = j+1, na + mm = occ(m,ispin) + hthree -= three_e_diag_parrallel_spin(mm,jj,iorb) + enddo + enddo + !! same-spin/oposite-spin + do j = 1, na + jj = occ(j,ispin) + do m = 1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree -= (direct_int - exchange_int) + enddo + enddo + !! oposite-spin/opposite-spin + do j = 1, nb + jj = occ(j,other_spin) + do m = j+1, nb + mm = occ(m,other_spin) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR + hthree -= (direct_int - exchange_int) + enddo + enddo + endif + end 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 069a1d53..dda4bd00 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -19,28 +19,36 @@ subroutine test_slater_tc_opt integer :: i,j double precision :: hmono, htwoe, htot, hthree double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot - double precision :: accu ,i_count + double precision :: accu_d ,i_count, accu accu = 0.d0 + accu_d = 0.d0 i_count = 0.d0 do i = 1, N_det ! do i = 14,14 call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot) + call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) +! print*,hthree,hnewthnewree +! print*,htot,hnewtot,dabs(hnewtot-htot) + accu_d += dabs(htot-hnewtot) +! if(dabs(htot-hnewtot).gt.1.d-8)then + print*,i + print*,htot,hnewtot,dabs(htot-hnewtot) +! endif do j = 1, N_det -! do j = 1, 1 if(i==j)cycle call single_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, htot) call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot) if(dabs(htot).gt.1.d-10)then -! if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then + if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then print*,j,i i_count += 1.D0 print*,htot,hnewtot,dabs(htot-hnewtot) accu += dabs(htot-hnewtot) -! endif + endif endif enddo enddo - print*,'accu = ',accu/i_count + print*,'accu_d = ',accu_d/N_det end