TABLE OF CONTENTS


ABINIT/m_xc_vdw [ Modules ]

[ Top ] [ Modules ]

NAME

  m_xc_vdw

FUNCTION

  Calculates van der Waals corrections to exchange-correlation.

COPYRIGHT

  Copyright (C) 2010-2018 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

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 #include "abi_xc_vdw.h"
27 
28 module m_xc_vdw
29 
30  use defs_basis
31  use iso_c_binding
32  use m_abicore
33  use m_errors
34  use libxc_functionals
35  use m_splines
36 #ifdef HAVE_NETCDF
37  use netcdf
38 #endif
39 
40  use m_io_tools,      only : flush_unit, open_file
41  use m_numeric_tools, only : simpson_int, cspint
42  use m_integrals,     only : radsintr
43 
44  implicit none
45 
46 #if defined DEV_YP_VDWXC
47 
48   public :: &
49 &   xc_vdw_type, &
50 &   xc_vdw_aggregate, &
51 &   xc_vdw_energy, &
52 &   xc_vdw_done, &
53 &   xc_vdw_get_params, &
54 &   xc_vdw_init, &
55 &   xc_vdw_libxc_init, &
56 &   xc_vdw_memcheck, &
57 &   xc_vdw_read, &
58 &   xc_vdw_set_functional, &
59 &   xc_vdw_show, &
60 &   xc_vdw_status, &
61 &   xc_vdw_trigger, &
62 &   xc_vdw_write
63 
64   private

ABINIT/m_xc_vdw/vdw_df_create_mesh [ Functions ]

[ Top ] [ 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

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

2826 subroutine vdw_df_create_mesh(mesh,npts,mesh_type,mesh_cutoff, &
2827 & mesh_inc,mesh_ratio,mesh_tolerance,mesh_file,avoid_zero,exact_max)
2828 
2829 
2830 !This section has been created automatically by the script Abilint (TD).
2831 !Do not modify the following lines by hand.
2832 #undef ABI_FUNC
2833 #define ABI_FUNC 'vdw_df_create_mesh'
2834 !End of the abilint section
2835 
2836   implicit none
2837 
2838 !Arguments ------------------------------------
2839   integer,intent(in) :: npts,mesh_type
2840   real(dp),intent(in) :: mesh_cutoff
2841   real(dp),intent(out) :: mesh(npts)
2842   logical,optional,intent(in) :: avoid_zero,exact_max
2843   real(dp),optional,intent(in) :: mesh_inc,mesh_ratio,mesh_tolerance
2844   character(len=*),optional,intent(in) :: mesh_file
2845 
2846 !Local variables ------------------------------
2847   logical :: my_avoid_zero,my_exact_max
2848   integer :: im,unt
2849   real(dp) :: exp_inc,exp_ofs,exp_tol,lin_inc,log_inc,log_ofs,log_tol
2850   character(len=500) :: msg
2851 
2852 ! *************************************************************************
2853 
2854 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2855   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2856 #endif
2857 
2858   if ( present(avoid_zero) ) then
2859     my_avoid_zero = avoid_zero
2860   else
2861     my_avoid_zero = .false.
2862   end if
2863 
2864   if ( present(exact_max) ) then
2865     my_exact_max = exact_max
2866   else
2867     my_exact_max = .false.
2868   end if
2869 
2870 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2871       write(msg,'(1x,a,a,1x,a,1x,i6,a,1x,a,i6,a,1x,a,1x,e10.3)') &
2872 &       "> Mesh parameters",ch10, &
2873 &       ">   npts    :",npts,ch10, &
2874 &       ">   type    :",mesh_type,ch10, &
2875 &       ">   cut-off :",mesh_cutoff
2876       call wrtout(std_out,msg,'COLL')
2877 #endif
2878 
2879   select case(mesh_type)
2880 
2881     case(0)
2882 
2883       if ( present(mesh_file) ) then
2884 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2885         write(msg,'(1x,a,1x,a)') ">   file    :",mesh_file
2886         call wrtout(std_out,msg,'COLL')
2887 #endif
2888         if (open_file(mesh_file,msg, newunit=unt,status="old",action="read") /= 0) then
2889           MSG_ERROR(msg)
2890         end if
2891         read(unt,'(e20.8)') mesh
2892         close(unit=unt)
2893       else
2894         MSG_ERROR("  You must specify a mesh file for this mesh type!")
2895       end if
2896 
2897     case(1)
2898 
2899       if ( present(mesh_ratio) ) then
2900         if ( present(mesh_tolerance) ) then
2901           log_tol = mesh_tolerance
2902         else
2903           log_tol = sqrt(my_vdw_params%tolerance)
2904         end if
2905 !DEBUG
2906 !        log_inc = mesh_ratio / (npts - 2)
2907         log_inc = log(mesh_ratio) / (npts - 1)
2908 !ENDDEBUG
2909         log_ofs = mesh_cutoff / (exp(log_inc*(npts - 1)) - 1)
2910 
2911 !!DEBUG
2912 !! Previous definition of log_ofs implies that the mesh
2913 !! points are to be calculated as:
2914 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 )
2915 !! which in turn means that x(0)=0 and is not consistent with the
2916 !! mesh(im) points given at line 3056.
2917 !! If the latter is adopted mesh_ratio above should be computed as:
2918 !! mesh_ratio=log ( (x(n)-x(n-1))/(x(2)-x(1)) ) and not merely as
2919 !! mesh_ratio= (x(n)-x(n-1))/(x(2)-x(1))
2920 !! in order to log_inc definition makes sense.
2921 !! The latter relation between log_inc and mesh_ratio is valid even if
2922 !! mesh points do not start at 0:
2923 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 ) + x_0
2924 !!ENDDEBUG
2925 
2926 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2927         write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') &
2928 &         ">   ratio   :",mesh_ratio,ch10, &
2929 &         ">   log_inc :",log_inc,ch10, &
2930 &         ">   log_ofs :",log_ofs,ch10, &
2931 &         ">   log_tol :",log_tol
2932         call wrtout(std_out,msg,'COLL')
2933 #endif
2934         if ( log_inc < log_tol ) then
2935           MSG_WARNING("  The logarithmic increment is very low!")
2936         end if
2937 
2938         do im = 1,npts
2939 !DEBUG
2940 !          mesh(im) = log_ofs * exp(log_inc * (im-1)) this is not consistent
2941 !       with definition of log_ofs at line 3025.
2942           mesh(im) = log_ofs * (exp(log_inc * (im-1)) - 1)
2943 !ENDDEBUG
2944         end do
2945 
2946         !FIXME: check the correctness of 0.2
2947         if ( my_avoid_zero ) then
2948           mesh(1) = exp(log_inc * 0.2) * mesh_cutoff / &
2949 &           (exp(log_inc*(npts - 1)) - one)
2950         end if
2951 
2952         if ( my_exact_max ) then
2953           mesh(npts) = mesh_cutoff
2954         end if
2955       else
2956         MSG_ERROR("  You must specify a mesh ratio for this mesh type!")
2957       end if
2958 
2959       ! Make sure the last point is qcut
2960       mesh(npts) = mesh_cutoff
2961 
2962     case(2)
2963 
2964       if ( present(mesh_inc) ) then
2965         lin_inc = mesh_inc
2966       else
2967         lin_inc = mesh_cutoff / npts
2968       end if
2969 
2970 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2971       write(msg,'(1x,a,1x,e10.3)') ">   lin_inc :",lin_inc
2972       call wrtout(std_out,msg,'COLL')
2973 #endif
2974       do im = 1,npts
2975         mesh(im) = lin_inc * (im - 1)
2976       end do
2977 
2978     case(3)
2979 
2980       if ( present(mesh_ratio) ) then
2981         if ( present(mesh_tolerance) ) then
2982           exp_tol = mesh_tolerance
2983         else
2984           exp_tol = sqrt(my_vdw_params%tolerance)
2985         end if
2986         exp_inc = mesh_ratio / (npts - 2)
2987         exp_ofs = mesh_cutoff * (log(exp_inc*(npts - 1)) + 1)
2988 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2989         write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') &
2990 &         ">   ratio   :",mesh_ratio,ch10, &
2991 &         ">   exp_inc :",exp_inc,ch10, &
2992 &         ">   exp_ofs :",exp_ofs,ch10, &
2993 &         ">   exp_tol :",exp_tol
2994         call wrtout(std_out,msg,'COLL')
2995 #endif
2996         if ( log_inc < log_tol ) then
2997           MSG_WARNING("  The logarithmic increment is very low!")
2998         end if
2999 
3000         do im = 1,npts
3001           mesh(im) = exp_ofs / log(exp_inc * (im-1))
3002         end do
3003       else
3004         MSG_ERROR("  You must specify a mesh ratio for this mesh type!")
3005       end if
3006 
3007     case default
3008 
3009       write(msg,'(2x,a,1x,i2)') "Unknown mesh type",mesh_type
3010       MSG_ERROR(msg)
3011 
3012   end select
3013 
3014 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
3015   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
3016 #endif
3017 
3018 end subroutine vdw_df_create_mesh

ABINIT/m_xc_vdw/vdw_df_filter [ Functions ]

