PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 15:12:19 +02:00
Code Releases Activity

Merge pull request #99 from QuantumPackage/cleaning_dft

Browse Source 
Cleaning dft
This commit is contained in:
Anthony Scemama 2020-04-21 15:07:36 +02:00 committed by GitHub
commit 1362042e23
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
69 changed files with 5818 additions and 164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/README.rst
View File

@ -22,8 +22,8 @@ The core modules of the QP
Ex : if "exchange_functional" == "sr_pbe", then energy_x will contain the exchange correlation functional defined in "functiona/sr_pbe.irp.f", which corresponds to the short-range PBE functional (at the value mu_erf for the range separation parameter)
*** How are handled the DFT functionals in QP2 ?
================================================
*** How to add a new functional in QP2
======================================
Creating a new functional and propagating it through the whole QP2 programs is easy as all dependencies are handled by a script.

63
src/basis_correction/51.basis_c.bats Normal file
View File

@ -0,0 +1,63 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
source $QP_ROOT/quantum_package.rc
function run() {
thresh=$2
test_exe fci || skip
qp edit --check
qp set perturbation do_pt2 False
qp set determinants n_det_max 8000
qp set determinants n_states 1
qp set davidson threshold_davidson 1.e-10
qp set davidson n_states_diag 8
qp run fci
energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
eq $energy1 $1 $thresh
}
function run_md() {
thresh=$2
qp set mu_of_r mu_of_r_potential cas_ful
file_out=${EZFIO_FILE}.basis_corr.out
qp run basis_correction | tee $file_out
energy1="$(grep 'ECMD SU-PBE-OT , state 1 =' ${file_out} | cut -d '=' -f 2)"
eq $energy1 $1 $thresh
}
function run_sd() {
thresh=$2
qp set mu_of_r mu_of_r_potential hf
qp set_frozen_core
file_out=${EZFIO_FILE}.basis_corr.out
qp run basis_correction | tee $file_out
energy1="$(grep 'ECMD PBE-UEG , state 1 =' ${file_out} | cut -d '=' -f 2)"
eq $energy1 $1 $thresh
}
@test "O2 CAS" {
qp set_file o2_cas.gms.ezfio
qp set_mo_class -c "[1-2]" -a "[3-10]" -d "[11-46]"
run -149.72435425 3.e-4 10000
qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
run_md -0.1160222327 1.e-6
}
@test "LiF RHF" {
qp set_file lif.ezfio
run_sd -0.0649431665 1.e-6
}
@test "F ROHF" {
qp set_file f.ezfio
run_sd -0.0355395041 1.e-6
}
@test "Be RHF" {
qp set_file be.ezfio
run_sd -0.0325139011 1.e-6
}

3
src/basis_correction/NEED Normal file
View File

@ -0,0 +1,3 @@
mu_of_r
dft_utils_func
dft_one_e

27
src/basis_correction/README.rst Normal file
View File

@ -0,0 +1,27 @@
================
basis_correction
================
This module proposes the various flavours of the DFT-based basis set correction originally proposed in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714.
This basis set correction relies mainy on :
+) The definition of a range-separation function \mu(r) varying in space to mimic the incompleteness of the basis set used to represent the coulomb interaction. This procedure needs a two-body rdm representing qualitatively the spacial distribution of the opposite spin electron pairs.
Two types of \mu(r) are proposed, according to the strength of correlation, through the keyword "mu_of_r_potential" in the module "mu_of_r":
a) "mu_of_r_potential = hf" uses the two-body rdm of a HF-like wave function (i.e. a single Slater determinant developped with the MOs stored in the EZFIO folder).
When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK.
See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results.
Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals.
b) "mu_of_r_potential = cas_ful" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder).
If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation.
+) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor. Chem. Acc.114, 305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density.
The two main flavours of Ecmd depends on the strength of correlation in the system:
a) for weakly correlated systems, the ECMD PBE-UEG functional based on the seminal work of in RSDFT (see JCP, 150, 084103 1-10 (2019)) and adapted for the basis set correction in JPCL, 10, 2931-2937 (2019) uses the exact on-top pair density of the UEG at large mu and the PBE correlation functional at mu = 0. As shown in JPCL, 10, 2931-2937 (2019), such a functional is more accurate than the ECMD LDA for weakly correlated systems.
b) for strongly correlated systems, the ECMD PBE-OT, which uses the extrapolated on-top pair density of the CAS wave function thanks to the large \mu behaviour of the on-top pair density, is accurate, but suffers from S_z dependence (i.e. is not invariant with respect to S_z) because of the spin-polarization dependence of the PBE correlation functional entering at mu=0.
An alternative is ECMD SU-PBE-OT which uses the same on-top pair density that ECMD PBE-OT but a ZERO spin-polarization to remove the S_z dependence. As shown in ???????????, this strategy is one of the more accurate and respects S_z invariance and size consistency if the CAS wave function is correctly chosen.

1
src/basis_correction/TODO Normal file
View File

@ -0,0 +1 @@
change all correlation functionals with the pbe_on_top general

27
src/basis_correction/basis_correction.irp.f Normal file
View File

@ -0,0 +1,27 @@
program basis_correction
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
read_wf = .True.
touch read_wf
no_core_density = .True.
touch no_core_density
provide mo_two_e_integrals_in_map
call print_basis_correction
! call print_e_b
end
subroutine print_e_b
implicit none
print *, 'Hello world'
print*,'ecmd_lda_mu_of_r = ',ecmd_lda_mu_of_r
print*,'ecmd_pbe_ueg_mu_of_r = ',ecmd_pbe_ueg_mu_of_r
print*,'ecmd_pbe_ueg_eff_xi_mu_of_r = ',ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'psi_energy + E^B_LDA = ',psi_energy + ecmd_lda_mu_of_r
print*,'psi_energy + E^B_PBE_UEG = ',psi_energy + ecmd_pbe_ueg_mu_of_r
print*,'psi_energy + E^B_PBE_UEG_Xi = ',psi_energy + ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'mu_average_prov = ',mu_average_prov
end

92
src/basis_correction/eff_xi_based_func.irp.f Normal file
View File

@ -0,0 +1,92 @@
BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_eff_xi_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_pbe_ueg_eff_xi_mu_of_r = multi-determinantal Ecmd within the PBE-UEG and effective spin polarization approximation with mu(r),
!
! see Eqs. 30 in ???????????
!
! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638)
!
! and replaces the approximation of the exact on-top pair density by the exact on-top of the UEG
!
! !!!! BUT !!!! with an EFFECTIVE SPIN POLARIZATION DEPENDING ON THE ON-TOP PAIR DENSITY
!
! See P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531 (1995). for original Ref., and Eq. 29 in ???????????
END_DOC
implicit none
double precision :: weight,density
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),g0_UEG_mu_inf,on_top
ecmd_pbe_ueg_eff_xi_mu_of_r = 0.d0
print*,'Providing ecmd_pbe_ueg_eff_xi_mu_of_r ...'
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight=final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! We use the effective spin density to define rho_a/rho_b
rho_a = 0.5d0 * (density + effective_spin_dm(ipoint,istate))
rho_b = 0.5d0 * (density - effective_spin_dm(ipoint,istate))
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
! with the effective rho_a and rho_b
on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE)
ecmd_pbe_ueg_eff_xi_mu_of_r(istate) += eps_c_md_PBE * weight
enddo
enddo
double precision :: wall1, wall0
call wall_time(wall1)
print*,'Time for the ecmd_pbe_ueg_eff_xi_mu_of_r:',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_lda_eff_xi_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_lda_eff_xi_mu_of_r = multi-determinantal Ecmd within the LDA and effective spin polarization approximation with mu(r),
!
! corresponds to equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714
!
! !!!! BUT !!!! with an EFFECTIVE SPIN POLARIZATION DEPENDING ON THE ON-TOP PAIR DENSITY
!
! See P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531 (1995). for original Ref., and Eq. 29 in ???????????
END_DOC
implicit none
integer :: ipoint,istate
double precision :: rho_a, rho_b, ec
logical :: dospin
double precision :: wall0,wall1,weight,mu,density
dospin = .true. ! JT dospin have to be set to true for open shell
print*,'Providing ecmd_lda_eff_xi_mu_of_r ...'
ecmd_lda_eff_xi_mu_of_r = 0.d0
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
mu = mu_of_r_prov(ipoint,istate)
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)
rho_a = 0.5d0 * (density + effective_spin_dm(ipoint,istate))
rho_b = 0.5d0 * (density - effective_spin_dm(ipoint,istate))
call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
if(isnan(ec))then
print*,'ec is nan'
stop
endif
ecmd_lda_eff_xi_mu_of_r(istate) += weight * ec
enddo
enddo
call wall_time(wall1)
print*,'Time for ecmd_lda_eff_xi_mu_of_r :',wall1-wall0
END_PROVIDER

183
src/basis_correction/pbe_on_top.irp.f Normal file
View File

@ -0,0 +1,183 @@
BEGIN_PROVIDER [double precision, ecmd_pbe_on_top_mu_of_r, (N_states)]
BEGIN_DOC
!
! Ecmd functional evaluated with mu(r) and depending on
! +) the on-top pair density
!
! +) the total density, density gradients
!
! +) the spin density
!
! Defined originally in Eq. (25) of JCP, 150, 084103 1-10 (2019) for RS-DFT calculations, but evaluated with mu(r).
!
! Such a functional is built by interpolating between two regimes :
!
! +) the large mu behaviour in cst/(\mu^3) \int dr on-top(r) where on-top(r) is supposed to be the exact on-top of the system
!
! +) mu= 0 with the usal ec_pbe(rho_a,rho_b,grad_rho_a,grad_rho_b)
!
! Here the approximation to the exact on-top is done through the assymptotic expansion (in \mu) of the exact on-top pair density (see Eq. 29) but with a mu(r) instead of a constant mu
!
! Such an asymptotic expansion was introduced in P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
END_DOC
implicit none
double precision :: weight
double precision :: eps_c_md_on_top_PBE,on_top_extrap,mu_correction_of_on_top
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
ecmd_pbe_on_top_mu_of_r = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
! depends on (rho_a, rho_b) <==> (rho_tot,spin_pol)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
ecmd_pbe_on_top_mu_of_r(istate) += eps_c_md_on_top_PBE * weight
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_on_top_su_mu_of_r, (N_states)]
BEGIN_DOC
!
! Ecmd functional evaluated with mu(r) and depending on
! +) the on-top pair density
!
! +) the total density, density gradients
!
! +) !!!!! NO SPIN POLAIRIZATION !!!!!
!
! Defined originally in Eq. (25) of JCP, 150, 084103 1-10 (2019) for RS-DFT calculations, but evaluated with mu(r).
!
! Such a functional is built by interpolating between two regimes :
!
! +) the large mu behaviour in cst/(\mu^3) \int dr on-top(r) where on-top(r) is supposed to be the exact on-top of the system
!
! +) mu= 0 with the usal ec_pbe(rho_a,rho_b,grad_rho_a,grad_rho_b)
!
! Here the approximation to the exact on-top is done through the assymptotic expansion (in \mu) of the exact on-top pair density (see Eq. 29) but with a mu(r) instead of a constant mu
!
! Such an asymptotic expansion was introduced in P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
END_DOC
implicit none
double precision :: weight
double precision :: eps_c_md_on_top_PBE,on_top_extrap,mu_correction_of_on_top
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top,density
ecmd_pbe_on_top_su_mu_of_r = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! rho_a = rho_b = rho_tot/2 ==> NO SPIN POLARIZATION
rho_a = 0.5d0 * density
rho_b = 0.5d0 * density
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
ecmd_pbe_on_top_su_mu_of_r(istate) += eps_c_md_on_top_PBE * weight
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_on_top_no_extrap_su_mu_of_r, (N_states)]
BEGIN_DOC
!
! Ecmd functional evaluated with mu(r) and depending on
! +) the on-top pair density
!
! +) the total density, density gradients
!
! +) !!!!! NO SPIN POLAIRIZATION !!!!!
!
! Defined originally in Eq. (25) of JCP, 150, 084103 1-10 (2019) for RS-DFT calculations, but evaluated with mu(r).
!
! Such a functional is built by interpolating between two regimes :
!
! +) the large mu behaviour in cst/(\mu^3) \int dr on-top(r) where on-top(r) is supposed to be the exact on-top of the system
!
! +) mu= 0 with the usal ec_pbe(rho_a,rho_b,grad_rho_a,grad_rho_b)
!
! Here the approximation to the exact on-top is done through the assymptotic expansion (in \mu) of the exact on-top pair density (see Eq. 29) but with a mu(r) instead of a constant mu
!
! Such an asymptotic expansion was introduced in P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
END_DOC
implicit none
double precision :: weight
double precision :: eps_c_md_on_top_PBE,on_top_extrap,mu_correction_of_on_top
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top,density
ecmd_pbe_on_top_no_extrap_su_mu_of_r = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mu = mu_of_r_prov(ipoint,istate)
density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! rho_a = rho_b = rho_tot/2 ==> NO SPIN POLARIZATION
rho_a = 0.5d0 * density
rho_b = 0.5d0 * density
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We DO NOT take the extrapolated on-top pair density, but there is * 2 because of normalization
on_top_extrap = 2.d0 * on_top
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
ecmd_pbe_on_top_no_extrap_su_mu_of_r(istate) += eps_c_md_on_top_PBE * weight
enddo
enddo
END_PROVIDER

