LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

added tc_bi_ortho

Browse Source
This commit is contained in:
eginer 2022-10-05 17:01:27 +02:00
parent 8c51af6c5a
commit 57e592870c
26 changed files with 3067 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/tc_bi_ortho/EZFIO.cfg Normal file
View File

@ -0,0 +1,11 @@
[psi_l_coef_bi_ortho]
interface: ezfio
doc: Coefficients for the left wave function
type: double precision
size: (determinants.n_det,determinants.n_states)
[psi_r_coef_bi_ortho]
interface: ezfio
doc: Coefficients for the right wave function
type: double precision
size: (determinants.n_det,determinants.n_states)

6
src/tc_bi_ortho/NEED Normal file
View File

@ -0,0 +1,6 @@
bi_ort_ints
bi_ortho_mos
tc_keywords
non_hermit_dav
dav_general_mat
tc_scf

104
src/tc_bi_ortho/e_corr_bi_ortho.irp.f Normal file
View File

@ -0,0 +1,104 @@
use bitmasks ! you need to include the bitmasks_module.f90 features
BEGIN_PROVIDER [ double precision, e_tilde_00]
implicit none
double precision :: hmono,htwoe,hthree,htot
call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot)
e_tilde_00 = htot
END_PROVIDER
BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth]
&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_single]
&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_double]
implicit none
integer :: i,degree
double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e
e_pt2_tc_bi_orth = 0.d0
e_pt2_tc_bi_orth_single = 0.d0
e_pt2_tc_bi_orth_double = 0.d0
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
e_pt2_tc_bi_orth += coef_pt1 * htilde_ij
if(degree == 1)then
e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij
else
! print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij
e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij
endif
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00]
implicit none
double precision :: hmono,htwoe,hthree,htilde_ij
call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00)
e_tilde_bi_orth_00 += nuclear_repulsion
END_PROVIDER
BEGIN_PROVIDER [ double precision, e_corr_bi_orth ]
&BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj ]
&BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth ]
&BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth ]
implicit none
integer :: i,degree
double precision :: hmono,htwoe,hthree,htilde_ij
e_corr_bi_orth = 0.d0
e_corr_single_bi_orth = 0.d0
e_corr_double_bi_orth = 0.d0
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
if(degree == 1)then
e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
else if(degree == 2)then
e_corr_double_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
! print*,'coef_wf , e_cor',reigvec_tc_bi_orth(i,1)/reigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)
endif
enddo
e_corr_bi_orth_proj = e_corr_single_bi_orth + e_corr_double_bi_orth
e_corr_bi_orth = eigval_right_tc_bi_orth(1) - e_tilde_bi_orth_00
END_PROVIDER
BEGIN_PROVIDER [ double precision, e_tc_left_right ]
implicit none
integer :: i,j
double precision :: hmono,htwoe,hthree,htilde_ij,accu
e_tc_left_right = 0.d0
accu = 0.d0
do i = 1, N_det
accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1)
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,htilde_ij)
e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1)
enddo
enddo
e_tc_left_right *= 1.d0/accu
e_tc_left_right += nuclear_repulsion
END_PROVIDER
BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)]
implicit none
integer :: i,degree
double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree==0)then
coef_pt1_bi_ortho(i) = 1.d0
else
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
coef_pt1_bi_ortho(i)= coef_pt1
endif
enddo
END_PROVIDER

92
src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f Normal file
View File

@ -0,0 +1,92 @@
subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze)
use bitmasks
BEGIN_DOC
! Application of H_TC on a vector
!
! v(i,istate) = \sum_j u(j,istate) H_TC(i,j), with:
! H_TC(i,j) = < Di | H_TC | Dj >
!
END_DOC
implicit none
integer, intent(in) :: N_st, sze
double precision, intent(in) :: u(sze,N_st)
double precision, intent(inout) :: v(sze,N_st)
integer :: i, j, istate
double precision :: htot
PROVIDE N_int
PROVIDE psi_det
! TODO : transform it with the bi-linear representation in terms of alpha-beta.
i = 1
j = 1
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
!$OMP SHARED(N_st, sze, N_int, psi_det, u, v) &
!$OMP PRIVATE(istate, i, j, htot)
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v(i,istate) = v(i,istate) + htot * u(j,istate)
enddo
enddo
enddo
!$OMP END PARALLEL DO
end
subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze)
use bitmasks
BEGIN_DOC
! Application of (H_TC)^dagger on a vector
!
! v(i,istate) = \sum_j u(j,istate) H_TC(j,i), with:
! H_TC(i,j) = < Di | H_TC | Dj >
!
END_DOC
implicit none
integer, intent(in) :: N_st, sze
double precision, intent(in) :: u(sze,N_st)
double precision, intent(inout) :: v(sze,N_st)
integer :: i, j, istate
double precision :: htot
PROVIDE N_int
PROVIDE psi_det
i = 1
j = 1
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot)
v = 0.d0
!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) &
!$OMP SHARED(N_st, sze, N_int, psi_det, u, v) &
!$OMP PRIVATE(istate, i, j, htot)
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
v(i,istate) = v(i,istate) + htot * u(j,istate)
enddo
enddo
enddo
!$OMP END PARALLEL DO
end

319
src/tc_bi_ortho/normal_ordered.irp.f Normal file
View File