[ Top ] [ 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

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

2238 subroutine vdw_df_filter(nqpts,nrpts,rcut,gcut,ngpts,sofswt)
2239 
2240 
2241 !This section has been created automatically by the script Abilint (TD).
2242 !Do not modify the following lines by hand.
2243 #undef ABI_FUNC
2244 #define ABI_FUNC 'vdw_df_filter'
2245 !End of the abilint section
2246 
2247   implicit none
2248 
2249 !Arguments ------------------------------------
2250   integer,intent(in) :: nqpts,nrpts,sofswt
2251   integer,intent(out) :: ngpts
2252   real(dp),intent(in) :: rcut,gcut
2253 
2254 !Local variables ------------------------------
2255   integer :: ig,iq1,iq2,ir
2256   real(dp) :: dg,dr,lstep,ptmp,q1,q2,xr,yq1,yqn,yr1,yrn
2257   real(dp),allocatable :: gmesh(:),rmesh(:),utmp(:)
2258 
2259 ! *************************************************************************
2260 
2261   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2262 
2263   ! Init
2264 
2265 
2266   dr = rcut / nrpts
2267   dg = pi / rcut
2268   ngpts = 1 + int(gcut/dg)
2269 
2270   ABI_ALLOCATE(utmp,(nrpts))
2271 
2272   ! Create radial mesh for radial FT: src/32_util/radsintr.F90
2273   ABI_ALLOCATE(rmesh,(nrpts))
2274   forall(ir=1:nrpts) rmesh(ir) = dr * dble(ir-1)
2275 
2276   if ( sofswt == 1 ) then
2277   ! Create reciprocal radial mesh
2278    ABI_ALLOCATE(gmesh,(ngpts))
2279    forall(ig=1:ngpts) gmesh(ig) = dg * dble(ig-1)
2280   end if
2281 
2282 
2283   ! Build filtered kernel for each (q1,q2) pair
2284   do iq1=1,nqpts
2285     do iq2=1,iq1
2286       q1 = qmesh(iq1)
2287       q2 = qmesh(iq2)
2288 
2289       if ( sofswt == 1 ) then
2290       ! Build kernel in real space
2291       ! Note: smoothly going to zero when approaching rcut
2292       ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values.
2293        do ir=1,nrpts
2294          xr=rmesh(ir)
2295          phir(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * &
2296 &          (one - ((ir - 1) / nrpts)**8)**4
2297        end do
2298 
2299       ! Obtain kernel in reciprocal space
2300        call radsintr(phir(:,iq1,iq2),phig(:,iq1,iq2),ngpts,nrpts,gmesh,rmesh,yq1,yqn)
2301 
2302       ! Filter in reciprocal space
2303        do ig=1,ngpts
2304          phig(ig,iq1,iq2) = phig(ig,iq1,iq2) * &
2305 &          (one - ((ig - 1) / ngpts)**8)**4
2306        end do
2307        phig(ngpts+1:nrpts,iq1,iq2) = zero
2308 
2309       ! Go back to real space
2310        call radsintr(phig(:,iq1,iq2),phir(:,iq1,iq2),nrpts,ngpts,rmesh,gmesh,yr1,yrn)
2311 
2312       ! Calculate second derivative in real space
2313        d2phidr2(1,iq1,iq2) = zero
2314        d2phidr2(nrpts,iq1,iq2) = zero
2315        utmp(1) = zero
2316        do ir=2,nrpts-1
2317          ptmp = half * d2phidr2(ir-1,iq1,iq2) + two
2318          d2phidr2(ir,iq1,iq2) = (half - one) / ptmp
2319          utmp(ir) = (three * (phir(ir+1,iq1,iq2) + phir(ir-1,iq1,iq2) - &
2320 &          two*phir(ir,iq1,iq2)) / (dr**2) - half * utmp(ir-1)) / ptmp
2321        end do
2322        do ir=nrpts-1,1,-1
2323          d2phidr2(ir,iq1,iq2) = d2phidr2(ir,iq1,iq2) * &
2324 &          d2phidr2(ir+1,iq1,iq2) + utmp(ir)
2325        end do
2326 
2327       ! Calculate second derivative in reciprocal space
2328        d2phidg2(1,iq1,iq2) = zero
2329        d2phidg2(ngpts,iq1,iq2) = zero
2330        utmp(1) = zero
2331        do ig=2,ngpts-1
2332          ptmp = half * d2phidg2(ig-1,iq1,iq2) + two
2333          d2phidg2(ig,iq1,iq2) = (half - one) / ptmp
2334          utmp(ig) = (three * (phig(ig+1,iq1,iq2) + phig(ig-1,iq1,iq2) - &
2335 &          two*phig(ig,iq1,iq2)) / (dr**2) - half * utmp(ig-1)) / ptmp
2336        end do
2337        do ig=ngpts-1,1,-1
2338          d2phidg2(ig,iq1,iq2) = d2phidg2(ig,iq1,iq2) * d2phidg2(ig+1,iq1,iq2) + &
2339 &          utmp(ig)
2340        end do
2341 
2342       ! Symmetrize kernels & derivatives
2343        phir(:,iq2,iq1) = phir(:,iq1,iq2)
2344        phig(:,iq2,iq1) = phig(:,iq1,iq2)
2345        d2phidr2(:,iq2,iq1) = d2phidr2(:,iq1,iq2)
2346        d2phidg2(:,iq2,iq1) = d2phidg2(:,iq1,iq2)
2347       end if ! sofswt == 1
2348 
2349       if ( sofswt == 0) then
2350       ! Build unsoftened kernel in real space
2351       ! Note: smoothly going to zero when approaching rcut
2352       ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values.
2353        do ir=1,nrpts
2354          xr=rmesh(ir)
2355          phir_u(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * &
2356 &          (one - ((ir - 1) / nrpts)**8)**4
2357        end do
2358       ! Symmetrize unsoftened kernel
2359        phir_u(:,iq2,iq1) = phir_u(:,iq1,iq2)
2360       end if ! sofswt == 0
2361 
2362      end do
2363    end do
2364 
2365   ABI_DEALLOCATE(utmp)
2366 
2367   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
2368 
2369 end subroutine vdw_df_filter

ABINIT/m_xc_vdw/vdw_df_indexof [ Functions ]

[ Top ] [ 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)

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

3044 function vdw_df_indexof(list_1d,npts,value)
3045 
3046 
3047 !This section has been created automatically by the script Abilint (TD).
3048 !Do not modify the following lines by hand.
3049 #undef ABI_FUNC
3050 #define ABI_FUNC 'vdw_df_indexof'
3051 !End of the abilint section
3052 
3053   implicit none
3054 
3055 !Arguments ------------------------------------
3056   integer,intent(in) :: npts
3057   real(dp),intent(in) :: value,list_1d(npts)
3058 
3059 !Local variables ------------------------------
3060   integer :: imin,imax,imid
3061   integer :: vdw_df_indexof
3062 
3063 ! *************************************************************************
3064 
3065   imin = 1
3066   imax = npts
3067 
3068   do
3069     imid = (imin + imax) / 2
3070     if ( value < list_1d(imid) ) then
3071       imax = imid
3072     else
3073       imin = imid
3074     end if
3075     if ( imin + 1 >= imax ) exit
3076   end do
3077 
3078   vdw_df_indexof = imin
3079 
3080 end function vdw_df_indexof

ABINIT/m_xc_vdw/vdw_df_internal_checks [ Functions ]

[ Top ] [ Functions ]

NAME

  vdw_df_internal_checks

FUNCTION

  Performs consistency checks of the vdW-DF kernel.

INPUTS

  test_mode= bitfield to enable/disable tests

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

3191 subroutine vdw_df_internal_checks(test_mode)
3192 
3193 
3194 !This section has been created automatically by the script Abilint (TD).
3195 !Do not modify the following lines by hand.
3196 #undef ABI_FUNC
3197 #define ABI_FUNC 'vdw_df_internal_checks'
3198 !End of the abilint section
3199 
3200   implicit none
3201 
3202 !Arguments ------------------------------------
3203   integer,intent(in) :: test_mode
3204 
3205 !Local variables ------------------------------
3206   character(len=512) :: msg
3207   integer :: ndpts,ntest,id1
3208   real(dp) :: acutmin,aratio,damax,dsoft,phisoft
3209   real(dp) :: delta_test,d1,d2
3210   real(dp),allocatable :: kern_test_c(:),kern_test_i(:)
3211   real(dp),allocatable :: intg_calc(:),intg_int(:)
3212 
3213 ! *************************************************************************
3214 
3215   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
3216 
3217   ! Here starts a test which will evaluate the integral over D for
3218   ! delta=0 (see Fig. 1 of DRSLL04).
3219   ! The integration will be performed over a finer and linear d mesh,
3220   ! not dmesh, for which the kernel values will be obtained by both
3221   ! interpolation and direct calculation. The range spanned by this mesh
3222   ! is half of dmesh range. Introduced variables: delta_test, kern_test_c,
3223   ! kern_test_i, ntest, intg_calc, intg_int.
3224   if ( iand(test_mode,1) /= 0 ) then
3225 
3226     write(msg,'(1x,a)') "[vdW-DF Check] Test #01: Kernel integral over D"
3227     call wrtout(std_out,msg,'COLL')
3228 
3229     acutmin = my_vdw_params%acutmin
3230     aratio = my_vdw_params%aratio
3231     damax = my_vdw_params%damax
3232     dsoft = my_vdw_params%dsoft
3233     ndpts = my_vdw_params%ndpts
3234     phisoft = my_vdw_params%phisoft
3235 
3236     ntest = 600
3237     delta_test = dmesh(ndpts) / (2 * ntest)
3238 
3239     ABI_ALLOCATE(kern_test_c,(ntest))
3240     ABI_ALLOCATE(kern_test_i,(ntest))
3241     ABI_ALLOCATE(intg_calc,(ntest))
3242     ABI_ALLOCATE(intg_int,(ntest))
3243 
3244 #if defined DEBUG_VERBOSE
3245     write(msg,'(a,4x,3(1x,a10),a)') ch10, "D","Calcul.","Interp.",ch10
3246     call wrtout(std_out,msg,'COLL')
3247 #endif
3248 
3249     do id1 = 1,ntest
3250       d1=id1*delta_test
3251       kern_test_c(id1) = (four_pi * d1**2) * &
3252 &     vdw_df_kernel_value(d1,d1,acutmin,aratio,damax)
3253       kern_test_i(id1) = (four_pi * d1**2) * vdw_df_interpolate(d1,d1,0)
3254 #if defined DEBUG_VERBOSE
3255       write(msg,'(4x,3(1x,e10.3))') id1*delta_test, &
3256 &       kern_test_c(id1), kern_test_i(id1)
3257       call wrtout(std_out,msg,'COLL')
3258 #endif
3259     end do
3260 
3261     call simpson_int(ntest,delta_test,kern_test_c,intg_calc)
3262     call simpson_int(ntest,delta_test,kern_test_i,intg_int)
3263 
3264     write(msg,'(a,4x,a,2(a,4x,a,e10.3))') ch10, &
3265 &     "> Integrals over D", ch10, ">   calculated   :",intg_calc(ntest),ch10, &
3266 &     ">   interpolated :",intg_int(ntest)
3267     call wrtout(std_out,msg,'COLL')
3268     if ( abs(intg_int(ntest) - intg_calc(ntest)) < tol3 ) then
3269       write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 PASSED"
3270     else
3271       write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 FAILED"
3272     end if
3273     call wrtout(std_out,msg,'COLL')
3274 
3275     ABI_DEALLOCATE(kern_test_c)
3276     ABI_DEALLOCATE(kern_test_i)
3277     ABI_DEALLOCATE(intg_calc)
3278     ABI_DEALLOCATE(intg_int)
3279 
3280   end if
3281 
3282   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
3283 
3284 end subroutine vdw_df_internal_checks

ABINIT/m_xc_vdw/vdw_df_interpolate [ Functions ]

[ Top ] [ 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

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

3103 function vdw_df_interpolate(d1,d2,sofswt)
3104 
3105   use m_errors
3106 
3107 !This section has been created automatically by the script Abilint (TD).
3108 !Do not modify the following lines by hand.
3109 #undef ABI_FUNC
3110 #define ABI_FUNC 'vdw_df_interpolate'
3111 !End of the abilint section
3112 
3113   implicit none
3114 
3115 !Arguments ------------------------------------
3116   real(dp),intent(in) :: d1,d2
3117   integer,intent(in) :: sofswt
3118 
3119 !Local variables ------------------------------
3120   real(dp) :: vdw_df_interpolate
3121   integer :: id1,id2,ix1,ix2,ndpts
3122   real(dp) :: dcut,xd(4,4)
3123 
3124 ! *************************************************************************
3125 
3126   vdw_df_interpolate = zero
3127 
3128   ndpts = size(dmesh(:))
3129   dcut = dmesh(ndpts)
3130 
3131   if ( (d1 >= dcut) .or. (d2 >= dcut) ) then
3132     return
3133   end if
3134 
3135   ! Find closest indices in dmesh
3136   id1 = vdw_df_indexof(dmesh,ndpts,d1)
3137   id2 = vdw_df_indexof(dmesh,ndpts,d2)
3138 
3139   ! Find intervals
3140   xd(1,:) = one
3141   xd(2,1) = (d1 - dmesh(id1)) / (dmesh(id1+1) - dmesh(id1))
3142   xd(2,2) = (d2 - dmesh(id2)) / (dmesh(id2+1) - dmesh(id2))
3143   xd(3,:) = xd(2,:)**2
3144   xd(4,:) = xd(2,:)**3
3145 
3146 
3147   select case ( sofswt )
3148 
3149     case (1)
3150       do ix1=1,4
3151         do ix2=1,4
3152           vdw_df_interpolate = vdw_df_interpolate + &
3153 &           phi_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2)
3154         end do
3155       end do
3156 
3157     case (0)
3158       do ix1=1,4
3159         do ix2=1,4
3160           vdw_df_interpolate = vdw_df_interpolate + &
3161 &           phi_u_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2)
3162         end do
3163       end do
3164 
3165     case default
3166       MSG_ERROR("  vdW-DF soft switch must be 0 or 1")
3167 
3168   end select
3169 
3170 end function vdw_df_interpolate

ABINIT/m_xc_vdw/vdw_df_kernel [ Functions ]

[ Top ] [ 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

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

2402 function vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax)
2403 
2404 
2405 !This section has been created automatically by the script Abilint (TD).
2406 !Do not modify the following lines by hand.
2407 #undef ABI_FUNC
2408 #define ABI_FUNC 'vdw_df_kernel'
2409 !End of the abilint section
2410 
2411   implicit none
2412 
2413 !Arguments ------------------------------------
2414   real(dp),intent(in) :: d1,d2,dsoft,phisoft,acutmin,aratio,damax
2415 
2416 !Local variables ------------------------------
2417   real(dp) :: vdw_df_kernel
2418   character(len=512) :: msg
2419   integer :: napts
2420   real(dp) :: deltad,dtol,dtmp,d1m,d1p,d2m,d2p,phid,phim,phip,phix
2421 
2422 ! *************************************************************************
2423 
2424 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2425   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2426 #endif
2427 
2428   ! Init
2429   dtol = my_vdw_params%tolerance
2430   dtmp = sqrt(d1**2 + d2**2)
2431 
2432   if ( dtmp < dsoft ) then
2433 
2434     ! Calculate a softened version of phi when (d1,d2) -> (0,0)
2435     ! using a polynomial expansion for nonzero distances and an
2436     ! exponential fit on dsoft when d -> 0 too.
2437     ! Note 1: the exponential fit is inspired from the definition
2438     !         of phi0_soft in RS09.
2439     ! Note 2: contrary to RS09, deltad is proportional to dsoft.
2440 
2441     if ( phisoft < 0 ) then
2442       MSG_ERROR('the softened kernel must be positive')
2443     end if
2444 
2445     vdw_df_kernel = phisoft
2446 
2447     if ( dtmp >= dtol ) then
2448       deltad = dsoft / 100.0_dp
2449       d1m = (dsoft - deltad) * d1 / dtmp
2450       d1p = (dsoft + deltad) * d1 / dtmp
2451       d2m = (dsoft - deltad) * d2 / dtmp
2452       d2p = (dsoft + deltad) * d2 / dtmp
2453 
2454       phim = vdw_df_kernel_value(d1m,d2m,acutmin,aratio,damax)
2455       phip = vdw_df_kernel_value(d1p,d2p,acutmin,aratio,damax)
2456       phix = (phim + phip) / two
2457       phid = (phip - phim) / (two * deltad)
2458 
2459       vdw_df_kernel = phisoft + &
2460 &       (four * (phix - phisoft) - phid * dsoft) &
2461 &       * dtmp**2 / (two * dsoft**2) + &
2462 &       (two  * (phisoft - phix) + phid * dsoft) &
2463 &       * dtmp**4 / (two * dsoft**4)
2464     end if
2465 
2466   else
2467 
2468     vdw_df_kernel = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)
2469 
2470   end if ! dtmp < dsoft
2471 
2472 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2473   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
2474 #endif
2475 
2476 end function vdw_df_kernel

ABINIT/m_xc_vdw/vdw_df_kernel_value [ Functions ]

[ Top ] [ 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

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

2502 function vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)
2503 
2504 
2505 !This section has been created automatically by the script Abilint (TD).
2506 !Do not modify the following lines by hand.
2507 #undef ABI_FUNC
2508 #define ABI_FUNC 'vdw_df_kernel_value'
2509 !End of the abilint section
2510 
2511   implicit none
2512 
2513 !Arguments ------------------------------------
2514   real(dp),intent(in) :: d1,d2,acutmin,aratio,damax
2515 
2516 !Local variables ------------------------------
2517   real(dp) :: vdw_df_kernel_value
2518   character(len=512) :: msg
2519   integer :: ia,ib,napts,metyp
2520   real(dp) :: aa,bb,atol,atmp,acut,btmp,dain,damin,tau,ttmp,wtmp
2521   real(dp),allocatable :: amesh(:),amesh_cos(:),amesh_sin(:),nu1(:),nu2(:)
2522   real(dp),allocatable :: dphida(:),dphidb(:),dphidd(:),ycoef(:,:),work(:)
2523   real(dp),allocatable :: eint(:)
2524 ! *************************************************************************
2525 
2526 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2527   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2528 #endif
2529 
2530   ! Create angular mesh
2531   !
2532   ! Note: Contrary to RS09, we use a linear mesh, but leave
2533   !       the door open to other kinds of meshes. Indeed, the use of a
2534   !       logarithmic mesh is interesting only when both d1 and d2 are
2535   !       small, in which case the kernel is usually softened.
2536 
2537   ! Init
2538 
2539 !!  Testing the same way in which amesh is obtained in siesta:
2540 !!  napts = nint(max(acutmin,aratio*sqrt(d1**2+d2**2))) This is
2541 !!  very different to the way the number of a-mesh points are determined
2542 !!  in Siesta.  Actually this almost corresponds to their definition of acut
2543 
2544   acut = max(acutmin,aratio*sqrt(d1**2+d2**2))
2545 
2546 !! Obtaining the number of amesh points, it is introduced damin
2547 !! and dain:
2548   damin = 0.01d0
2549   dain = min(damax,0.1*sqrt(d1**2+d2**2))
2550   dain = max(dain,damin)
2551 
2552   if ( dain == damax ) then !Linear mesh
2553    metyp = 2
2554    napts = nint( acut / damax )
2555   else
2556    bb = acut * dain / (damax-dain)
2557    aa = log( 1 + dain/bb )
2558    napts = 1 + nint( log(1+acut/bb) / aa )
2559    metyp = 1
2560   end if
2561 !! getting napts end
2562 
2563   atol = my_vdw_params%tolerance
2564 
2565 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2566   write(msg,'(1x,"vdw_df_kernel_value:",1x,"napts = ",i8)') napts
2567   call wrtout(std_out,msg,'COLL')
2568 #endif
2569 
2570   ABI_ALLOCATE(amesh,(napts))
2571   ABI_ALLOCATE(amesh_cos,(napts))
2572   ABI_ALLOCATE(amesh_sin,(napts))
2573   ABI_ALLOCATE(dphida,(napts))
2574   ABI_ALLOCATE(dphidb,(napts))
2575   ABI_ALLOCATE(nu1,(napts))
2576   ABI_ALLOCATE(nu2,(napts))
2577   ABI_ALLOCATE(ycoef,(3,napts))
2578   ABI_ALLOCATE(work,(napts))
2579   ABI_ALLOCATE(eint,(napts))
2580   ! Init amesh
2581   !FIXME: allow for other mesh types
2582   !Now intruducing metyp, setting mesh_cutoff to acut, removing
2583   !mesh_inc=damax and including mesh_ratio=
2584   call vdw_df_create_mesh(amesh,napts,metyp,acut,mesh_ratio=damax/dain)
2585 
2586   ! Init cos(a), sin(a), nu1(a), and nu2(a)
2587   tau = 4 * pi / 9
2588 !$OMP  PARALLEL DO &
2589 !$OMP& IF ( napts > 511 ) &
2590 !$OMP& DEFAULT(SHARED) PRIVATE(ia)
2591   do ia = 1,napts
2592     atmp = amesh(ia)
2593     amesh_cos(ia) = cos(atmp)
2594     amesh_sin(ia) = sin(atmp)
2595     if ( atmp == zero ) then  !it was changed the original condition: atmp < atol
2596       nu1(ia) = d1**2 / 2 / tau
2597       nu2(ia) = d2**2 / 2 / tau
2598     else
2599       if ( d1 == zero ) then  !it was changed the original condition: d1 < atol
2600         nu1(ia) = atmp*atmp / 2
2601       else
2602         nu1(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d1)**2))
2603       end if
2604       if ( d2 == zero ) then  !it was changed the original condition: d2 < atol
2605         nu2(ia) = atmp*atmp / 2
2606       else
2607         nu2(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d2)**2))
2608       end if
2609     end if
2610   end do
2611 !$OMP END PARALLEL DO
2612 
2613   ! Integrate
2614   dphida(1) = 0
2615 !$OMP PARALLEL DO &
2616 !$OMP& IF ( napts > 255 ) &
2617 !$OMP& DEFAULT(SHARED) PRIVATE(ia,ib,atmp,btmp,ttmp,wtmp)
2618   do ia = 2,napts
2619     atmp = amesh(ia)
2620     dphidb(1) = 0
2621 
2622     do ib = 2,napts
2623       btmp = amesh(ib)
2624 
2625       ! eq.(15) DRSLL04
2626       ttmp = half * ( 1 / (nu1(ia) + nu1(ib)) + 1 / (nu2(ia) + nu2(ib)) ) * &
2627 &       ( 1 / (nu1(ia) + nu2(ia)) / (nu1(ib) + nu2(ib)) + &
2628 &         1 / (nu1(ia) + nu2(ib)) / (nu2(ia) + nu1(ib)) )
2629 
2630       ! eq.(16) DRSLL04
2631       wtmp = two * ( (three - atmp*atmp) * btmp * amesh_cos(ib) * &
2632 &       amesh_sin(ia) + (three - btmp*btmp) * atmp * amesh_cos(ia) * &
2633 &       amesh_sin(ib) + (atmp*atmp + btmp*btmp - three) * &
2634 &       amesh_sin(ia) * amesh_sin(ib) - &
2635 &       three * atmp * btmp * amesh_cos(ia) * amesh_cos(ib) ) / &
2636 &       (atmp * btmp)**3
2637 
2638       dphidb(ib) = atmp*atmp * btmp*btmp * wtmp * ttmp
2639 
2640     end do ! ib = 2,napts
2641 
2642 !    call spline_integrate(dphida(ia),napts,damax,dphidb)
2643      call cspint(dphidb,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,dphida(ia))
2644   end do ! ia = 2,napts
2645 !$OMP END PARALLEL DO
2646 
2647   ! Final value of the kernel
2648 !  call spline_integrate(vdw_df_kernel_value,napts,damax,dphida)
2649      call cspint(dphida,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,vdw_df_kernel_value)
2650   vdw_df_kernel_value = vdw_df_kernel_value * 2 / pi**2
2651 
2652   ! Clean-up the mess
2653   ABI_DEALLOCATE(amesh)
2654   ABI_DEALLOCATE(amesh_cos)
2655   ABI_DEALLOCATE(amesh_sin)
2656   ABI_DEALLOCATE(nu1)
2657   ABI_DEALLOCATE(nu2)
2658   ABI_DEALLOCATE(dphida)
2659   ABI_DEALLOCATE(dphidb)
2660   ABI_DEALLOCATE(ycoef)
2661   ABI_DEALLOCATE(work)
2662   ABI_DEALLOCATE(eint)
2663 
2664 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2665   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
2666 #endif
2667 
2668 end function vdw_df_kernel_value

