LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-26 07:02:14 +02:00
Code Issues Releases Wiki Activity

second order for 1h2p with mono excitations

Browse Source
This commit is contained in:
Emmanuel Giner 2016-09-18 19:46:13 +02:00
parent c2e04c647c
commit fbb1409b35
5 changed files with 837 additions and 270 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

603
plugins/MRPT_Utils/new_way_second_order_coef.irp.f
View File

@ -19,6 +19,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
double precision :: get_mo_bielec_integral_schwartz
double precision :: active_int(n_act_orb,2)
double precision :: hij,phase
integer :: index_orb_act_mono(N_det,6)
!matrix_2h1p = 0.d0
elec_num_tab_local = 0
@ -47,10 +48,12 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
integer :: idx(0:N_det)
double precision :: delta_e(n_act_orb,2,N_states)
integer :: istate
integer :: index_orb_act_mono(N_det,3)
do idet = 1, N_det
call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
! if(idet == 81)then
! call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx)
! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r)
do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a)
@ -109,82 +112,52 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!! Second order coefficient : interactions between the perturbers throw the active space
do a = 1, n_act_orb
do jspin = 1, 2
do ispin = 1, 2
if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle
! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = iand(perturb_dets(inint,1,a,jspin,ispin),cas_bitmask(inint,1,1))
det_tmp(inint,2) = iand(perturb_dets(inint,2,a,jspin,ispin),cas_bitmask(inint,1,1))
enddo
do istate = 1, N_states
coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0
enddo
do b = 1, n_act_orb
do kspin = jspin , jspin
integer :: degree_scalar
if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle
do inint = 1, N_int
det_tmp_j(inint,1) = iand(perturb_dets(inint,1,b,kspin,ispin),cas_bitmask(inint,1,1))
det_tmp_j(inint,2) = iand(perturb_dets(inint,2,b,kspin,ispin),cas_bitmask(inint,1,1))
enddo
call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int)
if (degree_scalar > 2 .or. degree_scalar == 0)cycle
! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet >
! print*, '**********************'
! integer(bit_kind) :: det_bis(N_int,2)
! call debug_det(det_tmp,N_int)
! call debug_det(det_tmp_j,N_int)
! do inint = 1, N_int
! det_bis(inint,1) = perturb_dets(inint,1,b,kspin,ispin)
! det_bis(inint,2) = perturb_dets(inint,2,b,kspin,ispin)
! enddo
! call debug_det(det_bis,N_int)
call i_H_j_dyall(det_tmp,det_tmp_j,N_int,hij)
do istate = 1, N_states
coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) &
* hij / delta_e(a,jspin,istate)
if(dabs(hij).gt.0.01d0)then
print*,degree_scalar, hij
print*, coef_perturb_from_idet(b,kspin,ispin,istate,1)* hij / delta_e(a,jspin,istate),coef_perturb_from_idet(a,jspin,ispin,istate,1)
endif
enddo
enddo
enddo
enddo
enddo
enddo
do a = 1, n_act_orb
do jspin = 1, 2
do ispin = 1, 2
if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle
do istate = 1, N_states
! print*, coef_perturb_from_idet(a,jspin,ispin,istate,1),coef_perturb_from_idet(a,jspin,ispin,istate,2)
coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(a,jspin,ispin,istate,1)
enddo
enddo
enddo
enddo
! stop
!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a_{b} a^{\dagger}_a | Idet>
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a^{\dagger}_b a_{a} | Idet>
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | K_{ab} | Idet>
integer :: i_hole,i_part
double precision :: hij_test
double precision :: fock_operator_local(n_act_orb,n_act_orb,2)
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,1)) !!! a_{b}
index_orb_act_mono(idx(jdet),3) = 1
else
! Mono beta
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,2)) !!! a_{b}
index_orb_act_mono(idx(jdet),3) = 2
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
i_part = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b}
kspin = 1 !!! kspin
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
else
! Mono beta
i_hole = list_act_reverse(exc(1,1,2)) !!! a_a
i_part = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b}
kspin = 2 !!! kspin
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
else if(degree(jdet)==2)then
call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
! Mono alpha
index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a
index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b}
index_orb_act_mono(idx(jdet),3) = 1
! Mono beta
index_orb_act_mono(idx(jdet),4) = list_act_reverse(exc(1,1,2)) !!! a_a
index_orb_act_mono(idx(jdet),5) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b}
index_orb_act_mono(idx(jdet),6) = 2
endif
else
index_orb_act_mono(idx(jdet),1) = -1
@ -192,56 +165,181 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p)
enddo
integer :: kspin
integer :: corb,i_ok
integer(bit_kind) :: det_tmp_bis(N_int,2)
double precision :: hib , hab , hja
double precision :: delta_e_ab(N_states)
double precision :: hib_test,hja_test,hab_test
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator
! are connected by the presence of the perturbers determinants |det_tmp>
aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb}
borb = index_orb_act_mono(idx(jdet),2) ! a_{borb}
kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation
! the determinants Idet and Jdet interact throw the following operator
! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet >
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS
if(degree(jdet) == 1)then
! ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator
! ! are connected by the presence of the perturbers determinants |det_tmp>
aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb}
borb = index_orb_act_mono(idx(jdet),2) ! a_{borb}
kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
! ! the determinants Idet and Jdet interact throw the following operator
! ! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet >
do jspin = 1, 2
if (jspin .ne. kspin)then
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
do corb = 1, n_act_orb
if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle
! ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{corb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! ! < idet | H | det_tmp > = phase * (ir|cv)
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
hib= phase * active_int(corb,1)
endif
! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp >
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok)
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin)
enddo
double precision :: hja
! you determine the interaction between the excited determinant and the other parent | Jdet >
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet >
! hja = < det_tmp | H | Jdet >
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,2) - active_int(borb,1) )
else
hja = phase * active_int(borb,1)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = (fock_operator_local(aorb,borb,kspin) ) * phase
if(isnan(hab))then
print*, '1'
stop
endif
! < jdet | H | det_tmp_bis > = phase * (ir|cv)
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
hja= phase * (active_int(corb,1))
endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate)
matrix_2h1p(idx(jdet),idet,istate) = matrix_2h1p(idx(jdet),idet,istate) + &
hib / delta_e(corb,jspin,istate) * hab / delta_e_ab(istate) * hja
! ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp )
! ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis)
! ! < det_tmp_bis | H | Jdet >
enddo
enddo ! corb
else
if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
do corb = 1, n_act_orb
if(corb == aorb .or. corb == borb) cycle
if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle
! ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,jspin} a_{iorb,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
hib= phase * active_int(corb,1)
endif
! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp >
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok)
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = fock_operator_local(aorb,borb,kspin) * phase
if(isnan(hab))then
print*, '2'
stop
endif
! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) )
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
hja= phase * (active_int(corb,1))
endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate)
matrix_2h1p(idx(jdet),idet,istate) = matrix_2h1p(idx(jdet),idet,istate) + &
hib / delta_e(corb,jspin,istate) * hab / delta_e_ab(istate) * hja
! ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp )
! ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis)
! ! < det_tmp_bis | H | Jdet >
enddo
enddo ! corb
endif
enddo
enddo
!
else !! Double excitation operators
!
if (index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then !! spin exchange
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
!!! ! first combination of spin :: | det_tmp > = a^{\dagger}_{aorb,beta} | Idet >
jspin = 2
aorb = index_orb_act_mono(idx(jdet),1) ! hole of the alpha electron
borb = index_orb_act_mono(idx(jdet),2) ! particle of the alpha electron
if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
! | det_tmp > = a^{\dagger}_{aorb,beta} | Idet >
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(aorb,1) - active_int(aorb,2))
else
hib= phase * (active_int(aorb,1))
endif
if(hib .ne. perturb_dets_hij(aorb,jspin,ispin))then
print*, 'pb !!'
print*, 'hib .ne. perturb_dets_hij(aorb,jspin,ispin)'
stop
endif
enddo !! ispin
else if(index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),4))then !! closed shell double excitation
else
call get_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,degree_scalar,phase,N_int)
integer :: h1,h2,p1,p2,s1,s2 , degree_scalar
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*, h1,p1,h2,p2,s1,s2
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
print*, idet,idx(jdet)
print*, 'pb !!!!!!!!!!!!!'
call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx)
stop
endif
endif
do istate = 1, N_states
matrix_2h1p(idx(jdet),idet,istate) += hja * coef_perturb_from_idet(aorb,kspin,ispin,istate,2)
enddo
enddo ! ispin
else
! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations
!
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
do ispin = 1, 2
do kspin = 1, 2
if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
do a = 1, n_act_orb ! First active
do istate = 1, N_states
matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate,2) * perturb_dets_hij(a,kspin,ispin)
enddo
enddo
enddo
enddo
!! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations
!!
!! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet >
!!do ispin = 1, 2
!! do kspin = 1, 2
!! if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count
!! do a = 1, n_act_orb ! First active
!! do istate = 1, N_states
!! matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate,2) * perturb_dets_hij(a,kspin,ispin)
!! enddo
!! enddo
!! enddo
!!enddo
!
endif
enddo
@ -261,8 +359,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
use bitmasks
implicit none
double precision , intent(inout) :: matrix_1h2p(N_det,N_det,*)
integer :: i,v,r,a,b
integer :: iorb, vorb, rorb, aorb, borb
integer :: i,v,r,a,b,c
integer :: iorb, vorb, rorb, aorb, borb,corb
integer :: ispin,jspin
integer :: idet,jdet
integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2)
@ -292,6 +390,9 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
double precision :: hja
double precision :: contrib_hij
double precision :: fock_operator_local(n_act_orb,n_act_orb,2)
double precision :: fock_operator_from_core(n_act_orb,n_act_orb)
double precision :: fock_operator_from_virt(n_act_orb,n_act_orb)
double precision :: fock_operator_from_act(n_act_orb,n_act_orb,n_act_orb,2)
accu_contrib = 0.d0
!matrix_1h2p = 0.d0
@ -394,10 +495,11 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a^{\dagger}_b a_{a} | Idet>
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | K_{ab} | Idet>
integer :: i_hole,i_part
double precision :: hij_test
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
! print*, degree(jdet)
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
@ -408,7 +510,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),N_int,hij)
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
else
@ -419,7 +521,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),N_int,hij)
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
@ -439,11 +541,17 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
endif
enddo
integer ::dorb,i_ok
integer(bit_kind) :: det_tmp_bis(N_int,2)
double precision :: hib , hab
double precision :: delta_e_ab(N_states)
double precision :: hib_test,hja_test,hab_test
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
if(degree(jdet) == 1)then
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS
if(degree(jdet) == 1)then
! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator
! are connected by the presence of the perturbers determinants |det_tmp>
aorb = index_orb_act_mono(idx(jdet),1) ! a_{aorb}
@ -454,11 +562,6 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
integer ::corb,dorb,i_ok
integer(bit_kind) :: det_tmp_bis(N_int,2)
double precision :: hib , hab
double precision :: delta_e_ab(N_states)
double precision :: hib_test,hja_test,hab_test
if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count
do jspin = 1, 2
if (jspin .ne. kspin)then
@ -472,18 +575,12 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * (ir|cv)
! call i_H_j(det_tmp,psi_det(1,1,idet),N_int,hib)
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
hib= phase * active_int(corb,1)
endif
! if(hib_test .ne. hib)then
! print*, 'hib_test .ne. hib'
! print*, hib, hib_test
! stop
! endif
! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp >
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok)
@ -491,21 +588,14 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
! call i_H_j(det_tmp_bis,det_tmp,N_int,hab)
hab = fock_operator_local(aorb,borb,kspin) * phase
hab = (fock_operator_local(aorb,borb,kspin) ) * phase
! < jdet | H | det_tmp_bis > = phase * (ir|cv)
! call i_H_j(det_tmp_bis,psi_det(1,1,idx(jdet)),N_int,hja)
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
hja= phase * (active_int(corb,1))
endif
! if(hja_test .ne. hja)then
! print*, 'hja_test .ne. hja'
! print*, hja, hja_test
! stop
! endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate)
matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + &
@ -519,7 +609,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
do corb = 1, n_act_orb
if(corb == aorb .or. corb == borb) cycle
if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,kspin} a_{iorb,ispin} | Idet >
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,jspin} a_{iorb,ispin} | Idet >
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
@ -527,38 +617,25 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) )
! call i_H_j(det_tmp,psi_det(1,1,idet),N_int,hib)
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(corb,1) - active_int(corb,2))
else
hib= phase * active_int(corb,1)
endif
! if(hib_test .ne. hib)then
! print*, 'hib_test .ne. hib jspin == kspin'
! print*, hib, hib_test
! stop
! endif
! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp >
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok)
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = fock_operator_local(aorb,borb,kspin) * phase
! call i_H_j(det_tmp_bis,det_tmp,N_int,hab)
! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) )
! call i_H_j(det_tmp_bis,psi_det(1,1,idx(jdet)),N_int,hja)
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(corb,1) - active_int(corb,2))
else
hja= phase * (active_int(corb,1))
endif
! if(hja_test .ne. hja)then
! print*, 'hja_test .ne. hja'
! print*, hja, hja_test
! stop
! endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate)
matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + &
@ -574,14 +651,114 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
enddo ! ispin
else
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of double excitations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! call debug_det(psi_det(1,1,idet),N_int)
! call debug_det(psi_det(1,1,idx(jdet)),N_int)
! pause
! a^{\dagger}_r a_{i} (ispin)
aorb = index_orb_act_mono(idx(jdet),4) ! hole of a beta electron
borb = index_orb_act_mono(idx(jdet),5) ! propagation of the hole :: mono excitation of alpha spin
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
! ! first combination of spin :: | det_tmp > = a_{aorb,beta} | Idet >
jspin = 2
if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(borb,1) - active_int(borb,2))
else
hib= phase * (active_int(borb,1))
endif
if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then
call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),1,i_ok)
if(i_ok .ne. 1)then
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
print*, aorb, borb
call debug_det(det_tmp,N_int)
stop
endif
else
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),1,i_ok)
endif
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
if (aorb == borb)then
print*, 'iahaha'
stop
endif
hab = fock_operator_local(aorb,borb,1) * phase
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(borb,1) - active_int(borb,2))
else
hja= phase * (active_int(borb,1))
endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(aorb,jspin,istate) + one_anhil_one_creat(borb,aorb,1,1,istate)
matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + &
hib / delta_e(aorb,jspin,istate) * hab / delta_e_ab(istate) * hja
! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp )
! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis)
! < det_tmp_bis | H | Jdet >
enddo !! istate
! ! second combination of spin :: | det_tmp > = a_{aorb,alpha} | Idet >
jspin = 1
if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin)
enddo
call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int)
if(ispin == jspin)then
hib= phase * (active_int(borb,1) - active_int(borb,2))
else
hib= phase * (active_int(borb,1))
endif
if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then
call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),2,i_ok)
if(i_ok .ne. 1)then
call debug_det(psi_det(1,1,idet),N_int)
call debug_det(psi_det(1,1,idx(jdet)),N_int)
print*, aorb, borb
call debug_det(det_tmp,N_int)
stop
endif
else
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),2,i_ok)
endif
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
! < det_tmp | H | det_tmp_bis > = F_{aorb,borb}
hab = fock_operator_local(aorb,borb,2) * phase
call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int)
if(ispin == jspin)then
hja= phase * (active_int(borb,1) - active_int(borb,2))
else
hja= phase * (active_int(borb,1))
endif
do istate = 1, N_states
delta_e_ab(istate) = delta_e(aorb,jspin,istate) + one_anhil_one_creat(borb,aorb,1,1,istate)
matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + &
hib / delta_e(aorb,jspin,istate) * hab / delta_e_ab(istate) * hja
! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp )
! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis)
! < det_tmp_bis | H | Jdet >
enddo !! istate
enddo !! ispin
endif
endif !! en of test if jdet is a single or a double excitation of type K_ab
else
else !! jdet is idet
! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations
!
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet >
@ -603,14 +780,13 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
endif
enddo
enddo !! jdet
enddo
enddo
enddo
enddo
print* , 'accu_contrib = ',accu_contrib
@ -618,104 +794,3 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p)
end
! do a = 1, n_act_orb
! do jspin = 1, 2
! do ispin = 1, 2
! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle
! ! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet >
! do inint = 1, N_int
! det_tmp(inint,1) = perturb_dets(inint,1,a,jspin,ispin)
! det_tmp(inint,2) = perturb_dets(inint,2,a,jspin,ispin)
! enddo
! do istate = 1, N_states
! coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0
! enddo
! do b = 1, n_act_orb
! do kspin = jspin , jspin
! integer :: degree_scalar
! if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle
! do inint = 1, N_int
! det_tmp_j(inint,1) = perturb_dets(inint,1,b,kspin,ispin)
! det_tmp_j(inint,2) = perturb_dets(inint,2,b,kspin,ispin)
! enddo
! call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int)
! if (degree_scalar > 2 .or. degree_scalar == 0)cycle
! ! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet >
! call i_H_j(det_tmp,det_tmp_j,N_int,hij)
! do istate = 1, N_states
! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) &
! * hij / delta_e(a,jspin,istate)
! endif
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
! do a = 1, n_act_orb
! do jspin = 1, 2
! do ispin = 1, 2
! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle
! ! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet >
! do inint = 1, N_int
! det_tmp(inint,1) = iand(perturb_dets(inint,1,a,jspin,ispin),cas_bitmask(inint,1,1))
! det_tmp(inint,2) = iand(perturb_dets(inint,2,a,jspin,ispin),cas_bitmask(inint,1,1))
! enddo
! do istate = 1, N_states
! coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0
! enddo
! do b = 1, n_act_orb
! do kspin = jspin , jspin
! integer :: degree_scalar
! if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle
! do inint = 1, N_int
! det_tmp_j(inint,1) = iand(perturb_dets(inint,1,b,kspin,ispin),cas_bitmask(inint,1,1))
! det_tmp_j(inint,2) = iand(perturb_dets(inint,2,b,kspin,ispin),cas_bitmask(inint,1,1))
! enddo
! call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int)
! if (degree_scalar > 2 .or. degree_scalar == 0)cycle
! ! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet >
!! print*, '**********************'
!! integer(bit_kind) :: det_bis(N_int,2)
!! call debug_det(det_tmp,N_int)
!! call debug_det(det_tmp_j,N_int)
!! do inint = 1, N_int
!! det_bis(inint,1) = perturb_dets(inint,1,b,kspin,ispin)
!! det_bis(inint,2) = perturb_dets(inint,2,b,kspin,ispin)
!! enddo
!! call debug_det(det_bis,N_int)
! call i_H_j_dyall(det_tmp,det_tmp_j,N_int,hij)
! do istate = 1, N_states
! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) &
! * hij / delta_e(a,jspin,istate)
! if(dabs(hij).gt.0.01d0)then
! print*,degree_scalar, hij
! print*, coef_perturb_from_idet(b,kspin,ispin,istate,1)* hij / delta_e(a,jspin,istate),coef_perturb_from_idet(a,jspin,ispin,istate,1)
!
! endif
! enddo
! enddo
! enddo
! enddo
! enddo
! enddo
! do a = 1, n_act_orb
! do jspin = 1, 2
! do ispin = 1, 2
! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle
! do istate = 1, N_states
! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(a,jspin,ispin,istate,1)
! enddo
! enddo
! enddo
! enddo