@ -0,0 +1,319 @@
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
! Normal ordering of the three body interaction on the HF density
END_DOC
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i,h1,p1,h2,p2
integer :: hh1,hh2,pp1,pp2
integer :: Ne(2)
integer, allocatable :: occ(:,:)
integer(bit_kind), allocatable :: key_i_core(:,:)
double precision :: hthree_aba,hthree_aaa,hthree_aab
double precision :: wall0,wall1
PROVIDE N_int
allocate( occ(N_int*bit_kind_size,2) )
allocate( key_i_core(N_int,2) )
if(core_tc_op) then
do i = 1, N_int
key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
else
call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
endif
normal_two_body_bi_orth = 0.d0
print*,'Providing normal_two_body_bi_orth ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aba, hthree_aab, hthree_aaa) &
!$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth)
!$OMP DO SCHEDULE (static)
do hh1 = 1, n_act_orb
h1 = list_act(hh1)
do pp1 = 1, n_act_orb
p1 = list_act(pp1)
do hh2 = 1, n_act_orb
h2 = list_act(hh2)
do pp2 = 1, n_act_orb
p2 = list_act(pp2)
! opposite spin double excitations
call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aba)
! same spin double excitations with opposite spin contributions
if(h1<h2.and.p1.gt.p2)then
call give_aab_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aab) ! exchange h1<->h2
! same spin double excitations with same spin contributions
if(Ne(2).ge.3)then
call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
else
hthree_aaa = 0.d0
endif
else
call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
if(Ne(2).ge.3)then
call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
else
hthree_aaa = 0.d0
endif
endif
normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'Wall time for normal_two_body_bi_orth ',wall1-wall0
deallocate( occ )
deallocate( key_i_core )
END_PROVIDER
subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer, intent(in) :: Nint, h1, h2, p1, p2
integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
double precision, intent(out) :: hthree
integer :: ii, i
double precision :: int_direct, int_exc_12, int_exc_13, integral
!!!! double alpha/beta
hthree = 0.d0
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
int_direct = -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
int_exc_13 = -1.d0 * integral
call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
int_exc_12 = -1.d0 * integral
hthree += 2.d0 * int_direct - 1.d0 * ( int_exc_13 + int_exc_12)
enddo
do ii = Ne(2) + 1, Ne(1) ! purely open-shell part
i = occ(ii,1)
call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
int_direct = -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
int_exc_13 = -1.d0 * integral
call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral)
int_exc_12 = -1.d0 * integral
hthree += 1.d0 * int_direct - 0.5d0* ( int_exc_13 + int_exc_12)
enddo
end subroutine give_aba_contraction
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
! Normal ordered two-body sector of the three-body terms for opposite spin double excitations
END_DOC
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: h1, p1, h2, p2, i
integer :: hh1, hh2, pp1, pp2
integer :: Ne(2)
integer, allocatable :: occ(:,:)
integer(bit_kind), allocatable :: key_i_core(:,:)
double precision :: hthree
PROVIDE N_int
allocate( key_i_core(N_int,2) )
allocate( occ(N_int*bit_kind_size,2) )
if(core_tc_op)then
do i = 1, N_int
key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core,occ,Ne,N_int)
else
call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
endif
normal_two_body_bi_orth_ab = 0.d0
do hh1 = 1, n_act_orb
h1 = list_act(hh1)
do pp1 = 1, n_act_orb
p1 = list_act(pp1)
do hh2 = 1, n_act_orb
h2 = list_act(hh2)
do pp2 = 1, n_act_orb
p2 = list_act(pp2)
call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree)
normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree
enddo
enddo
enddo
enddo
deallocate( key_i_core )
deallocate( occ )
END_PROVIDER
BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_act_orb, n_act_orb, n_act_orb)]
BEGIN_DOC
! Normal ordered two-body sector of the three-body terms for same spin double excitations
END_DOC
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i,ii,j,h1,p1,h2,p2
integer :: hh1,hh2,pp1,pp2
integer :: Ne(2)
integer, allocatable :: occ(:,:)
integer(bit_kind), allocatable :: key_i_core(:,:)
double precision :: hthree_aab, hthree_aaa
PROVIDE N_int
allocate( key_i_core(N_int,2) )
allocate( occ(N_int*bit_kind_size,2) )
if(core_tc_op)then
do i = 1, N_int
key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
else
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
endif
normal_two_body_bi_orth_aa_bb = 0.d0
do hh1 = 1, n_act_orb
h1 = list_act(hh1)
do pp1 = 1 , n_act_orb
p1 = list_act(pp1)
do hh2 = 1, n_act_orb
h2 = list_act(hh2)
do pp2 = 1 , n_act_orb
p2 = list_act(pp2)
if(h1<h2.and.p1.gt.p2)then
call give_aab_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aab) ! exchange h1<->h2
if(Ne(2).ge.3)then
call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2
else
hthree_aaa = 0.d0
endif
else
call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab)
if(Ne(2).ge.3)then
call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa)
else
hthree_aaa = 0.d0
endif
endif
normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) = hthree_aab + hthree_aaa
enddo
enddo
enddo
enddo
deallocate( key_i_core )
deallocate( occ )
END_PROVIDER
subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer, intent(in) :: Nint, h1, h2, p1, p2
integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
double precision, intent(out) :: hthree
integer :: ii,i
double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23
double precision :: integral,int_exc_l,int_exc_ll
hthree = 0.d0
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
int_direct = -1.d0 * integral
call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
int_exc_l = -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
int_exc_ll= -1.d0 * integral
call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
int_exc_12= -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
int_exc_13= -1.d0 * integral
call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
int_exc_23= -1.d0 * integral
hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 )
enddo
do ii = Ne(2)+1,Ne(1) ! purely open-shell part
i = occ(ii,1)
call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral)
int_direct = -1.d0 * integral
call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral)
int_exc_l = -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral)
int_exc_ll= -1.d0 * integral
call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral)
int_exc_12= -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral)
int_exc_13= -1.d0 * integral
call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral)
int_exc_23= -1.d0 * integral
hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 ))
enddo
end subroutine give_aaa_contraction
subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
implicit none
use bitmasks ! you need to include the bitmasks_module.f90 features
integer, intent(in) :: Nint, h1, h2, p1, p2
integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
double precision, intent(out) :: hthree
integer :: ii, i
double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23
double precision :: integral, int_exc_l, int_exc_ll
hthree = 0.d0
do ii = 1, Ne(2) ! purely closed shell part
i = occ(ii,2)
call give_integrals_3_body_bi_ort(p2,p1,i,h2,h1,i,integral)
int_direct = -1.d0 * integral
call give_integrals_3_body_bi_ort(p1,p2,i,h2,h1,i,integral)
int_exc_23= -1.d0 * integral
hthree += 1.d0 * int_direct - int_exc_23
enddo
end subroutine give_aab_contraction

104
src/tc_bi_ortho/print_tc_wf.irp.f Normal file
View File

@ -0,0 +1,104 @@
program print_tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
! if(three_body_h_tc)then
! call provide_all_three_ints_bi_ortho
! endif
! call routine
call write_l_r_wf
end
subroutine write_l_r_wf
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
output=trim(ezfio_filename)//'.tc_wf'
i_unit_output = getUnitAndOpen(output,'w')
integer :: i
print*,'Writing the left-right wf'
do i = 1, N_det
write(i_unit_output,*)i,psi_l_coef_sorted_bi_ortho_left(i),psi_r_coef_sorted_bi_ortho_right(i)
enddo
end
subroutine routine
implicit none
integer :: i,degree
integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2
double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e,e_pt2
double precision :: contrib_pt,e_corr,coef,contrib,phase
double precision :: accu_positive,accu_positive_pt, accu_positive_core,accu_positive_core_pt
e_pt2 = 0.d0
accu_positive = 0.D0
accu_positive_pt = 0.D0
accu_positive_core = 0.d0
accu_positive_core_pt = 0.d0
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0)
delta_e = e_tilde_00 - e_i0
coef_pt1 = htilde_ij / delta_e
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij)
contrib_pt = coef_pt1 * htilde_ij
e_pt2 += contrib_pt
coef = psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1)
contrib = coef * htilde_ij
e_corr += contrib
call get_excitation(HF_bitmask,psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'*********'
if(degree==1)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
else if(degree ==2)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
endif
print*,'coef_pt1 = ',coef_pt1
print*,'coef = ',coef
print*,'contrib_pt ',contrib_pt
print*,'contrib = ',contrib
if(contrib.gt.0.d0)then
accu_positive += contrib
if(h1==1.or.h2==1)then
accu_positive_core += contrib
endif
if(dabs(contrib).gt.1.d-5)then
print*,'Found a positive contribution to correlation energy !!'
endif
endif
if(contrib_pt.gt.0.d0)then
accu_positive_pt += contrib_pt
if(h2==1.or.h1==1)then
accu_positive_core_pt += contrib_pt
endif
endif
endif
enddo
print*,''
print*,''
print*,'Total correlation energy = ',e_corr
print*,'Total correlation energy PT = ',e_pt2
print*,'Positive contribution to ecorr = ',accu_positive
print*,'Positive contribution to ecorr PT = ',accu_positive_pt
print*,'Pure core contribution = ',accu_positive_core
print*,'Pure core contribution PT = ',accu_positive_core_pt
end

157
src/tc_bi_ortho/psi_det_tc_sorted.irp.f Normal file
View File

@ -0,0 +1,157 @@
use bitmasks
BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ]
implicit none
BEGIN_DOC
! Contribution of determinants to the state-averaged density.
END_DOC
integer :: i,j,k
double precision :: f
psi_average_norm_contrib_tc(:) = 0.d0
do k=1,N_states
do i=1,N_det
psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) + &
dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k)
enddo
enddo
f = 1.d0/sum(psi_average_norm_contrib_tc(1:N_det))
do i=1,N_det
psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i)*f
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ]
&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ]
&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ]
implicit none
BEGIN_DOC
! Wave function sorted by determinants contribution to the norm (state-averaged)
!
! psi_det_sorted_tc_order(i) -> k : index in psi_det
END_DOC
integer :: i,j,k
integer, allocatable :: iorder(:)
allocate ( iorder(N_det) )
do i=1,N_det
psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i)
iorder(i) = i
enddo
call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det)
do i=1,N_det
do j=1,N_int
psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i))
psi_det_sorted_tc(j,2,i) = psi_det(j,2,iorder(i))
enddo
psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_sorted_tc(i)
psi_det_sorted_tc_order(iorder(i)) = i
enddo
double precision :: accu
do k=1,N_states
accu = 0.d0
do i=1,N_det
psi_coef_sorted_tc(i,k) = dsqrt(dabs(psi_l_coef_bi_ortho(iorder(i),k)*psi_r_coef_bi_ortho(iorder(i),k)))
accu += psi_coef_sorted_tc(i,k)**2
enddo
accu = 1.d0/dsqrt(accu)
do i=1,N_det
psi_coef_sorted_tc(i,k) *= accu
enddo
enddo
psi_det_sorted_tc(:,:,N_det+1:psi_det_size) = 0_bit_kind
psi_coef_sorted_tc(N_det+1:psi_det_size,:) = 0.d0
psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0
psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0
deallocate(iorder)
END_PROVIDER
BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)]
&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)]
BEGIN_DOC
! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc
! psi_l_coef_sorted_bi_ortho : left coefficients corresponding to psi_det_sorted_tc
END_DOC
implicit none
integer :: i, j, k
psi_r_coef_sorted_bi_ortho = 0.d0
psi_l_coef_sorted_bi_ortho = 0.d0
do i = 1, N_det
psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_bit, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_bit, (psi_det_size,N_states) ]
implicit none
BEGIN_DOC
! Determinants on which we apply $\langle i|H|psi \rangle$ for perturbation.
! They are sorted by determinants interpreted as integers. Useful
! to accelerate the search of a random determinant in the wave
! function.
END_DOC
call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, size(psi_coef,1), &
psi_det_sorted_tc_bit, psi_coef_sorted_tc_bit, N_states)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_right, (N_int,2,N_det) ]
&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho_right, (N_det)]
implicit none
BEGIN_DOC
! psi_det_sorted_tc_right : Slater determinants sorted by decreasing value of |right- coefficients|
!
! psi_r_coef_sorted_bi_ortho_right : right wave function according to psi_det_sorted_tc_right
END_DOC
integer, allocatable :: iorder(:)
double precision, allocatable :: coef(:)
integer :: i,j
allocate ( iorder(N_det) , coef(N_det))
do i=1,N_det
coef(i) = -dabs(psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1))
iorder(i) = i
enddo
call dsort(coef,iorder,N_det)
do i=1,N_det
do j=1,N_int
psi_det_sorted_tc_right(j,1,i) = psi_det(j,1,iorder(i))
psi_det_sorted_tc_right(j,2,i) = psi_det(j,2,iorder(i))
enddo
psi_r_coef_sorted_bi_ortho_right(i) = psi_r_coef_bi_ortho(iorder(i),1)/psi_r_coef_bi_ortho(iorder(1),1)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_left, (N_int,2,N_det) ]
&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho_left, (N_det)]
implicit none
BEGIN_DOC
! psi_det_sorted_tc_left : Slater determinants sorted by decreasing value of |LEFTt- coefficients|
!
! psi_r_coef_sorted_bi_ortho_left : LEFT wave function according to psi_det_sorted_tc_left
END_DOC
integer, allocatable :: iorder(:)
double precision, allocatable :: coef(:)
integer :: i,j
allocate ( iorder(N_det) , coef(N_det))
do i=1,N_det
coef(i) = -dabs(psi_l_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1))
iorder(i) = i
enddo
call dsort(coef,iorder,N_det)
do i=1,N_det
do j=1,N_int
psi_det_sorted_tc_left(j,1,i) = psi_det(j,1,iorder(i))
psi_det_sorted_tc_left(j,2,i) = psi_det(j,2,iorder(i))
enddo
psi_l_coef_sorted_bi_ortho_left(i) = psi_l_coef_bi_ortho(iorder(i),1)/psi_l_coef_bi_ortho(iorder(1),1)
enddo
END_PROVIDER

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

