PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 04:43:45 +01:00
Code Releases Activity

added mu_of_r

Browse Source
This commit is contained in:
Emmanuel Giner 2020-04-07 11:01:24 +02:00
parent c4ded9dcfb
commit 51cf96a506
12 changed files with 1233 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/mu_of_r/EZFIO.cfg Normal file
View File

@ -0,0 +1,18 @@
[mu_of_r_disk]
type: double precision
doc: array of the values of mu(r)
interface: ezfio
size: (becke_numerical_grid.n_points_final_grid,determinants.n_states)
[mu_of_r_potential]
type: character*(32)
doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated]
interface: ezfio, provider, ocaml
default: hf
[io_mu_of_r]
type: Disk_access
doc: Read/Write the mu(r) from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None

1
src/mu_of_r/NEED Normal file
View File

@ -0,0 +1 @@
cas_based_on_top

4
src/mu_of_r/README.rst Normal file
View File

@ -0,0 +1,4 @@
==================
mu_of_r_definition
==================

106
src/mu_of_r/basis_def.irp.f Normal file
View File

@ -0,0 +1,106 @@
BEGIN_PROVIDER [integer, n_occ_val_orb_for_hf,(2)]
&BEGIN_PROVIDER [integer, n_max_occ_val_orb_for_hf]
implicit none
BEGIN_DOC
! Number of OCCUPIED VALENCE ORBITALS for each spin to build the f_{HF}(r_1,r_2) function
!
! This is typically elec_alpha_num - n_core_orb for alpha electrons and elec_beta_num - n_core_orb for beta electrons
!
! This determines the size of the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937
END_DOC
integer :: i
n_occ_val_orb_for_hf = 0
! You browse the ALPHA ELECTRONS and check if its not a CORE ORBITAL
do i = 1, elec_alpha_num
if( trim(mo_class(i))=="Inactive" &
.or. trim(mo_class(i))=="Active" &
.or. trim(mo_class(i))=="Virtual" )then
n_occ_val_orb_for_hf(1) +=1
endif
enddo
! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL
do i = 1, elec_beta_num
if( trim(mo_class(i))=="Inactive" &
.or. trim(mo_class(i))=="Active" &
.or. trim(mo_class(i))=="Virtual" )then
n_occ_val_orb_for_hf(2) +=1
endif
enddo
n_max_occ_val_orb_for_hf = maxval(n_occ_val_orb_for_hf)
END_PROVIDER
BEGIN_PROVIDER [integer, list_valence_orb_for_hf, (n_max_occ_val_orb_for_hf,2)]
implicit none
BEGIN_DOC
! List of OCCUPIED valence orbitals for each spin to build the f_{HF}(r_1,r_2) function
!
! This corresponds to ALL OCCUPIED orbitals in the HF wave function, except those defined as "core"
!
! This determines the space \mathcal{A} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937
END_DOC
integer :: i,j
j = 0
! You browse the ALPHA ELECTRONS and check if its not a CORE ORBITAL
do i = 1, elec_alpha_num
if( trim(mo_class(i))=="Inactive" &
.or. trim(mo_class(i))=="Active" &
.or. trim(mo_class(i))=="Virtual" )then
j +=1
list_valence_orb_for_hf(j,1) = i
endif
enddo
j = 0
! You browse the BETA ELECTRONS and check if its not a CORE ORBITAL
do i = 1, elec_beta_num
if( trim(mo_class(i))=="Inactive" &
.or. trim(mo_class(i))=="Active" &
.or. trim(mo_class(i))=="Virtual" )then
j +=1
list_valence_orb_for_hf(j,2) = i
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [integer, n_basis_orb]
implicit none
BEGIN_DOC
! Defines the number of orbitals you will use to explore the basis set
!
! This determines the size of the space \mathcal{B} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! It corresponds to all MOs except those defined as "deleted"
END_DOC
n_basis_orb = n_all_but_del_orb
END_PROVIDER
BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)]
implicit none
BEGIN_DOC
! Defines the set of orbitals you will use to explore the basis set
!
! This determines the space \mathcal{B} of Eqs. (15-16) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! It corresponds to all MOs except those defined as "deleted"
END_DOC
integer :: i
do i = 1, n_all_but_del_orb
list_basis(i) = list_all_but_del_orb(i)
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)]
implicit none
integer :: ipoint,i,ii
do ipoint = 1, n_points_final_grid
do i = 1, n_basis_orb
ii = list_basis(i)
basis_mos_in_r_array(i,ipoint) = mos_in_r_array(ii,ipoint)
enddo
enddo
END_PROVIDER

