diff --git a/plugins/local/basis_correction/weak_corr_func.irp.f b/plugins/local/basis_correction/weak_corr_func.irp.f index a31ff3af..99ea504b 100644 --- a/plugins/local/basis_correction/weak_corr_func.irp.f +++ b/plugins/local/basis_correction/weak_corr_func.irp.f @@ -10,21 +10,28 @@ double precision :: rho_a, rho_b, ec, rho, p2 double precision :: wall0,wall1,weight,mu logical :: dospin + double precision, external :: g0_UEG_mu_inf 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) - if (mu_of_r_potential.EQ."proj")then + if (mu_of_r_potential.EQ."proj_DUMMY")then 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 = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) - p2 = !TODO - rho_a = 0.5d0*(rho - dsqrt(-p2 + rho*rho)) - rho_b = 0.5d0*(rho + dsqrt(-p2 + rho*rho)) + 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) + rho = rho_a + rho_b +! 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) call ESRC_MD_LDAERF (mu,rho_a,rho_b,dospin,ec) if(isnan(ec))then @@ -71,7 +78,7 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)] 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 + double precision, external :: g0_UEG_mu_inf ecmd_pbe_ueg_mu_of_r = 0.d0