82
src/basis_correction/print_routine.irp.f Normal file
View File

@ -0,0 +1,82 @@
subroutine print_basis_correction
implicit none
integer :: istate
provide mu_average_prov
if(mu_of_r_potential.EQ."hf")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r
endif
print*, ''
print*, ''
print*, '****************************************'
print*, '****************************************'
print*, 'Basis set correction for WFT using DFT Ecmd functionals'
print*, 'These functionals are accurate for short-range correlation'
print*, ''
print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018) '
print*, ' Journal of Physical Chemistry Letters 10, 2931-2937 (2019) '
print*, ' ???REF SC?'
print*, '****************************************'
print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential
if(mu_of_r_potential.EQ."hf")then
print*, ''
print*,'Using a HF-like two-body density to define mu(r)'
print*,'This assumes that HF is a qualitative representation of the wave function '
print*,'********************************************'
print*,'Functionals more suited for weak correlation'
print*,'********************************************'
print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
enddo
print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
else if(mu_of_r_potential.EQ."cas_ful")then
print*, ''
print*,'Using a CAS-like two-body density to define mu(r)'
print*,'This assumes that the CAS is a qualitative representation of the wave function '
print*,'********************************************'
print*,'Functionals more suited for weak correlation'
print*,'********************************************'
print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
enddo
print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'********************************************'
print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
enddo
print*,''
endif
print*,''
print*,'**************'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate)
enddo
end

83
src/basis_correction/weak_corr_func.irp.f Normal file
View File

@ -0,0 +1,83 @@
BEGIN_PROVIDER [double precision, ecmd_lda_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_lda_mu_of_r = multi-determinantal Ecmd within the LDA approximation with mu(r) ,
!
! see equation 40 in J. Chem. Phys. 149, 194301 (2018); https://doi.org/10.1063/1.5052714
END_DOC
implicit none
integer :: ipoint,istate
double precision :: rho_a, rho_b, ec
double precision :: wall0,wall1,weight,mu
logical :: dospin
dospin = .true. ! JT dospin have to be set to true for open shell
print*,'Providing ecmd_lda_mu_of_r ...'
ecmd_lda_mu_of_r = 0.d0
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
mu = mu_of_r_prov(ipoint,istate)
weight = final_weight_at_r_vector(ipoint)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
if(isnan(ec))then
print*,'ec is nan'
stop
endif
ecmd_lda_mu_of_r(istate) += weight * ec
enddo
enddo
call wall_time(wall1)
print*,'Time for ecmd_lda_mu_of_r :',wall1-wall0
END_PROVIDER
BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
BEGIN_DOC
! ecmd_pbe_ueg_mu_of_r = multi-determinantal Ecmd within the PBE-UEG approximation with mu(r) ,
!
! see Eqs. 13-14b in Phys.Chem.Lett.2019, 10, 2931 2937; https://pubs.acs.org/doi/10.1021/acs.jpclett.9b01176
!
! Based on the PBE-on-top functional (see Eqs. 26, 27 of J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.5082638)
!
! but it the on-top pair density of the UEG as an approximation of the exact on-top pair density
END_DOC
implicit none
double precision :: weight
integer :: ipoint,istate
double precision :: eps_c_md_PBE,mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
double precision :: g0_UEG_mu_inf
ecmd_pbe_ueg_mu_of_r = 0.d0
print*,'Providing ecmd_pbe_ueg_mu_of_r ...'
call wall_time(wall0)
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
weight=final_weight_at_r_vector(ipoint)
! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
mu = mu_of_r_prov(ipoint,istate)
rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
on_top = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
! The form of interpolated (mu=0 ---> mu=infinity) functional originally introduced in JCP, 150, 084103 1-10 (2019)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_PBE)
ecmd_pbe_ueg_mu_of_r(istate) += eps_c_md_PBE * weight
enddo
enddo
double precision :: wall1, wall0
call wall_time(wall1)
print*,'Time for the ecmd_pbe_ueg_mu_of_r:',wall1-wall0
END_PROVIDER

1
src/becke_numerical_grid/grid_becke_vector.irp.f
View File

@ -18,6 +18,7 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
enddo
print*,'n_points_final_grid = ',n_points_final_grid
print*,'n max point = ',n_points_integration_angular*(n_points_radial_grid*nucl_num - 1)
call ezfio_set_becke_numerical_grid_n_points_final_grid(n_points_final_grid)
END_PROVIDER
BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)]

2
src/cas_based_on_top/NEED Normal file
View File

@ -0,0 +1,2 @@
two_body_rdm
dft_utils_in_r

9
src/cas_based_on_top/README.rst Normal file
View File

@ -0,0 +1,9 @@
==============
on_top_density
==============
This plugin proposes different routines/providers to compute the on-top pair density of a CAS-based wave function.
This means that all determinants in psi_det must belong to an active-space.
As usual, see the file "example.irp.f" to get introduced to the main providers/routines.

125
src/cas_based_on_top/c_i_a_v_mos.irp.f Normal file
View File