@ -0,0 +1,44 @@
! ---
BEGIN_PROVIDER [ double precision, psi_bitcleft_bilinear_matrix_values, (N_det,N_states) ]
BEGIN_DOC
! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
! $D_\alpha^\dagger.C.D_\beta$
!
! Rows are $\alpha$ determinants and columns are $\beta$.
!
! Order refers to psi_det
END_DOC
use bitmasks
implicit none
integer :: k, l
if(N_det .eq. 1) then
do l = 1, N_states
psi_bitcleft_bilinear_matrix_values(1,l) = 1.d0
enddo
else
do l = 1, N_states
do k = 1, N_det
psi_bitcleft_bilinear_matrix_values(k,l) = psi_l_coef_bi_ortho(k,l)
enddo
enddo
PROVIDE psi_bilinear_matrix_order
do l = 1, N_states
call dset_order(psi_bitcleft_bilinear_matrix_values(1,l), psi_bilinear_matrix_order, N_det)
enddo
endif
END_PROVIDER
! ---

217
src/tc_bi_ortho/psi_r_l_prov.irp.f Normal file
View File

@ -0,0 +1,217 @@
use bitmasks
BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) ]
implicit none
BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file
! is empty.
END_DOC
integer :: i,k, N_int2
logical :: exists
character*(64) :: label
PROVIDE read_wf N_det mo_label ezfio_filename nproc
psi_l_coef_bi_ortho = 0.d0
do i=1,min(N_states,N_det)
psi_l_coef_bi_ortho(i,i) = 1.d0
enddo
if (mpi_master) then
if (read_wf) then
call ezfio_has_tc_bi_ortho_psi_l_coef_bi_ortho(exists)
! if (exists) then
! call ezfio_has_tc_bi_ortho_mo_label(exists)
! if (exists) then
! call ezfio_get_tc_bi_ortho_mo_label(label)
! exists = (label == mo_label)
! endif
! endif
if (exists) then
double precision, allocatable :: psi_l_coef_bi_ortho_read(:,:)
allocate (psi_l_coef_bi_ortho_read(N_det,N_states))
print *, 'Read psi_l_coef_bi_ortho', N_det, N_states
call ezfio_get_tc_bi_ortho_psi_l_coef_bi_ortho(psi_l_coef_bi_ortho_read)
do k=1,N_states
do i=1,N_det
psi_l_coef_bi_ortho(i,k) = psi_l_coef_bi_ortho_read(i,k)
enddo
enddo
deallocate(psi_l_coef_bi_ortho_read)
endif
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( psi_l_coef_bi_ortho, size(psi_l_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read psi_l_coef_bi_ortho with MPI'
endif
IRP_ENDIF
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) ]
implicit none
BEGIN_DOC
! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file
! is empty.
END_DOC
integer :: i,k, N_int2
logical :: exists
character*(64) :: label
PROVIDE read_wf N_det mo_label ezfio_filename nproc
psi_r_coef_bi_ortho = 0.d0
do i=1,min(N_states,N_det)
psi_r_coef_bi_ortho(i,i) = 1.d0
enddo
if (mpi_master) then
if (read_wf) then
call ezfio_has_tc_bi_ortho_psi_r_coef_bi_ortho(exists)
! if (exists) then
! call ezfio_has_tc_bi_ortho_mo_label(exists)
! if (exists) then
! call ezfio_get_tc_bi_ortho_mo_label(label)
! exists = (label == mo_label)
! endif
! endif
if (exists) then
double precision, allocatable :: psi_r_coef_bi_ortho_read(:,:)
allocate (psi_r_coef_bi_ortho_read(N_det,N_states))
print *, 'Read psi_r_coef_bi_ortho', N_det, N_states
call ezfio_get_tc_bi_ortho_psi_r_coef_bi_ortho(psi_r_coef_bi_ortho_read)
do k=1,N_states
do i=1,N_det
psi_r_coef_bi_ortho(i,k) = psi_r_coef_bi_ortho_read(i,k)
enddo
enddo
deallocate(psi_r_coef_bi_ortho_read)
endif
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( psi_r_coef_bi_ortho, size(psi_r_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read psi_r_coef_bi_ortho with MPI'
endif
IRP_ENDIF
END_PROVIDER
subroutine save_tc_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psilcoef,psircoef)
implicit none
BEGIN_DOC
! Save the wave function into the |EZFIO| file
END_DOC
use bitmasks
include 'constants.include.F'
integer, intent(in) :: ndet,nstates,dim_psicoef
integer(bit_kind), intent(in) :: psidet(N_int,2,ndet)
double precision, intent(in) :: psilcoef(dim_psicoef,nstates)
double precision, intent(in) :: psircoef(dim_psicoef,nstates)
integer*8, allocatable :: psi_det_save(:,:,:)
double precision, allocatable :: psil_coef_save(:,:)
double precision, allocatable :: psir_coef_save(:,:)
double precision :: accu_norm
integer :: i,j,k, ndet_qp_edit
if (mpi_master) then
ndet_qp_edit = min(ndet,N_det_qp_edit)
call ezfio_set_determinants_N_int(N_int)
call ezfio_set_determinants_bit_kind(bit_kind)
call ezfio_set_determinants_N_det(ndet)
call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit)
call ezfio_set_determinants_n_states(nstates)
call ezfio_set_determinants_mo_label(mo_label)
allocate (psi_det_save(N_int,2,ndet))
do i=1,ndet
do j=1,2
do k=1,N_int
psi_det_save(k,j,i) = transfer(psidet(k,j,i),1_8)
enddo
enddo
enddo
call ezfio_set_determinants_psi_det(psi_det_save)
call ezfio_set_determinants_psi_det_qp_edit(psi_det_save)
deallocate (psi_det_save)
allocate (psil_coef_save(ndet,nstates),psir_coef_save(ndet,nstates))
do k=1,nstates
do i=1,ndet
psil_coef_save(i,k) = psilcoef(i,k)
psir_coef_save(i,k) = psircoef(i,k)
enddo
enddo
call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(psil_coef_save)
call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(psir_coef_save)
deallocate (psil_coef_save,psir_coef_save)
! allocate (psi_coef_save(ndet_qp_edit,nstates))
! do k=1,nstates
! do i=1,ndet_qp_edit
! psi_coef_save(i,k) = psicoef(i,k)
! enddo
! enddo
!
! call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save)
! deallocate (psi_coef_save)
call write_int(6,ndet,'Saved determinantsi and psi_r/psi_l coef')
endif
end
subroutine save_tc_bi_ortho_wavefunction
implicit none
call save_tc_wavefunction_general(N_det,N_states,psi_det,size(psi_l_coef_bi_ortho, 1),psi_l_coef_bi_ortho,psi_r_coef_bi_ortho)
call routine_save_right_bi_ortho
end
subroutine routine_save_right_bi_ortho
implicit none
double precision, allocatable :: coef_tmp(:,:)
integer :: i
N_states = 1
allocate(coef_tmp(N_det, N_states))
do i = 1, N_det
coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
enddo
call save_wavefunction_general_unormalized(N_det,N_states,psi_det,size(coef_tmp,1),coef_tmp(1,1))
end
subroutine routine_save_left_right_bi_ortho
implicit none
double precision, allocatable :: coef_tmp(:,:)
integer :: i,n_states_tmp
n_states_tmp = 2
allocate(coef_tmp(N_det, n_states_tmp))
do i = 1, N_det
coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1)
coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1)
enddo
call save_wavefunction_general_unormalized(N_det,n_states_tmp,psi_det,size(coef_tmp,1),coef_tmp(1,1))
end

