LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

Fix basis correction

Browse Source
This commit is contained in:
Anthony Scemama 2025-03-20 11:49:21 +01:00
parent 854b917735
commit 8ca5bc68da
1 changed files with 13 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
plugins/local/basis_correction/weak_corr_func.irp.f
View File

@ -10,21 +10,28 @@
double precision :: rho_a, rho_b, ec, rho, p2 double precision :: rho_a, rho_b, ec, rho, p2
double precision :: wall0,wall1,weight,mu double precision :: wall0,wall1,weight,mu
logical :: dospin logical :: dospin
double precision, external :: g0_UEG_mu_inf
dospin = .true. ! JT dospin have to be set to true for open shell dospin = .true. ! JT dospin have to be set to true for open shell
print*,'Providing ecmd_lda_mu_of_r ...' print*,'Providing ecmd_lda_mu_of_r ...'
ecmd_lda_mu_of_r = 0.d0 ecmd_lda_mu_of_r = 0.d0
call wall_time(wall0) call wall_time(wall0)
if (mu_of_r_potential.EQ."proj")then if (mu_of_r_potential.EQ."proj_DUMMY")then
do istate = 1, N_states do istate = 1, N_states
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018) ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
mu = mu_of_r_prov(ipoint,istate) mu = mu_of_r_prov(ipoint,istate)
weight = final_weight_at_r_vector(ipoint) weight = final_weight_at_r_vector(ipoint)
rho = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
p2 = !TODO rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
rho_a = 0.5d0*(rho - dsqrt(-p2 + rho*rho)) rho = rho_a + rho_b
rho_b = 0.5d0*(rho + dsqrt(-p2 + rho*rho)) ! p2 = rho_a*rho_b
! We take the on-top pair density of the UEG which is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
p2 = 4.d0 * rho_a * rho_b * g0_UEG_mu_inf(rho_a,rho_b)
rho_a = 0.5d0*(rho + dsqrt(rho*rho - 4.d0*p2))
rho_b = 0.5d0*(rho - dsqrt(rho*rho - 4.d0*p2))
! rho_a = 0.5d0*rho
! rho_b = 0.5d0*rho
! Ecmd within the LDA approximation of PRB 73, 155111 (2006) ! Ecmd within the LDA approximation of PRB 73, 155111 (2006)
call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec) call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec)
if(isnan(ec))then if(isnan(ec))then
@ -71,7 +78,7 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
double precision :: weight double precision :: weight
integer :: ipoint,istate 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 :: 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 double precision, external :: g0_UEG_mu_inf
ecmd_pbe_ueg_mu_of_r = 0.d0 ecmd_pbe_ueg_mu_of_r = 0.d0