202
src/mu_of_r/example.irp.f Normal file
View File

@ -0,0 +1,202 @@
subroutine test_f_HF_valence_ab
implicit none
BEGIN_DOC
! routine to test the function f_HF(r1,r2)
!
! the integral over r1,r2 should be equal to the alpha/beta interaction of HF determinant
END_DOC
integer :: ipoint,i,j,i_i,j_j,jpoint
double precision :: accu_val,accu_ful, weight1,weight2, r1(3),integral_psi_val,integral_psi,r2(3),two_bod
accu_2 = 0.d0
! You compute the coulomb repulsion between alpha-beta electrons for HF
do i = 1, n_occ_val_orb_for_hf(1)
i_i = list_valence_orb_for_hf(i,1)
do j = 1, n_occ_val_orb_for_hf(2)
j_j = list_valence_orb_for_hf(j,2)
accu_2 += mo_two_e_integrals_jj(j_j,i_i)
enddo
enddo
print*,''
print*,''
print*,''
print*,'**************************'
print*,'**************************'
print*,'Routine to test the f_HF(r1,r2) function'
print*,'**************************'
print*,''
print*,''
print*,''
print*,'**************************'
print*,'<HF| We_ee^{ab}|HF> = ',accu_2
print*,'**************************'
print*,'semi analytical form '
accu_val = 0.d0
! You integrate on r2 the analytical integral over r1 of f_HF(r1,r2)
do ipoint = 1, n_points_final_grid
weight1 =final_weight_at_r_vector(ipoint)
r2(1) = final_grid_points(1,ipoint)
r2(2) = final_grid_points(2,ipoint)
r2(3) = final_grid_points(3,ipoint)
call integral_f_HF_valence_ab(r2,integral_psi_val)
accu_val += integral_psi_val * weight1
enddo
print*,'**************************'
! Should give you the alpha-beta repulsion of HF, excluding core contributions,
print*,'int dr1 dr2 f_HF(r1,r2) = ',accu_val
double precision :: accu_2
print*,'pure numerical form (might take quite some time as it grows as N_g^2 * N_e^2 * N_b^2 ...)'
! You integrate brut force on r1 and r2
accu_val = 0.d0
do jpoint = 1, n_points_final_grid
weight1 =final_weight_at_r_vector(jpoint)
r1(1) = final_grid_points(1,jpoint)
r1(2) = final_grid_points(2,jpoint)
r1(3) = final_grid_points(3,jpoint)
do ipoint = 1, n_points_final_grid
weight2 =final_weight_at_r_vector(ipoint)
r2(1) = final_grid_points(1,ipoint)
r2(2) = final_grid_points(2,ipoint)
r2(3) = final_grid_points(3,ipoint)
call f_HF_valence_ab(r1,r2,integral_psi_val,two_bod)
accu_val += integral_psi_val * weight1 * weight2
enddo
enddo
print*,'int dr1 dr2 f_HF(r1,r2) = ',accu_val
print*,'**************************'
print*,'**************************'
print*,'**************************'
accu_val = 0.d0
r1 = 0.d0
r1(1) = 0.5d0
print*,'r1 = ',r1
! You compute the integral over r2 of f_HF(r1,r2)
call integral_f_HF_valence_ab(r1,integral_psi)
do ipoint = 1, n_points_final_grid
weight1 =final_weight_at_r_vector(ipoint)
r2(1) = final_grid_points(1,ipoint)
r2(2) = final_grid_points(2,ipoint)
r2(3) = final_grid_points(3,ipoint)
call f_HF_valence_ab(r1,r2,integral_psi_val,two_bod)
accu_val += integral_psi_val * weight1
enddo
print*,'int dr2 f_HF(r1,r2) = ',integral_psi
print*,'analytical form = ',accu_val
print*,'**************************'
end
subroutine test_f_ii_valence_ab
implicit none
BEGIN_DOC
! routine to test the function f_ii(r1,r2)
!
! it should be the same that f_HF(r1,r2) only for inactive orbitals
END_DOC
integer :: ipoint
double precision :: accu_f, accu_n2, weight, r1(3),r2(3)
double precision :: accu_f_on_top
double precision :: f_HF_val_ab,two_bod_dens_hf,f_ii_val_ab,two_bod_dens_ii
accu_f = 0.d0
accu_n2 = 0.d0
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
r2 = r1
call f_HF_valence_ab(r1,r2,f_HF_val_ab,two_bod_dens_hf)
call give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens_ii)
accu_f += dabs(f_HF_val_ab - f_ii_val_ab) * weight
accu_n2+= dabs(two_bod_dens_hf - two_bod_dens_ii) * weight
accu_f_on_top += dabs(two_bod_dens_hf) * weight
enddo
print*,'**************************'
print*,''
print*,'accu_f = ',accu_f
print*,'accu_n2 = ',accu_n2
print*,''
print*,'accu_f_on_top = ',accu_f_on_top
end
subroutine test_f_ia_valence_ab
implicit none
BEGIN_DOC
! routine to test the function f_ii(r1,r2), f_ia(r1,r2) and f_aa(r1,r2)
END_DOC
integer :: ipoint,istate
double precision :: accu_f, accu_n2, weight, r1(3),r2(3)
double precision :: accu_f_on_top
double precision :: f_ref,f_comp,on_top_ref,on_top_comp
double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa
double precision :: accu
accu_f = 0.d0
accu_n2 = 0.d0
accu = 0.d0
istate = 1
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
r1(1) = final_grid_points(1,ipoint)
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
r2 = r1
call give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens_ii)
call give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens_ia,istate)
call give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens_aa,istate)
f_ref = f_psi_cas_ab_old(ipoint,istate)
f_comp = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab
on_top_ref = total_cas_on_top_density(ipoint,istate)
on_top_comp= two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
accu_f += dabs(f_ref - f_comp) * weight
accu_n2+= dabs(on_top_ref - on_top_comp) * weight
accu += f_ref * weight
enddo
print*,'**************************'
print*,''
print*,'accu_f = ',accu_f
print*,'accu_n2 = ',accu_n2
print*,''
print*,'accu = ',accu
end
subroutine test_f_ii_ia_aa_valence_ab
implicit none
BEGIN_DOC
! routine to test the function f_Psi(r1,r2) based on core/inactive/active orbitals
END_DOC
integer :: ipoint,istate
double precision :: accu_f, accu_n2, weight, r1(3),r2(3)
double precision :: accu_f_on_top
double precision :: f_ref,f_comp,on_top_ref,on_top_comp
double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa
double precision :: accu
accu_f = 0.d0
accu_n2 = 0.d0
accu = 0.d0
istate = 1
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
f_ref = f_psi_cas_ab(ipoint,istate)
f_comp = f_psi_cas_ab_old(ipoint,istate)
on_top_ref = total_cas_on_top_density(ipoint,istate)
on_top_comp= on_top_cas_mu_r(ipoint,istate)
accu_f += dabs(f_ref - f_comp) * weight
accu_n2+= dabs(on_top_ref - on_top_comp) * weight
accu += f_ref * weight
enddo
print*,'**************************'
print*,''
print*,'accu_f = ',accu_f
print*,'accu_n2 = ',accu_n2
print*,''
print*,'accu = ',accu
end