61
src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f Normal file
View File

@ -0,0 +1,61 @@
program save_bitcpsileft_for_qmcchem
integer :: iunit
logical :: exists
double precision :: e_ref
print *, ' '
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
print *, ' call save_for_qmcchem before '
print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
print *, ' '
call write_lr_spindeterminants()
e_ref = 0.d0
iunit = 13
open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
call ezfio_has_fci_energy_pt2(exists)
if(.not.exists) then
call ezfio_has_fci_energy(exists)
if(.not.exists) then
call ezfio_has_tc_scf_bitc_energy(exists)
if(exists) then
call ezfio_get_tc_scf_bitc_energy(e_ref)
endif
endif
endif
write(iunit,*) e_ref
close(iunit)
end
! --
subroutine write_lr_spindeterminants()
use bitmasks
implicit none
integer :: k, l
double precision, allocatable :: buffer(:,:)
PROVIDE psi_bitcleft_bilinear_matrix_values
allocate(buffer(N_det,N_states))
do l = 1, N_states
do k = 1, N_det
buffer(k,l) = psi_bitcleft_bilinear_matrix_values(k,l)
enddo
enddo
call ezfio_set_spindeterminants_psi_left_coef_matrix_values(buffer)
deallocate(buffer)
end subroutine write_lr_spindeterminants
! ---

15
src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f Normal file
View File

@ -0,0 +1,15 @@
program tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call routine_save_left_right_bi_ortho
! call test
end

33
src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f Normal file
View File

@ -0,0 +1,33 @@
program tc_natorb_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call save_tc_natorb
end
subroutine save_tc_natorb
implicit none
print*,'Saving the natorbs '
provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
call save_ref_determinant_nstates_1
call ezfio_set_determinants_read_wf(.False.)
end
subroutine save_ref_determinant_nstates_1
implicit none
use bitmasks
double precision :: buffer(1,N_states)
buffer = 0.d0
buffer(1,1) = 1.d0
call save_wavefunction_general(1,1,ref_bitmask,1,buffer)
end

61
src/tc_bi_ortho/select_dets_bi_ortho.irp.f Normal file
View File

@ -0,0 +1,61 @@
program tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
!!!!!!!!!!!!!!! WARNING NO 3-BODY
!!!!!!!!!!!!!!! WARNING NO 3-BODY
three_body_h_tc = .False.
touch three_body_h_tc
!!!!!!!!!!!!!!! WARNING NO 3-BODY
!!!!!!!!!!!!!!! WARNING NO 3-BODY
call routine_test
! call test
end
subroutine routine_test
implicit none
use bitmasks ! you need to include the bitmasks_module.f90 features
integer :: i,n_good,degree
integer(bit_kind), allocatable :: dets(:,:,:)
integer, allocatable :: iorder(:)
double precision, allocatable :: coef(:),coef_new(:,:)
double precision :: thr
allocate(coef(N_det), iorder(N_det))
do i = 1, N_det
iorder(i) = i
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree==1)then
coef(i) = -0.5d0
else
coef(i) = -dabs(coef_pt1_bi_ortho(i))
endif
enddo
call dsort(coef,iorder,N_det)
!thr = save_threshold
thr = 1d-15
n_good = 0
do i = 1, N_det
if(dabs(coef(i)).gt.thr)then
n_good += 1
endif
enddo
print*,'n_good = ',n_good
allocate(dets(N_int,2,n_good),coef_new(n_good,n_states))
do i = 1, n_good
dets(:,:,i) = psi_det(:,:,iorder(i))
coef_new(i,:) = psi_coef(iorder(i),:)
enddo
call save_wavefunction_general(n_good,n_states,dets,n_good,coef_new)
end

359
src/tc_bi_ortho/slater_tc.irp.f Normal file
View File

@ -0,0 +1,359 @@
!!!!!!
subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, 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
!!
!! 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) :: htot
double precision :: hmono,htwoe,hthree
integer :: degree
call get_excitation_degree(key_j, key_i, degree, Nint)
if(degree.gt.2)then
htot = 0.d0
else
call htilde_mu_mat_bi_ortho(key_j,key_i, Nint, hmono,htwoe,hthree,htot)
endif
end subroutine htilde_mu_mat_tot
! --
subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
implicit none
use bitmasks
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
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_bi_ortho(Nint, key_i, hmono, htwoe, htot)
else if (degree == 1)then
call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
else if(degree == 2)then
call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
endif
if(three_body_h_tc) then
if(degree == 2) then
if(.not.double_normal_ord) then
call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
endif
else if(degree == 1)then
call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
else if(degree == 0)then
call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
endif
endif
htot = hmono + htwoe + hthree
if(degree==0)then
htot += nuclear_repulsion
endif
end
subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot)
BEGIN_DOC
! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hmono,htwoe,htot
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
double precision :: get_mo_two_e_integral_tc_int
integer(bit_kind) :: key_i_core(Nint,2)
! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
!
! PROVIDE mo_integrals_erf_map core_energy nuclear_repulsion core_bitmask
! PROVIDE core_fock_operator
!
! PROVIDE j1b_gauss
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! do i = 1, Nint
! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
! enddo
! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
! hmono = core_energy - nuclear_repulsion
! else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
hmono = 0.d0
! endif
htwoe= 0.d0
htot = 0.d0
do ispin = 1, 2
do i = 1, Ne(ispin) !
ii = occ(i,ispin)
hmono += mo_bi_ortho_tc_one_e(ii,ii)
! if(j1b_gauss .eq. 1) then
! print*,'j1b not implemented for bi ortho TC'
! print*,'stopping ....'
! stop
! !hmono += mo_j1b_gauss_hermI (ii,ii) &
! ! + mo_j1b_gauss_hermII (ii,ii) &
! ! + mo_j1b_gauss_nonherm(ii,ii)
! endif
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! hmono += core_fock_operator(ii,ii) ! add the usual Coulomb - Exchange from the core
! endif
enddo
enddo
! alpha/beta two-body
ispin = 1
jspin = 2
do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1))
ii = occ(i,ispin)
do j = 1, Ne(jspin) ! electron 2
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii)
enddo
enddo
! alpha/alpha two-body
do i = 1, Ne(ispin)
ii = occ(i,ispin)
do j = i+1, Ne(ispin)
jj = occ(j,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
! beta/beta two-body
do i = 1, Ne(jspin)
ii = occ(i,jspin)
do j = i+1, Ne(jspin)
jj = occ(j,jspin)
htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii)
enddo
enddo
htot = hmono + htwoe
end
subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, 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, 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
integer :: other_spin(2)
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: get_mo_two_e_integral_tc_int,phase
! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
other_spin(1) = 2
other_spin(2) = 1
call get_excitation_degree(key_i, key_j, degree, Nint)
hmono = 0.d0
htwoe= 0.d0
htot = 0.d0
if(degree.ne.2)then
return
endif
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! do i = 1, Nint
! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
! enddo
! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
! else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
! endif
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
if(s1.ne.s2)then
! opposite spin two-body
! key_j, key_i
htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1)
if(double_normal_ord.and.+Ne(1).gt.2)then
htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ???
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(double_normal_ord.and.+Ne(1).gt.2)then
htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ???
htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ???
endif
endif
htwoe *= phase
htot = htwoe
end
subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, 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, 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)
! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e
!
! PROVIDE core_bitmask core_fock_operator mo_integrals_erf_map
! PROVIDE j1b_gauss
other_spin(1) = 2
other_spin(2) = 1
hmono = 0.d0
htwoe= 0.d0
htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.ne.1)then
return
endif
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! do i = 1, Nint
! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
! key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1))
! key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2))
! enddo
! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
! else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
! endif
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
hmono = mo_bi_ortho_tc_one_e(p1,h1) * phase
! if(j1b_gauss .eq. 1) then
! print*,'j1b not implemented for bi ortho TC'
! print*,'stopping ....'
! stop
! !hmono += ( mo_j1b_gauss_hermI (h1,p1) &
! ! + mo_j1b_gauss_hermII (h1,p1) &
! ! + mo_j1b_gauss_nonherm(h1,p1) ) * phase
! endif
! if(core_tc_op)then
! print*,'core_tc_op not already taken into account for bi ortho'
! print*,'stopping ...'
! stop
! hmono += phase * core_fock_operator(h1,p1)
! endif
! alpha/beta two-body
ispin = other_spin(s1)
if(s1==1)then
! single alpha
do i = 1, Ne(ispin) ! electron 2
ii = occ(i,ispin)
htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1)
enddo
else
! single beta
do i = 1, Ne(ispin) ! electron 1
ii = occ(i,ispin)
htwoe += mo_bi_ortho_tc_two_e(p1,ii,h1,ii)
enddo
endif
! ! same spin two-body
do i = 1, Ne(s1)
ii = occ(i,s1)
! (h1p1|ii ii) - (h1 ii| p1 ii)
htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1) - mo_bi_ortho_tc_two_e(p1,ii,ii,h1)
enddo
htwoe *= phase
htot = hmono + htwoe
end