23
plugins/MRPT_Utils/print_1h2p.irp.f
View File

@ -17,7 +17,8 @@ subroutine routine
enddo
enddo
enddo
call give_1h2p_contrib_sec_order(matrix_1h2p)
if(.False.)then
call give_1h2p_contrib(matrix_1h2p)
double precision :: accu
accu = 0.d0
do i = 1, N_det
@ -25,7 +26,25 @@ subroutine routine
accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
print*, 'accu', accu
print*, 'second order ', accu
endif
if(.True.)then
do i = 1, N_det
do j = 1, N_det
do istate = 1, N_states
matrix_1h2p(i,j,istate) = 0.d0
enddo
enddo
enddo
call give_1h2p_new(matrix_1h2p)
accu = 0.d0
do i = 1, N_det
do j = 1, N_det
accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1)
enddo
enddo
print*, 'third order ', accu
deallocate (matrix_1h2p)
end

303
plugins/MRPT_Utils/second_order_new.irp.f Normal file
View File

@ -0,0 +1,303 @@
subroutine give_1h2p_new(matrix_1h2p)
use bitmasks
implicit none
double precision , intent(inout) :: matrix_1h2p(N_det,N_det,*)
integer :: i,v,r,a,b,c
integer :: iorb, vorb, rorb, aorb, borb,corb
integer :: ispin,jspin
integer :: idet,jdet
integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2)
double precision :: perturb_dets_phase(n_act_orb,2,2)
double precision :: perturb_dets_hij(n_act_orb,2,2)
double precision :: perturb_dets_hpsi0(n_act_orb,2,2,N_states)
double precision :: coef_perturb_from_idet(n_act_orb,2,2,N_states,2)
logical :: already_generated(n_act_orb,2,2)
integer :: inint
integer :: elec_num_tab_local(2),acu_elec
integer(bit_kind) :: det_tmp(N_int,2)
integer(bit_kind) :: det_tmp_j(N_int,2)
integer :: exc(0:2,2,2)
integer :: accu_elec
double precision :: get_mo_bielec_integral_schwartz
double precision :: active_int(n_act_orb,2)
double precision :: hij,phase
double precision :: accu_contrib(N_states)
integer :: degree(N_det)
integer :: idx(0:N_det)
double precision :: delta_e(n_act_orb,2,N_states)
double precision :: delta_e_inv(n_act_orb,2,N_states)
integer :: istate
integer :: index_orb_act_mono(N_det,6)
double precision :: delta_e_inactive_virt(N_states)
integer :: kspin
double precision :: delta_e_ja(N_states)
double precision :: hja
double precision :: contrib_hij
double precision :: fock_operator_local(n_act_orb,n_act_orb,2)
double precision :: hij_test
integer ::i_ok
integer(bit_kind) :: det_tmp_bis(N_int,2)
double precision :: hib , hab
double precision :: delta_e_ab(N_states)
double precision :: hib_test,hja_test,hab_test
integer :: i_hole,i_part
double precision :: hia,hjb
integer :: other_spin(2)
other_spin(1) = 2
other_spin(2) = 1
accu_contrib = 0.d0
!matrix_1h2p = 0.d0
elec_num_tab_local = 0
do inint = 1, N_int
elec_num_tab_local(1) += popcnt(psi_det(inint,1,1))
elec_num_tab_local(2) += popcnt(psi_det(inint,2,1))
enddo
do i = 1, n_inact_orb ! First inactive
iorb = list_inact(i)
do v = 1, n_virt_orb ! First virtual
vorb = list_virt(v)
do r = 1, n_virt_orb ! Second virtual
rorb = list_virt(r)
! take all the integral you will need for i,j,r fixed
do a = 1, n_act_orb
aorb = list_act(a)
active_int(a,1) = get_mo_bielec_integral_schwartz(iorb,aorb,rorb,vorb,mo_integrals_map) ! direct
active_int(a,2) = get_mo_bielec_integral_schwartz(iorb,aorb,vorb,rorb,mo_integrals_map) ! exchange
perturb_dets_phase(a,1,1) = -1000.d0
perturb_dets_phase(a,1,2) = -1000.d0
perturb_dets_phase(a,2,2) = -1000.d0
perturb_dets_phase(a,2,1) = -1000.d0
enddo
do istate = 1, N_states
delta_e_inactive_virt(istate) = &
- fock_virt_total_spin_trace(rorb,istate) &
- fock_virt_total_spin_trace(vorb,istate) &
+ fock_core_inactive_total_spin_trace(iorb,istate)
enddo
do idet = 1, N_det
call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements
do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb)
do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb)
do a = 1, n_act_orb ! First active
aorb = list_act(a)
do istate = 1, N_states
perturb_dets_hpsi0(a,jspin,ispin,istate) = 0.d0
coef_perturb_from_idet(a,jspin,ispin,istate,1) = 0.d0
coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0
enddo
if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count
do inint = 1, N_int
det_tmp(inint,1) = psi_det(inint,1,idet)
det_tmp(inint,2) = psi_det(inint,2,idet)
enddo
! Do the excitation inactive -- > virtual
call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin
call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin
! Do the excitation active -- > virtual
call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin
call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin
! Check if the excitation is possible or not on psi_det(idet)
accu_elec= 0
do inint = 1, N_int
accu_elec+= popcnt(det_tmp(inint,jspin))
enddo
if(accu_elec .ne. elec_num_tab_local(jspin))then
perturb_dets_phase(a,jspin,ispin) = -1000.0d0
perturb_dets_hij(a,jspin,ispin) = 0.d0
do istate = 1, N_states
coef_perturb_from_idet(a,jspin,ispin,istate,1) = 0.d0
coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0
enddo
cycle
endif
do inint = 1, N_int
perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1)
perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
perturb_dets_phase(a,jspin,ispin) = phase
do istate = 1, N_states
delta_e(a,jspin,istate) = one_anhil(a,jspin,istate) + delta_e_inactive_virt(istate)
delta_e_inv(a,jspin,istate) = 1.d0 / delta_e(a,jspin,istate)
enddo
if(ispin == jspin)then
perturb_dets_hij(a,jspin,ispin) = phase * (active_int(a,1) - active_int(a,2) )
else
perturb_dets_hij(a,jspin,ispin) = phase * active_int(a,1)
endif
enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS
!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | a^{\dagger}_b a_{a} | Idet>
!!!!!!!!!!!!!!!!!!!!!!!!!!!! <Jdet | K_{ab} | Idet>
do jdet = 1, idx(0)
if(idx(jdet).ne.idet)then
if(degree(jdet)==1)then
call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Mono alpha
i_hole = list_act_reverse(exc(1,1,1)) !!! a_a
i_part = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b}
kspin = 1 !!! kspin
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
else
! Mono beta
i_hole = list_act_reverse(exc(1,1,2)) !!! a_a
i_part = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b}
kspin = 2 !!! kspin
index_orb_act_mono(idx(jdet),1) = i_hole
index_orb_act_mono(idx(jdet),2) = i_part
index_orb_act_mono(idx(jdet),3) = kspin
call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij)
fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator
fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator
endif
endif
else
index_orb_act_mono(idx(jdet),1) = -1
endif
enddo
do jdet = 1, idx(0)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS
if(degree(jdet) == 1)then
! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator
! are connected by the presence of the perturbers determinants |det_tmp>
aorb = index_orb_act_mono(idx(jdet),1) ! a_{aorb}
borb = index_orb_act_mono(idx(jdet),2) ! a^{\dagger}_{borb}
kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation
! the determinants Idet and Jdet interact throw the following operator
! | Jdet > = a^{\dagger}_{borb,kspin} a_{aorb, kspin} | Idet >
accu_contrib = 0.d0
do ispin = 1, 2 ! you loop on all possible spin for the excitation
! a^{\dagger}_r a_{i} (ispin)
if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count
! FIRST ORDER CONTRIBUTION
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet >
if(perturb_dets_phase(aorb,kspin,ispin) .le. -10.d0)cycle
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin)
enddo
call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hia = phase * (active_int(aorb,1) - active_int(aorb,2) )
else
hia = phase * active_int(aorb,1)
endif
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int)
if(kspin == ispin)then
hja = phase * (active_int(borb,1) - active_int(borb,2) )
else
hja = phase * active_int(borb,1)
endif
call i_H_j(det_tmp,psi_det(1,1,idx(jdet)),N_int,hij_test)
if(hij_test .ne. hja)then
print*, 'hij_test .ne. hja'
stop
endif
contrib_hij = hja * hia
do istate = 1, N_states
accu_contrib(istate) += contrib_hij * delta_e_inv(aorb,kspin,istate)
enddo
!!!! SECOND ORDER CONTRIBTIONS
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,jspin} a_{corb,jspin} a_{iorb,ispin} | Idet >
do jspin = 1, 2
do corb = 1, n_act_orb
if(perturb_dets_phase(corb,jspin,ispin) .le. -10.d0)cycle
do inint = 1, N_int
det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin)
det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin)
enddo
! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet >
call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok)
if(i_ok .ne. 1)cycle
call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int)
hia = perturb_dets_hij(corb,jspin,ispin)
hab = fock_operator_local(aorb,borb,kspin) * phase
call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int)
if(jspin == ispin)then
hjb = phase * (active_int(corb,1) - active_int(corb,2) )
else
hjb = phase * active_int(corb,1)
endif
do istate = 1, N_states
accu_contrib(istate)+=hia * delta_e_inv(corb,jspin,istate) & ! | Idet > --> | det_tmp >
! | det_tmp > --> | det_tmp_bis >
*hab / (delta_e(corb,jspin,istate) + one_anhil_one_creat(aorb,borb,kspin,kspin,istate)) &
*hjb
enddo
enddo
enddo
enddo ! ispin
do istate = 1, N_states
matrix_1h2p(idet,idx(jdet),istate) += accu_contrib(istate)
enddo
else if (degree(jdet) == 0)then
! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations
!
! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet >
accu_contrib = 0.d0
do ispin = 1, 2
do kspin = 1, 2
do a = 1, n_act_orb ! First active
if( perturb_dets_phase(a,kspin,ispin) .le. -10.d0)cycle
if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count
contrib_hij = perturb_dets_hij(a,kspin,ispin) * perturb_dets_hij(a,kspin,ispin)
do istate = 1, N_states
accu_contrib(istate) += contrib_hij * delta_e_inv(a,kspin,istate)
enddo
enddo
enddo
enddo
do istate = 1, N_states
matrix_1h2p(idet,idet,istate) += accu_contrib(istate)
enddo
endif
enddo !! jdet
enddo
enddo
enddo
enddo
end