138
src/mu_of_r/f_hf_utils.irp.f Normal file
View File

@ -0,0 +1,138 @@
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
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
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
end

313
src/mu_of_r/f_psi_i_a_v_utils.irp.f Normal file
View File

@ -0,0 +1,313 @@
subroutine give_f_ii_val_ab(r1,r2,f_ii_val_ab,two_bod_dens)
implicit none
BEGIN_DOC
! contribution from purely inactive orbitals to f_{\Psi^B}(r_1,r_2) for a CAS wave function
END_DOC
double precision, intent(in) :: r1(3),r2(3)
double precision, intent(out):: f_ii_val_ab,two_bod_dens
integer :: i,j,m,n,i_m,i_n
integer :: i_i,i_j
double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
double precision :: get_two_e_integral
! You get all orbitals in r_1 and r_2
allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
call give_all_mos_at_r(r1,mos_array_r1)
call give_all_mos_at_r(r2,mos_array_r2)
! You extract the inactive orbitals
allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
do i_m = 1, n_inact_orb
mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m))
enddo
do i_m = 1, n_inact_orb
mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m))
enddo
! You extract the orbitals belonging to the space \mathcal{B}
allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) )
do i_m = 1, n_basis_orb
mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m))
mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m))
enddo
f_ii_val_ab = 0.d0
two_bod_dens = 0.d0
! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space
do m = 1, n_inact_orb ! electron 1
! You browse all OCCUPIED BETA electrons in the \mathcal{A} space
do n = 1, n_inact_orb ! electron 2
! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2)
two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_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_ii_val_ab += two_e_int_ii_f(j,i,n,m) * mos_array_inact_r1(m) * mos_array_basis_r1(i) &
* mos_array_inact_r2(n) * mos_array_basis_r2(j)
enddo
enddo
enddo
enddo
end
subroutine give_f_ia_val_ab(r1,r2,f_ia_val_ab,two_bod_dens,istate)
BEGIN_DOC
! contribution from inactive and active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function
END_DOC
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: r1(3),r2(3)
double precision, intent(out):: f_ia_val_ab,two_bod_dens
integer :: i,orb_i,a,orb_a,n,m,b
double precision :: rho
double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:)
double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:)
f_ia_val_ab = 0.d0
two_bod_dens = 0.d0
! You get all orbitals in r_1 and r_2
allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
call give_all_mos_at_r(r1,mos_array_r1)
call give_all_mos_at_r(r2,mos_array_r2)
! You extract the inactive orbitals
allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) )
do i = 1, n_inact_orb
mos_array_inact_r1(i) = mos_array_r1(list_inact(i))
enddo
do i= 1, n_inact_orb
mos_array_inact_r2(i) = mos_array_r2(list_inact(i))
enddo
! You extract the active orbitals
allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) )
do i= 1, n_act_orb
mos_array_act_r1(i) = mos_array_r1(list_act(i))
enddo
do i= 1, n_act_orb
mos_array_act_r2(i) = mos_array_r2(list_act(i))
enddo
! You extract the orbitals belonging to the space \mathcal{B}
allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) )
do i= 1, n_basis_orb
mos_array_basis_r1(i) = mos_array_r1(list_basis(i))
enddo
do i= 1, n_basis_orb
mos_array_basis_r2(i) = mos_array_r2(list_basis(i))
enddo
! Contracted density : intermediate quantity
! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2)
allocate(rho_tilde(n_inact_orb,n_act_orb))
two_bod_dens = 0.d0
do a = 1, n_act_orb
do i = 1, n_inact_orb
rho_tilde(i,a) = 0.d0
do b = 1, n_act_orb
rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate)
two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho
rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b)
enddo
enddo
enddo
! Contracted two-e integrals : intermediate quantity
! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n >
allocate( v_tilde(n_act_orb,n_act_orb) )
allocate( integrals_array(mo_num,mo_num) )
v_tilde = 0.d0
do a = 1, n_act_orb
orb_a = list_act(a)
do i = 1, n_inact_orb
v_tilde(i,a) = 0.d0
orb_i = list_inact(i)
! call get_mo_two_e_integrals_ij(orb_i,orb_a,mo_num,integrals_array,mo_integrals_map)
do m = 1, n_basis_orb
do n = 1, n_basis_orb
! v_tilde(i,a) += integrals_array(n,m) * mos_array_basis_r2(n) * mos_array_basis_r1(m)
v_tilde(i,a) += two_e_int_ia_f(n,m,i,a) * mos_array_basis_r2(n) * mos_array_basis_r1(m)
enddo
enddo
enddo
enddo
do a = 1, n_act_orb
do i = 1, n_inact_orb
f_ia_val_ab += v_tilde(i,a) * rho_tilde(i,a)
enddo
enddo
end
subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
BEGIN_DOC
! contribution from purely active orbitals to f_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function
END_DOC
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: r1(3),r2(3)
double precision, intent(out):: f_aa_val_ab,two_bod_dens
integer :: i,orb_i,a,orb_a,n,m,b,c,d
double precision :: rho
double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:)
double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:)
double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:)
double precision, allocatable :: integrals_array(:,:),rho_tilde(:,:),v_tilde(:,:)
f_aa_val_ab = 0.d0
two_bod_dens = 0.d0
! You get all orbitals in r_1 and r_2
allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) )
call give_all_mos_at_r(r1,mos_array_r1)
call give_all_mos_at_r(r2,mos_array_r2)
! You extract the active orbitals
allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) )
do i= 1, n_act_orb
mos_array_act_r1(i) = mos_array_r1(list_act(i))
enddo
do i= 1, n_act_orb
mos_array_act_r2(i) = mos_array_r2(list_act(i))
enddo
! You extract the orbitals belonging to the space \mathcal{B}
allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) )
do i= 1, n_basis_orb
mos_array_basis_r1(i) = mos_array_r1(list_basis(i))
enddo
do i= 1, n_basis_orb
mos_array_basis_r2(i) = mos_array_r2(list_basis(i))
enddo
! Contracted density : intermediate quantity
! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2)
allocate(rho_tilde(n_act_orb,n_act_orb))
two_bod_dens = 0.d0
rho_tilde = 0.d0
do a = 1, n_act_orb ! 1
do b = 1, n_act_orb ! 2
do c = 1, n_act_orb ! 1
do d = 1, n_act_orb ! 2
rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
rho_tilde(b,a) += rho
two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b)
enddo
enddo
enddo
enddo
! Contracted two-e integrals : intermediate quantity
! v_tilde(i,a) = \sum_{m,n} phi_m(1) * phi_n(2) < i a | m n >
allocate( v_tilde(n_act_orb,n_act_orb) )
v_tilde = 0.d0
do a = 1, n_act_orb
do b = 1, n_act_orb
v_tilde(b,a) = 0.d0
do m = 1, n_basis_orb
do n = 1, n_basis_orb
v_tilde(b,a) += two_e_int_aa_f(n,m,b,a) * mos_array_basis_r2(n) * mos_array_basis_r1(m)
enddo
enddo
enddo
enddo
do a = 1, n_act_orb
do b = 1, n_act_orb
f_aa_val_ab += v_tilde(b,a) * rho_tilde(b,a)
enddo
enddo
end
BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act_orb,n_act_orb)]
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_{ii}(r_1,r_2)
!
! two_e_int_aa_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 :: integrals_array(mo_num,mo_num),get_two_e_integral
do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m)
do orb_n = 1, n_act_orb ! electron 2
n = list_act(orb_n)
call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
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_aa_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
! two_e_int_aa_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_act_orb)]
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_{ia}(r_1,r_2)
!
! two_e_int_aa_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 :: integrals_array(mo_num,mo_num),get_two_e_integral
do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m)
do orb_n = 1, n_inact_orb ! electron 2
n = list_inact(orb_n)
call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
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_ia_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
two_e_int_ia_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_inact_orb,n_inact_orb)]
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_{ii}(r_1,r_2)
!
! two_e_int_ii_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,integrals_array(mo_num,mo_num)
do orb_m = 1, n_inact_orb ! electron 1
m = list_inact(orb_m)
do orb_n = 1, n_inact_orb ! electron 2
n = list_inact(orb_n)
call get_mo_two_e_integrals_ij(m,n,mo_num,integrals_array,mo_integrals_map)
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_ii_f(orb_j,orb_i,orb_n,orb_m) = get_two_e_integral(m,n,i,j,mo_integrals_map)
two_e_int_ii_f(orb_j,orb_i,orb_n,orb_m) = integrals_array(j,i)
enddo
enddo
enddo
enddo
END_PROVIDER

