9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-13 14:34:30 +02:00
qp2/src/mu_of_r/f_hf_utils.irp.f
eginer 46dfb551a8
Some checks failed
continuous-integration/drone/push Build is failing
modified the factor two in rdm
2023-03-04 13:07:12 +01:00

145 lines
5.9 KiB
Fortran

BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max_occ_val_orb_for_hf,n_max_occ_val_orb_for_hf)]
implicit none
BEGIN_DOC
! list of two-electron integrals (built with the MOs belonging to the \mathcal{B} space)
!
! needed to compute the function f_{HF}(r_1,r_2)
!
! two_e_int_hf_f(j,i,n,m) = < j i | n m > where all orbitals belong to "list_basis"
END_DOC
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1
m = list_valence_orb_for_hf(orb_m,1)
do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2
n = list_valence_orb_for_hf(orb_n,1)
do orb_i = 1, n_basis_orb ! electron 1
i = list_basis(orb_i)
do orb_j = 1, n_basis_orb ! electron 2
j = list_basis(orb_j)
! 2 1 2 1
two_e_int_hf_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
enddo
enddo
enddo
enddo
END_PROVIDER
subroutine f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens)
implicit none
BEGIN_DOC
! f_HF_val_ab(r1,r2) = function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937)
!
! two_bod_dens = on-top pair density of the HF wave function
!
! < HF | wee_{\alpha\beta} | HF > = \int (r1,r2) f_HF_ab(r1,r2) excluding all contributions from "core" "electrons"
END_DOC
double precision, intent(in) :: r1(3), r2(3)
double precision, intent(out):: f_HF_val_ab,two_bod_dens
integer :: i,j,m,n,i_m,i_n
integer :: i_i,i_j
double precision :: mo_two_e_integral
double precision, allocatable :: mos_array_r1(:)
double precision, allocatable :: mos_array_r2(:)
double precision, allocatable :: mos_array_valence_r1(:),mos_array_valence_r2(:)
double precision, allocatable :: mos_array_valence_hf_r1(:),mos_array_valence_hf_r2(:)
double precision :: get_two_e_integral
allocate(mos_array_valence_r1(n_basis_orb) , mos_array_valence_r2(n_basis_orb), mos_array_r1(mo_num), mos_array_r2(mo_num))
allocate(mos_array_valence_hf_r1(n_occ_val_orb_for_hf(1)) , mos_array_valence_hf_r2(n_occ_val_orb_for_hf(2)) )
! You get all orbitals in r_1 and r_2
call give_all_mos_at_r(r1,mos_array_r1)
call give_all_mos_at_r(r2,mos_array_r2)
! You extract the occupied ALPHA/BETA orbitals belonging to the space \mathcal{A}
do i_m = 1, n_occ_val_orb_for_hf(1)
mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1))
enddo
do i_m = 1, n_occ_val_orb_for_hf(2)
mos_array_valence_hf_r2(i_m) = mos_array_r2(list_valence_orb_for_hf(i_m,2))
enddo
! You extract the orbitals belonging to the space \mathcal{B}
do i_m = 1, n_basis_orb
mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m))
mos_array_valence_r2(i_m) = mos_array_r2(list_basis(i_m))
enddo
f_HF_val_ab = 0.d0
two_bod_dens = 0.d0
! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
do m = 1, n_occ_val_orb_for_hf(1)! electron 1
! You browse all OCCUPIED BETA electrons in the \mathcal{A} space
do n = 1, n_occ_val_orb_for_hf(2)! electron 2
! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
two_bod_dens += mos_array_valence_hf_r1(m) * mos_array_valence_hf_r1(m) * mos_array_valence_hf_r2(n) * mos_array_valence_hf_r2(n)
! You browse all COUPLE OF ORBITALS in the \mathacal{B} space
do i = 1, n_basis_orb
do j = 1, n_basis_orb
! 2 1 2 1
f_HF_val_ab += two_e_int_hf_f(j,i,n,m) &
* mos_array_valence_r1(i) * mos_array_valence_hf_r1(m) &
* mos_array_valence_r2(j) * mos_array_valence_hf_r2(n)
enddo
enddo
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
f_HF_val_ab *= 2.d0
two_bod_dens *= 2.d0
end
subroutine integral_f_HF_valence_ab(r1,int_f_HF_val_ab)
implicit none
BEGIN_DOC
! in_f_HF_val_ab(r_1) = \int dr_2 f_{\Psi^B}(r_1,r_2)
!
! where f_{\Psi^B}(r_1,r_2) is defined by Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
! for alpha beta spins and an HF wave function and excluding the "core" orbitals (see Eq. 16a of Phys.Chem.Lett.2019, 10, 2931 2937)
!
! Such function can be used to test if the f_HF_val_ab(r_1,r_2) is correctly built.
!
! < HF | wee_{\alpha\beta} | HF > = \int (r1) int_f_HF_val_ab(r_1)
END_DOC
double precision, intent(in) :: r1(3)
double precision, intent(out):: int_f_HF_val_ab
integer :: i,j,m,n,i_m,i_n
integer :: i_i,i_j
double precision :: mo_two_e_integral
double precision :: mos_array_r1(mo_num)
double precision, allocatable :: mos_array_valence_r1(:)
double precision, allocatable :: mos_array_valence_hf_r1(:)
double precision :: get_two_e_integral
call give_all_mos_at_r(r1,mos_array_r1)
allocate(mos_array_valence_r1( n_basis_orb ))
allocate(mos_array_valence_hf_r1( n_occ_val_orb_for_hf(1) ) )
do i_m = 1, n_occ_val_orb_for_hf(1)
mos_array_valence_hf_r1(i_m) = mos_array_r1(list_valence_orb_for_hf(i_m,1))
enddo
do i_m = 1, n_basis_orb
mos_array_valence_r1(i_m) = mos_array_r1(list_basis(i_m))
enddo
int_f_HF_val_ab = 0.d0
! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
do m = 1, n_occ_val_orb_for_hf(1)! electron 1
! You browse all OCCUPIED BETA electrons in the \mathcal{A} space
do n = 1, n_occ_val_orb_for_hf(2)! electron 2
! You browse all ORBITALS in the \mathacal{B} space
do i = 1, n_basis_orb
! due to integration in real-space and the use of orthonormal MOs, a Kronecker delta_jn shoes up
j = n
! 2 1 2 1
int_f_HF_val_ab += two_e_int_hf_f(j,i,n,m) &
* mos_array_valence_r1(i) * mos_array_valence_hf_r1(m)
enddo
enddo
enddo
! multiply by two to adapt to the N(N-1) normalization condition of the active two-rdm
int_f_HF_val_ab *= 2.d0
end