From f877b41ef545b56bcdf85184531e2e9934175c44 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner LCT Date: Thu, 13 Jun 2019 19:07:01 +0200 Subject: [PATCH] added kinetic density --- src/dft_utils_in_r/kin_dens.irp.f | 15 +++++++++++++++ src/dft_utils_one_e/ec_scan.irp.f | 32 +++++++++++++++---------------- 2 files changed, 30 insertions(+), 17 deletions(-) create mode 100644 src/dft_utils_in_r/kin_dens.irp.f diff --git a/src/dft_utils_in_r/kin_dens.irp.f b/src/dft_utils_in_r/kin_dens.irp.f new file mode 100644 index 00000000..008d6b2d --- /dev/null +++ b/src/dft_utils_in_r/kin_dens.irp.f @@ -0,0 +1,15 @@ +BEGIN_PROVIDER [double precision, kinetic_density_generalized, (n_points_final_grid)] + implicit none + integer :: i,j,m,i_point + kinetic_density_generalized = 0.d0 + do i_point = 1, n_points_final_grid + do i = 1, mo_num + do j = 1, mo_num + do m = 1, 3 + kinetic_density_generalized(i_point) += 0.5d0 * mos_grad_in_r_array_tranp(m,j,i_point) * mos_grad_in_r_array_tranp(m,i,i_point) * one_e_dm_mo_for_dft(j,i,1) + enddo + enddo + enddo + enddo + +END_PROVIDER diff --git a/src/dft_utils_one_e/ec_scan.irp.f b/src/dft_utils_one_e/ec_scan.irp.f index 0aa4b5ca..7a4b587b 100644 --- a/src/dft_utils_one_e/ec_scan.irp.f +++ b/src/dft_utils_one_e/ec_scan.irp.f @@ -23,7 +23,11 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2) cst_3pi2 = 3.d0 * pi*pi drho2 = max(grad_rho_2,thr) drho = dsqrt(drho2) - spin_d = max(nup-ndo,thr) + if((nup-ndo).gt.0.d0)then + spin_d = max(nup-ndo,thr) + else + spin_d = min(nup-ndo,-thr) + endif c_1c = 0.64d0 c_2c = 1.5d0 d_c = 0.7d0 @@ -37,20 +41,24 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2) xi = spin_d/rho rs = (cst_43 * pi * rho)**(-cst_13) s = drho/( 2.d0 * cst_3pi2**(cst_13) * rho**cst_43 ) - t_w = drho2 * cst_18 + t_w = drho2 * cst_18 * rho_inv ds_xi = 0.5d0 * ( (1.d0+xi)**cst_53 + (1.d0 - xi)**cst_53) t_unif = 0.3d0 * (cst_3pi2)**cst_23 * rho**cst_53*ds_xi t_unif = max(t_unif,thr) alpha = (tau - t_w)/t_unif cst_1alph= 1.d0 - alpha - cst_1alph= max(cst_1alph,thr) + if(cst_1alph.gt.0.d0)then + cst_1alph= max(cst_1alph,thr) + else + cst_1alph= min(cst_1alph,-thr) + endif inv_1alph= 1.d0/cst_1alph phi = 0.5d0 * ( (1.d0+xi)**cst_23 + (1.d0 - xi)**cst_23) phi_3 = phi*phi*phi t = (cst_3pi2/16.d0)**cst_13 * s / (phi * rs**0.5d0) - w_1 = dexp(-e_c_lsda1/(gama * phi_3)) + w_1 = dexp(-e_c_lsda1/(gama * phi_3)) - 1.d0 a = beta_rs(rs) /(gama * w_1) - g_at2 = 1.d0/(1.d0 + 4.d0 * a*t*t) + g_at2 = 1.d0/(1.d0 + 4.d0 * a*t*t)**0.25d0 h1 = gama * phi_3 * dlog(1.d0 + w_1 * (1.d0 - g_at2)) ! interpolation function fc_alpha = dexp(-c_1c * alpha * inv_1alph) * step_f(cst_1alph) - d_c * dexp(c_2c * inv_1alph) * step_f(-cst_1alph) @@ -64,24 +72,14 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2) beta_inf = 0.066725d0 * 0.1d0 / 0.1778d0 cx_xi = -3.d0/(4.d0*pi) * (9.d0 * pi/4.d0)**cst_13 * dx_xi - if(dabs(xi).lt.thr)then - x_inf = 0.128026d0 - else - x_inf = (cst_3pi2/16.d0)**cst_23 * beta_inf * phi/(cx_xi - f0) - endif + x_inf = 0.128026d0 f0 = -0.9d0 g_inf = 1.d0/(1.d0 + 4.d0 * x_inf * s*s)**0.25d0 - h0 = b_1c * dlog(1.d0 + w_0 * (1.d0 * g_inf)) + h0 = b_1c * dlog(1.d0 + w_0 * (1.d0 - g_inf)) e_c_0 = (e_c_lsda0 + h0) * gc_xi ec_scan = e_c_1 + fc_alpha * (e_c_0 - e_c_1) - if(isnan(ec_scan))then - print*,'isnan(ec_scan)' - print*,'e_c_1 + fc_alpha * e_c_0' - print*, e_c_1 , fc_alpha , e_c_0 - stop - endif end double precision function step_f(x)