39
src/mu_of_r/f_psi_old.irp.f Normal file
View File

@ -0,0 +1,39 @@
BEGIN_PROVIDER [double precision, f_psi_cas_ab_old, (n_points_final_grid,N_states)]
implicit none
BEGIN_DOC
!
! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid and for all states
!
! Assumes that the wave function in psi_det is developped within an active space defined
!
END_DOC
integer :: ipoint,k,l,istate
double precision :: wall0,wall1
print*,'Providing f_psi_cas_ab_old ..... '
provide full_occ_2_rdm_ab_mo
call wall_time(wall0)
provide core_inact_act_V_kl_contracted full_occ_2_rdm_cntrctd
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,k,l,istate) &
!$OMP SHARED (n_core_inact_act_orb, n_points_final_grid, full_occ_2_rdm_cntrctd, core_inact_act_V_kl_contracted, f_psi_cas_ab_old,N_states)
!$OMP DO
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
f_psi_cas_ab_old(ipoint,istate) = 0.d0
do l = 1, n_core_inact_act_orb ! 2
do k = 1, n_core_inact_act_orb ! 1
f_psi_cas_ab_old(ipoint,istate) += core_inact_act_V_kl_contracted(k,l,ipoint) * full_occ_2_rdm_cntrctd(k,l,ipoint,istate)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'Time to provide f_psi_cas_ab_old = ',wall1 - wall0
END_PROVIDER