ABINIT/m_xc_vdw/vdw_df_ldaxc [ Functions ]

[ Top ] [ 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

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

2699 subroutine vdw_df_ldaxc(npts_rho,nspden,ngrad,rho_grho, &
2700 &                       exc_lda,vxc_lda,vxcg_lda)
2701 
2702 
2703 !This section has been created automatically by the script Abilint (TD).
2704 !Do not modify the following lines by hand.
2705 #undef ABI_FUNC
2706 #define ABI_FUNC 'vdw_df_ldaxc'
2707 !End of the abilint section
2708 
2709   implicit none
2710 
2711 !Arguments ------------------------------------
2712   integer,intent(in) :: nspden,npts_rho,ngrad
2713   real(dp),intent(in) :: rho_grho(npts_rho,nspden,ngrad)
2714   real(dp),intent(out) :: exc_lda(2,npts_rho),vxc_lda(2,npts_rho,nspden)
2715   real(dp),intent(out) :: vxcg_lda(2,3,npts_rho)
2716 
2717 !Local variables ------------------------------
2718   integer :: ii
2719   logical :: is_gga
2720   real(dp),allocatable :: rho_tmp(:,:),grho2(:,:)
2721 
2722 ! *************************************************************************
2723 
2724   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2725 
2726   is_gga=libxc_functionals_isgga(vdw_funcs)
2727 
2728   ABI_ALLOCATE(rho_tmp,(npts_rho,nspden))
2729   if (is_gga) then
2730     ABI_ALLOCATE(grho2,(npts_rho,2*min(nspden,2)-1))
2731   end if
2732 
2733   ! Convert the quantities provided by ABINIT to the ones needed by LibXC
2734   !
2735   ! Notes:
2736   !
2737   !   * Abinit passes rho_up in the spin-unpolarized case, while LibXC
2738   !     expects the total density, but as we use rhonow, we don't need
2739   !     to convert anything.
2740   !
2741   !   * In the case off gradients:
2742   !
2743   !       - Abinit passes |grho_up|^2 while LibXC needs |grho_tot|^2;
2744   !       - Abinit passes |grho_up|^2, |grho_dn|^2, and |grho_tot|^2,
2745   !         while LibXC needs |grho_up|^2, grho_up.grho_dn,
2746   !         and |grho_dn|^2;
2747   !
2748   !     but again the use of rho_grho lets us build these quantities
2749   !     directly.
2750   !
2751   ! See ~abinit/src/56_xc/rhotoxc.F90 for details.
2752   if (nspden==1) then
2753     do ii=1,npts_rho
2754       rho_tmp(ii,1)=half*rho_grho(ii,1,1)
2755       if (is_gga) grho2(ii,1)=quarter*(rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2 &
2756 &                                     +rho_grho(ii,1,4)**2)
2757     end do
2758   else
2759     do ii=1,npts_rho
2760       rho_tmp(ii,1)=rho_grho(ii,2,1)
2761       rho_tmp(ii,2)=rho_grho(ii,1,1)-rho_grho(ii,2,1)
2762       if (is_gga) then
2763         grho2(ii,1)=rho_grho(ii,2,2)**2+rho_grho(ii,2,3)**2+rho_grho(ii,2,4)**2
2764         grho2(ii,2)=(rho_grho(ii,1,2)-rho_grho(ii,2,2))**2 &
2765 &                  +(rho_grho(ii,1,3)-rho_grho(ii,2,3))**2 &
2766 &                  +(rho_grho(ii,1,4)-rho_grho(ii,2,4))**2
2767         grho2(ii,3)=rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2+rho_grho(ii,1,4)**2
2768       end if
2769     end do
2770   end if
2771 
2772   if (is_gga) then
2773     call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp,exc_lda,vxc_lda,&
2774 &                                 grho2=grho2,vxcgr=vxcg_lda,xc_functionals=vdw_funcs)
2775   else
2776     call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp,exc_lda,vxc_lda,&
2777 &                                 xc_functionals=vdw_funcs)
2778   end if
2779 
2780   ABI_DEALLOCATE(rho_tmp)
2781   if (is_gga) then
2782     ABI_DEALLOCATE(grho2)
2783   end if
2784 
2785   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
2786 
2787 end subroutine vdw_df_ldaxc

ABINIT/m_xc_vdw/vdw_df_set_tweaks [ Functions ]

[ Top ] [ Functions ]

NAME

  vdw_df_set_tweaks

FUNCTION

  Expands user-specified tweaks for internal use.

INPUTS

  input_tweaks= condensed tweaks

OUTPUT

  output_tweaks= expanded tweaks

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

3371 subroutine vdw_df_set_tweaks(input_tweaks,output_tweaks)
3372 
3373 
3374 !This section has been created automatically by the script Abilint (TD).
3375 !Do not modify the following lines by hand.
3376 #undef ABI_FUNC
3377 #define ABI_FUNC 'vdw_df_set_tweaks'
3378 !End of the abilint section
3379 
3380   implicit none
3381 
3382 !Arguments ------------------------------------
3383   integer,intent(in) :: input_tweaks
3384   type(vdw_df_tweaks_type),intent(out) :: output_tweaks
3385 
3386 !Local variables-------------------------------
3387 
3388 ! *************************************************************************
3389 
3390   output_tweaks%amesh_type = iand(ishft(input_tweaks,-10),3)
3391   output_tweaks%deriv_type = iand(ishft(input_tweaks,-2),3)
3392   output_tweaks%dmesh_type = iand(ishft(input_tweaks,-6),3)
3393   output_tweaks%interp_type = iand(ishft(input_tweaks,-4),3)
3394   output_tweaks%qmesh_type = iand(ishft(input_tweaks,-8),3)
3395   output_tweaks%run_type = iand(input_tweaks,3)
3396 
3397 end subroutine vdw_df_set_tweaks

ABINIT/m_xc_vdw/vdw_df_tweaks_type [ Types ]

[ Top ] [ Types ]

NAME

  vdw_df_tweaks_type

FUNCTION

  Tweaks for van der Waals XC calculations (use at your own risks).

SOURCE

159   type vdw_df_tweaks_type
160 
161     integer :: amesh_type
162     ! A-mesh type
163 
164     integer :: deriv_type
165     ! Derivation mode
166 
167     integer :: dmesh_type
168     ! D-mesh type
169 
170     integer :: interp_type
171     ! Interpolation mode
172 
173     integer :: qmesh_type
174     ! Q-mesh type
175 
176     integer :: run_type
177     ! Run mode
178 
179   end type vdw_df_tweaks_type

ABINIT/m_xc_vdw/vdw_df_write_func [ Functions ]

[ Top ] [ 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)

PARENTS

CHILDREN

      wrtout

SOURCE

3418 subroutine vdw_df_write_func(func_name,mode)
3419 
3420 
3421 !This section has been created automatically by the script Abilint (TD).
3422 !Do not modify the following lines by hand.
3423 #undef ABI_FUNC
3424 #define ABI_FUNC 'vdw_df_write_func'
3425 !End of the abilint section
3426 
3427   implicit none
3428 
3429 !Arguments ------------------------------------
3430   character(len=*),intent(in) :: func_name,mode
3431 
3432 !Local variables-------------------------------
3433   character(len=512) :: msg
3434 
3435 ! *************************************************************************
3436 
3437   write(msg,'(a,1x,a,1x,a,1x,a)') ch10,"=== subprogram :",func_name,"==="
3438   call wrtout(std_out,msg,mode)
3439 
3440 end subroutine vdw_df_write_func

ABINIT/m_xc_vdw/xc_vdw_aggregate [ Functions ]