293
src/tc_bi_ortho/slater_tc_3e.irp.f Normal file
View File

@ -0,0 +1,293 @@
subroutine provide_all_three_ints_bi_ortho
implicit none
BEGIN_DOC
! routine that provides all necessary three-electron integrals
END_DOC
if(three_body_h_tc)then
PROVIDE three_e_3_idx_direct_bi_ort three_e_3_idx_cycle_1_bi_ort three_e_3_idx_cycle_2_bi_ort
PROVIDE three_e_3_idx_exch23_bi_ort three_e_3_idx_exch13_bi_ort three_e_3_idx_exch12_bi_ort
PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_cycle_2_bi_ort
PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort
endif
if(.not.double_normal_ord)then
PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort
PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort
else
PROVIDE normal_two_body_bi_orth
endif
end
subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
BEGIN_DOC
! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
END_DOC
use bitmasks
implicit none
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2)
double precision, intent(out) :: hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: direct_int, exchange_int
double precision :: sym_3_e_int_from_6_idx_tensor
double precision :: three_e_diag_parrallel_spin
if(core_tc_op)then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core,occ,Ne,Nint)
else
call bitstring_to_list_ab(key_i,occ,Ne,Nint)
endif
hthree = 0.d0
if(Ne(1)+Ne(2).ge.3)then
!! ! alpha/alpha/beta three-body
do i = 1, Ne(1)
ii = occ(i,1)
do j = i+1, Ne(1)
jj = occ(j,1)
do m = 1, Ne(2)
mm = occ(m,2)
! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR
! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR
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
enddo
enddo
enddo
! beta/beta/alpha three-body
do i = 1, Ne(2)
ii = occ(i,2)
do j = i+1, Ne(2)
jj = occ(j,2)
do m = 1, Ne(1)
mm = occ(m,1)
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
enddo
enddo
enddo
! alpha/alpha/alpha three-body
do i = 1, Ne(1)
ii = occ(i,1) ! 1
do j = i+1, Ne(1)
jj = occ(j,1) ! 2
do m = j+1, Ne(1)
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
enddo
enddo
enddo
! beta/beta/beta three-body
do i = 1, Ne(2)
ii = occ(i,2) ! 1
do j = i+1, Ne(2)
jj = occ(j,2) ! 2
do m = j+1, Ne(2)
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
enddo
enddo
enddo
endif
end
subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
BEGIN_DOC
! <key_j | H_tilde | key_i> for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
!!
!! 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) :: hthree
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 :: direct_int,phase,exchange_int,three_e_single_parrallel_spin
double precision :: sym_3_e_int_from_6_idx_tensor
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
hthree = 0.d0
call get_excitation_degree(key_i,key_j,degree,Nint)
if(degree.ne.1)then
return
endif
if(core_tc_op)then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1))
key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
endif
call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2)
! alpha/alpha/beta three-body
! print*,'IN SLAT RULES'
if(Ne(1)+Ne(2).ge.3)then
! hole of spin s1 :: contribution from purely other spin
ispin = other_spin(s1) ! ispin is the other spin than s1
do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1
ii = occ(i,ispin)
do j = i+1, Ne(ispin) ! j has the same spin than s1
jj = occ(j,ispin)
! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is
! < h1 j i | p1 j i > - < h1 j i | p1 i j >
!
! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
! exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor
direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1)
hthree += direct_int - exchange_int
enddo
enddo
! hole of spin s1 :: contribution from mixed other spin / same spin
do i = 1, Ne(ispin) ! other spin
ii = occ(i,ispin) ! other spin
do j = 1, Ne(s1) ! same spin
jj = occ(j,s1) ! same spin
! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor
! exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor
direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1)
! < h1 j i | p1 j i > - < h1 j i | j p1 i >
hthree += direct_int - exchange_int
! print*,'h1,p1,ii,jj = ',h1,p1,ii,jj
! print*,direct_int, exchange_int
enddo
enddo
!
! hole of spin s1 :: PURE SAME SPIN CONTRIBUTIONS !!!
do i = 1, Ne(s1)
ii = occ(i,s1)
do j = i+1, Ne(s1)
jj = occ(j,s1)
! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1)
hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR
enddo
enddo
endif
hthree *= phase
end
subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
BEGIN_DOC
! <key_j | H_tilde | key_i> for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS
!!
!! 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) :: hthree
integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm
integer :: degree,exc(0:2,2,2)
integer :: h1, p1, h2, p2, s1, s2
double precision :: phase
integer :: other_spin(2)
integer(bit_kind) :: key_i_core(Nint,2)
double precision :: direct_int,exchange_int,sym_3_e_int_from_6_idx_tensor
double precision :: three_e_double_parrallel_spin
other_spin(1) = 2
other_spin(2) = 1
call get_excitation_degree(key_i, key_j, degree, Nint)
hthree = 0.d0
if(degree.ne.2)then
return
endif
if(core_tc_op)then
do i = 1, Nint
key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
enddo
call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
else
call bitstring_to_list_ab(key_i, occ, Ne, Nint)
endif
call get_double_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2)
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
else ! different spin excitation
do m = 1, Ne(s1)
mm = occ(m,s1) !
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
hthree += direct_int - exchange_int
enddo
do m = 1, Ne(s2)
mm = occ(m,s2) !
direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1)
exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
hthree += direct_int - exchange_int
enddo
endif
endif
hthree *= phase
end

111
src/tc_bi_ortho/symmetrized_3_e_int.irp.f Normal file
View File

