PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-19 04:22:32 +01:00
Code Releases Activity

Merge pull request #235 from eginer/dev

Browse Source 
Optimisation HTC matrix elements
This commit is contained in:
Anthony Scemama 2023-01-23 16:21:37 +01:00 committed by GitHub
commit 600ef80784
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 1693 additions and 167 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/bi_ort_ints/three_body_ijm.irp.f
View File

@ -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) = <mji|-L|jim>
! three_e_3_idx_cycle_1_bi_ort(m,j,i) = <mji|-L|jim>
!
! 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
!

4
src/bi_ort_ints/three_body_ijmk.irp.f
View File

@ -195,7 +195,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num,
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch13_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_4_idx_exch13_bi_ort(m,j,k,i) = <mjk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
END_DOC
@ -241,7 +241,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num,
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch12_bi_ort(m,j,k,i) = <mjk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_4_idx_exch12_bi_ort(m,j,k,i) = <mjk|-L|mij> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
!

6
src/bi_ort_ints/three_body_ijmkl.irp.f
View File

@ -7,7 +7,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
END_DOC
@ -202,7 +202,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num,
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
!
@ -251,7 +251,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num,
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|mij> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
!

49
src/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -199,3 +199,52 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num) ]
implicit none
BEGIN_DOC
! mo_bi_ortho_tc_two_e_jj(i,j) = J_ij = <ji|W-K|ji>
! mo_bi_ortho_tc_two_e_jj_exchange(i,j) = K_ij = <ij|W-K|ji>
! mo_bi_ortho_tc_two_e_jj_anti(i,j) = J_ij - K_ij
END_DOC
integer :: i,j
double precision :: get_two_e_integral
mo_bi_ortho_tc_two_e_jj = 0.d0
mo_bi_ortho_tc_two_e_jj_exchange = 0.d0
do i=1,mo_num
do j=1,mo_num
mo_bi_ortho_tc_two_e_jj(i,j) = mo_bi_ortho_tc_two_e(j,i,j,i)
mo_bi_ortho_tc_two_e_jj_exchange(i,j) = mo_bi_ortho_tc_two_e(i,j,j,i)
mo_bi_ortho_tc_two_e_jj_anti(i,j) = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)]
&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)]
implicit none
BEGIN_DOC
! tc_2e_3idx_coulomb_integrals(j,k,i) = <jk|ji>
!
! tc_2e_3idx_exchange_integrals(j,k,i) = <kj|ji>
END_DOC
integer :: i,j,k,l
double precision :: get_two_e_integral
double precision :: integral
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i )
tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i )
enddo
enddo
enddo
END_PROVIDER

2
src/determinants/single_excitations.irp.f
View File

@ -28,7 +28,7 @@ BEGIN_PROVIDER [double precision, fock_operator_closed_shell_ref_bitmask, (mo_nu
integer :: occ_virt(N_int*bit_kind_size,2)
integer(bit_kind) :: key_test(N_int)
integer(bit_kind) :: key_virt(N_int,2)
fock_operator_closed_shell_ref_bitmask = 0.d0
call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
do i = 1, N_int
key_virt(i,1) = full_ijkl_bitmask(i)

8
src/determinants/slater_rules.irp.f
View File

@ -1811,12 +1811,12 @@ double precision function diag_H_mat_elem(det_in,Nint)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
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, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta
det_tmp = ref_bitmask
do ispin=1,2

3
src/tc_bi_ortho/slater_tc.irp.f
View File

@ -324,6 +324,9 @@ subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
! if(h1==14.and.p1==2)then
! print*,'h1,p1 old = ',h1,p1
! endif
hmono = mo_bi_ortho_tc_one_e(p1,h1) * phase

37
src/tc_bi_ortho/slater_tc_3e.irp.f
View File

@ -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
@ -269,20 +256,16 @@ subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
if(Ne(1)+Ne(2).ge.3)then
if(s1==s2)then ! same spin excitation
ispin = other_spin(s1)
! print*,'htilde ij'
do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
mm = occ(m,ispin)
!! direct_int = three_body_ints_bi_ort(mm,p2,p1,mm,h2,h1)
!! exchange_int = three_body_ints_bi_ort(mm,p2,p1,mm,h1,h2)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
! print*,direct_int,exchange_int
hthree += direct_int - exchange_int
enddo
do m = 1, Ne(s1) ! pure contribution from s1
mm = occ(m,s1)
hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1)
enddo
do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1)
mm = occ(m,ispin)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
hthree += direct_int - exchange_int
enddo
do m = 1, Ne(s1) ! pure contribution from s1
mm = occ(m,s1)
hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1)
enddo
else ! different spin excitation
do m = 1, Ne(s1)
mm = occ(m,s1) !