[ Top ] [ 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.

PARENTS

      rhotoxc

CHILDREN

      wrtout

SOURCE

280 subroutine xc_vdw_aggregate(volume,gprimd,npts_rho,nspden,ngrad,nr1,nr2,nr3, &
281 & rho_grho,deltae_vdw,exc_vdw,decdrho_vdw,decdgrho_vdw,stress_vdw)
282 
283 
284 !This section has been created automatically by the script Abilint (TD).
285 !Do not modify the following lines by hand.
286 #undef ABI_FUNC
287 #define ABI_FUNC 'xc_vdw_aggregate'
288 !End of the abilint section
289 
290   implicit none
291 
292 #if defined HAVE_FFT_FFTW3
293   include "fftw3.f03"
294 #endif
295 
296 !Arguments ------------------------------------
297   integer,intent(in) :: npts_rho,nspden,ngrad,nr1,nr2,nr3
298   real(dp),intent(in) :: volume,gprimd(3,3)
299   real(dp),intent(in) :: rho_grho(npts_rho,nspden,ngrad)
300   real(dp),intent(out) :: exc_vdw,deltae_vdw
301   real(dp),intent(out) :: decdrho_vdw(nspden),decdgrho_vdw(3,nspden)
302   real(dp),intent(out),optional :: stress_vdw(3,3)
303 
304 !Local variables ------------------------------
305   character(len=512) :: msg
306   integer :: ig,ip1,ip2,iq1,iq2,ir1,ir2,ir3,is1,is2,ngpts,nqpts
307   integer(kind=SIZEOF_PTRDIFF_T) :: fftw3_plan
308   real(dp) :: a1,a2,a3,b1,b2,b3,gtmp,gcut,sg
309   real(dp) :: ex,ec,vx,vc
310   real(dp) :: exc_nl,eps_vdw,deltae_uns,dexc,dexcg(3)
311   real(dp) :: dg,dvol,rho_tmp,gtol
312   real(dp) :: exc_tmp,decdrho_tmp(nspden),decdgrho_tmp(3,nspden)
313   real(dp),allocatable :: exc_lda(:,:),vxc_lda(:,:),vxcg_lda(:,:,:)
314   real(dp),allocatable :: gvec(:,:),theta(:,:,:)
315   real(dp),allocatable :: t3dr(:,:,:,:)
316   complex(dp),allocatable :: t3dg(:,:,:,:)
317   real(dp),allocatable :: ptmp(:,:,:),ttmp(:,:),utmp(:),wtmp(:)
318 
319 ! *************************************************************************
320 
321   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
322 
323   ! Init
324   ngpts = my_vdw_params%ngpts
325   nqpts = my_vdw_params%nqpts
326   gcut = my_vdw_params%gcut
327   gtol = my_vdw_params%tolerance
328   exc_vdw = zero
329   deltae_uns = zero
330   deltae_vdw = zero
331   if ( present(stress_vdw) ) stress_vdw(:,:) = zero
332   decdgrho_tmp(:,:) = zero
333   decdrho_tmp(:) = zero
334   dvol = volume / npts_rho
335 
336   ! Check the status of the switch
337   if ( .not. vdw_switch ) return
338 
339   ABI_ALLOCATE(exc_lda,(2,npts_rho))
340   ABI_ALLOCATE(vxc_lda,(2,npts_rho))
341   ABI_ALLOCATE(vxcg_lda,(2,3,npts_rho))
342   ABI_ALLOCATE(gvec,(3,npts_rho))
343   ABI_ALLOCATE(theta,(nqpts,nspden,5))
344   ABI_ALLOCATE(t3dr,(nr1,nr2,nr3,nqpts))
345   ABI_ALLOCATE(t3dg,(nr1,nr2,nr3,nqpts))
346   ABI_ALLOCATE(ttmp,(nqpts,nr1*nr2*nr3))
347   ABI_ALLOCATE(ptmp,(nqpts,nqpts,2))
348   ABI_ALLOCATE(utmp,(nqpts))
349   ABI_ALLOCATE(wtmp,(nqpts))
350 
351   ! Calculate XC energy density from LDA / GGA
352   call vdw_df_ldaxc(npts_rho,nspden,ngrad,rho_grho,exc_lda,vxc_lda,vxcg_lda)
353 
354   ! Build theta in 3D and g-vectors
355   if ( npts_rho /= nr1*nr2*nr3 ) then
356     MSG_WARNING('The 3D reconstruction of the density might be wrong (npts /= nr1*nr2*nr3)')
357   end if
358 #if defined DEBUG_VERBOSE
359   write(msg,'(a,1x,i8,1x,a)') "Will now call xc_vdw_energy", &
360 &   nr1 * nr2 * nr3,"times"
361   MSG_COMMENT(msg)
362 #endif
363 
364   ip1 = 1
365   ip2 = 1
366   do ir3=1,nr3
367     do ir2=1,nr2
368       do ir1=1,nr1
369         gvec(:,ip1) = gprimd(:,1) * ir1 + gprimd(:,2) * ir2 + gprimd(:,3) * ir3
370         ex = exc_lda(1,ip1)
371         ec = exc_lda(2,ip1)
372         vx = vxc_lda(1,ip1)
373         vc = vxc_lda(2,ip1)
374         theta(:,:,:) = zero
375         call xc_vdw_energy(nspden,rho_grho(ip1,1:nspden,1), &
376 &         rho_grho(ip1,1:nspden,2:ngrad), &
377 &         ex,ec,vx,vc,theta,eps_vdw)
378         t3dr(ir1,ir2,ir3,1:nqpts) = theta(1:nqpts,1,1)
379         rho_tmp = sum(rho_grho(ip1,1:nspden,1))
380         deltae_uns = deltae_uns + rho_tmp * eps_vdw * dvol
381         ip1 = ip1 + 1
382       end do
383     end do
384 #if defined DEBUG_VERBOSE
385     if ( (ip1 - ip2) > 100 ) then
386       write(msg,'(1x,a,1x,i3,"% complete")') &
387 &       '[vdW-DF Energy]',int(ip1*100.0/(nr1*nr2*nr3))
388       call wrtout(std_out,msg,'COLL')
389       ip2 = ip1
390     end if
391 #endif
392   end do
393 
394   ! Fourier-transform theta
395 #if defined HAVE_FFT_FFTW3
396   do iq1=1,nqpts
397     call dfftw_plan_dft_r2c_3d(fftw3_plan,nr1,nr2,nr3, &
398 &     t3dr(:,:,:,iq1),t3dg(:,:,:,iq1),FFTW_ESTIMATE)
399     call dfftw_execute(fftw3_plan)
400     call dfftw_destroy_plan(fftw3_plan)
401     t3dg(:,:,:,iq1) = t3dg(:,:,:,iq1) / (nr1 * nr2 * nr3)
402   end do
403 #else
404   MSG_ERROR('vdW-DF calculations require FFTW3 support')
405 #endif
406 
407   ! Repack theta
408   ip1 = 1
409   do ir3=1,nr3
410     do ir2=1,nr2
411       do ir1=1,nr1
412         ttmp(:,ip1) = t3dg(ir1,ir2,ir3,:)
413         ip1 = ip1 + 1
414       end do
415     end do
416   end do
417 
418   ! Reset 3D counters, since theta will be reconstructed in 3D on the fly
419   ir1 = 1
420   ir2 = 1
421   ir3 = 1
422 
423   ! Go through reciprocal vectors
424   ! FIXME: get values of G-vectors
425   do ip1=1,npts_rho
426     if ( (ir3 > nr3) .and. (ir1 /= 1) .and. (ir2 /= 1) ) then
427       MSG_WARNING('buffer overflow reconstructing theta in 3D')
428     end if
429 
430     gtmp = sqrt(sum(gvec(:,ip1)**2))
431 
432     if ( gtmp < gcut ) then
433 
434       ! Interpolate phi in reciprocal space:
435       !   * ptmp(:,:,1) = phi(g)
436       !   * ptmp(:,:,2) = dphi(g)/dg
437       ! Note: do this on the fly, to go as fast as possible (this is equivalent
438       !       to a call of the 'splint' routine)
439       ptmp(:,:,:) = zero
440       dg = pi / my_vdw_params%rcut
441       ig = int(gtmp * ngpts / gcut)
442       a1 = ((ig + 1) * dg - gtmp) / dg
443       b1 = one - a1
444       a2 = (3 * a1**2 - one) * dg / six
445       b2 = (3 * b1**2 - one) * dg / six
446       a3 = (a1**3 - a1) * dg**2 / six
447       b3 = (b1**3 - b1) * dg**2 / six
448       do iq2 = 1,nqpts
449         do iq1 = 1,iq2
450           ptmp(iq1,iq2,1) = a1 * phig(ig,iq1,iq2) + b1 * phig(ig+1,iq1,iq2) &
451 &           + a3 * d2phidg2(ig,iq1,iq2) + b3 * d2phidg2(ig+1,iq1,iq2)
452           ptmp(iq1,iq2,2) = (phig(ig+1,iq1,iq2) - phig(ig,iq1,iq2)) / dg &
453 &           - a2 * d2phidg2(ig,iq1,iq2) + b2 * d2phidg2(ig+1,iq1,iq2)
454         end do
455       end do
456 
457       do iq2 = 1,nqpts
458         do iq1 = 1,iq2-1
459           ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:)
460         end do
461       end do
462 
463       ! Calculate contributions to integral in Fourier space:
464       ! FIXME: find back # of integral in paper
465       utmp(:) = matmul(ttmp(:,ip1),ptmp(:,:,1))
466 
467       ! Calculate contribution to stress in reciprocal space
468       ! Note: contribution of g=0 is zero
469       if ( present(stress_vdw) ) then
470         if ( gtmp > gtol ) then
471           wtmp(:) = matmul(ttmp(:,ip1),ptmp(:,:,2))
472           sg = sum(wtmp(:)*ttmp(:,ip1)) * volume / gtmp
473           do is2=1,3
474             do is1=1,3
475               stress_vdw(is1,is2) = stress_vdw(is1,is2) - &
476 &               sg * gvec(is1,ip1) * gvec(is2,ip1)
477             end do
478           end do
479         end if ! gtmp > gtol
480       end if ! present(stress_vdw)
481 
482     else
483 
484       utmp(:) = zero
485 
486     end if ! gtmp < gcut
487 
488     ! Reconstruct the integrand in 3D
489     t3dg(ir1,ir2,ir3,:) = utmp(:)
490 
491     ir1 = ir1 + 1
492     if ( ir1 > nr1 ) then
493       ir2 = ir2 + 1
494       ir1 = 1
495     end if
496     if ( ir2 > nr2 ) then
497       ir3 = ir3 + 1
498       ir2 = 1
499     end if
500   end do !ip1=1,npts_rho
501 
502   ! Fourier-transform back the integrand
503 #if defined HAVE_FFT_FFTW3
504   !write(msg,'(a)') ch10
505   !call wrtout(std_out,msg,'COLL')
506   do iq1=1,nqpts
507     !write(msg,'(1x,a,i4.4)') "xc_vdw_aggregate: backward FFT #",iq1
508     !call wrtout(std_out,msg,'COLL')
509 
510     call dfftw_plan_dft_c2r_3d(fftw3_plan,nr1,nr2,nr3, &
511 &     t3dg(:,:,:,iq1),t3dr(:,:,:,iq1),FFTW_ESTIMATE)
512     call dfftw_execute(fftw3_plan)
513     call dfftw_destroy_plan(fftw3_plan)
514   end do
515 #else
516   MSG_ERROR('vdW-DF calculations require FFTW3 support')
517 #endif
518 
519   ! Repack the integrand
520   ip1 = 1
521 !$OMP  PARALLEL DO COLLAPSE(3) &
522 !$OMP& DEFAULT(SHARED) PRIVATE(ir1,ir2,ir3)
523   do ir3=1,nr3
524     do ir2=1,nr2
525       do ir1=1,nr1
526         ttmp(:,ip1) = t3dr(ir1,ir2,ir3,1:nqpts)
527         ip1 = ip1 + 1
528       end do
529     end do
530   end do
531 !$OMP END PARALLEL DO
532 
533 #if defined DEBUG_VERBOSE
534   write(msg,'(a,1x,i8.8,1x,a)') "Will now call xc_vdw_energy",npts_rho,"times"
535   MSG_COMMENT(msg)
536 #endif
537 
538   ! Calculate and integrate vdW corrections at each point
539   do ip1=1,npts_rho
540 
541     ! Get local contributions
542     ex = exc_lda(1,ip1)
543     ec = exc_lda(2,ip1)
544     vx = vxc_lda(1,ip1)
545     vc = vxc_lda(2,ip1)
546     theta(:,:,:) = zero
547     call xc_vdw_energy(nspden,rho_grho(ip1,1:nspden,1), &
548 &     rho_grho(ip1,1:nspden,2:ngrad), &
549 &     ex,ec,vx,vc,theta,exc_tmp,decdrho_tmp,decdgrho_tmp)
550 
551     ! Get nonlocal contributons
552     ! Note: is2 represents cartesian coordinates here.
553     rho_tmp = sum(rho_grho(ip1,1:nspden,1))
554     deltae_vdw = deltae_vdw + rho_tmp * exc_tmp * dvol
555     exc_vdw = exc_vdw + rho_tmp * &
556 &     ( half * ( sum(ttmp(:,ip1) * theta(:,1,1)) / &
557 &       (rho_tmp + tiny(rho_tmp)) ) * dvol ) * dvol
558     !Correctness of multiplication by dvol above depends on how fftw3 deals
559     !with the volume element in direct or inverse FFT. Here it was assumed
560     !that fftw3 does not multiply by the volume upon FFT^-1
561     do is1=1,nspden
562       decdrho_vdw(is1) = decdrho_vdw(is1) + decdrho_tmp(is1) + &
563 &       sum(ttmp(:,ip1) * theta(:,is1,2))
564       do is2=1,3
565         decdgrho_vdw(is2,is1) = decdgrho_vdw(is2,is1) + decdgrho_tmp(is2,is1) + &
566 &         sum(ttmp(:,ip1) * theta(:,is1,is2+2))
567       end do
568     end do
569 
570 #if defined DEBUG_VERBOSE
571     if ( modulo(ip1,int(npts_rho/100)) == 0 ) then
572       write(msg,'(1x,a,1x,i3,"% complete")') &
573 &       "[vdW-DF Energy]",int(ip1*100.0/npts_rho)
574       call wrtout(std_out,msg,'COLL')
575     end if
576 #endif
577 
578   end do ! Loop on npts_rho
579 
580 
581   deltae_vdw = deltae_uns + deltae_vdw
582 
583 #if defined DEBUG_VERBOSE
584   write(msg,'(1x,a)') "[vdW-DF Enrgy] 100% complete"
585 #endif
586 
587   ! Display results
588   write(msg,'(a,1x,3a,2(3x,a,1x,e12.5,a),3x,a,1(1x,e12.5),a, &
589 &   3x,a,3(1x,e12.5),a)') &
590 &   ch10,"[vdW-DF] xc_vdw_aggregate: reporting vdW-DF contributions", &
591 &   ch10,ch10, &
592 &   "DeltaE_vdw = ",deltae_vdw,ch10, &
593 &   "Exc_vdw = ",exc_vdw,ch10, &
594 &   "dExc_vdw/drho = ",decdrho_vdw(:),ch10, &
595 &   "dExc_vdw/dgrho = ",decdgrho_vdw(:,:),ch10
596   call wrtout(std_out,msg,'COLL')
597 
598   ! Final adjustments of stress
599   if ( present(stress_vdw) ) then
600     stress_vdw(:,:) = stress_vdw(:,:) / volume
601 
602     write(msg,'(3x,a,a,3(5x,3(1x,e12.5),a))') &
603 &   "Stress_vdw = ",ch10, &
604 &      stress_vdw(1,:),ch10, &
605 &      stress_vdw(2,:),ch10, &
606 &      stress_vdw(3,:),ch10
607     call wrtout(std_out,msg,'COLL')
608   end if
609 
610   ! Clean-up the mess
611   ABI_DEALLOCATE(exc_lda)
612   ABI_DEALLOCATE(vxc_lda)
613   ABI_DEALLOCATE(vxcg_lda)
614   ABI_DEALLOCATE(gvec)
615   ABI_DEALLOCATE(theta)
616   ABI_DEALLOCATE(t3dr)
617   ABI_DEALLOCATE(t3dg)
618   ABI_DEALLOCATE(ptmp)
619   ABI_DEALLOCATE(ttmp)
620   ABI_DEALLOCATE(utmp)
621   ABI_DEALLOCATE(wtmp)
622 
623   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
624 
625 end subroutine xc_vdw_aggregate

ABINIT/m_xc_vdw/xc_vdw_done [ Functions ]