@ -0,0 +1,111 @@
subroutine give_all_perm_for_three_e(n,l,k,m,j,i,idx_list,phase)
implicit none
BEGIN_DOC
! returns all the list of permutting indices for the antimmetrization of
!
! (k^dagger l^dagger n^dagger m j i) <nlk|L|mji> when all indices have the same spins
!
! idx_list(:,i) == list of the 6 indices corresponding the permutation "i"
!
! phase(i) == phase of the permutation "i"
!
! there are in total 6 permutations with different indices
END_DOC
integer, intent(in) :: n,l,k,m,j,i
integer, intent(out) :: idx_list(6,6)
double precision :: phase(6)
integer :: list(6)
!!! CYCLIC PERMUTATIONS
phase(1:3) = 1.d0
!!! IDENTITY PERMUTATION
list = (/n,l,k,m,j,i/)
idx_list(:,1) = list(:)
!!! FIRST CYCLIC PERMUTATION
list = (/n,l,k,j,i,m/)
idx_list(:,2) = list(:)
!!! FIRST CYCLIC PERMUTATION
list = (/n,l,k,i,m,j/)
idx_list(:,3) = list(:)
!!! NON CYCLIC PERMUTATIONS
phase(1:3) = -1.d0
!!! PARTICLE 1 is FIXED
list = (/n,l,k,j,m,i/)
idx_list(:,4) = list(:)
!!! PARTICLE 2 is FIXED
list = (/n,l,k,i,j,m/)
idx_list(:,5) = list(:)
!!! PARTICLE 3 is FIXED
list = (/n,l,k,m,i,j/)
idx_list(:,6) = list(:)
end
double precision function sym_3_e_int_from_6_idx_tensor(n,l,k,m,j,i)
implicit none
BEGIN_DOC
! returns all good combinations of permutations of integrals with the good signs
!
! for a given (k^dagger l^dagger n^dagger m j i) <nlk|L|mji> when all indices have the same spins
END_DOC
integer, intent(in) :: n,l,k,m,j,i
sym_3_e_int_from_6_idx_tensor = three_body_ints_bi_ort(n,l,k,m,j,i) & ! direct
+ three_body_ints_bi_ort(n,l,k,j,i,m) & ! 1st cyclic permutation
+ three_body_ints_bi_ort(n,l,k,i,m,j) & ! 2nd cyclic permutation
- three_body_ints_bi_ort(n,l,k,j,m,i) & ! elec 1 is kept fixed
- three_body_ints_bi_ort(n,l,k,i,j,m) & ! elec 2 is kept fixed
- three_body_ints_bi_ort(n,l,k,m,i,j) ! elec 3 is kept fixed
end
double precision function direct_sym_3_e_int(n,l,k,m,j,i)
implicit none
BEGIN_DOC
! returns all good combinations of permutations of integrals with the good signs
!
! for a given (k^dagger l^dagger n^dagger m j i) <nlk|L|mji> when all indices have the same spins
END_DOC
integer, intent(in) :: n,l,k,m,j,i
double precision :: integral
direct_sym_3_e_int = 0.d0
call give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral) ! direct
direct_sym_3_e_int += integral
call give_integrals_3_body_bi_ort(n,l,k,j,i,m,integral) ! 1st cyclic permutation
direct_sym_3_e_int += integral
call give_integrals_3_body_bi_ort(n,l,k,i,m,j,integral) ! 2nd cyclic permutation
direct_sym_3_e_int += integral
call give_integrals_3_body_bi_ort(n,l,k,j,m,i,integral) ! elec 1 is kept fixed
direct_sym_3_e_int += -integral
call give_integrals_3_body_bi_ort(n,l,k,i,j,m,integral) ! elec 2 is kept fixed
direct_sym_3_e_int += -integral
call give_integrals_3_body_bi_ort(n,l,k,m,i,j,integral) ! elec 3 is kept fixed
direct_sym_3_e_int += -integral
end
double precision function three_e_diag_parrallel_spin(m,j,i)
implicit none
integer, intent(in) :: i,j,m
three_e_diag_parrallel_spin = three_e_3_idx_direct_bi_ort(m,j,i) ! direct
three_e_diag_parrallel_spin += three_e_3_idx_cycle_1_bi_ort(m,j,i) + three_e_3_idx_cycle_2_bi_ort(m,j,i) & ! two cyclic permutations
- three_e_3_idx_exch23_bi_ort(m,j,i) - three_e_3_idx_exch13_bi_ort(m,j,i) & ! two first exchange
- three_e_3_idx_exch12_bi_ort(m,j,i) ! last exchange
end
double precision function three_e_single_parrallel_spin(m,j,k,i)
implicit none
integer, intent(in) :: i,k,j,m
three_e_single_parrallel_spin = three_e_4_idx_direct_bi_ort(m,j,k,i) ! direct
three_e_single_parrallel_spin += three_e_4_idx_cycle_1_bi_ort(m,j,k,i) + three_e_4_idx_cycle_2_bi_ort(m,j,k,i) & ! two cyclic permutations
- three_e_4_idx_exch23_bi_ort(m,j,k,i) - three_e_4_idx_exch13_bi_ort(m,j,k,i) & ! two first exchange
- three_e_4_idx_exch12_bi_ort(m,j,k,i) ! last exchange
end
double precision function three_e_double_parrallel_spin(m,l,j,k,i)
implicit none
integer, intent(in) :: i,k,j,m,l
three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! direct
three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations
- three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) & ! two first exchange
- three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange
end

61
src/tc_bi_ortho/tc_bi_ortho.irp.f Normal file
View File

@ -0,0 +1,61 @@
program tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end.
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call routine_diag
! call test
end
subroutine test
implicit none
integer :: i,j
double precision :: hmono,htwoe,hthree,htot
use bitmasks
print*,'test'
! call htilde_mu_mat_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot)
call double_htilde_mu_mat_bi_ortho(N_int,psi_det(1,1,1), psi_det(1,1,2), hmono, htwoe, htot)
print*,hmono, htwoe, htot
end
subroutine routine_diag
implicit none
! provide eigval_right_tc_bi_orth
provide overlap_bi_ortho
! provide htilde_matrix_elmt_bi_ortho
integer ::i,j
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)
print*,'e_tc_left_right = ',e_tc_left_right
print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00
print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth
print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single
print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double
print*,'***'
print*,'e_corr_bi_orth = ',e_corr_bi_orth
print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj
print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth
print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth
print*,'Left/right eigenvectors'
do i = 1,N_det
write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1)
enddo
do j=1,N_states
do i=1,N_det
psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
enddo
enddo
SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
call save_tc_bi_ortho_wavefunction
! call routine_save_left_right_bi_ortho
end

24
src/tc_bi_ortho/tc_bi_ortho_prop.irp.f Normal file
View File

@ -0,0 +1,24 @@
program tc_bi_ortho_prop
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
! call routine_diag
call test
end
subroutine test
implicit none
integer :: i
print*,'TC Dipole components'
do i= 1, 3
print*,tc_bi_ortho_dipole(i,1)
enddo
end

24
src/tc_bi_ortho/tc_cisd_sc2.irp.f Normal file
View File

@ -0,0 +1,24 @@
program tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call test
end
subroutine test
implicit none
! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
! allocate(dressing_dets(N_det),e_corr_dets(N_det))
! e_corr_dets = 0.d0
! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
provide eigval_tc_cisd_sc2_bi_ortho
end

110
src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f Normal file
View File

@ -0,0 +1,110 @@
BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)]
&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)]
implicit none
integer :: it,n_real,degree,i
double precision :: e_before, e_current,thr, hmono,htwoe,hthree
double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:)
double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det))
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 1 .or. degree == 2)then
call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i))
endif
enddo
do i = 1, N_det
e_corr_dets(i) = reigvec_tc_bi_orth(i,1) * h0j(i)/reigvec_tc_bi_orth(1,1)
enddo
print*,'Starting from ',eigval_right_tc_bi_orth(1)
e_before = 0.d0
e_current = 10.d0
thr = 1.d-5
it = 0
dressing_dets = 0.d0
do while (dabs(E_before-E_current).gt.thr)
it += 1
E_before = E_current
h_sc2 = htilde_matrix_elmt_bi_ortho
call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets)
do i = 1, N_det
print*,'dressing_dets(i) = ',dressing_dets(i)
h_sc2(i,i) += dressing_dets(i)
enddo
call non_hrmt_real_diag(N_det,h_sc2,&
leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
n_real,eigval_right_tmp)
do i = 1, N_det
e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1)
enddo
E_current = eigval_right_tmp(1)
print*,'it, E(SC)^2 = ',it,E_current
enddo
eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states)
reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states)
END_PROVIDER
subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets)
implicit none
use bitmasks
integer, intent(in) :: ndet
integer(bit_kind), intent(in) :: dets(N_int,2,ndet)
double precision, intent(in) :: e_corr_dets(ndet)
double precision, intent(out) :: dressing_dets(ndet)
integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:)
integer(bit_kind), allocatable :: hole_part(:,:,:)
integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2
integer(bit_kind) :: xorvec(2,N_int)
double precision :: phase
dressing_dets = 0.d0
allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet))
do i = 2, ndet
call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int)
do j = 1, N_int
hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i))
hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i))
enddo
if(degree(i) == 1)then
call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
else if(degree(i) == 2)then
call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int)
endif
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
hole(1,i) = h1
hole(2,i) = h2
part(1,i) = p1
part(2,i) = p2
spin(1,i) = s1
spin(2,i) = s2
enddo
integer :: same
if(elec_alpha_num+elec_beta_num<3)return
do i = 2, ndet
do j = i+1, ndet
same = 0
if(degree(i) == degree(j) .and. degree(i)==1)cycle
do k = 1, N_int
xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j))
xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j))
same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2))
enddo
! print*,'i,j',i,j
! call debug_det(dets(1,1,i),N_int)
! call debug_det(hole_part(1,1,i),N_int)
! call debug_det(dets(1,1,j),N_int)
! call debug_det(hole_part(1,1,j),N_int)
! print*,'same = ',same
if(same.eq.0)then
dressing_dets(i) += e_corr_dets(j)
dressing_dets(j) += e_corr_dets(i)
endif
enddo
enddo
end

179
src/tc_bi_ortho/tc_h_eigvectors.irp.f Normal file
View File