151
src/mu_of_r/f_psi_utils.irp.f Normal file
View File

@ -0,0 +1,151 @@
BEGIN_PROVIDER [double precision, core_inact_act_V_kl_contracted, (n_core_inact_act_orb,n_core_inact_act_orb,n_points_final_grid)]
implicit none
BEGIN_DOC
! core_inact_act_V_kl_contracted(k,l,ipoint) = \sum_{ij} V_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint)
!
! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
END_DOC
integer :: ipoint,k,l
do k = 1, n_core_inact_act_orb
do l = 1, n_core_inact_act_orb
do ipoint = 1, n_points_final_grid
core_inact_act_V_kl_contracted(k,l,ipoint) = full_occ_v_kl_cntrctd(ipoint,k,l)
enddo
enddo
enddo
free full_occ_v_kl_cntrctd
END_PROVIDER
BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd, (n_core_inact_act_orb,n_core_inact_act_orb,n_points_final_grid,N_states)]
implicit none
BEGIN_DOC
! full_occ_2_rdm_cntrctd(k,l,ipoint,istate) = \sum_{ij} \Gamma_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint)
!
! where \Gamma_{ij}^{kl}(istate) = <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
!
! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
END_DOC
integer :: ipoint,k,l,istate
do istate = 1, N_states
do k = 1, n_core_inact_act_orb
do l = 1, n_core_inact_act_orb
do ipoint = 1, n_points_final_grid
full_occ_2_rdm_cntrctd(k,l,ipoint,istate) = full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate)
enddo
enddo
enddo
enddo
free full_occ_2_rdm_cntrctd_trans
END_PROVIDER
BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_grid,n_core_inact_act_orb,n_core_inact_act_orb,N_states)]
implicit none
BEGIN_DOC
! full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) = \sum_{ij} \Gamma_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint)
!
! where \Gamma_{ij}^{kl}(istate) = <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
!
! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
END_DOC
integer :: i,j,k,l,istate
integer :: ipoint
double precision, allocatable :: mos_array_r(:),r(:)
provide full_occ_2_rdm_ab_mo
double precision :: wall0,wall1
print*,'Providing full_occ_2_rdm_cntrctd_trans ..... '
call wall_time(wall0)
full_occ_2_rdm_cntrctd_trans = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,k,l,i,j,istate) &
!$OMP SHARED (n_core_inact_act_orb, n_points_final_grid, full_occ_2_rdm_cntrctd_trans, final_grid_points,full_occ_2_rdm_ab_mo,core_inact_act_mos_in_r_array,N_states )
!$OMP DO
do istate = 1, N_states
do l = 1, n_core_inact_act_orb ! 2
do k = 1, n_core_inact_act_orb ! 1
do ipoint = 1, n_points_final_grid
do j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb
! 1 2 1 2
full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'Time to provide full_occ_2_rdm_cntrctd_trans = ',wall1 - wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, full_occ_v_kl_cntrctd, (n_points_final_grid,n_core_inact_act_orb,n_core_inact_act_orb)]
implicit none
BEGIN_DOC
! full_occ_v_kl_cntrctd(ipoint,k,l) = \sum_{ij} V_{ij}^{kl} phi_i(r_ipoint) phi_j(r_ipoint)
!
! This is needed to build the function f_{\Psi^B}(X_1,X_2) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018)
!
END_DOC
integer :: i,j,k,l,kk,ll,ii,jj
integer :: ipoint
double precision, allocatable :: integrals_array(:,:), mos_array_r(:),r(:), integrals_basis(:,:)
! just not to mess with parallelization
allocate(integrals_array(mo_num,mo_num))
k = 1
l = 1
call get_mo_two_e_integrals_ij(k,l,mo_num,integrals_array,mo_integrals_map)
provide basis_mos_in_r_array
deallocate(integrals_array)
double precision :: wall0,wall1
call wall_time(wall0)
full_occ_v_kl_cntrctd = 0.d0
print*,'Providing full_occ_v_kl_cntrctd ..... '
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,kk,ll,k,l,i,j,ii,jj,integrals_array,integrals_basis) &
!$OMP SHARED (mo_num, n_points_final_grid, n_basis_orb, list_basis, full_occ_v_kl_cntrctd, mo_integrals_map,final_grid_points,basis_mos_in_r_array, n_core_inact_act_orb, list_core_inact_act)
allocate(integrals_array(mo_num,mo_num), integrals_basis(n_basis_orb,n_basis_orb))
!$OMP DO
do l = 1, n_core_inact_act_orb! 2
ll = list_core_inact_act(l)
do k = 1, n_core_inact_act_orb ! 1
kk = list_core_inact_act(k)
call get_mo_two_e_integrals_ij(kk,ll,mo_num,integrals_array,mo_integrals_map)
do j = 1, n_basis_orb
jj = list_basis(j)
do i = 1, n_basis_orb
ii = list_basis(i)
integrals_basis(i,j) = integrals_array(ii,jj)
enddo
enddo
do ipoint = 1, n_points_final_grid
do j = 1, n_basis_orb ! condition on mo_num in order to ensure the correct CBS limit
do i = 1, n_basis_orb !
!1 2 1 2
full_occ_v_kl_cntrctd(ipoint,k,l) += integrals_basis(i,j) * basis_mos_in_r_array(j,ipoint) * basis_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
deallocate(integrals_array,integrals_basis)
!$OMP END PARALLEL
call wall_time(wall1)
print*,'Time to provide full_occ_v_kl_cntrctd = ',wall1 - wall0
END_PROVIDER