[ Top ] [ 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

PARENTS

      driver,vdw_kernelgen

CHILDREN

      wrtout

SOURCE

869 subroutine xc_vdw_done(vdw_params)
870 
871 
872 !This section has been created automatically by the script Abilint (TD).
873 !Do not modify the following lines by hand.
874 #undef ABI_FUNC
875 #define ABI_FUNC 'xc_vdw_done'
876 !End of the abilint section
877 
878   implicit none
879 
880 !Arguments ------------------------------------
881   type(xc_vdw_type), intent(in) :: vdw_params
882 
883 !Local variables ------------------------------
884   character(len=512) :: msg
885   integer :: i
886 
887 ! *************************************************************************
888 
889   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
890 
891   call libxc_functionals_end(xc_functionals=vdw_funcs)
892 
893   ABI_DEALLOCATE(dmesh)
894   ABI_DEALLOCATE(qmesh)
895   ABI_DEALLOCATE(qpoly_basis)
896   ABI_DEALLOCATE(phi)
897   ABI_DEALLOCATE(phi_bicubic)
898   ABI_DEALLOCATE(phi_u_bicubic)
899   ABI_DEALLOCATE(phir)
900   ABI_DEALLOCATE(phir_u)
901   ABI_DEALLOCATE(d2phidr2)
902   ABI_DEALLOCATE(phig)
903   ABI_DEALLOCATE(d2phidg2)
904 
905   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
906 
907 end subroutine xc_vdw_done

ABINIT/m_xc_vdw/xc_vdw_energy [ Functions ]

[ Top ] [ 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.

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

663 subroutine xc_vdw_energy(nspden,rho,grho,ex_lda,ec_lda,vx_lda,vc_lda, &
664 & theta,eps,dexc,dexcg)
665 
666 
667 !This section has been created automatically by the script Abilint (TD).
668 !Do not modify the following lines by hand.
669 #undef ABI_FUNC
670 #define ABI_FUNC 'xc_vdw_energy'
671 !End of the abilint section
672 
673   implicit none
674 
675 !Arguments ------------------------------------
676   integer,intent(in) :: nspden
677   real(dp),intent(in) :: ex_lda,ec_lda,vx_lda,vc_lda
678   real(dp),intent(in) :: rho(nspden),grho(nspden,3)
679   real(dp),intent(inout) :: theta(my_vdw_params%nqpts,nspden,5)
680   real(dp),intent(out),optional :: eps,dexc(2),dexcg(3,nspden)
681 
682 !Local variables ------------------------------
683   character(len=512) :: msg
684   logical :: calc_corrections
685   integer :: iq0,iq1,iq2,ir,is,ix,nqpts,nrpts,ns
686   real(dp) :: decdrho,dexdrho,dvcdrho,dvxdrho,ec,ex,vc,vx
687   real(dp) :: dr,kappa,rho_tmp,grho_tmp(3),grho2_tmp,rtol,qcut,zab
688   real(dp) :: dq(3),ptmp(2),q0(5),qtmp(2)
689   real(dp) :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts)
690   real(dp) :: qpoly(my_vdw_params%nqpts,3)
691 
692 ! *************************************************************************
693 
694   ! Init
695   nqpts = my_vdw_params%nqpts
696   nrpts = my_vdw_params%nrpts
697   ns = my_vdw_params%nsmooth
698   qcut = my_vdw_params%qcut
699   rtol = my_vdw_params%tolerance
700   zab = my_vdw_params%zab
701 
702   calc_corrections = .false.
703   if ( present(eps) .and. present(dexc) .and. present(dexcg) ) then
704     calc_corrections = .true.
705   end if
706 
707   ! Sum density over spin
708   rho_tmp = sum(rho(1:nspden))
709   kappa = (three * pi**2 * rho_tmp)**third
710   forall(ix=1:3) grho_tmp(ix) = sum(grho(1:nspden,ix))
711   grho2_tmp = sum(grho_tmp**2)
712 
713   ! Calculate local wavevector q0 of eqs. 11-12 from DRSLL04
714   !   * q0(1)   = q0
715   !   * q0(2)   = dq0 / drho
716   !   * q0(3:5) = dq0 / dgrho2
717   ! Notes:
718   !   * treating rho->0 separately (divide by 0)
719   !   * treating ex_lda->0 separately (divide by 0)
720   q0(:) = zero
721   if ( (rho_tmp < rtol) .or. (ex_lda < rtol) ) then
722     q0(1) = qcut
723   else
724     q0(1) = kappa * (one + ec_lda / ex_lda - zab / nine * grho2_tmp / &
725 &     (two * kappa * rho_tmp)**2)
726 
727     if ( calc_corrections ) then
728      q0(2) = ( (vc_lda - ec_lda) / (rho_tmp * ex_lda) - ec_lda * &
729 &      (vx_lda - ex_lda) / &
730 &      (rho_tmp * ex_lda**2) + two * zab / nine * grho2_tmp / &
731 &      (two * kappa * rho_tmp)**3 * &
732 &      eight * kappa / three ) * kappa + q0(1) / three * rho_tmp
733      q0(3:5) = -two * (zab / nine) / (two * kappa * rho_tmp)**2 * &
734 &      grho_tmp(:)
735     end if
736   end if
737 
738   if ( q0(1) > qcut ) then
739     q0(1) = qcut
740     q0(2:5)= zero
741   end if
742 
743   ! Smoothen q0 near qcut exponentially  eq. 5 RS09. (Saturation)
744   !   * qtmp(1) = qs = qcut * (1 - exp(-Sum_i=1,ns 1/i (x/xc)**i))
745   !   * qtmp(2) = dqs / dq |q=q0
746   qtmp(1) = (q0(1) / qcut) / ns
747   qtmp(2) = one / qcut
748   do is=ns-1,1,-1
749     qtmp(1) = (qtmp(1) + one / is) * q0(1) / qcut
750     qtmp(2) = (qtmp(2) * q0(1) + one) / qcut
751   end do
752   qtmp(2) = qtmp(2) * qcut * exp(-qtmp(1))
753 
754   q0(1) = qcut * (one - exp(-qtmp(1)))
755   q0(2:5) = q0(2:5) * qtmp(2)
756 
757   ! Calculate polynomial coefficients for cubic-spline interpolation at q0
758   !   * qpoly(:,1) = coefficients
759   !   * qpoly(:,2) = first derivatives
760   qpoly(:,:) = zero
761   iq0 = vdw_df_indexof(qmesh(:),nqpts,q0(1))
762   dq(1) = qmesh(iq0+1) - qmesh(iq0)
763   dq(2) = (qmesh(iq0+1) - q0(1)) / dq(1)
764   dq(3) = (q0(1) - qmesh(iq0)) / dq(1)
765   do iq1=1,nqpts
766    qpoly(iq1,1) = dq(2) * qpoly_basis(iq0,iq1,1) + dq(3) * &
767 &    qpoly_basis(iq0+1,iq1,1) + &
768 &    ((dq(2)**3 - dq(2)) * qpoly_basis(iq0,iq1,2) + &
769 &    (dq(3)**3 - dq(3)) * qpoly_basis(iq0+1,iq1,2)) * dq(1)**2 / six
770 
771    if ( calc_corrections ) then
772     qpoly(iq1,2) = -(qpoly_basis(iq0,iq1,1) - &
773 &     qpoly_basis(iq0+1,iq1,1)) / dq(1) - &
774 &     ((three * dq(2)**2 - one) * qpoly_basis(iq0,iq1,2) - &
775 &     (three * dq(3)**2 - one) * qpoly_basis(iq0+1,iq1,2)) * dq(1) / six
776    end if
777   end do
778 
779   ! Pre-compute integrand for vdW energy correction
780   ! by cubic-spline interpolation and store it in
781   ! qpoly(:,3). This is done twice: one for the
782   ! unsoftened kernel and other one for the
783   ! softened kernel
784   dr = my_vdw_params%rcut / nrpts
785   ztmp(:,:) = zero
786 
787   if ( calc_corrections ) then
788    do ir=1,nrpts
789      do iq1=1,nqpts
790        do iq2=1,iq1
791          ztmp(iq1,iq2) = ztmp(iq1,iq2) - two * pi * phir(ir,iq1,iq2) * dr
792        end do
793      end do
794    end do
795   else
796    do ir=1,nrpts
797      do iq1=1,nqpts
798        do iq2=1,iq1
799          ztmp(iq1,iq2) = ztmp(iq1,iq2) + two * pi * phir_u(ir,iq1,iq2) * dr
800        end do
801      end do
802    end do
803   end if ! calc_corrections
804 
805 
806   do iq1=1,nqpts
807     do iq2=1,iq1-1
808       ztmp(iq2,iq1) = ztmp(iq1,iq2)
809     end do
810   end do
811   qpoly(:,3) = matmul(qpoly(:,1),ztmp(:,:))
812 
813   ! Calculate theta and its derivatives (see RS09)
814   !   * theta(:,:,1)   = theta
815   !   * theta(:,:,2)   = dtheta / drho
816   !   * theta(:,:,3:5) = dtheta / dgrho2
817   ! Note: trick to go from Abinit to LibXC conventions
818 
819   do is=1,nspden
820     theta(:,is,1) = qpoly(:,1) * (2 - nspden) * rho(is)
821   end do
822 
823   ! Calculate theta derivatives
824   if ( calc_corrections ) then
825     do is=1,nspden
826       theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(is)
827       do ix=3,5
828         theta(:,is,ix) = qpoly(:,2) * q0(ix) * (2 - nspden) * rho(is)
829       end do
830     end do
831   end if
832 
833   ! Calculate energy corrections
834     eps = zero
835     ptmp(1) = sum(qpoly(:,3) * qpoly(:,1))
836     eps = ptmp(1) * rho_tmp
837   if ( calc_corrections ) then
838      dexc(:) = zero
839      dexcg(:,:) = zero
840      ptmp(2) = two * sum(qpoly(:,2) * qpoly(:,3))
841      dexc(:) = two * ptmp(1) * rho_tmp + ptmp(2) * qtmp(2) * rho_tmp**2
842     do ix=3,5
843       dexcg(ix-2,:) = ptmp(2) * q0(ix) * rho_tmp**2
844     end do
845   end if ! calc_corrections
846 
847 end subroutine xc_vdw_energy

ABINIT/m_xc_vdw/xc_vdw_get_params [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_get_params

FUNCTION

  Exports internal vdW-DF parameters.

OUTPUT

  vdw_params= van der Waals parameters

PARENTS

      vdw_kernelgen

CHILDREN

      wrtout

SOURCE

928 subroutine xc_vdw_get_params(vdw_params)
929 
930 
931 !This section has been created automatically by the script Abilint (TD).
932 !Do not modify the following lines by hand.
933 #undef ABI_FUNC
934 #define ABI_FUNC 'xc_vdw_get_params'
935 !End of the abilint section
936 
937   implicit none
938 
939 !Arguments ------------------------------------
940   type(xc_vdw_type), intent(out)  :: vdw_params
941 
942 !Local variables ------------------------------
943   character(len=512) :: msg
944 
945 ! *************************************************************************
946 
947   vdw_params%functional = my_vdw_params%functional
948   vdw_params%zab        = my_vdw_params%zab
949   vdw_params%ndpts      = my_vdw_params%ndpts
950   vdw_params%dcut       = my_vdw_params%dcut
951   vdw_params%dratio     = my_vdw_params%dratio
952   vdw_params%dsoft      = my_vdw_params%dsoft
953   vdw_params%phisoft    = my_vdw_params%phisoft
954   vdw_params%nqpts      = my_vdw_params%nqpts
955   vdw_params%qcut       = my_vdw_params%qcut
956   vdw_params%qratio     = my_vdw_params%qratio
957   vdw_params%nrpts      = my_vdw_params%nrpts
958   vdw_params%rcut       = my_vdw_params%rcut
959   vdw_params%rsoft      = my_vdw_params%rsoft
960   vdw_params%ngpts      = my_vdw_params%ngpts
961   vdw_params%gcut       = my_vdw_params%gcut
962   vdw_params%acutmin    = my_vdw_params%acutmin
963   vdw_params%aratio     = my_vdw_params%aratio
964   vdw_params%damax      = my_vdw_params%damax
965   vdw_params%damin      = my_vdw_params%damin
966   vdw_params%nsmooth    = my_vdw_params%nsmooth
967   vdw_params%tolerance  = my_vdw_params%tolerance
968   vdw_params%tweaks     = my_vdw_params%tweaks
969 
970 end subroutine xc_vdw_get_params

ABINIT/m_xc_vdw/xc_vdw_init [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_init

FUNCTION

  Calculates the van der Waals kernel.

INPUTS

  vdw_params= parameters for the van der Waals calculations

PARENTS

      driver,vdw_kernelgen

CHILDREN

      wrtout

SOURCE

 991 subroutine xc_vdw_init(vdw_params)
 992 
 993 
 994 !This section has been created automatically by the script Abilint (TD).
 995 !Do not modify the following lines by hand.
 996 #undef ABI_FUNC
 997 #define ABI_FUNC 'xc_vdw_init'
 998 !End of the abilint section
 999 
1000   implicit none
1001 
1002 !Arguments ------------------------------------
1003   type(xc_vdw_type), intent(in)  :: vdw_params
1004 
1005 !Local variables ------------------------------
1006   character(len=*),parameter :: dmesh_file = "vdw_df_dmesh.pts"
1007   character(len=*),parameter :: qmesh_file = "vdw_df_qmesh.pts"
1008 
1009   character(len=512) :: msg
1010   integer :: id1,id2,iq1,iq2,ir,ir0,ixc,jd1,jd2,ndpts,ngpts,nqpts,nrpts,ntest,sofswt
1011   real(dp) :: d1,d2,dcut,dd,delta_test,dsoft,phisoft,acutmin,aratio,damax
1012   real(dp) :: phimm,phimp,phipm,phipp,phipn,phi_tmp,ptmp(2)
1013   real(dp),allocatable :: deltd(:),kern_spline_der(:,:,:)
1014   real(dp),allocatable :: btmp(:,:),utmp(:),xde(:)
1015 
1016 ! *************************************************************************
1017 
1018   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
1019 
1020   ! Global init
1021   call xc_vdw_set_params(vdw_params)
1022   call vdw_df_set_tweaks(my_vdw_params%tweaks,my_vdw_tweaks)
1023 
1024   ! Local init
1025   ndpts   = my_vdw_params%ndpts
1026   nqpts   = my_vdw_params%nqpts
1027   nrpts   = my_vdw_params%nrpts
1028   acutmin = my_vdw_params%acutmin
1029   aratio  = my_vdw_params%aratio
1030   damax   = my_vdw_params%damax
1031   dsoft   = my_vdw_params%dsoft
1032   phisoft = my_vdw_params%phisoft
1033 
1034   ! Create d-mesh
1035   ! Note: avoid zero and make sure the last point is exactly dcut
1036   ABI_ALLOCATE(dmesh,(ndpts))
1037 #if defined DEBUG_VERBOSE
1038   write(msg,'(1x,a)') "[vdW-DF Build] Generating D-mesh"
1039   call wrtout(std_out,msg,'COLL')
1040 #endif
1041   call vdw_df_create_mesh(dmesh,ndpts,my_vdw_tweaks%dmesh_type, &
1042 &   my_vdw_params%dcut,mesh_ratio=my_vdw_params%dratio, &
1043 &   mesh_tolerance=my_vdw_params%tolerance,mesh_file=dmesh_file, &
1044 &   avoid_zero=.true.,exact_max=.true.)
1045 
1046   ! Create q-mesh
1047   ABI_ALLOCATE(qmesh,(nqpts))
1048 #if defined DEBUG_VERBOSE
1049   write(msg,'(1x,a)') "[vdW-DF Build] Generating Q-mesh"
1050   call wrtout(std_out,msg,'COLL')
1051 #endif
1052   call vdw_df_create_mesh(qmesh,nqpts,my_vdw_tweaks%qmesh_type, &
1053 &   my_vdw_params%qcut,mesh_ratio=my_vdw_params%qratio, &
1054 &   mesh_tolerance=my_vdw_params%tolerance, mesh_file=qmesh_file)
1055 
1056   ! Build polynomial basis for cubic-spline interpolation
1057   !   * qpoly_basis(:,1) = polynomial basis
1058   !   * qpoly_basis(:,2) = second derivative
1059   ABI_ALLOCATE(qpoly_basis,(nqpts,nqpts,2))
1060   ABI_ALLOCATE(btmp,(nqpts,nqpts))
1061   ABI_ALLOCATE(utmp,(nqpts))
1062 #if defined DEBUG_VERBOSE
1063   write(msg,'(1x,a)') "[vdW-DF Build] Generating polynomial basis"
1064   call wrtout(std_out,msg,'COLL')
1065 #endif
1066   btmp(:,:) = zero
1067   forall(iq1=1:nqpts) btmp(iq1,iq1) = one
1068   utmp(:) = zero
1069   qpoly_basis(:,:,:) = zero
1070   do iq1=1,nqpts
1071     qpoly_basis(iq1,iq1,1) = one
1072     qpoly_basis(iq1,1,2) = zero
1073     utmp(1) = zero
1074     do iq2=2,nqpts-1
1075       ptmp(1) = (qmesh(iq2) - qmesh(iq2-1)) / &
1076 &     (qmesh(iq2+1) - qmesh(iq2-1))
1077       ptmp(2) = ptmp(1) * qpoly_basis(iq1,iq2-1,2) + two
1078       qpoly_basis(iq1,iq2,2) = (ptmp(1) - one) / ptmp(2)
1079       utmp(iq2) = ( six * ((btmp(iq1,iq2+1) - btmp(iq1,iq2)) / &
1080 &     (qmesh(iq2+1) - qmesh(iq2)) - (btmp(iq1,iq2) - &
1081 &     btmp(iq1,iq2-1)) / (qmesh(iq2) - qmesh(iq2-1))) / &
1082 &     (qmesh(iq2+1) - qmesh(iq2-1)) - ptmp(1) * utmp(iq2-1)) / &
1083 &     ptmp(2)
1084     end do
1085     utmp(nqpts) = zero
1086     qpoly_basis(iq1,nqpts,2) = zero
1087     do iq2=nqpts-1,1,-1
1088       qpoly_basis(iq1,iq2,2) = qpoly_basis(iq1,iq2,2) * &
1089 &       qpoly_basis(iq1,iq2+1,2) + utmp(iq2)
1090     end do
1091   end do
1092   ABI_DEALLOCATE(btmp)
1093   ABI_DEALLOCATE(utmp)
1094 
1095   ! Create kernel and its derivatives
1096   ! Note: using 4 neighbours for derivatives
1097   ABI_ALLOCATE(phi,(ndpts,ndpts,4))
1098 #if defined DEBUG_VERBOSE
1099   write(msg,'(1x,a)') "[vdW-DF Build] Building kernel and its derivatives"
1100   call wrtout(std_out,msg,'COLL')
1101 #endif
1102 
1103   if ( phisoft < zero ) then
1104     phisoft = three * half * exp(-dsoft * two / sqrt(two))
1105     my_vdw_params%phisoft = phisoft
1106   end if
1107 
1108   ! Building of the softened kernel
1109 
1110   sofswt = 1
1111   do id1=1,ndpts
1112     d1 = dmesh(id1)
1113     phi(id1,id1,1) = vdw_df_kernel(d1,d1,dsoft,phisoft,acutmin,aratio,damax)
1114 
1115     ! Delta(d) should be at least 10^-6
1116     ! WARNING: changed tol3 to tol6 and changed max to min
1117     dd = min(my_vdw_params%dratio*d1,tol6)
1118 
1119       phimm = vdw_df_kernel(d1-dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax)
1120       phipm = vdw_df_kernel(d1+dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax)
1121       phipp = vdw_df_kernel(d1+dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax)
1122       phimp = vdw_df_kernel(d1-dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax)
1123 
1124       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1125       phi(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1126       phi(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1127       phi(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1128 
1129     do id2=1,id1-1
1130       d2 = dmesh(id2)
1131       phi(id1,id2,1) = vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax)
1132       phi(id2,id1,1) = phi(id1,id2,1)
1133 
1134       phimm = vdw_df_kernel(d1-dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax)
1135       phipm = vdw_df_kernel(d1+dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax)
1136       phipp = vdw_df_kernel(d1+dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax)
1137       phimp = vdw_df_kernel(d1-dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax)
1138 
1139       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1140       phi(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1141       phi(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1142       phi(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1143 
1144       phi(id2,id1,2) = phi(id1,id2,2)
1145       phi(id2,id1,3) = phi(id1,id2,3)
1146       phi(id2,id1,4) = phi(id1,id2,4)
1147     end do
1148 
1149     write(msg,'(1x,a,1x,i3,"% complete")') &
1150 &     '[vdW-DF Build]',int(id1*100.0/ndpts)
1151     call wrtout(std_out,msg,'COLL')
1152 #if defined DEBUG_VERBOSE
1153     call flush_unit(std_out)
1154 #endif
1155   end do
1156   write(msg,'(a)') " "
1157   call wrtout(std_out,msg,'COLL')
1158 
1159   ! Set boundary conditions
1160   ! Note: will have to check that the borders are smooth enough
1161   ! These boundary conditions produce large discontinuities, now testing its removal
1162   ! phi(ndpts,:,:) = zero
1163   ! phi(:,ndpts,:) = zero
1164   ! phi(1,:,2:4)   = zero
1165   ! phi(:,1,2:4)   = zero
1166   ! Ar2 results show better convergence of E_vdW than when boundary conditions were used.
1167   ! Kernel plot show that there are not dicontinuities as well.
1168 
1169 #if defined DEBUG_VERBOSE
1170   write(msg,'(1x,a)') "[vdW-DF Build] Building filtered kernel"
1171   call wrtout(std_out,msg,'COLL')
1172 #endif
1173 
1174   ! Calculate coefficients for bicubic interpolation
1175   ABI_ALLOCATE(phi_bicubic,(4,4,ndpts,ndpts))
1176   call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi(:,:,1),phi(:,:,2),&
1177 &   phi(:,:,3),phi(:,:,4),phi_bicubic)
1178 
1179   ! Build filtered kernel
1180   ABI_ALLOCATE(phir,(nrpts,nqpts,nqpts))
1181   ABI_ALLOCATE(d2phidr2,(nrpts,nqpts,nqpts))
1182   ABI_ALLOCATE(phig,(nrpts,nqpts,nqpts))
1183   ABI_ALLOCATE(d2phidg2,(nrpts,nqpts,nqpts))
1184   call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt)
1185   my_vdw_params%ngpts = ngpts
1186 
1187   ! Find closest indices in dmesh
1188   ! FIXME: something is missing or should be removed here
1189   ABI_ALLOCATE(kern_spline_der,(ndpts,ndpts,3))
1190   ABI_ALLOCATE(deltd,(2))
1191   ABI_ALLOCATE(xde,(2))
1192 
1193   kern_spline_der(:,:,:) = zero
1194 
1195   do jd1=1,ndpts
1196     do jd2=1,ndpts
1197 
1198       d1 = dmesh(jd1)
1199       d2 = dmesh(jd2)
1200 
1201       if ( (d1 < dcut) .or. (d2 < dcut) ) then
1202         id1 = vdw_df_indexof(dmesh,ndpts,d1)
1203         id2 = vdw_df_indexof(dmesh,ndpts,d2)
1204 
1205         deltd(1) = dmesh(id1+1) - dmesh(id1)
1206         deltd(2) = dmesh(id2+1) - dmesh(id2)
1207 
1208         xde(1) = (d1 - dmesh(id1)) / deltd(1)
1209         xde(2) = (d2 - dmesh(id2)) / deltd(2)
1210 
1211         kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1)
1212         kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2)
1213         kern_spline_der(jd1,jd2,3) = &
1214 &         phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) )
1215       end if
1216 
1217     end do
1218   end do
1219 
1220   ABI_DEALLOCATE(kern_spline_der)
1221   ABI_DEALLOCATE(deltd)
1222   ABI_DEALLOCATE(xde)
1223 
1224   ! Building of the unsoftened kernel
1225 
1226   ! Create unsoftened kernel and its derivatives
1227   ! Note: using 4 neighbours for derivatives
1228   ABI_ALLOCATE(phi_u,(ndpts,ndpts,4))
1229 #if defined DEBUG_VERBOSE
1230   write(msg,'(1x,a)') "[vdW-DF Build] Building unsoftened kernel and its derivatives"
1231   call wrtout(std_out,msg,'COLL')
1232 #endif
1233 
1234   do id1=1,ndpts
1235     d1 = dmesh(id1)
1236     phi_u(id1,id1,1) = vdw_df_kernel_value(d1,d1,acutmin,aratio,damax)
1237 
1238     ! Delta(d) should be at least 10^-6
1239     ! WARNING: changed tol3 to tol6 and changed max to min
1240     dd = min(my_vdw_params%dratio*d1,tol6)
1241 
1242       phimm = vdw_df_kernel_value(d1-dd,d1-dd,acutmin,aratio,damax)
1243       phipm = vdw_df_kernel_value(d1+dd,d1-dd,acutmin,aratio,damax)
1244       phipp = vdw_df_kernel_value(d1+dd,d1+dd,acutmin,aratio,damax)
1245       phimp = vdw_df_kernel_value(d1-dd,d1+dd,acutmin,aratio,damax)
1246 
1247       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1248       phi_u(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1249       phi_u(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1250       phi_u(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1251 
1252     do id2=1,id1-1
1253       d2 = dmesh(id2)
1254       phi_u(id1,id2,1) = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)
1255       phi_u(id2,id1,1) = phi_u(id1,id2,1)
1256 
1257       phimm = vdw_df_kernel_value(d1-dd,d2-dd,acutmin,aratio,damax)
1258       phipm = vdw_df_kernel_value(d1+dd,d2-dd,acutmin,aratio,damax)
1259       phipp = vdw_df_kernel_value(d1+dd,d2+dd,acutmin,aratio,damax)
1260       phimp = vdw_df_kernel_value(d1-dd,d2+dd,acutmin,aratio,damax)
1261 
1262       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1263       phi_u(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1264       phi_u(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1265       phi_u(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1266 
1267       phi_u(id2,id1,2) = phi_u(id1,id2,2)
1268       phi_u(id2,id1,3) = phi_u(id1,id2,3)
1269       phi_u(id2,id1,4) = phi_u(id1,id2,4)
1270     end do
1271 
1272     write(msg,'(1x,a,1x,i3,"% complete")') &
1273 &     '[vdW-DF-Uns Build]',int(id1*100.0/ndpts)
1274     call wrtout(std_out,msg,'COLL')
1275 #if defined DEBUG_VERBOSE
1276     call flush_unit(std_out)
1277 #endif
1278   end do
1279   write(msg,'(a)') " "
1280   call wrtout(std_out,msg,'COLL')
1281 
1282   ! Set boundary conditions
1283   ! Note: will have to check that the borders are smooth enough
1284   ! These boundary conditions produce large discontinuities, now testing its removal
1285   ! phi(ndpts,:,:) = zero
1286   ! phi(:,ndpts,:) = zero
1287   ! phi(1,:,2:4)   = zero
1288   ! phi(:,1,2:4)   = zero
1289   ! Ar2 results show better convergence of E_vdW than when boundary conditions were used.
1290   ! Kernel plot show that there are not dicontinuities as well.
1291 
1292 #if defined DEBUG_VERBOSE
1293   write(msg,'(1x,a)') "[vdW-DF Build] Building filtered kernel"
1294   call wrtout(std_out,msg,'COLL')
1295 #endif
1296 
1297   ! Calculate coefficients for bicubic interpolation
1298   ABI_ALLOCATE(phi_u_bicubic,(4,4,ndpts,ndpts))
1299   call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi_u(:,:,1),phi_u(:,:,2),&
1300 &   phi_u(:,:,3),phi_u(:,:,4),phi_u_bicubic)
1301 
1302   ! Build filtered kernel
1303   ABI_ALLOCATE(phir_u,(nrpts,nqpts,nqpts))
1304   ! For the local correction we just need the kernel not its
1305   ! derivatives.
1306   sofswt = 0
1307   !ABI_ALLOCATE(d2phidr2,(nrpts,nqpts,nqpts))
1308   !ABI_ALLOCATE(phig,(nrpts,nqpts,nqpts))
1309   !ABI_ALLOCATE(d2phidg2,(nrpts,nqpts,nqpts))
1310   call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt)
1311   !  my_vdw_params%ngpts = ngpts already defined above
1312 
1313   ! The following is not needed for the local correction:
1314   ! Find closest indices in dmesh
1315   ! FIXME: something is missing or should be removed here
1316   ! ABI_ALLOCATE(kern_spline_der,(ndpts,ndpts,3))
1317   ! ABI_ALLOCATE(deltd,(2))
1318   ! ABI_ALLOCATE(xde,(2))
1319 
1320   ! kern_spline_der(:,:,:) = zero
1321 
1322   ! do jd1=1,ndpts
1323   !   do jd2=1,ndpts
1324 
1325   !     d1 = dmesh(jd1)
1326   !     d2 = dmesh(jd2)
1327 
1328   !     if ( (d1 < dcut) .or. (d2 < dcut) ) then
1329   !       id1 = vdw_df_indexof(dmesh,ndpts,d1)
1330   !       id2 = vdw_df_indexof(dmesh,ndpts,d2)
1331 
1332   !       deltd(1) = dmesh(id1+1) - dmesh(id1)
1333   !       deltd(2) = dmesh(id2+1) - dmesh(id2)
1334 
1335   !       xde(1) = (d1 - dmesh(id1)) / deltd(1)
1336   !       xde(2) = (d2 - dmesh(id2)) / deltd(2)
1337 
1338   !       kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1)
1339   !       kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2)
1340   !       kern_spline_der(jd1,jd2,3) = &
1341   !&         phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) )
1342   !     end if
1343 
1344   !   end do
1345   ! end do
1346 
1347   ! ABI_DEALLOCATE(kern_spline_der)
1348   ! ABI_DEALLOCATE(deltd)
1349   ! ABI_DEALLOCATE(xde)
1350 
1351   call vdw_df_internal_checks(1)
1352 
1353   vdw_switch = .true.
1354 
1355   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
1356 
1357 end subroutine xc_vdw_init

ABINIT/m_xc_vdw/xc_vdw_libxc_init [ Functions ]

[ Top ] [ 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.

PARENTS

      driver

CHILDREN

      wrtout

SOURCE

1386 subroutine xc_vdw_libxc_init(ixc_vdw)
1387 
1388 
1389 !This section has been created automatically by the script Abilint (TD).
1390 !Do not modify the following lines by hand.
1391 #undef ABI_FUNC
1392 #define ABI_FUNC 'xc_vdw_libxc_init'
1393 !End of the abilint section
1394 
1395   implicit none
1396 
1397 !Arguments ------------------------------------
1398   integer,intent(in) :: ixc_vdw
1399 
1400 !Local variables ------------------------------
1401   character(len=*),parameter :: c_vdw1 = "GGA_C_PBE"
1402   character(len=*),parameter :: x_vdw1 = "GGA_X_PBE_R"
1403   character(len=*),parameter :: x_vdw2 = "GGA_X_PW86"
1404   character(len=*),parameter :: x_vdw3 = "GGA_X_C09X"
1405   integer :: i,ii,ic,ix,ixc
1406   character(len=500) :: msg
1407 
1408 ! *************************************************************************
1409 
1410   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
1411 
1412   ! Select LibXC functionals
1413   select case (ixc_vdw)
1414     case (1)
1415       ix = libxc_functionals_getid(x_vdw1)
1416       ic = libxc_functionals_getid(c_vdw1)
1417     case (2)
1418       ix = libxc_functionals_getid(x_vdw2)
1419       ic = libxc_functionals_getid(c_vdw1)
1420     case (3)
1421       ix = libxc_functionals_getid(x_vdw3)
1422       ic = libxc_functionals_getid(c_vdw1)
1423       MSG_ERROR('[vdW-DF] C09 not available for now')
1424     case default
1425       MSG_ERROR('[vdW-DF] invalid setting of vdw_xc')
1426   end select
1427   if ( (ix == -1) .or. (ic == -1) ) then
1428     MSG_ERROR('[vdW-DF] unable to set LibXC parameters')
1429   end if
1430   ixc = -(ix * 1000 + ic)
1431 
1432   ! Propagate to internal parameters
1433   my_vdw_params%functional = ixc_vdw
1434 
1435   ! XC functional init
1436   call libxc_functionals_init(ixc,1,xc_functionals=vdw_funcs)
1437 
1438   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
1439 
1440 end subroutine xc_vdw_libxc_init

ABINIT/m_xc_vdw/xc_vdw_memcheck [ Functions ]

[ Top ] [ 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

PARENTS

      driver,vdw_kernelgen

CHILDREN

      wrtout

SOURCE

1462 subroutine xc_vdw_memcheck(unt,vp)
1463 
1464 
1465 !This section has been created automatically by the script Abilint (TD).
1466 !Do not modify the following lines by hand.
1467 #undef ABI_FUNC
1468 #define ABI_FUNC 'xc_vdw_memcheck'
1469 !End of the abilint section
1470 
1471   implicit none
1472 
1473 !Arguments ------------------------------------
1474   integer,intent(in) :: unt
1475   type(xc_vdw_type), intent(in), optional :: vp
1476 
1477 !Local variables ------------------------------
1478   character(len=1536) :: msg
1479   integer :: napts,ndpts,nqpts,nrpts,id1,id2
1480   real(dp) :: dmax,dtmp,mem_perm,mem_temp
1481   type(xc_vdw_type) :: my_vp
1482 
1483 ! *************************************************************************
1484 
1485   if ( present(vp) ) then
1486     my_vp = vp
1487   else
1488     my_vp = my_vdw_params
1489   end if
1490 
1491   ndpts = my_vp%ndpts
1492   nqpts = my_vp%nqpts
1493   nrpts = my_vp%nrpts
1494 
1495   if ( .not. allocated(dmesh) ) return
1496 
1497   dmax = zero
1498   do id1=1,ndpts
1499     do id2=1,id1
1500       dtmp = sqrt(dmesh(id1)**2+dmesh(id2)**2)
1501       if ( dtmp > dmax ) dmax = dtmp
1502     end do
1503   end do
1504   napts = nint(max(my_vp%acutmin,my_vp%aratio*dmax))
1505 
1506   mem_perm = ( &
1507 &   4 * ndpts * (1 + 5 * ndpts ) + &
1508 &   2 * nqpts + (2 + 4 * nrpts) * nqpts * nqpts &
1509 &   ) * sizeof(one) / 1048576.0_dp
1510 
1511   mem_temp = ( &
1512 &   7 * napts + &
1513 &   4 * nqpts + 3* nqpts * nqpts + &
1514 &   5* nrpts + 2 &
1515 &   ) * sizeof(one) / 1048576.0_dp
1516 
1517   write(msg,'(a,1x,3(a),12(3x,a,1x,i8,1x,"elts",a), &
1518 &   a,10x,a,1x,f10.3,1x,"Mb",a,a,13(3x,a,1x,i8,1x,"elts",a), &
1519 &   a,10x,a,1x,f10.3,1x,"Mb",a,a,14x,a,1x,f10.3,1x,"Mb",2(a),1x,2(a))') &
1520 &   ch10, &
1521 &   '==== vdW-DF memory estimation ====',ch10,ch10, &
1522 &   'd2phidg2 .........',nrpts*nqpts*nqpts,ch10, &
1523 &   'd2phidr2 .........',nrpts*nqpts*nqpts,ch10, &
1524 &   'dmesh ............',ndpts * 4,ch10, &
1525 &   'phi ..............',ndpts*ndpts*4,ch10, &
1526 &   'phi_u ............',ndpts*ndpts*4,ch10, &
1527 &   'phi_bicubic ......',4*4*ndpts*ndpts,ch10, &
1528 &   'phi_u_bicubic ....',4*4*ndpts*ndpts,ch10, &
1529 &   'phig .............',nrpts*nqpts*nqpts,ch10, &
1530 &   'phir .............',nrpts*nqpts*nqpts,ch10, &
1531 &   'phir_u ...........',nrpts*nqpts*nqpts,ch10, &
1532 &   'qmesh ............',nqpts * 4,ch10, &
1533 &   'qpoly_basis.......',nqpts * nqpts * 2,ch10,ch10, &
1534 &   'Permanent =',mem_perm,ch10,ch10, &
1535 &   'amesh ............',napts,ch10, &
1536 &   'amesh_cos ........',napts,ch10, &
1537 &   'amesh_sin ........',napts,ch10, &
1538 &   'btmp .............',nqpts*nqpts,ch10, &
1539 &   'dphida ...........',napts,ch10, &
1540 &   'dphidb ...........',napts,ch10, &
1541 &   'ftmp .............',2*(2*nrpts+1),ch10, &
1542 &   'nu1 ..............',napts,ch10, &
1543 &   'nu2 ..............',napts,ch10, &
1544 &   'qpoly ............',nqpts*3,ch10, &
1545 &   'utmp(q) ..........',nqpts,ch10, &
1546 &   'utmp(r) ..........',nrpts,ch10, &
1547 &   'wtmp .............',nqpts*nqpts*2,ch10,ch10, &
1548 &   'Temporary =',mem_temp,ch10,ch10, &
1549 &   'TOTAL =',mem_perm+mem_temp,ch10, &
1550 &   ch10, &
1551 &   '==================================',ch10
1552 
1553   call wrtout(unt,msg,'COLL')
1554 
1555 end subroutine xc_vdw_memcheck

ABINIT/m_xc_vdw/xc_vdw_read [ Functions ]

[ Top ] [ 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

PARENTS

      driver

CHILDREN

      wrtout

SOURCE

1579 subroutine xc_vdw_read(filename)
1580 
1581 
1582 !This section has been created automatically by the script Abilint (TD).
1583 !Do not modify the following lines by hand.
1584 #undef ABI_FUNC
1585 #define ABI_FUNC 'xc_vdw_read'
1586 !End of the abilint section
1587 
1588   implicit none
1589 
1590 !Arguments ------------------------------------
1591   character(len=*), intent(in)  :: filename
1592 
1593 !Local variables ------------------------------
1594   character(len=512) :: msg
1595   character(len=64) :: format_string
1596   integer :: ncid,dimids(4),varids(13)
1597   integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id
1598   integer :: nderivs,npoly,ndpts,ngpts,nqpts,nrpts
1599 
1600 ! *************************************************************************
1601 
1602   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
1603 
1604 #if defined HAVE_NETCDF
1605   write(msg,'("Reading vdW-DF data from",1x,a)') trim(filename)
1606   call wrtout(std_out,msg,'COLL')
1607 
1608   ! Open file for reading
1609   NETCDF_VDWXC_CHECK(nf90_open(trim(filename),NF90_NOWRITE,ncid=ncid))
1610 
1611   ! Get file format and generator
1612   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'file_format',format_string))
1613   write(msg,'(3x,"File format:",1x,a)') trim(format_string)
1614   call wrtout(std_out,msg,'COLL')
1615   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'generator',msg))
1616   write(msg,'(3x,"Generator:",1x,a)') trim(msg)
1617   call wrtout(std_out,msg,'COLL')
1618 
1619   ! Get attributes
1620   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin))
1621   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio))
1622   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax))
1623   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin))
1624   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut))
1625   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio))
1626   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft))
1627   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional))
1628   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut))
1629   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth))
1630   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft))
1631   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut))
1632   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio))
1633   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut))
1634   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft))
1635   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance))
1636   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks))
1637   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab))
1638 
1639   ! Get dimensions
1640   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nderivs',nderivs_id))
1641   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'npoly',npoly_id))
1642   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ndpts',ndpts_id))
1643   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ngpts',ngpts_id))
1644   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nqpts',nqpts_id))
1645   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nrpts',nrpts_id))
1646 
1647   ! Get varids
1648   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qmesh',varids(2)))
1649   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'dmesh',varids(3)))
1650   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi',varids(4)))
1651   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u',varids(5)))
1652   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_bicubic',varids(6)))
1653   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u_bicubic',varids(7)))
1654   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir',varids(8)))
1655   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir_u',varids(9)))
1656   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidr2',varids(10)))
1657   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phig',varids(11)))
1658   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidg2',varids(12)))
1659   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qpoly_basis',varids(13)))
1660 
1661   ! Get dimensions
1662   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nderivs_id,len=nderivs))
1663   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,npoly_id,len=npoly))
1664   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ndpts_id,len=ndpts))
1665   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ngpts_id,len=ngpts))
1666   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nqpts_id,len=nqpts))
1667   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nrpts_id,len=nrpts))
1668 
1669   ! Get qmesh
1670   ABI_ALLOCATE(qmesh,(nqpts))
1671   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(2),qmesh))
1672 
1673   ! Get dmesh
1674   ABI_ALLOCATE(dmesh,(ndpts))
1675   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(3),dmesh))
1676 
1677   ! Get phi
1678   ABI_ALLOCATE(phi,(ndpts,ndpts,nderivs))
1679   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(4),phi))
1680 
1681   ! Get phi_u
1682   ABI_ALLOCATE(phi_u,(ndpts,ndpts,nderivs))
1683   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(5),phi_u))
1684 
1685   ! Get phi_bicubic
1686   ABI_ALLOCATE(phi_bicubic,(nderivs,nderivs,ndpts,ndpts))
1687   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(6),phi_bicubic))
1688 
1689   ! Get phi_u_bicubic
1690   ABI_ALLOCATE(phi_u_bicubic,(nderivs,nderivs,ndpts,ndpts))
1691   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(7),phi_u_bicubic))
1692 
1693   ! Get phir
1694   ABI_ALLOCATE(phir,(nrpts,nqpts,nqpts))
1695   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(8),phir))
1696 
1697   ! Get phir_u
1698   ABI_ALLOCATE(phir_u,(nrpts,nqpts,nqpts))
1699   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(9),phir_u))
1700 
1701   ! Get d2phidr2
1702   ABI_ALLOCATE(d2phidr2,(nrpts,nqpts,nqpts))
1703   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(10),d2phidr2))
1704 
1705   ! Get phig
1706   ABI_ALLOCATE(phig,(nrpts,nqpts,nqpts))
1707   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(11),phig))
1708 
1709   ! Get d2phidg2
1710   ABI_ALLOCATE(d2phidg2,(nrpts,nqpts,nqpts))
1711   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(12),d2phidg2))
1712 
1713   ! Get qpoly_basis
1714   ABI_ALLOCATE(qpoly_basis,(nqpts,nqpts,npoly))
1715   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(13),qpoly_basis))
1716 
1717   ! Close file
1718   NETCDF_VDWXC_CHECK(nf90_close(ncid))
1719 
1720   ! Update my_vdw_params
1721   my_vdw_params%ndpts = ndpts
1722   my_vdw_params%ngpts = ngpts
1723   my_vdw_params%nqpts = nqpts
1724   my_vdw_params%nrpts = nrpts
1725 #else
1726   MSG_ERROR('reading vdW-DF variables requires NetCDF')
1727 #endif
1728 
1729   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
1730 
1731 end subroutine xc_vdw_read

