mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-09-12 02:38:31 +02:00
Merge branch 'dev' of https://github.com/QuantumPackage/qp2 into dev
This commit is contained in:
commit
cbaeecb11f
@ -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)
|
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.
|
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
63
src/basis_correction/51.basis_c.bats
Normal 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
3
src/basis_correction/NEED
Normal file
@ -0,0 +1,3 @@
|
|||||||
|
mu_of_r
|
||||||
|
dft_utils_func
|
||||||
|
dft_one_e
|
27
src/basis_correction/README.rst
Normal file
27
src/basis_correction/README.rst
Normal 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
1
src/basis_correction/TODO
Normal file
@ -0,0 +1 @@
|
|||||||
|
change all correlation functionals with the pbe_on_top general
|
27
src/basis_correction/basis_correction.irp.f
Normal file
27
src/basis_correction/basis_correction.irp.f
Normal 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
92
src/basis_correction/eff_xi_based_func.irp.f
Normal 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
183
src/basis_correction/pbe_on_top.irp.f
Normal 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
82
src/basis_correction/print_routine.irp.f
Normal 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
83
src/basis_correction/weak_corr_func.irp.f
Normal 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
|
@ -18,6 +18,7 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
|
|||||||
enddo
|
enddo
|
||||||
print*,'n_points_final_grid = ',n_points_final_grid
|
print*,'n_points_final_grid = ',n_points_final_grid
|
||||||
print*,'n max point = ',n_points_integration_angular*(n_points_radial_grid*nucl_num - 1)
|
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
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)]
|
BEGIN_PROVIDER [double precision, final_grid_points, (3,n_points_final_grid)]
|
||||||
|
2
src/cas_based_on_top/NEED
Normal file
2
src/cas_based_on_top/NEED
Normal file
@ -0,0 +1,2 @@
|
|||||||
|
two_body_rdm
|
||||||
|
dft_utils_in_r
|
9
src/cas_based_on_top/README.rst
Normal file
9
src/cas_based_on_top/README.rst
Normal 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
125
src/cas_based_on_top/c_i_a_v_mos.irp.f
Normal 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
47
src/cas_based_on_top/cas_based_density.irp.f
Normal 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
19
src/cas_based_on_top/cas_based_on_top.irp.f
Normal 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
99
src/cas_based_on_top/cas_dens_prov.irp.f
Normal 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
65
src/cas_based_on_top/cas_dens_rout.irp.f
Normal 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
37
src/cas_based_on_top/cas_one_e_rdm.irp.f
Normal 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
68
src/cas_based_on_top/eff_spin_dens.irp.f
Normal 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
41
src/cas_based_on_top/example.irp.f
Normal 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
76
src/cas_based_on_top/on_top_cas_prov.irp.f
Normal 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
118
src/cas_based_on_top/on_top_cas_rout.irp.f
Normal 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
73
src/cas_based_on_top/on_top_grad.irp.f
Normal 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
|
||||||
|
|
||||||
|
|
@ -1,2 +1,8 @@
|
|||||||
dft_utils_one_e
|
dft_utils_func
|
||||||
functionals
|
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
|
||||||
|
2
src/dft_utils_func/NEED
Normal file
2
src/dft_utils_func/NEED
Normal file
@ -0,0 +1,2 @@
|
|||||||
|
density_for_dft
|
||||||
|
dft_utils_in_r
|
152
src/dft_utils_func/ecmd_lda.irp.f
Normal file
152
src/dft_utils_func/ecmd_lda.irp.f
Normal 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
53
src/dft_utils_func/ecmd_pbe_general.irp.f
Normal 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
72
src/dft_utils_func/ecmd_pbe_on_top.irp.f
Normal 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
194
src/dft_utils_func/ecmd_pbe_ueg.irp.f
Normal 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
|
||||||
|
|
121
src/dft_utils_func/on_top_from_ueg.irp.f
Normal file
121
src/dft_utils_func/on_top_from_ueg.irp.f
Normal 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
|
||||||
|
|
@ -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, (ao_num,n_points_final_grid)]
|
||||||
&BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_num)]
|
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! aos_in_r_array(i,j) = value of the ith ao on the jth grid point
|
! 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
|
END_DOC
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
double precision :: aos_array(ao_num), r(3)
|
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
|
do i = 1, n_points_final_grid
|
||||||
r(1) = final_grid_points(1,i)
|
r(1) = final_grid_points(1,i)
|
||||||
r(2) = final_grid_points(2,i)
|
r(2) = final_grid_points(2,i)
|
||||||
@ -15,11 +17,30 @@
|
|||||||
call give_all_aos_at_r(r,aos_array)
|
call give_all_aos_at_r(r,aos_array)
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
aos_in_r_array(j,i) = aos_array(j)
|
aos_in_r_array(j,i) = aos_array(j)
|
||||||
aos_in_r_array_transp(i,j) = aos_array(j)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
END_PROVIDER
|
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)]
|
BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_grid,3)]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
@ -30,6 +51,10 @@
|
|||||||
integer :: i,j,m
|
integer :: i,j,m
|
||||||
double precision :: aos_array(ao_num), r(3)
|
double precision :: aos_array(ao_num), r(3)
|
||||||
double precision :: aos_grad_array(3,ao_num)
|
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
|
do i = 1, n_points_final_grid
|
||||||
r(1) = final_grid_points(1,i)
|
r(1) = final_grid_points(1,i)
|
||||||
r(2) = final_grid_points(2,i)
|
r(2) = final_grid_points(2,i)
|
||||||
@ -41,15 +66,16 @@
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
implicit none
|
||||||
BEGIN_DOC
|
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
|
! k = 1 : x, k= 2, y, k 3, z
|
||||||
END_DOC
|
END_DOC
|
||||||
@ -57,49 +83,18 @@
|
|||||||
double precision :: aos_array(ao_num), r(3)
|
double precision :: aos_array(ao_num), r(3)
|
||||||
double precision :: aos_grad_array(3,ao_num)
|
double precision :: aos_grad_array(3,ao_num)
|
||||||
do i = 1, n_points_final_grid
|
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 m = 1, 3
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
aos_grad_in_r_array_transp(i,j,m) = 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_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)
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
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, (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
|
implicit none
|
||||||
BEGIN_DOC
|
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(i,j,k) = value of the kth component of the laplacian of jth ao on the ith 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
|
|
||||||
!
|
!
|
||||||
! k = 1 : x, k= 2, y, k 3, z
|
! k = 1 : x, k= 2, y, k 3, z
|
||||||
END_DOC
|
END_DOC
|
||||||
@ -107,6 +102,10 @@
|
|||||||
double precision :: aos_array(ao_num), r(3)
|
double precision :: aos_array(ao_num), r(3)
|
||||||
double precision :: aos_grad_array(ao_num,3)
|
double precision :: aos_grad_array(ao_num,3)
|
||||||
double precision :: aos_lapl_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 m = 1, 3
|
||||||
do i = 1, n_points_final_grid
|
do i = 1, n_points_final_grid
|
||||||
r(1) = final_grid_points(1,i)
|
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)
|
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
|
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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
@ -29,11 +29,18 @@
|
|||||||
double precision, allocatable :: dm_a(:),dm_b(:), dm_a_grad(:,:), dm_b_grad(:,:)
|
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(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))
|
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 istate = 1, N_states
|
||||||
do i = 1, n_points_final_grid
|
do i = 1, n_points_final_grid
|
||||||
r(1) = final_grid_points(1,i)
|
r(1) = final_grid_points(1,i)
|
||||||
r(2) = final_grid_points(2,i)
|
r(2) = final_grid_points(2,i)
|
||||||
r(3) = final_grid_points(3,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)
|
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))
|
* (dm_a(istate) + dm_b(istate))
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -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, (mo_num,n_points_final_grid)]
|
||||||
&BEGIN_PROVIDER[double precision, mos_in_r_array_transp,(n_points_final_grid,mo_num)]
|
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! mos_in_r_array(i,j) = value of the ith mo on the jth grid point
|
! 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
|
END_DOC
|
||||||
integer :: i,j
|
integer :: i,j
|
||||||
double precision :: mos_array(mo_num), r(3)
|
double precision :: mos_array(mo_num), r(3)
|
||||||
@ -15,14 +12,49 @@
|
|||||||
call give_all_mos_at_r(r,mos_array)
|
call give_all_mos_at_r(r,mos_array)
|
||||||
do j = 1, mo_num
|
do j = 1, mo_num
|
||||||
mos_in_r_array(j,i) = mos_array(j)
|
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
|
||||||
enddo
|
enddo
|
||||||
END_PROVIDER
|
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,(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
|
implicit none
|
||||||
BEGIN_DOC
|
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
|
! 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
|
! k = 1 : x, k= 2, y, k 3, z
|
||||||
END_DOC
|
END_DOC
|
||||||
integer :: m
|
integer :: m
|
||||||
print*,'mo_num,n_points_final_grid',mo_num,n_points_final_grid
|
|
||||||
mos_grad_in_r_array = 0.d0
|
mos_grad_in_r_array = 0.d0
|
||||||
do m=1,3
|
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)
|
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
|
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
|
integer :: i,j
|
||||||
|
mos_grad_in_r_array = 0.d0
|
||||||
do i = 1, n_points_final_grid
|
do i = 1, n_points_final_grid
|
||||||
do j = 1, mo_num
|
do j = 1, mo_num
|
||||||
do m = 1, 3
|
do m = 1, 3
|
||||||
|
@ -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
|
|
@ -15,7 +15,7 @@ density_for_dft
|
|||||||
determinants
|
determinants
|
||||||
dft_keywords
|
dft_keywords
|
||||||
dft_utils_in_r
|
dft_utils_in_r
|
||||||
dft_utils_one_e
|
dft_utils_func
|
||||||
dressing
|
dressing
|
||||||
electrons
|
electrons
|
||||||
ezfio_files
|
ezfio_files
|
||||||
|
@ -49,6 +49,18 @@ function run {
|
|||||||
run hcn.xyz 1 0 aug-cc-pvdz
|
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" {
|
@test "N2" {
|
||||||
run n2.xyz 1 0 cc-pvtz
|
run n2.xyz 1 0 cc-pvtz
|
||||||
}
|
}
|
||||||
|
@ -25,3 +25,9 @@ function run {
|
|||||||
run cu_nh3_4_2plus.gms.out cu_nh3_4_2plus.ezfio
|
run cu_nh3_4_2plus.gms.out cu_nh3_4_2plus.ezfio
|
||||||
qp set scf_utils thresh_scf 1.e-10
|
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]"
|
||||||
|
}
|
||||||
|
@ -1 +1 @@
|
|||||||
dft_utils_one_e
|
dft_utils_func
|
||||||
|
@ -159,10 +159,10 @@ END_PROVIDER
|
|||||||
enddo
|
enddo
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
do m = 1,3
|
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_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_xyz(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_xyz(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_xyz(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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -315,8 +315,8 @@ END_PROVIDER
|
|||||||
enddo
|
enddo
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
do m = 1,3
|
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_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_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(m,j,i)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
@ -62,11 +62,11 @@ END_PROVIDER
|
|||||||
do istate = 1, n_states
|
do istate = 1, n_states
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
do j = 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_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_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_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_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_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_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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -83,8 +83,8 @@ END_PROVIDER
|
|||||||
do istate = 1, n_states
|
do istate = 1, n_states
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
do j = 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_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_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_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
|
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_vc_alpha_sr_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_vc_beta_sr_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_vx_alpha_sr_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_vx_beta_sr_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_d_vc_alpha_sr_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_d_vc_beta_sr_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_d_vx_alpha_sr_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_d_vx_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! intermediates to compute the sr_pbe potentials
|
! 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
|
END_DOC
|
||||||
integer :: istate,i,j,m
|
integer :: istate,i,j,m
|
||||||
double precision :: mu,weight
|
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 :: 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 :: 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
|
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_d_vc_alpha_sr_pbe_w= 0.d0
|
||||||
aos_dsr_vc_beta_pbe_w = 0.d0
|
aos_d_vc_beta_sr_pbe_w = 0.d0
|
||||||
aos_dsr_vx_alpha_pbe_w= 0.d0
|
aos_d_vx_alpha_sr_pbe_w= 0.d0
|
||||||
aos_dsr_vx_beta_pbe_w = 0.d0
|
aos_d_vx_beta_sr_pbe_w = 0.d0
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
do i = 1, n_points_final_grid
|
do i = 1, n_points_final_grid
|
||||||
weight = final_weight_at_r_vector(i)
|
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) )
|
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
|
enddo
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
aos_sr_vc_alpha_pbe_w(j,i,istate) = vc_rho_a * 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_sr_vc_beta_pbe_w (j,i,istate) = vc_rho_b * 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_sr_vx_alpha_pbe_w(j,i,istate) = vx_rho_a * 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_sr_vx_beta_pbe_w (j,i,istate) = vx_rho_b * 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
|
enddo
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
do m = 1,3
|
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_d_vc_alpha_sr_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp(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_d_vc_beta_sr_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp(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_d_vx_alpha_sr_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp(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_vx_beta_sr_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp(m,j,i)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -169,10 +169,10 @@ END_PROVIDER
|
|||||||
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_scal_x_alpha_ao_sr_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_scal_c_alpha_ao_sr_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_scal_x_beta_ao_sr_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_c_beta_ao_sr_pbe, (ao_num,ao_num,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
! intermediates to compute the sr_pbe potentials
|
! intermediates to compute the sr_pbe potentials
|
||||||
!
|
!
|
||||||
@ -180,33 +180,33 @@ END_PROVIDER
|
|||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
|
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
|
||||||
END_DOC
|
END_DOC
|
||||||
pot_sr_scal_c_alpha_ao_pbe = 0.d0
|
pot_scal_c_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_scal_x_alpha_ao_pbe = 0.d0
|
pot_scal_x_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_scal_c_beta_ao_pbe = 0.d0
|
pot_scal_c_beta_ao_sr_pbe = 0.d0
|
||||||
pot_sr_scal_x_beta_ao_pbe = 0.d0
|
pot_scal_x_beta_ao_sr_pbe = 0.d0
|
||||||
double precision :: wall_1,wall_2
|
double precision :: wall_1,wall_2
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
! correlation alpha
|
! correlation alpha
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! correlation beta
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange alpha
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange beta
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
enddo
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
@ -214,10 +214,10 @@ END_PROVIDER
|
|||||||
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_grad_x_alpha_ao_sr_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_grad_x_beta_ao_sr_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_grad_c_alpha_ao_sr_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_c_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals
|
! 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
|
integer :: istate
|
||||||
double precision :: wall_1,wall_2
|
double precision :: wall_1,wall_2
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
pot_sr_grad_c_alpha_ao_pbe = 0.d0
|
pot_grad_c_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_grad_x_alpha_ao_pbe = 0.d0
|
pot_grad_x_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_grad_c_beta_ao_pbe = 0.d0
|
pot_grad_c_beta_ao_sr_pbe = 0.d0
|
||||||
pot_sr_grad_x_beta_ao_pbe = 0.d0
|
pot_grad_x_beta_ao_sr_pbe = 0.d0
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
! correlation alpha
|
! correlation alpha
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! correlation beta
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange alpha
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange beta
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
enddo
|
||||||
|
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
@ -257,13 +257,13 @@ END_PROVIDER
|
|||||||
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_vxc_alpha_sr_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_vxc_beta_sr_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_d_vxc_alpha_sr_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_d_vxc_beta_sr_pbe_w , (ao_num,n_points_final_grid,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
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
|
END_DOC
|
||||||
integer :: istate,i,j,m
|
integer :: istate,i,j,m
|
||||||
double precision :: mu,weight
|
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 :: 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
|
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_d_vxc_alpha_sr_pbe_w = 0.d0
|
||||||
aos_dsr_vxc_beta_pbe_w = 0.d0
|
aos_d_vxc_beta_sr_pbe_w = 0.d0
|
||||||
|
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
do i = 1, n_points_final_grid
|
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) )
|
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
|
enddo
|
||||||
do j = 1, ao_num
|
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_vxc_alpha_sr_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_beta_sr_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i)
|
||||||
enddo
|
enddo
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
do m = 1,3
|
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_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_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_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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -322,36 +322,36 @@ END_PROVIDER
|
|||||||
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_scal_xc_alpha_ao_sr_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_beta_ao_sr_pbe, (ao_num,ao_num,N_states)]
|
||||||
implicit none
|
implicit none
|
||||||
integer :: istate
|
integer :: istate
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
|
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the scalar part of the potential
|
||||||
END_DOC
|
END_DOC
|
||||||
pot_sr_scal_xc_alpha_ao_pbe = 0.d0
|
pot_scal_xc_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_scal_xc_beta_ao_pbe = 0.d0
|
pot_scal_xc_beta_ao_sr_pbe = 0.d0
|
||||||
double precision :: wall_1,wall_2
|
double precision :: wall_1,wall_2
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
! exchange - correlation alpha
|
! exchange - correlation alpha
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange - correlation beta
|
||||||
call dgemm('N','T',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
enddo
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, pot_sr_grad_xc_alpha_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_sr_grad_xc_beta_ao_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
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
! intermediate quantity for the calculation of the vxc potentials for the GGA functionals related to the gradienst of the density and orbitals
|
! 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
|
integer :: istate
|
||||||
double precision :: wall_1,wall_2
|
double precision :: wall_1,wall_2
|
||||||
call wall_time(wall_1)
|
call wall_time(wall_1)
|
||||||
pot_sr_grad_xc_alpha_ao_pbe = 0.d0
|
pot_grad_xc_alpha_ao_sr_pbe = 0.d0
|
||||||
pot_sr_grad_xc_beta_ao_pbe = 0.d0
|
pot_grad_xc_beta_ao_sr_pbe = 0.d0
|
||||||
do istate = 1, N_states
|
do istate = 1, N_states
|
||||||
! exchange - correlation alpha
|
! exchange - correlation alpha
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
! exchange - correlation beta
|
||||||
call dgemm('N','N',ao_num,ao_num,n_points_final_grid,1.d0, &
|
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, &
|
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
|
enddo
|
||||||
|
|
||||||
call wall_time(wall_2)
|
call wall_time(wall_2)
|
||||||
|
@ -21,6 +21,19 @@ function run() {
|
|||||||
run b2_stretched.ezfio -48.9950585434279
|
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
|
@test "SiH2_3B1" { # 0.539000 1.51094s
|
||||||
run sih2_3b1.ezfio -289.9654718453571
|
run sih2_3b1.ezfio -289.9654718453571
|
||||||
}
|
}
|
||||||
|
@ -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
|
integer :: i,j,k
|
||||||
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
|
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)
|
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
|
||||||
mos_array=0d0
|
mos_array = 0.d0
|
||||||
mos_grad_array=0d0
|
mos_grad_array = 0.d0
|
||||||
mos_lapl_array=0d0
|
mos_lapl_array = 0.d0
|
||||||
do j = 1, mo_num
|
do j = 1, mo_num
|
||||||
do k=1, ao_num
|
do k=1, ao_num
|
||||||
mos_array(j) += mo_coef(k,j)*aos_array(k)
|
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,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,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_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,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,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_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
end
|
end
|
||||||
|
18
src/mu_of_r/EZFIO.cfg
Normal file
18
src/mu_of_r/EZFIO.cfg
Normal 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
1
src/mu_of_r/NEED
Normal file
@ -0,0 +1 @@
|
|||||||
|
cas_based_on_top
|
4
src/mu_of_r/README.rst
Normal file
4
src/mu_of_r/README.rst
Normal file
@ -0,0 +1,4 @@
|
|||||||
|
==================
|
||||||
|
mu_of_r_definition
|
||||||
|
==================
|
||||||
|
|
106
src/mu_of_r/basis_def.irp.f
Normal file
106
src/mu_of_r/basis_def.irp.f
Normal 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
202
src/mu_of_r/example.irp.f
Normal 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
138
src/mu_of_r/f_hf_utils.irp.f
Normal 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
313
src/mu_of_r/f_psi_i_a_v_utils.irp.f
Normal 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
39
src/mu_of_r/f_psi_old.irp.f
Normal 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
151
src/mu_of_r/f_psi_utils.irp.f
Normal 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
90
src/mu_of_r/f_val_general.irp.f
Normal 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
150
src/mu_of_r/mu_of_r_conditions.irp.f
Normal 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
21
src/mu_of_r/test_proj_op.irp.f
Normal 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
1
tests/input/be.xyz
Normal file
@ -0,0 +1 @@
|
|||||||
|
Be
|
1
tests/input/f.xyz
Normal file
1
tests/input/f.xyz
Normal file
@ -0,0 +1 @@
|
|||||||
|
F
|
4
tests/input/lif.xyz
Normal file
4
tests/input/lif.xyz
Normal file
@ -0,0 +1,4 @@
|
|||||||
|
2
|
||||||
|
|
||||||
|
Li 0. 0. 1.56359565
|
||||||
|
F 0. 0. 0.
|
2385
tests/input/o2_cas.gms.out
Normal file
2385
tests/input/o2_cas.gms.out
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user