@ -0,0 +1,125 @@
BEGIN_PROVIDER[double precision, core_mos_in_r_array, (n_core_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, core_mos_in_r_array_transp,(n_points_final_grid,n_core_orb)]
implicit none
BEGIN_DOC
! all COREE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_core_orb
j = list_core(i)
do k = 1, n_points_final_grid
core_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_core_orb
core_mos_in_r_array(i,k) = core_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, inact_mos_in_r_array, (n_inact_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, inact_mos_in_r_array_transp,(n_points_final_grid,n_inact_orb)]
implicit none
BEGIN_DOC
! all INACTIVE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_inact_orb
j = list_inact(i)
do k = 1, n_points_final_grid
inact_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_inact_orb
inact_mos_in_r_array(i,k) = inact_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, act_mos_in_r_array, (n_act_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, act_mos_in_r_array_transp,(n_points_final_grid,n_act_orb)]
implicit none
BEGIN_DOC
! all ACTIVE MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_act_orb
j = list_act(i)
do k = 1, n_points_final_grid
act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_act_orb
act_mos_in_r_array(i,k) = act_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, virt_mos_in_r_array, (n_virt_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, virt_mos_in_r_array_transp,(n_points_final_grid,n_virt_orb)]
implicit none
BEGIN_DOC
! all VIRTUAL MOs on the grid points, arranged in two different ways
END_DOC
integer :: i,j,k
do i = 1, n_virt_orb
j = list_virt(i)
do k = 1, n_points_final_grid
virt_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_virt_orb
virt_mos_in_r_array(i,k) = virt_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array, (n_core_inact_act_orb,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, core_inact_act_mos_in_r_array_transp,(n_points_final_grid,n_core_inact_act_orb)]
implicit none
integer :: i,j,k
do i = 1, n_core_inact_act_orb
j = list_core_inact_act(i)
do k = 1, n_points_final_grid
core_inact_act_mos_in_r_array_transp(k,i) = mos_in_r_array_transp(k,j)
enddo
enddo
do k = 1, n_points_final_grid
do i = 1, n_core_inact_act_orb
core_inact_act_mos_in_r_array(i,k) = core_inact_act_mos_in_r_array_transp(k,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, core_inact_act_mos_grad_in_r_array, (3,n_core_inact_act_orb,n_points_final_grid)]
implicit none
integer :: i,j,k,l
do i = 1, n_core_inact_act_orb
j = list_core_inact_act(i)
do k = 1, n_points_final_grid
do l = 1, 3
core_inact_act_mos_grad_in_r_array(l,i,k) = mos_grad_in_r_array(j,k,l)
enddo
enddo
enddo
END_PROVIDER

47
src/cas_based_on_top/cas_based_density.irp.f Normal file
View File

@ -0,0 +1,47 @@
program cas_based_density
implicit none
BEGIN_DOC
! TODO : Small example to use the different quantities in this plugin
END_DOC
!! You force QP2 to read the wave function in the EZFIO folder
!! It is assumed that all Slater determinants in the wave function
!! belongs to an active space defined by core, inactive and active list of orbitals
read_wf = .True.
touch read_wf
call routine_test_cas_based_density
end
subroutine routine_test_cas_based_density
implicit none
integer :: ipoint, istate
double precision :: accu_n_elec(N_states),accu_n_elec_2(N_states)
! PROVIDERS
accu_n_elec = 0.d0
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
accu_n_elec(istate) += one_e_cas_total_density(ipoint,istate) * final_weight_at_r_vector(ipoint)
enddo
print*,'istate = ',istate
print*,'accu_n_elec = ',accu_n_elec(istate)
enddo
! ROUTINES
double precision :: r(3),core_dens,inact_dens,act_dens(2,N_states),total_cas_dens(N_states)
accu_n_elec_2 = 0.d0
do ipoint = 1, n_points_final_grid
r(:) = final_grid_points(:,ipoint)
call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
do istate = 1, N_states
accu_n_elec_2(istate) += total_cas_dens(istate) * final_weight_at_r_vector(ipoint)
enddo
enddo
do istate = 1, N_states
print*,'istate = ',istate
print*,'accu_n_elec = ',accu_n_elec_2(istate)
enddo
end

19
src/cas_based_on_top/cas_based_on_top.irp.f Normal file
View File

@ -0,0 +1,19 @@
program cas_based_on_top_density
implicit none
BEGIN_DOC
! TODO : Small example to use the different quantities in this plugin
END_DOC
!! You force QP2 to read the wave function in the EZFIO folder
!! It is assumed that all Slater determinants in the wave function
!! belongs to an active space defined by core, inactive and active list of orbitals
read_wf = .True.
touch read_wf
! call routine_test_cas_based_on_top_density
call routine
end
subroutine routine
implicit none
provide total_cas_on_top_density
end

99
src/cas_based_on_top/cas_dens_prov.irp.f Normal file
View File

@ -0,0 +1,99 @@
BEGIN_PROVIDER [double precision, one_e_cas_total_density ,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_cas_total_density = TOTAL DENSITY FOR a CAS wave function
!
! WARNING : if "no_core_density" == .True. then the core part of density is ignored
END_DOC
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
one_e_cas_total_density(ipoint,istate) = one_e_act_density_alpha(ipoint,istate) + one_e_act_density_beta(ipoint,istate) &
+ 2.d0 * inact_density(ipoint)
if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY
one_e_cas_total_density(ipoint,istate) += 2.d0 * core_density(ipoint)
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_density_alpha,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_act_density_alpha = pure ACTIVE part of the DENSITY for ALPHA ELECTRONS
END_DOC
one_e_act_density_alpha = 0.d0
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
do i = 1, n_act_orb
do j = 1, n_act_orb
one_e_act_density_alpha(ipoint,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_density_beta,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! one_e_act_density_beta = pure ACTIVE part of the DENSITY for BETA ELECTRONS
END_DOC
one_e_act_density_beta = 0.d0
integer :: ipoint,i,j,istate
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
do i = 1, n_act_orb
do j = 1, n_act_orb
one_e_act_density_beta(ipoint,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos_in_r_array(j,ipoint) * act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, inact_density, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! INACTIVE part of the density for alpha/beta.
!
! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE INACTIVE,
!
! YOU MUST MULTIPLY BY TWO
END_DOC
integer :: i,j
inact_density = 0.d0
do i = 1, n_points_final_grid
do j = 1, n_inact_orb
inact_density(i) += inact_mos_in_r_array(j,i) **2
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, core_density, (n_points_final_grid) ]
implicit none
BEGIN_DOC
! CORE part of the density for alpha/beta.
!
! WARNING :: IF YOU NEED THE TOTAL DENSITY COMING FROM THE CORE,
!
! YOU MUST MULTIPLY BY TWO
END_DOC
integer :: i,j
core_density = 0.d0
do i = 1, n_points_final_grid
do j = 1, n_core_orb
core_density(i) += core_mos_in_r_array(j,i) **2
enddo
enddo
END_PROVIDER

65
src/cas_based_on_top/cas_dens_rout.irp.f Normal file
View File

@ -0,0 +1,65 @@
subroutine give_cas_density_in_r(core_dens,inact_dens,act_dens,total_cas_dens,r)
implicit none
BEGIN_DOC
! returns the different component of the density at grid point r(3) for a CAS wave function
!
! core_dens : density coming from the CORE orbitals
!
! inact_dens : density coming from the INACT orbitals
!
! act_dens(1/2,1:N_states) : active part of the alpha/beta electrons for all states
!
! total_cas_dens : total density of the cas wave function
!
! WARNING : if "no_core_density" == .True. then the core part of density is ignored in total_cas_dens
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: core_dens, inact_dens, act_dens(2,N_states), total_cas_dens(N_states)
double precision, allocatable :: mos_array(:),act_mos(:)
allocate(mos_array(mo_num))
call give_all_mos_at_r(r,mos_array)
integer :: i,iorb,j,jorb,istate
! core part of the density
core_dens = 0.d0
do i = 1, n_core_orb
iorb = list_core(i)
core_dens += mos_array(iorb)*mos_array(iorb)
enddo
core_dens = core_dens * 2.d0
! inactive part of the density
inact_dens = 0.d0
do i = 1, n_inact_orb
iorb = list_inact(i)
inact_dens += mos_array(iorb)*mos_array(iorb)
enddo
inact_dens = inact_dens * 2.d0
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
iorb = list_act(i)
act_mos(i) = mos_array(iorb)
enddo
! active part of the density for alpha/beta and all states
act_dens = 0.d0
do istate = 1, N_states
do i = 1, n_act_orb
do j = 1, n_act_orb
act_dens(1,istate) += one_e_act_dm_alpha_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
act_dens(2,istate) += one_e_act_dm_beta_mo_for_dft(j,i,istate) * act_mos(j) * act_mos(i)
enddo
enddo
enddo
! TOTAL density for all states
do istate = 1, N_states
total_cas_dens(istate) = inact_dens + act_dens(1,istate) + act_dens(2,istate)
if(.not.no_core_density)then !!! YOU ADD THE CORE DENSITY
total_cas_dens(istate) += core_dens
endif
enddo
end

37
src/cas_based_on_top/cas_one_e_rdm.irp.f Normal file
View File

@ -0,0 +1,37 @@
BEGIN_PROVIDER [double precision, one_e_act_dm_beta_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
implicit none
BEGIN_DOC
! one_e_act_dm_beta_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the BETA ELECTRONS
END_DOC
integer :: i,j,ii,jj,istate
do istate = 1, N_states
do ii = 1, n_act_orb
i = list_act(ii)
do jj = 1, n_act_orb
j = list_act(jj)
one_e_act_dm_beta_mo_for_dft(jj,ii,istate) = one_e_dm_mo_beta_for_dft(j,i,istate)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, one_e_act_dm_alpha_mo_for_dft, (n_act_orb,n_act_orb,N_states)]
implicit none
BEGIN_DOC
! one_e_act_dm_alpha_mo_for_dft = pure ACTIVE part of the ONE ELECTRON REDUCED DENSITY MATRIX for the ALPHA ELECTRONS
END_DOC
integer :: i,j,ii,jj,istate
do istate = 1, N_states
do ii = 1, n_act_orb
i = list_act(ii)
do jj = 1, n_act_orb
j = list_act(jj)
one_e_act_dm_alpha_mo_for_dft(jj,ii,istate) = one_e_dm_mo_alpha_for_dft(j,i,istate)
enddo
enddo
enddo
END_PROVIDER

68
src/cas_based_on_top/eff_spin_dens.irp.f Normal file
View File

@ -0,0 +1,68 @@
BEGIN_PROVIDER [double precision, effective_spin_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_spin_dm, (3,n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! effective_spin_dm(r_i) = \sqrt( n(r)^2 - 4 * ontop(r) )
! effective spin density obtained from the total density and on-top pair density
! see equation (6) of Phys. Chem. Chem. Phys., 2015, 17, 22412--22422 | 22413
END_DOC
provide total_cas_on_top_density
integer :: i_point,i_state,i
double precision :: n2,m2,thr
thr = 1.d-14
effective_spin_dm = 0.d0
grad_effective_spin_dm = 0.d0
do i_state = 1, N_states
do i_point = 1, n_points_final_grid
n2 = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state) + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
! density squared
n2 = n2 * n2
if(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state).gt.thr)then
effective_spin_dm(i_point,i_state) = dsqrt(n2 - 4.D0 * total_cas_on_top_density(i_point,i_state))
if(isnan(effective_spin_dm(i_point,i_state)))then
print*,'isnan(effective_spin_dm(i_point,i_state)'
stop
endif
m2 = effective_spin_dm(i_point,i_state)
m2 = 0.5d0 / m2 ! 1/(2 * sqrt(n(r)^2 - 4 * ontop(r)) )
do i = 1, 3
grad_effective_spin_dm(i,i_point,i_state) = m2 * ( one_e_stuff_for_pbe(i,i_point,i_state) - 4.d0 * grad_total_cas_on_top_density(i,i_point,i_state) )
enddo
else
effective_spin_dm(i_point,i_state) = 0.d0
grad_effective_spin_dm(:,i_point,i_state) = 0.d0
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, effective_alpha_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, effective_beta_dm, (n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_alpha_dm, (3,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, grad_effective_beta_dm, (3,n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! effective_alpha_dm(r_i) = 1/2 * (effective_spin_dm(r_i) + n(r_i))
! effective_beta_dm(r_i) = 1/2 * (-effective_spin_dm(r_i) + n(r_i))
END_DOC
provide total_cas_on_top_density
integer :: i_point,i_state,i
double precision :: n,grad_n
do i_state = 1, N_states
do i_point = 1, n_points_final_grid
n = (one_e_dm_and_grad_alpha_in_r(4,i_point,i_state) + one_e_dm_and_grad_beta_in_r(4,i_point,i_state))
effective_alpha_dm(i_point,i_state) = 0.5d0 * (n + effective_spin_dm(i_point,i_state))
effective_beta_dm(i_point,i_state) = 0.5d0 * (n - effective_spin_dm(i_point,i_state))
do i = 1, 3
grad_n = (one_e_dm_and_grad_alpha_in_r(i,i_point,i_state) + one_e_dm_and_grad_beta_in_r(i,i_point,i_state))
grad_effective_alpha_dm(i,i_point,i_state) = 0.5d0 * (grad_n + grad_effective_spin_dm(i,i_point,i_state) )
grad_effective_beta_dm(i,i_point,i_state) = 0.5d0 * (grad_n - grad_effective_spin_dm(i,i_point,i_state) )
enddo
enddo
enddo
END_PROVIDER

41
src/cas_based_on_top/example.irp.f Normal file
View File

@ -0,0 +1,41 @@
subroutine write_on_top_in_real_space
implicit none
BEGIN_DOC
! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid
END_DOC
double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist
integer :: nz,i,istate
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'.on_top'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
zmax = 2.0d0
print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
print*,'nucl_coord(2,3) = ',nucl_coord(2,3)
dist = dabs(nucl_coord(1,3) - nucl_coord(2,3))
zmax += dist
zcenter = (nucl_coord(1,3) + nucl_coord(2,3))*0.5d0
print*,'zcenter = ',zcenter
print*,'zmax = ',zmax
nz = 1000
dz = zmax / dble(nz)
r(:) = 0.d0
r(3) = zcenter -zmax * 0.5d0
print*,'r(3) = ',r(3)
istate = 1
do i = 1, nz
call give_on_top_in_r_one_state(r,istate,on_top_in_r)
call give_cas_density_in_r(r,total_density)
write(i_unit_output,*)r(3),on_top_in_r,total_density
r(3) += dz
enddo
end

76
src/cas_based_on_top/on_top_cas_prov.irp.f Normal file
View File

@ -0,0 +1,76 @@
subroutine act_on_top_on_grid_pt(ipoint,istate,pure_act_on_top_of_r)
implicit none
BEGIN_DOC
! act_on_top_on_grid_pt returns the purely ACTIVE part of the on top pair density
!
! at the grid point ipoint, for the state istate
END_DOC
integer, intent(in) :: ipoint,istate
double precision, intent(out) :: pure_act_on_top_of_r
double precision :: phi_i,phi_j,phi_k,phi_l
integer :: i,j,k,l
ASSERT (istate <= N_states)
pure_act_on_top_of_r = 0.d0
do l = 1, n_act_orb
phi_l = act_mos_in_r_array(l,ipoint)
do k = 1, n_act_orb
phi_k = act_mos_in_r_array(k,ipoint)
do j = 1, n_act_orb
phi_j = act_mos_in_r_array(j,ipoint)
do i = 1, n_act_orb
phi_i = act_mos_in_r_array(i,ipoint)
! 1 2 1 2
pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l
enddo
enddo
enddo
enddo
end
BEGIN_PROVIDER [double precision, total_cas_on_top_density,(n_points_final_grid,N_states) ]
implicit none
BEGIN_DOC
! on top pair density :: n2(r,r) at each of the Becke's grid point of a CAS-BASED wf
!
! Contains all core/inact/act contribution.
!
! !!!!! WARNING !!!!! If no_core_density then you REMOVE ALL CONTRIBUTIONS COMING FROM THE CORE ORBITALS
END_DOC
integer :: i_point,istate
double precision :: wall_0,wall_1,core_inact_dm,pure_act_on_top_of_r
logical :: no_core
print*,'providing the total_cas_on_top_density'
! for parallelization
provide inact_density core_density one_e_act_density_beta one_e_act_density_alpha act_mos_in_r_array
i_point = 1
istate = 1
call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
call wall_time(wall_0)
if(no_core_density)then
print*,'USING THE VALENCE ONLY TWO BODY DENSITY'
endif
do istate = 1, N_states
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_point,core_inact_dm,pure_act_on_top_of_r) &
!$OMP SHARED(total_cas_on_top_density,n_points_final_grid,inact_density,core_density,one_e_act_density_beta,one_e_act_density_alpha,no_core_density,istate)
do i_point = 1, n_points_final_grid
call act_on_top_on_grid_pt(i_point,istate,pure_act_on_top_of_r)
if(no_core_density) then
core_inact_dm = inact_density(i_point)
else
core_inact_dm = (inact_density(i_point) + core_density(i_point))
endif
total_cas_on_top_density(i_point,istate) = pure_act_on_top_of_r + core_inact_dm * (one_e_act_density_beta(i_point,istate) + one_e_act_density_alpha(i_point,istate)) + core_inact_dm*core_inact_dm
enddo
!$OMP END PARALLEL DO
enddo
call wall_time(wall_1)
print*,'provided the total_cas_on_top_density'
print*,'Time to provide :',wall_1 - wall_0
END_PROVIDER

118
src/cas_based_on_top/on_top_cas_rout.irp.f Normal file
View File

@ -0,0 +1,118 @@
subroutine give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
implicit none
double precision, intent(in) :: r(3),mos_array(mo_num)
double precision, intent(out):: core_density_in_r,inact_density_in_r,act_density_in_r(2,N_states),total_density(N_states)
BEGIN_DOC
! core, inactive and active part of the density for alpha/beta electrons
!
! the density coming from the core and inactive are the same for alpha/beta electrons
!
! act_density(1/2, i) = alpha/beta density for the ith state
!
! total_density(i) = 2 * (core_density_in_r+inact_density_in_r) + act_density_in_r(1,i) + act_density_in_r(2,i)
END_DOC
integer :: i,j,istate
core_density_in_r = 0.d0
do i = 1, n_core_orb
j = list_core(i)
core_density_in_r += mos_array(j) * mos_array(j)
enddo
inact_density_in_r = 0.d0
do i = 1, n_inact_orb
j = list_inact(i)
inact_density_in_r += mos_array(j) * mos_array(j)
enddo
double precision, allocatable :: act_mos(:)
double precision :: tmp
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
j = list_act(i)
act_mos(i) = mos_array(j)
enddo
act_density_in_r = 0.d0
do istate = 1, N_states
do i = 1, n_act_orb
do j = 1, n_act_orb
tmp = act_mos(i) * act_mos(j)
act_density_in_r(1,istate) += tmp * one_e_act_dm_alpha_mo_for_dft(j,i,istate)
act_density_in_r(2,istate) += tmp * one_e_act_dm_beta_mo_for_dft(j,i,istate)
enddo
enddo
total_density(istate) = 2.d0 * (core_density_in_r + inact_density_in_r) + act_density_in_r(1,istate) + act_density_in_r(2,istate)
enddo
end
subroutine give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
implicit none
BEGIN_DOC
! gives the purely active on-top pair density for a given state
END_DOC
integer, intent(in) :: istate
double precision, intent(in) :: r(3),mos_array(mo_num)
double precision, intent(out) :: act_on_top
double precision :: phi_i,phi_j,phi_k,phi_l
integer :: i,j,k,l
double precision, allocatable :: act_mos(:)
double precision :: tmp
allocate(act_mos(n_act_orb))
do i = 1, n_act_orb
j = list_act(i)
act_mos(i) = mos_array(j)
enddo
act_on_top = 0.d0
do l = 1, n_act_orb
phi_l = act_mos(l)
do k = 1, n_act_orb
phi_k = act_mos(k)
do j = 1, n_act_orb
phi_j = act_mos(j)
tmp = phi_l * phi_k * phi_j
do i = 1, n_act_orb
phi_i = act_mos(i)
! 1 2 1 2
act_on_top += act_2_rdm_ab_mo(i,j,k,l,istate) * tmp * phi_i
enddo
enddo
enddo
enddo
end
subroutine give_on_top_in_r_one_state(r,istate,on_top_in_r)
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: r(3)
double precision, intent(out) :: on_top_in_r
BEGIN_DOC
! on top pair density in r for the state istate a CAS-BASED wf
!
! note that if no_core_density .EQ. .True., all core contributions are excluded
END_DOC
double precision, allocatable :: mos_array(:)
provide act_2_rdm_ab_mo one_e_act_dm_alpha_mo_for_dft one_e_act_dm_beta_mo_for_dft
allocate(mos_array(mo_num))
call give_all_mos_at_r(r,mos_array)
double precision :: core_density_in_r, inact_density_in_r, act_density_in_r(2,N_states), total_density(N_states)
double precision :: act_on_top,core_inact_dm
! getting the different part of the density in r
call give_core_inact_act_density_in_r(r,mos_array,core_density_in_r,inact_density_in_r,act_density_in_r, total_density)
! getting the purely active part of the density in r
call give_active_on_top_in_r_one_state(r,istate,mos_array,act_on_top)
if(no_core_density) then
core_inact_dm = inact_density_in_r
else
core_inact_dm = core_density_in_r + inact_density_in_r
endif
on_top_in_r = act_on_top + core_inact_dm * (act_density_in_r(1,istate) + act_density_in_r(2,istate)) + core_inact_dm*core_inact_dm
end

73
src/cas_based_on_top/on_top_grad.irp.f Normal file
View File

@ -0,0 +1,73 @@
subroutine give_on_top_gradient(ipoint,istate,ontop_grad)
implicit none
BEGIN_DOC
! on top pair density and its gradient evaluated at a given point of the grid
! ontop_grad(1:3) :: gradients of the on-top pair density
! ontop_grad(4) :: on-top pair density
END_DOC
double precision, intent(out) :: ontop_grad(4)
integer, intent(in) :: ipoint,istate
double precision :: phi_jkl,phi_ikl,phi_ijl,phi_ijk
integer :: i,j,k,l,m
ontop_grad = 0.d0
do l = 1, n_core_inact_act_orb
do k = 1, n_core_inact_act_orb
do j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb
phi_jkl = core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ikl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(k,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ijl = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(l,ipoint)
phi_ijk = core_inact_act_mos_in_r_array(i,ipoint) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(k,ipoint)
! 1 2 1 2
ontop_grad(4) += phi_ijk * core_inact_act_mos_in_r_array(l,ipoint) * full_occ_2_rdm_ab_mo(i,j,k,l,istate)
do m = 1,3
ontop_grad (m) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * &
( core_inact_act_mos_grad_in_r_array(m,i,ipoint) * phi_jkl + core_inact_act_mos_grad_in_r_array(m,j,ipoint) * phi_ikl + &
core_inact_act_mos_grad_in_r_array(m,k,ipoint) * phi_ijl + core_inact_act_mos_grad_in_r_array(m,l,ipoint) * phi_ijk )
enddo
enddo
enddo
enddo
enddo
end
BEGIN_PROVIDER [double precision, grad_total_cas_on_top_density,(4,n_points_final_grid,N_states) ]
&BEGIN_PROVIDER [double precision, wall_time_core_inact_act_on_top_of_r ]
implicit none
BEGIN_DOC
! grad_total_cas_on_top_density(1:3,ipoint,istate) : provider for the on top pair density gradient (x,y,z) for the point 'ipoint' and state 'istate'
!
! grad_total_cas_on_top_density(4,ipoint,istate) : on top pair density for the point 'ipoint' and state 'istate'
END_DOC
integer :: i_point,i_state,i
double precision :: wall_0,wall_1
double precision :: core_inact_act_on_top_of_r_from_provider,ontop_grad(4)
print*,'providing the core_inact_act_on_top_of_r'
i_point = 1
provide full_occ_2_rdm_ab_mo
i_state = 1
call give_on_top_gradient(i_point,i_state,ontop_grad)
call wall_time(wall_0)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_point,i_state,ontop_grad) &
!$OMP SHARED(grad_total_cas_on_top_density,n_points_final_grid,N_states)
do i_point = 1, n_points_final_grid
do i_state = 1, N_states
call give_on_top_gradient(i_point,i_state,ontop_grad)
do i = 1, 4
grad_total_cas_on_top_density(i,i_point,i_state) = ontop_grad(i)
enddo
enddo
enddo
!$OMP END PARALLEL DO
call wall_time(wall_1)
print*,'provided the core_inact_act_on_top_of_r'
print*,'Time to provide :',wall_1 - wall_0
wall_time_core_inact_act_on_top_of_r = wall_1 - wall_0
END_PROVIDER

8
src/dft_one_e/NEED
View File

@ -1,2 +1,8 @@
dft_utils_one_e
dft_utils_func
functionals
mo_one_e_ints
mo_two_e_ints
ao_one_e_ints
ao_two_e_ints
mo_two_e_erf_ints
ao_two_e_erf_ints

0
src/dft_utils_one_e/effective_pot.irp.f → src/dft_one_e/effective_pot.irp.f
View File

0
src/dft_utils_one_e/mu_erf_dft.irp.f → src/dft_one_e/mu_erf_dft.irp.f
View File

0
src/dft_utils_one_e/one_e_energy_dft.irp.f → src/dft_one_e/one_e_energy_dft.irp.f
View File

0
src/dft_utils_one_e/sr_coulomb.irp.f → src/dft_one_e/sr_coulomb.irp.f
View File

2
src/dft_utils_func/NEED Normal file
View File

@ -0,0 +1,2 @@
density_for_dft
dft_utils_in_r

0
src/dft_utils_one_e/README.rst → src/dft_utils_func/README.rst
View File

152
src/dft_utils_func/ecmd_lda.irp.f Normal file
View File

@ -0,0 +1,152 @@
!****************************************************************************
subroutine ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,e)
!*****************************************************************************
! Short-range spin-dependent LDA correlation functional with multideterminant reference
! for OEP calculations from Section V of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Input: rhot : total density
! rhos : spin density
! mu : Interation parameter
! dospin : use spin density
!
! Ouput: e : energy
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision, intent(in) :: rho_a,rho_b,mu
logical, intent(in) :: dospin
double precision, intent(out):: e
double precision :: e1
double precision :: rhoa,rhob
double precision :: rhot, rhos
rhoa=max(rho_a,1.0d-15)
rhob=max(rho_b,1.0d-15)
rhot = rhoa + rhob
rhos = rhoa - rhob
call ec_only_lda_sr(mu,rho_a,rho_b,e1)
if(isnan(e1))then
print*,'e1 is NaN'
print*,mu,rho_a,rho_b
stop
endif
call DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e)
if(isnan(e))then
print*,'e is NaN'
print*,mu,rhot,rhos
stop
endif
e = e1 + e
end
!****************************************************************************
subroutine DELTA_LRSR_LDAERF (rhot,rhos,mu,dospin,e)
!*****************************************************************************
! LDA approximation to term Delta_(LR-SR) from Eq. 42 of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Input: rhot : total density
! rhos : spin density
! mu : Interation parameter
! dospin : use spin density
!
! Ouput: e : energy
!
! Warning: not tested for z != 0
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision rhot, rhos, mu
logical dospin
double precision e
double precision f13, f83, pi, rsfac, alpha2
double precision rs, rs2, rs3
double precision rhoa, rhob, z, z2, onepz, onemz, zp, zm, phi8
double precision g0f, g0s
double precision bd2, bd3
double precision c45, c4, c5
double precision bc2, bc4, bc3t, bc5t, d0
double precision delta2,delta3,delta4,delta5,delta6
double precision delta
parameter(f13 = 0.333333333333333d0)
parameter(f83 = 2.6666666666666665d0)
parameter(pi = 3.141592653589793d0)
parameter(rsfac = 0.620350490899400d0)
parameter(alpha2 = 0.2715053589826032d0)
rs = rsfac/(rhot**f13)
rs2 = rs*rs
rs3 = rs2*rs
! Spin-unpolarized case
if (.not.dospin) then
z = 0.d0
! Spin-polarized case
else
rhoa=max((rhot+rhos)*.5d0,1.0d-15)
rhob=max((rhot-rhos)*.5d0,1.0d-15)
z=min((rhoa-rhob)/(rhoa+rhob),0.9999999999d0)
endif
z2=z*z
bd2=dexp(-0.547d0*rs)*(-0.388d0*rs+0.676*rs2)/rs2
bd3=dexp(-0.31d0*rs)*(-4.95d0*rs+rs2)/rs3
onepz=1.d0+z
onemz=1.d0-z
phi8=0.5d0*(onepz**f83+onemz**f83)
zp=onepz/2.d0
zm=onemz/2.d0
c45=(zp**2)*g0s(rs*zp**(-f13))+(zm**2)*g0s(rs*zm**(-f13))
c4=c45+(1.d0-z2)*bd2-phi8/(5.d0*alpha2*rs2)
c5=c45+(1.d0-z2)*bd3
bc2=-3.d0*(1-z2)*(g0f(rs)-0.5d0)/(8.d0*rs3)
bc4=-9.d0*c4/(64.d0*rs3)
bc3t=-(1-z2)*g0f(rs)*(2.d0*dsqrt(2.d0)-1)/(2.d0*dsqrt(pi)*rs3)
bc5t = -3.d0*c5*(3.d0-dsqrt(2.d0))/(20.d0*dsqrt(2.d0*pi)*rs3)
d0=(0.70605d0+0.12927d0*z2)*rs
delta2=0.073867d0*(rs**(1.5d0))
delta3=4*(d0**6)*bc3t+(d0**8)*bc5t;
delta4=4*(d0**6)*bc2+(d0**8)*bc4;
delta5=(d0**8)*bc3t;
delta6=(d0**8)*bc2;
delta=(delta2*(mu**2)+delta3*(mu**3)+delta4*(mu**4)+delta5*(mu**5)+delta6*(mu**6))/((1+(d0**2)*(mu**2))**4)
! multiply by rhot to get energy density
e=delta*rhot
end
!*****************************************************************************
double precision function g0s(rs)
!*****************************************************************************
! g"(0,rs,z=1) from Eq. 32 of
! Paziani, Moroni, Gori-Giorgi and Bachelet, PRB 73, 155111 (2006)
!
! Created: 26-08-11, J. Toulouse
!*****************************************************************************
implicit none
double precision rs
double precision rs2, f53, alpha2
parameter(f53 = 1.6666666666666667d0)
parameter(alpha2 = 0.2715053589826032d0)
rs2=rs*rs
g0s=(2.d0**f53)*(1.d0-0.02267d0*rs)/((5.d0*alpha2*rs2)*(1.d0+0.4319d0*rs+0.04d0*rs2))
end

53
src/dft_utils_func/ecmd_pbe_general.irp.f Normal file
View File

@ -0,0 +1,53 @@
subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,eps_c_md_on_top_PBE)
implicit none
BEGIN_DOC
!
! General e_cmd functional interpolating between :
!
! +) the large mu behaviour in cst/(\mu^3) on-top
!
! +) mu= 0 with the usal ec_pbe at
!
! Depends on : mu, the density (rho_a,rho_b), the square of the density gradient (grad_rho_a,grad_rho_b)
!
! the flavour of on-top densiyt (on_top) you fill in: in principle it should be the exact on-top
!
! The form of the functional was originally introduced in JCP, 150, 084103 1-10 (2019)
!
END_DOC
double precision, intent(in) :: mu,rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),on_top
double precision, intent(out) :: eps_c_md_on_top_PBE
double precision :: pi, e_pbe,beta,denom
double precision :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
integer :: m
pi = 4.d0 * datan(1.d0)
eps_c_md_on_top_PBE = 0.d0
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do m = 1, 3
grad_rho_a_2 += grad_rho_a(m)*grad_rho_a(m)
grad_rho_b_2 += grad_rho_b(m)*grad_rho_b(m)
grad_rho_a_b += grad_rho_a(m)*grad_rho_b(m)
enddo
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi)* on_top
if (dabs(denom) > 1.d-12) then
! quantity of Eq. (26)
beta = (3.d0*e_PBE)/denom
eps_c_md_on_top_PBE = e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_on_top_PBE =0.d0
endif
end