@ -0,0 +1,179 @@
use bitmasks
BEGIN_PROVIDER [ integer, index_HF_psi_det]
implicit none
integer :: i,degree
do i = 1, N_det
call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
if(degree == 0)then
index_HF_psi_det = i
exit
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)]
&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)]
&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)]
&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)]
&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ]
BEGIN_DOC
! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis
END_DOC
implicit none
integer :: i, idx_dress, j, istate
logical :: converged, dagger
integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
PROVIDE N_det N_int
if(n_det.le.N_det_max_full)then
allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det))
call non_hrmt_real_diag(N_det,htilde_matrix_elmt_bi_ortho,&
leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,&
n_real_tc_bi_orth_eigval_right,eigval_right_tmp)
double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
integer, allocatable :: iorder(:)
allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det))
do i = 1,N_det
iorder(i) = i
coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
enddo
call dsort(coef_hf_r,iorder,N_det)
igood_r = iorder(1)
print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1)
do i = 1,N_det
iorder(i) = i
coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i))
enddo
call dsort(coef_hf_l,iorder,N_det)
igood_l = iorder(1)
print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1)
if(igood_r.ne.igood_l.and.igood_r.ne.1)then
print *,''
print *,'Warning, the left and right eigenvectors are "not the same" '
print *,'Warning, the ground state is not dominated by HF...'
print *,'State with largest RIGHT coefficient of HF ',igood_r
print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r)
print *,'State with largest LEFT coefficient of HF ',igood_l
print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l)
endif
if(state_following_tc)then
print *,'Following the states with the largest coef on HF'
print *,'igood_r,igood_l',igood_r,igood_l
i= igood_r
eigval_right_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det
reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i)
! print*,reigvec_tc_bi_orth(j,1)
enddo
i= igood_l
eigval_left_tc_bi_orth(1) = eigval_right_tmp(i)
do j = 1, N_det
leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i)
enddo
else
do i = 1, N_states
eigval_right_tc_bi_orth(i) = eigval_right_tmp(i)
eigval_left_tc_bi_orth(i) = eigval_right_tmp(i)
do j = 1, N_det
reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i)
leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i)
enddo
enddo
endif
else
double precision, allocatable :: H_jj(:),vec_tmp(:,:)
external htc_bi_ortho_calc_tdav
external htcdag_bi_ortho_calc_tdav
allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag))
do i = 1, N_det
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i))
enddo
!!!! Preparing the left-eigenvector
print*,'Computing the left-eigenvector '
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(:,istate) = psi_l_coef_bi_ortho(:,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav)
do istate = 1, N_states
leigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
enddo
print*,'Computing the right-eigenvector '
!!!! Preparing the right-eigenvector
vec_tmp = 0.d0
do istate = 1, N_states
vec_tmp(:,istate) = psi_r_coef_bi_ortho(:,istate)
enddo
do istate = N_states+1, n_states_diag
vec_tmp(istate,istate) = 1.d0
enddo
call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav)
do istate = 1, N_states
reigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate)
enddo
deallocate(H_jj)
endif
call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,N_det,N_det,N_states)
print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1)
norm_ground_left_right_bi_orth = 0.d0
do j = 1, N_det
norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1)
enddo
print*,'norm l/r = ',norm_ground_left_right_bi_orth
END_PROVIDER
subroutine bi_normalize(u_l,u_r,n,ld,nstates)
!!!! Normalization of the scalar product of the left/right eigenvectors
double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates)
integer, intent(in) :: n,ld,nstates
integer :: i
double precision :: accu, tmp
do i = 1, nstates
!!!! Normalization of right eigenvectors |Phi>
accu = 0.d0
do j = 1, n
accu += u_r(j,i) * u_r(j,i)
enddo
accu = 1.d0/dsqrt(accu)
print*,'accu_r = ',accu
do j = 1, n
u_r(j,i) *= accu
enddo
tmp = u_r(1,i) / dabs(u_r(1,i))
do j = 1, n
u_r(j,i) *= tmp
enddo
!!!! Adaptation of the norm of the left eigenvector such that <chi|Phi> = 1
accu = 0.d0
do j = 1, n
accu += u_l(j,i) * u_r(j,i)
! print*,j, u_l(j,i) , u_r(j,i)
enddo
if(accu.gt.0.d0)then
accu = 1.d0/dsqrt(accu)
else
accu = 1.d0/dsqrt(-accu)
endif
tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i))
do j = 1, n
u_l(j,i) *= accu * tmp
u_r(j,i) *= accu
enddo
enddo
end

41
src/tc_bi_ortho/tc_hmat.irp.f Normal file
View File

@ -0,0 +1,41 @@
BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
BEGIN_DOC
! htilde_matrix_elmt_bi_ortho(j,i) = <J| H^tilde |I>
!
! WARNING !!!!!!!!! IT IS NOT HERMITIAN !!!!!!!!!
END_DOC
implicit none
integer :: i, j
double precision :: hmono,htwoe,hthree,htot
PROVIDE N_int
!$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hmono, htwoe, hthree, htot) &
!$OMP SHARED (N_det, psi_det, N_int,htilde_matrix_elmt_bi_ortho)
do i = 1, N_det
do j = 1, N_det
! < J | Htilde | I >
call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
htilde_matrix_elmt_bi_ortho(j,i) = htot
enddo
enddo
!$OMP END PARALLEL DO
! print*,'htilde_matrix_elmt_bi_ortho = '
! do i = 1, min(100,N_det)
! write(*,'(100(F16.10,X))')htilde_matrix_elmt_bi_ortho(1:min(100,N_det),i)
! enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)]
implicit none
integer ::i,j
do i = 1, N_det
do j = 1, N_det
htilde_matrix_elmt_bi_ortho_tranp(j,i) = htilde_matrix_elmt_bi_ortho(i,j)
enddo
enddo
END_PROVIDER

268
src/tc_bi_ortho/tc_prop.irp.f Normal file
View File