ABINIT/m_xc_vdw/xc_vdw_set_functional [ Functions ]

[ Top ] [ 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

PARENTS

CHILDREN

      wrtout

SOURCE

1752 subroutine xc_vdw_set_functional(vdw_func,vdw_zab)
1753 
1754 
1755 !This section has been created automatically by the script Abilint (TD).
1756 !Do not modify the following lines by hand.
1757 #undef ABI_FUNC
1758 #define ABI_FUNC 'xc_vdw_set_functional'
1759 !End of the abilint section
1760 
1761   implicit none
1762 
1763 !Arguments ------------------------------------
1764   integer, intent(in)  :: vdw_func
1765   real(dp), optional, intent(in) :: vdw_zab
1766 
1767 !Local variables ------------------------------
1768   character(len=512) :: msg
1769 
1770 ! *************************************************************************
1771 
1772   my_vdw_params%functional = vdw_func
1773   if ( present(vdw_zab) ) then
1774     my_vdw_params%zab = vdw_zab
1775   end if
1776 
1777 end subroutine xc_vdw_set_functional

ABINIT/m_xc_vdw/xc_vdw_set_params [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_set_params

FUNCTION

  Imports external vdW-DF parameters.

 INPUT
  vdw_params= van der Waals parameters

PARENTS

      m_xc_vdw

CHILDREN

      wrtout

SOURCE

1798 subroutine xc_vdw_set_params(vdw_params)
1799 
1800 
1801 !This section has been created automatically by the script Abilint (TD).
1802 !Do not modify the following lines by hand.
1803 #undef ABI_FUNC
1804 #define ABI_FUNC 'xc_vdw_set_params'
1805 !End of the abilint section
1806 
1807   implicit none
1808 
1809 !Arguments ------------------------------------
1810   type(xc_vdw_type), intent(in)  :: vdw_params
1811 
1812 !Local variables ------------------------------
1813   character(len=512) :: msg
1814 
1815 ! *************************************************************************
1816 
1817   my_vdw_params%functional = vdw_params%functional
1818   my_vdw_params%zab        = vdw_params%zab
1819   my_vdw_params%ndpts      = vdw_params%ndpts
1820   my_vdw_params%dcut       = vdw_params%dcut
1821   my_vdw_params%dratio     = vdw_params%dratio
1822   my_vdw_params%dsoft      = vdw_params%dsoft
1823   my_vdw_params%phisoft    = vdw_params%phisoft
1824   my_vdw_params%nqpts      = vdw_params%nqpts
1825   my_vdw_params%qcut       = vdw_params%qcut
1826   my_vdw_params%qratio     = vdw_params%qratio
1827   my_vdw_params%nrpts      = vdw_params%nrpts
1828   my_vdw_params%rcut       = vdw_params%rcut
1829   my_vdw_params%rsoft      = vdw_params%rsoft
1830   my_vdw_params%ngpts      = vdw_params%ngpts
1831   my_vdw_params%gcut       = vdw_params%gcut
1832   my_vdw_params%acutmin    = vdw_params%acutmin
1833   my_vdw_params%aratio     = vdw_params%aratio
1834   my_vdw_params%damax      = vdw_params%damax
1835   my_vdw_params%damin      = vdw_params%damin
1836   my_vdw_params%nsmooth    = vdw_params%nsmooth
1837   my_vdw_params%tolerance  = vdw_params%tolerance
1838   my_vdw_params%tweaks     = vdw_params%tweaks
1839 
1840 end subroutine xc_vdw_set_params

ABINIT/m_xc_vdw/xc_vdw_show [ Functions ]

[ Top ] [ 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

PARENTS

      driver,vdw_kernelgen

CHILDREN

      wrtout

SOURCE

1862 subroutine xc_vdw_show(unt,vp)
1863 
1864 
1865 !This section has been created automatically by the script Abilint (TD).
1866 !Do not modify the following lines by hand.
1867 #undef ABI_FUNC
1868 #define ABI_FUNC 'xc_vdw_show'
1869 !End of the abilint section
1870 
1871   implicit none
1872 
1873 !Arguments ------------------------------------
1874   integer,intent(in) :: unt
1875   type(xc_vdw_type), intent(in), optional :: vp
1876 
1877 !Local variables ------------------------------
1878   character(len=1536) :: msg
1879   type(xc_vdw_type) :: my_vp
1880 
1881 ! *************************************************************************
1882 
1883   if ( present(vp) ) then
1884     my_vp = vp
1885   else
1886     my_vp = my_vdw_params
1887   end if
1888 
1889   write(msg,'(a,1x,3(a),3x,a,1x,i9,a,3x,a,1x,f9.3,a, &
1890 &   5(3x,a,1x,i9,a),14(3x,a,1x,f9.3,a),1(3x,a,1x,i9,a),a,1x,a,a)') &
1891 &   ch10, &
1892 &   '==== vdW-DF internal parameters ====',ch10,ch10, &
1893 &   'functional .............',my_vp%functional,ch10, &
1894 &   'zab ....................',my_vp%zab,ch10, &
1895 &   'd-mesh points ..........',my_vp%ndpts,ch10, &
1896 &   'q-mesh points ..........',my_vp%nqpts,ch10, &
1897 &   'r-mesh points ..........',my_vp%nrpts,ch10, &
1898 &   'g-mesh points ..........',my_vp%ngpts,ch10, &
1899 &   'smoothing iterations ...',my_vp%nsmooth,ch10, &
1900 &   'd-mesh cut-off .........',my_vp%dcut,ch10, &
1901 &   'd-mesh ratio ...........',my_vp%dratio,ch10, &
1902 &   'd-mesh softening .......',my_vp%dsoft,ch10, &
1903 &   'phi softening ..........',my_vp%phisoft,ch10, &
1904 &   'q-mesh cut-off .........',my_vp%qcut,ch10, &
1905 &   'q-mesh ratio ...........',my_vp%qratio,ch10, &
1906 &   'r-mesh cut-off .........',my_vp%rcut,ch10, &
1907 &   'r-mesh softening .......',my_vp%rsoft,ch10, &
1908 &   'g-mesh cut-off .........',my_vp%gcut,ch10, &
1909 &   'a-mesh min cut-off .....',my_vp%acutmin,ch10, &
1910 &   'a-mesh ratio ...........',my_vp%aratio,ch10, &
1911 &   'a-mesh delta max .......',my_vp%damax,ch10, &
1912 &   'a-mesh delta min .......',my_vp%damin,ch10, &
1913 &   'global tolerance .......',my_vp%tolerance,ch10, &
1914 &   'tweaks .................',my_vp%tweaks,ch10, &
1915 &   ch10, &
1916 &   '====================================',ch10
1917 
1918   call wrtout(unt,msg,'COLL')
1919 
1920 end subroutine xc_vdw_show

ABINIT/m_xc_vdw/xc_vdw_status [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_status

FUNCTION

  Returns the status of the main vdW-DF switch.

INPUTS

  None

PARENTS

  Will be filled automatically by the parent script

CHILDREN

  Will be filled automatically by the parent script

SOURCE

1941 function xc_vdw_status()
1942 
1943 
1944 !This section has been created automatically by the script Abilint (TD).
1945 !Do not modify the following lines by hand.
1946 #undef ABI_FUNC
1947 #define ABI_FUNC 'xc_vdw_status'
1948 !End of the abilint section
1949 
1950   implicit none
1951 
1952 !Arguments ------------------------------------
1953 
1954 !Local variables ------------------------------
1955   logical :: xc_vdw_status
1956 
1957 ! *************************************************************************
1958 
1959   xc_vdw_status = vdw_switch
1960 
1961 end function xc_vdw_status

ABINIT/m_xc_vdw/xc_vdw_trigger [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_trigger

FUNCTION

  Switches on and off the calculation of vdW interactions.

INPUTS

  condition= boolean condition to trigger the calculations

PARENTS

      driver,scprqt

CHILDREN

      wrtout

SOURCE

1982 subroutine xc_vdw_trigger(condition)
1983 
1984 
1985 !This section has been created automatically by the script Abilint (TD).
1986 !Do not modify the following lines by hand.
1987 #undef ABI_FUNC
1988 #define ABI_FUNC 'xc_vdw_trigger'
1989 !End of the abilint section
1990 
1991   implicit none
1992 
1993 !Arguments ------------------------------------
1994   logical, intent(in) :: condition
1995 
1996 !Local variables ------------------------------
1997   character(len=512) :: msg
1998 
1999 ! *************************************************************************
2000 
2001   if ( my_vdw_params%functional /= 0 ) then
2002     if ( vdw_switch ) then
2003       write (msg,'(1x,a,a)') &
2004 &       "[vdW-DF] xc_vdw_trigger: keeping xc_vdw_aggregate enabled",ch10
2005     else
2006       if ( condition ) then
2007         write (msg,'(1x,a,a)') &
2008 &         "[vdW-DF] xc_vdw_trigger: enabling xc_vdw_aggregate",ch10
2009       else
2010         write (msg,'(1x,a,a)') &
2011 &         "[vdW-DF] xc_vdw_trigger: disabling xc_vdw_aggregate",ch10
2012       end if
2013 
2014       vdw_switch = condition
2015     end if
2016 
2017     call wrtout(std_out,msg,'COLL')
2018   end if
2019 
2020 end subroutine xc_vdw_trigger

ABINIT/m_xc_vdw/xc_vdw_type [ Types ]

[ Top ] [ Types ]

NAME

  xc_vdw_type

FUNCTION

  Tunable parameters for calculations relying upon van der Waals XC.

SOURCE

 76   type xc_vdw_type
 77 
 78     integer :: functional = 0
 79     ! Functional to use for the van der Waals correction
 80     ! (see the description of the vdw_xc input variable for possible values)
 81 
 82     real(dp) :: zab = -0.8491_dp
 83     ! Parameter of the vdW functional, introduced by DRSLL04
 84 
 85     integer :: ndpts = 20
 86     ! Number of points for the d-mesh
 87 
 88     real(dp) :: dcut = 30.0_dp
 89     ! Cut-off for the d-mesh
 90 
 91     real(dp) :: dratio = 20.0_dp
 92     ! Ratio between highest and lowest d
 93 
 94     real(dp) :: dsoft = 1.0_dp
 95     ! Distance within the d-mesh below which the kernel will be softened
 96 
 97     real(dp) :: phisoft = -1.0_dp
 98     ! Value of the softened kernel for d=0
 99     ! Will be set automatically if negative (default behaviour)
100 
101     integer :: nqpts = 30
102     ! Number of points of the q-mesh
103 
104     real(dp) :: qcut = 5.0_dp
105     ! Cut-off for the q-mesh
106 
107     real(dp) :: qratio = 20.0_dp
108     ! Ratio between highest and lowest q
109 
110     integer :: nrpts = 2048
111     ! Number of points in real space for the kernel
112 
113     real(dp) :: rcut = 100.0_dp
114     ! Real-space cut-off for the kernel
115 
116     real(dp) :: rsoft = 0.0_dp
117     ! Radius below which the kernel will be softened
118 
119     integer :: ngpts = -1
120     ! Number of wavevectors for the kernel
121 
122     real(dp) :: gcut = 5.0_dp
123     ! Wavevector cut-off for the kernel (=sqrt(ecut))
124 
125     real(dp) :: acutmin = 10.0_dp
126     ! Minimum cut-off for the angular mesh
127 
128     real(dp) :: aratio = 30.0_dp
129     ! Ratio between highest and lowest a !DEBUG this definition has to
130     ! be corrected since its use in the code is different.
131 
132     real(dp) :: damax = 0.5_dp
133     ! Maximum variation in the angular mesh
134 
135     real(dp) :: damin = 1.0e-2_dp
136     ! Minimum variation in the angular mesh
137 
138     integer :: nsmooth = 12
139     ! Saturation level to smoothen q0 near qcut
140 
141     real(dp) :: tolerance = 1.0e-13_dp
142     ! Global numerical tolerance for the boundary values of the kernel
143 
144     integer :: tweaks = 0
145     ! Tweaks of the implementation (modify with extreme care)
146 
147   end type xc_vdw_type

ABINIT/m_xc_vdw/xc_vdw_write [ Functions ]

[ Top ] [ 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

PARENTS

      driver,vdw_kernelgen

CHILDREN

      wrtout

SOURCE

2044 subroutine xc_vdw_write(filename)
2045 
2046 
2047 !This section has been created automatically by the script Abilint (TD).
2048 !Do not modify the following lines by hand.
2049 #undef ABI_FUNC
2050 #define ABI_FUNC 'xc_vdw_write'
2051 !End of the abilint section
2052 
2053   implicit none
2054 
2055 !Arguments ------------------------------------
2056   character(len=*), intent(in)  :: filename
2057 
2058 !Local variables ------------------------------
2059   character(len=512) :: msg
2060   integer :: ncid,dimids(4),varids(13)
2061   integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id
2062 
2063 ! *************************************************************************
2064 
2065   VDWXC_DBG_ENTER("COLL",ABI_FUNC)
2066 
2067 #if defined HAVE_NETCDF
2068   write(msg,'(a,1x,"Writing vdW-DF data to",1x,a)') ch10,trim(filename)
2069   call wrtout(std_out,msg,'COLL')
2070 
2071   ! Create file (overwriting any existing data)
2072   NETCDF_VDWXC_CHECK(nf90_create(trim(filename),NF90_CLOBBER,ncid))
2073 
2074   ! Write attributes
2075   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'file_format',vdw_format_string))
2076   write(msg,'(a,1x,a,1x,"for",1x,a,1x,"(build",1x,a,")")') &
2077 &   PACKAGE_NAME,ABINIT_VERSION,ABINIT_TARGET,ABINIT_VERSION_BUILD
2078   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'generator',msg))
2079   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin))
2080   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio))
2081   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax))
2082   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin))
2083   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut))
2084   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio))
2085   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft))
2086   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional))
2087   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut))
2088   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth))
2089   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft))
2090   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut))
2091   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio))
2092   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut))
2093   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft))
2094   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance))
2095   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks))
2096   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab))
2097 
2098   ! Define dimensions
2099   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nderivs',4,nderivs_id))
2100   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'npoly',2,npoly_id))
2101   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ndpts',my_vdw_params%ndpts,ndpts_id))
2102   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ngpts',my_vdw_params%ngpts,ngpts_id))
2103   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nqpts',my_vdw_params%nqpts,nqpts_id))
2104   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nrpts',my_vdw_params%nrpts,nrpts_id))
2105 
2106   ! Define qmesh
2107   dimids(1) = nqpts_id
2108   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qmesh',NF90_DOUBLE,dimids(1),varids(2)))
2109 
2110   ! Define dmesh
2111   dimids(1) = ndpts_id
2112   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'dmesh',NF90_DOUBLE,dimids(1),varids(3)))
2113 
2114   ! Define phi
2115   dimids(1) = ndpts_id
2116   dimids(2) = ndpts_id
2117   dimids(3) = nderivs_id
2118   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi',NF90_DOUBLE,dimids(1:3),varids(4)))
2119 
2120   ! Define phi_u
2121   dimids(1) = ndpts_id
2122   dimids(2) = ndpts_id
2123   dimids(3) = nderivs_id
2124   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u',NF90_DOUBLE,dimids(1:3),varids(5)))
2125 
2126   ! Define phi_bicubic
2127   dimids(1) = nderivs_id
2128   dimids(2) = nderivs_id
2129   dimids(3) = ndpts_id
2130   dimids(4) = ndpts_id
2131   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_bicubic',NF90_DOUBLE,dimids(1:4),varids(6)))
2132 
2133   ! Define phi_u_bicubic
2134   dimids(1) = nderivs_id
2135   dimids(2) = nderivs_id
2136   dimids(3) = ndpts_id
2137   dimids(4) = ndpts_id
2138   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u_bicubic',NF90_DOUBLE,dimids(1:4),varids(7)))
2139 
2140   ! Define phir
2141   dimids(1) = nrpts_id
2142   dimids(2) = nqpts_id
2143   dimids(3) = nqpts_id
2144   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir',NF90_DOUBLE,dimids(1:3),varids(8)))
2145 
2146   ! Define phir_u
2147   dimids(1) = nrpts_id
2148   dimids(2) = nqpts_id
2149   dimids(3) = nqpts_id
2150   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir_u',NF90_DOUBLE,dimids(1:3),varids(9)))
2151 
2152   ! Define d2phidr2
2153   dimids(1) = nrpts_id
2154   dimids(2) = nqpts_id
2155   dimids(3) = nqpts_id
2156   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidr2',NF90_DOUBLE,dimids(1:3),varids(10)))
2157 
2158   ! Define phig
2159   dimids(1) = nrpts_id
2160   dimids(2) = nqpts_id
2161   dimids(3) = nqpts_id
2162   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phig',NF90_DOUBLE,dimids(1:3),varids(11)))
2163 
2164   ! Define d2phidg2
2165   dimids(1) = nrpts_id
2166   dimids(2) = nqpts_id
2167   dimids(3) = nqpts_id
2168   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidg2',NF90_DOUBLE,dimids(1:3),varids(12)))
2169 
2170   ! Define qpoly_basis
2171   dimids(1) = nqpts_id
2172   dimids(2) = nqpts_id
2173   dimids(3) = npoly_id
2174   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qpoly_basis',NF90_DOUBLE,dimids(1:3),varids(13)))
2175 
2176   ! Stop definitions
2177   NETCDF_VDWXC_CHECK(nf90_enddef(ncid))
2178 
2179   ! Write variables
2180   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(2),qmesh))
2181   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(3),dmesh))
2182   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(4),phi))
2183   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(5),phi_u))
2184   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(6),phi_bicubic))
2185   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(7),phi_u_bicubic))
2186   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(8),phir))
2187   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(9),phir_u))
2188   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(10),d2phidr2))
2189   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(11),phig))
2190   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(12),d2phidg2))
2191   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(13),qpoly_basis))
2192 
2193   ! Close file
2194   NETCDF_VDWXC_CHECK(nf90_close(ncid))
2195 #else
2196   MSG_WARNING('writing vdW-DF variables requires NetCDF - skipping')
2197 #endif
2198 
2199   VDWXC_DBG_EXIT("COLL",ABI_FUNC)
2200 
2201 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.