72
src/dft_utils_func/ecmd_pbe_on_top.irp.f Normal file
View File

@ -0,0 +1,72 @@
subroutine ec_md_on_top_PBE_mu_corrected(mu,r,two_dm,eps_c_md_on_top_PBE)
implicit none
BEGIN_DOC
! enter with "r(3)", and "two_dm(N_states)" which is the on-top pair density at "r" for each states
!
! you get out with the energy density defined in J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.508263
!
! by Eq. (26), which includes the correction of the on-top pair density of Eq. (29).
END_DOC
double precision, intent(in) :: mu , r(3), two_dm
double precision, intent(out) :: eps_c_md_on_top_PBE(N_states)
double precision :: two_dm_in_r, pi, e_pbe(N_states),beta(N_states),mu_correction_of_on_top
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo,on_top_corrected
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_on_top_PBE = 0.d0
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE(istate))
! correction of the on-top pair density according to Eq. (29)
on_top_corrected = mu_correction_of_on_top(mu,two_dm)
! quantity of Eq. (27) with a factor two according to the difference of normalization
! between the on-top of the JCP paper and that of QP2
beta(istate) = (3.d0*e_PBE(istate))/( (-2.d0+sqrt(2d0))*sqrt(2.d0*pi)*2.d0* on_top_corrected)
! quantity of Eq. (26)
eps_c_md_on_top_PBE(istate)=e_PBE(istate)/(1.d0+beta(istate)*mu**3)
enddo
end
double precision function mu_correction_of_on_top(mu,on_top)
implicit none
BEGIN_DOC
! mu-based correction to the on-top pair density provided by the assymptotic expansion of
!
! P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
!
! This is used in J. Chem. Phys.150, 084103 (2019); Eq. (29).
END_DOC
double precision, intent(in) :: mu,on_top
double precision :: pi
pi = 4.d0 * datan(1.d0)
mu_correction_of_on_top = on_top / ( 1.d0 + 2.d0/(dsqrt(pi)*mu) )
mu_correction_of_on_top = max(mu_correction_of_on_top ,1.d-15)
end