44
src/tc_bi_ortho/slater_tc_opt.irp.f Normal file
View File

@ -0,0 +1,44 @@
subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
BEGIN_DOC
!
! <key_j | H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element in terms of single, two and three electron contribution.
!! WARNING !!
!
! Non hermitian !!
!
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree
hmono = 0.d0
htwoe = 0.d0
htot = 0.d0
hthree = 0.D0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2) return
if(degree == 0)then
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)
else if(degree == 2)then
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
endif
if(degree==0) then
htot += nuclear_repulsion
endif
end
! ---

279
src/tc_bi_ortho/slater_tc_opt_diag.irp.f Normal file
View File

@ -0,0 +1,279 @@
BEGIN_PROVIDER [ double precision, ref_tc_energy_tot]
&BEGIN_PROVIDER [ double precision, ref_tc_energy_1e]
&BEGIN_PROVIDER [ double precision, ref_tc_energy_2e]
&BEGIN_PROVIDER [ double precision, ref_tc_energy_3e]
implicit none
BEGIN_DOC
! Various component of the TC energy for the reference "HF" Slater determinant
END_DOC
double precision :: hmono, htwoe, htot, hthree
call diag_htilde_mu_mat_bi_ortho(N_int,HF_bitmask , hmono, htwoe, htot)
ref_tc_energy_1e = hmono
ref_tc_energy_2e = htwoe
if(three_body_h_tc)then
call diag_htilde_three_body_ints_bi_ort(N_int, HF_bitmask, hthree)
ref_tc_energy_3e = hthree
else
ref_tc_energy_3e = 0.d0
endif
ref_tc_energy_tot = ref_tc_energy_1e + ref_tc_energy_2e + ref_tc_energy_3e
END_PROVIDER
subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree, htot)
implicit none
BEGIN_DOC
! Computes $\langle i|H|i \rangle$.
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot,hthree
integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: particle(Nint,2)
integer :: i, nexc(2), ispin
integer :: occ_particle(Nint*bit_kind_size,2)
integer :: occ_hole(Nint*bit_kind_size,2)
integer(bit_kind) :: det_tmp(Nint,2)
integer :: na, nb
ASSERT (Nint > 0)
ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num)
ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num)
nexc(1) = 0
nexc(2) = 0
do i=1,Nint
hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1))
hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2))
particle(i,1) = iand(hole(i,1),det_in(i,1))
particle(i,2) = iand(hole(i,2),det_in(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
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
!call debug_det(det_in,Nint)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
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, Nint)
ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta
det_tmp = ref_bitmask
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)
do i=1,nexc(ispin)
!DIR$ FORCEINLINE
call ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, hmono,htwoe,hthree, Nint,na,nb)
!DIR$ FORCEINLINE
call a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, hmono,htwoe,hthree, Nint,na,nb)
enddo
enddo
htot = hmono+htwoe+hthree
end
subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes one- and two-body energy corresponding
!
! to the ADDITION of an electron in an orbital 'iorb' of spin 'ispin'
!
! onto a determinant 'key'.
!
! in output, the determinant key is changed by the ADDITION of that electron
!
! and the quantities hmono,htwoe,hthree are INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hmono,htwoe,hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,mm,j,m
double precision :: direct_int, exchange_int
if (iorb < 1) then
print *, irp_here, ': iorb < 1'
print *, iorb, mo_num
stop -1
endif
if (iorb > mo_num) then
print *, irp_here, ': iorb > mo_num'
print *, iorb, mo_num
stop -1
endif
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
ASSERT (tmp(1) == elec_alpha_num)
ASSERT (tmp(2) == elec_beta_num)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k >0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb)
! Same spin
do i=1,na
htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
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_prov(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
subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes one- and two-body energy corresponding
!
! to the REMOVAL of an electron in an orbital 'iorb' of spin 'ispin'
!
! onto a determinant 'key'.
!
! in output, the determinant key is changed by the REMOVAL of that electron
!
! and the quantities hmono,htwoe,hthree are INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint
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
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,mm,j,m
integer :: tmp(2)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k>0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
key(k,ispin) = ibclr(key(k,ispin),l)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
na = na-1
hmono = hmono - mo_bi_ortho_tc_one_e(iorb,iorb)
! Same spin
do i=1,na
htwoe= htwoe- mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
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_prov(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

421
src/tc_bi_ortho/slater_tc_opt_double.irp.f Normal file
View File

@ -0,0 +1,421 @@
subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
BEGIN_DOC
! <key_j | H_tilde | key_i> 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
integer :: degree_i,degree_j
call get_excitation_degree(ref_bitmask,key_i,degree_i,N_int)
call get_excitation_degree(ref_bitmask,key_j,degree_j,N_int)
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
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(three_body_h_tc)then
if(.not.double_normal_ord)then
if(degree_i>degree_j)then
call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
else
call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
endif
elseif(double_normal_ord.and.elec_num+elec_num.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
if(degree_i>degree_j)then
call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree)
else
call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
endif
elseif(double_normal_ord.and.elec_num+elec_num.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
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_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(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_prov(ipart,p2,h2,p1,h1)
ihole=occ_hole(i,ispin)
hthree -= three_e_double_parrallel_spin_prov(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
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/beta double excitations
!
! from contraction with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_beta a_h2_beta a_h1_alpha
END_DOC
integer :: i,h1,p1,h2,p2
integer :: hh1,hh2,pp1,pp2,m,mm
integer :: Ne(2)
integer, allocatable :: occ(:,:)
double precision :: contrib
allocate( occ(N_int*bit_kind_size,2) )
call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
call give_contrib_for_abab(1,1,1,1,occ,Ne,contrib)
eff_2_e_from_3_e_ab = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_ab)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
do hh2 = 1, n_act_orb !! beta
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1)
do pp2 = 1, n_act_orb !! beta
p2 = list_act(pp2)
call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
eff_2_e_from_3_e_ab(p2,p1,h2,h1) = contrib
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_beta
!
! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == alpha
!! h2,p2 == beta
contrib = 0.d0
do mm = 1, Ne(1) !! alpha
m = occ(mm,1)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and m
exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int
enddo
do mm = 1, Ne(2) !! beta
m = occ(mm,2)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h2,p2) and m
exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int
enddo
end
BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_aa, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/alpha double excitations
!
! from contractionelec_alpha_num with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_alpha a_h2_alpha a_h1_alpha
!
! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill
!
! |||| h2>h1, p2>p1 ||||
END_DOC
integer :: i,h1,p1,h2,p2
integer :: hh1,hh2,pp1,pp2,m,mm
integer :: Ne(2)
integer, allocatable :: occ(:,:)
double precision :: contrib
allocate( occ(N_int*bit_kind_size,2) )
call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
call give_contrib_for_aaaa(1 ,1 ,1 ,1 ,occ,Ne,contrib)
eff_2_e_from_3_e_aa = 100000000.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_aa)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! alpha
h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! alpha
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! alpha
p1 = list_act(pp1)
do pp2 = pp1+1, n_act_orb !! alpha
p2 = list_act(pp2)
call give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
eff_2_e_from_3_e_aa(p2,p1,h2,h1) = contrib
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
subroutine give_contrib_for_aaaa(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_alpha (h2,p2)_alpha
!
! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == alpha
!! h2,p2 == alpha
contrib = 0.d0
do mm = 1, Ne(1) !! alpha ==> pure parallele spin contribution
m = occ(mm,1)
contrib += three_e_double_parrallel_spin_prov(m,p2,h2,p1,h1)
enddo
do mm = 1, Ne(2) !! beta
m = occ(mm,2)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and (h2,p2)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int
enddo
end
BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_bb, (mo_num, mo_num, mo_num, mo_num)]
implicit none
BEGIN_DOC
! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for beta/beta double excitations
!
! from contractionelec_beta_num with HF density = a^{dagger}_p1_beta a^{dagger}_p2_beta a_h2_beta a_h1_beta
!
! WARNING :: to be coherent with the phase convention used in the Hamiltonian matrix elements, you must fulfill
!
! |||| h2>h1, p2>p1 ||||
END_DOC
integer :: i,h1,p1,h2,p2
integer :: hh1,hh2,pp1,pp2,m,mm
integer :: Ne(2)
integer, allocatable :: occ(:,:)
double precision :: contrib
allocate( occ(N_int*bit_kind_size,2) )
call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
call give_contrib_for_bbbb(1,1 ,1 ,1 ,occ,Ne,contrib)
eff_2_e_from_3_e_bb = 100000000.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, contrib) &
!$OMP SHARED (n_act_orb, list_act, Ne,occ, eff_2_e_from_3_e_bb)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb !! beta
h1 = list_act(hh1)
do hh2 = hh1+1, n_act_orb !! beta
h2 = list_act(hh2)
do pp1 = 1, n_act_orb !! beta
p1 = list_act(pp1)
do pp2 = pp1+1, n_act_orb !! beta
p2 = list_act(pp2)
call give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
eff_2_e_from_3_e_bb(p2,p1,h2,h1) = contrib
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
subroutine give_contrib_for_bbbb(h1,h2,p1,p2,occ,Ne,contrib)
implicit none
BEGIN_DOC
! gives the contribution for a double excitation (h1,p1)_beta (h2,p2)_beta
!
! on top of a determinant whose occupied orbitals is in (occ, Ne)
END_DOC
integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
double precision, intent(out) :: contrib
integer :: mm,m
double precision :: direct_int, exchange_int
!! h1,p1 == beta
!! h2,p2 == beta
contrib = 0.d0
do mm = 1, Ne(2) !! beta ==> pure parallele spin contribution
m = occ(mm,1)
contrib += three_e_double_parrallel_spin_prov(m,p2,h2,p1,h1)
enddo
do mm = 1, Ne(1) !! alpha
m = occ(mm,1)
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
! exchange between (h1,p1) and (h2,p2)
exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1)
contrib += direct_int - exchange_int
enddo
end