90
src/mu_of_r/f_val_general.irp.f Normal file
View File

@ -0,0 +1,90 @@
BEGIN_PROVIDER [double precision, f_psi_cas_ab, (n_points_final_grid,N_states)]
&BEGIN_PROVIDER [double precision, on_top_cas_mu_r, (n_points_final_grid,N_states)]
implicit none
BEGIN_DOC
!
! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid and for all states and for a CAS wave function
!
! Assumes that the wave function in psi_det is developped within an active space defined
!
END_DOC
integer :: ipoint,istate
double precision :: wall0,wall1,r(3)
double precision :: f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa
double precision :: accu
accu = 0.d0
r = 0.d0
istate = 1
! To initialize parallelization
call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii)
call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate)
call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate)
provide final_grid_points act_2_rdm_ab_mo
print*,'Providing f_psi_cas_ab..... '
call wall_time(wall0)
do istate = 1, N_states
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,r,f_ii_val_ab,two_bod_dens_ii,f_ia_val_ab,two_bod_dens_ia,f_aa_val_ab,two_bod_dens_aa) &
!$OMP SHARED (n_points_final_grid,f_psi_cas_ab,on_top_cas_mu_r,final_grid_points,istate)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
! inactive-inactive part of f_psi(r1,r2)
call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii)
! inactive-active part of f_psi(r1,r2)
call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate)
! active-active part of f_psi(r1,r2)
call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate)
f_psi_cas_ab(ipoint,istate) = f_ii_val_ab + f_ia_val_ab + f_aa_val_ab
on_top_cas_mu_r(ipoint,istate) = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
enddo
!$OMP END DO
!$OMP END PARALLEL
enddo
call wall_time(wall1)
print*,'Time to provide f_psi_cas_ab = ',wall1 - wall0
print*,'accu = ',accu
END_PROVIDER
BEGIN_PROVIDER [double precision, f_psi_hf_ab, (n_points_final_grid)]
&BEGIN_PROVIDER [double precision, on_top_hf_mu_r, (n_points_final_grid)]
implicit none
BEGIN_DOC
!
! Function f_{\Psi^B}(r,r) of Eq. (22) of J. Chem. Phys. 149, 194301 (2018) on each point of the grid for a HF wave function
!
END_DOC
integer :: ipoint
double precision :: wall0,wall1,r(3),f_HF_val_ab,two_bod_dens
f_psi_hf_ab = 0.d0
r = 0.d0
! To initialize parallelization
call f_HF_valence_ab(r,r,f_HF_val_ab,two_bod_dens)
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,r,f_HF_val_ab,two_bod_dens) &
!$OMP SHARED (n_points_final_grid,f_psi_hf_ab,on_top_hf_mu_r,final_grid_points)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
call f_HF_valence_ab(r,r,f_HF_val_ab,two_bod_dens)
f_psi_hf_ab(ipoint) = f_HF_val_ab
on_top_hf_mu_r(ipoint) = two_bod_dens
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'Time to provide f_psi_hf_ab = ',wall1 - wall0
END_PROVIDER