194
src/dft_utils_func/ecmd_pbe_ueg.irp.f Normal file
View File

@ -0,0 +1,194 @@
subroutine ecmd_pbe_ueg_at_r(mu,r,eps_c_md_PBE)
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with r(3), you get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu , r(3)
double precision, intent(out) :: eps_c_md_PBE(N_states)
double precision :: pi, e_PBE, beta
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
double precision :: g0_UEG_mu_inf, denom
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(mu == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a(istate)*rho_b(istate)*g0_UEG_mu_inf(rho_a(istate),rho_b(istate))
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end
subroutine eps_c_md_PBE_from_density(mu,rho_a,rho_b, grad_rho_a, grad_rho_b,eps_c_md_PBE) ! EG
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with the alpha/beta density and density gradients
!
! You get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu(N_states) , rho_a(N_states),rho_b(N_states), grad_rho_a(3,N_states),grad_rho_b(3,N_states)
double precision, intent(out) :: eps_c_md_PBE(N_states)
double precision :: pi, e_PBE, beta
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
double precision :: g0_UEG_mu_inf, denom
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
grad_rho_a_2 = 0.d0
grad_rho_b_2 = 0.d0
grad_rho_a_b = 0.d0
do istate = 1, N_states
do m = 1, 3
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
enddo
enddo
do istate = 1, N_states
! convertion from (alpha,beta) formalism to (closed, open) formalism
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(mu(istate) == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a(istate)*rho_b(istate)*g0_UEG_mu_inf(rho_a(istate),rho_b(istate))
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu(istate)**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end
subroutine eps_c_md_PBE_at_grid_pt(mu,i_point,eps_c_md_PBE)
implicit none
BEGIN_DOC
! provides the integrand of Eq. (13) of Phys.Chem.Lett.2019, 10, 2931 2937
!
! !!! WARNING !!! This is the total integrand of Eq. (13), which is e_cmd * n
!
! such a function is based on the exact behaviour of the Ecmd at large mu
!
! but with the exact on-top estimated with that of the UEG
!
! You enter with the alpha/beta density and density gradients
!
! You get out with eps_c_md_PBE(1:N_states)
END_DOC
double precision, intent(in) :: mu
double precision, intent(out) :: eps_c_md_PBE(N_states)
integer, intent(in) :: i_point
double precision :: two_dm, pi, e_pbe,beta,mu_correction_of_on_top
double precision :: grad_rho_a(3),grad_rho_b(3)
double precision :: grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
double precision :: rhoc,rhoo,ec_pbe_88
double precision :: delta,two_dm_corr,rho_a,rho_b
double precision :: grad_rho_2,denom,g0_UEG_mu_inf
double precision :: sigmacc,sigmaco,sigmaoo
integer :: m, istate
pi = 4.d0 * datan(1.d0)
eps_c_md_PBE = 0.d0
do istate = 1, N_states
! total and spin density
rhoc = one_e_dm_and_grad_alpha_in_r(4,i_point,istate) + one_e_dm_and_grad_beta_in_r(4,i_point,istate)
rhoo = one_e_dm_and_grad_alpha_in_r(4,i_point,istate) - one_e_dm_and_grad_beta_in_r(4,i_point,istate)
! gradients of the effective spin density
grad_rho_a_2 = 0.D0
grad_rho_b_2 = 0.D0
grad_rho_a_b = 0.D0
do m = 1, 3
grad_rho_a_2 += one_e_dm_and_grad_alpha_in_r(m,i_point,istate)**2.d0
grad_rho_b_2 += one_e_dm_and_grad_beta_in_r(m,i_point,istate) **2.d0
grad_rho_a_b += one_e_dm_and_grad_alpha_in_r(m,i_point,istate) * one_e_dm_and_grad_beta_in_r(m,i_point,istate)
enddo
sigmacc = grad_rho_a_2 + grad_rho_b_2 + 2.d0 * grad_rho_a_b
sigmaco = 0.d0
sigmaoo = 0.d0
rho_a = one_e_dm_and_grad_alpha_in_r(4,i_point,istate)
rho_b = one_e_dm_and_grad_beta_in_r(4,i_point,istate)
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE)
if(e_PBE.gt.0.d0)then
print*,'PBE gt 0 with regular dens'
endif
if(mu == 0.d0) then
eps_c_md_PBE(istate)=e_PBE
else
! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
denom = (-2.d0+dsqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a*rho_b*g0_UEG_mu_inf(rho_a,rho_b)
if (dabs(denom) > 1.d-12) then
beta = (3.d0*e_PBE)/denom
! Ecmd functional with the UEG ontop pair density when mu -> infty
! and the usual PBE correlation energy when mu = 0
eps_c_md_PBE(istate)=e_PBE/(1.d0+beta*mu**3)
else
eps_c_md_PBE(istate)=0.d0
endif
endif
enddo
end

0
src/dft_utils_one_e/garbage_func.irp.f → src/dft_utils_func/garbage_func.irp.f
View File

121
src/dft_utils_func/on_top_from_ueg.irp.f Normal file
View File

@ -0,0 +1,121 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
double precision function correction_to_on_top_from_UEG(mu,r,istate)
implicit none
integer, intent(in) :: istate
double precision, intent(in) :: mu,r(3)
double precision :: rho_a(N_states),rho_b(N_states)
double precision :: g0_UEG_mu_inf, g0_UEG_mu
call dm_dft_alpha_beta_at_r(r,rho_a,rho_b)
correction_to_on_top_from_UEG = g0_UEG_mu_inf(rho_a(istate),rho_b(istate)) / g0_UEG_mu(mu,rho_a(istate),rho_b(istate))
end
double precision function g0_UEG_mu_inf(rho_a,rho_b)
BEGIN_DOC
! Pair distribution function g0(n_alpha,n_beta) of the Colombic UEG
!
! Taken from Eq. (46) P. Gori-Giorgi and A. Savin, Phys. Rev. A 73, 032506 (2006).
END_DOC
implicit none
double precision, intent(in) :: rho_a,rho_b
double precision :: rho,pi,x
double precision :: B, C, D, E, d2, rs, ahd
rho = rho_a+rho_b
pi = 4d0 * datan(1d0)
ahd = -0.36583d0
d2 = 0.7524d0
B = -2d0 * ahd - d2
C = 0.08193d0
D = -0.01277d0
E = 0.001859d0
if (dabs(rho) > 1.d-12) then
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs
g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*exp(x)
else
g0_UEG_mu_inf= 0.d0
endif
end
double precision function g0_UEG_mu(mu,rho_a,rho_b)
implicit none
BEGIN_DOC
! Pair distribution function g0(n_alpha,n_beta) of the UEG interacting with the long range interaction erf(mu r12)/r12
!
! Taken from P. Gori-Giorgi and A. Savin, Phys. Rev. A 73, 032506 (2006).
END_DOC
double precision, intent(in) :: rho_a,rho_b,mu
double precision :: zeta,pi,rho,x,alpha
double precision :: B, C, D, E, d2, rs, ahd, h_func, kf
pi = 4d0 * datan(1d0)
rho = rho_a+rho_b
alpha = (4d0/(9d0*pi))**(1d0/3d0)
ahd = -0.36583d0
d2 = 0.7524d0
B = -2d0 * ahd - d2
C = 0.08193d0
D = -0.01277d0
E = 0.001859d0
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
kf = (alpha*rs)**(-1d0)
zeta = mu / kf
x = -d2*rs*h_func(zeta)/ahd
g0_UEG_mu = (exp(x)/2d0) * (1d0- B*(h_func(zeta)/ahd)*rs + C*((h_func(zeta)**2d0)/(ahd**2d0))*(rs**2d0) + D*((h_func(zeta)**3d0)/(ahd**3d0))*(rs**3d0) + E*((h_func(zeta)**4d0)/(ahd**4d0))*(rs**4d0) )
end
double precision function h_func(zeta)
implicit none
double precision, intent(in) :: zeta
double precision :: pi
double precision :: a1, a2, b1, b2, b3, ahd, alpha
pi = 4d0 * datan(1d0)
ahd = -0.36583d0
alpha = (4d0/(9d0*pi))**(1d0/3d0)
a1 = -(6d0*alpha/pi)*(1d0-log(2d0))
b1 = 1.4919d0
b3 = 1.91528d0
a2 = ahd * b3
b2 = (a1 - (b3*alpha/sqrt(pi)))/ahd
h_func = (a1*zeta**2d0 + a2*zeta**3d0) / (1d0 + b1*zeta + b2*zeta**2d0 + b3*zeta**3d0)
end
!-------------------------------------------------------------------------------------------------------------------------------------------
subroutine g0_dg0(rho, rho_a, rho_b, g0, dg0drho)
implicit none
BEGIN_DOC
! Give the on-top pair distribution function g0 and its derivative according to rho dg0drho
END_DOC
double precision, intent (in) :: rho, rho_a, rho_b
double precision, intent (out) :: g0, dg0drho
double precision :: pi
double precision :: g0_UEG_mu_inf, dg0drs
double precision :: C1, F1, D1, E1, B1, rs
pi = dacos(-1.d0)
C1 = 0.0819306d0
F1 = 0.752411d0
D1 = -0.0127713d0
E1 = 0.00185898d0
B1 = 0.7317d0 - F1
rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0)
g0 = g0_UEG_mu_inf(rho_a, rho_b)
dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*exp(-F1*rs)
dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
end subroutine g0_dg0

0
src/dft_utils_one_e/rho_ab_to_rho_tot.irp.f → src/dft_utils_func/rho_ab_to_rho_tot.irp.f
View File

0
src/dft_utils_one_e/routines_exc_sr_lda.irp.f → src/dft_utils_func/routines_exc_sr_lda.irp.f
View File

0
src/dft_utils_one_e/routines_exc_sr_pbe.irp.f → src/dft_utils_func/routines_exc_sr_pbe.irp.f
View File

0
src/dft_utils_one_e/utils.irp.f → src/dft_utils_func/utils.irp.f
View File

96
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -1,13 +1,15 @@
BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
implicit none
BEGIN_DOC
! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
!
! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,j) &
!$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
@ -15,11 +17,30 @@
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
aos_in_r_array(j,i) = aos_array(j)
aos_in_r_array_transp(i,j) = aos_array(j)
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
implicit none
BEGIN_DOC
! aos_in_r_array_transp(i,j) = value of the jth ao on the ith grid point
END_DOC
integer :: i,j
double precision :: aos_array(ao_num), r(3)
do i = 1, n_points_final_grid
do j = 1, ao_num
aos_in_r_array_transp(i,j) = aos_in_r_array(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
implicit none
BEGIN_DOC
@ -30,6 +51,10 @@
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,aos_grad_array,j,m) &
!$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
@ -41,15 +66,16 @@
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp, (n_points_final_grid,ao_num,3)]
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp, (3,ao_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! aos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth ao on the ith grid point
! aos_grad_in_r_array_transp(k,i,j) = value of the kth component of the gradient of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
@ -57,49 +83,18 @@
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
do m = 1, 3
do j = 1, ao_num
aos_grad_in_r_array_transp(i,j,m) = aos_grad_array(m,j)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_xyz, (3,ao_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! aos_grad_in_r_array_transp_xyz(k,i,j) = value of the kth component of the gradient of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
do m = 1, 3
do j = 1, ao_num
aos_grad_in_r_array_transp_xyz(m,j,i) = aos_grad_array(m,j)
aos_grad_in_r_array_transp(m,j,i) = aos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
&BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
implicit none
BEGIN_DOC
! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith ao on the jth grid point
!
! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
@ -107,6 +102,10 @@
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(ao_num,3)
double precision :: aos_lapl_array(ao_num,3)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) &
!$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points)
do m = 1, 3
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
@ -115,7 +114,24 @@
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
do j = 1, ao_num
aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_array(j,m)
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
implicit none
!
! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
integer :: i,j,m
do m = 1, 3
do i = 1, n_points_final_grid
do j = 1, ao_num
aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_in_r_array(j,i,m)
enddo
enddo
enddo

