mirror of
https://github.com/QuantumPackage/qp2.git
synced 2025-01-03 09:05:39 +01:00
beginning to work on double exc with optimization
This commit is contained in:
eginer
parent
f1137bc883
commit
ac2ebda9ce
@ -256,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) !
|
||||
|
||||
50
src/tc_bi_ortho/slater_tc_opt.irp.f
Normal file
50
src/tc_bi_ortho/slater_tc_opt.irp.f
Normal file
@ -0,0 +1,50 @@
|
||||
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_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
|
||||
endif
|
||||
|
||||
htot = hmono + htwoe + hthree
|
||||
if(degree==0) then
|
||||
htot += nuclear_repulsion
|
||||
endif
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
@ -32,6 +32,11 @@ 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
|
||||
@ -40,7 +45,7 @@ subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
|
||||
!! h2,p2 == beta
|
||||
contrib = 0.d0
|
||||
do mm = 1, Ne(1) !! alpha
|
||||
m = occ(m,1)
|
||||
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)
|
||||
@ -48,10 +53,142 @@ subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
|
||||
enddo
|
||||
|
||||
do mm = 1, Ne(2) !! beta
|
||||
m = occ(m,2)
|
||||
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)
|
||||
eff_2_e_from_3_e_aa = 100000000.d0
|
||||
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
|
||||
|
||||
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
|
||||
double precision :: three_e_double_parrallel_spin
|
||||
!! 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(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)
|
||||
eff_2_e_from_3_e_bb = 100000000.d0
|
||||
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
|
||||
|
||||
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
|
||||
double precision :: three_e_double_parrallel_spin
|
||||
!! 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(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
|
||||
|
||||
|
||||
@ -39,7 +39,7 @@ 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)
|
||||
normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
|
||||
accu += dabs(hthree-normal)
|
||||
enddo
|
||||
@ -48,28 +48,82 @@ subroutine test
|
||||
enddo
|
||||
print*,'accu opposite spin = ',accu
|
||||
|
||||
!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
|
||||
! 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)
|
||||
! accu += dabs(hthree-normal)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
!enddo
|
||||
!print*,'accu same spin = ',accu
|
||||
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
|
||||
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
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user