PARENTS

CHILDREN

      wrtout

SOURCE

3311 subroutine vdw_df_netcdf_ioerr(ncerr,file_name,file_line)
3312 
3313 !Arguments ------------------------------------
3314 
3315 !This section has been created automatically by the script Abilint (TD).
3316 !Do not modify the following lines by hand.
3317 #undef ABI_FUNC
3318 #define ABI_FUNC 'vdw_df_netcdf_ioerr'
3319 !End of the abilint section
3320 
3321  integer,intent(in) :: ncerr
3322  character(len=*),optional,intent(in) :: file_name
3323  integer,optional,intent(in) :: file_line
3324 
3325 !Local variables-------------------------------
3326  character(len=1024) :: msg
3327  character(len=fnlen) :: my_file_name = "N/D"
3328  integer :: my_file_line = -1
3329 
3330 ! *************************************************************************
3331 
3332   if ( present(file_name) ) then
3333     my_file_name = trim(file_name)
3334   end if
3335   if ( present(file_line) ) then
3336     my_file_line = file_line
3337   end if
3338 
3339 #if defined HAVE_NETCDF
3340   write(msg,'(a,a,3x,a)') &
3341 &  'NetCDF returned the error:',ch10,trim(nf90_strerror(ncerr))
3342   if ( ncerr /= NF90_NOERR ) then
3343     call die(msg,trim(my_file_name),my_file_line)
3344   end if
3345 #endif
3346 
3347 end subroutine vdw_df_netcdf_ioerr