14
src/dft_utils_in_r/dm_in_r.irp.f
View File

@ -29,11 +29,18 @@
double precision, allocatable :: dm_a(:),dm_b(:), dm_a_grad(:,:), dm_b_grad(:,:)
allocate(dm_a(N_states),dm_b(N_states), dm_a_grad(3,N_states), dm_b_grad(3,N_states))
allocate(aos_array(ao_num),grad_aos_array(3,ao_num))
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP SHARED(n_points_final_grid,final_grid_points,N_states, &
!$OMP one_e_dm_and_grad_alpha_in_r,one_e_dm_and_grad_beta_in_r, &
!$OMP one_e_grad_2_dm_alpha_at_r,one_e_grad_2_dm_beta_at_r, &
!$OMP scal_prod_grad_one_e_dm_ab,one_e_stuff_for_pbe) &
!$OMP PRIVATE (istate,i,r,dm_a,dm_b,dm_a_grad,dm_b_grad,aos_array, grad_aos_array)
do istate = 1, N_states
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, dm_a_grad, dm_b_grad, aos_array, grad_aos_array)
@ -72,6 +79,7 @@
* (dm_a(istate) + dm_b(istate))
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

54
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -1,10 +1,7 @@
BEGIN_PROVIDER[double precision, mos_in_r_array, (mo_num,n_points_final_grid)]
&BEGIN_PROVIDER[double precision, mos_in_r_array_transp,(n_points_final_grid,mo_num)]
implicit none
BEGIN_DOC
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
!
! mos_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
END_DOC
integer :: i,j
double precision :: mos_array(mo_num), r(3)
@ -15,14 +12,49 @@
call give_all_mos_at_r(r,mos_array)
do j = 1, mo_num
mos_in_r_array(j,i) = mos_array(j)
mos_in_r_array_transp(i,j) = mos_array(j)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_in_r_array_omp, (mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
END_DOC
integer :: i,j
double precision :: mos_array(mo_num), r(3)
!$OMP PARALLEL DO &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,r,mos_array,j) &
!$OMP SHARED(mos_in_r_array_omp,n_points_final_grid,mo_num,final_grid_points)
do i = 1, n_points_final_grid
r(1) = final_grid_points(1,i)
r(2) = final_grid_points(2,i)
r(3) = final_grid_points(3,i)
call give_all_mos_at_r(r,mos_array)
do j = 1, mo_num
mos_in_r_array_omp(j,i) = mos_array(j)
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_in_r_array_transp,(n_points_final_grid,mo_num)]
implicit none
BEGIN_DOC
! mos_in_r_array_transp(i,j) = value of the jth mo on the ith grid point
END_DOC
integer :: i,j
do i = 1, n_points_final_grid
do j = 1, mo_num
mos_in_r_array_transp(i,j) = mos_in_r_array(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_grad_in_r_array,(mo_num,n_points_final_grid,3)]
&BEGIN_PROVIDER[double precision, mos_grad_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! mos_grad_in_r_array(i,j,k) = value of the kth component of the gradient of ith mo on the jth grid point
@ -32,12 +64,22 @@
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: m
print*,'mo_num,n_points_final_grid',mo_num,n_points_final_grid
mos_grad_in_r_array = 0.d0
do m=1,3
call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_grad_in_r_array(1,1,m),ao_num,0.d0,mos_grad_in_r_array(1,1,m),mo_num)
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_grad_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! mos_grad_in_r_array_transp(i,j,k) = value of the kth component of the gradient of jth mo on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: m
integer :: i,j
mos_grad_in_r_array = 0.d0
do i = 1, n_points_final_grid
do j = 1, mo_num
do m = 1, 3