460
src/tc_bi_ortho/slater_tc_opt_single.irp.f Normal file
View File

@ -0,0 +1,460 @@
subroutine single_htilde_mu_mat_fock_bi_ortho (Nint, key_j, key_i, hmono, htwoe, hthree, htot)
BEGIN_DOC
! <key_j | H_tilde | key_i> for single 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
double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13
integer :: other_spin(2)
integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2)
other_spin(1) = 2
other_spin(2) = 1
hmono = 0.d0
htwoe = 0.d0
hthree = 0.d0
htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.ne.1)then
return
endif
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
call get_single_excitation_from_fock_tc(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,hthree,htot)
end
subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,htwoe,hthree,htot)
use bitmasks
implicit none
integer,intent(in) :: h,p,spin
double precision, intent(in) :: phase
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
double precision, intent(out) :: hmono,htwoe,hthree,htot
integer(bit_kind) :: differences(N_int,2)
integer(bit_kind) :: hole(N_int,2)
integer(bit_kind) :: partcl(N_int,2)
integer :: occ_hole(N_int*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i
double precision :: buffer_c(mo_num),buffer_x(mo_num)
do i=1, mo_num
buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
enddo
do i = 1, N_int
differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
partcl(i,1) = iand(differences(i,1),key_i(i,1))
partcl(i,2) = iand(differences(i,2),key_i(i,2))
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
hmono = mo_bi_ortho_tc_one_e(p,h)
htwoe = fock_op_2_e_tc_closed_shell(p,h)
! holes :: direct terms
do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,1)
htwoe -= buffer_c(i)
enddo
do i0 = 1, n_occ_ab_hole(2)
i = occ_hole(i0,2)
htwoe -= buffer_c(i)
enddo
! holes :: exchange terms
do i0 = 1, n_occ_ab_hole(spin)
i = occ_hole(i0,spin)
htwoe += buffer_x(i)
enddo
! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,1)
htwoe += buffer_c(i)
enddo
do i0 = 1, n_occ_ab_partcl(2)
i = occ_partcl(i0,2)
htwoe += buffer_c(i)
enddo
! particles :: exchange terms
do i0 = 1, n_occ_ab_partcl(spin)
i = occ_partcl(i0,spin)
htwoe -= buffer_x(i)
enddo
hthree = 0.d0
if (three_body_h_tc)then
call three_comp_fock_elem(key_i,h,p,spin,hthree)
endif
htwoe = htwoe * phase
hmono = hmono * phase
hthree = hthree * phase
htot = htwoe + hmono + hthree
end
subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
implicit none
integer,intent(in) :: h_fock,p_fock,ispin_fock
integer(bit_kind), intent(in) :: key_i(N_int,2)
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)
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
!! Initialize the matrix element with the reference ROHF Slater determinant Fock element
if(ispin_fock==1)then
hthree = fock_a_tot_3e_bi_orth(p_fock,h_fock)
else
hthree = fock_b_tot_3e_bi_orth(p_fock,h_fock)
endif
det_tmp = ref_bitmask
do ispin=1,2
na = elec_num_tab(ispin)
nb = elec_num_tab(iand(ispin,1)+1)
do i=1,nexc(ispin)
!DIR$ FORCEINLINE
call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
!DIR$ FORCEINLINE
call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
enddo
enddo
end
subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes the contribution to the three-electron part of the Fock operator
!
! a^dagger_{p_fock} a_{h_fock} of spin ispin_fock
!
! on top of a determinant 'key' on which you ADD an electron of spin ispin in orbital iorb
!
! in output, the determinant key is changed by the ADDITION of that electron
!
! the output hthree is INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint, h_fock,p_fock, ispin_fock
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,j
double precision :: direct_int, exchange_int
if (iorb < 1) then
print *, irp_here, ': iorb < 1'
print *, iorb, mo_num
stop -1
endif
if (iorb > mo_num) then
print *, irp_here, ': iorb > mo_num'
print *, iorb, mo_num
stop -1
endif
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
ASSERT (tmp(1) == elec_alpha_num)
ASSERT (tmp(2) == elec_beta_num)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k >0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
ASSERT (l >= 0)
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
!! spin of other electrons == ispin
if(ispin == ispin_fock)then
!! in what follows :: jj == other electrons in the determinant
!! :: iorb == electron that has been added of spin ispin
!! :: p_fock, h_fock == hole particle of spin ispin_fock
!! jj = ispin = ispin_fock >> pure parallel spin
do j = 1, na
jj = occ(j,ispin)
hthree += three_e_single_parrallel_spin_prov(jj,iorb,p_fock,h_fock)
enddo
!! spin of jj == other spin than ispin AND ispin_fock
!! exchange between the iorb and (h_fock, p_fock)
do j = 1, nb
jj = occ(j,other_spin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch12_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree += direct_int - exchange_int
enddo
else !! ispin NE to ispin_fock
!! jj = ispin BUT NON EQUAL TO ispin_fock
!! exchange between the jj and iorb
do j = 1, na
jj = occ(j,ispin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch23_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree += direct_int - exchange_int
enddo
!! jj = other_spin than ispin BUT jj == ispin_fock
!! exchange between jj and (h_fock,p_fock)
do j = 1, nb
jj = occ(j,other_spin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch13_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree += direct_int - exchange_int
enddo
endif
na = na+1
end
subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb)
use bitmasks
implicit none
BEGIN_DOC
! Routine that computes the contribution to the three-electron part of the Fock operator
!
! a^dagger_{p_fock} a_{h_fock} of spin ispin_fock
!
! on top of a determinant 'key' on which you REMOVE an electron of spin ispin in orbital iorb
!
! in output, the determinant key is changed by the REMOVAL of that electron
!
! the output hthree is INCREMENTED
END_DOC
integer, intent(in) :: iorb, ispin, Nint, h_fock,p_fock, ispin_fock
integer, intent(inout) :: na, nb
integer(bit_kind), intent(inout) :: key(Nint,2)
double precision, intent(inout) :: hthree
double precision :: direct_int, exchange_int
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i,jj,mm,j,m
integer :: tmp(2)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
k = shiftr(iorb-1,bit_kind_shift)+1
ASSERT (k>0)
l = iorb - shiftl(k-1,bit_kind_shift)-1
key(k,ispin) = ibclr(key(k,ispin),l)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
na = na-1
!! spin of other electrons == ispin
if(ispin == ispin_fock)then
!! in what follows :: jj == other electrons in the determinant
!! :: iorb == electron that has been added of spin ispin
!! :: p_fock, h_fock == hole particle of spin ispin_fock
!! jj = ispin = ispin_fock >> pure parallel spin
do j = 1, na
jj = occ(j,ispin)
hthree -= three_e_single_parrallel_spin_prov(jj,iorb,p_fock,h_fock)
enddo
!! spin of jj == other spin than ispin AND ispin_fock
!! exchange between the iorb and (h_fock, p_fock)
do j = 1, nb
jj = occ(j,other_spin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch12_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree -= direct_int - exchange_int
enddo
else !! ispin NE to ispin_fock
!! jj = ispin BUT NON EQUAL TO ispin_fock
!! exchange between the jj and iorb
do j = 1, na
jj = occ(j,ispin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch23_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree -= direct_int - exchange_int
enddo
!! jj = other_spin than ispin BUT jj == ispin_fock
!! exchange between jj and (h_fock,p_fock)
do j = 1, nb
jj = occ(j,other_spin)
direct_int = three_e_4_idx_direct_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
exchange_int = three_e_4_idx_exch13_bi_ort(jj,iorb,p_fock,h_fock) ! USES 4-IDX TENSOR
hthree -= direct_int - exchange_int
enddo
endif
end
BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ]
implicit none
BEGIN_DOC
! Closed-shell part of the Fock operator for the TC operator
END_DOC
integer :: h0,p0,h,p,k0,k,i
integer :: n_occ_ab(2)
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab_virt(2)
integer :: occ_virt(N_int*bit_kind_size,2)
integer(bit_kind) :: key_test(N_int)
integer(bit_kind) :: key_virt(N_int,2)
double precision :: accu
fock_op_2_e_tc_closed_shell = -1000.d0
call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
do i = 1, N_int
key_virt(i,1) = full_ijkl_bitmask(i)
key_virt(i,2) = full_ijkl_bitmask(i)
key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
! docc ---> virt single excitations
do h0 = 1, n_occ_ab(1)
h=occ(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab(1)
p=occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
! virt ---> virt single excitations
do h0 = 1, n_occ_ab_virt(1)
h=occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p = occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1)
p=occ_virt(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
! docc ---> docc single excitations
do h0 = 1, n_occ_ab(1)
h=occ(h0,1)
do p0 = 1, n_occ_ab(1)
p = occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
do h0 = 1, n_occ_ab(1)
h = occ(h0,1)
do p0 = 1, n_occ_ab(1)
p=occ(p0,1)
accu = 0.d0
do k0 = 1, n_occ_ab(1)
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo
! do i = 1, mo_num
! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i)
! enddo
END_PROVIDER

140
src/tc_bi_ortho/symmetrized_3_e_int_prov.irp.f Normal file
View File

@ -0,0 +1,140 @@
BEGIN_PROVIDER [ double precision, three_e_diag_parrallel_spin_prov, (mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator ON A BI ORTHONORMAL BASIS
!
! three_e_diag_parrallel_spin_prov(m,j,i) = All combinations of the form <mji|-L|mji> for same spin matrix elements
!
! notice the -1 sign: in this way three_e_diag_parrallel_spin_prov can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
integer :: i, j, m
double precision :: integral, wall1, wall0, three_e_diag_parrallel_spin
three_e_diag_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_diag_parrallel_spin_prov ...'
integral = three_e_diag_parrallel_spin(1,1,1) ! to provide all stuffs
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
!$OMP SHARED (mo_num,three_e_diag_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do j = 1, mo_num
do m = j, mo_num
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin(m,j,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do i = 1, mo_num
do j = 1, mo_num
do m = 1, j
three_e_diag_parrallel_spin_prov(m,j,i) = three_e_diag_parrallel_spin_prov(j,m,i)
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for three_e_diag_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [ double precision, three_e_single_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_single_parrallel_spin_prov(m,j,k,i) = All combination of <mjk|-L|mji> for same spin matrix elements
!
! 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
!
END_DOC
implicit none
integer :: i, j, k, m
double precision :: integral, wall1, wall0, three_e_single_parrallel_spin
three_e_single_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_single_parrallel_spin_prov ...'
integral = three_e_single_parrallel_spin(1,1,1,1)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,integral) &
!$OMP SHARED (mo_num,three_e_single_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_single_parrallel_spin_prov(m,j,k,i) = three_e_single_parrallel_spin(m,j,k,i)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_single_parrallel_spin_prov', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_double_parrallel_spin_prov, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_double_parrallel_spin_prov(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! 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
END_DOC
implicit none
integer :: i, j, k, m, l
double precision :: integral, wall1, wall0, three_e_double_parrallel_spin
three_e_double_parrallel_spin_prov = 0.d0
print *, ' Providing the three_e_double_parrallel_spin_prov ...'
call wall_time(wall0)
integral = three_e_double_parrallel_spin(1,1,1,1,1)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_double_parrallel_spin_prov)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
three_e_double_parrallel_spin_prov(m,l,j,k,i) = three_e_double_parrallel_spin(m,l,j,k,i)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_double_parrallel_spin_prov', wall1 - wall0
END_PROVIDER

102
src/tc_bi_ortho/test_normal_order.irp.f
View File

@ -10,6 +10,7 @@ program test_normal_order
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call provide_all_three_ints_bi_ortho
call test
end
@ -28,7 +29,7 @@ subroutine test
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_beta_num+1, mo_num
do p1 = elec_alpha_num+1, mo_num
do h2 = 1, elec_beta_num
do p2 = elec_beta_num+1, mo_num
det_i = ref_bitmask
@ -38,36 +39,93 @@ subroutine test
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
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 = normal_two_body_bi_orth_ab(p2,h2,p1,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
print*,'accu opposite spin = ',accu
stop
s1 = 2
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_beta_num+1, mo_num
do h2 = h1+1, elec_beta_num
do p2 = elec_beta_num+1, mo_num
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
hthree *= phase
normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
accu += dabs(hthree-normal)
enddo
! p2=6
! p1=5
! h2=2
! h1=1
s1 = 1
s2 = 1
accu = 0.d0
do h1 = 1, elec_alpha_num
do p1 = elec_alpha_num+1, mo_num
do p2 = p1+1, mo_num
do h2 = h1+1, elec_alpha_num
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
integer :: hh1, pp1, hh2, pp2, ss1, ss2
call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
hthree *= phase
! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
if(dabs(hthree).lt.1.d-10)cycle
if(dabs(hthree-normal).gt.1.d-10)then
print*,pp2,pp1,hh2,hh1
print*,p2,p1,h2,h1
print*,hthree,normal,dabs(hthree-normal)
stop
endif
! print*,hthree,normal,dabs(hthree-normal)
accu += dabs(hthree-normal)
enddo
enddo
enddo
print*,'accu same spin = ',accu
enddo
print*,'accu same spin alpha = ',accu
s1 = 2
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_beta_num+1, mo_num
do p2 = p1+1, mo_num
do h2 = h1+1, elec_beta_num
det_i = ref_bitmask
call do_single_excitation(det_i,h1,p1,s1,i_ok)
if(i_ok.ne.1)cycle
call do_single_excitation(det_i,h2,p2,s2,i_ok)
if(i_ok.ne.1)cycle
call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call get_excitation_degree(ref_bitmask,det_i,degree,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
hthree *= phase
! normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
normal = eff_2_e_from_3_e_bb(p2,p1,h2,h1)
if(dabs(hthree).lt.1.d-10)cycle
if(dabs(hthree-normal).gt.1.d-10)then
print*,pp2,pp1,hh2,hh1
print*,p2,p1,h2,h1
print*,hthree,normal,dabs(hthree-normal)
stop
endif
! print*,hthree,normal,dabs(hthree-normal)
accu += dabs(hthree-normal)
enddo
enddo
enddo
enddo
print*,'accu same spin beta = ',accu
end

303
src/tc_bi_ortho/test_tc_bi_ortho.irp.f
View File

@ -11,121 +11,210 @@ program tc_bi_ortho
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
! call routine_2
call test_rout
! call test_slater_tc_opt
call timing_tot
! call timing_diag
! call timing_single
! call timing_double
end
subroutine test_rout
subroutine test_slater_tc_opt
implicit none
integer :: i,j,ii,jj
use bitmasks ! you need to include the bitmasks_module.f90 features
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
det_i(:,:)= psi_det(:,:,1)
call debug_det(det_i,N_int)
integer, allocatable :: occ(:,:)
integer :: n_occ_ab(2)
allocate(occ(N_int*bit_kind_size,2))
call bitstring_to_list_ab(det_i, occ, n_occ_ab, N_int)
double precision :: hmono, htwoe, htot
call diag_htilde_mu_mat_bi_ortho(N_int, det_i, hmono, htwoe, htot)
print*,'hmono, htwoe, htot'
print*, hmono, htwoe, htot
print*,'alpha electrons orbital occupancy'
do i = 1, n_occ_ab(1) ! browsing the alpha electrons
j = occ(i,1)
print*,j,mo_bi_ortho_tc_one_e(j,j)
enddo
print*,'beta electrons orbital occupancy'
do i = 1, n_occ_ab(2) ! browsing the beta electrons
j = occ(i,2)
print*,j,mo_bi_ortho_tc_one_e(j,j)
enddo
print*,'alpha beta'
do i = 1, n_occ_ab(1)
ii = occ(i,1)
do j = 1, n_occ_ab(2)
jj = occ(j,2)
print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
print*,'alpha alpha'
do i = 1, n_occ_ab(1)
ii = occ(i,1)
do j = 1, n_occ_ab(1)
jj = occ(j,1)
print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
print*,'beta beta'
do i = 1, n_occ_ab(2)
ii = occ(i,2)
do j = 1, n_occ_ab(2)
jj = occ(j,2)
print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
end
subroutine routine_2
implicit none
integer :: i
double precision :: bi_ortho_mo_ints
print*,'H matrix'
integer :: i,j,degree
double precision :: hmono, htwoe, htot, hthree
double precision :: hnewmono, hnewtwoe, hnewthree, hnewtot
double precision :: accu_d ,i_count, accu
accu = 0.d0
accu_d = 0.d0
i_count = 0.d0
do i = 1, N_det
write(*,'(1000(F16.5,X))')htilde_matrix_elmt_bi_ortho(:,i)
enddo
i = 1
double precision :: phase
integer :: degree,h1, p1, h2, p2, s1, s2, exc(0:2,2,2)
call get_excitation_degree(ref_bitmask, psi_det(1,1,i), degree, N_int)
if(degree==2)then
call get_double_excitation(ref_bitmask, psi_det(1,1,i), exc, phase, N_int)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
print*,'h1,h2,p1,p2'
print*, h1,h2,p1,p2
print*,mo_bi_ortho_tc_two_e(p1,p2,h1,h2),mo_bi_ortho_tc_two_e(h1,h2,p1,p2)
endif
print*,'coef'
do i = 1, ao_num
print*,i,mo_l_coef(i,8),mo_r_coef(i,8)
enddo
! print*,'mdlqfmlqgmqglj'
! print*,'mo_bi_ortho_tc_two_e()',mo_bi_ortho_tc_two_e(2,2,3,3)
! print*,'bi_ortho_mo_ints ',bi_ortho_mo_ints(2,2,3,3)
print*,'Overlap'
do i = 1, mo_num
write(*,'(100(F16.10,X))')overlap_bi_ortho(:,i)
do j = 1,N_det
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hnewmono, hnewtwoe, hnewthree, hnewtot)
if(dabs(htot).gt.1.d-15)then
i_count += 1.D0
accu += dabs(htot-hnewtot)
if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then
call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
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
endif
enddo
enddo
print*,'accu = ',accu/i_count
end
subroutine routine
subroutine timing_tot
implicit none
double precision :: hmono,htwoe,hthree,htot
integer(bit_kind), allocatable :: key1(:,:)
integer(bit_kind), allocatable :: key2(:,:)
allocate(key1(N_int,2),key2(N_int,2))
use bitmasks
key1 = ref_bitmask
call htilde_mu_mat_bi_ortho(key1,key1, N_int, hmono,htwoe,hthree,htot)
key2 = key1
integer :: h,p,i_ok
h = 1
p = 8
call do_single_excitation(key2,h,p,1,i_ok)
call debug_det(key2,N_int)
call htilde_mu_mat_bi_ortho(key2,key1, N_int, hmono,htwoe,hthree,htot)
! print*,'fock_matrix_tc_mo_alpha(p,h) = ',fock_matrix_tc_mo_alpha(p,h)
print*,'htot = ',htot
print*,'hmono = ',hmono
print*,'htwoe = ',htwoe
double precision :: bi_ortho_mo_ints
print*,'bi_ortho_mo_ints(1,p,1,h)',bi_ortho_mo_ints(1,p,1,h)
integer :: i,j
double precision :: wall0, wall1
double precision, allocatable :: mat_old(:,:),mat_new(:,:)
double precision :: hmono, htwoe, hthree, htot, i_count
integer :: degree
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot)
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall0)
i_count = 0.d0
do i = 1, N_det
do j = 1, N_det
! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
i_count += 1.d0
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for old hij for total = ',wall1 - wall0
call wall_time(wall0)
i_count = 0.d0
do i = 1, N_det
do j = 1, N_det
! call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
i_count += 1.d0
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for new hij for total = ',wall1 - wall0
call i_H_j(psi_det(1,1,1), psi_det(1,1,2),N_int,htot)
call wall_time(wall0)
i_count = 0.d0
do i = 1, N_det
do j = 1, N_det
call i_H_j(psi_det(1,1,j), psi_det(1,1,i),N_int,htot)
i_count += 1.d0
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for new hij STANDARD = ',wall1 - wall0
end
subroutine timing_diag
implicit none
integer :: i,j
double precision :: wall0, wall1
double precision, allocatable :: mat_old(:,:),mat_new(:,:)
double precision :: hmono, htwoe, hthree, htot, i_count
integer :: degree
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall0)
i_count = 0.d0
do i = 1, N_det
do j = i,i
i_count += 1.d0
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for old hij for diagonal= ',wall1 - wall0
call wall_time(wall0)
i_count = 0.d0
do i = 1, N_det
do j = i,i
i_count += 1.d0
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for new hij for diagonal= ',wall1 - wall0
end
subroutine timing_single
implicit none
integer :: i,j
double precision :: wall0, wall1,accu
double precision, allocatable :: mat_old(:,:),mat_new(:,:)
double precision :: hmono, htwoe, hthree, htot, i_count
integer :: degree
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
i_count = 0.d0
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
if(degree.ne.1)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo
enddo
print*,'i_count = ',i_count
print*,'time for old hij for singles = ',accu
i_count = 0.d0
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
if(degree.ne.1)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo
enddo
print*,'i_count = ',i_count
print*,'time for new hij for singles = ',accu
end
subroutine timing_double
implicit none
integer :: i,j
double precision :: wall0, wall1,accu
double precision, allocatable :: mat_old(:,:),mat_new(:,:)
double precision :: hmono, htwoe, hthree, htot, i_count
integer :: degree
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,1), psi_det(1,1,1), N_int, hmono, htwoe, hthree, htot)
i_count = 0.d0
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
if(degree.ne.2)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo
enddo
print*,'i_count = ',i_count
print*,'time for old hij for doubles = ',accu
i_count = 0.d0
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,j), psi_det(1,1,i),degree,N_int)
if(degree.ne.2)cycle
i_count += 1.d0
call wall_time(wall0)
call htilde_mu_mat_opt_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
call wall_time(wall1)
accu += wall1 - wall0
enddo
enddo
call wall_time(wall1)
print*,'i_count = ',i_count
print*,'time for new hij for doubles = ',accu
end