150
src/mu_of_r/mu_of_r_conditions.irp.f Normal file
View File

@ -0,0 +1,150 @@
BEGIN_PROVIDER [double precision, mu_of_r_prov, (n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! general variable for mu(r)
!
! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the two-body density matrix are excluded
END_DOC
integer :: ipoint,istate
double precision :: wall0,wall1
print*,'providing mu_of_r ...'
! PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
call wall_time(wall0)
if (read_mu_of_r) then
print*,'Reading mu(r) from disk ...'
call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_prov)
return
endif
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
if(mu_of_r_potential.EQ."hf")then
mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint)
else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate)
else
print*,'you requested the following mu_of_r_potential'
print*,mu_of_r_potential
print*,'which does not correspond to any of the options for such keyword'
stop
endif
enddo
enddo
if (write_mu_of_r) then
print*,'Writing mu(r) on disk ...'
call ezfio_set_mu_of_r_io_mu_of_r('Read')
call ezfio_set_mu_of_r_mu_of_r_disk(mu_of_r_prov)
endif
call wall_time(wall1)
print*,'Time to provide mu_of_r = ',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_of_r_hf, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO)
!
! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the two-body density matrix are excluded
END_DOC
integer :: ipoint
double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
print*,'providing mu_of_r_hf ...'
call wall_time(wall0)
sqpi = dsqrt(dacos(-1.d0))
provide f_psi_hf_ab
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) &
!$OMP ShARED (n_points_final_grid,mu_of_r_hf,f_psi_hf_ab,on_top_hf_mu_r,sqpi)
do ipoint = 1, n_points_final_grid
f_hf = f_psi_hf_ab(ipoint)
on_top = on_top_hf_mu_r(ipoint)
if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
w_hf = 1.d+10
else
w_hf = f_hf / on_top
endif
mu_of_r_hf(ipoint) = w_hf * sqpi * 0.5d0
enddo
!$OMP END PARALLEL DO
call wall_time(wall1)
print*,'Time to provide mu_of_r_hf = ',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_of_r_psi_cas, (n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! mu(r) computed with a wave function developped in an active space
!
! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the one- and two-body density matrix are excluded
END_DOC
integer :: ipoint,istate
double precision :: wall0,wall1,f_psi,on_top,w_psi,sqpi
print*,'providing mu_of_r_psi_cas ...'
call wall_time(wall0)
sqpi = dsqrt(dacos(-1.d0))
provide f_psi_cas_ab
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) &
!$OMP SHARED (n_points_final_grid,mu_of_r_psi_cas,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
f_psi = f_psi_cas_ab(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
w_psi = 1.d+10
else
w_psi = f_psi / on_top
endif
mu_of_r_psi_cas(ipoint,istate) = w_psi * sqpi * 0.5d0
enddo
enddo
!$OMP END PARALLEL DO
call wall_time(wall1)
print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
implicit none
BEGIN_DOC
! average value of mu(r) weighted with the total one-e density and divised by the number of electrons
!
! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals
!
! in the one- and two-body density matrix are excluded
END_DOC
integer :: ipoint,istate
double precision :: weight,density
mu_average_prov = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight =final_weight_at_r_vector(ipoint)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) &
+ one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle
mu_average_prov(istate) += mu_of_r_prov(ipoint,istate) * weight * density
enddo
mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
enddo
END_PROVIDER

21
src/mu_of_r/test_proj_op.irp.f Normal file
View File

@ -0,0 +1,21 @@
program projected_operators
implicit none
BEGIN_DOC
! TODO
END_DOC
read_wf = .True.
touch read_wf
! You specify that you want to avoid any contribution from
! orbitals coming from core
no_core_density = .True.
touch no_core_density
mu_of_r_potential = "cas_ful"
touch mu_of_r_potential
print*,'Using Valence Only functions'
! call test_f_HF_valence_ab
! call routine_full_mos
! call test_f_ii_valence_ab
call test_f_ia_valence_ab
call test_f_ii_ia_aa_valence_ab
end