8
src/dft_utils_one_e/NEED
View File

@ -1,8 +0,0 @@
density_for_dft
dft_utils_in_r
mo_one_e_ints
mo_two_e_ints
ao_one_e_ints
ao_two_e_ints
mo_two_e_erf_ints
ao_two_e_erf_ints

2
src/dummy/NEED
View File

@ -15,7 +15,7 @@ density_for_dft
determinants
dft_keywords
dft_utils_in_r
dft_utils_one_e
dft_utils_func
dressing
electrons
ezfio_files

12
src/ezfio_files/00.create.bats
View File

@ -49,6 +49,18 @@ function run {
run hcn.xyz 1 0 aug-cc-pvdz
}
@test "LiF" {
run lif.xyz 1 0 cc-pvtz
}
@test "F" {
run f.xyz 2 0 cc-pvtz
}
@test "Be" {
run be.xyz 1 0 cc-pvtz
}
@test "N2" {
run n2.xyz 1 0 cc-pvtz
}

6
src/ezfio_files/01.convert.bats
View File

@ -25,3 +25,9 @@ function run {
run cu_nh3_4_2plus.gms.out cu_nh3_4_2plus.ezfio
qp set scf_utils thresh_scf 1.e-10
}
@test "O2 CAS GAMESS" { # 1.38541s
run o2_cas.gms.out o2_cas.gms.ezfio
qp set scf_utils thresh_scf 1.e-10
qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
}

2
src/functionals/NEED
View File

@ -1 +1 @@
dft_utils_one_e
dft_utils_func

12
src/functionals/pbe.irp.f
View File

@ -159,10 +159,10 @@ END_PROVIDER
enddo
do j = 1, ao_num
do m = 1,3
aos_d_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp(m,j,i)
enddo
enddo
enddo
@ -315,8 +315,8 @@ END_PROVIDER
enddo
do j = 1, ao_num
do m = 1,3
aos_d_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp(m,j,i)
enddo
enddo
enddo

172
src/functionals/sr_pbe.irp.f
View File