@ -0,0 +1,268 @@
BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ]
implicit none
BEGIN_DOC
! tc_transition_matrix(p,h,istate,jstate) = <Chi_istate| a^\dagger_p a_h |Phi_jstate>
!
! where <Chi_istate| and |Phi_jstate> are the left/right eigenvectors on a bi-ortho basis
END_DOC
integer :: i,j,istate,jstate,m,n,p,h
double precision :: phase
integer, allocatable :: occ(:,:)
integer :: n_occ_ab(2),degree,exc(0:2,2,2)
allocate(occ(N_int*bit_kind_size,2))
tc_transition_matrix = 0.d0
do istate = 1, N_states
do jstate = 1, N_states
do i = 1, N_det
do j = 1, N_det
call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int)
if(degree.gt.1)then
cycle
else if (degree == 0)then
call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int)
do p = 1, n_occ_ab(1) ! browsing the alpha electrons
m = occ(p,1)
tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
do p = 1, n_occ_ab(2) ! browsing the beta electrons
m = occ(p,1)
tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
enddo
else
call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int)
if (exc(0,1,1) == 1) then
! Single alpha
h = exc(1,1,1) ! hole in psi_det(1,1,j)
p = exc(1,2,1) ! particle in psi_det(1,1,j)
else
! Single beta
h = exc(1,1,2) ! hole in psi_det(1,1,j)
p = exc(1,2,2) ! particle in psi_det(1,1,j)
endif
tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate)
endif
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)]
implicit none
BEGIN_DOC
! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals)
! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !!
END_DOC
double precision, allocatable :: dm_tmp(:,:)
integer :: i,j,k,n_real
allocate( dm_tmp(mo_num,mo_num))
dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1)
print*,'dm_tmp'
do i = 1, mo_num
write(*,'(100(F16.10,X))')-dm_tmp(:,i)
enddo
! call non_hrmt_diag_split_degen( mo_num, dm_tmp&
call non_hrmt_fock_mat( mo_num, dm_tmp&
! call non_hrmt_bieig( mo_num, dm_tmp&
, natorb_tc_leigvec_mo, natorb_tc_reigvec_mo&
, n_real, natorb_tc_eigval )
double precision :: accu
accu = 0.d0
do i = 1, n_real
print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i)
accu += -natorb_tc_eigval(i)
enddo
print*,'accu = ',accu
dm_tmp = 0.d0
do i = 1, n_real
accu = 0.d0
do k = 1, mo_num
accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i)
enddo
accu = 1.d0/dsqrt(dabs(accu))
natorb_tc_reigvec_mo(:,i) *= accu
natorb_tc_leigvec_mo(:,i) *= accu
do j = 1, n_real
do k = 1, mo_num
dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j)
enddo
enddo
enddo
double precision :: accu_d, accu_nd
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, mo_num
accu_d += dm_tmp(i,i)
! write(*,'(100(F16.10,X))')dm_tmp(:,i)
do j = 1, mo_num
if(i==j)cycle
accu_nd += dabs(dm_tmp(j,i))
enddo
enddo
print*,'Trace of the overlap between TC natural orbitals ',accu_d
print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)]
implicit none
integer ::i,j,k
print*,'Diagonal elements of the Fock matrix before '
do i = 1, mo_num
write(*,*)i,Fock_matrix_tc_mo_tot(i,i)
enddo
double precision, allocatable :: fock_diag(:)
allocate(fock_diag(mo_num))
fock_diag = 0.d0
do i = 1, mo_num
fock_diag(i) = 0.d0
do j = 1, mo_num
do k = 1, mo_num
fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i)
enddo
enddo
enddo
integer, allocatable :: iorder(:)
allocate(iorder(mo_num))
do i = 1, mo_num
iorder(i) = i
enddo
call dsort(fock_diag,iorder,mo_num)
print*,'Diagonal elements of the Fock matrix after '
do i = 1, mo_num
write(*,*)i,fock_diag(i)
enddo
do i = 1, mo_num
fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i))
do j = 1, mo_num
fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i))
fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i))
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ]
BEGIN_DOC
! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
!
! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo
END_DOC
implicit none
integer :: i, j, k, q, p
double precision :: accu, accu_d
double precision, allocatable :: tmp(:,:)
! ! MO_R x R
call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
, mo_r_coef, size(mo_r_coef, 1) &
, fock_diag_sorted_r_natorb, size(fock_diag_sorted_r_natorb, 1) &
, 0.d0, natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) )
!
! MO_L x L
call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
, mo_l_coef, size(mo_l_coef, 1) &
, fock_diag_sorted_l_natorb, size(fock_diag_sorted_l_natorb, 1) &
, 0.d0, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1) )
allocate( tmp(mo_num,ao_num) )
! tmp <-- L.T x S_ao
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
, 0.d0, tmp, size(tmp, 1) )
! S <-- tmp x R
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
, tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) &
, 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) )
deallocate( tmp )
! ---
double precision :: norm
do i = 1, mo_num
norm = 1.d0/dsqrt(dabs(overlap_natorb_tc_eigvec_ao(i,i)))
do j = 1, mo_num
natorb_tc_reigvec_ao(j,i) *= norm
natorb_tc_leigvec_ao(j,i) *= norm
enddo
enddo
allocate( tmp(mo_num,ao_num) )
! tmp <-- L.T x S_ao
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
, 0.d0, tmp, size(tmp, 1) )
! S <-- tmp x R
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
, tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) &
, 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) )
deallocate( tmp )
accu_d = 0.d0
accu = 0.d0
do i = 1, mo_num
accu_d += overlap_natorb_tc_eigvec_ao(i,i)
do j = 1, mo_num
if(i==j)cycle
accu += dabs(overlap_natorb_tc_eigvec_ao(j,i))
enddo
enddo
print*,'Trace of the overlap_natorb_tc_eigvec_ao = ',accu_d
print*,'mo_num = ',mo_num
print*,'L1 norm of extra diagonal elements of overlap matrix ',accu
accu = accu / dble(mo_num**2)
END_PROVIDER
BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)]
implicit none
integer :: i,j,istate,m
double precision :: nuclei_part(3)
tc_bi_ortho_dipole = 0.d0
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i)
tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i)
tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i)
enddo
enddo
enddo
nuclei_part = 0.d0
do m = 1, 3
do i = 1,nucl_num
nuclei_part(m) += nucl_charge(i) * nucl_coord(i,m)
enddo
enddo
!
do istate = 1, N_states
do m = 1, 3
tc_bi_ortho_dipole(m,istate) += nuclei_part(m)
enddo
enddo
END_PROVIDER

73
src/tc_bi_ortho/test_normal_order.irp.f Normal file
View File

@ -0,0 +1,73 @@
program test_normal_order
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call test
end
subroutine test
implicit none
use bitmasks ! you need to include the bitmasks_module.f90 features
integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2)
integer :: exc(0:2,2,2)
integer(bit_kind), allocatable :: det_i(:,:)
double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal
integer, allocatable :: occ(:,:)
allocate( occ(N_int*bit_kind_size,2) )
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
allocate(det_i(N_int,2))
s1 = 1
s2 = 2
accu = 0.d0
do h1 = 1, elec_beta_num
do p1 = elec_beta_num+1, mo_num
do h2 = 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)
call htilde_mu_mat_bi_ortho(det_i,HF_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_ab(p2,h2,p1,h1)
accu += dabs(hthree-normal)
enddo
enddo
enddo
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
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
end

131
src/tc_bi_ortho/test_tc_bi_ortho.irp.f Normal file
View File

@ -0,0 +1,131 @@
program tc_bi_ortho
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
! call routine_2
call test_rout
end
subroutine test_rout
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'
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)
enddo
end
subroutine routine
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)
end

169
src/tc_bi_ortho/test_tc_fock.irp.f Normal file
View File

@ -0,0 +1,169 @@
program test_tc_fock
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
read_wf = .True.
touch read_wf
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
!call routine_1
!call routine_2
call routine_3()
end
! ---
subroutine routine_0
implicit none
use bitmasks ! you need to include the bitmasks_module.f90 features
integer :: i,a,j,m,i_ok
integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2,degree
integer(bit_kind), allocatable :: det_i(:,:)
double precision :: hmono,htwoe,hthree,htilde_ij,phase
double precision :: same, op, tot, accu
allocate(det_i(N_int,2))
s1 = 1
accu = 0.d0
do i = 1, elec_alpha_num ! occupied
do a = elec_alpha_num+1, mo_num ! virtual
det_i = ref_bitmask
call do_single_excitation(det_i,i,a,s1,i_ok)
if(i_ok == -1)then
print*,'PB !!'
print*,i,a
stop
endif
! call debug_det(det_i,N_int)
call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
op = fock_3_mat_a_op_sh_bi_orth(a,i)
same = fock_3_mat_a_sa_sh_bi_orth(a,i)
! same = 0.d0
tot = same + op
if(dabs(tot - phase*hthree).gt.1.d-10)then
print*,'------'
print*,i,a,phase
print*,'hthree = ',phase*hthree
print*,'fock = ',tot
print*,'same,op= ',same,op
print*,dabs(tot - phase*hthree)
stop
endif
accu += dabs(tot - phase*hthree)
enddo
enddo
print*,'accu = ',accu
end subroutine routine_0
! ---
subroutine routine_1
implicit none
integer :: i, a
double precision :: accu
accu = 0.d0
do i = 1, mo_num
do a = 1, mo_num
accu += dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) )
!if(dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) .gt. 1.d-10)then
print*, i, a
print*, dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) &
, fock_3_mat_a_op_sh_bi_orth_old(a,i), fock_3_mat_a_op_sh_bi_orth(a,i)
!endif
enddo
enddo
print *, 'accu = ', accu
end subroutine routine_1
! ---
subroutine routine_2
implicit none
integer :: i, a
double precision :: accu
accu = 0.d0
do i = 1, mo_num
do a = 1, mo_num
accu += dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) )
!if(dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) .gt. 1.d-10)then
print*, i, a
print*, dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) &
, fock_3_mat_a_sa_sh_bi_orth_old(a,i), fock_3_mat_a_sa_sh_bi_orth(a,i)
!endif
enddo
enddo
print *, 'accu = ', accu
end subroutine routine_2
! ---
subroutine routine_3()
use bitmasks ! you need to include the bitmasks_module.f90 features
implicit none
integer :: i, a, i_ok, s1
double precision :: hmono, htwoe, hthree, htilde_ij
double precision :: err_ai, err_tot
integer(bit_kind), allocatable :: det_i(:,:)
allocate(det_i(N_int,2))
err_tot = 0.d0
s1 = 1
det_i = ref_bitmask
call debug_det(det_i, N_int)
print*, ' HF det'
call debug_det(det_i, N_int)
do i = 1, elec_alpha_num ! occupied
do a = elec_alpha_num+1, mo_num ! virtual
det_i = ref_bitmask
call do_single_excitation(det_i, i, a, s1, i_ok)
if(i_ok == -1) then
print*, 'PB !!'
print*, i, a
stop
endif
!print*, ' excited det'
!call debug_det(det_i, N_int)
call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
err_ai = dabs(htilde_ij)
if(err_ai .gt. 1d-7) then
print*, ' warning on', i, a
print*, hmono, htwoe, htilde_ij
endif
err_tot += err_ai
write(22, *) htilde_ij
enddo
enddo
print *, ' err_tot = ', err_tot
deallocate(det_i)
end subroutine routine_3
! ---