TABLE OF CONTENTS
- ABINIT/m_xc_vdw
- ABINIT/m_xc_vdw/vdw_df_create_mesh
- ABINIT/m_xc_vdw/vdw_df_filter
- ABINIT/m_xc_vdw/vdw_df_indexof
- ABINIT/m_xc_vdw/vdw_df_internal_checks
- ABINIT/m_xc_vdw/vdw_df_interpolate
- ABINIT/m_xc_vdw/vdw_df_kernel
- ABINIT/m_xc_vdw/vdw_df_kernel_value
- ABINIT/m_xc_vdw/vdw_df_ldaxc
- ABINIT/m_xc_vdw/vdw_df_set_tweaks
- ABINIT/m_xc_vdw/vdw_df_tweaks_type
- ABINIT/m_xc_vdw/vdw_df_write_func
- ABINIT/m_xc_vdw/xc_vdw_aggregate
- ABINIT/m_xc_vdw/xc_vdw_done
- ABINIT/m_xc_vdw/xc_vdw_energy
- ABINIT/m_xc_vdw/xc_vdw_get_params
- ABINIT/m_xc_vdw/xc_vdw_init
- ABINIT/m_xc_vdw/xc_vdw_libxc_init
- ABINIT/m_xc_vdw/xc_vdw_memcheck
- ABINIT/m_xc_vdw/xc_vdw_read
- ABINIT/m_xc_vdw/xc_vdw_set_functional
- ABINIT/m_xc_vdw/xc_vdw_set_params
- ABINIT/m_xc_vdw/xc_vdw_show
- ABINIT/m_xc_vdw/xc_vdw_status
- ABINIT/m_xc_vdw/xc_vdw_trigger
- ABINIT/m_xc_vdw/xc_vdw_type
- ABINIT/m_xc_vdw/xc_vdw_write
- m_xc_vdw/vdw_df_netcdf_ioerr
- m_xc_vdw/vdw_df_saturation
ABINIT/m_xc_vdw [ Modules ]
NAME
m_xc_vdw
FUNCTION
Calculates van der Waals corrections to exchange-correlation.
COPYRIGHT
Copyright (C) 2010-2024 ABINIT group (Yann Pouillon, Camilo Espejo) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
DRSLL04 = doi:10.1103/PhysRevLett.92.246401 [[cite:Dion2004]] RS09 = doi:10.1103/PhysRevLett.103.096102 [[cite:Romanperez2009]]
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 #include "abi_xc_vdw.h" 26 27 module m_xc_vdw 28 29 use defs_basis 30 use, intrinsic :: iso_c_binding 31 use m_abicore 32 use m_errors 33 use libxc_functionals 34 use m_splines 35 #ifdef HAVE_NETCDF 36 use netcdf 37 #endif 38 39 use m_io_tools, only : flush_unit, open_file 40 use m_numeric_tools, only : simpson_int, cspint 41 use m_integrals, only : radsintr 42 43 implicit none 44 45 #if defined DEV_YP_VDWXC 46 47 public :: & 48 & xc_vdw_type, & 49 & xc_vdw_aggregate, & 50 & xc_vdw_energy, & 51 & xc_vdw_done, & 52 & xc_vdw_get_params, & 53 & xc_vdw_init, & 54 & xc_vdw_libxc_init, & 55 & xc_vdw_memcheck, & 56 & xc_vdw_read, & 57 & xc_vdw_set_functional, & 58 & xc_vdw_show, & 59 & xc_vdw_status, & 60 & xc_vdw_trigger, & 61 & xc_vdw_write 62 63 private
ABINIT/m_xc_vdw/vdw_df_create_mesh [ Functions ]
NAME
vdw_df_create_mesh
FUNCTION
Creates a 1D mesh following user specifications.
INPUTS
npts= length of the mesh mesh_type= integer representing the type of the mesh mesh_cutoff= where to cut the mesh mesh_inc(optional)= mesh increment, if linear avoid_zero(optional)= whether to avoid zero in the mesh exact_max(optional)= whether to force the last mesh element to be strictly equal to the cutoff mesh_ratio(optional)= ratio between the first and last segment (or inverse), when applies mesh_tolerance(optional)= tolerance to apply for the mesh mesh_file(optional)= formatted file (e20.8) where to read the mesh, when applies
OUTPUT
mesh= the desired mesh
SOURCE
2860 subroutine vdw_df_create_mesh(mesh,npts,mesh_type,mesh_cutoff, & 2861 & mesh_inc,mesh_ratio,mesh_tolerance,mesh_file,avoid_zero,exact_max) 2862 2863 !Arguments ------------------------------------ 2864 integer,intent(in) :: npts,mesh_type 2865 real(dp),intent(in) :: mesh_cutoff 2866 real(dp),intent(out) :: mesh(npts) 2867 logical,optional,intent(in) :: avoid_zero,exact_max 2868 real(dp),optional,intent(in) :: mesh_inc,mesh_ratio,mesh_tolerance 2869 character(len=*),optional,intent(in) :: mesh_file 2870 2871 !Local variables ------------------------------ 2872 logical :: my_avoid_zero,my_exact_max 2873 integer :: im,unt 2874 real(dp) :: exp_inc,exp_ofs,exp_tol,lin_inc,log_inc,log_ofs,log_tol 2875 character(len=500) :: msg 2876 2877 ! ************************************************************************* 2878 2879 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2880 DBG_ENTER("COLL") 2881 #endif 2882 2883 if ( present(avoid_zero) ) then 2884 my_avoid_zero = avoid_zero 2885 else 2886 my_avoid_zero = .false. 2887 end if 2888 2889 if ( present(exact_max) ) then 2890 my_exact_max = exact_max 2891 else 2892 my_exact_max = .false. 2893 end if 2894 2895 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2896 write(msg,'(1x,a,a,1x,a,1x,i6,a,1x,a,i6,a,1x,a,1x,e10.3)') & 2897 & "> Mesh parameters",ch10, & 2898 & "> npts :",npts,ch10, & 2899 & "> type :",mesh_type,ch10, & 2900 & "> cut-off :",mesh_cutoff 2901 call wrtout(std_out,msg,'COLL') 2902 #endif 2903 2904 select case(mesh_type) 2905 2906 case(0) 2907 2908 if ( present(mesh_file) ) then 2909 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2910 write(msg,'(1x,a,1x,a)') "> file :",mesh_file 2911 call wrtout(std_out,msg,'COLL') 2912 #endif 2913 if (open_file(mesh_file,msg, newunit=unt,status="old",action="read") /= 0) then 2914 ABI_ERROR(msg) 2915 end if 2916 read(unt,'(e20.8)') mesh 2917 close(unit=unt) 2918 else 2919 ABI_ERROR(" You must specify a mesh file for this mesh type!") 2920 end if 2921 2922 case(1) 2923 2924 if ( present(mesh_ratio) ) then 2925 if ( present(mesh_tolerance) ) then 2926 log_tol = mesh_tolerance 2927 else 2928 log_tol = sqrt(my_vdw_params%tolerance) 2929 end if 2930 !DEBUG 2931 ! log_inc = mesh_ratio / (npts - 2) 2932 log_inc = log(mesh_ratio) / (npts - 1) 2933 !ENDDEBUG 2934 log_ofs = mesh_cutoff / (exp(log_inc*(npts - 1)) - 1) 2935 2936 !!DEBUG 2937 !! Previous definition of log_ofs implies that the mesh 2938 !! points are to be calculated as: 2939 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 ) 2940 !! which in turn means that x(0)=0 and is not consistent with the 2941 !! mesh(im) points given at line 3056. 2942 !! If the latter is adopted mesh_ratio above should be computed as: 2943 !! mesh_ratio=log ( (x(n)-x(n-1))/(x(2)-x(1)) ) and not merely as 2944 !! mesh_ratio= (x(n)-x(n-1))/(x(2)-x(1)) 2945 !! in order to log_inc definition makes sense. 2946 !! The latter relation between log_inc and mesh_ratio is valid even if 2947 !! mesh points do not start at 0: 2948 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 ) + x_0 2949 !!ENDDEBUG 2950 2951 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2952 write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') & 2953 & "> ratio :",mesh_ratio,ch10, & 2954 & "> log_inc :",log_inc,ch10, & 2955 & "> log_ofs :",log_ofs,ch10, & 2956 & "> log_tol :",log_tol 2957 call wrtout(std_out,msg,'COLL') 2958 #endif 2959 if ( log_inc < log_tol ) then 2960 ABI_WARNING(" The logarithmic increment is very low!") 2961 end if 2962 2963 do im = 1,npts 2964 !DEBUG 2965 ! mesh(im) = log_ofs * exp(log_inc * (im-1)) this is not consistent 2966 ! with definition of log_ofs at line 3025. 2967 mesh(im) = log_ofs * (exp(log_inc * (im-1)) - 1) 2968 !ENDDEBUG 2969 end do 2970 2971 !FIXME: check the correctness of 0.2 2972 if ( my_avoid_zero ) then 2973 mesh(1) = exp(log_inc * 0.2) * mesh_cutoff / & 2974 & (exp(log_inc*(npts - 1)) - one) 2975 end if 2976 2977 if ( my_exact_max ) then 2978 mesh(npts) = mesh_cutoff 2979 end if 2980 else 2981 ABI_ERROR(" You must specify a mesh ratio for this mesh type!") 2982 end if 2983 2984 ! Make sure the last point is qcut 2985 mesh(npts) = mesh_cutoff 2986 2987 case(2) 2988 2989 if ( present(mesh_inc) ) then 2990 lin_inc = mesh_inc 2991 else 2992 lin_inc = mesh_cutoff / npts 2993 end if 2994 2995 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2996 write(msg,'(1x,a,1x,e10.3)') "> lin_inc :",lin_inc 2997 call wrtout(std_out,msg,'COLL') 2998 #endif 2999 do im = 1,npts 3000 mesh(im) = lin_inc * (im - 1) 3001 end do 3002 3003 case(3) 3004 3005 if ( present(mesh_ratio) ) then 3006 if ( present(mesh_tolerance) ) then 3007 exp_tol = mesh_tolerance 3008 else 3009 exp_tol = sqrt(my_vdw_params%tolerance) 3010 end if 3011 exp_inc = mesh_ratio / (npts - 2) 3012 exp_ofs = mesh_cutoff * (log(exp_inc*(npts - 1)) + 1) 3013 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 3014 write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') & 3015 & "> ratio :",mesh_ratio,ch10, & 3016 & "> exp_inc :",exp_inc,ch10, & 3017 & "> exp_ofs :",exp_ofs,ch10, & 3018 & "> exp_tol :",exp_tol 3019 call wrtout(std_out,msg,'COLL') 3020 #endif 3021 if ( log_inc < log_tol ) then 3022 ABI_WARNING(" The logarithmic increment is very low!") 3023 end if 3024 3025 do im = 1,npts 3026 mesh(im) = exp_ofs / log(exp_inc * (im-1)) 3027 end do 3028 else 3029 ABI_ERROR(" You must specify a mesh ratio for this mesh type!") 3030 end if 3031 3032 case default 3033 3034 write(msg,'(2x,a,1x,i2)') "Unknown mesh type",mesh_type 3035 ABI_ERROR(msg) 3036 3037 end select 3038 3039 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 3040 DBG_EXIT("COLL") 3041 #endif 3042 3043 end subroutine vdw_df_create_mesh
ABINIT/m_xc_vdw/vdw_df_filter [ Functions ]
NAME
vdw_df_filter
FUNCTION
Softens the kernel by applying a filter in reciprocal space.
INPUTS
nqpts= number of q-mesh points nrpts= number of mesh points in real space rcut= cut-off in real space gcut= cut-off in reciprocal space sofswt= Driver of the softening in real space OUTPUTS ngpts= number of mesh points in reciprocal space
SIDE EFFECTS
phir= the softened kernel in real space phir_u= the unsoftened kernel in real space phig= the softened kernel in reciprocal space
SOURCE
2200 subroutine vdw_df_filter(nqpts,nrpts,rcut,gcut,ngpts,sofswt) 2201 2202 !Arguments ------------------------------------ 2203 integer,intent(in) :: nqpts,nrpts,sofswt 2204 integer,intent(out) :: ngpts 2205 real(dp),intent(in) :: rcut 2206 real(dp),intent(inout) :: gcut 2207 !Local variables ------------------------------ 2208 integer :: ig,iq1,iq2,ir,rfftt,id1,ndpts,ng 2209 real(dp) :: dg,dr,lstep,ptmp,q1,q2,qcut,qtol,un,xr 2210 real(dp) :: x1,x2,yq1,yqn,yr1,yrn,gmax 2211 real(dp),allocatable :: gmesh(:),rmesh(:),utmp(:),phiraux(:,:,:) 2212 real(dp),allocatable :: q1sat(:),q2sat(:) 2213 ! ************************************************************************* 2214 2215 DBG_ENTER("COLL") 2216 2217 ! Init 2218 qtol = my_vdw_params%tolerance 2219 qcut = my_vdw_params%qcut 2220 write(std_out,*) 'gcut before =', gcut 2221 dr = rcut / nrpts 2222 ngpts = nrpts 2223 dg = pi / rcut 2224 gcut = pi / dr !kmax in siesta 2225 gmax = 10.639734883681651 !kcut in siesta for a particular 2226 ! system (Ar). This should be computed each time in terms of 2227 ! the FFT and BZ parameters, not hard wired like here. 2228 ng = 1 + int(gmax/dg) ! nk in siesta 2229 rfftt = 2 !Type of radial sin transform 2230 write(std_out,*) 'gmax from siesta value for kmax=', gmax 2231 write(std_out,*) 'ng from siesta gmax=', ng 2232 write(std_out,*) 'gcut after =', gcut 2233 write(std_out,*) 'ngpts after =', ngpts 2234 2235 ABI_MALLOC(utmp,(0:nrpts)) 2236 2237 ! Create radial mesh for radial FT: 2238 ! shared/common/src/32_util/radsintr.F90 2239 ABI_MALLOC(rmesh,(0:nrpts)) 2240 forall(ir=1:nrpts) rmesh(ir) = dr * dble(ir) 2241 2242 if ( sofswt == 1 ) then 2243 ! Create reciprocal radial mesh 2244 ABI_MALLOC(gmesh,(0:ngpts)) 2245 !----my debug------------------------------- 2246 ABI_MALLOC(phiraux,(0:nrpts,nqpts,nqpts)) 2247 !----end my debug--------------------------- 2248 forall(ig=0:ngpts) gmesh(ig) = dg * dble(ig) 2249 end if 2250 2251 ABI_MALLOC(q1sat,(2)) 2252 ABI_MALLOC(q2sat,(2)) 2253 2254 ! Build filtered kernel for each (q1,q2) pair 2255 do iq1=1,nqpts 2256 !Inverse saturation of q1-mesh value 2257 x1 = 0 2258 x2 = 1.5_dp * qcut 2259 do 2260 q1 = (x1+x2)/2 2261 call vdw_df_saturation( q1, qcut, q1sat ) 2262 if (abs(qmesh(iq1)-q1sat(1))<qtol) then 2263 exit 2264 else if (q1sat(1) < qmesh(iq1)) then 2265 x1 = q1 2266 else 2267 x2 = q1 2268 end if 2269 end do ! end calculation of unsaturated qmesh(iq1) 2270 2271 do iq2=1,iq1 2272 !Rather than this q2 = qmesh(iq2), the values of q-mesh 2273 !are unsaturated. 2274 x1 = 0 2275 x2 = 1.5_dp * qcut 2276 do 2277 q2 = (x1+x2)/2 2278 call vdw_df_saturation( q2, qcut, q2sat ) 2279 if (abs(qmesh(iq2)-q2sat(1))<qtol) then 2280 exit 2281 else if (q2sat(1) < qmesh(iq2)) then 2282 x1 = q2 2283 else 2284 x2 = q2 2285 end if 2286 end do ! end calculation of unsaturated qmesh(iq2) 2287 2288 if ( sofswt == 1 ) then 2289 ! Build kernel in real space 2290 ! Note: smoothly going to zero when approaching rcut 2291 ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values. 2292 do ir=1,nrpts 2293 xr=rmesh(ir) 2294 phir(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * & 2295 & (one - ( xr / rmesh(nrpts))**8)**4 2296 end do 2297 !-----debug------------------------------------- 2298 phiraux(:,iq1,iq2) = phir(:,iq1,iq2) 2299 !-----end debug--------------------------------- 2300 2301 ! Obtain kernel in reciprocal space 2302 call radsintr(phir(:,iq1,iq2),phig(:,iq1,iq2),ngpts,nrpts,gmesh,rmesh,yq1,yqn,rfftt) 2303 if (rfftt == 2) then 2304 phig(:,iq1,iq2) = phig(:,iq1,iq2) * (2*pi)**1.5_dp 2305 end if 2306 2307 ! Filter in reciprocal space 2308 phig(ng:ngpts,iq1,iq2) = zero 2309 do ig=1,ng 2310 phig(ig,iq1,iq2) = phig(ig,iq1,iq2) * & 2311 & (one - ( gmesh(ig) / gmesh(ng))**8)**4 2312 end do 2313 2314 ! Go back to real space 2315 call radsintr(phig(:,iq1,iq2),phir(:,iq1,iq2),nrpts,ngpts,rmesh,gmesh,yr1,yrn,rfftt) 2316 phir(:,iq1,iq2) = phir(:,iq1,iq2) / (2*pi)**1.5_dp 2317 2318 ! Calculate second derivative in real space 2319 d2phidr2(0,iq1,iq2) = -half 2320 utmp(0) = (three / dr) * (phir(1,iq1,iq2)-phir(0,iq1,iq2)) / dr 2321 do ir=1,nrpts-1 2322 ptmp = half * d2phidr2(ir-1,iq1,iq2) + two 2323 d2phidr2(ir,iq1,iq2) = (half - one) / ptmp 2324 utmp(ir) = (three * (phir(ir+1,iq1,iq2) + phir(ir-1,iq1,iq2) - & 2325 & two*phir(ir,iq1,iq2)) / (dr**2) - half * utmp(ir-1)) / ptmp 2326 end do 2327 un = (three / dr) * (phir(nrpts-1,iq1,iq2)-phir(nrpts,iq1,iq2)) / dr 2328 d2phidr2(nrpts,iq1,iq2) = (un - half * utmp(nrpts-1)) / & 2329 (half * d2phidr2(nrpts-1,iq1,iq2) + one) 2330 do ir=nrpts-1,0,-1 2331 d2phidr2(ir,iq1,iq2) = d2phidr2(ir,iq1,iq2) * & 2332 & d2phidr2(ir+1,iq1,iq2) + utmp(ir) 2333 end do 2334 2335 ! Calculate second derivative in reciprocal space 2336 d2phidg2(0,iq1,iq2) = -half 2337 utmp(0) = (three / dg) * (phig(1,iq1,iq2)-phig(0,iq1,iq2)) / dg 2338 do ig=1,ngpts-1 2339 ptmp = half * d2phidg2(ig-1,iq1,iq2) + two 2340 d2phidg2(ig,iq1,iq2) = (half - one) / ptmp 2341 utmp(ig) = (three * (phig(ig+1,iq1,iq2) + phig(ig-1,iq1,iq2) - & 2342 & two*phig(ig,iq1,iq2)) / (dg**2) - half * utmp(ig-1)) / ptmp 2343 end do 2344 un = (three / dg) * (phig(ngpts-1,iq1,iq2)-phig(ngpts,iq1,iq2)) / dg 2345 d2phidg2(ngpts,iq1,iq2) = (un - half * utmp(ngpts-1)) / & 2346 (half * d2phidg2(ngpts-1,iq1,iq2) + one) 2347 do ig=ngpts-1,0,-1 2348 d2phidg2(ig,iq1,iq2) = d2phidg2(ig,iq1,iq2) * & 2349 & d2phidg2(ig+1,iq1,iq2) + utmp(ig) 2350 end do 2351 2352 ! Symmetrize kernels & derivatives 2353 phir(:,iq2,iq1) = phir(:,iq1,iq2) 2354 phig(:,iq2,iq1) = phig(:,iq1,iq2) 2355 d2phidr2(:,iq2,iq1) = d2phidr2(:,iq1,iq2) 2356 d2phidg2(:,iq2,iq1) = d2phidg2(:,iq1,iq2) 2357 !---------my debug---------------------- 2358 phiraux(:,iq2,iq1) = phiraux(:,iq1,iq2) 2359 !----end my debug---------------- 2360 else! sofswt == 1 2361 2362 ! Build unsoftened kernel in real space 2363 ! Note: smoothly going to zero when approaching rcut 2364 ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values. 2365 do ir=0,nrpts 2366 xr=rmesh(ir) 2367 phir_u(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * & 2368 & (one - ( xr / rmesh(nrpts))**8)**4 2369 end do 2370 2371 ! Symmetrize unsoftened kernel 2372 phir_u(:,iq2,iq1) = phir_u(:,iq1,iq2) 2373 end if ! sofswt == 1 2374 2375 end do !loop on iq2 2376 end do !loop on iq1 2377 2378 !---------my debug---------------------- 2379 if ( sofswt == 1 ) then 2380 ! open(unit=68,file='phir-gfil-RADSINTR-A.dat') 2381 ! open(unit=67,file='d2phidr2-gfil-RADSINTR-A.dat') 2382 ! do ir=0,nrpts 2383 ! write(68,*) dble(ir)*dr, ((phiraux(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2384 ! write(67,*) dble(ir)*dr, ((d2phidr2(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2385 ! end do 2386 ! close(unit=67) 2387 ! close(unit=68) 2388 2389 ! open(unit=78,file='phig-gfil-RADSINTR-A.dat') 2390 ! open(unit=77,file='d2phidg2-gfil-RADSINTR-A.dat') 2391 ! do ir=0,nrpts 2392 ! write(78,*) dble(ir)*dg, ((phig(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2393 ! write(77,*) dble(ir)*dg, ((d2phidg2(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2394 ! end do 2395 ! close(unit=77) 2396 ! close(unit=78) 2397 2398 ! open(unit=87,file='phirfil-gfil-RADSINTR-A.dat') 2399 ! do ir=0,nrpts 2400 ! write(87,*) dble(ir)*dr, ((phir(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2401 ! end do 2402 ! close(unit=87) 2403 2404 ABI_FREE(phiraux) 2405 2406 end if 2407 2408 ! if (sofswt == 0 ) then 2409 ! open(unit=88,file='phir-uns-gfil-RADSINTR-A.dat') 2410 ! do ir=0,nrpts 2411 ! write(88,*) dble(ir)*dr, ((phir_u(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 2412 ! end do 2413 ! close(unit=88) 2414 ! end if 2415 !----end my debug---------------- 2416 2417 ABI_FREE(utmp) 2418 ABI_FREE(q1sat) 2419 ABI_FREE(q2sat) 2420 2421 DBG_EXIT("COLL") 2422 2423 end subroutine vdw_df_filter
ABINIT/m_xc_vdw/vdw_df_indexof [ Functions ]
NAME
vdw_df_indexof
FUNCTION
Returns the index of a specified value in a sorted array.
INPUTS
list_1d= sorted array npts= length of the sorted array value= value the index of which has to be found
NOTES
The returned index is such that list_1d(index) < value < list_1d(index+1)
SOURCE
3063 function vdw_df_indexof(list_1d,npts,value) 3064 3065 !Arguments ------------------------------------ 3066 integer,intent(in) :: npts 3067 real(dp),intent(in) :: value,list_1d(npts) 3068 3069 !Local variables ------------------------------ 3070 integer :: imin,imax,imid 3071 integer :: vdw_df_indexof 3072 3073 ! ************************************************************************* 3074 3075 imin = 1 3076 imax = npts 3077 3078 do 3079 imid = (imin + imax) / 2 3080 if ( value < list_1d(imid) ) then 3081 imax = imid 3082 else 3083 imin = imid 3084 end if 3085 if ( imin + 1 >= imax ) exit 3086 end do 3087 3088 vdw_df_indexof = imin 3089 3090 end function vdw_df_indexof
ABINIT/m_xc_vdw/vdw_df_internal_checks [ Functions ]
NAME
vdw_df_internal_checks
FUNCTION
Performs consistency checks of the vdW-DF kernel.
INPUTS
test_mode= bitfield to enable/disable tests
SOURCE
3181 subroutine vdw_df_internal_checks(test_mode) 3182 3183 #if defined HAVE_FFTW3 3184 include "fftw3.f03" 3185 #endif 3186 3187 !Arguments ------------------------------------ 3188 integer,intent(in) :: test_mode 3189 3190 !Local variables ------------------------------ 3191 character(len=512) :: msg 3192 integer :: ndpts,ngpts,nqpts,nrpts,ntest,id1,ig,iq1,iq2,jj,ir,nn,ii 3193 integer :: iostatus,ip1,ir1,ir2,ir3,nr1,nr2,nr3,i1,i2,i3 3194 integer(kind=SIZEOF_PTRDIFF_T) :: plan 3195 real(dp) :: acutmin,aratio,damax,dsoft,phisoft 3196 real(dp) :: a1,a2,a3,delta_test,d1,d2i,gcut,gtmp,grho2,grhotmp2(3) 3197 real(dp) :: b1,b2,b3,dg,dr,rcut,rtmp,dvol,exc_tmp,deltae_vdw 3198 real(dp),allocatable :: gprimdt(:,:),kern_test_c(:),kern_test_i(:) 3199 real(dp),allocatable :: intg_calc(:),intg_int(:),ptmp(:,:,:) 3200 real(dp),allocatable :: extin(:,:), grhotmp(:,:),ztmp(:,:) 3201 real(dp),allocatable :: decdrho_tmp(:),decdgrho_tmp(:,:) 3202 real(dp),allocatable :: t3dr(:,:,:,:),t3dr1d(:,:),gvecs(:,:),t3dr2(:,:,:,:) 3203 real(dp),allocatable :: u1dtmp(:,:) 3204 real(dp) :: theta(my_vdw_params%nqpts,1,5),dummyv 3205 complex(dp),allocatable :: t3dg(:,:,:,:),t3dg1d(:,:),utmp(:) 3206 ! real(dp),allocatable :: t3dg(:,:,:,:),t3dg1d(:,:),utmp(:),u1dtmp(:,:) 3207 ! ************************************************************************* 3208 3209 DBG_ENTER("COLL") 3210 3211 ! Here starts a test which will evaluate the integral over D for 3212 ! delta=0 (see Fig. 1 of DRSLL04). 3213 ! The integration will be performed over a finer and linear d mesh, 3214 ! not dmesh, for which the kernel values will be obtained by both 3215 ! interpolation and direct calculation. The range spanned by this mesh 3216 ! is half of dmesh range. Introduced variables: delta_test, kern_test_c, 3217 ! kern_test_i, ntest, intg_calc, intg_int. 3218 if ( iand(test_mode,1) /= 0 ) then 3219 3220 write(msg,'(1x,a)') "[vdW-DF Check] Test #01: Kernel integral over D" 3221 call wrtout(std_out,msg,'COLL') 3222 3223 acutmin = my_vdw_params%acutmin 3224 aratio = my_vdw_params%aratio 3225 damax = my_vdw_params%damax 3226 dsoft = my_vdw_params%dsoft 3227 ndpts = my_vdw_params%ndpts 3228 phisoft = my_vdw_params%phisoft 3229 3230 ntest = 600 3231 delta_test = dmesh(ndpts) / (2 * ntest) 3232 3233 ABI_MALLOC(kern_test_c,(ntest)) 3234 ABI_MALLOC(kern_test_i,(ntest)) 3235 ABI_MALLOC(intg_calc,(ntest)) 3236 ABI_MALLOC(intg_int,(ntest)) 3237 3238 #if defined DEBUG_VERBOSE 3239 write(msg,'(a,4x,3(1x,a10),a)') ch10, "D","Calcul.","Interp.",ch10 3240 call wrtout(std_out,msg,'COLL') 3241 #endif 3242 3243 do id1 = 1,ntest 3244 d1=id1*delta_test 3245 kern_test_c(id1) = (four_pi * d1**2) * & 3246 & vdw_df_kernel_value(d1,d1,acutmin,aratio,damax) 3247 kern_test_i(id1) = (four_pi * d1**2) * vdw_df_interpolate(d1,d1,0) 3248 #if defined DEBUG_VERBOSE 3249 write(msg,'(4x,3(1x,e10.3))') id1*delta_test, & 3250 & kern_test_c(id1), kern_test_i(id1) 3251 call wrtout(std_out,msg,'COLL') 3252 #endif 3253 end do 3254 3255 call simpson_int(ntest,delta_test,kern_test_c,intg_calc) 3256 call simpson_int(ntest,delta_test,kern_test_i,intg_int) 3257 3258 write(msg,'(a,4x,a,2(a,4x,a,e10.3))') ch10, & 3259 & "> Integrals over D", ch10, "> calculated :",intg_calc(ntest),ch10, & 3260 & "> interpolated :",intg_int(ntest) 3261 call wrtout(std_out,msg,'COLL') 3262 if ( abs(intg_int(ntest) - intg_calc(ntest)) < tol3 ) then 3263 write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 PASSED" 3264 else 3265 write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 FAILED" 3266 end if 3267 call wrtout(std_out,msg,'COLL') 3268 3269 ABI_FREE(kern_test_c) 3270 ABI_FREE(kern_test_i) 3271 ABI_FREE(intg_calc) 3272 ABI_FREE(intg_int) 3273 3274 ! The following test aims at checking the interpolation in both 3275 ! real and reciprocal space of the radial representation of the 3276 ! Kernel 3277 3278 write(msg,'(1x,a)') "[vdW-DF Check] Test #02: Interpolation & 3279 & of radial reciprocal space representation of the Kernel" 3280 call wrtout(std_out,msg,'COLL') 3281 3282 ! Interpolate phi in reciprocal space: 3283 ! * ptmp(:,:,1) = phi(g) 3284 ! * ptmp(:,:,2) = dphi(g)/dg 3285 ! The interolation will be performed in a three times 3286 ! finer mesh than that used for the spline, and for 3287 ! iq1=1q2=1,2 only. 3288 3289 rcut = my_vdw_params%rcut 3290 nrpts = my_vdw_params%nrpts 3291 gcut = my_vdw_params%gcut 3292 ngpts = my_vdw_params%ngpts 3293 ! begin my debug 3294 write(std_out,*) 'ngpts after init: ',ngpts 3295 ! end my debug 3296 ABI_MALLOC(ptmp,(2,2,2)) 3297 3298 dg = pi / rcut !gcut / (ngpts-one) 3299 #if defined DEBUG_VERBOSE 3300 write(msg,'(a,5x,5(1x,a10),a)') ch10, "g","ig","phi(g,1,1)","phi(g,1,2)" & 3301 & ,"phi(g,2,2)", ch10 3302 call wrtout(std_out,msg,'COLL') 3303 #endif 3304 do jj = 0,3*ngpts-1 !3*(ngpts-1) 3305 ptmp(:,:,:) = zero 3306 gtmp = jj * dg / three 3307 ig = int(gtmp / dg) + 1 3308 a1 = ig - gtmp / dg !((ig + 1) * dg - gtmp) / dg 3309 b1 = one - a1 3310 a2 = (3 * a1**2 - one) * dg / six 3311 b2 = (3 * b1**2 - one) * dg / six 3312 a3 = (a1**3 - a1) * dg**2 / six 3313 b3 = (b1**3 - b1) * dg**2 / six 3314 do iq2 = 1,2 3315 do iq1 = 1,iq2 3316 ! ptmp(iq1,iq2,1) = a1 * phig(ig,iq1,iq2) + b1 * phig(ig+1,iq1,iq2) & 3317 !& + a3 * d2phidg2(ig,iq1,iq2) + b3 * d2phidg2(ig+1,iq1,iq2) 3318 ! ptmp(iq1,iq2,2) = (phig(ig+1,iq1,iq2) - phig(ig,iq1,iq2)) / dg & 3319 !& - a2 * d2phidg2(ig,iq1,iq2) + b2 * d2phidg2(ig+1,iq1,iq2) 3320 ptmp(iq1,iq2,1) = a1 * phig(ig-1,iq1,iq2) + b1 * phig(ig,iq1,iq2) & 3321 & + a3 * d2phidg2(ig-1,iq1,iq2) + b3 * d2phidg2(ig,iq1,iq2) 3322 ptmp(iq1,iq2,2) = (phig(ig,iq1,iq2) - phig(ig-1,iq1,iq2)) / dg & 3323 & - a2 * d2phidg2(ig-1,iq1,iq2) + b2 * d2phidg2(ig,iq1,iq2) 3324 end do 3325 end do 3326 3327 do iq2 = 1,2 3328 do iq1 = 1,iq2-1 3329 ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:) 3330 end do 3331 end do 3332 #if defined DEBUG_VERBOSE 3333 write(msg,'(5x,5(1x,f16.8))') gtmp, real(ig), & 3334 & ptmp(1,1,1), ptmp(1,2,1), ptmp(2,2,1) 3335 call wrtout(std_out,msg,'COLL') 3336 #endif 3337 end do 3338 3339 write(msg,'(1x,a)') "[vdW-DF Check] Test #02: Interpolation & 3340 & of radial real space representation of the Kernel" 3341 call wrtout(std_out,msg,'COLL') 3342 3343 ! Interpolate phi in real space: 3344 ! * ptmp(:,:,1) = phi(r) 3345 ! * ptmp(:,:,2) = dphi(r)/dr 3346 ! The interolation will be performed in a ten times 3347 ! finer mesh than that used for the spline and for 3348 ! iq1=1q2=1,2. 3349 3350 dr = rcut / nrpts !(nrpts-one) 3351 #if defined DEBUG_VERBOSE 3352 write(msg,'(a,4x,4(1x,a10),a)') ch10, "r","phi(r,1,1)","phi(r,1,2)" & 3353 & ,"phi(r,2,2)", ch10 3354 call wrtout(std_out,msg,'COLL') 3355 #endif 3356 do jj = 0,3*nrpts-1 !-1 3357 ptmp(:,:,:) = zero 3358 rtmp = jj * dr / three 3359 ir = int(rtmp / dr) + 1 3360 a1 = ir - ( rtmp / dr ) 3361 b1 = one - a1 3362 a2 = (3 * a1**2 - one) * dr / six 3363 b2 = (3 * b1**2 - one) * dr / six 3364 a3 = (a1**3 - a1) * dr**2 / six 3365 b3 = (b1**3 - b1) * dr**2 / six 3366 do iq2 = 1,2 3367 do iq1 = 1,iq2 3368 ! ptmp(iq1,iq2,1) = a1 * phir(ir,iq1,iq2) + b1 * phir(ir+1,iq1,iq2) & 3369 !& + a3 * d2phidr2(ir,iq1,iq2) + b3 * d2phidr2(ir+1,iq1,iq2) 3370 ! ptmp(iq1,iq2,2) = (phir(ir+1,iq1,iq2) - phir(ir,iq1,iq2)) / dr & 3371 !& - a2 * d2phidr2(ir,iq1,iq2) + b2 * d2phidr2(ir+1,iq1,iq2) 3372 ptmp(iq1,iq2,1) = a1 * phir(ir-1,iq1,iq2) + b1 * phir(ir,iq1,iq2) & 3373 & + a3 * d2phidr2(ir-1,iq1,iq2) + b3 * d2phidr2(ir,iq1,iq2) 3374 ptmp(iq1,iq2,2) = (phir(ir,iq1,iq2) - phir(ir-1,iq1,iq2)) / dr & 3375 & - a2 * d2phidr2(ir-1,iq1,iq2) + b2 * d2phidr2(ir,iq1,iq2) 3376 end do 3377 end do 3378 3379 do iq2 = 1,2 3380 do iq1 = 1,iq2-1 3381 ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:) 3382 end do 3383 end do 3384 #if defined DEBUG_VERBOSE 3385 write(msg,'(4x,4(1x,f16.8))') rtmp, & 3386 & ptmp(1,1,1), ptmp(1,2,1), ptmp(2,2,1) 3387 call wrtout(std_out,msg,'COLL') 3388 #endif 3389 end do 3390 3391 ABI_FREE(ptmp) 3392 3393 !-------my debug------------------------------------------------ 3394 ! The following test aims at checking the calculation of q0, 3395 ! the interpolation polynomials at q0, and the local cusp 3396 ! correction from precomputed electronic density and its 3397 ! gradient. Also theta and the nonlocal correction is 3398 ! computed. 3399 3400 write(msg,'(1x,a)') "[vdW-DF Check] Test #03: q0 and & 3401 & interpolation polynomials from external rho and grho" 3402 call wrtout(std_out,msg,'COLL') 3403 3404 nqpts = my_vdw_params%nqpts 3405 3406 open (unit=98, file='rhogrhoext.dat', status='old',action='read') 3407 open (unit=59, file='rho-grho2-q0-qpoly-check.dat',status='replace') 3408 open (unit=55, file='rho-epscusp-check.dat',status='replace') 3409 open (unit=60, file='rho-grho2-theta-check.dat',status='replace') 3410 open (unit=61, file='rho-grho2-dtdd-check.dat',status='replace') 3411 ! open (unit=62, file='rho-grho2-dtdgd-check.dat',status='replace') 3412 3413 nn=0 3414 do 3415 read(98,*, iostat=iostatus) dummyv 3416 if(iostatus/=0) then ! to avoid end of file error. 3417 exit 3418 else 3419 nn=nn+1 3420 end if 3421 end do 3422 write(std_out,*) 'Number of lines in rhogrhoext.dat:', nn 3423 3424 ABI_MALLOC(extin,(nn,9)) 3425 ABI_MALLOC(grhotmp,(1,3)) 3426 ABI_MALLOC(ztmp,(nqpts,nqpts)) 3427 ABI_MALLOC(decdrho_tmp,(1)) 3428 ABI_MALLOC(decdgrho_tmp,(3,1)) 3429 ABI_MALLOC(gprimdt,(3,3)) 3430 3431 write(std_out,*) 'Allocated arrays...' 3432 3433 ! extin contains rho, grho(1), grho(2), grho(3), ex, ec, vx, vc, q0 3434 rewind(98) 3435 3436 write(std_out,*) 'Rewinded file' 3437 3438 do ii=1,nn 3439 read(98,*) (extin(ii,jj),jj=1,9) 3440 end do 3441 3442 write(std_out,*) 'Readed external file...' 3443 3444 ! Pre-compute integrand for vdW energy correction 3445 ztmp(:,:) = zero 3446 do ir=1,nrpts 3447 do iq1=1,nqpts 3448 do iq2=1,iq1 3449 ztmp(iq1,iq2) = ztmp(iq1,iq2) - & 3450 & two * pi * dr * phir(ir,iq1,iq2) * (ir*dr)**2 3451 !& two * pi * ( phir_u(ir,iq1,iq2) - phir(ir,iq1,iq2) ) * dr 3452 end do 3453 end do 3454 end do 3455 do iq1=2,nqpts 3456 do iq2=1,iq1-1 3457 ztmp(iq2,iq1) = ztmp(iq1,iq2) 3458 end do 3459 end do 3460 !---------my debug---------------------------------- 3461 open(34,file='table-localcorr.dat',status='replace') 3462 do iq2=1,nqpts 3463 do iq1=1,nqpts 3464 write(34,*) iq1, iq2, ztmp(iq1,iq2) 3465 enddo 3466 enddo 3467 close(34) 3468 write(std_out,*) 'Table file written...' 3469 3470 !--------end my debug------------------------------ 3471 3472 ! Compute q0, polynomials, theta from Abinit routines 3473 ! Also get nonlocal correction to the energy 3474 open(unit=81,file='gvecs.dat',status='replace') 3475 theta(:,:,:) = zero 3476 deltae_vdw = zero 3477 dvol = 0.025742973 ! this is a particular value for the 3478 ! precomputed density. Same for gprimd 3479 3480 gprimdt(1,1)= 0.22166114341982476 3481 gprimdt(2,1)= 0.0d0 3482 gprimdt(3,1)= 0.0d0 3483 gprimdt(1,2)= 0.0d0 3484 gprimdt(2,2)= 0.33249171512973713 3485 gprimdt(3,2)= 0.0d0 3486 gprimdt(1,3)= 0.0d0 3487 gprimdt(2,3)= 0.0d0 3488 gprimdt(3,3)= 0.33249171512973713 3489 3490 nr1=96 !particular values for this precomputed density 3491 nr2=64 3492 nr3=64 3493 3494 ABI_MALLOC(t3dr,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts)) 3495 ABI_MALLOC(t3dr2,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts)) 3496 ABI_MALLOC(t3dg,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts)) 3497 ABI_MALLOC(t3dr1d,(nqpts,nr1*nr2*nr3)) 3498 ABI_MALLOC(t3dg1d,(nqpts,nr1*nr2*nr3)) 3499 ABI_MALLOC(gvecs,(3,nr1*nr2*nr3)) 3500 ABI_MALLOC(utmp,(nqpts)) 3501 ABI_MALLOC(u1dtmp,(nqpts,nr1*nr2*nr3)) 3502 3503 ip1 = 1 3504 do ir3=0,nr3-1 !1,nr3 3505 do ir2=0,nr2-1 !1,nr2 3506 do ir1=0,nr1-1 !1,nr1 3507 decdrho_tmp(:) = zero 3508 decdgrho_tmp(:,:) = zero 3509 grhotmp(1,1:3) = extin(ip1,2:4) 3510 forall(ii=1:3) grhotmp2(ii) = grhotmp(1,ii) 3511 grho2 = sum(grhotmp2**2) 3512 !write(*,*) 'Call to xc_vdw_energy number:', ip1 3513 call xc_vdw_energy(1,extin(ip1,1),grhotmp,extin(ip1,5),& 3514 & extin(ip1,6),extin(ip1,7),extin(ip1,8),theta,& 3515 & exc_tmp,decdrho_tmp,decdgrho_tmp,ztmp) 3516 t3dr(ir1,ir2,ir3,1:nqpts) = theta(1:nqpts,1,1) 3517 deltae_vdw = deltae_vdw + extin(ip1,1) * exc_tmp * dvol 3518 ! forall(ii=1:3) grhotmp2(ii) = grhotmp(1,ii) 3519 ! grho2 = sum(grhotmp2**2) 3520 ! exc_vdw = exc_vdw + rho_tmp * & 3521 !& ( half * ( sum(ttmp2(:,ip1) * theta(:,1,1)) / & 3522 !& (rho_tmp + tiny(rho_tmp)) ) * dvol ) 3523 if (ir1 > nr1/2) then 3524 i1 = ir1 - nr1 3525 else 3526 i1 = ir1 3527 end if 3528 if (ir2 > nr2/2) then 3529 i2 = ir2 - nr2 3530 else 3531 i2 = ir2 3532 end if 3533 if (ir3 > nr3/2) then 3534 i3 = ir3 - nr3 3535 else 3536 i3 = ir3 3537 end if 3538 gvecs(:,ip1) = gprimdt(:,1) * i1 + gprimdt(:,2)& 3539 & * i2 + gprimdt(:,3) * i3 3540 write(81,*) ip1, gvecs(:,ip1) 3541 write(55,*) extin(ip1,1), exc_tmp 3542 ! write(60,*) extin(ip1,1), grho2, (theta(iq1,1,1),iq1=1,nqpts) 3543 ! write(61,*) extin(ip1,1), grho2, (theta(iq1,1,2),iq1=1,nqpts) 3544 ! write(62,*) extin(ip1,1), grho2, & 3545 !& (theta(iq1,1,3),theta(iq1,1,4),theta(iq1,1,5),iq1=1,nqpts) 3546 ip1 = ip1 +1 3547 end do 3548 end do 3549 end do 3550 write(std_out,*) 'Cusp correction for this density:' 3551 write(std_out,*) 'deltae_vdw=',deltae_vdw 3552 close(unit=81) 3553 close(unit=98) 3554 close(unit=59) 3555 close(unit=55) 3556 close(unit=60) 3557 close(unit=61) 3558 ! close(unit=62) 3559 ! Evaluation of 3D FFT of a component of theta 3560 ! open (unit=63, file='thetar-3D-iq20.dat', status='old',action='read') 3561 ! open (unit=63, file='R-thetag-1D-check.dat', status='replace') 3562 ! open (unit=64, file='I-thetag-1D-check.dat', status='replace') 3563 open (unit=65, file='thetar-1D-check.dat', status='replace') 3564 ! open (unit=66, file='thetar-1D-after-check.dat', status='replace') 3565 ! nn=0 3566 ! do 3567 ! read(63,*, iostat=iostatus) dummyv 3568 ! if(iostatus/=0) then ! to avoid end of file error. 3569 ! exit 3570 ! else 3571 ! nn=nn+1 3572 ! end if 3573 ! end do 3574 ! rewind(63) 3575 ! write(*,*) 'Rewinded file' 3576 3577 ! write(*,*) 'Number of lines in thetar-3D-iq20.dat:', nn 3578 ! nr1=96 !particular values for this precomputed density 3579 ! nr2=64 3580 ! nr3=64 3581 ! ABI_MALLOC(t3dr,(nr1,nr2,nr3,1)) 3582 ! ABI_MALLOC(t3dg,(nr1,nr2,nr3,1)) 3583 ! ABI_MALLOC(t3dg1d,(nr1*nr2*nr3)) 3584 ! ABI_MALLOC(gvecs,(3,nr1*nr2*nr3)) 3585 ! do ir3=1,nr3 3586 ! do ir2=1,nr2 3587 ! do ir1=1,nr1 3588 ! read(63,*) t3dr(ir1,ir2,ir3,1) 3589 ! end do 3590 ! end do 3591 ! end do 3592 3593 ! Fourier-transform theta 3594 #if defined HAVE_FFTW3 3595 do iq1 = 1,nqpts 3596 call dfftw_plan_dft_r2c_3d(plan,nr1,nr2,nr3, & 3597 & t3dr(:,:,:,iq1),t3dg(:,:,:,iq1),FFTW_ESTIMATE) 3598 call dfftw_execute(plan) 3599 call dfftw_destroy_plan(plan) 3600 t3dg(:,:,:,iq1) = t3dg(:,:,:,iq1) / (nr1 * nr2 * nr3) 3601 end do 3602 ! write(*,*) 'FFTW3 call number', iq1 3603 !fftw_plan_r2r_3d 3604 #else 3605 ABI_ERROR('vdW-DF calculations require FFTW3 support') 3606 #endif 3607 ! Repack theta 3608 ip1 = 1 3609 do ir3=0,nr3-1 !1,nr3 3610 do ir2=0,nr2-1 !1,nr2 3611 do ir1=0,nr1-1 !1,nr1 3612 t3dg1d(1:nqpts,ip1) = t3dg(ir1,ir2,ir3,1:nqpts) 3613 t3dr1d(1:nqpts,ip1) = t3dr(ir1,ir2,ir3,1:nqpts) 3614 ! write(63,*) ip1, (real(t3dg1d(iq1,ip1)),iq1=1,nqpts) 3615 ! write(64,*) ip1, (aimag(t3dg1d(iq1,ip1)),iq1=1,nqpts) 3616 ! write(65,*) ip1, (t3dr1d(iq1,ip1),iq1=1,nqpts) 3617 ip1 = ip1 + 1 3618 end do 3619 end do 3620 end do 3621 write(std_out,*) 'Finished 1D repack of theta.' 3622 3623 ! Fourier-transform back theta just to check fftw3 3624 !#if defined HAVE_FFT_FFTW3 3625 ! write(*,*) "Backward FFT of theta - test" 3626 ! do iq1=1,nqpts 3627 ! call dfftw_plan_dft_c2r_3d(plan,nr1,nr2,nr3, & 3628 !& t3dg(:,:,:,iq1),t3dr2(:,:,:,iq1),FFTW_ESTIMATE) 3629 ! call dfftw_execute(plan) 3630 ! call dfftw_destroy_plan(plan) 3631 ! end do 3632 ! ip1 = 1 3633 ! do ir3=0,nr3-1 !1,nr3 3634 ! do ir2=0,nr2-1 !1,nr2 3635 ! do ir1=0,nr1-1 !1,nr1 3636 ! write(66,*) ip1, (t3dr2(ir1,ir2,ir3,iq1),iq1=1,nqpts) 3637 ! ip1 = ip1 + 1 3638 ! end do 3639 ! end do 3640 ! end do 3641 !#else 3642 ! MSG_ERROR('vdW-DF calculations require FFTW3 support') 3643 !#endif 3644 3645 open(unit=82,file='Interp-phi-check.dat',status='replace') 3646 ! Reset 3D counters, since theta will be reconstructed in 3D on the fly 3647 ir1 = 0 !1 3648 ir2 = 0 !1 3649 ir3 = 0 !1 3650 ! Go through reciprocal vectors to build u 3651 ABI_MALLOC(ptmp,(nqpts,nqpts,2)) 3652 do ip1=1,nr1*nr2*nr3 3653 gtmp = sqrt(sum(gvecs(:,ip1)**2)) 3654 ! Interpolate kernel in reciprocal space 3655 if ( gtmp < gcut ) then 3656 3657 ! Interpolate phi in reciprocal space: 3658 ! * ptmp(:,:,1) = phi(g) 3659 ! * ptmp(:,:,2) = dphi(g)/dg 3660 ! Note: do this on the fly, to go as fast as possible (this is equivalent 3661 ! to a call of the 'splint' routine) 3662 ptmp(:,:,:) = zero 3663 dg = pi / my_vdw_params%rcut 3664 ig = int(gtmp / dg) + 1 3665 a1 = ig - gtmp / dg 3666 b1 = one - a1 3667 a2 = (3 * a1**2 - one) * dg / six 3668 b2 = (3 * b1**2 - one) * dg / six 3669 a3 = (a1**3 - a1) * dg**2 / six 3670 b3 = (b1**3 - b1) * dg**2 / six 3671 do iq2 = 1,nqpts 3672 do iq1 = 1,iq2 3673 ptmp(iq1,iq2,1) = a1 * phig(ig-1,iq1,iq2) + b1 * phig(ig,iq1,iq2) & 3674 & + a3 * d2phidg2(ig-1,iq1,iq2) + b3 * d2phidg2(ig,iq1,iq2) 3675 ptmp(iq1,iq2,2) = (phig(ig,iq1,iq2) - phig(ig-1,iq1,iq2)) / dg & 3676 & - a2 * d2phidg2(ig-1,iq1,iq2) + b2 * d2phidg2(ig,iq1,iq2) 3677 end do 3678 end do 3679 3680 do iq2 = 1,nqpts 3681 do iq1 = 1,iq2-1 3682 ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:) 3683 end do 3684 end do 3685 !---------my debug------------------------------------------------------------ 3686 ! write diagonal components of interpolated kernel at each point 3687 write(82,*) ip1, (ptmp(iq1,iq1,1),iq1=1,nqpts) 3688 !--------end my debug-------------------------------------------------------- 3689 ! Calculate contributions to integral in Fourier space: 3690 ! Eq(12) from RS09 3691 utmp(:) = matmul(t3dg1d(:,ip1),ptmp(:,:,1)) 3692 else 3693 3694 write(std_out,*) 'WARNING: gtmp > gcut!!! This is at' 3695 write(std_out,*) 'ip1=',ip1 3696 utmp(:) = (zero,zero) 3697 3698 end if ! gtmp < gcut 3699 3700 ! Reconstruct the integrand in 3D 3701 t3dg(ir1,ir2,ir3,:) = utmp(:) 3702 3703 ir1 = ir1 + 1 3704 if ( ir1 > nr1-1 ) then 3705 ir2 = ir2 + 1 3706 ir1 = 0 3707 end if 3708 if ( ir2 > nr2-1 ) then 3709 ir3 = ir3 + 1 3710 ir2 = 0 3711 end if 3712 end do !ip1=1,npts_rho 3713 3714 ! Fourier-transform back the integrand 3715 #if defined HAVE_FFTW3 3716 do iq1=1,nqpts 3717 call dfftw_plan_dft_c2r_3d(plan,nr1,nr2,nr3, & 3718 t3dg(:,:,:,iq1),t3dr(:,:,:,iq1),FFTW_ESTIMATE) 3719 call dfftw_execute(plan) 3720 call dfftw_destroy_plan(plan) 3721 end do 3722 #else 3723 ABI_ERROR('vdW-DF calculations require FFTW3 support') 3724 #endif 3725 ! Write u for comparison 3726 open (unit=80, file='ureal1d-check.dat', status='replace') 3727 ! Repack the integrand 3728 ip1 = 1 3729 !$OMP PARALLEL DO COLLAPSE(3) & 3730 !$OMP& DEFAULT(SHARED) PRIVATE(ir1,ir2,ir3) 3731 do ir3=0,nr3-1 3732 do ir2=0,nr2-1 3733 do ir1=0,nr1-1 3734 u1dtmp(:,ip1) = t3dr(ir1,ir2,ir3,1:nqpts) !real space u 3735 write(80,*) ip1, (u1dtmp(iq1,ip1),iq1=1,nqpts) 3736 ip1 = ip1 + 1 3737 end do 3738 end do 3739 end do 3740 ! ttmp2(:,ip1) corresponds to u_alpha,i at eq(11) 3741 ! from RS09. 3742 !$OMP END PARALLEL DO 3743 3744 3745 3746 ! close(unit=63) 3747 ! close(unit=64) 3748 close(unit=65) 3749 ! close(unit=66) 3750 close(unit=80) 3751 close(unit=82) 3752 ABI_FREE(extin) 3753 ABI_FREE(grhotmp) 3754 ABI_FREE(ztmp) 3755 ABI_FREE(ptmp) 3756 ABI_FREE(decdrho_tmp) 3757 ABI_FREE(decdgrho_tmp) 3758 ABI_FREE(t3dr) 3759 ABI_FREE(t3dr2) 3760 ABI_FREE(t3dg) 3761 ABI_FREE(t3dr1d) 3762 ABI_FREE(t3dg1d) 3763 ABI_FREE(gprimdt) 3764 ABI_FREE(gvecs) 3765 ABI_FREE(utmp) 3766 ABI_FREE(u1dtmp) 3767 !--------end my debug-------------------------------------------- 3768 3769 end if 3770 3771 DBG_EXIT("COLL") 3772 3773 end subroutine vdw_df_internal_checks
ABINIT/m_xc_vdw/vdw_df_interpolate [ Functions ]
NAME
vdw_df_interpolate
FUNCTION
Calculates values of Phi(d1,d2) by interpolating the kernel.
INPUTS
d1= first coordinate d2= second coordinate sofswt= switch for the kernel softening
SOURCE
3107 function vdw_df_interpolate(d1,d2,sofswt) 3108 3109 use m_errors 3110 3111 !Arguments ------------------------------------ 3112 real(dp),intent(in) :: d1,d2 3113 integer,intent(in) :: sofswt 3114 3115 !Local variables ------------------------------ 3116 real(dp) :: vdw_df_interpolate 3117 integer :: id1,id2,ix1,ix2,ndpts 3118 real(dp) :: dcut,xd(4,4) 3119 3120 ! ************************************************************************* 3121 3122 vdw_df_interpolate = zero 3123 3124 ndpts = size(dmesh(:)) 3125 dcut = dmesh(ndpts) 3126 3127 if ( (d1 >= dcut) .or. (d2 >= dcut) ) then 3128 return 3129 end if 3130 3131 ! Find closest indices in dmesh 3132 id1 = vdw_df_indexof(dmesh,ndpts,d1) 3133 id2 = vdw_df_indexof(dmesh,ndpts,d2) 3134 3135 ! Find intervals 3136 xd(1,:) = one 3137 xd(2,1) = (d1 - dmesh(id1)) / (dmesh(id1+1) - dmesh(id1)) 3138 xd(2,2) = (d2 - dmesh(id2)) / (dmesh(id2+1) - dmesh(id2)) 3139 xd(3,:) = xd(2,:)**2 3140 xd(4,:) = xd(2,:)**3 3141 3142 3143 select case ( sofswt ) 3144 3145 case (1) 3146 do ix2=1,4 3147 do ix1=1,4 3148 vdw_df_interpolate = vdw_df_interpolate + & 3149 & phi_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2) 3150 end do 3151 end do 3152 3153 case (0) 3154 do ix1=1,4 3155 do ix2=1,4 3156 vdw_df_interpolate = vdw_df_interpolate + & 3157 & phi_u_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2) 3158 end do 3159 end do 3160 3161 case default 3162 ABI_ERROR(" vdW-DF soft switch must be 0 or 1") 3163 3164 end select 3165 3166 end function vdw_df_interpolate
ABINIT/m_xc_vdw/vdw_df_kernel [ Functions ]
NAME
vdw_df_kernel
FUNCTION
Calculates the van der Waals kernel for specified d-coordinates. Decides whether to use direct integration of Eq.(14) of Dion et al., PRL 92, 246401 (2004) [[cite:Dion2004]], or to return a 4th-order polynomial for small distances.
INPUTS
d1= first coordinate d2= second coordinate dsoft= distance below which the kernel will be softened acutmin= minimum angular cut-off aratio= ratio between highest and lowest angular delta damax= maximum angular delta
SIDE EFFECTS
phisoft= value of the softened kernel at d=0 will be automatically set if negative
SOURCE
2450 function vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax) 2451 2452 !Arguments ------------------------------------ 2453 real(dp),intent(in) :: d1,d2,dsoft,phisoft,acutmin,aratio,damax 2454 2455 !Local variables ------------------------------ 2456 real(dp) :: vdw_df_kernel 2457 character(len=512) :: msg 2458 integer :: napts 2459 real(dp) :: deltad,dtol,dtmp,d1m,d1p,d2m,d2p,phid,phim,phip,phix 2460 2461 ! ************************************************************************* 2462 2463 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2464 DBG_ENTER("COLL") 2465 #endif 2466 2467 ! Init 2468 dtol = my_vdw_params%tolerance 2469 dtmp = sqrt(d1**2 + d2**2) 2470 2471 if ( dtmp < dsoft ) then 2472 2473 ! Calculate a softened version of phi when (d1,d2) -> (0,0) 2474 ! using a polynomial expansion for nonzero distances and an 2475 ! exponential fit on dsoft when d -> 0 too. 2476 ! Note 1: the exponential fit is inspired from the definition 2477 ! of phi0_soft in RS09. 2478 ! Note 2: contrary to RS09, deltad is proportional to dsoft. 2479 2480 if ( phisoft < 0 ) then 2481 ABI_ERROR('the softened kernel must be positive') 2482 end if 2483 2484 vdw_df_kernel = phisoft 2485 2486 if ( dtmp > zero ) then 2487 deltad = dsoft / 100.0_dp 2488 d1m = (dsoft - deltad) * d1 / dtmp 2489 d1p = (dsoft + deltad) * d1 / dtmp 2490 d2m = (dsoft - deltad) * d2 / dtmp 2491 d2p = (dsoft + deltad) * d2 / dtmp 2492 2493 phim = vdw_df_kernel_value(d1m,d2m,acutmin,aratio,damax) 2494 phip = vdw_df_kernel_value(d1p,d2p,acutmin,aratio,damax) 2495 phix = (phim + phip) / two 2496 phid = (phip - phim) / (two * deltad) 2497 2498 vdw_df_kernel = phisoft + & 2499 & (four * (phix - phisoft) - phid * dsoft) & 2500 & * dtmp**2 / (two * dsoft**2) + & 2501 & (two * (phisoft - phix) + phid * dsoft) & 2502 & * dtmp**4 / (two * dsoft**4) 2503 end if 2504 2505 else 2506 2507 vdw_df_kernel = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax) 2508 2509 end if ! dtmp < dsoft 2510 2511 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2512 DBG_EXIT("COLL") 2513 #endif 2514 2515 end function vdw_df_kernel
ABINIT/m_xc_vdw/vdw_df_kernel_value [ Functions ]
NAME
vdw_df_kernel_value
FUNCTION
Calculates the van der Waals kernel for specified d-coordinates. Uses direct integration of Eq.(14) of Dion et al., PRL 92, 246401 (2004) [[cite:Dion2004]].
INPUTS
d1= first coordinate d2= second coordinate acutmin= minimum angular cut-off aratio= ratio between highest and lowest angular delta damax= maximum angular delta
SOURCE
2535 function vdw_df_kernel_value(d1,d2,acutmin,aratio,damax) 2536 2537 !Arguments ------------------------------------ 2538 real(dp),intent(in) :: d1,d2,acutmin,aratio,damax 2539 2540 !Local variables ------------------------------ 2541 real(dp) :: vdw_df_kernel_value 2542 character(len=512) :: msg 2543 integer :: ia,ib,napts,metyp 2544 real(dp) :: aa,bb,atol,atmp,acut,btmp,dain,damin,tau,ttmp,wtmp 2545 real(dp),allocatable :: amesh(:),amesh_cos(:),amesh_sin(:),nu1(:),nu2(:) 2546 real(dp),allocatable :: dphida(:),dphidb(:),dphidd(:),ycoef(:,:),work(:) 2547 real(dp),allocatable :: eint(:) 2548 ! ************************************************************************* 2549 2550 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2551 DBG_ENTER("COLL") 2552 #endif 2553 2554 ! Create angular mesh 2555 ! 2556 ! Note: Contrary to RS09, we use a linear mesh, but leave 2557 ! the door open to other kinds of meshes. Indeed, the use of a 2558 ! logarithmic mesh is interesting only when both d1 and d2 are 2559 ! small, in which case the kernel is usually softened. 2560 2561 ! Init 2562 2563 !! Testing the same way in which amesh is obtained in siesta: 2564 !! napts = nint(max(acutmin,aratio*sqrt(d1**2+d2**2))) This is 2565 !! very different to the way the number of a-mesh points are determined 2566 !! in Siesta. Actually this almost corresponds to their definition of acut 2567 2568 acut = max(acutmin,aratio*sqrt(d1**2+d2**2)) 2569 2570 !! Obtaining the number of amesh points, it is introduced damin 2571 !! and dain: 2572 damin = 0.01d0 2573 dain = min(damax,0.1*sqrt(d1**2+d2**2)) 2574 dain = max(dain,damin) 2575 2576 if ( dain == damax ) then !Linear mesh 2577 metyp = 2 2578 napts = nint( acut / damax ) 2579 else 2580 bb = acut * dain / (damax-dain) 2581 aa = log( 1 + dain/bb ) 2582 napts = 1 + nint( log(1+acut/bb) / aa ) 2583 metyp = 1 2584 end if 2585 !! getting napts end 2586 2587 atol = my_vdw_params%tolerance 2588 2589 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2590 write(msg,'(1x,"vdw_df_kernel_value:",1x,"napts = ",i8)') napts 2591 call wrtout(std_out,msg,'COLL') 2592 #endif 2593 2594 ABI_MALLOC(amesh,(napts)) 2595 ABI_MALLOC(amesh_cos,(napts)) 2596 ABI_MALLOC(amesh_sin,(napts)) 2597 ABI_MALLOC(dphida,(napts)) 2598 ABI_MALLOC(dphidb,(napts)) 2599 ABI_MALLOC(nu1,(napts)) 2600 ABI_MALLOC(nu2,(napts)) 2601 ABI_MALLOC(ycoef,(3,napts)) 2602 ABI_MALLOC(work,(napts)) 2603 ABI_MALLOC(eint,(napts)) 2604 ! Init amesh 2605 !FIXME: allow for other mesh types 2606 !Now intruducing metyp, setting mesh_cutoff to acut, removing 2607 !mesh_inc=damax and including mesh_ratio= 2608 call vdw_df_create_mesh(amesh,napts,metyp,acut,mesh_ratio=damax/dain) 2609 2610 ! Init cos(a), sin(a), nu1(a), and nu2(a) 2611 tau = 4 * pi / 9 2612 !$OMP PARALLEL DO & 2613 !$OMP& IF ( napts > 511 ) & 2614 !$OMP& DEFAULT(SHARED) PRIVATE(ia) 2615 do ia = 1,napts 2616 atmp = amesh(ia) 2617 amesh_cos(ia) = cos(atmp) 2618 amesh_sin(ia) = sin(atmp) 2619 if ( atmp == zero ) then !it was changed the original condition: atmp < atol 2620 nu1(ia) = d1**2 / 2 / tau 2621 nu2(ia) = d2**2 / 2 / tau 2622 else 2623 if ( d1 == zero ) then !it was changed the original condition: d1 < atol 2624 nu1(ia) = atmp*atmp / 2 2625 else 2626 nu1(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d1)**2)) 2627 end if 2628 if ( d2 == zero ) then !it was changed the original condition: d2 < atol 2629 nu2(ia) = atmp*atmp / 2 2630 else 2631 nu2(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d2)**2)) 2632 end if 2633 end if 2634 end do 2635 !$OMP END PARALLEL DO 2636 2637 ! Integrate 2638 dphida(1) = 0 2639 !$OMP PARALLEL DO & 2640 !$OMP& IF ( napts > 255 ) & 2641 !$OMP& DEFAULT(SHARED) PRIVATE(ia,ib,atmp,btmp,ttmp,wtmp) 2642 do ia = 2,napts 2643 atmp = amesh(ia) 2644 dphidb(1) = 0 2645 2646 do ib = 2,napts 2647 btmp = amesh(ib) 2648 2649 ! eq.(15) DRSLL04 2650 ttmp = half * ( 1 / (nu1(ia) + nu1(ib)) + 1 / (nu2(ia) + nu2(ib)) ) * & 2651 & ( 1 / (nu1(ia) + nu2(ia)) / (nu1(ib) + nu2(ib)) + & 2652 & 1 / (nu1(ia) + nu2(ib)) / (nu2(ia) + nu1(ib)) ) 2653 2654 ! eq.(16) DRSLL04 2655 wtmp = two * ( (three - atmp*atmp) * btmp * amesh_cos(ib) * & 2656 & amesh_sin(ia) + (three - btmp*btmp) * atmp * amesh_cos(ia) * & 2657 & amesh_sin(ib) + (atmp*atmp + btmp*btmp - three) * & 2658 & amesh_sin(ia) * amesh_sin(ib) - & 2659 & three * atmp * btmp * amesh_cos(ia) * amesh_cos(ib) ) / & 2660 & (atmp * btmp)**3 2661 2662 dphidb(ib) = atmp*atmp * btmp*btmp * wtmp * ttmp 2663 2664 end do ! ib = 2,napts 2665 2666 ! call spline_integrate(dphida(ia),napts,damax,dphidb) 2667 call cspint(dphidb,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,dphida(ia)) 2668 end do ! ia = 2,napts 2669 !$OMP END PARALLEL DO 2670 2671 ! Final value of the kernel 2672 ! call spline_integrate(vdw_df_kernel_value,napts,damax,dphida) 2673 call cspint(dphida,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,vdw_df_kernel_value) 2674 vdw_df_kernel_value = vdw_df_kernel_value * 2 / pi**2 2675 2676 ! Clean-up the mess 2677 ABI_FREE(amesh) 2678 ABI_FREE(amesh_cos) 2679 ABI_FREE(amesh_sin) 2680 ABI_FREE(nu1) 2681 ABI_FREE(nu2) 2682 ABI_FREE(dphida) 2683 ABI_FREE(dphidb) 2684 ABI_FREE(ycoef) 2685 ABI_FREE(work) 2686 ABI_FREE(eint) 2687 2688 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 ) 2689 DBG_EXIT("COLL") 2690 #endif 2691 2692 end function vdw_df_kernel_value
ABINIT/m_xc_vdw/vdw_df_ldaxc [ Functions ]
NAME
vdw_df_ldaxc
FUNCTION
Calculates LDA-based XC energy density for other vdW routines.
INPUTS
ixc= functional to apply nspden= number of spin components npts_rho= number of density points rho= density grho2= density gradient
OUTPUT
ex= exchange energy ec= correlation energy vx= exchange potential vc= correlation potential
SOURCE
2717 subroutine vdw_df_ldaxc(npts_rho,nspden,ngrad,rho_grho, & 2718 & exc_lda,vxc_lda,vxcg_lda) 2719 2720 !Arguments ------------------------------------ 2721 integer,intent(in) :: nspden,npts_rho,ngrad 2722 real(dp),intent(in) :: rho_grho(npts_rho,nspden,ngrad) 2723 real(dp),intent(out) :: exc_lda(2,npts_rho),vxc_lda(2,npts_rho,nspden) 2724 real(dp),intent(out) :: vxcg_lda(2,3,npts_rho) 2725 2726 !Local variables ------------------------------ 2727 integer :: ii 2728 logical :: is_gga 2729 type(libxc_functional_type) :: aux_funcx(2),aux_funcc(2) 2730 real(dp),allocatable :: rho_tmp(:,:)!,grho2(:,:) 2731 2732 ! ************************************************************************* 2733 2734 DBG_ENTER("COLL") 2735 2736 ! is_gga=libxc_functionals_isgga(vdw_funcs) 2737 2738 ABI_MALLOC(rho_tmp,(npts_rho,nspden)) 2739 ! if (is_gga) then 2740 ! ABI_MALLOC(grho2,(npts_rho,2*min(nspden,2)-1)) 2741 ! end if 2742 2743 ! Convert the quantities provided by ABINIT to the ones needed by LibXC 2744 ! 2745 ! Notes: 2746 ! 2747 ! * Abinit passes rho_up in the spin-unpolarized case, while LibXC 2748 ! expects the total density, but as we use rhonow, we don't need 2749 ! to convert anything. 2750 ! 2751 ! * In the case off gradients: 2752 ! 2753 ! - Abinit passes |grho_up|^2 while LibXC needs |grho_tot|^2; 2754 ! - Abinit passes |grho_up|^2, |grho_dn|^2, and |grho_tot|^2, 2755 ! while LibXC needs |grho_up|^2, grho_up.grho_dn, 2756 ! and |grho_dn|^2; 2757 ! 2758 ! but again the use of rho_grho lets us build these quantities 2759 ! directly. 2760 ! 2761 ! See ~abinit/src/56_xc/rhotoxc.F90 for details. 2762 ! However in routine libxc_functionals_getvxc of m_libxc_functionals the transformation 2763 ! to libxc input takes place, this means that we do not need to provide here the 2764 ! quatities for libxc. 2765 2766 if (nspden==1) then 2767 do ii=1,npts_rho 2768 rho_tmp(ii,1)=half*rho_grho(ii,1,1) 2769 ! if (is_gga) grho2(ii,1)=quarter*(rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2 & 2770 !& +rho_grho(ii,1,4)**2) 2771 end do 2772 else 2773 do ii=1,npts_rho 2774 rho_tmp(ii,1)=rho_grho(ii,2,1) 2775 rho_tmp(ii,2)=rho_grho(ii,1,1)-rho_grho(ii,2,1) 2776 ! if (is_gga) then 2777 ! grho2(ii,1)=rho_grho(ii,2,2)**2+rho_grho(ii,2,3)**2+rho_grho(ii,2,4)**2 2778 ! grho2(ii,2)=(rho_grho(ii,1,2)-rho_grho(ii,2,2))**2 & 2779 ! +(rho_grho(ii,1,3)-rho_grho(ii,2,3))**2 & 2780 ! +(rho_grho(ii,1,4)-rho_grho(ii,2,4))**2 2781 ! grho2(ii,3)=rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2+rho_grho(ii,1,4)**2 2782 ! end if 2783 end do 2784 end if 2785 2786 ! if (is_gga) then 2787 ! call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp,exc_lda,vxc_lda,& 2788 !& grho2=grho2,vxcgr=vxcg_lda,xc_functionals=vdw_funcs) 2789 ! else 2790 2791 ! Get LDA exchange energy and potential 2792 2793 call libxc_functionals_init(-001000,1,xc_functionals=aux_funcx) 2794 2795 write(std_out,*) 'aux_funcx=', aux_funcx(1)%id, aux_funcx(2)%id 2796 is_gga=libxc_functionals_isgga(aux_funcx) 2797 write(std_out,*) 'LDA exchange is GGA?', is_gga 2798 2799 call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp(:,nspden),exc_lda(1,:),vxc_lda(1,:,nspden),& 2800 & xc_functionals=aux_funcx) 2801 2802 call libxc_functionals_end(xc_functionals=aux_funcx) 2803 2804 ! Get LDA correlation energy and potential 2805 2806 call libxc_functionals_init(-000012,1,xc_functionals=aux_funcc) 2807 2808 write(std_out,*) 'aux_funcc=', aux_funcc(1)%id, aux_funcc(2)%id 2809 is_gga=libxc_functionals_isgga(aux_funcc) 2810 write(std_out,*) 'LDA correlation is GGA?', is_gga 2811 2812 call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp(:,nspden),exc_lda(2,:),vxc_lda(2,:,nspden),& 2813 & xc_functionals=aux_funcc) 2814 2815 2816 call libxc_functionals_end(xc_functionals=aux_funcc) 2817 2818 ! end if 2819 2820 ABI_FREE(rho_tmp) 2821 ! if (is_gga) then 2822 ! ABI_FREE(grho2) 2823 ! end if 2824 2825 DBG_EXIT("COLL") 2826 2827 end subroutine vdw_df_ldaxc
ABINIT/m_xc_vdw/vdw_df_set_tweaks [ Functions ]
NAME
vdw_df_set_tweaks
FUNCTION
Expands user-specified tweaks for internal use.
INPUTS
input_tweaks= condensed tweaks
OUTPUT
output_tweaks= expanded tweaks
SOURCE
3894 subroutine vdw_df_set_tweaks(input_tweaks,output_tweaks) 3895 3896 !Arguments ------------------------------------ 3897 integer,intent(in) :: input_tweaks 3898 type(vdw_df_tweaks_type),intent(out) :: output_tweaks 3899 3900 !Local variables------------------------------- 3901 3902 ! ************************************************************************* 3903 3904 output_tweaks%amesh_type = iand(ishft(input_tweaks,-10),3) 3905 output_tweaks%deriv_type = iand(ishft(input_tweaks,-2),3) 3906 output_tweaks%dmesh_type = iand(ishft(input_tweaks,-6),3) 3907 output_tweaks%interp_type = iand(ishft(input_tweaks,-4),3) 3908 output_tweaks%qmesh_type = iand(ishft(input_tweaks,-8),3) 3909 output_tweaks%run_type = iand(input_tweaks,3) 3910 3911 end subroutine vdw_df_set_tweaks
ABINIT/m_xc_vdw/vdw_df_tweaks_type [ Types ]
NAME
vdw_df_tweaks_type
FUNCTION
Tweaks for van der Waals XC calculations (use at your own risks).
SOURCE
158 type vdw_df_tweaks_type 159 160 integer :: amesh_type 161 ! A-mesh type 162 163 integer :: deriv_type 164 ! Derivation mode 165 166 integer :: dmesh_type 167 ! D-mesh type 168 169 integer :: interp_type 170 ! Interpolation mode 171 172 integer :: qmesh_type 173 ! Q-mesh type 174 175 integer :: run_type 176 ! Run mode 177 178 end type vdw_df_tweaks_type
ABINIT/m_xc_vdw/vdw_df_write_func [ Functions ]
NAME
vdw_df_write_func
FUNCTION
Writes function names for debugging purposes.
INPUTS
func_name= name of the function to print mode= write mode (see wrtout)
SOURCE
3927 subroutine vdw_df_write_func(func_name,mode) 3928 3929 !Arguments ------------------------------------ 3930 character(len=*),intent(in) :: func_name,mode 3931 3932 !Local variables------------------------------- 3933 character(len=512) :: msg 3934 3935 ! ************************************************************************* 3936 3937 write(msg,'(a,1x,a,1x,a,1x,a)') ch10,"=== subprogram :",func_name,"===" 3938 call wrtout(std_out,msg,mode) 3939 3940 end subroutine vdw_df_write_func
ABINIT/m_xc_vdw/xc_vdw_aggregate [ Functions ]
NAME
xc_vdw_aggregate
FUNCTION
Aggregates the various quantities calculated by the other vdW routines and produces energy, potential and stress tensor.
INPUTS
volume= the volume of the cell gprimd= unit vectors in reciprocal space npts_rho= number of density points to treat ngr2= number of components of grho2 nspden= number of spin components nr1= 1st dimension of real-space mesh nr2= 2nd dimension of real-space mesh nr3= 3rd dimension of real-space mesh rho= electronic density grho2= gradient of the density lda_ex= exchange energy computed at the LDA level lda_ec= correlation energy computed at the LDA level lda_vx= exchange potential computed at the LDA level lda_vc= correlation potential computed at the LDA level OUTPUTS dexc= vdW correction to the exchange-correlation energy dexcg= vdW correction to the exchange-correlation potential theta= see RS09
NOTES
exc_vdw includes deltae_vdw.
SOURCE
273 subroutine xc_vdw_aggregate(volume,gprimd,npts_rho,nspden,ngrad,nr1,nr2,nr3, & 274 & rho_grho,deltae_vdw,exc_vdw,decdrho_vdw,decdgrho_vdw,stress_vdw) 275 276 #if defined HAVE_FFTW3 277 include "fftw3.f03" 278 #endif 279 280 !Arguments ------------------------------------ 281 integer,intent(in) :: npts_rho,nspden,ngrad,nr1,nr2,nr3 282 real(dp),intent(in) :: volume,gprimd(3,3) 283 real(dp),intent(in) :: rho_grho(npts_rho,nspden,ngrad) 284 real(dp),intent(out) :: exc_vdw,deltae_vdw 285 real(dp),intent(out) :: decdrho_vdw(npts_rho,nspden),decdgrho_vdw(npts_rho,3,nspden) 286 real(dp),intent(out),optional :: stress_vdw(3,3) 287 288 !Local variables ------------------------------ 289 character(len=512) :: msg 290 integer :: ig,ip1,ip2,iq1,iq2,ir,ir1,ir2,ir3,is1,is2,ix,ngpts,nqpts,nrpts 291 integer(kind=SIZEOF_PTRDIFF_T) :: fftw3_plan 292 real(dp) :: a1,a2,a3,b1,b2,b3,gtmp,gcut,sg,dr 293 real(dp) :: ex,ec,vx,vc 294 real(dp) :: exc_nl,eps_vdw,deltae_uns,dexc,dexcg(3) 295 real(dp) :: dg,dvol,rho_tmp,gtol,grho_tmp(3),ngrho !last two for debug 296 real(dp) :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts) 297 real(dp) :: exc_tmp,decdrho_tmp(nspden),decdgrho_tmp(3,nspden) 298 real(dp),allocatable :: exc_lda(:,:),vxc_lda(:,:),vxcg_lda(:,:,:) 299 real(dp),allocatable :: gvec(:,:),theta(:,:,:) 300 real(dp),allocatable :: t3dr(:,:,:,:) 301 complex(dp),allocatable :: t3dg(:,:,:,:),ttmp(:,:),utmp(:) 302 real(dp),allocatable :: ptmp(:,:,:),wtmp(:),ttmp2(:,:) !,utmp(:) 303 304 ! ************************************************************************* 305 306 DBG_ENTER("COLL") 307 308 ! Init 309 nrpts = my_vdw_params%nrpts 310 ngpts = my_vdw_params%ngpts 311 nqpts = my_vdw_params%nqpts 312 gcut = my_vdw_params%gcut 313 gtol = my_vdw_params%tolerance 314 exc_vdw = zero 315 deltae_uns = zero 316 deltae_vdw = zero 317 if ( present(stress_vdw) ) stress_vdw(:,:) = zero 318 dr = my_vdw_params%rcut / nrpts 319 ztmp(:,:) = zero 320 decdgrho_tmp(:,:) = zero 321 decdrho_tmp(:) = zero 322 dvol = volume / npts_rho 323 324 ! Check the status of the switch 325 if ( .not. vdw_switch ) return 326 327 ABI_MALLOC(exc_lda,(2,npts_rho)) 328 ABI_MALLOC(vxc_lda,(2,npts_rho)) 329 ABI_MALLOC(vxcg_lda,(2,3,npts_rho)) 330 ABI_MALLOC(gvec,(3,npts_rho)) 331 ABI_MALLOC(theta,(nqpts,nspden,5)) 332 ABI_MALLOC(t3dr,(nr1,nr2,nr3,nqpts)) 333 ABI_MALLOC(t3dg,(nr1,nr2,nr3,nqpts)) 334 ABI_MALLOC(ttmp,(nqpts,nr1*nr2*nr3)) 335 ABI_MALLOC(ttmp2,(nqpts,nr1*nr2*nr3)) 336 ABI_MALLOC(ptmp,(nqpts,nqpts,2)) 337 ABI_MALLOC(utmp,(nqpts)) 338 ABI_MALLOC(wtmp,(nqpts)) 339 340 ! Calculate XC energy density from LDA / GGA 341 call vdw_df_ldaxc(npts_rho,nspden,ngrad,rho_grho,exc_lda,vxc_lda,vxcg_lda) 342 343 ! Pre-compute integrand for vdW energy correction 344 345 ztmp(:,:) = zero 346 do ir=1,nrpts 347 do iq1=1,nqpts 348 do iq2=1,iq1 349 ztmp(iq1,iq2) = ztmp(iq1,iq2) - & 350 & two * pi * dr * phir(ir,iq1,iq2) * (ir*dr)**2 351 !& two * pi * ( phir_u(ir,iq1,iq2) - phir(ir,iq1,iq2) ) * dr 352 end do 353 end do 354 end do 355 do iq1=2,nqpts 356 do iq2=1,iq1-1 357 ztmp(iq2,iq1) = ztmp(iq1,iq2) 358 end do 359 end do 360 361 !---------my debug---------------------------------- 362 ! open(34,file='table-localcorr.dat',status='replace') 363 ! do iq2=1,nqpts 364 ! do iq1=1,nqpts 365 ! write(34,*) iq1, iq2, ztmp(iq1,iq2) 366 ! enddo 367 ! enddo 368 ! close(34) 369 !--------end my debug------------------------------ 370 371 ! Build theta in 3D and g-vectors 372 if ( npts_rho /= nr1*nr2*nr3 ) then 373 ABI_WARNING('The 3D reconstruction of the density might be wrong (npts /= nr1*nr2*nr3)') 374 end if 375 #if defined DEBUG_VERBOSE 376 write(msg,'(a,1x,i8,1x,a)') "Will now call xc_vdw_energy", & 377 & nr1 * nr2 * nr3,"times" 378 ABI_COMMENT(msg) 379 #endif 380 !---my debug---------------------- 381 ! write(*,*) 'gprimd(i,1)=',(gprimd(ix,1),ix=1,3) 382 ! write(*,*) 'gprimd(i,2)=',(gprimd(ix,2),ix=1,3) 383 ! write(*,*) 'gprimd(i,3)=',(gprimd(ix,3),ix=1,3) 384 !----end my debug---------------- 385 ip1 = 1 386 ip2 = 1 387 do ir3=1,nr3 388 do ir2=1,nr2 389 do ir1=1,nr1 390 gvec(:,ip1) = two * pi * (gprimd(:,1) * ir1 + gprimd(:,2)& 391 & * ir2 + gprimd(:,3) * ir3) 392 ex = exc_lda(1,ip1) 393 ec = exc_lda(2,ip1) 394 vx = vxc_lda(1,ip1) 395 vc = vxc_lda(2,ip1) 396 theta(:,:,:) = zero 397 call xc_vdw_energy(nspden,rho_grho(ip1,1:nspden,1), & 398 & rho_grho(ip1,1:nspden,2:ngrad), & 399 & ex,ec,vx,vc,theta) 400 t3dr(ir1,ir2,ir3,1:nqpts) = theta(1:nqpts,1,1) 401 ! rho_tmp = sum(rho_grho(ip1,1:nspden,1)) 402 ! deltae_uns = deltae_uns + rho_tmp * eps_vdw * dvol 403 ip1 = ip1 + 1 404 end do 405 end do 406 #if defined DEBUG_VERBOSE 407 if ( (ip1 - ip2) > 100 ) then 408 write(msg,'(1x,a,1x,i3,"% complete")') & 409 & '[vdW-DF Energy]',int(ip1*100.0/(nr1*nr2*nr3)) 410 call wrtout(std_out,msg,'COLL') 411 ip2 = ip1 412 end if 413 #endif 414 end do 415 416 ! Fourier-transform theta 417 #if defined HAVE_FFTW3 418 do iq1=1,nqpts 419 call dfftw_plan_dft_r2c_3d(fftw3_plan,nr1,nr2,nr3, & 420 & t3dr(:,:,:,iq1),t3dg(:,:,:,iq1),FFTW_ESTIMATE) 421 call dfftw_execute(fftw3_plan) 422 call dfftw_destroy_plan(fftw3_plan) 423 t3dg(:,:,:,iq1) = t3dg(:,:,:,iq1) / (nr1 * nr2 * nr3) 424 end do 425 #else 426 ABI_ERROR('vdW-DF calculations require FFTW3 support') 427 #endif 428 429 ! Repack theta 430 ip1 = 1 431 do ir3=1,nr3 432 do ir2=1,nr2 433 do ir1=1,nr1 434 ttmp(:,ip1) = t3dg(ir1,ir2,ir3,:) 435 ip1 = ip1 + 1 436 end do 437 end do 438 end do 439 440 ! Reset 3D counters, since theta will be reconstructed in 3D on the fly 441 ir1 = 1 442 ir2 = 1 443 ir3 = 1 444 445 ! Go through reciprocal vectors 446 ! FIXME: get values of G-vectors 447 do ip1=1,npts_rho 448 if ( (ir3 > nr3) .and. (ir1 /= 1) .and. (ir2 /= 1) ) then 449 ABI_WARNING('buffer overflow reconstructing theta in 3D') 450 end if 451 452 gtmp = sqrt(sum(gvec(:,ip1)**2)) 453 454 if ( gtmp < gcut ) then 455 456 ! Interpolate phi in reciprocal space: 457 ! * ptmp(:,:,1) = phi(g) 458 ! * ptmp(:,:,2) = dphi(g)/dg 459 ! Note: do this on the fly, to go as fast as possible (this is equivalent 460 ! to a call of the 'splint' routine) 461 ptmp(:,:,:) = zero 462 dg = pi / my_vdw_params%rcut 463 ig = int(gtmp / dg) + 1 464 a1 = ig - gtmp / dg 465 b1 = one - a1 466 a2 = (3 * a1**2 - one) * dg / six 467 b2 = (3 * b1**2 - one) * dg / six 468 a3 = (a1**3 - a1) * dg**2 / six 469 b3 = (b1**3 - b1) * dg**2 / six 470 do iq2 = 1,nqpts 471 do iq1 = 1,iq2 472 ptmp(iq1,iq2,1) = a1 * phig(ig-1,iq1,iq2) + b1 * phig(ig,iq1,iq2) & 473 & + a3 * d2phidg2(ig-1,iq1,iq2) + b3 * d2phidg2(ig,iq1,iq2) 474 ptmp(iq1,iq2,2) = (phig(ig,iq1,iq2) - phig(ig-1,iq1,iq2)) / dg & 475 & - a2 * d2phidg2(ig-1,iq1,iq2) + b2 * d2phidg2(ig,iq1,iq2) 476 end do 477 end do 478 479 do iq2 = 1,nqpts 480 do iq1 = 1,iq2-1 481 ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:) 482 end do 483 end do 484 485 ! Calculate contributions to integral in Fourier space: 486 ! Eq(12) from RS09. 487 utmp(:) = matmul(ttmp(:,ip1),ptmp(:,:,1)) 488 489 ! Calculate contribution to stress in reciprocal space 490 ! Note: contribution of g=0 is zero 491 if ( present(stress_vdw) ) then 492 if ( gtmp > gtol ) then 493 wtmp(:) = matmul(ttmp(:,ip1),ptmp(:,:,2)) 494 sg = sum(wtmp(:)*ttmp(:,ip1)) * volume / gtmp 495 do is2=1,3 496 do is1=1,3 497 stress_vdw(is1,is2) = stress_vdw(is1,is2) - & 498 & sg * gvec(is1,ip1) * gvec(is2,ip1) 499 end do 500 end do 501 end if ! gtmp > gtol 502 end if ! present(stress_vdw) 503 504 else 505 506 utmp(:) = (zero,zero) 507 508 end if ! gtmp < gcut 509 510 ! Reconstruct the integrand in 3D 511 t3dg(ir1,ir2,ir3,:) = utmp(:) 512 513 ir1 = ir1 + 1 514 if ( ir1 > nr1 ) then 515 ir2 = ir2 + 1 516 ir1 = 1 517 end if 518 if ( ir2 > nr2 ) then 519 ir3 = ir3 + 1 520 ir2 = 1 521 end if 522 end do !ip1=1,npts_rho 523 524 ! Fourier-transform back the integrand 525 #if defined HAVE_FFTW3 526 !write(msg,'(a)') ch10 527 !call wrtout(std_out,msg,'COLL') 528 do iq1=1,nqpts 529 !write(msg,'(1x,a,i4.4)') "xc_vdw_aggregate: backward FFT #",iq1 530 !call wrtout(std_out,msg,'COLL') 531 532 call dfftw_plan_dft_c2r_3d(fftw3_plan,nr1,nr2,nr3, & 533 & t3dg(:,:,:,iq1),t3dr(:,:,:,iq1),FFTW_ESTIMATE) 534 call dfftw_execute(fftw3_plan) 535 call dfftw_destroy_plan(fftw3_plan) 536 end do 537 #else 538 ABI_ERROR('vdW-DF calculations require FFTW3 support') 539 #endif 540 541 ! Repack the integrand 542 ip1 = 1 543 !$OMP PARALLEL DO COLLAPSE(3) & 544 !$OMP& DEFAULT(SHARED) PRIVATE(ir1,ir2,ir3) 545 do ir3=1,nr3 546 do ir2=1,nr2 547 do ir1=1,nr1 548 ttmp2(:,ip1) = t3dr(ir1,ir2,ir3,1:nqpts) 549 ip1 = ip1 + 1 550 end do 551 end do 552 end do 553 ! ttmp2(:,ip1) corresponds to u_alpha,i at eq(11) 554 ! from RS09. 555 !$OMP END PARALLEL DO 556 557 #if defined DEBUG_VERBOSE 558 write(msg,'(a,1x,i8.8,1x,a)') "Will now call xc_vdw_energy",npts_rho,"times" 559 ABI_COMMENT(msg) 560 #endif 561 !--------my debug------------------------------------------ 562 open(unit=56,file='rho-LDA-xc.dat',status='replace') 563 open(unit=57,file='rho-eps.dat',status='replace') 564 open(unit=58,file='rho-grho-qpoly.dat',status='replace') 565 !--------end debug---------------------------------------- 566 ! Calculate and integrate vdW corrections at each point 567 do ip1=1,npts_rho 568 569 ! Get local contributions 570 ex = exc_lda(1,ip1) 571 ec = exc_lda(2,ip1) 572 vx = vxc_lda(1,ip1) 573 vc = vxc_lda(2,ip1) 574 theta(:,:,:) = zero 575 call xc_vdw_energy(nspden,rho_grho(ip1,1:nspden,1), & 576 & rho_grho(ip1,1:nspden,2:ngrad), & 577 & ex,ec,vx,vc,theta,exc_tmp,decdrho_tmp,decdgrho_tmp,ztmp) 578 579 ! Get nonlocal contributons 580 ! Note: is2 represents cartesian coordinates here. 581 if (nspden==1) then 582 rho_tmp = rho_grho(ip1,nspden,1) 583 else 584 rho_tmp = rho_grho(ip1,1,1) 585 end if 586 !rho_tmp = sum(rho_grho(ip1,1:nspden,1)) 587 deltae_vdw = deltae_vdw + rho_tmp * exc_tmp * dvol 588 exc_vdw = exc_vdw + rho_tmp * & 589 & ( half * ( sum(ttmp2(:,ip1) * theta(:,1,1)) / & 590 & (rho_tmp + tiny(rho_tmp)) ) * dvol ) 591 592 !----------------my debug-------------------------------------------------- 593 594 ! write(56,'(5e15.6)') rho_tmp, ex, ec, vx, vc 595 ! write(57,'(2e15.6)') rho_tmp, exc_tmp !ex, ec, vx, vc 596 597 !---------------end my debug----------------------------------------------- 598 599 !Correctness of multiplication by dvol above depends on how fftw3 deals 600 !with the volume element in direct or inverse FFT. Here it was assumed 601 !that fftw3 does not multiply by the volume upon FFT^-1 602 do is1=1,nspden 603 decdrho_vdw(ip1,is1) = decdrho_tmp(is1) + & 604 & sum(ttmp2(:,ip1) * theta(:,is1,2)) !CHECK how it is defined 605 !array decdrho_tmp(is1) 606 ! decdrho_vdw(is1) = decdrho_vdw(is1) + decdrho_tmp(is1) + & 607 !& sum(ttmp(:,ip1) * theta(:,is1,2)) 608 do is2=1,3 609 decdgrho_vdw(ip1,is2,is1) = decdgrho_tmp(is2,is1) + & 610 & sum(ttmp2(:,ip1) * theta(:,is1,is2+2)) 611 ! decdgrho_vdw(is2,is1) = decdgrho_vdw(is2,is1) + decdgrho_tmp(is2,is1) + & 612 !& sum(ttmp(:,ip1) * theta(:,is1,is2+2)) 613 end do 614 end do 615 616 #if defined DEBUG_VERBOSE 617 if ( modulo(ip1,int(npts_rho/100)) == 0 ) then 618 write(msg,'(1x,a,1x,i3,"% complete")') & 619 & "[vdW-DF Energy]",int(ip1*100.0/npts_rho) 620 call wrtout(std_out,msg,'COLL') 621 end if 622 #endif 623 624 end do ! Loop on npts_rho 625 !------my debug-------------------------------------------- 626 close (unit=56) 627 close (unit=57) 628 close (unit=58) 629 !------end debug------------------------------------------- 630 631 ! deltae_vdw = deltae_uns + deltae_vdw 632 633 #if defined DEBUG_VERBOSE 634 write(msg,'(1x,a)') "[vdW-DF Enrgy] 100% complete" 635 #endif 636 637 ! Display results 638 write(msg,'(a,1x,3a,2(3x,a,1x,e12.5,a),3x,a,1(1x,e12.5),a, & 639 & 3x,a,1(1x,e12.5),a)') & 640 & ch10,"[vdW-DF] xc_vdw_aggregate: reporting vdW-DF contributions", & 641 & ch10,ch10, & 642 & "DeltaE_vdw = ",deltae_vdw,ch10, & 643 & "Exc_vdw = ",exc_vdw,ch10 !, & 644 !& "dExc_vdw/drho = ",decdrho_vdw(:,:),ch10, & 645 !& "dExc_vdw/dgrho = ",decdgrho_vdw(:,:,:),ch10 646 call wrtout(std_out,msg,'COLL') 647 648 ! Final adjustments of stress 649 if ( present(stress_vdw) ) then 650 stress_vdw(:,:) = stress_vdw(:,:) / volume 651 652 write(msg,'(3x,a,a,3(5x,3(1x,e12.5),a))') & 653 & "Stress_vdw = ",ch10, & 654 & stress_vdw(1,:),ch10, & 655 & stress_vdw(2,:),ch10, & 656 & stress_vdw(3,:),ch10 657 call wrtout(std_out,msg,'COLL') 658 end if 659 660 ! Clean-up the mess 661 ABI_FREE(exc_lda) 662 ABI_FREE(vxc_lda) 663 ABI_FREE(vxcg_lda) 664 ABI_FREE(gvec) 665 ABI_FREE(theta) 666 ABI_FREE(t3dr) 667 ABI_FREE(t3dg) 668 ABI_FREE(ptmp) 669 ABI_FREE(ttmp) 670 ABI_FREE(ttmp2) 671 ABI_FREE(utmp) 672 ABI_FREE(wtmp) 673 674 DBG_EXIT("COLL") 675 676 end subroutine xc_vdw_aggregate
ABINIT/m_xc_vdw/xc_vdw_done [ Functions ]
NAME
xc_vdw_done
FUNCTION
Cleans-up the mess once van der Waals corrections to the energy are not needed anymore.
INPUTS
vdw_params= van der Waals parameters
SOURCE
925 subroutine xc_vdw_done(vdw_params) 926 927 !Arguments ------------------------------------ 928 type(xc_vdw_type), intent(in) :: vdw_params 929 930 !Local variables ------------------------------ 931 character(len=512) :: msg 932 integer :: i 933 934 ! ************************************************************************* 935 936 DBG_ENTER("COLL") 937 938 !--------my debug---------------------------- 939 !call libxc_functionals_end(xc_functionals=vdw_funcs) 940 !--------end my debug------------------------ 941 ABI_FREE(dmesh) 942 ABI_FREE(qmesh) 943 ABI_FREE(qpoly_basis) 944 ABI_FREE(phi) 945 ABI_FREE(phi_bicubic) 946 ABI_FREE(phi_u_bicubic) 947 ABI_FREE(phir) 948 ABI_FREE(phir_u) 949 ABI_FREE(d2phidr2) 950 ABI_FREE(phig) 951 ABI_FREE(d2phidg2) 952 953 DBG_EXIT("COLL") 954 955 end subroutine xc_vdw_done
ABINIT/m_xc_vdw/xc_vdw_energy [ Functions ]
NAME
xc_vdw_energy
FUNCTION
Calculates the exchange-correlation energy correction due to van der Waals interactions at one point.
INPUTS
nspden= number of spin components rho= electronic density grho= gradient of the density lda_ex= exchange energy computed at the LDA level lda_ec= correlation energy computed at the LDA level lda_vx= exchange potential computed at the LDA level lda_vc= correlation potential computed at the LDA level OUTPUTS theta= theta and its derivatives (see RS09) dexc= vdW correction to the exchange-correlation energy dexcg= vdW correction to the exchange-correlation potential
NOTES
The derivatives of theta are calculated only if the optional arguments are present. For performance reasons, this routine should not allocate any variable.
SOURCE
708 subroutine xc_vdw_energy(nspden,rho,grho,ex_lda,ec_lda,vx_lda,vc_lda, & 709 & theta,eps,dexc,dexcg,ztmp) 710 711 !Arguments ------------------------------------ 712 integer,intent(in) :: nspden 713 real(dp),intent(in) :: ex_lda,ec_lda,vx_lda,vc_lda 714 real(dp),intent(in) :: rho(nspden),grho(nspden,3) 715 real(dp),intent(inout) :: theta(my_vdw_params%nqpts,nspden,5) 716 real(dp),intent(in),optional :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts) 717 real(dp),intent(out),optional :: eps,dexc(nspden),dexcg(3,nspden) 718 719 !Local variables ------------------------------ 720 character(len=512) :: msg 721 logical :: calc_corrections 722 integer :: iq0,iq1,iq2,ir,is,ix,nqpts,nrpts,ns 723 real(dp) :: decdrho,dexdrho,dvcdrho,dvxdrho,ec,ex,vc,vx 724 real(dp) :: dr,kappa,rho_tmp,grho_tmp(3),grho2_tmp,rtol,qcut,zab 725 real(dp) :: dq(3),ptmp(2),q0(5),qtmp(2) 726 ! real(dp) :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts) 727 real(dp) :: qpoly(my_vdw_params%nqpts,3) 728 729 ! ************************************************************************* 730 731 ! Init 732 nqpts = my_vdw_params%nqpts 733 ! nrpts = my_vdw_params%nrpts 734 ns = my_vdw_params%nsmooth 735 qcut = my_vdw_params%qcut 736 rtol = my_vdw_params%tolerance 737 zab = my_vdw_params%zab 738 739 calc_corrections = .false. 740 if ( present(eps) .and. present(dexc) .and. present(dexcg) ) then 741 calc_corrections = .true. 742 end if 743 744 ! Sum density over spin 745 !rho_tmp = sum(rho(1:nspden)) 746 ! Get total density 747 if (nspden==1) then 748 rho_tmp = rho(nspden) 749 else 750 rho_tmp = rho(1) 751 end if 752 753 kappa = (three * pi**2 * rho_tmp)**third 754 forall(ix=1:3) grho_tmp(ix) = sum(grho(1:nspden,ix)) 755 grho2_tmp = sum(grho_tmp**2) 756 757 ! Calculate local wavevector q0 of eqs. 11-12 from DRSLL04 758 ! * q0(1) = q0 759 ! * q0(2) = dq0 / drho 760 ! * q0(3:5) = dq0 / dgrho 761 ! Notes: 762 ! * treating rho->0 separately (divide by 0) 763 ! * treating ex_lda->0 separately (divide by 0) 764 q0(:) = zero 765 ! if ( rho_tmp < rtol ) then 766 if ( rho_tmp < rtol ) then ! .or. (ex_lda < rtol) ) then 767 q0(1) = qcut 768 q0(2:5) = zero 769 else 770 q0(1) = kappa * (one + ec_lda/ex_lda - zab/nine * grho2_tmp / & 771 & (two * kappa * rho_tmp)**2) 772 773 if ( calc_corrections ) then 774 q0(2) = ( (vc_lda - ec_lda) / (rho_tmp * ex_lda) - ec_lda * & 775 & (vx_lda - ex_lda) / & 776 & (rho_tmp * ex_lda**2) + two * zab / nine * grho2_tmp / & 777 & (two * kappa * rho_tmp)**3 * & 778 & eight * kappa / three ) * kappa + q0(1) / (three * rho_tmp) 779 q0(3:5) = -two * kappa * (zab / nine) / (two * kappa * rho_tmp)**2 * & 780 & grho_tmp(1:3) 781 end if 782 ! Smoothen q0 near qcut exponentially eq. 5 RS09. (Saturation) 783 ! * q0s(1) = q0c * (1 - exp(-Sum_i=1,ns 1/i (q0/q0c)**i)) 784 ! * q0s(2) = dq0s / dq |q=q0 785 call vdw_df_saturation(q0(1),qcut,qtmp) 786 q0(1) = qtmp(1) 787 q0(2:5) = q0(2:5) * qtmp(2) 788 end if 789 790 !if ( q0(1) > qcut ) then 791 ! q0(1) = qcut 792 ! q0(2:5)= zero 793 !end if 794 795 !!qtmp(1) = (q0(1) / qcut) / ns 796 !!qtmp(2) = one / qcut 797 !!do is=ns-1,1,-1 798 !! qtmp(1) = (qtmp(1) + one / is) * q0(1) / qcut 799 !! qtmp(2) = (qtmp(2) * q0(1) + one) / qcut 800 !!end do 801 !!qtmp(2) = qtmp(2) * qcut * exp(-qtmp(1)) 802 803 !!q0(1) = qcut * (one - exp(-qtmp(1))) 804 !!q0(2:5) = q0(2:5) * qtmp(2) 805 806 ! Calculate polynomial coefficients for cubic-spline interpolation at q0 807 ! * qpoly(:,1) = coefficients 808 ! * qpoly(:,2) = first derivatives 809 qpoly(:,:) = zero 810 iq0 = vdw_df_indexof(qmesh(:),nqpts,q0(1)) 811 dq(1) = qmesh(iq0+1) - qmesh(iq0) 812 dq(2) = (qmesh(iq0+1) - q0(1)) / dq(1) 813 dq(3) = (q0(1) - qmesh(iq0)) / dq(1) 814 do iq1=1,nqpts 815 qpoly(iq1,1) = dq(2) * qpoly_basis(iq0,iq1,1) + dq(3) * & 816 & qpoly_basis(iq0+1,iq1,1) + & 817 & ((dq(2)**3 - dq(2)) * qpoly_basis(iq0,iq1,2) + & 818 & (dq(3)**3 - dq(3)) * qpoly_basis(iq0+1,iq1,2)) * dq(1)**2 / six 819 820 if ( calc_corrections ) then 821 qpoly(iq1,2) = -(qpoly_basis(iq0,iq1,1) - & 822 & qpoly_basis(iq0+1,iq1,1)) / dq(1) - & 823 & ((three * dq(2)**2 - one) * qpoly_basis(iq0,iq1,2) - & 824 & (three * dq(3)**2 - one) * qpoly_basis(iq0+1,iq1,2)) * dq(1) / six 825 end if 826 end do 827 828 ! Pre-compute integrand for vdW energy correction 829 ! by cubic-spline interpolation and store it in 830 ! qpoly(:,3). This is done twice: one for the 831 ! unsoftened kernel and other one for the 832 ! softened kernel 833 !dr = my_vdw_params%rcut / nrpts 834 !ztmp(:,:) = zero 835 836 ! if ( calc_corrections ) then 837 ! do ir=1,nrpts 838 ! do iq1=1,nqpts 839 ! do iq2=1,iq1 840 ! ztmp(iq1,iq2) = ztmp(iq1,iq2) - two * pi * phir(ir,iq1,iq2) * dr 841 ! end do 842 ! end do 843 ! end do 844 ! else 845 ! do ir=1,nrpts 846 ! do iq1=1,nqpts 847 ! do iq2=1,iq1 848 ! ztmp(iq1,iq2) = ztmp(iq1,iq2) + two * pi * phir_u(ir,iq1,iq2) * dr 849 ! end do 850 ! end do 851 ! end do 852 ! end if ! calc_corrections 853 854 ! do iq1=1,nqpts 855 ! do iq2=1,iq1-1 856 ! ztmp(iq2,iq1) = ztmp(iq1,iq2) 857 ! end do 858 ! end do 859 ! qpoly(:,3) = matmul(qpoly(:,1),ztmp(:,:)) 860 861 ! Calculate theta and its derivatives (see RS09) 862 ! * theta(:,:,1) = theta 863 ! * theta(:,:,2) = dtheta / drho 864 ! * theta(:,:,3:5) = dtheta / dgrho 865 ! Note: trick to go from Abinit to LibXC conventions 866 867 do is=1,nspden 868 !theta(:,is,1) = qpoly(:,1) * (2 - nspden) * rho(is) 869 theta(:,is,1) = qpoly(:,1) * rho(1) !CHECK THIS! 870 end do 871 !-------my debug-------------------------------------------------- 872 ! write(59,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts) 873 !-------end my debug---------------------------------------------- 874 ! Calculate theta derivatives 875 if ( calc_corrections ) then 876 do is=1,nspden 877 ! theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(is) 878 theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(1) 879 ! theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(is) 880 do ix=3,5 881 ! theta(:,is,ix) = qpoly(:,2) * q0(ix) * (2 - nspden) * rho(is) 882 theta(:,is,ix) = qpoly(:,2) * q0(ix) * rho(1) 883 end do 884 end do 885 end if 886 887 ! Calculate energy corrections 888 889 if ( calc_corrections ) then 890 !-------my debug---------------------------------------------- 891 ! write(58,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts) 892 !! write(59,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts) 893 !! write(60,*) rho_tmp, grho2_tmp, (theta(iq1,1,1),iq1=1,nqpts) 894 !! write(61,*) rho_tmp, grho2_tmp, (theta(iq1,1,2),iq1=1,nqpts) 895 !-------end my debug----------------------------------------- 896 qpoly(:,3) = matmul(qpoly(:,1),ztmp(:,:)) 897 eps = zero 898 ptmp(1) = sum(qpoly(:,3) * qpoly(:,1)) 899 eps = ptmp(1) * rho_tmp 900 dexc(:) = zero 901 dexcg(:,:) = zero 902 ptmp(2) = two * sum(qpoly(:,2) * qpoly(:,3)) 903 dexc(:) = two * ptmp(1) * rho_tmp + ptmp(2) * q0(2) * rho_tmp**2 904 do ix=3,5 905 dexcg(ix-2,:) = ptmp(2) * q0(ix) * rho_tmp**2 906 end do 907 end if ! calc_corrections 908 909 end subroutine xc_vdw_energy
ABINIT/m_xc_vdw/xc_vdw_get_params [ Functions ]
NAME
xc_vdw_get_params
FUNCTION
Exports internal vdW-DF parameters.
OUTPUT
vdw_params= van der Waals parameters
SOURCE
970 subroutine xc_vdw_get_params(vdw_params) 971 972 !Arguments ------------------------------------ 973 type(xc_vdw_type), intent(out) :: vdw_params 974 975 !Local variables ------------------------------ 976 character(len=512) :: msg 977 978 ! ************************************************************************* 979 980 vdw_params%functional = my_vdw_params%functional 981 vdw_params%zab = my_vdw_params%zab 982 vdw_params%ndpts = my_vdw_params%ndpts 983 vdw_params%dcut = my_vdw_params%dcut 984 vdw_params%dratio = my_vdw_params%dratio 985 vdw_params%dsoft = my_vdw_params%dsoft 986 vdw_params%phisoft = my_vdw_params%phisoft 987 vdw_params%nqpts = my_vdw_params%nqpts 988 vdw_params%qcut = my_vdw_params%qcut 989 vdw_params%qratio = my_vdw_params%qratio 990 vdw_params%nrpts = my_vdw_params%nrpts 991 vdw_params%rcut = my_vdw_params%rcut 992 vdw_params%rsoft = my_vdw_params%rsoft 993 vdw_params%ngpts = my_vdw_params%ngpts 994 vdw_params%gcut = my_vdw_params%gcut 995 vdw_params%acutmin = my_vdw_params%acutmin 996 vdw_params%aratio = my_vdw_params%aratio 997 vdw_params%damax = my_vdw_params%damax 998 vdw_params%damin = my_vdw_params%damin 999 vdw_params%nsmooth = my_vdw_params%nsmooth 1000 vdw_params%tolerance = my_vdw_params%tolerance 1001 vdw_params%tweaks = my_vdw_params%tweaks 1002 1003 end subroutine xc_vdw_get_params
ABINIT/m_xc_vdw/xc_vdw_init [ Functions ]
NAME
xc_vdw_init
FUNCTION
Calculates the van der Waals kernel.
INPUTS
vdw_params= parameters for the van der Waals calculations
SOURCE
1018 subroutine xc_vdw_init(vdw_params) 1019 1020 !Arguments ------------------------------------ 1021 type(xc_vdw_type), intent(in) :: vdw_params 1022 1023 !Local variables ------------------------------ 1024 character(len=*),parameter :: dmesh_file = "vdw_df_dmesh.pts" 1025 character(len=*),parameter :: qmesh_file = "vdw_df_qmesh.pts" 1026 1027 character(len=512) :: msg 1028 integer :: id1,id2,iq1,iq2,ir,ir0,ixc,jd1,jd2,ndpts,ngpts,nqpts,nrpts,ntest,sofswt 1029 integer :: ii,jj ! DEBUG 1030 real(dp) :: d1,d2,dcut,dd,delta_test,dsoft,phisoft,acutmin,aratio,damax 1031 real(dp) :: phimm,phimp,phipm,phipp,phipn,phi_tmp,ptmp(2) 1032 real(dp),allocatable :: deltd(:),kern_spline_der(:,:,:) 1033 real(dp),allocatable :: btmp(:,:),utmp(:),xde(:) 1034 1035 ! ************************************************************************* 1036 1037 DBG_ENTER("COLL") 1038 1039 ! Global init 1040 call xc_vdw_set_params(vdw_params) 1041 call vdw_df_set_tweaks(my_vdw_params%tweaks,my_vdw_tweaks) 1042 1043 ! Local init 1044 ndpts = my_vdw_params%ndpts 1045 nqpts = my_vdw_params%nqpts 1046 nrpts = my_vdw_params%nrpts 1047 acutmin = my_vdw_params%acutmin 1048 aratio = my_vdw_params%aratio 1049 damax = my_vdw_params%damax 1050 dsoft = my_vdw_params%dsoft 1051 phisoft = my_vdw_params%phisoft 1052 1053 ! Create d-mesh 1054 ! Note: Not avoid zero and make sure the last point is exactly dcut 1055 ABI_MALLOC(dmesh,(ndpts)) 1056 #if defined DEBUG_VERBOSE 1057 write(msg,'(1x,a)') "[vdW-DF Build] Generating D-mesh" 1058 call wrtout(std_out,msg,'COLL') 1059 #endif 1060 call vdw_df_create_mesh(dmesh,ndpts,my_vdw_tweaks%dmesh_type, & 1061 & my_vdw_params%dcut,mesh_ratio=my_vdw_params%dratio, & 1062 & mesh_tolerance=my_vdw_params%tolerance,mesh_file=dmesh_file, & 1063 & avoid_zero=.false.,exact_max=.true.) 1064 !DEBUG 1065 ! write(*,*) '-----D-mesh:-------' 1066 ! do id1 = 1,ndpts 1067 ! write(*,*) id1, dmesh(id1) 1068 ! end do 1069 !END DEBUG 1070 ! Create q-mesh 1071 ABI_MALLOC(qmesh,(nqpts)) 1072 #if defined DEBUG_VERBOSE 1073 write(msg,'(1x,a)') "[vdW-DF Build] Generating Q-mesh" 1074 call wrtout(std_out,msg,'COLL') 1075 #endif 1076 call vdw_df_create_mesh(qmesh,nqpts,my_vdw_tweaks%qmesh_type, & 1077 & my_vdw_params%qcut,mesh_ratio=my_vdw_params%qratio, & 1078 & mesh_tolerance=my_vdw_params%tolerance, mesh_file=qmesh_file) 1079 1080 ! Build polynomial basis for cubic-spline interpolation 1081 ! * qpoly_basis(:,1) = polynomial basis 1082 ! * qpoly_basis(:,2) = second derivative 1083 ABI_MALLOC(qpoly_basis,(nqpts,nqpts,2)) 1084 ! ABI_MALLOC(btmp,(nqpts,nqpts)) 1085 ABI_MALLOC(utmp,(nqpts)) 1086 #if defined DEBUG_VERBOSE 1087 write(msg,'(1x,a)') "[vdW-DF Build] Generating polynomial basis" 1088 call wrtout(std_out,msg,'COLL') 1089 #endif 1090 ! btmp(:,:) = zero 1091 ! forall(iq1=1:nqpts) btmp(iq1,iq1) = one 1092 utmp(:) = zero 1093 qpoly_basis(:,:,:) = zero 1094 do iq1=1,nqpts 1095 qpoly_basis(iq1,iq1,1) = one 1096 qpoly_basis(1,iq1,2) = zero 1097 utmp(1) = zero 1098 do iq2=2,nqpts-1 1099 ptmp(1) = (qmesh(iq2)-qmesh(iq2-1)) / (qmesh(iq2+1)-qmesh(iq2-1)) 1100 ptmp(2) = ptmp(1) * qpoly_basis(iq2-1,iq1,2) + two 1101 qpoly_basis(iq2,iq1,2) = (ptmp(1) - one) / ptmp(2) 1102 utmp(iq2) = ( 6*( (qpoly_basis(iq2+1,iq1,1)-qpoly_basis(iq2,iq1,1))& 1103 & /(qmesh(iq2+1)-qmesh(iq2)) - (qpoly_basis(iq2,iq1,1)-& 1104 & qpoly_basis(iq2-1,iq1,1))/(qmesh(iq2)-qmesh(iq2-1)) ) /& 1105 & (qmesh(iq2+1)-qmesh(iq2-1)) - ptmp(1)*utmp(iq2-1) ) / ptmp(2) 1106 1107 ! ptmp(1) = (qmesh(iq2) - qmesh(iq2-1)) / & 1108 !& (qmesh(iq2+1) - qmesh(iq2-1)) 1109 ! ptmp(2) = ptmp(1) * qpoly_basis(iq1,iq2-1,2) + two 1110 ! qpoly_basis(iq1,iq2,2) = (ptmp(1) - one) / ptmp(2) 1111 ! utmp(iq2) = ( six * ((btmp(iq1,iq2+1) - btmp(iq1,iq2)) / & 1112 !& (qmesh(iq2+1) - qmesh(iq2)) - (btmp(iq1,iq2) - & 1113 !& btmp(iq1,iq2-1)) / (qmesh(iq2) - qmesh(iq2-1))) / & 1114 !& (qmesh(iq2+1) - qmesh(iq2-1)) - ptmp(1) * utmp(iq2-1)) / & 1115 !& ptmp(2) 1116 end do 1117 ! utmp(nqpts) = zero 1118 qpoly_basis(nqpts,iq1,2) = zero 1119 do iq2=nqpts-1,1,-1 1120 qpoly_basis(iq2,iq1,2) = qpoly_basis(iq2,iq1,2) * & 1121 & qpoly_basis(iq2+1,iq1,2) + utmp(iq2) 1122 end do 1123 end do 1124 ! ABI_DEALLOCATE(btmp) 1125 ABI_FREE(utmp) 1126 1127 !DEBUG----------------------- 1128 ! open(37,file='qmesh-pofq-d2pdq2-check.dat',status='replace') 1129 ! do iq1=1,nqpts 1130 ! do iq2=1,nqpts 1131 ! write(37,'(4e15.6)') qmesh(iq1), qmesh(iq2), qpoly_basis(iq1,iq2,1) & 1132 ! , qpoly_basis(iq1,iq2,2) 1133 ! end do 1134 ! end do 1135 !END DEBUG------------------ 1136 1137 ! Create kernel and its derivatives 1138 ! Note: using 4 neighbours for derivatives 1139 ABI_MALLOC(phi,(ndpts,ndpts,4)) 1140 #if defined DEBUG_VERBOSE 1141 write(msg,'(1x,a)') "[vdW-DF Build] Building kernel and its derivatives" 1142 call wrtout(std_out,msg,'COLL') 1143 #endif 1144 1145 if ( phisoft < zero ) then 1146 phisoft = 0.4_dp !three * half * exp(-dsoft * two / sqrt(two)) 1147 my_vdw_params%phisoft = phisoft 1148 end if 1149 1150 ! Building of the softened kernel 1151 1152 sofswt = 1 1153 do id1=1,ndpts 1154 d1 = dmesh(id1) 1155 phi(id1,id1,1) = vdw_df_kernel(d1,d1,dsoft,phisoft,acutmin,aratio,damax) 1156 1157 ! Delta(d) should be at least 10^-6 1158 ! WARNING: changed tol3 to tol6 and changed max to min 1159 dd = max(0.01_dp*d1,tol3) !min(my_vdw_params%dratio*d1,tol6) 1160 1161 phimm = vdw_df_kernel(d1-dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax) 1162 phipm = vdw_df_kernel(d1+dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax) 1163 phipp = vdw_df_kernel(d1+dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax) 1164 phimp = vdw_df_kernel(d1-dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax) 1165 1166 ! Using five point stencil formula for crossed derivative phi(id1,id2,4) 1167 phi(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd) 1168 phi(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd) 1169 phi(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2) 1170 1171 do id2=1,id1-1 1172 d2 = dmesh(id2) 1173 phi(id1,id2,1) = vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax) 1174 phi(id2,id1,1) = phi(id1,id2,1) 1175 1176 phimm = vdw_df_kernel(d1-dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax) 1177 phipm = vdw_df_kernel(d1+dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax) 1178 phipp = vdw_df_kernel(d1+dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax) 1179 phimp = vdw_df_kernel(d1-dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax) 1180 1181 ! Using five point stencil formula for crossed derivative phi(id1,id2,4) 1182 phi(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd) 1183 phi(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd) 1184 phi(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2) 1185 1186 phi(id2,id1,2) = phi(id1,id2,3) 1187 phi(id2,id1,3) = phi(id1,id2,2) 1188 phi(id2,id1,4) = phi(id1,id2,4) 1189 end do 1190 1191 write(msg,'(1x,a,1x,i3,"% complete")') & 1192 & '[vdW-DF Build]',int(id1*100.0/ndpts) 1193 call wrtout(std_out,msg,'COLL') 1194 #if defined DEBUG_VERBOSE 1195 call flush_unit(std_out) 1196 #endif 1197 end do 1198 write(msg,'(a)') " " 1199 call wrtout(std_out,msg,'COLL') 1200 1201 ! Set boundary conditions 1202 ! Note: will have to check that the borders are smooth enough 1203 ! These boundary conditions produce large discontinuities, now testing its removal 1204 phi(ndpts,:,:) = zero 1205 phi(:,ndpts,:) = zero 1206 phi(1,:,2) = zero 1207 phi(:,1,3) = zero 1208 phi(1,:,4) = zero 1209 phi(:,1,4) = zero 1210 1211 !DEBUG------ 1212 ! write(*,*) 'd1 -- d2 -- phi -- dphi/dd1 -- dphi/dd2 -- d2phi/dd1dd2' 1213 ! do id1 = 1,ndpts 1214 ! do id2 = 1,ndpts 1215 ! d1 = dmesh(id1) 1216 ! d2 = dmesh(id2) 1217 ! write(*,'(2f10.6,4e15.6)') d1, d2, phi(id1,id2,1), & 1218 ! phi(id1,id2,2), phi(id1,id2,3), phi(id1,id2,4) 1219 ! end do 1220 ! end do 1221 !END DEBUG------- 1222 ! Ar2 results show better convergence of E_vdW than when boundary conditions were used. 1223 ! Kernel plot show that there are not dicontinuities as well. 1224 1225 #if defined DEBUG_VERBOSE 1226 write(msg,'(1x,a)') "[vdW-DF Build] Building filtered kernel" 1227 call wrtout(std_out,msg,'COLL') 1228 #endif 1229 1230 ! Calculate coefficients for bicubic interpolation 1231 ABI_MALLOC(phi_bicubic,(4,4,ndpts,ndpts)) 1232 call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi(:,:,1),phi(:,:,2),& 1233 & phi(:,:,3),phi(:,:,4),phi_bicubic) 1234 1235 !DEBUG------------------------------ 1236 ! open (unit=41, file='phi-bicubic-check.dat',status='replace') 1237 ! do id1=1,ndpts 1238 ! do id2=1,ndpts 1239 ! write(41,*) id1, id2, ((phi_bicubic(ii,jj,id1,id2),ii=1,4),jj=1,4) 1240 ! end do 1241 ! end do 1242 ! close(unit=41) 1243 !END DEBUG------------------------ 1244 ! Build filtered kernel 1245 ABI_MALLOC(phir,(0:nrpts,nqpts,nqpts)) 1246 ABI_MALLOC(d2phidr2,(0:nrpts,nqpts,nqpts)) 1247 ABI_MALLOC(phig,(0:nrpts,nqpts,nqpts)) 1248 ABI_MALLOC(d2phidg2,(0:nrpts,nqpts,nqpts)) 1249 call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt) 1250 my_vdw_params%ngpts = ngpts 1251 1252 ! Find closest indices in dmesh 1253 ! FIXME: something is missing or should be removed here 1254 ! ABI_MALLOC(kern_spline_der,(ndpts,ndpts,3)) 1255 ! ABI_MALLOC(deltd,(2)) 1256 ! ABI_MALLOC(xde,(2)) 1257 1258 ! kern_spline_der(:,:,:) = zero 1259 1260 ! do jd1=1,ndpts 1261 ! do jd2=1,ndpts 1262 1263 ! d1 = dmesh(jd1) 1264 ! d2 = dmesh(jd2) 1265 1266 ! if ( (d1 < dcut) .or. (d2 < dcut) ) then 1267 ! id1 = vdw_df_indexof(dmesh,ndpts,d1) 1268 ! id2 = vdw_df_indexof(dmesh,ndpts,d2) 1269 1270 ! deltd(1) = dmesh(id1+1) - dmesh(id1) 1271 ! deltd(2) = dmesh(id2+1) - dmesh(id2) 1272 1273 ! xde(1) = (d1 - dmesh(id1)) / deltd(1) 1274 ! xde(2) = (d2 - dmesh(id2)) / deltd(2) 1275 1276 ! kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1) 1277 ! kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2) 1278 ! kern_spline_der(jd1,jd2,3) = & 1279 ! phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) ) 1280 ! end if 1281 1282 ! end do 1283 ! end do 1284 1285 ! ABI_FREE(kern_spline_der) 1286 ! ABI_FREE(deltd) 1287 ! ABI_FREE(xde) 1288 1289 ! Building of the unsoftened kernel 1290 1291 ! Create unsoftened kernel and its derivatives 1292 ! Note: using 4 neighbours for derivatives 1293 ABI_MALLOC(phi_u,(ndpts,ndpts,4)) 1294 #if defined DEBUG_VERBOSE 1295 write(msg,'(1x,a)') "[vdW-DF Build] Building unsoftened kernel and its derivatives" 1296 call wrtout(std_out,msg,'COLL') 1297 #endif 1298 !------my debug-------------------------------------------------------- 1299 phi_u(:,:,:) = zero 1300 1301 ! do id1=1,ndpts 1302 ! d1 = dmesh(id1) 1303 ! phi_u(id1,id1,1) = vdw_df_kernel_value(d1,d1,acutmin,aratio,damax) 1304 ! 1305 ! ! Delta(d) should be at least 10^-6 1306 ! ! WARNING: changed tol3 to tol6 and changed max to min 1307 ! dd = max(0.01_dp*d1,tol3) !min(my_vdw_params%dratio*d1,tol6) 1308 ! 1309 ! phimm = vdw_df_kernel_value(d1-dd,d1-dd,acutmin,aratio,damax) 1310 ! phipm = vdw_df_kernel_value(d1+dd,d1-dd,acutmin,aratio,damax) 1311 ! phipp = vdw_df_kernel_value(d1+dd,d1+dd,acutmin,aratio,damax) 1312 ! phimp = vdw_df_kernel_value(d1-dd,d1+dd,acutmin,aratio,damax) 1313 ! 1314 ! ! Using five point stencil formula for crossed derivative phi(id1,id2,4) 1315 ! phi_u(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd) 1316 ! phi_u(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd) 1317 ! phi_u(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2) 1318 ! 1319 ! do id2=1,id1-1 1320 ! d2 = dmesh(id2) 1321 ! phi_u(id1,id2,1) = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax) 1322 ! phi_u(id2,id1,1) = phi_u(id1,id2,1) 1323 ! 1324 ! phimm = vdw_df_kernel_value(d1-dd,d2-dd,acutmin,aratio,damax) 1325 ! phipm = vdw_df_kernel_value(d1+dd,d2-dd,acutmin,aratio,damax) 1326 ! phipp = vdw_df_kernel_value(d1+dd,d2+dd,acutmin,aratio,damax) 1327 ! phimp = vdw_df_kernel_value(d1-dd,d2+dd,acutmin,aratio,damax) 1328 ! 1329 ! ! Using five point stencil formula for crossed derivative phi(id1,id2,4) 1330 ! phi_u(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd) 1331 ! phi_u(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd) 1332 ! phi_u(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2) 1333 ! 1334 ! phi_u(id2,id1,2) = phi_u(id1,id2,2) 1335 ! phi_u(id2,id1,3) = phi_u(id1,id2,3) 1336 ! phi_u(id2,id1,4) = phi_u(id1,id2,4) 1337 ! end do 1338 ! 1339 ! write(msg,'(1x,a,1x,i3,"% complete")') & 1340 ! & '[vdW-DF-Uns Build]',int(id1*100.0/ndpts) 1341 ! call wrtout(std_out,msg,'COLL') 1342 ! #if defined DEBUG_VERBOSE 1343 ! call flush_unit(std_out) 1344 ! #endif 1345 ! end do 1346 ! write(msg,'(a)') " " 1347 ! call wrtout(std_out,msg,'COLL') 1348 ! 1349 ! ! Set boundary conditions 1350 ! ! Note: will have to check that the borders are smooth enough 1351 ! ! These boundary conditions produce large discontinuities, now testing its removal 1352 ! phi_u(ndpts,:,:) = zero 1353 ! phi_u(:,ndpts,:) = zero 1354 ! phi_u(1,:,2) = zero 1355 ! phi_u(:,1,3) = zero 1356 ! phi_u(1,:,4) = zero 1357 ! phi_u(:,1,4) = zero 1358 ! 1359 ! 1360 ! ! Ar2 results show better convergence of E_vdW than when boundary conditions were used. 1361 ! ! Kernel plot show that there are not dicontinuities as well. 1362 !-------end my debug------------------------------------------------------- 1363 #if defined DEBUG_VERBOSE 1364 write(msg,'(1x,a)') "[vdW-DF Build] Building filtered unsoftened kernel" 1365 call wrtout(std_out,msg,'COLL') 1366 #endif 1367 1368 ! Calculate coefficients for bicubic interpolation 1369 ABI_MALLOC(phi_u_bicubic,(4,4,ndpts,ndpts)) 1370 call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi_u(:,:,1),phi_u(:,:,2),& 1371 & phi_u(:,:,3),phi_u(:,:,4),phi_u_bicubic) 1372 1373 ! Build filtered kernel 1374 ABI_MALLOC(phir_u,(0:nrpts,nqpts,nqpts)) 1375 ! For the local correction we just need the kernel not its 1376 ! derivatives. 1377 sofswt = 0 1378 !ABI_MALLOC(d2phidr2,(nrpts,nqpts,nqpts)) 1379 !ABI_MALLOC(phig,(nrpts,nqpts,nqpts)) 1380 !ABI_MALLOC(d2phidg2,(nrpts,nqpts,nqpts)) 1381 call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt) 1382 ! my_vdw_params%ngpts = ngpts already defined above 1383 1384 ! The following is not needed for the local correction: 1385 ! Find closest indices in dmesh 1386 ! FIXME: something is missing or should be removed here 1387 ! ABI_MALLOC(kern_spline_der,(ndpts,ndpts,3)) 1388 ! ABI_MALLOC(deltd,(2)) 1389 ! ABI_MALLOC(xde,(2)) 1390 1391 ! kern_spline_der(:,:,:) = zero 1392 1393 ! do jd1=1,ndpts 1394 ! do jd2=1,ndpts 1395 1396 ! d1 = dmesh(jd1) 1397 ! d2 = dmesh(jd2) 1398 1399 ! if ( (d1 < dcut) .or. (d2 < dcut) ) then 1400 ! id1 = vdw_df_indexof(dmesh,ndpts,d1) 1401 ! id2 = vdw_df_indexof(dmesh,ndpts,d2) 1402 1403 ! deltd(1) = dmesh(id1+1) - dmesh(id1) 1404 ! deltd(2) = dmesh(id2+1) - dmesh(id2) 1405 1406 ! xde(1) = (d1 - dmesh(id1)) / deltd(1) 1407 ! xde(2) = (d2 - dmesh(id2)) / deltd(2) 1408 1409 ! kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1) 1410 ! kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2) 1411 ! kern_spline_der(jd1,jd2,3) = & 1412 !& phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) ) 1413 ! end if 1414 1415 ! end do 1416 ! end do 1417 1418 ! ABI_FREE(kern_spline_der) 1419 ! ABI_FREE(deltd) 1420 ! ABI_FREE(xde) 1421 1422 call vdw_df_internal_checks(1) 1423 1424 vdw_switch = .true. 1425 1426 DBG_EXIT("COLL") 1427 1428 end subroutine xc_vdw_init
ABINIT/m_xc_vdw/xc_vdw_libxc_init [ Functions ]
NAME
xc_vdw_libxc_init
FUNCTION
Initializes LibXC parameters for the LDA-based part of vdW-DF calculations.
INPUTS
ixc_vdw= vdW-DF functional to apply
OUTPUT
(only writing)
SIDE EFFECTS
Internal variable vdw_funcs is set according to specified ixc_vdw. Internal variable my_vdw_params receives the selected functional.
SOURCE
1451 subroutine xc_vdw_libxc_init(ixc_vdw) 1452 1453 !Arguments ------------------------------------ 1454 integer,intent(in) :: ixc_vdw 1455 1456 !Local variables ------------------------------ 1457 character(len=*),parameter :: c_vdw1 = "LDA_C_PW" 1458 character(len=*),parameter :: x_vdw1 = "LDA_X" !"GGA_X_PBE_R" 1459 character(len=*),parameter :: x_vdw2 = "GGA_X_PW86" 1460 character(len=*),parameter :: x_vdw3 = "GGA_X_C09X" 1461 integer :: i,ii,ic,ix,ixc 1462 character(len=500) :: msg 1463 1464 ! ************************************************************************* 1465 1466 DBG_ENTER("COLL") 1467 1468 ! Select LibXC functionals 1469 select case (ixc_vdw) 1470 case (1) 1471 ix = libxc_functionals_getid(x_vdw1) 1472 ic = libxc_functionals_getid(c_vdw1) 1473 case (2) 1474 ix = libxc_functionals_getid(x_vdw2) 1475 ic = libxc_functionals_getid(c_vdw1) 1476 case (3) 1477 ix = libxc_functionals_getid(x_vdw3) 1478 ic = libxc_functionals_getid(c_vdw1) 1479 ABI_ERROR('[vdW-DF] C09 not available for now') 1480 case default 1481 ABI_ERROR('[vdW-DF] invalid setting of vdw_xc') 1482 end select 1483 if ( (ix == -1) .or. (ic == -1) ) then 1484 ABI_ERROR('[vdW-DF] unable to set LibXC parameters') 1485 end if 1486 ixc = -(ix * 1000 + ic) 1487 1488 ! Propagate to internal parameters 1489 my_vdw_params%functional = ixc_vdw !explore what happens with this variable 1490 1491 ! XC functional init 1492 call libxc_functionals_init(ixc,1,xc_functionals=vdw_funcs) 1493 1494 !------------my debug------------ 1495 ! write(*,*) 'vdw_funcs=', vdw_funcs(1)%id, vdw_funcs(2)%id 1496 ! call libxc_functionals_end(xc_functionals=vdw_funcs) 1497 !-----------end my debug------- 1498 1499 DBG_EXIT("COLL") 1500 1501 end subroutine xc_vdw_libxc_init
ABINIT/m_xc_vdw/xc_vdw_memcheck [ Functions ]
NAME
xc_vdw_memcheck
FUNCTION
Estimates the memory to be used by the vdW-DF method.
INPUTS
unt= unit to write the data to vp= van der Waals parameters
SOURCE
1517 subroutine xc_vdw_memcheck(unt,vp) 1518 1519 !Arguments ------------------------------------ 1520 integer,intent(in) :: unt 1521 type(xc_vdw_type), intent(in), optional :: vp 1522 1523 !Local variables ------------------------------ 1524 character(len=1536) :: msg 1525 integer :: napts,ndpts,nqpts,nrpts,id1,id2 1526 real(dp) :: dmax,dtmp,mem_perm,mem_temp 1527 type(xc_vdw_type) :: my_vp 1528 1529 ! ************************************************************************* 1530 1531 if ( present(vp) ) then 1532 my_vp = vp 1533 else 1534 my_vp = my_vdw_params 1535 end if 1536 1537 ndpts = my_vp%ndpts 1538 nqpts = my_vp%nqpts 1539 nrpts = my_vp%nrpts 1540 1541 if ( .not. allocated(dmesh) ) return 1542 1543 dmax = zero 1544 do id1=1,ndpts 1545 do id2=1,id1 1546 dtmp = sqrt(dmesh(id1)**2+dmesh(id2)**2) 1547 if ( dtmp > dmax ) dmax = dtmp 1548 end do 1549 end do 1550 napts = nint(max(my_vp%acutmin,my_vp%aratio*dmax)) 1551 1552 mem_perm = ( & 1553 & 4 * ndpts * (1 + 5 * ndpts ) + & 1554 & 2 * nqpts + (2 + 4 * nrpts) * nqpts * nqpts & 1555 & ) * sizeof(one) / 1048576.0_dp 1556 1557 mem_temp = ( & 1558 & 7 * napts + & 1559 & 4 * nqpts + 3* nqpts * nqpts + & 1560 & 5* nrpts + 2 & 1561 & ) * sizeof(one) / 1048576.0_dp 1562 1563 write(msg,'(a,1x,3(a),12(3x,a,1x,i8,1x,"elts",a), & 1564 & a,10x,a,1x,f10.3,1x,"Mb",a,a,13(3x,a,1x,i8,1x,"elts",a), & 1565 & a,10x,a,1x,f10.3,1x,"Mb",a,a,14x,a,1x,f10.3,1x,"Mb",2(a),1x,2(a))') & 1566 & ch10, & 1567 & '==== vdW-DF memory estimation ====',ch10,ch10, & 1568 & 'd2phidg2 .........',nrpts*nqpts*nqpts,ch10, & 1569 & 'd2phidr2 .........',nrpts*nqpts*nqpts,ch10, & 1570 & 'dmesh ............',ndpts * 4,ch10, & 1571 & 'phi ..............',ndpts*ndpts*4,ch10, & 1572 & 'phi_u ............',ndpts*ndpts*4,ch10, & 1573 & 'phi_bicubic ......',4*4*ndpts*ndpts,ch10, & 1574 & 'phi_u_bicubic ....',4*4*ndpts*ndpts,ch10, & 1575 & 'phig .............',nrpts*nqpts*nqpts,ch10, & 1576 & 'phir .............',nrpts*nqpts*nqpts,ch10, & 1577 & 'phir_u ...........',nrpts*nqpts*nqpts,ch10, & 1578 & 'qmesh ............',nqpts * 4,ch10, & 1579 & 'qpoly_basis.......',nqpts * nqpts * 2,ch10,ch10, & 1580 & 'Permanent =',mem_perm,ch10,ch10, & 1581 & 'amesh ............',napts,ch10, & 1582 & 'amesh_cos ........',napts,ch10, & 1583 & 'amesh_sin ........',napts,ch10, & 1584 & 'btmp .............',nqpts*nqpts,ch10, & 1585 & 'dphida ...........',napts,ch10, & 1586 & 'dphidb ...........',napts,ch10, & 1587 & 'ftmp .............',2*(2*nrpts+1),ch10, & 1588 & 'nu1 ..............',napts,ch10, & 1589 & 'nu2 ..............',napts,ch10, & 1590 & 'qpoly ............',nqpts*3,ch10, & 1591 & 'utmp(q) ..........',nqpts,ch10, & 1592 & 'utmp(r) ..........',nrpts,ch10, & 1593 & 'wtmp .............',nqpts*nqpts*2,ch10,ch10, & 1594 & 'Temporary =',mem_temp,ch10,ch10, & 1595 & 'TOTAL =',mem_perm+mem_temp,ch10, & 1596 & ch10, & 1597 & '==================================',ch10 1598 1599 call wrtout(unt,msg,'COLL') 1600 1601 end subroutine xc_vdw_memcheck
ABINIT/m_xc_vdw/xc_vdw_read [ Functions ]
NAME
xc_vdw_read
FUNCTION
Reads vdW-DF variables from disk.
INPUTS
filename= file to read data from
TODO
design an extension for ETSF_IO
SOURCE
1619 subroutine xc_vdw_read(filename) 1620 1621 !Arguments ------------------------------------ 1622 character(len=*), intent(in) :: filename 1623 1624 !Local variables ------------------------------ 1625 character(len=512) :: msg 1626 character(len=64) :: format_string 1627 integer :: ncid,dimids(4),varids(13) 1628 integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id 1629 integer :: nderivs,npoly,ndpts,ngpts,nqpts,nrpts 1630 1631 ! ************************************************************************* 1632 1633 DBG_ENTER("COLL") 1634 1635 #if defined HAVE_NETCDF 1636 write(msg,'("Reading vdW-DF data from",1x,a)') trim(filename) 1637 call wrtout(std_out,msg,'COLL') 1638 1639 ! Open file for reading 1640 NETCDF_VDWXC_CHECK(nf90_open(trim(filename),NF90_NOWRITE,ncid=ncid)) 1641 1642 ! Get file format and generator 1643 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'file_format',format_string)) 1644 write(msg,'(3x,"File format:",1x,a)') trim(format_string) 1645 call wrtout(std_out,msg,'COLL') 1646 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'generator',msg)) 1647 write(msg,'(3x,"Generator:",1x,a)') trim(msg) 1648 call wrtout(std_out,msg,'COLL') 1649 1650 ! Get attributes 1651 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin)) 1652 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio)) 1653 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax)) 1654 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin)) 1655 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut)) 1656 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio)) 1657 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft)) 1658 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional)) 1659 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut)) 1660 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth)) 1661 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft)) 1662 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut)) 1663 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio)) 1664 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut)) 1665 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft)) 1666 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance)) 1667 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks)) 1668 NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab)) 1669 1670 ! Get dimensions 1671 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nderivs',nderivs_id)) 1672 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'npoly',npoly_id)) 1673 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ndpts',ndpts_id)) 1674 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ngpts',ngpts_id)) 1675 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nqpts',nqpts_id)) 1676 NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nrpts',nrpts_id)) 1677 1678 ! Get varids 1679 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qmesh',varids(2))) 1680 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'dmesh',varids(3))) 1681 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi',varids(4))) 1682 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u',varids(5))) 1683 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_bicubic',varids(6))) 1684 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u_bicubic',varids(7))) 1685 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir',varids(8))) 1686 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir_u',varids(9))) 1687 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidr2',varids(10))) 1688 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phig',varids(11))) 1689 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidg2',varids(12))) 1690 NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qpoly_basis',varids(13))) 1691 1692 ! DEBUG 1693 ! write(*,*) ch10,'Inside xc_vdw_read routine',ch10 1694 ! write(*,*) 'nrpts before =',nrpts,ch10 1695 ! write(*,*) 'varids(8)=',varids(8) 1696 ! write(*,*) 'ndpts before =',ndpts,ch10 1697 ! nrpts = nrpts + 1 1698 ! ENDDEBUG 1699 ! Get dimensions 1700 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nderivs_id,len=nderivs)) 1701 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,npoly_id,len=npoly)) 1702 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ndpts_id,len=ndpts)) 1703 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ngpts_id,len=ngpts)) 1704 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nqpts_id,len=nqpts)) 1705 NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nrpts_id,len=nrpts)) 1706 1707 ! DEBUG 1708 ! write(*,*) 'nrpts after =',nrpts,ch10 1709 ! write(*,*) 'ndpts after =',ndpts,ch10 1710 ! nrpts = nrpts + 1 1711 ! ENDDEBUG 1712 1713 ! Get qmesh 1714 ABI_MALLOC(qmesh,(nqpts)) 1715 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(2),qmesh)) 1716 1717 ! Get dmesh 1718 ABI_MALLOC(dmesh,(ndpts)) 1719 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(3),dmesh)) 1720 1721 ! Get phi 1722 ABI_MALLOC(phi,(ndpts,ndpts,nderivs)) 1723 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(4),phi)) 1724 1725 ! Get phi_u 1726 ABI_MALLOC(phi_u,(ndpts,ndpts,nderivs)) 1727 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(5),phi_u)) 1728 1729 ! Get phi_bicubic 1730 ABI_MALLOC(phi_bicubic,(nderivs,nderivs,ndpts,ndpts)) 1731 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(6),phi_bicubic)) 1732 1733 ! Get phi_u_bicubic 1734 ABI_MALLOC(phi_u_bicubic,(nderivs,nderivs,ndpts,ndpts)) 1735 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(7),phi_u_bicubic)) 1736 1737 ! Get phir 1738 ABI_MALLOC(phir,(0:nrpts-1,nqpts,nqpts)) 1739 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(8),phir)) 1740 1741 ! Get phir_u 1742 ABI_MALLOC(phir_u,(0:nrpts-1,nqpts,nqpts)) 1743 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(9),phir_u)) 1744 1745 ! Get d2phidr2 1746 ABI_MALLOC(d2phidr2,(0:nrpts-1,nqpts,nqpts)) 1747 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(10),d2phidr2)) 1748 1749 ! Get phig 1750 ABI_MALLOC(phig,(0:nrpts-1,nqpts,nqpts)) 1751 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(11),phig)) 1752 1753 ! Get d2phidg2 1754 ABI_MALLOC(d2phidg2,(0:nrpts-1,nqpts,nqpts)) 1755 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(12),d2phidg2)) 1756 1757 ! Get qpoly_basis 1758 ABI_MALLOC(qpoly_basis,(nqpts,nqpts,npoly)) 1759 NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(13),qpoly_basis)) 1760 1761 ! Close file 1762 NETCDF_VDWXC_CHECK(nf90_close(ncid)) 1763 1764 ! Update my_vdw_params 1765 my_vdw_params%ndpts = ndpts 1766 my_vdw_params%ngpts = ngpts 1767 my_vdw_params%nqpts = nqpts 1768 my_vdw_params%nrpts = nrpts-1 1769 1770 ! DEBUG 1771 ! write(*,*) ch10,'At the end of xc_vdw_read',ch10 1772 ! write(*,*) 'nrpts= ',nrpts 1773 ! END DEBUG 1774 1775 #else 1776 ABI_ERROR('reading vdW-DF variables requires NetCDF') 1777 #endif 1778 1779 DBG_EXIT("COLL") 1780 1781 end subroutine xc_vdw_read
ABINIT/m_xc_vdw/xc_vdw_set_functional [ Functions ]
NAME
xc_vdw_set_functional
FUNCTION
Sets the vdW-DF parameters directly related to the functional.
INPUTS
vdw_func= van der Waals functional vdw_zab= Zab parameter to use for the calculations
SOURCE
1797 subroutine xc_vdw_set_functional(vdw_func,vdw_zab) 1798 1799 !Arguments ------------------------------------ 1800 integer, intent(in) :: vdw_func 1801 real(dp), optional, intent(in) :: vdw_zab 1802 1803 !Local variables ------------------------------ 1804 character(len=512) :: msg 1805 1806 ! ************************************************************************* 1807 1808 my_vdw_params%functional = vdw_func 1809 if ( present(vdw_zab) ) then 1810 my_vdw_params%zab = vdw_zab 1811 end if 1812 1813 end subroutine xc_vdw_set_functional
ABINIT/m_xc_vdw/xc_vdw_set_params [ Functions ]
NAME
xc_vdw_set_params
FUNCTION
Imports external vdW-DF parameters. INPUT vdw_params= van der Waals parameters
SOURCE
1828 subroutine xc_vdw_set_params(vdw_params) 1829 1830 !Arguments ------------------------------------ 1831 type(xc_vdw_type), intent(in) :: vdw_params 1832 1833 !Local variables ------------------------------ 1834 character(len=512) :: msg 1835 1836 ! ************************************************************************* 1837 1838 my_vdw_params%functional = vdw_params%functional 1839 my_vdw_params%zab = vdw_params%zab 1840 my_vdw_params%ndpts = vdw_params%ndpts 1841 my_vdw_params%dcut = vdw_params%dcut 1842 my_vdw_params%dratio = vdw_params%dratio 1843 my_vdw_params%dsoft = vdw_params%dsoft 1844 my_vdw_params%phisoft = vdw_params%phisoft 1845 my_vdw_params%nqpts = vdw_params%nqpts 1846 my_vdw_params%qcut = vdw_params%qcut 1847 my_vdw_params%qratio = vdw_params%qratio 1848 my_vdw_params%nrpts = vdw_params%nrpts 1849 my_vdw_params%rcut = vdw_params%rcut 1850 my_vdw_params%rsoft = vdw_params%rsoft 1851 my_vdw_params%ngpts = vdw_params%ngpts 1852 my_vdw_params%gcut = vdw_params%gcut 1853 my_vdw_params%acutmin = vdw_params%acutmin 1854 my_vdw_params%aratio = vdw_params%aratio 1855 my_vdw_params%damax = vdw_params%damax 1856 my_vdw_params%damin = vdw_params%damin 1857 my_vdw_params%nsmooth = vdw_params%nsmooth 1858 my_vdw_params%tolerance = vdw_params%tolerance 1859 my_vdw_params%tweaks = vdw_params%tweaks 1860 1861 end subroutine xc_vdw_set_params
ABINIT/m_xc_vdw/xc_vdw_show [ Functions ]
NAME
xc_vdw_show
FUNCTION
Displays the parameters in use for the vdW-DF corrections.
INPUTS
unt= unit to write the data to vp= van der Waals parameters
SOURCE
1877 subroutine xc_vdw_show(unt,vp) 1878 1879 !Arguments ------------------------------------ 1880 integer,intent(in) :: unt 1881 type(xc_vdw_type), intent(in), optional :: vp 1882 1883 !Local variables ------------------------------ 1884 character(len=1536) :: msg 1885 type(xc_vdw_type) :: my_vp 1886 1887 ! ************************************************************************* 1888 1889 if ( present(vp) ) then 1890 my_vp = vp 1891 else 1892 my_vp = my_vdw_params 1893 end if 1894 1895 write(msg,'(a,1x,3(a),3x,a,1x,i9,a,3x,a,1x,f9.3,a, & 1896 & 5(3x,a,1x,i9,a),14(3x,a,1x,f9.3,a),1(3x,a,1x,i9,a),a,1x,a,a)') & 1897 & ch10, & 1898 & '==== vdW-DF internal parameters ====',ch10,ch10, & 1899 & 'functional .............',my_vp%functional,ch10, & 1900 & 'zab ....................',my_vp%zab,ch10, & 1901 & 'd-mesh points ..........',my_vp%ndpts,ch10, & 1902 & 'q-mesh points ..........',my_vp%nqpts,ch10, & 1903 & 'r-mesh points ..........',my_vp%nrpts,ch10, & 1904 & 'g-mesh points ..........',my_vp%ngpts,ch10, & 1905 & 'smoothing iterations ...',my_vp%nsmooth,ch10, & 1906 & 'd-mesh cut-off .........',my_vp%dcut,ch10, & 1907 & 'd-mesh ratio ...........',my_vp%dratio,ch10, & 1908 & 'd-mesh softening .......',my_vp%dsoft,ch10, & 1909 & 'phi softening ..........',my_vp%phisoft,ch10, & 1910 & 'q-mesh cut-off .........',my_vp%qcut,ch10, & 1911 & 'q-mesh ratio ...........',my_vp%qratio,ch10, & 1912 & 'r-mesh cut-off .........',my_vp%rcut,ch10, & 1913 & 'r-mesh softening .......',my_vp%rsoft,ch10, & 1914 & 'g-mesh cut-off .........',my_vp%gcut,ch10, & 1915 & 'a-mesh min cut-off .....',my_vp%acutmin,ch10, & 1916 & 'a-mesh ratio ...........',my_vp%aratio,ch10, & 1917 & 'a-mesh delta max .......',my_vp%damax,ch10, & 1918 & 'a-mesh delta min .......',my_vp%damin,ch10, & 1919 & 'global tolerance .......',my_vp%tolerance,ch10, & 1920 & 'tweaks .................',my_vp%tweaks,ch10, & 1921 & ch10, & 1922 & '====================================',ch10 1923 1924 call wrtout(unt,msg,'COLL') 1925 1926 end subroutine xc_vdw_show
ABINIT/m_xc_vdw/xc_vdw_status [ Functions ]
NAME
xc_vdw_status
FUNCTION
Returns the status of the main vdW-DF switch.
INPUTS
None
SOURCE
1941 function xc_vdw_status() 1942 1943 !Arguments ------------------------------------ 1944 1945 !Local variables ------------------------------ 1946 logical :: xc_vdw_status 1947 1948 ! ************************************************************************* 1949 1950 xc_vdw_status = vdw_switch 1951 1952 end function xc_vdw_status
ABINIT/m_xc_vdw/xc_vdw_trigger [ Functions ]
NAME
xc_vdw_trigger
FUNCTION
Switches on and off the calculation of vdW interactions.
INPUTS
condition= boolean condition to trigger the calculations
SOURCE
1967 subroutine xc_vdw_trigger(condition) 1968 1969 !Arguments ------------------------------------ 1970 logical, intent(in) :: condition 1971 1972 !Local variables ------------------------------ 1973 character(len=512) :: msg 1974 1975 ! ************************************************************************* 1976 1977 if ( my_vdw_params%functional /= 0 ) then 1978 1979 if ( vdw_switch ) then 1980 if (.not. condition) then 1981 write (msg,'(1x,a,a)') & 1982 & "[vdW-DF] xc_vdw_trigger: disabling xc_vdw_aggregate",ch10 1983 vdw_switch = condition 1984 else 1985 write (msg,'(1x,a,a)') & 1986 & "[vdW-DF] xc_vdw_trigger: keeping xc_vdw_aggregate enabled",ch10 1987 end if 1988 else 1989 if ( condition ) then 1990 write (msg,'(1x,a,a)') & 1991 & "[vdW-DF] xc_vdw_trigger: enabling xc_vdw_aggregate",ch10 1992 else 1993 write (msg,'(1x,a,a)') & 1994 & "[vdW-DF] xc_vdw_trigger: disabling xc_vdw_aggregate",ch10 1995 end if 1996 1997 vdw_switch = condition 1998 end if 1999 2000 call wrtout(std_out,msg,'COLL') 2001 end if 2002 2003 end subroutine xc_vdw_trigger
ABINIT/m_xc_vdw/xc_vdw_type [ Types ]
NAME
xc_vdw_type
FUNCTION
Tunable parameters for calculations relying upon van der Waals XC.
SOURCE
75 type xc_vdw_type 76 77 integer :: functional = 0 78 ! Functional to use for the van der Waals correction 79 ! (see the description of the vdw_xc input variable for possible values) 80 81 real(dp) :: zab = -0.8491_dp 82 ! Parameter of the vdW functional, introduced by DRSLL04 83 84 integer :: ndpts = 20 85 ! Number of points for the d-mesh 86 87 real(dp) :: dcut = 30.0_dp 88 ! Cut-off for the d-mesh 89 90 real(dp) :: dratio = 20.0_dp 91 ! Ratio between highest and lowest d 92 93 real(dp) :: dsoft = 1.0_dp 94 ! Distance within the d-mesh below which the kernel will be softened 95 96 real(dp) :: phisoft = -1.0_dp 97 ! Value of the softened kernel for d=0 98 ! Will be set automatically if negative (default behaviour) 99 100 integer :: nqpts = 30 101 ! Number of points of the q-mesh 102 103 real(dp) :: qcut = 5.0_dp 104 ! Cut-off for the q-mesh 105 106 real(dp) :: qratio = 20.0_dp 107 ! Ratio between highest and lowest q 108 109 integer :: nrpts = 2048 110 ! Number of points in real space for the kernel 111 112 real(dp) :: rcut = 100.0_dp 113 ! Real-space cut-off for the kernel 114 115 real(dp) :: rsoft = 0.0_dp 116 ! Radius below which the kernel will be softened 117 118 integer :: ngpts = -1 119 ! Number of wavevectors for the kernel 120 121 real(dp) :: gcut = 5.0_dp 122 ! Wavevector cut-off for the kernel (=sqrt(ecut)) 123 124 real(dp) :: acutmin = 10.0_dp 125 ! Minimum cut-off for the angular mesh 126 127 real(dp) :: aratio = 30.0_dp 128 ! Ratio between highest and lowest a !DEBUG this definition has to 129 ! be corrected since its use in the code is different. 130 131 real(dp) :: damax = 0.5_dp 132 ! Maximum variation in the angular mesh 133 134 real(dp) :: damin = 1.0e-2_dp 135 ! Minimum variation in the angular mesh 136 137 integer :: nsmooth = 12 138 ! Saturation level to smoothen q0 near qcut 139 140 real(dp) :: tolerance = 1.0e-15_dp 141 ! Global numerical tolerance for the boundary values of the kernel 142 143 integer :: tweaks = 0 144 ! Tweaks of the implementation (modify with extreme care) 145 146 end type xc_vdw_type
ABINIT/m_xc_vdw/xc_vdw_write [ Functions ]
NAME
xc_vdw_write
FUNCTION
Writes vdW-DF variables to disk.
INPUTS
filename= file to write data to
TODO
FIXME: design an extension for ETSF_IO
SOURCE
2021 subroutine xc_vdw_write(filename) 2022 2023 !Arguments ------------------------------------ 2024 character(len=*), intent(in) :: filename 2025 2026 !Local variables ------------------------------ 2027 character(len=512) :: msg 2028 integer :: ncid,dimids(4),varids(13) 2029 integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id 2030 2031 ! ************************************************************************* 2032 2033 DBG_ENTER("COLL") 2034 2035 #if defined HAVE_NETCDF 2036 write(msg,'(a,1x,"Writing vdW-DF data to",1x,a)') ch10,trim(filename) 2037 call wrtout(std_out,msg,'COLL') 2038 2039 ! Create file (overwriting any existing data) 2040 NETCDF_VDWXC_CHECK(nf90_create(trim(filename),NF90_CLOBBER,ncid)) 2041 2042 ! Write attributes 2043 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'file_format',vdw_format_string)) 2044 write(msg,'(a,1x,a,1x,"for",1x,a,1x,"(build",1x,a,")")') & 2045 & PACKAGE_NAME,ABINIT_VERSION,ABINIT_TARGET,ABINIT_VERSION_BUILD 2046 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'generator',msg)) 2047 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin)) 2048 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio)) 2049 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax)) 2050 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin)) 2051 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut)) 2052 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio)) 2053 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft)) 2054 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional)) 2055 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut)) 2056 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth)) 2057 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft)) 2058 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut)) 2059 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio)) 2060 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut)) 2061 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft)) 2062 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance)) 2063 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks)) 2064 NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab)) 2065 2066 ! Define dimensions 2067 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nderivs',4,nderivs_id)) 2068 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'npoly',2,npoly_id)) 2069 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ndpts',my_vdw_params%ndpts,ndpts_id)) 2070 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ngpts',my_vdw_params%ngpts,ngpts_id)) 2071 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nqpts',my_vdw_params%nqpts,nqpts_id)) 2072 NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nrpts',my_vdw_params%nrpts+1,nrpts_id)) 2073 2074 ! Define qmesh 2075 dimids(1) = nqpts_id 2076 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qmesh',NF90_DOUBLE,dimids(1),varids(2))) 2077 2078 ! Define dmesh 2079 dimids(1) = ndpts_id 2080 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'dmesh',NF90_DOUBLE,dimids(1),varids(3))) 2081 2082 ! Define phi 2083 dimids(1) = ndpts_id 2084 dimids(2) = ndpts_id 2085 dimids(3) = nderivs_id 2086 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi',NF90_DOUBLE,dimids(1:3),varids(4))) 2087 2088 ! Define phi_u 2089 dimids(1) = ndpts_id 2090 dimids(2) = ndpts_id 2091 dimids(3) = nderivs_id 2092 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u',NF90_DOUBLE,dimids(1:3),varids(5))) 2093 2094 ! Define phi_bicubic 2095 dimids(1) = nderivs_id 2096 dimids(2) = nderivs_id 2097 dimids(3) = ndpts_id 2098 dimids(4) = ndpts_id 2099 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_bicubic',NF90_DOUBLE,dimids(1:4),varids(6))) 2100 2101 ! Define phi_u_bicubic 2102 dimids(1) = nderivs_id 2103 dimids(2) = nderivs_id 2104 dimids(3) = ndpts_id 2105 dimids(4) = ndpts_id 2106 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u_bicubic',NF90_DOUBLE,dimids(1:4),varids(7))) 2107 2108 ! Define phir 2109 dimids(1) = nrpts_id 2110 dimids(2) = nqpts_id 2111 dimids(3) = nqpts_id 2112 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir',NF90_DOUBLE,dimids(1:3),varids(8))) 2113 2114 ! Define phir_u 2115 dimids(1) = nrpts_id 2116 dimids(2) = nqpts_id 2117 dimids(3) = nqpts_id 2118 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir_u',NF90_DOUBLE,dimids(1:3),varids(9))) 2119 2120 ! Define d2phidr2 2121 dimids(1) = nrpts_id 2122 dimids(2) = nqpts_id 2123 dimids(3) = nqpts_id 2124 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidr2',NF90_DOUBLE,dimids(1:3),varids(10))) 2125 2126 ! Define phig 2127 dimids(1) = nrpts_id 2128 dimids(2) = nqpts_id 2129 dimids(3) = nqpts_id 2130 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phig',NF90_DOUBLE,dimids(1:3),varids(11))) 2131 2132 ! Define d2phidg2 2133 dimids(1) = nrpts_id 2134 dimids(2) = nqpts_id 2135 dimids(3) = nqpts_id 2136 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidg2',NF90_DOUBLE,dimids(1:3),varids(12))) 2137 2138 ! Define qpoly_basis 2139 dimids(1) = nqpts_id 2140 dimids(2) = nqpts_id 2141 dimids(3) = npoly_id 2142 NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qpoly_basis',NF90_DOUBLE,dimids(1:3),varids(13))) 2143 2144 ! Stop definitions 2145 NETCDF_VDWXC_CHECK(nf90_enddef(ncid)) 2146 2147 ! Write variables 2148 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(2),qmesh)) 2149 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(3),dmesh)) 2150 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(4),phi)) 2151 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(5),phi_u)) 2152 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(6),phi_bicubic)) 2153 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(7),phi_u_bicubic)) 2154 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(8),phir)) 2155 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(9),phir_u)) 2156 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(10),d2phidr2)) 2157 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(11),phig)) 2158 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(12),d2phidg2)) 2159 NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(13),qpoly_basis)) 2160 2161 ! Close file 2162 NETCDF_VDWXC_CHECK(nf90_close(ncid)) 2163 #else 2164 ABI_WARNING('writing vdW-DF variables requires NetCDF - skipping') 2165 #endif 2166 2167 DBG_EXIT("COLL") 2168 2169 end subroutine xc_vdw_write
m_xc_vdw/vdw_df_netcdf_ioerr [ Functions ]
[ Top ] [ m_xc_vdw ] [ Functions ]
NAME
vdw_df_netcdf_ioerr
FUNCTION
Reports errors occurring when allocating memory.
INPUTS
array_name= name of the array not properly allocated. array_size= size of the array. file_name= name of the file where the allocation was performed. file_line= line number in the file where the allocation was performed.
NOTES
This routine is usually interfaced with the macros defined in abi_common.h and uses this information to define a line offset.
SOURCE
3847 subroutine vdw_df_netcdf_ioerr(ncerr,file_name,file_line) 3848 3849 !Arguments ------------------------------------ 3850 integer,intent(in) :: ncerr 3851 character(len=*),optional,intent(in) :: file_name 3852 integer,optional,intent(in) :: file_line 3853 3854 !Local variables------------------------------- 3855 character(len=1024) :: msg 3856 character(len=fnlen) :: my_file_name = "N/D" 3857 integer :: my_file_line = -1 3858 3859 ! ************************************************************************* 3860 3861 if ( present(file_name) ) then 3862 my_file_name = trim(file_name) 3863 end if 3864 if ( present(file_line) ) then 3865 my_file_line = file_line 3866 end if 3867 3868 #if defined HAVE_NETCDF 3869 write(msg,'(a,a,3x,a)') & 3870 & 'NetCDF returned the error:',ch10,trim(nf90_strerror(ncerr)) 3871 if ( ncerr /= NF90_NOERR ) then 3872 call die(msg,trim(my_file_name),my_file_line) 3873 end if 3874 #endif 3875 3876 end subroutine vdw_df_netcdf_ioerr
m_xc_vdw/vdw_df_saturation [ Functions ]
[ Top ] [ m_xc_vdw ] [ Functions ]
NAME
vdw_df_saturation
FUNCTION
Obtain a saturated value of q0 which is always inside the interpolation range (eq. 5 RS09.) * q0s = q0c * (1 - exp(-Sum_i=1,ns 1/i (q0/q0c)**i)) * dq0sdq0 = dq0s / dq |q=q0
INPUTS
q0 = q0 value to be saturated q0c = Saturation value, max q0 value OUTPUTS q0s(1) = saturated q0 value q0s(2) = derivative of q0s(1) wrt q at q0
NOTES
This routine is usually interfaced with the macros defined in abi_common.h and uses this information to define a line offset.
SOURCE
3800 subroutine vdw_df_saturation(q0,q0c,q0s) 3801 3802 !Arguments ------------------------------------ 3803 real(dp),intent(in) :: q0, q0c 3804 real(dp),intent(out) :: q0s(2) 3805 !integer,intent(in) :: ncerr 3806 !character(len=*),optional,intent(in) :: file_name 3807 !integer,optional,intent(in) :: file_line 3808 3809 !Local variables------------------------------- 3810 integer :: is,ns 3811 real(dp) :: ptm, dptmdq 3812 ! ************************************************************************* 3813 3814 ns = my_vdw_params%nsmooth 3815 3816 ptm = (q0 / q0c) / ns 3817 dptmdq = 1._dp / q0c 3818 do is=ns-1,1,-1 3819 ptm = (ptm + 1._dp / is) * q0 / q0c 3820 dptmdq = (dptmdq * q0 + 1) / q0c ! (dptmdq * q0 + one) / q0c 3821 end do 3822 q0s(1) = q0c * (1 - exp(-ptm)) !q0c * (one - exp(-ptm)) 3823 q0s(2) = q0c * dptmdq * exp(-ptm) 3824 3825 end subroutine vdw_df_saturation