@ -62,11 +62,11 @@ END_PROVIDER
do istate = 1, n_states
do i = 1, ao_num
do j = 1, ao_num
potential_x_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(i,j,istate)
potential_x_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(i,j,istate)
potential_x_alpha_ao_sr_pbe(j,i,istate) = pot_scal_x_alpha_ao_sr_pbe(j,i,istate) + pot_grad_x_alpha_ao_sr_pbe(j,i,istate) + pot_grad_x_alpha_ao_sr_pbe(i,j,istate)
potential_x_beta_ao_sr_pbe(j,i,istate) = pot_scal_x_beta_ao_sr_pbe(j,i,istate) + pot_grad_x_beta_ao_sr_pbe(j,i,istate) + pot_grad_x_beta_ao_sr_pbe(i,j,istate)
potential_c_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(i,j,istate)
potential_c_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(i,j,istate)
potential_c_alpha_ao_sr_pbe(j,i,istate) = pot_scal_c_alpha_ao_sr_pbe(j,i,istate) + pot_grad_c_alpha_ao_sr_pbe(j,i,istate) + pot_grad_c_alpha_ao_sr_pbe(i,j,istate)
potential_c_beta_ao_sr_pbe(j,i,istate) = pot_scal_c_beta_ao_sr_pbe(j,i,istate) + pot_grad_c_beta_ao_sr_pbe(j,i,istate) + pot_grad_c_beta_ao_sr_pbe(i,j,istate)
enddo
enddo
enddo
@ -83,8 +83,8 @@ END_PROVIDER
do istate = 1, n_states
do i = 1, ao_num
do j = 1, ao_num
potential_xc_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(i,j,istate)
potential_xc_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_xc_beta_ao_pbe(j,i,istate) + pot_sr_grad_xc_beta_ao_pbe(j,i,istate) + pot_sr_grad_xc_beta_ao_pbe(i,j,istate)
potential_xc_alpha_ao_sr_pbe(j,i,istate) = pot_scal_xc_alpha_ao_sr_pbe(j,i,istate) + pot_grad_xc_alpha_ao_sr_pbe(j,i,istate) + pot_grad_xc_alpha_ao_sr_pbe(i,j,istate)
potential_xc_beta_ao_sr_pbe(j,i,istate) = pot_scal_xc_beta_ao_sr_pbe(j,i,istate) + pot_grad_xc_beta_ao_sr_pbe(j,i,istate) + pot_grad_xc_beta_ao_sr_pbe(i,j,istate)
enddo
enddo
enddo
@ -93,19 +93,19 @@ END_PROVIDER
BEGIN_PROVIDER[double precision, aos_sr_vc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_sr_vc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_sr_vx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_sr_vx_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vx_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vx_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
BEGIN_PROVIDER[double precision, aos_vc_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_vc_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_vx_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vc_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vc_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vx_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vx_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
implicit none
BEGIN_DOC
! intermediates to compute the sr_pbe potentials
!
! aos_sr_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
! aos_vxc_alpha_sr_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
END_DOC
integer :: istate,i,j,m
double precision :: mu,weight
@ -114,10 +114,10 @@ END_PROVIDER
double precision :: contrib_grad_xa(3),contrib_grad_xb(3),contrib_grad_ca(3),contrib_grad_cb(3)
double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
aos_dsr_vc_alpha_pbe_w= 0.d0
aos_dsr_vc_beta_pbe_w = 0.d0
aos_dsr_vx_alpha_pbe_w= 0.d0
aos_dsr_vx_beta_pbe_w = 0.d0
aos_d_vc_alpha_sr_pbe_w= 0.d0
aos_d_vc_beta_sr_pbe_w = 0.d0
aos_d_vx_alpha_sr_pbe_w= 0.d0
aos_d_vx_beta_sr_pbe_w = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
weight = final_weight_at_r_vector(i)
@ -150,17 +150,17 @@ END_PROVIDER
contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 * grad_rho_b(m) + vx_grad_rho_a_b * grad_rho_a(m) )
enddo
do j = 1, ao_num
aos_sr_vc_alpha_pbe_w(j,i,istate) = vc_rho_a * aos_in_r_array(j,i)
aos_sr_vc_beta_pbe_w (j,i,istate) = vc_rho_b * aos_in_r_array(j,i)
aos_sr_vx_alpha_pbe_w(j,i,istate) = vx_rho_a * aos_in_r_array(j,i)
aos_sr_vx_beta_pbe_w (j,i,istate) = vx_rho_b * aos_in_r_array(j,i)
aos_vc_alpha_sr_pbe_w(j,i,istate) = vc_rho_a * aos_in_r_array(j,i)
aos_vc_beta_sr_pbe_w (j,i,istate) = vc_rho_b * aos_in_r_array(j,i)
aos_vx_alpha_sr_pbe_w(j,i,istate) = vx_rho_a * aos_in_r_array(j,i)
aos_vx_beta_sr_pbe_w (j,i,istate) = vx_rho_b * aos_in_r_array(j,i)
enddo
do j = 1, ao_num
do m = 1,3
aos_dsr_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_dsr_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_dsr_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_dsr_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vc_alpha_sr_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vc_beta_sr_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vx_alpha_sr_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vx_beta_sr_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp(m,j,i)
enddo
enddo
enddo
@ -169,10 +169,10 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [double precision, pot_sr_scal_x_alpha_ao_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_scal_c_alpha_ao_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_scal_x_beta_ao_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_scal_c_beta_ao_pbe, (ao_num,ao_num,N_states)]
BEGIN_PROVIDER [double precision, pot_scal_x_alpha_ao_sr_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_scal_c_alpha_ao_sr_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_scal_x_beta_ao_sr_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_scal_c_beta_ao_sr_pbe, (ao_num,ao_num,N_states)]
implicit none
! intermediates to compute the sr_pbe potentials
!
@ -180,33 +180,33 @@ END_PROVIDER
BEGIN_DOC
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
END_DOC
pot_sr_scal_c_alpha_ao_pbe = 0.d0
pot_sr_scal_x_alpha_ao_pbe = 0.d0
pot_sr_scal_c_beta_ao_pbe = 0.d0
pot_sr_scal_x_beta_ao_pbe = 0.d0
pot_scal_c_alpha_ao_sr_pbe = 0.d0
pot_scal_x_alpha_ao_sr_pbe = 0.d0
pot_scal_c_beta_ao_sr_pbe = 0.d0
pot_scal_x_beta_ao_sr_pbe = 0.d0
double precision :: wall_1,wall_2
call wall_time(wall_1)
do istate = 1, N_states
! correlation alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vc_alpha_pbe_w(1,1,istate),size(aos_sr_vc_alpha_pbe_w,1), &
aos_vc_alpha_sr_pbe_w(1,1,istate),size(aos_vc_alpha_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_c_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_c_alpha_ao_pbe,1))
pot_scal_c_alpha_ao_sr_pbe(1,1,istate),size(pot_scal_c_alpha_ao_sr_pbe,1))
! correlation beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vc_beta_pbe_w(1,1,istate),size(aos_sr_vc_beta_pbe_w,1), &
aos_vc_beta_sr_pbe_w(1,1,istate),size(aos_vc_beta_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_c_beta_ao_pbe(1,1,istate),size(pot_sr_scal_c_beta_ao_pbe,1))
pot_scal_c_beta_ao_sr_pbe(1,1,istate),size(pot_scal_c_beta_ao_sr_pbe,1))
! exchange alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vx_alpha_pbe_w(1,1,istate),size(aos_sr_vx_alpha_pbe_w,1), &
aos_vx_alpha_sr_pbe_w(1,1,istate),size(aos_vx_alpha_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_x_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_x_alpha_ao_pbe,1))
pot_scal_x_alpha_ao_sr_pbe(1,1,istate),size(pot_scal_x_alpha_ao_sr_pbe,1))
! exchange beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vx_beta_pbe_w(1,1,istate),size(aos_sr_vx_beta_pbe_w,1), &
aos_vx_beta_sr_pbe_w(1,1,istate),size(aos_vx_beta_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_x_beta_ao_pbe(1,1,istate), size(pot_sr_scal_x_beta_ao_pbe,1))
pot_scal_x_beta_ao_sr_pbe(1,1,istate), size(pot_scal_x_beta_ao_sr_pbe,1))
enddo
call wall_time(wall_2)
@ -214,10 +214,10 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [double precision, pot_sr_grad_x_alpha_ao_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_grad_x_beta_ao_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_grad_c_alpha_ao_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_grad_c_beta_ao_pbe,(ao_num,ao_num,N_states)]
BEGIN_PROVIDER [double precision, pot_grad_x_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_grad_x_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_grad_c_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_grad_c_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals
@ -225,31 +225,31 @@ END_PROVIDER
integer :: istate
double precision :: wall_1,wall_2
call wall_time(wall_1)
pot_sr_grad_c_alpha_ao_pbe = 0.d0
pot_sr_grad_x_alpha_ao_pbe = 0.d0
pot_sr_grad_c_beta_ao_pbe = 0.d0
pot_sr_grad_x_beta_ao_pbe = 0.d0
pot_grad_c_alpha_ao_sr_pbe = 0.d0
pot_grad_x_alpha_ao_sr_pbe = 0.d0
pot_grad_c_beta_ao_sr_pbe = 0.d0
pot_grad_x_beta_ao_sr_pbe = 0.d0
do istate = 1, N_states
! correlation alpha
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vc_alpha_pbe_w(1,1,istate),size(aos_dsr_vc_alpha_pbe_w,1), &
aos_d_vc_alpha_sr_pbe_w(1,1,istate),size(aos_d_vc_alpha_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_c_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_c_alpha_ao_pbe,1))
pot_grad_c_alpha_ao_sr_pbe(1,1,istate),size(pot_grad_c_alpha_ao_sr_pbe,1))
! correlation beta
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vc_beta_pbe_w(1,1,istate),size(aos_dsr_vc_beta_pbe_w,1), &
aos_d_vc_beta_sr_pbe_w(1,1,istate),size(aos_d_vc_beta_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_c_beta_ao_pbe(1,1,istate),size(pot_sr_grad_c_beta_ao_pbe,1))
pot_grad_c_beta_ao_sr_pbe(1,1,istate),size(pot_grad_c_beta_ao_sr_pbe,1))
! exchange alpha
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vx_alpha_pbe_w(1,1,istate),size(aos_dsr_vx_alpha_pbe_w,1), &
aos_d_vx_alpha_sr_pbe_w(1,1,istate),size(aos_d_vx_alpha_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_x_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_x_alpha_ao_pbe,1))
pot_grad_x_alpha_ao_sr_pbe(1,1,istate),size(pot_grad_x_alpha_ao_sr_pbe,1))
! exchange beta
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vx_beta_pbe_w(1,1,istate),size(aos_dsr_vx_beta_pbe_w,1), &
aos_d_vx_beta_sr_pbe_w(1,1,istate),size(aos_d_vx_beta_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_x_beta_ao_pbe(1,1,istate),size(pot_sr_grad_x_beta_ao_pbe,1))
pot_grad_x_beta_ao_sr_pbe(1,1,istate),size(pot_grad_x_beta_ao_sr_pbe,1))
enddo
call wall_time(wall_2)
@ -257,13 +257,13 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_sr_vxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_sr_vxc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vxc_alpha_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_dsr_vxc_beta_pbe_w , (ao_num,n_points_final_grid,N_states)]
BEGIN_PROVIDER[double precision, aos_vxc_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_vxc_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vxc_alpha_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
&BEGIN_PROVIDER[double precision, aos_d_vxc_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
implicit none
BEGIN_DOC
! aos_sr_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
! aos_vxc_alpha_sr_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
END_DOC
integer :: istate,i,j,m
double precision :: mu,weight
@ -273,8 +273,8 @@ END_PROVIDER
double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
aos_dsr_vxc_alpha_pbe_w = 0.d0
aos_dsr_vxc_beta_pbe_w = 0.d0
aos_d_vxc_alpha_sr_pbe_w = 0.d0
aos_d_vxc_beta_sr_pbe_w = 0.d0
do istate = 1, N_states
do i = 1, n_points_final_grid
@ -307,13 +307,13 @@ END_PROVIDER
contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 * grad_rho_b(m) + vx_grad_rho_a_b * grad_rho_a(m) )
enddo
do j = 1, ao_num
aos_sr_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i)
aos_sr_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i)
aos_vxc_alpha_sr_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i)
aos_vxc_beta_sr_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i)
enddo
do j = 1, ao_num
do m = 1,3
aos_dsr_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_dsr_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
aos_d_vxc_alpha_sr_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp(m,j,i)
aos_d_vxc_beta_sr_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp(m,j,i)
enddo
enddo
enddo
@ -322,36 +322,36 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [double precision, pot_sr_scal_xc_alpha_ao_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_scal_xc_beta_ao_pbe, (ao_num,ao_num,N_states)]
BEGIN_PROVIDER [double precision, pot_scal_xc_alpha_ao_sr_pbe, (ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_scal_xc_beta_ao_sr_pbe, (ao_num,ao_num,N_states)]
implicit none
integer :: istate
BEGIN_DOC
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
END_DOC
pot_sr_scal_xc_alpha_ao_pbe = 0.d0
pot_sr_scal_xc_beta_ao_pbe = 0.d0
pot_scal_xc_alpha_ao_sr_pbe = 0.d0
pot_scal_xc_beta_ao_sr_pbe = 0.d0
double precision :: wall_1,wall_2
call wall_time(wall_1)
do istate = 1, N_states
! exchange - correlation alpha
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vxc_alpha_pbe_w(1,1,istate),size(aos_sr_vxc_alpha_pbe_w,1), &
aos_vxc_alpha_sr_pbe_w(1,1,istate),size(aos_vxc_alpha_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_xc_alpha_ao_pbe(1,1,istate),size(pot_sr_scal_xc_alpha_ao_pbe,1))
pot_scal_xc_alpha_ao_sr_pbe(1,1,istate),size(pot_scal_xc_alpha_ao_sr_pbe,1))
! exchange - correlation beta
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_sr_vxc_beta_pbe_w(1,1,istate),size(aos_sr_vxc_beta_pbe_w,1), &
aos_vxc_beta_sr_pbe_w(1,1,istate),size(aos_vxc_beta_sr_pbe_w,1), &
aos_in_r_array,size(aos_in_r_array,1),1.d0, &
pot_sr_scal_xc_beta_ao_pbe(1,1,istate),size(pot_sr_scal_xc_beta_ao_pbe,1))
pot_scal_xc_beta_ao_sr_pbe(1,1,istate),size(pot_scal_xc_beta_ao_sr_pbe,1))
enddo
call wall_time(wall_2)
END_PROVIDER
BEGIN_PROVIDER [double precision, pot_sr_grad_xc_alpha_ao_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_sr_grad_xc_beta_ao_pbe,(ao_num,ao_num,N_states)]
BEGIN_PROVIDER [double precision, pot_grad_xc_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
&BEGIN_PROVIDER [double precision, pot_grad_xc_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
implicit none
BEGIN_DOC
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals
@ -359,19 +359,19 @@ END_PROVIDER
integer :: istate
double precision :: wall_1,wall_2
call wall_time(wall_1)
pot_sr_grad_xc_alpha_ao_pbe = 0.d0
pot_sr_grad_xc_beta_ao_pbe = 0.d0
pot_grad_xc_alpha_ao_sr_pbe = 0.d0
pot_grad_xc_beta_ao_sr_pbe = 0.d0
do istate = 1, N_states
! exchange - correlation alpha
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vxc_alpha_pbe_w(1,1,istate),size(aos_dsr_vxc_alpha_pbe_w,1), &
aos_d_vxc_alpha_sr_pbe_w(1,1,istate),size(aos_d_vxc_alpha_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_xc_alpha_ao_pbe(1,1,istate),size(pot_sr_grad_xc_alpha_ao_pbe,1))
pot_grad_xc_alpha_ao_sr_pbe(1,1,istate),size(pot_grad_xc_alpha_ao_sr_pbe,1))
! exchange - correlation beta
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
aos_dsr_vxc_beta_pbe_w(1,1,istate),size(aos_dsr_vxc_beta_pbe_w,1), &
aos_d_vxc_beta_sr_pbe_w(1,1,istate),size(aos_d_vxc_beta_sr_pbe_w,1), &
aos_in_r_array_transp,size(aos_in_r_array_transp,1),1.d0, &
pot_sr_grad_xc_beta_ao_pbe(1,1,istate),size(pot_sr_grad_xc_beta_ao_pbe,1))
pot_grad_xc_beta_ao_sr_pbe(1,1,istate),size(pot_grad_xc_beta_ao_sr_pbe,1))
enddo
call wall_time(wall_2)

13
src/hartree_fock/10.hf.bats
View File

@ -21,6 +21,19 @@ function run() {
run b2_stretched.ezfio -48.9950585434279
}
@test "LiF" { # 3 s
run lif.ezfio -106.9801081911955
}
@test "Be" { # 3 s
run be.ezfio -14.57287346825270
}
@test "F" { # 3 s
run f.ezfio -99.40093527229389
}
@test "SiH2_3B1" { # 0.539000 1.51094s
run sih2_3b1.ezfio -289.9654718453571
}

20
src/mo_basis/mos_in_r.irp.f
View File

@ -37,18 +37,18 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
integer :: i,j,k
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
mos_array=0d0
mos_grad_array=0d0
mos_lapl_array=0d0
mos_array = 0.d0
mos_grad_array = 0.d0
mos_lapl_array = 0.d0
do j = 1, mo_num
do k=1, ao_num
mos_array(j) += mo_coef(k,j)*aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
mos_lapl_array(j,1) += mo_coef(k,j)*aos_lapl_array(k,1)
mos_lapl_array(j,2) += mo_coef(k,j)*aos_lapl_array(k,2)
mos_lapl_array(j,3) += mo_coef(k,j)*aos_lapl_array(k,3)
mos_array(j) += mo_coef(k,j) * aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j) * aos_grad_array(k,1)
mos_grad_array(j,2) += mo_coef(k,j) * aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3)
mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1)
mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2)
mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3)
enddo
enddo
end

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

1
tests/input/be.xyz Normal file
View File

@ -0,0 +1 @@
Be

1
tests/input/f.xyz Normal file
View File

@ -0,0 +1 @@
F

4
tests/input/lif.xyz Normal file
View File

@ -0,0 +1,4 @@
2
Li 0. 0. 1.56359565
F 0. 0. 0.

2385
tests/input/o2_cas.gms.out Normal file
View File

File diff suppressed because it is too large Load Diff