4
src/Determinants/create_excitations.irp.f
View File

@ -31,8 +31,8 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok)
n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
enddo
if(n_elec_tmp .ne. elec_num)then
print*, n_elec_tmp,elec_num
call debug_det(key_in,N_int)
!print*, n_elec_tmp,elec_num
!call debug_det(key_in,N_int)
i_ok = -1
endif
end

174
src/Determinants/slater_rules.irp.f
View File

@ -1367,8 +1367,178 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
exchange_1 = popcnt(xor(iand(key1(1,1,i),key1(1,2,i)),iand(key2(1,1),key2(1,2))))
exchange_2 = popcnt(iand(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2))))
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2))))
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2))))
if (d > 4)cycle
if (d ==4)then
if(exchange_1 .eq. 0 ) then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else
cycle
endif
! pause
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else if (Nint==2) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2)))
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,2),key2(1,2)))) + &
popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,2),key2(2,2))))
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + &
popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2))))
if (d > 4)cycle
if (d ==4)then
if(exchange_1 .eq. 0 ) then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else
cycle
endif
! pause
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else if (Nint==3) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2))) + &
popcnt(xor( key1(2,1,i), key2(2,1))) + &
popcnt(xor( key1(2,2,i), key2(2,2))) + &
popcnt(xor( key1(3,1,i), key2(3,1))) + &
popcnt(xor( key1(3,2,i), key2(3,2)))
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + &
popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,1),key2(2,2)))) + &
popcnt(xor(ior(key1(3,1,i),key1(3,2,i)),ior(key2(3,1),key2(3,2))))
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + &
popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) + &
popcnt(ior(xor(key1(3,1,i),key2(3,1)),xor(key1(3,2,i),key2(3,2))))
if (d > 4)cycle
if (d ==4)then
if(exchange_1 .eq. 0 ) then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else
cycle
endif
! pause
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = 0
exchange_1 = 0
!DIR$ LOOP COUNT MIN(4)
do m=1,Nint
d = d + popcnt(xor( key1(m,1,i), key2(m,1))) &
+ popcnt(xor( key1(m,2,i), key2(m,2)))
exchange_1 = popcnt(xor(ior(key1(m,1,i),key1(m,2,i)),ior(key2(m,1),key2(m,2))))
exchange_2 = popcnt(ior(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2))))
enddo
if (d > 4)cycle
if (d ==4)then
if(exchange_1 .eq. 0 ) then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
else
cycle
endif
! pause
else
degree(l) = ishft(d,-1)
idx(l) = i
l = l+1
endif
enddo
endif
idx(0) = l-1
end
subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degree,Nint,sze,idx)
use bitmasks
implicit none
BEGIN_DOC
! Applies get_excitation_degree to an array of determinants and return only the mono excitations
! and the connections through exchange integrals
END_DOC
integer, intent(in) :: Nint, sze
integer(bit_kind), intent(in) :: key1(Nint,2,sze)
integer(bit_kind), intent(in) :: key2(Nint,2)
integer, intent(out) :: degree(sze)
integer, intent(out) :: idx(0:sze)
integer :: i,l,d,m
integer :: exchange_1,exchange_2
ASSERT (Nint > 0)
ASSERT (sze > 0)
l=1
if (Nint==1) then
!DIR$ LOOP COUNT (1000)
do i=1,sze
d = popcnt(xor( key1(1,1,i), key2(1,1))) + &
popcnt(xor( key1(1,2,i), key2(1,2)))
exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2))))
exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2))))
if(i==99)then
integer(bit_kind) :: key_test(N_int,2)
key_test(1,2) = 0_bit_kind
call debug_det(key2,N_int)
key_test(1,1) = ior(key2(1,1),key2(1,2))
call debug_det(key_test,N_int)
key_test(1,1) = ior(key1(1,1,i),key1(1,2,i))
call debug_det(key1(1,1,i),N_int)
call debug_det(key_test,N_int)
key_test(1,1) = xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))
call debug_det(key_test,N_int)
print*, exchange_1 , exchange_2
stop
endif
if (d > 4)cycle
if (d ==4)then
if(exchange_1 .eq. 0 ) then