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-2022 ABINIT group (Yann Pouillon, Camilo Espejo)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  DRSLL04 = doi:10.1103/PhysRevLett.92.246401 [[cite:Dion2004]]
  RS09 = doi:10.1103/PhysRevLett.103.096102 [[cite:Romanperez2009]]

SOURCE

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

ABINIT/m_xc_vdw/vdw_df_create_mesh [ Functions ]

[ 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

SOURCE

2860 subroutine vdw_df_create_mesh(mesh,npts,mesh_type,mesh_cutoff, &
2861 & mesh_inc,mesh_ratio,mesh_tolerance,mesh_file,avoid_zero,exact_max)
2862 
2863 !Arguments ------------------------------------
2864   integer,intent(in) :: npts,mesh_type
2865   real(dp),intent(in) :: mesh_cutoff
2866   real(dp),intent(out) :: mesh(npts)
2867   logical,optional,intent(in) :: avoid_zero,exact_max
2868   real(dp),optional,intent(in) :: mesh_inc,mesh_ratio,mesh_tolerance
2869   character(len=*),optional,intent(in) :: mesh_file
2870 
2871 !Local variables ------------------------------
2872   logical :: my_avoid_zero,my_exact_max
2873   integer :: im,unt
2874   real(dp) :: exp_inc,exp_ofs,exp_tol,lin_inc,log_inc,log_ofs,log_tol
2875   character(len=500) :: msg
2876 
2877 ! *************************************************************************
2878 
2879 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2880   DBG_ENTER("COLL")
2881 #endif
2882 
2883   if ( present(avoid_zero) ) then
2884     my_avoid_zero = avoid_zero
2885   else
2886     my_avoid_zero = .false.
2887   end if
2888 
2889   if ( present(exact_max) ) then
2890     my_exact_max = exact_max
2891   else
2892     my_exact_max = .false.
2893   end if
2894 
2895 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2896       write(msg,'(1x,a,a,1x,a,1x,i6,a,1x,a,i6,a,1x,a,1x,e10.3)') &
2897 &       "> Mesh parameters",ch10, &
2898 &       ">   npts    :",npts,ch10, &
2899 &       ">   type    :",mesh_type,ch10, &
2900 &       ">   cut-off :",mesh_cutoff
2901       call wrtout(std_out,msg,'COLL')
2902 #endif
2903 
2904   select case(mesh_type)
2905 
2906     case(0)
2907 
2908       if ( present(mesh_file) ) then
2909 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2910         write(msg,'(1x,a,1x,a)') ">   file    :",mesh_file
2911         call wrtout(std_out,msg,'COLL')
2912 #endif
2913         if (open_file(mesh_file,msg, newunit=unt,status="old",action="read") /= 0) then
2914           ABI_ERROR(msg)
2915         end if
2916         read(unt,'(e20.8)') mesh
2917         close(unit=unt)
2918       else
2919         ABI_ERROR("  You must specify a mesh file for this mesh type!")
2920       end if
2921 
2922     case(1)
2923 
2924       if ( present(mesh_ratio) ) then
2925         if ( present(mesh_tolerance) ) then
2926           log_tol = mesh_tolerance
2927         else
2928           log_tol = sqrt(my_vdw_params%tolerance)
2929         end if
2930 !DEBUG
2931 !        log_inc = mesh_ratio / (npts - 2)
2932         log_inc = log(mesh_ratio) / (npts - 1)
2933 !ENDDEBUG
2934         log_ofs = mesh_cutoff / (exp(log_inc*(npts - 1)) - 1)
2935 
2936 !!DEBUG
2937 !! Previous definition of log_ofs implies that the mesh
2938 !! points are to be calculated as:
2939 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 )
2940 !! which in turn means that x(0)=0 and is not consistent with the
2941 !! mesh(im) points given at line 3056.
2942 !! If the latter is adopted mesh_ratio above should be computed as:
2943 !! mesh_ratio=log ( (x(n)-x(n-1))/(x(2)-x(1)) ) and not merely as
2944 !! mesh_ratio= (x(n)-x(n-1))/(x(2)-x(1))
2945 !! in order to log_inc definition makes sense.
2946 !! The latter relation between log_inc and mesh_ratio is valid even if
2947 !! mesh points do not start at 0:
2948 !! x(i)=log_ofs( exp(log_inc*(npts - 1)) - 1 ) + x_0
2949 !!ENDDEBUG
2950 
2951 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2952         write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') &
2953 &         ">   ratio   :",mesh_ratio,ch10, &
2954 &         ">   log_inc :",log_inc,ch10, &
2955 &         ">   log_ofs :",log_ofs,ch10, &
2956 &         ">   log_tol :",log_tol
2957         call wrtout(std_out,msg,'COLL')
2958 #endif
2959         if ( log_inc < log_tol ) then
2960           ABI_WARNING("  The logarithmic increment is very low!")
2961         end if
2962 
2963         do im = 1,npts
2964 !DEBUG
2965 !          mesh(im) = log_ofs * exp(log_inc * (im-1)) this is not consistent
2966 !       with definition of log_ofs at line 3025.
2967           mesh(im) = log_ofs * (exp(log_inc * (im-1)) - 1)
2968 !ENDDEBUG
2969         end do
2970 
2971         !FIXME: check the correctness of 0.2
2972         if ( my_avoid_zero ) then
2973           mesh(1) = exp(log_inc * 0.2) * mesh_cutoff / &
2974 &           (exp(log_inc*(npts - 1)) - one)
2975         end if
2976 
2977         if ( my_exact_max ) then
2978           mesh(npts) = mesh_cutoff
2979         end if
2980       else
2981         ABI_ERROR("  You must specify a mesh ratio for this mesh type!")
2982       end if
2983 
2984       ! Make sure the last point is qcut
2985       mesh(npts) = mesh_cutoff
2986 
2987     case(2)
2988 
2989       if ( present(mesh_inc) ) then
2990         lin_inc = mesh_inc
2991       else
2992         lin_inc = mesh_cutoff / npts
2993       end if
2994 
2995 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2996       write(msg,'(1x,a,1x,e10.3)') ">   lin_inc :",lin_inc
2997       call wrtout(std_out,msg,'COLL')
2998 #endif
2999       do im = 1,npts
3000         mesh(im) = lin_inc * (im - 1)
3001       end do
3002 
3003     case(3)
3004 
3005       if ( present(mesh_ratio) ) then
3006         if ( present(mesh_tolerance) ) then
3007           exp_tol = mesh_tolerance
3008         else
3009           exp_tol = sqrt(my_vdw_params%tolerance)
3010         end if
3011         exp_inc = mesh_ratio / (npts - 2)
3012         exp_ofs = mesh_cutoff * (log(exp_inc*(npts - 1)) + 1)
3013 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
3014         write(msg,'(1x,a,1x,e10.3,3(a,1x,a,1x,e10.3))') &
3015 &         ">   ratio   :",mesh_ratio,ch10, &
3016 &         ">   exp_inc :",exp_inc,ch10, &
3017 &         ">   exp_ofs :",exp_ofs,ch10, &
3018 &         ">   exp_tol :",exp_tol
3019         call wrtout(std_out,msg,'COLL')
3020 #endif
3021         if ( log_inc < log_tol ) then
3022           ABI_WARNING("  The logarithmic increment is very low!")
3023         end if
3024 
3025         do im = 1,npts
3026           mesh(im) = exp_ofs / log(exp_inc * (im-1))
3027         end do
3028       else
3029         ABI_ERROR("  You must specify a mesh ratio for this mesh type!")
3030       end if
3031 
3032     case default
3033 
3034       write(msg,'(2x,a,1x,i2)') "Unknown mesh type",mesh_type
3035       ABI_ERROR(msg)
3036 
3037   end select
3038 
3039 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
3040   DBG_EXIT("COLL")
3041 #endif
3042 
3043 end subroutine vdw_df_create_mesh

ABINIT/m_xc_vdw/vdw_df_filter [ Functions ]

[ 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

SOURCE

2200 subroutine vdw_df_filter(nqpts,nrpts,rcut,gcut,ngpts,sofswt)
2201 
2202 !Arguments ------------------------------------
2203   integer,intent(in) :: nqpts,nrpts,sofswt
2204   integer,intent(out) :: ngpts
2205   real(dp),intent(in) :: rcut
2206   real(dp),intent(inout) :: gcut
2207 !Local variables ------------------------------
2208   integer :: ig,iq1,iq2,ir,rfftt,id1,ndpts,ng
2209   real(dp) :: dg,dr,lstep,ptmp,q1,q2,qcut,qtol,un,xr
2210   real(dp) :: x1,x2,yq1,yqn,yr1,yrn,gmax
2211   real(dp),allocatable :: gmesh(:),rmesh(:),utmp(:),phiraux(:,:,:)
2212   real(dp),allocatable :: q1sat(:),q2sat(:)
2213 ! *************************************************************************
2214 
2215   DBG_ENTER("COLL")
2216 
2217   ! Init
2218   qtol = my_vdw_params%tolerance
2219   qcut = my_vdw_params%qcut
2220   write(std_out,*) 'gcut before =', gcut  
2221   dr = rcut / nrpts
2222   ngpts = nrpts
2223   dg = pi / rcut
2224   gcut =  pi / dr  !kmax in siesta
2225   gmax = 10.639734883681651 !kcut in siesta for a particular 
2226 ! system (Ar). This should be computed each time in terms of 
2227 ! the FFT and BZ parameters, not hard wired like here. 
2228   ng = 1 + int(gmax/dg) ! nk in siesta
2229   rfftt = 2 !Type of radial sin transform
2230   write(std_out,*) 'gmax from siesta value for kmax=', gmax
2231   write(std_out,*) 'ng from siesta gmax=', ng
2232   write(std_out,*) 'gcut after =', gcut
2233   write(std_out,*) 'ngpts after =', ngpts
2234 
2235   ABI_MALLOC(utmp,(0:nrpts))
2236 
2237   ! Create radial mesh for radial FT: 
2238   ! shared/common/src/32_util/radsintr.F90
2239   ABI_MALLOC(rmesh,(0:nrpts))
2240   forall(ir=1:nrpts) rmesh(ir) = dr * dble(ir)
2241 
2242   if ( sofswt == 1 ) then
2243   ! Create reciprocal radial mesh
2244    ABI_MALLOC(gmesh,(0:ngpts))
2245 !----my debug-------------------------------   
2246    ABI_MALLOC(phiraux,(0:nrpts,nqpts,nqpts))
2247 !----end my debug---------------------------
2248    forall(ig=0:ngpts) gmesh(ig) = dg * dble(ig)
2249   end if
2250 
2251   ABI_MALLOC(q1sat,(2))
2252   ABI_MALLOC(q2sat,(2))
2253 
2254   ! Build filtered kernel for each (q1,q2) pair
2255   do iq1=1,nqpts
2256     !Inverse saturation of q1-mesh value 
2257     x1 = 0
2258     x2 = 1.5_dp * qcut
2259     do
2260       q1 = (x1+x2)/2
2261       call vdw_df_saturation( q1, qcut, q1sat )
2262       if (abs(qmesh(iq1)-q1sat(1))<qtol) then
2263         exit
2264       else if (q1sat(1) < qmesh(iq1)) then
2265         x1 = q1
2266       else
2267         x2 = q1
2268       end if
2269     end do ! end calculation of unsaturated qmesh(iq1)
2270 
2271     do iq2=1,iq1      
2272     !Rather than this q2 = qmesh(iq2), the values of q-mesh
2273     !are unsaturated. 
2274       x1 = 0
2275       x2 = 1.5_dp * qcut
2276       do
2277         q2 = (x1+x2)/2
2278         call vdw_df_saturation( q2, qcut, q2sat )
2279         if (abs(qmesh(iq2)-q2sat(1))<qtol) then
2280           exit
2281         else if (q2sat(1) < qmesh(iq2)) then
2282           x1 = q2
2283         else
2284           x2 = q2
2285         end if                           
2286       end do ! end calculation of unsaturated qmesh(iq2)
2287 
2288       if ( sofswt == 1 ) then
2289       ! Build kernel in real space
2290       ! Note: smoothly going to zero when approaching rcut
2291       ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values.
2292        do ir=1,nrpts
2293          xr=rmesh(ir)
2294          phir(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * &
2295 &          (one - ( xr / rmesh(nrpts))**8)**4
2296        end do
2297 !-----debug-------------------------------------     
2298        phiraux(:,iq1,iq2) = phir(:,iq1,iq2)
2299 !-----end debug---------------------------------
2300 
2301       ! Obtain kernel in reciprocal space
2302        call radsintr(phir(:,iq1,iq2),phig(:,iq1,iq2),ngpts,nrpts,gmesh,rmesh,yq1,yqn,rfftt)
2303        if (rfftt == 2) then
2304          phig(:,iq1,iq2) = phig(:,iq1,iq2) * (2*pi)**1.5_dp
2305        end if
2306 
2307       ! Filter in reciprocal space
2308        phig(ng:ngpts,iq1,iq2) = zero 
2309        do ig=1,ng
2310          phig(ig,iq1,iq2) = phig(ig,iq1,iq2) * &
2311 &          (one - ( gmesh(ig) / gmesh(ng))**8)**4
2312        end do
2313 
2314       ! Go back to real space
2315        call radsintr(phig(:,iq1,iq2),phir(:,iq1,iq2),nrpts,ngpts,rmesh,gmesh,yr1,yrn,rfftt)
2316        phir(:,iq1,iq2) = phir(:,iq1,iq2) / (2*pi)**1.5_dp
2317        
2318       ! Calculate second derivative in real space
2319        d2phidr2(0,iq1,iq2) = -half
2320        utmp(0) = (three / dr) * (phir(1,iq1,iq2)-phir(0,iq1,iq2)) / dr
2321        do ir=1,nrpts-1
2322          ptmp = half * d2phidr2(ir-1,iq1,iq2) + two
2323          d2phidr2(ir,iq1,iq2) = (half - one) / ptmp
2324          utmp(ir) = (three * (phir(ir+1,iq1,iq2) + phir(ir-1,iq1,iq2) - &
2325 &          two*phir(ir,iq1,iq2)) / (dr**2) - half * utmp(ir-1)) / ptmp
2326        end do
2327        un = (three / dr) * (phir(nrpts-1,iq1,iq2)-phir(nrpts,iq1,iq2)) / dr
2328        d2phidr2(nrpts,iq1,iq2) = (un - half * utmp(nrpts-1)) /  &
2329          (half * d2phidr2(nrpts-1,iq1,iq2) + one) 
2330        do ir=nrpts-1,0,-1
2331          d2phidr2(ir,iq1,iq2) = d2phidr2(ir,iq1,iq2) * &
2332 &          d2phidr2(ir+1,iq1,iq2) + utmp(ir)
2333        end do
2334 
2335       ! Calculate second derivative in reciprocal space
2336        d2phidg2(0,iq1,iq2) = -half 
2337        utmp(0) = (three / dg) * (phig(1,iq1,iq2)-phig(0,iq1,iq2)) / dg
2338        do ig=1,ngpts-1
2339          ptmp = half * d2phidg2(ig-1,iq1,iq2) + two
2340          d2phidg2(ig,iq1,iq2) = (half - one) / ptmp
2341          utmp(ig) = (three * (phig(ig+1,iq1,iq2) + phig(ig-1,iq1,iq2) - &
2342 &          two*phig(ig,iq1,iq2)) / (dg**2) - half * utmp(ig-1)) / ptmp
2343        end do
2344        un = (three / dg) * (phig(ngpts-1,iq1,iq2)-phig(ngpts,iq1,iq2)) / dg
2345        d2phidg2(ngpts,iq1,iq2) = (un - half * utmp(ngpts-1)) /  &
2346          (half * d2phidg2(ngpts-1,iq1,iq2) + one)
2347        do ig=ngpts-1,0,-1
2348          d2phidg2(ig,iq1,iq2) = d2phidg2(ig,iq1,iq2) * &                          
2349 &         d2phidg2(ig+1,iq1,iq2) + utmp(ig)
2350        end do
2351 
2352       ! Symmetrize kernels & derivatives
2353        phir(:,iq2,iq1) = phir(:,iq1,iq2)
2354        phig(:,iq2,iq1) = phig(:,iq1,iq2)
2355        d2phidr2(:,iq2,iq1) = d2phidr2(:,iq1,iq2)
2356        d2phidg2(:,iq2,iq1) = d2phidg2(:,iq1,iq2)
2357 !---------my debug----------------------
2358        phiraux(:,iq2,iq1) = phiraux(:,iq1,iq2)
2359 !----end my debug----------------
2360       else! sofswt == 1
2361 
2362       ! Build unsoftened kernel in real space
2363       ! Note: smoothly going to zero when approaching rcut
2364       ! radial mesh is indeed a mesh of |\vec{r1}-\vec{r2}| values.
2365        do ir=0,nrpts
2366          xr=rmesh(ir)
2367          phir_u(ir,iq1,iq2) = vdw_df_interpolate(q1*xr,q2*xr,sofswt) * &
2368 &          (one - ( xr / rmesh(nrpts))**8)**4
2369        end do
2370 
2371       ! Symmetrize unsoftened kernel
2372        phir_u(:,iq2,iq1) = phir_u(:,iq1,iq2)
2373       end if ! sofswt == 1
2374 
2375     end do !loop on iq2
2376   end do !loop on iq1
2377 
2378 !---------my debug----------------------
2379    if ( sofswt == 1 ) then  
2380 !    open(unit=68,file='phir-gfil-RADSINTR-A.dat')
2381 !    open(unit=67,file='d2phidr2-gfil-RADSINTR-A.dat')
2382 !    do ir=0,nrpts
2383 !      write(68,*) dble(ir)*dr, ((phiraux(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts)
2384 !      write(67,*) dble(ir)*dr, ((d2phidr2(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts) 
2385 !    end do                                                                  
2386 !    close(unit=67)
2387 !    close(unit=68)
2388 
2389 !    open(unit=78,file='phig-gfil-RADSINTR-A.dat')                                 
2390 !    open(unit=77,file='d2phidg2-gfil-RADSINTR-A.dat')
2391 !    do ir=0,nrpts
2392 !      write(78,*) dble(ir)*dg, ((phig(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts)
2393 !      write(77,*) dble(ir)*dg, ((d2phidg2(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts)
2394 !    end do
2395 !    close(unit=77)
2396 !    close(unit=78)
2397 
2398 !    open(unit=87,file='phirfil-gfil-RADSINTR-A.dat')
2399 !    do ir=0,nrpts
2400 !      write(87,*) dble(ir)*dr, ((phir(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts)
2401 !    end do
2402 !    close(unit=87)
2403      
2404      ABI_FREE(phiraux)     
2405 
2406    end if
2407    
2408 !  if (sofswt == 0 ) then 
2409 !    open(unit=88,file='phir-uns-gfil-RADSINTR-A.dat')
2410 !    do ir=0,nrpts
2411 !       write(88,*) dble(ir)*dr, ((phir_u(ir,iq1,iq2),iq1=1,iq2),iq2=1,nqpts)
2412 !    end do
2413 !    close(unit=88)
2414 !  end if
2415 !----end my debug----------------
2416 
2417   ABI_FREE(utmp)
2418   ABI_FREE(q1sat)
2419   ABI_FREE(q2sat)
2420 
2421   DBG_EXIT("COLL")
2422 
2423 end subroutine vdw_df_filter

ABINIT/m_xc_vdw/vdw_df_indexof [ Functions ]

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

SOURCE

3063 function vdw_df_indexof(list_1d,npts,value)
3064 
3065 !Arguments ------------------------------------
3066   integer,intent(in) :: npts
3067   real(dp),intent(in) :: value,list_1d(npts)
3068 
3069 !Local variables ------------------------------
3070   integer :: imin,imax,imid
3071   integer :: vdw_df_indexof
3072 
3073 ! *************************************************************************
3074 
3075   imin = 1
3076   imax = npts
3077 
3078   do
3079     imid = (imin + imax) / 2
3080     if ( value < list_1d(imid) ) then
3081       imax = imid
3082     else
3083       imin = imid
3084     end if
3085     if ( imin + 1 >= imax ) exit
3086   end do
3087 
3088   vdw_df_indexof = imin
3089 
3090 end function vdw_df_indexof

ABINIT/m_xc_vdw/vdw_df_internal_checks [ Functions ]

[ Top ] [ Functions ]

NAME

  vdw_df_internal_checks

FUNCTION

  Performs consistency checks of the vdW-DF kernel.

INPUTS

  test_mode= bitfield to enable/disable tests

SOURCE

3181 subroutine vdw_df_internal_checks(test_mode)
3182 
3183 #if defined HAVE_FFTW3
3184   include "fftw3.f03"
3185 #endif
3186 
3187 !Arguments ------------------------------------
3188   integer,intent(in) :: test_mode
3189 
3190 !Local variables ------------------------------
3191   character(len=512) :: msg
3192   integer :: ndpts,ngpts,nqpts,nrpts,ntest,id1,ig,iq1,iq2,jj,ir,nn,ii
3193   integer :: iostatus,ip1,ir1,ir2,ir3,nr1,nr2,nr3,i1,i2,i3 
3194   integer(kind=SIZEOF_PTRDIFF_T) :: plan 
3195   real(dp) :: acutmin,aratio,damax,dsoft,phisoft
3196   real(dp) :: a1,a2,a3,delta_test,d1,d2i,gcut,gtmp,grho2,grhotmp2(3)
3197   real(dp) :: b1,b2,b3,dg,dr,rcut,rtmp,dvol,exc_tmp,deltae_vdw
3198   real(dp),allocatable :: gprimdt(:,:),kern_test_c(:),kern_test_i(:)
3199   real(dp),allocatable :: intg_calc(:),intg_int(:),ptmp(:,:,:) 
3200   real(dp),allocatable :: extin(:,:), grhotmp(:,:),ztmp(:,:) 
3201   real(dp),allocatable :: decdrho_tmp(:),decdgrho_tmp(:,:)
3202   real(dp),allocatable :: t3dr(:,:,:,:),t3dr1d(:,:),gvecs(:,:),t3dr2(:,:,:,:)
3203   real(dp),allocatable :: u1dtmp(:,:)
3204   real(dp) :: theta(my_vdw_params%nqpts,1,5),dummyv 
3205   complex(dp),allocatable :: t3dg(:,:,:,:),t3dg1d(:,:),utmp(:) 
3206 ! real(dp),allocatable :: t3dg(:,:,:,:),t3dg1d(:,:),utmp(:),u1dtmp(:,:)
3207 ! *************************************************************************
3208 
3209   DBG_ENTER("COLL")
3210 
3211   ! Here starts a test which will evaluate the integral over D for
3212   ! delta=0 (see Fig. 1 of DRSLL04).
3213   ! The integration will be performed over a finer and linear d mesh,
3214   ! not dmesh, for which the kernel values will be obtained by both
3215   ! interpolation and direct calculation. The range spanned by this mesh
3216   ! is half of dmesh range. Introduced variables: delta_test, kern_test_c,
3217   ! kern_test_i, ntest, intg_calc, intg_int.
3218   if ( iand(test_mode,1) /= 0 ) then
3219 
3220     write(msg,'(1x,a)') "[vdW-DF Check] Test #01: Kernel integral over D"
3221     call wrtout(std_out,msg,'COLL')
3222 
3223     acutmin = my_vdw_params%acutmin
3224     aratio = my_vdw_params%aratio
3225     damax = my_vdw_params%damax
3226     dsoft = my_vdw_params%dsoft
3227     ndpts = my_vdw_params%ndpts
3228     phisoft = my_vdw_params%phisoft
3229 
3230     ntest = 600
3231     delta_test = dmesh(ndpts) / (2 * ntest)
3232 
3233     ABI_MALLOC(kern_test_c,(ntest))
3234     ABI_MALLOC(kern_test_i,(ntest))
3235     ABI_MALLOC(intg_calc,(ntest))
3236     ABI_MALLOC(intg_int,(ntest))
3237 
3238 #if defined DEBUG_VERBOSE
3239     write(msg,'(a,4x,3(1x,a10),a)') ch10, "D","Calcul.","Interp.",ch10
3240     call wrtout(std_out,msg,'COLL')
3241 #endif
3242 
3243     do id1 = 1,ntest
3244       d1=id1*delta_test
3245       kern_test_c(id1) = (four_pi * d1**2) * &
3246 &     vdw_df_kernel_value(d1,d1,acutmin,aratio,damax)
3247       kern_test_i(id1) = (four_pi * d1**2) * vdw_df_interpolate(d1,d1,0)
3248 #if defined DEBUG_VERBOSE
3249       write(msg,'(4x,3(1x,e10.3))') id1*delta_test, &
3250 &       kern_test_c(id1), kern_test_i(id1)
3251       call wrtout(std_out,msg,'COLL')
3252 #endif
3253     end do
3254 
3255     call simpson_int(ntest,delta_test,kern_test_c,intg_calc)
3256     call simpson_int(ntest,delta_test,kern_test_i,intg_int)
3257 
3258     write(msg,'(a,4x,a,2(a,4x,a,e10.3))') ch10, &
3259 &     "> Integrals over D", ch10, ">   calculated   :",intg_calc(ntest),ch10, &
3260 &     ">   interpolated :",intg_int(ntest)
3261     call wrtout(std_out,msg,'COLL')
3262     if ( abs(intg_int(ntest) - intg_calc(ntest)) < tol3 ) then
3263       write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 PASSED"
3264     else
3265       write(msg,'(a,1x,a)') ch10,"[vdW-DF Check] Test #01 FAILED"
3266     end if
3267     call wrtout(std_out,msg,'COLL')
3268 
3269     ABI_FREE(kern_test_c)
3270     ABI_FREE(kern_test_i)
3271     ABI_FREE(intg_calc)
3272     ABI_FREE(intg_int)
3273 
3274 !   The following test aims at checking the interpolation in both 
3275 !   real and reciprocal space of the radial representation of the
3276 !   Kernel
3277 
3278     write(msg,'(1x,a)') "[vdW-DF Check] Test #02: Interpolation &
3279 &   of radial reciprocal space representation of the Kernel"   
3280     call wrtout(std_out,msg,'COLL')
3281 
3282 !   Interpolate phi in reciprocal space:
3283 !    * ptmp(:,:,1) = phi(g)
3284 !    * ptmp(:,:,2) = dphi(g)/dg
3285 !   The interolation will be performed in a three times 
3286 !   finer mesh than that used for the spline, and for 
3287 !   iq1=1q2=1,2 only.
3288 
3289     rcut = my_vdw_params%rcut
3290     nrpts = my_vdw_params%nrpts
3291     gcut = my_vdw_params%gcut 
3292     ngpts = my_vdw_params%ngpts
3293 !   begin my debug
3294     write(std_out,*) 'ngpts after init: ',ngpts
3295 !   end my debug
3296     ABI_MALLOC(ptmp,(2,2,2))
3297 
3298     dg = pi / rcut !gcut / (ngpts-one)
3299 #if defined DEBUG_VERBOSE
3300     write(msg,'(a,5x,5(1x,a10),a)') ch10, "g","ig","phi(g,1,1)","phi(g,1,2)" &
3301 &     ,"phi(g,2,2)", ch10
3302     call wrtout(std_out,msg,'COLL')
3303 #endif
3304     do jj = 0,3*ngpts-1 !3*(ngpts-1)
3305       ptmp(:,:,:) = zero
3306       gtmp = jj * dg / three
3307       ig = int(gtmp / dg) + 1
3308       a1 = ig - gtmp / dg  !((ig + 1) * dg - gtmp) / dg
3309       b1 = one - a1
3310       a2 = (3 * a1**2 - one) * dg / six
3311       b2 = (3 * b1**2 - one) * dg / six
3312       a3 = (a1**3 - a1) * dg**2 / six
3313       b3 = (b1**3 - b1) * dg**2 / six
3314       do iq2 = 1,2
3315         do iq1 = 1,iq2
3316 !          ptmp(iq1,iq2,1) = a1 * phig(ig,iq1,iq2) + b1 * phig(ig+1,iq1,iq2) &
3317 !&           + a3 * d2phidg2(ig,iq1,iq2) + b3 * d2phidg2(ig+1,iq1,iq2)
3318 !          ptmp(iq1,iq2,2) = (phig(ig+1,iq1,iq2) - phig(ig,iq1,iq2)) / dg &
3319 !&           - a2 * d2phidg2(ig,iq1,iq2) + b2 * d2phidg2(ig+1,iq1,iq2)
3320           ptmp(iq1,iq2,1) = a1 * phig(ig-1,iq1,iq2) + b1 * phig(ig,iq1,iq2) &
3321 &           + a3 * d2phidg2(ig-1,iq1,iq2) + b3 * d2phidg2(ig,iq1,iq2)
3322           ptmp(iq1,iq2,2) = (phig(ig,iq1,iq2) - phig(ig-1,iq1,iq2)) / dg &
3323 &           - a2 * d2phidg2(ig-1,iq1,iq2) + b2 * d2phidg2(ig,iq1,iq2)
3324         end do
3325       end do
3326   
3327       do iq2 = 1,2
3328         do iq1 = 1,iq2-1
3329           ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:)
3330         end do
3331       end do
3332 #if defined DEBUG_VERBOSE
3333       write(msg,'(5x,5(1x,f16.8))') gtmp, real(ig), &
3334 &       ptmp(1,1,1), ptmp(1,2,1), ptmp(2,2,1)
3335       call wrtout(std_out,msg,'COLL')
3336 #endif
3337     end do
3338 
3339     write(msg,'(1x,a)') "[vdW-DF Check] Test #02: Interpolation &
3340 &   of radial real space representation of the Kernel"   
3341     call wrtout(std_out,msg,'COLL')
3342 
3343 !     Interpolate phi in real space:
3344 !       * ptmp(:,:,1) = phi(r)
3345 !       * ptmp(:,:,2) = dphi(r)/dr
3346 !     The interolation will be performed in a ten times 
3347 !     finer mesh than that used for the spline and for 
3348 !     iq1=1q2=1,2.
3349 
3350     dr = rcut / nrpts !(nrpts-one)
3351 #if defined DEBUG_VERBOSE
3352     write(msg,'(a,4x,4(1x,a10),a)') ch10, "r","phi(r,1,1)","phi(r,1,2)" &
3353 &     ,"phi(r,2,2)", ch10
3354     call wrtout(std_out,msg,'COLL')
3355 #endif
3356     do jj = 0,3*nrpts-1 !-1
3357       ptmp(:,:,:) = zero
3358       rtmp = jj * dr / three
3359       ir = int(rtmp / dr) + 1
3360       a1 = ir  - ( rtmp / dr ) 
3361       b1 = one - a1
3362       a2 = (3 * a1**2 - one) * dr / six
3363       b2 = (3 * b1**2 - one) * dr / six
3364       a3 = (a1**3 - a1) * dr**2 / six
3365       b3 = (b1**3 - b1) * dr**2 / six
3366       do iq2 = 1,2
3367         do iq1 = 1,iq2
3368 !          ptmp(iq1,iq2,1) = a1 * phir(ir,iq1,iq2) + b1 * phir(ir+1,iq1,iq2) &
3369 !&           + a3 * d2phidr2(ir,iq1,iq2) + b3 * d2phidr2(ir+1,iq1,iq2)
3370 !          ptmp(iq1,iq2,2) = (phir(ir+1,iq1,iq2) - phir(ir,iq1,iq2)) / dr &
3371 !&           - a2 * d2phidr2(ir,iq1,iq2) + b2 * d2phidr2(ir+1,iq1,iq2)
3372           ptmp(iq1,iq2,1) = a1 * phir(ir-1,iq1,iq2) + b1 * phir(ir,iq1,iq2) &
3373 &           + a3 * d2phidr2(ir-1,iq1,iq2) + b3 * d2phidr2(ir,iq1,iq2)
3374           ptmp(iq1,iq2,2) = (phir(ir,iq1,iq2) - phir(ir-1,iq1,iq2)) / dr &
3375 &           - a2 * d2phidr2(ir-1,iq1,iq2) + b2 * d2phidr2(ir,iq1,iq2)
3376         end do
3377       end do
3378   
3379       do iq2 = 1,2
3380         do iq1 = 1,iq2-1
3381           ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:)
3382         end do
3383       end do
3384 #if defined DEBUG_VERBOSE
3385       write(msg,'(4x,4(1x,f16.8))') rtmp, &
3386 &       ptmp(1,1,1), ptmp(1,2,1), ptmp(2,2,1)
3387       call wrtout(std_out,msg,'COLL')
3388 #endif
3389     end do
3390 
3391     ABI_FREE(ptmp)
3392 
3393 !-------my debug------------------------------------------------
3394 !   The following test aims at checking the calculation of q0, 
3395 !   the interpolation  polynomials at q0, and the local cusp 
3396 !   correction from precomputed electronic density and its 
3397 !   gradient. Also theta and the nonlocal correction is
3398 !   computed. 
3399 
3400     write(msg,'(1x,a)') "[vdW-DF Check] Test #03: q0 and &
3401 &   interpolation polynomials from external rho and grho"   
3402     call wrtout(std_out,msg,'COLL')
3403 
3404     nqpts = my_vdw_params%nqpts
3405 
3406     open (unit=98, file='rhogrhoext.dat', status='old',action='read')
3407     open (unit=59, file='rho-grho2-q0-qpoly-check.dat',status='replace')
3408     open (unit=55, file='rho-epscusp-check.dat',status='replace')
3409     open (unit=60, file='rho-grho2-theta-check.dat',status='replace')
3410     open (unit=61, file='rho-grho2-dtdd-check.dat',status='replace')
3411 !    open (unit=62, file='rho-grho2-dtdgd-check.dat',status='replace')
3412  
3413     nn=0
3414     do
3415       read(98,*, iostat=iostatus) dummyv
3416       if(iostatus/=0) then ! to avoid end of file error.
3417         exit
3418       else
3419         nn=nn+1
3420       end if
3421     end do
3422     write(std_out,*) 'Number of lines in rhogrhoext.dat:', nn
3423 
3424     ABI_MALLOC(extin,(nn,9))
3425     ABI_MALLOC(grhotmp,(1,3))
3426     ABI_MALLOC(ztmp,(nqpts,nqpts)) 
3427     ABI_MALLOC(decdrho_tmp,(1))
3428     ABI_MALLOC(decdgrho_tmp,(3,1))
3429     ABI_MALLOC(gprimdt,(3,3))
3430 
3431     write(std_out,*) 'Allocated arrays...'
3432 
3433     ! extin contains rho, grho(1), grho(2), grho(3), ex, ec, vx, vc, q0
3434     rewind(98)
3435 
3436     write(std_out,*) 'Rewinded file'
3437 
3438     do ii=1,nn
3439       read(98,*) (extin(ii,jj),jj=1,9)
3440     end do
3441 
3442     write(std_out,*) 'Readed external file...'
3443     
3444     ! Pre-compute integrand for vdW energy correction 
3445     ztmp(:,:) = zero
3446     do ir=1,nrpts
3447       do iq1=1,nqpts
3448         do iq2=1,iq1
3449           ztmp(iq1,iq2) = ztmp(iq1,iq2) - &
3450   &       two * pi * dr * phir(ir,iq1,iq2) * (ir*dr)**2 
3451   !&       two * pi * ( phir_u(ir,iq1,iq2) - phir(ir,iq1,iq2) ) * dr
3452         end do
3453       end do
3454     end do
3455     do iq1=2,nqpts
3456       do iq2=1,iq1-1
3457          ztmp(iq2,iq1) = ztmp(iq1,iq2)
3458       end do
3459     end do
3460 !---------my debug----------------------------------
3461     open(34,file='table-localcorr.dat',status='replace')
3462     do iq2=1,nqpts
3463        do iq1=1,nqpts
3464          write(34,*) iq1, iq2, ztmp(iq1,iq2)
3465        enddo
3466     enddo
3467     close(34)
3468     write(std_out,*) 'Table file written...'
3469 
3470 !--------end my debug------------------------------  
3471 
3472     ! Compute q0, polynomials, theta from Abinit routines
3473     ! Also get nonlocal correction to the energy
3474     open(unit=81,file='gvecs.dat',status='replace')
3475     theta(:,:,:) = zero
3476     deltae_vdw = zero
3477     dvol = 0.025742973 ! this is a particular value for the 
3478                        ! precomputed density. Same for gprimd
3479     
3480     gprimdt(1,1)= 0.22166114341982476
3481     gprimdt(2,1)= 0.0d0
3482     gprimdt(3,1)= 0.0d0
3483     gprimdt(1,2)= 0.0d0
3484     gprimdt(2,2)= 0.33249171512973713
3485     gprimdt(3,2)= 0.0d0
3486     gprimdt(1,3)= 0.0d0
3487     gprimdt(2,3)= 0.0d0
3488     gprimdt(3,3)= 0.33249171512973713
3489 
3490     nr1=96  !particular values for this precomputed density
3491     nr2=64
3492     nr3=64
3493     
3494     ABI_MALLOC(t3dr,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts))
3495     ABI_MALLOC(t3dr2,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts))
3496     ABI_MALLOC(t3dg,(0:nr1-1,0:nr2-1,0:nr3-1,nqpts)) 
3497     ABI_MALLOC(t3dr1d,(nqpts,nr1*nr2*nr3))
3498     ABI_MALLOC(t3dg1d,(nqpts,nr1*nr2*nr3))
3499     ABI_MALLOC(gvecs,(3,nr1*nr2*nr3))
3500     ABI_MALLOC(utmp,(nqpts))
3501     ABI_MALLOC(u1dtmp,(nqpts,nr1*nr2*nr3))
3502      
3503     ip1 = 1 
3504     do ir3=0,nr3-1 !1,nr3
3505       do ir2=0,nr2-1 !1,nr2
3506         do ir1=0,nr1-1 !1,nr1        
3507           decdrho_tmp(:) = zero                                    
3508           decdgrho_tmp(:,:) = zero
3509           grhotmp(1,1:3) = extin(ip1,2:4)
3510           forall(ii=1:3) grhotmp2(ii) = grhotmp(1,ii) 
3511           grho2 = sum(grhotmp2**2) 
3512           !write(*,*) 'Call to xc_vdw_energy number:', ip1
3513           call xc_vdw_energy(1,extin(ip1,1),grhotmp,extin(ip1,5),&
3514    &      extin(ip1,6),extin(ip1,7),extin(ip1,8),theta,&
3515    &      exc_tmp,decdrho_tmp,decdgrho_tmp,ztmp)
3516           t3dr(ir1,ir2,ir3,1:nqpts) = theta(1:nqpts,1,1)
3517           deltae_vdw = deltae_vdw + extin(ip1,1) * exc_tmp * dvol
3518 !          forall(ii=1:3) grhotmp2(ii) = grhotmp(1,ii) 
3519 !          grho2 = sum(grhotmp2**2) 
3520 !      exc_vdw = exc_vdw + rho_tmp * &
3521 !&     ( half * ( sum(ttmp2(:,ip1) * theta(:,1,1)) / &
3522 !&       (rho_tmp + tiny(rho_tmp)) ) * dvol )  
3523           if (ir1 > nr1/2) then 
3524             i1 = ir1 - nr1
3525           else
3526             i1 = ir1 
3527           end if
3528           if (ir2 > nr2/2) then
3529             i2 = ir2 - nr2
3530           else
3531             i2 = ir2
3532           end if  
3533           if (ir3 > nr3/2) then
3534             i3 = ir3 - nr3
3535           else
3536             i3 = ir3
3537           end if 
3538           gvecs(:,ip1) =  gprimdt(:,1) * i1 + gprimdt(:,2)&
3539 &                         * i2 + gprimdt(:,3) * i3
3540           write(81,*) ip1, gvecs(:,ip1)
3541           write(55,*) extin(ip1,1), exc_tmp
3542 !          write(60,*) extin(ip1,1), grho2, (theta(iq1,1,1),iq1=1,nqpts)
3543 !          write(61,*) extin(ip1,1), grho2, (theta(iq1,1,2),iq1=1,nqpts)
3544 !          write(62,*) extin(ip1,1), grho2, &
3545 !&         (theta(iq1,1,3),theta(iq1,1,4),theta(iq1,1,5),iq1=1,nqpts)
3546           ip1 = ip1 +1
3547         end do 
3548       end do
3549     end do
3550     write(std_out,*) 'Cusp correction for this density:'
3551     write(std_out,*) 'deltae_vdw=',deltae_vdw
3552     close(unit=81)
3553     close(unit=98)
3554     close(unit=59)
3555     close(unit=55) 
3556     close(unit=60)
3557     close(unit=61) 
3558 !    close(unit=62)
3559 !   Evaluation of 3D FFT of a component of theta
3560 !   open (unit=63, file='thetar-3D-iq20.dat', status='old',action='read')
3561 !   open (unit=63, file='R-thetag-1D-check.dat', status='replace')
3562 !   open (unit=64, file='I-thetag-1D-check.dat', status='replace')
3563     open (unit=65, file='thetar-1D-check.dat', status='replace')
3564 !    open (unit=66, file='thetar-1D-after-check.dat', status='replace')
3565 !   nn=0
3566 !   do
3567 !     read(63,*, iostat=iostatus) dummyv
3568 !     if(iostatus/=0) then ! to avoid end of file error.
3569 !       exit
3570 !     else
3571 !       nn=nn+1
3572 !     end if
3573 !   end do
3574 !   rewind(63)
3575 !   write(*,*) 'Rewinded file'
3576 
3577 !   write(*,*) 'Number of lines in thetar-3D-iq20.dat:', nn
3578 !   nr1=96  !particular values for this precomputed density
3579 !   nr2=64
3580 !   nr3=64
3581 !   ABI_MALLOC(t3dr,(nr1,nr2,nr3,1))
3582 !   ABI_MALLOC(t3dg,(nr1,nr2,nr3,1)) 
3583 !   ABI_MALLOC(t3dg1d,(nr1*nr2*nr3))
3584 !   ABI_MALLOC(gvecs,(3,nr1*nr2*nr3))
3585 !   do ir3=1,nr3
3586 !     do ir2=1,nr2
3587 !       do ir1=1,nr1
3588 !        read(63,*) t3dr(ir1,ir2,ir3,1)
3589 !       end do 
3590 !     end do 
3591 !   end do
3592 
3593   ! Fourier-transform theta
3594 #if defined HAVE_FFTW3
3595     do iq1 = 1,nqpts
3596       call dfftw_plan_dft_r2c_3d(plan,nr1,nr2,nr3, &
3597   &   t3dr(:,:,:,iq1),t3dg(:,:,:,iq1),FFTW_ESTIMATE)
3598       call dfftw_execute(plan)
3599       call dfftw_destroy_plan(plan)
3600       t3dg(:,:,:,iq1) = t3dg(:,:,:,iq1) / (nr1 * nr2 * nr3)
3601     end do
3602 !      write(*,*) 'FFTW3 call number', iq1
3603 !fftw_plan_r2r_3d
3604 #else
3605     ABI_ERROR('vdW-DF calculations require FFTW3 support')
3606 #endif
3607     ! Repack theta
3608     ip1 = 1
3609     do ir3=0,nr3-1 !1,nr3
3610       do ir2=0,nr2-1 !1,nr2
3611         do ir1=0,nr1-1 !1,nr1
3612           t3dg1d(1:nqpts,ip1) = t3dg(ir1,ir2,ir3,1:nqpts)
3613           t3dr1d(1:nqpts,ip1) = t3dr(ir1,ir2,ir3,1:nqpts) 
3614 !          write(63,*) ip1, (real(t3dg1d(iq1,ip1)),iq1=1,nqpts)
3615 !          write(64,*) ip1, (aimag(t3dg1d(iq1,ip1)),iq1=1,nqpts)
3616 !         write(65,*) ip1, (t3dr1d(iq1,ip1),iq1=1,nqpts)
3617           ip1 = ip1 + 1
3618         end do
3619       end do
3620     end do
3621     write(std_out,*) 'Finished 1D repack of theta.'
3622 
3623     ! Fourier-transform back theta just to check fftw3
3624 !#if defined HAVE_FFT_FFTW3
3625 !    write(*,*) "Backward FFT of theta - test"
3626 !    do iq1=1,nqpts
3627 !     call dfftw_plan_dft_c2r_3d(plan,nr1,nr2,nr3, &
3628 !&        t3dg(:,:,:,iq1),t3dr2(:,:,:,iq1),FFTW_ESTIMATE)
3629 !     call dfftw_execute(plan)
3630 !     call dfftw_destroy_plan(plan)
3631 !    end do
3632 !    ip1 = 1
3633 !    do ir3=0,nr3-1 !1,nr3
3634 !      do ir2=0,nr2-1 !1,nr2
3635 !        do ir1=0,nr1-1 !1,nr1
3636 !          write(66,*) ip1, (t3dr2(ir1,ir2,ir3,iq1),iq1=1,nqpts)
3637 !          ip1 = ip1 + 1
3638 !        end do
3639 !      end do
3640 !    end do
3641 !#else
3642 !    MSG_ERROR('vdW-DF calculations require FFTW3 support')
3643 !#endif
3644 
3645     open(unit=82,file='Interp-phi-check.dat',status='replace') 
3646    ! Reset 3D counters, since theta will be reconstructed in 3D on the fly
3647     ir1 = 0 !1
3648     ir2 = 0 !1
3649     ir3 = 0 !1
3650     ! Go through reciprocal vectors to build u
3651     ABI_MALLOC(ptmp,(nqpts,nqpts,2))
3652     do ip1=1,nr1*nr2*nr3
3653       gtmp = sqrt(sum(gvecs(:,ip1)**2))  
3654     ! Interpolate kernel in reciprocal space
3655       if ( gtmp < gcut ) then
3656 
3657       ! Interpolate phi in reciprocal space:
3658       !   * ptmp(:,:,1) = phi(g)
3659       !   * ptmp(:,:,2) = dphi(g)/dg
3660       ! Note: do this on the fly, to go as fast as possible (this is equivalent
3661       !       to a call of the 'splint' routine)
3662         ptmp(:,:,:) = zero
3663         dg = pi / my_vdw_params%rcut
3664         ig = int(gtmp / dg) + 1
3665         a1 = ig  - gtmp / dg
3666         b1 = one - a1
3667         a2 = (3 * a1**2 - one) * dg / six
3668         b2 = (3 * b1**2 - one) * dg / six
3669         a3 = (a1**3 - a1) * dg**2 / six
3670         b3 = (b1**3 - b1) * dg**2 / six
3671         do iq2 = 1,nqpts
3672           do iq1 = 1,iq2
3673             ptmp(iq1,iq2,1) = a1 * phig(ig-1,iq1,iq2) + b1 * phig(ig,iq1,iq2) &
3674 &           + a3 * d2phidg2(ig-1,iq1,iq2) + b3 * d2phidg2(ig,iq1,iq2)
3675             ptmp(iq1,iq2,2) = (phig(ig,iq1,iq2) - phig(ig-1,iq1,iq2)) / dg &
3676 &           - a2 * d2phidg2(ig-1,iq1,iq2) + b2 * d2phidg2(ig,iq1,iq2)
3677           end do
3678         end do
3679 
3680         do iq2 = 1,nqpts
3681           do iq1 = 1,iq2-1
3682             ptmp(iq2,iq1,:) = ptmp(iq1,iq2,:)
3683           end do
3684         end do
3685 !---------my debug------------------------------------------------------------
3686 ! write diagonal components of interpolated kernel at each point
3687         write(82,*) ip1, (ptmp(iq1,iq1,1),iq1=1,nqpts)
3688 !--------end my debug--------------------------------------------------------
3689         ! Calculate contributions to integral in Fourier space:
3690         ! Eq(12) from RS09 
3691         utmp(:) = matmul(t3dg1d(:,ip1),ptmp(:,:,1))                      
3692       else       
3693 
3694         write(std_out,*) 'WARNING: gtmp > gcut!!! This is at'
3695         write(std_out,*) 'ip1=',ip1
3696         utmp(:) = (zero,zero)
3697 
3698       end if ! gtmp < gcut
3699 
3700       ! Reconstruct the integrand in 3D
3701       t3dg(ir1,ir2,ir3,:) = utmp(:)
3702                                         
3703       ir1 = ir1 + 1
3704       if ( ir1 > nr1-1 ) then
3705         ir2 = ir2 + 1
3706         ir1 = 0
3707       end if
3708       if ( ir2 > nr2-1 ) then
3709         ir3 = ir3 + 1
3710         ir2 = 0
3711       end if
3712     end do !ip1=1,npts_rho
3713 
3714     ! Fourier-transform back the integrand
3715 #if defined HAVE_FFTW3
3716     do iq1=1,nqpts
3717       call dfftw_plan_dft_c2r_3d(plan,nr1,nr2,nr3, &
3718        t3dg(:,:,:,iq1),t3dr(:,:,:,iq1),FFTW_ESTIMATE)
3719       call dfftw_execute(plan)
3720       call dfftw_destroy_plan(plan)
3721     end do
3722 #else
3723     ABI_ERROR('vdW-DF calculations require FFTW3 support')
3724 #endif
3725    ! Write u for comparison
3726     open (unit=80, file='ureal1d-check.dat', status='replace')
3727    ! Repack the integrand
3728     ip1 = 1
3729 !$OMP  PARALLEL DO COLLAPSE(3) &
3730 !$OMP& DEFAULT(SHARED) PRIVATE(ir1,ir2,ir3)
3731     do ir3=0,nr3-1
3732       do ir2=0,nr2-1
3733         do ir1=0,nr1-1
3734           u1dtmp(:,ip1) = t3dr(ir1,ir2,ir3,1:nqpts) !real space u
3735           write(80,*) ip1, (u1dtmp(iq1,ip1),iq1=1,nqpts)
3736           ip1 = ip1 + 1
3737         end do
3738       end do
3739     end do
3740    ! ttmp2(:,ip1) corresponds to u_alpha,i at eq(11)
3741    ! from RS09. 
3742 !$OMP END PARALLEL DO
3743 
3744 
3745 
3746 !  close(unit=63)
3747 !  close(unit=64)
3748    close(unit=65)
3749 !    close(unit=66)
3750     close(unit=80)
3751     close(unit=82)
3752     ABI_FREE(extin) 
3753     ABI_FREE(grhotmp) 
3754     ABI_FREE(ztmp)
3755     ABI_FREE(ptmp)
3756     ABI_FREE(decdrho_tmp)
3757     ABI_FREE(decdgrho_tmp)
3758     ABI_FREE(t3dr)
3759     ABI_FREE(t3dr2)
3760     ABI_FREE(t3dg)
3761     ABI_FREE(t3dr1d)
3762     ABI_FREE(t3dg1d)
3763     ABI_FREE(gprimdt) 
3764     ABI_FREE(gvecs)
3765     ABI_FREE(utmp)
3766     ABI_FREE(u1dtmp)
3767 !--------end my debug--------------------------------------------
3768 
3769   end if
3770 
3771   DBG_EXIT("COLL")
3772 
3773 end subroutine vdw_df_internal_checks

ABINIT/m_xc_vdw/vdw_df_interpolate [ Functions ]

[ 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

SOURCE

3107 function vdw_df_interpolate(d1,d2,sofswt)
3108 
3109   use m_errors
3110 
3111 !Arguments ------------------------------------
3112   real(dp),intent(in) :: d1,d2
3113   integer,intent(in) :: sofswt
3114 
3115 !Local variables ------------------------------
3116   real(dp) :: vdw_df_interpolate
3117   integer :: id1,id2,ix1,ix2,ndpts
3118   real(dp) :: dcut,xd(4,4)
3119 
3120 ! *************************************************************************
3121 
3122   vdw_df_interpolate = zero
3123 
3124   ndpts = size(dmesh(:))
3125   dcut = dmesh(ndpts)
3126 
3127   if ( (d1 >= dcut) .or. (d2 >= dcut) ) then
3128     return
3129   end if
3130 
3131   ! Find closest indices in dmesh
3132   id1 = vdw_df_indexof(dmesh,ndpts,d1)
3133   id2 = vdw_df_indexof(dmesh,ndpts,d2)
3134 
3135   ! Find intervals
3136   xd(1,:) = one
3137   xd(2,1) = (d1 - dmesh(id1)) / (dmesh(id1+1) - dmesh(id1))
3138   xd(2,2) = (d2 - dmesh(id2)) / (dmesh(id2+1) - dmesh(id2))
3139   xd(3,:) = xd(2,:)**2
3140   xd(4,:) = xd(2,:)**3
3141 
3142 
3143   select case ( sofswt )
3144 
3145     case (1)
3146       do ix2=1,4
3147         do ix1=1,4
3148           vdw_df_interpolate = vdw_df_interpolate + &
3149 &           phi_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2)
3150         end do
3151       end do
3152 
3153     case (0)
3154       do ix1=1,4
3155         do ix2=1,4
3156           vdw_df_interpolate = vdw_df_interpolate + &
3157 &           phi_u_bicubic(ix1,ix2,id1,id2) * xd(ix1,1) * xd(ix2,2)
3158         end do
3159       end do
3160 
3161     case default
3162       ABI_ERROR("  vdW-DF soft switch must be 0 or 1")
3163 
3164   end select
3165 
3166 end function vdw_df_interpolate

ABINIT/m_xc_vdw/vdw_df_kernel [ Functions ]

[ 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

SOURCE

2450 function vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax)
2451 
2452 !Arguments ------------------------------------
2453   real(dp),intent(in) :: d1,d2,dsoft,phisoft,acutmin,aratio,damax
2454 
2455 !Local variables ------------------------------
2456   real(dp) :: vdw_df_kernel
2457   character(len=512) :: msg
2458   integer :: napts
2459   real(dp) :: deltad,dtol,dtmp,d1m,d1p,d2m,d2p,phid,phim,phip,phix
2460 
2461 ! *************************************************************************
2462 
2463 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2464   DBG_ENTER("COLL")
2465 #endif
2466 
2467   ! Init
2468   dtol = my_vdw_params%tolerance
2469   dtmp = sqrt(d1**2 + d2**2)
2470 
2471   if ( dtmp < dsoft ) then
2472 
2473     ! Calculate a softened version of phi when (d1,d2) -> (0,0)
2474     ! using a polynomial expansion for nonzero distances and an
2475     ! exponential fit on dsoft when d -> 0 too.
2476     ! Note 1: the exponential fit is inspired from the definition
2477     !         of phi0_soft in RS09.
2478     ! Note 2: contrary to RS09, deltad is proportional to dsoft.
2479 
2480     if ( phisoft < 0 ) then
2481       ABI_ERROR('the softened kernel must be positive')
2482     end if
2483 
2484     vdw_df_kernel = phisoft
2485 
2486     if ( dtmp > zero ) then
2487       deltad = dsoft / 100.0_dp
2488       d1m = (dsoft - deltad) * d1 / dtmp
2489       d1p = (dsoft + deltad) * d1 / dtmp
2490       d2m = (dsoft - deltad) * d2 / dtmp
2491       d2p = (dsoft + deltad) * d2 / dtmp
2492 
2493       phim = vdw_df_kernel_value(d1m,d2m,acutmin,aratio,damax)
2494       phip = vdw_df_kernel_value(d1p,d2p,acutmin,aratio,damax)
2495       phix = (phim + phip) / two
2496       phid = (phip - phim) / (two * deltad)
2497 
2498       vdw_df_kernel = phisoft + &
2499 &       (four * (phix - phisoft) - phid * dsoft) &
2500 &       * dtmp**2 / (two * dsoft**2) + &
2501 &       (two  * (phisoft - phix) + phid * dsoft) &
2502 &       * dtmp**4 / (two * dsoft**4)
2503     end if
2504 
2505   else
2506 
2507     vdw_df_kernel = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)
2508 
2509   end if ! dtmp < dsoft
2510 
2511 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2512   DBG_EXIT("COLL")
2513 #endif
2514 
2515 end function vdw_df_kernel

ABINIT/m_xc_vdw/vdw_df_kernel_value [ Functions ]

[ 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

SOURCE

2535 function vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)
2536 
2537 !Arguments ------------------------------------
2538   real(dp),intent(in) :: d1,d2,acutmin,aratio,damax
2539 
2540 !Local variables ------------------------------
2541   real(dp) :: vdw_df_kernel_value
2542   character(len=512) :: msg
2543   integer :: ia,ib,napts,metyp
2544   real(dp) :: aa,bb,atol,atmp,acut,btmp,dain,damin,tau,ttmp,wtmp
2545   real(dp),allocatable :: amesh(:),amesh_cos(:),amesh_sin(:),nu1(:),nu2(:)
2546   real(dp),allocatable :: dphida(:),dphidb(:),dphidd(:),ycoef(:,:),work(:)
2547   real(dp),allocatable :: eint(:)
2548 ! *************************************************************************
2549 
2550 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2551   DBG_ENTER("COLL")
2552 #endif
2553 
2554   ! Create angular mesh
2555   !
2556   ! Note: Contrary to RS09, we use a linear mesh, but leave
2557   !       the door open to other kinds of meshes. Indeed, the use of a
2558   !       logarithmic mesh is interesting only when both d1 and d2 are
2559   !       small, in which case the kernel is usually softened.
2560 
2561   ! Init
2562 
2563 !!  Testing the same way in which amesh is obtained in siesta:
2564 !!  napts = nint(max(acutmin,aratio*sqrt(d1**2+d2**2))) This is
2565 !!  very different to the way the number of a-mesh points are determined
2566 !!  in Siesta.  Actually this almost corresponds to their definition of acut
2567 
2568   acut = max(acutmin,aratio*sqrt(d1**2+d2**2))
2569 
2570 !! Obtaining the number of amesh points, it is introduced damin
2571 !! and dain:
2572   damin = 0.01d0
2573   dain = min(damax,0.1*sqrt(d1**2+d2**2))
2574   dain = max(dain,damin)
2575 
2576   if ( dain == damax ) then !Linear mesh
2577    metyp = 2
2578    napts = nint( acut / damax )
2579   else
2580    bb = acut * dain / (damax-dain)
2581    aa = log( 1 + dain/bb )
2582    napts = 1 + nint( log(1+acut/bb) / aa )
2583    metyp = 1
2584   end if
2585 !! getting napts end
2586 
2587   atol = my_vdw_params%tolerance
2588 
2589 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2590   write(msg,'(1x,"vdw_df_kernel_value:",1x,"napts = ",i8)') napts
2591   call wrtout(std_out,msg,'COLL')
2592 #endif
2593 
2594   ABI_MALLOC(amesh,(napts))
2595   ABI_MALLOC(amesh_cos,(napts))
2596   ABI_MALLOC(amesh_sin,(napts))
2597   ABI_MALLOC(dphida,(napts))
2598   ABI_MALLOC(dphidb,(napts))
2599   ABI_MALLOC(nu1,(napts))
2600   ABI_MALLOC(nu2,(napts))
2601   ABI_MALLOC(ycoef,(3,napts))
2602   ABI_MALLOC(work,(napts))
2603   ABI_MALLOC(eint,(napts))
2604   ! Init amesh
2605   !FIXME: allow for other mesh types
2606   !Now intruducing metyp, setting mesh_cutoff to acut, removing
2607   !mesh_inc=damax and including mesh_ratio=
2608   call vdw_df_create_mesh(amesh,napts,metyp,acut,mesh_ratio=damax/dain)
2609 
2610   ! Init cos(a), sin(a), nu1(a), and nu2(a)
2611   tau = 4 * pi / 9
2612 !$OMP  PARALLEL DO &
2613 !$OMP& IF ( napts > 511 ) &
2614 !$OMP& DEFAULT(SHARED) PRIVATE(ia)
2615   do ia = 1,napts
2616     atmp = amesh(ia)
2617     amesh_cos(ia) = cos(atmp)
2618     amesh_sin(ia) = sin(atmp)
2619     if ( atmp == zero ) then  !it was changed the original condition: atmp < atol
2620       nu1(ia) = d1**2 / 2 / tau
2621       nu2(ia) = d2**2 / 2 / tau
2622     else
2623       if ( d1 == zero ) then  !it was changed the original condition: d1 < atol
2624         nu1(ia) = atmp*atmp / 2
2625       else
2626         nu1(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d1)**2))
2627       end if
2628       if ( d2 == zero ) then  !it was changed the original condition: d2 < atol
2629         nu2(ia) = atmp*atmp / 2
2630       else
2631         nu2(ia) = atmp*atmp / 2 / (1 - exp(-tau*(atmp/d2)**2))
2632       end if
2633     end if
2634   end do
2635 !$OMP END PARALLEL DO
2636 
2637   ! Integrate
2638   dphida(1) = 0
2639 !$OMP PARALLEL DO &
2640 !$OMP& IF ( napts > 255 ) &
2641 !$OMP& DEFAULT(SHARED) PRIVATE(ia,ib,atmp,btmp,ttmp,wtmp)
2642   do ia = 2,napts
2643     atmp = amesh(ia)
2644     dphidb(1) = 0
2645 
2646     do ib = 2,napts
2647       btmp = amesh(ib)
2648 
2649       ! eq.(15) DRSLL04
2650       ttmp = half * ( 1 / (nu1(ia) + nu1(ib)) + 1 / (nu2(ia) + nu2(ib)) ) * &
2651 &       ( 1 / (nu1(ia) + nu2(ia)) / (nu1(ib) + nu2(ib)) + &
2652 &         1 / (nu1(ia) + nu2(ib)) / (nu2(ia) + nu1(ib)) )
2653 
2654       ! eq.(16) DRSLL04
2655       wtmp = two * ( (three - atmp*atmp) * btmp * amesh_cos(ib) * &
2656 &       amesh_sin(ia) + (three - btmp*btmp) * atmp * amesh_cos(ia) * &
2657 &       amesh_sin(ib) + (atmp*atmp + btmp*btmp - three) * &
2658 &       amesh_sin(ia) * amesh_sin(ib) - &
2659 &       three * atmp * btmp * amesh_cos(ia) * amesh_cos(ib) ) / &
2660 &       (atmp * btmp)**3
2661 
2662       dphidb(ib) = atmp*atmp * btmp*btmp * wtmp * ttmp
2663 
2664     end do ! ib = 2,napts
2665 
2666 !    call spline_integrate(dphida(ia),napts,damax,dphidb)
2667      call cspint(dphidb,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,dphida(ia))
2668   end do ! ia = 2,napts
2669 !$OMP END PARALLEL DO
2670 
2671   ! Final value of the kernel
2672 !  call spline_integrate(vdw_df_kernel_value,napts,damax,dphida)
2673      call cspint(dphida,amesh,napts,amesh(1),amesh(napts),ycoef,eint,work,vdw_df_kernel_value)
2674   vdw_df_kernel_value = vdw_df_kernel_value * 2 / pi**2
2675 
2676   ! Clean-up the mess
2677   ABI_FREE(amesh)
2678   ABI_FREE(amesh_cos)
2679   ABI_FREE(amesh_sin)
2680   ABI_FREE(nu1)
2681   ABI_FREE(nu2)
2682   ABI_FREE(dphida)
2683   ABI_FREE(dphidb)
2684   ABI_FREE(ycoef)
2685   ABI_FREE(work)
2686   ABI_FREE(eint)
2687 
2688 #if ( defined DEBUG_VERBOSE ) && ( DEBUG_VERBOSE > 1 )
2689   DBG_EXIT("COLL")
2690 #endif
2691 
2692 end function vdw_df_kernel_value

ABINIT/m_xc_vdw/vdw_df_ldaxc [ Functions ]

[ 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

SOURCE

2717 subroutine vdw_df_ldaxc(npts_rho,nspden,ngrad,rho_grho, &
2718 &                       exc_lda,vxc_lda,vxcg_lda)
2719 
2720 !Arguments ------------------------------------
2721   integer,intent(in) :: nspden,npts_rho,ngrad
2722   real(dp),intent(in) :: rho_grho(npts_rho,nspden,ngrad)
2723   real(dp),intent(out) :: exc_lda(2,npts_rho),vxc_lda(2,npts_rho,nspden)
2724   real(dp),intent(out) :: vxcg_lda(2,3,npts_rho)
2725 
2726 !Local variables ------------------------------
2727   integer :: ii
2728   logical :: is_gga
2729   type(libxc_functional_type) :: aux_funcx(2),aux_funcc(2)
2730   real(dp),allocatable :: rho_tmp(:,:)!,grho2(:,:)
2731 
2732 ! *************************************************************************
2733 
2734   DBG_ENTER("COLL")
2735 
2736 ! is_gga=libxc_functionals_isgga(vdw_funcs)
2737 
2738   ABI_MALLOC(rho_tmp,(npts_rho,nspden))
2739 ! if (is_gga) then
2740 !   ABI_MALLOC(grho2,(npts_rho,2*min(nspden,2)-1))
2741 ! end if
2742 
2743   ! Convert the quantities provided by ABINIT to the ones needed by LibXC
2744   !
2745   ! Notes:
2746   !
2747   !   * Abinit passes rho_up in the spin-unpolarized case, while LibXC
2748   !     expects the total density, but as we use rhonow, we don't need
2749   !     to convert anything.
2750   !
2751   !   * In the case off gradients:
2752   !
2753   !       - Abinit passes |grho_up|^2 while LibXC needs |grho_tot|^2;
2754   !       - Abinit passes |grho_up|^2, |grho_dn|^2, and |grho_tot|^2,
2755   !         while LibXC needs |grho_up|^2, grho_up.grho_dn,
2756   !         and |grho_dn|^2;
2757   !
2758   !     but again the use of rho_grho lets us build these quantities
2759   !     directly.
2760   !
2761   ! See ~abinit/src/56_xc/rhotoxc.F90 for details.
2762   ! However in routine libxc_functionals_getvxc of m_libxc_functionals the transformation
2763   ! to libxc input takes place, this means that we do not need to provide here the
2764   ! quatities for libxc.   
2765 
2766   if (nspden==1) then
2767     do ii=1,npts_rho
2768       rho_tmp(ii,1)=half*rho_grho(ii,1,1)
2769 !     if (is_gga) grho2(ii,1)=quarter*(rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2 &
2770 !&                                     +rho_grho(ii,1,4)**2)
2771     end do
2772   else
2773     do ii=1,npts_rho
2774       rho_tmp(ii,1)=rho_grho(ii,2,1)
2775       rho_tmp(ii,2)=rho_grho(ii,1,1)-rho_grho(ii,2,1)
2776 !     if (is_gga) then
2777 !       grho2(ii,1)=rho_grho(ii,2,2)**2+rho_grho(ii,2,3)**2+rho_grho(ii,2,4)**2
2778 !       grho2(ii,2)=(rho_grho(ii,1,2)-rho_grho(ii,2,2))**2 &
2779 !                  +(rho_grho(ii,1,3)-rho_grho(ii,2,3))**2 &
2780 !                  +(rho_grho(ii,1,4)-rho_grho(ii,2,4))**2
2781 !       grho2(ii,3)=rho_grho(ii,1,2)**2+rho_grho(ii,1,3)**2+rho_grho(ii,1,4)**2
2782 !     end if
2783     end do
2784   end if
2785 
2786 !  if (is_gga) then
2787 !    call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp,exc_lda,vxc_lda,&
2788 !&                                 grho2=grho2,vxcgr=vxcg_lda,xc_functionals=vdw_funcs)
2789 !  else
2790 
2791 !  Get LDA exchange energy and potential
2792 
2793     call libxc_functionals_init(-001000,1,xc_functionals=aux_funcx)
2794 
2795     write(std_out,*) 'aux_funcx=', aux_funcx(1)%id, aux_funcx(2)%id
2796     is_gga=libxc_functionals_isgga(aux_funcx)
2797     write(std_out,*) 'LDA exchange is GGA?', is_gga
2798 
2799     call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp(:,nspden),exc_lda(1,:),vxc_lda(1,:,nspden),&
2800 &                                 xc_functionals=aux_funcx)
2801 
2802     call libxc_functionals_end(xc_functionals=aux_funcx)
2803 
2804 !  Get LDA correlation energy and potential
2805 
2806     call libxc_functionals_init(-000012,1,xc_functionals=aux_funcc)
2807 
2808     write(std_out,*) 'aux_funcc=', aux_funcc(1)%id, aux_funcc(2)%id
2809     is_gga=libxc_functionals_isgga(aux_funcc)
2810     write(std_out,*) 'LDA correlation is GGA?', is_gga
2811 
2812     call libxc_functionals_getvxc(1,1,npts_rho,nspden,1,rho_tmp(:,nspden),exc_lda(2,:),vxc_lda(2,:,nspden),&
2813 &                                 xc_functionals=aux_funcc)
2814 
2815 
2816     call libxc_functionals_end(xc_functionals=aux_funcc)
2817 
2818 !  end if
2819 
2820   ABI_FREE(rho_tmp)
2821 !  if (is_gga) then
2822 !    ABI_FREE(grho2)
2823 !  end if
2824 
2825   DBG_EXIT("COLL")
2826 
2827 end subroutine vdw_df_ldaxc

ABINIT/m_xc_vdw/vdw_df_set_tweaks [ Functions ]

[ 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

SOURCE

3894 subroutine vdw_df_set_tweaks(input_tweaks,output_tweaks)
3895 
3896 !Arguments ------------------------------------
3897   integer,intent(in) :: input_tweaks
3898   type(vdw_df_tweaks_type),intent(out) :: output_tweaks
3899 
3900 !Local variables-------------------------------
3901 
3902 ! *************************************************************************
3903 
3904   output_tweaks%amesh_type = iand(ishft(input_tweaks,-10),3)
3905   output_tweaks%deriv_type = iand(ishft(input_tweaks,-2),3)
3906   output_tweaks%dmesh_type = iand(ishft(input_tweaks,-6),3)
3907   output_tweaks%interp_type = iand(ishft(input_tweaks,-4),3)
3908   output_tweaks%qmesh_type = iand(ishft(input_tweaks,-8),3)
3909   output_tweaks%run_type = iand(input_tweaks,3)
3910 
3911 end subroutine vdw_df_set_tweaks

ABINIT/m_xc_vdw/vdw_df_tweaks_type [ Types ]

[ Top ] [ Types ]

NAME

  vdw_df_tweaks_type

FUNCTION

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

SOURCE

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

ABINIT/m_xc_vdw/vdw_df_write_func [ Functions ]

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

SOURCE

3927 subroutine vdw_df_write_func(func_name,mode)
3928 
3929 !Arguments ------------------------------------
3930   character(len=*),intent(in) :: func_name,mode
3931 
3932 !Local variables-------------------------------
3933   character(len=512) :: msg
3934 
3935 ! *************************************************************************
3936 
3937   write(msg,'(a,1x,a,1x,a,1x,a)') ch10,"=== subprogram :",func_name,"==="
3938   call wrtout(std_out,msg,mode)
3939 
3940 end subroutine vdw_df_write_func

ABINIT/m_xc_vdw/xc_vdw_aggregate [ Functions ]

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

SOURCE

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

ABINIT/m_xc_vdw/xc_vdw_done [ Functions ]

[ 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

SOURCE

925 subroutine xc_vdw_done(vdw_params)
926 
927 !Arguments ------------------------------------
928   type(xc_vdw_type), intent(in) :: vdw_params
929 
930 !Local variables ------------------------------
931   character(len=512) :: msg
932   integer :: i
933 
934 ! *************************************************************************
935 
936   DBG_ENTER("COLL")
937 
938 !--------my debug----------------------------
939   !call libxc_functionals_end(xc_functionals=vdw_funcs)
940 !--------end my debug------------------------
941   ABI_FREE(dmesh)
942   ABI_FREE(qmesh)
943   ABI_FREE(qpoly_basis)
944   ABI_FREE(phi)
945   ABI_FREE(phi_bicubic)
946   ABI_FREE(phi_u_bicubic)
947   ABI_FREE(phir)
948   ABI_FREE(phir_u)
949   ABI_FREE(d2phidr2)
950   ABI_FREE(phig)
951   ABI_FREE(d2phidg2)
952 
953   DBG_EXIT("COLL")
954 
955 end subroutine xc_vdw_done

ABINIT/m_xc_vdw/xc_vdw_energy [ Functions ]

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

SOURCE

708 subroutine xc_vdw_energy(nspden,rho,grho,ex_lda,ec_lda,vx_lda,vc_lda, &
709 & theta,eps,dexc,dexcg,ztmp)
710 
711 !Arguments ------------------------------------
712   integer,intent(in) :: nspden
713   real(dp),intent(in) :: ex_lda,ec_lda,vx_lda,vc_lda
714   real(dp),intent(in) :: rho(nspden),grho(nspden,3)
715   real(dp),intent(inout) :: theta(my_vdw_params%nqpts,nspden,5)
716   real(dp),intent(in),optional :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts)
717   real(dp),intent(out),optional :: eps,dexc(nspden),dexcg(3,nspden)
718 
719 !Local variables ------------------------------
720   character(len=512) :: msg
721   logical :: calc_corrections
722   integer :: iq0,iq1,iq2,ir,is,ix,nqpts,nrpts,ns
723   real(dp) :: decdrho,dexdrho,dvcdrho,dvxdrho,ec,ex,vc,vx
724   real(dp) :: dr,kappa,rho_tmp,grho_tmp(3),grho2_tmp,rtol,qcut,zab
725   real(dp) :: dq(3),ptmp(2),q0(5),qtmp(2)
726 !  real(dp) :: ztmp(my_vdw_params%nqpts,my_vdw_params%nqpts)
727   real(dp) :: qpoly(my_vdw_params%nqpts,3)
728 
729 ! *************************************************************************
730 
731   ! Init
732   nqpts = my_vdw_params%nqpts
733 !  nrpts = my_vdw_params%nrpts
734   ns = my_vdw_params%nsmooth
735   qcut = my_vdw_params%qcut
736   rtol = my_vdw_params%tolerance
737   zab = my_vdw_params%zab
738 
739   calc_corrections = .false.
740   if ( present(eps) .and. present(dexc) .and. present(dexcg) ) then
741     calc_corrections = .true.
742   end if
743 
744   ! Sum density over spin
745   !rho_tmp = sum(rho(1:nspden))
746   ! Get total density
747   if (nspden==1) then
748     rho_tmp = rho(nspden)
749   else
750     rho_tmp = rho(1)
751   end if
752 
753   kappa = (three * pi**2 * rho_tmp)**third
754   forall(ix=1:3) grho_tmp(ix) = sum(grho(1:nspden,ix))
755   grho2_tmp = sum(grho_tmp**2)
756 
757   ! Calculate local wavevector q0 of eqs. 11-12 from DRSLL04
758   !   * q0(1)   = q0
759   !   * q0(2)   = dq0 / drho
760   !   * q0(3:5) = dq0 / dgrho
761   ! Notes:
762   !   * treating rho->0 separately (divide by 0)
763   !   * treating ex_lda->0 separately (divide by 0)
764   q0(:) = zero
765 !  if ( rho_tmp < rtol ) then
766   if ( rho_tmp < rtol ) then ! .or. (ex_lda < rtol) ) then
767     q0(1) = qcut
768     q0(2:5) = zero
769   else
770     q0(1) = kappa * (one + ec_lda/ex_lda - zab/nine * grho2_tmp / &
771 &     (two * kappa * rho_tmp)**2)
772 
773     if ( calc_corrections ) then
774      q0(2) = ( (vc_lda - ec_lda) / (rho_tmp * ex_lda) - ec_lda * &
775 &      (vx_lda - ex_lda) / &
776 &      (rho_tmp * ex_lda**2) + two * zab / nine * grho2_tmp / &
777 &      (two * kappa * rho_tmp)**3 * &
778 &      eight * kappa / three ) * kappa + q0(1) / (three * rho_tmp)
779      q0(3:5) = -two * kappa * (zab / nine) / (two * kappa * rho_tmp)**2 * &
780 &      grho_tmp(1:3)
781     end if
782     ! Smoothen q0 near qcut exponentially  eq. 5 RS09. (Saturation)
783     !   * q0s(1) = q0c * (1 - exp(-Sum_i=1,ns 1/i (q0/q0c)**i)) 
784     !   * q0s(2) = dq0s / dq |q=q0
785     call vdw_df_saturation(q0(1),qcut,qtmp)
786     q0(1) = qtmp(1)
787     q0(2:5) = q0(2:5) * qtmp(2)
788   end if
789 
790   !if ( q0(1) > qcut ) then
791   !  q0(1) = qcut
792   !  q0(2:5)= zero
793   !end if
794 
795 !!qtmp(1) = (q0(1) / qcut) / ns
796 !!qtmp(2) = one / qcut
797 !!do is=ns-1,1,-1
798 !!  qtmp(1) = (qtmp(1) + one / is) * q0(1) / qcut
799 !!  qtmp(2) = (qtmp(2) * q0(1) + one) / qcut
800 !!end do
801 !!qtmp(2) = qtmp(2) * qcut * exp(-qtmp(1))
802 
803 !!q0(1) = qcut * (one - exp(-qtmp(1)))
804 !!q0(2:5) = q0(2:5) * qtmp(2)
805 
806   ! Calculate polynomial coefficients for cubic-spline interpolation at q0
807   !   * qpoly(:,1) = coefficients
808   !   * qpoly(:,2) = first derivatives
809   qpoly(:,:) = zero
810   iq0 = vdw_df_indexof(qmesh(:),nqpts,q0(1))
811   dq(1) = qmesh(iq0+1) - qmesh(iq0)
812   dq(2) = (qmesh(iq0+1) - q0(1)) / dq(1)
813   dq(3) = (q0(1) - qmesh(iq0)) / dq(1)
814   do iq1=1,nqpts
815    qpoly(iq1,1) = dq(2) * qpoly_basis(iq0,iq1,1) + dq(3) * &
816 &    qpoly_basis(iq0+1,iq1,1) + &
817 &    ((dq(2)**3 - dq(2)) * qpoly_basis(iq0,iq1,2) + &
818 &    (dq(3)**3 - dq(3)) * qpoly_basis(iq0+1,iq1,2)) * dq(1)**2 / six
819 
820    if ( calc_corrections ) then
821     qpoly(iq1,2) = -(qpoly_basis(iq0,iq1,1) - &
822 &     qpoly_basis(iq0+1,iq1,1)) / dq(1) - &
823 &     ((three * dq(2)**2 - one) * qpoly_basis(iq0,iq1,2) - &
824 &     (three * dq(3)**2 - one) * qpoly_basis(iq0+1,iq1,2)) * dq(1) / six
825    end if
826   end do
827 
828   ! Pre-compute integrand for vdW energy correction
829   ! by cubic-spline interpolation and store it in
830   ! qpoly(:,3). This is done twice: one for the
831   ! unsoftened kernel and other one for the
832   ! softened kernel
833   !dr = my_vdw_params%rcut / nrpts
834   !ztmp(:,:) = zero
835 
836 ! if ( calc_corrections ) then
837 !  do ir=1,nrpts
838 !    do iq1=1,nqpts
839 !      do iq2=1,iq1
840 !        ztmp(iq1,iq2) = ztmp(iq1,iq2) - two * pi * phir(ir,iq1,iq2) * dr
841 !      end do
842 !    end do
843 !  end do
844 ! else
845 !  do ir=1,nrpts
846 !    do iq1=1,nqpts
847 !      do iq2=1,iq1
848 !        ztmp(iq1,iq2) = ztmp(iq1,iq2) + two * pi * phir_u(ir,iq1,iq2) * dr
849 !      end do
850 !    end do
851 !  end do
852 ! end if ! calc_corrections
853 
854 ! do iq1=1,nqpts
855 !   do iq2=1,iq1-1
856 !     ztmp(iq2,iq1) = ztmp(iq1,iq2)
857 !   end do
858 ! end do
859 ! qpoly(:,3) = matmul(qpoly(:,1),ztmp(:,:))
860 
861   ! Calculate theta and its derivatives (see RS09)
862   !   * theta(:,:,1)   = theta
863   !   * theta(:,:,2)   = dtheta / drho
864   !   * theta(:,:,3:5) = dtheta / dgrho
865   ! Note: trick to go from Abinit to LibXC conventions
866 
867   do is=1,nspden
868     !theta(:,is,1) = qpoly(:,1) * (2 - nspden) * rho(is)
869     theta(:,is,1) = qpoly(:,1) * rho(1)  !CHECK THIS!
870   end do
871 !-------my debug--------------------------------------------------
872 !  write(59,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts) 
873 !-------end my debug----------------------------------------------
874   ! Calculate theta derivatives 
875   if ( calc_corrections ) then
876     do is=1,nspden
877 !      theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(is) 
878       theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(1)
879 !      theta(:,is,2) = qpoly(:,1) + qpoly(:,2) * q0(2) * rho(is) 
880       do ix=3,5
881 !        theta(:,is,ix) = qpoly(:,2) * q0(ix) * (2 - nspden) * rho(is)
882         theta(:,is,ix) = qpoly(:,2) * q0(ix) * rho(1)
883       end do
884     end do
885   end if
886 
887   ! Calculate energy corrections
888 
889   if ( calc_corrections ) then
890 !-------my debug----------------------------------------------
891 !  write(58,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts) 
892 !!     write(59,*) rho_tmp, grho2_tmp, q0(1),(qpoly(iq1,1),iq1=1,nqpts)
893 !!     write(60,*) rho_tmp, grho2_tmp, (theta(iq1,1,1),iq1=1,nqpts)
894 !!     write(61,*) rho_tmp, grho2_tmp, (theta(iq1,1,2),iq1=1,nqpts)
895 !-------end my debug-----------------------------------------
896      qpoly(:,3) = matmul(qpoly(:,1),ztmp(:,:))
897      eps = zero
898      ptmp(1) = sum(qpoly(:,3) * qpoly(:,1))
899      eps = ptmp(1) * rho_tmp
900      dexc(:) = zero
901      dexcg(:,:) = zero
902      ptmp(2) = two * sum(qpoly(:,2) * qpoly(:,3))
903      dexc(:) = two * ptmp(1) * rho_tmp + ptmp(2) * q0(2) * rho_tmp**2
904     do ix=3,5
905       dexcg(ix-2,:) = ptmp(2) * q0(ix) * rho_tmp**2
906     end do
907   end if ! calc_corrections
908 
909 end subroutine xc_vdw_energy

ABINIT/m_xc_vdw/xc_vdw_get_params [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_get_params

FUNCTION

  Exports internal vdW-DF parameters.

OUTPUT

  vdw_params= van der Waals parameters

SOURCE

 970 subroutine xc_vdw_get_params(vdw_params)
 971 
 972 !Arguments ------------------------------------
 973   type(xc_vdw_type), intent(out)  :: vdw_params
 974 
 975 !Local variables ------------------------------
 976   character(len=512) :: msg
 977 
 978 ! *************************************************************************
 979 
 980   vdw_params%functional = my_vdw_params%functional
 981   vdw_params%zab        = my_vdw_params%zab
 982   vdw_params%ndpts      = my_vdw_params%ndpts
 983   vdw_params%dcut       = my_vdw_params%dcut
 984   vdw_params%dratio     = my_vdw_params%dratio
 985   vdw_params%dsoft      = my_vdw_params%dsoft
 986   vdw_params%phisoft    = my_vdw_params%phisoft
 987   vdw_params%nqpts      = my_vdw_params%nqpts
 988   vdw_params%qcut       = my_vdw_params%qcut
 989   vdw_params%qratio     = my_vdw_params%qratio
 990   vdw_params%nrpts      = my_vdw_params%nrpts
 991   vdw_params%rcut       = my_vdw_params%rcut
 992   vdw_params%rsoft      = my_vdw_params%rsoft
 993   vdw_params%ngpts      = my_vdw_params%ngpts
 994   vdw_params%gcut       = my_vdw_params%gcut
 995   vdw_params%acutmin    = my_vdw_params%acutmin
 996   vdw_params%aratio     = my_vdw_params%aratio
 997   vdw_params%damax      = my_vdw_params%damax
 998   vdw_params%damin      = my_vdw_params%damin
 999   vdw_params%nsmooth    = my_vdw_params%nsmooth
1000   vdw_params%tolerance  = my_vdw_params%tolerance
1001   vdw_params%tweaks     = my_vdw_params%tweaks
1002 
1003 end subroutine xc_vdw_get_params

ABINIT/m_xc_vdw/xc_vdw_init [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_init

FUNCTION

  Calculates the van der Waals kernel.

INPUTS

  vdw_params= parameters for the van der Waals calculations

SOURCE

1018 subroutine xc_vdw_init(vdw_params)
1019 
1020 !Arguments ------------------------------------
1021   type(xc_vdw_type), intent(in)  :: vdw_params
1022 
1023 !Local variables ------------------------------
1024   character(len=*),parameter :: dmesh_file = "vdw_df_dmesh.pts"
1025   character(len=*),parameter :: qmesh_file = "vdw_df_qmesh.pts"
1026 
1027   character(len=512) :: msg
1028   integer :: id1,id2,iq1,iq2,ir,ir0,ixc,jd1,jd2,ndpts,ngpts,nqpts,nrpts,ntest,sofswt
1029   integer :: ii,jj ! DEBUG
1030   real(dp) :: d1,d2,dcut,dd,delta_test,dsoft,phisoft,acutmin,aratio,damax
1031   real(dp) :: phimm,phimp,phipm,phipp,phipn,phi_tmp,ptmp(2)
1032   real(dp),allocatable :: deltd(:),kern_spline_der(:,:,:)
1033   real(dp),allocatable :: btmp(:,:),utmp(:),xde(:)
1034 
1035 ! *************************************************************************
1036 
1037   DBG_ENTER("COLL")
1038 
1039   ! Global init
1040   call xc_vdw_set_params(vdw_params)
1041   call vdw_df_set_tweaks(my_vdw_params%tweaks,my_vdw_tweaks)
1042 
1043   ! Local init
1044   ndpts   = my_vdw_params%ndpts
1045   nqpts   = my_vdw_params%nqpts
1046   nrpts   = my_vdw_params%nrpts
1047   acutmin = my_vdw_params%acutmin
1048   aratio  = my_vdw_params%aratio
1049   damax   = my_vdw_params%damax
1050   dsoft   = my_vdw_params%dsoft
1051   phisoft = my_vdw_params%phisoft
1052 
1053   ! Create d-mesh
1054   ! Note: Not avoid zero and make sure the last point is exactly dcut
1055   ABI_MALLOC(dmesh,(ndpts))
1056 #if defined DEBUG_VERBOSE
1057   write(msg,'(1x,a)') "[vdW-DF Build] Generating D-mesh"
1058   call wrtout(std_out,msg,'COLL')
1059 #endif
1060   call vdw_df_create_mesh(dmesh,ndpts,my_vdw_tweaks%dmesh_type, &
1061 &   my_vdw_params%dcut,mesh_ratio=my_vdw_params%dratio, &
1062 &   mesh_tolerance=my_vdw_params%tolerance,mesh_file=dmesh_file, &
1063 &   avoid_zero=.false.,exact_max=.true.)
1064 !DEBUG
1065 !  write(*,*) '-----D-mesh:-------'
1066 !  do id1 = 1,ndpts
1067 !    write(*,*) id1, dmesh(id1)
1068 !  end do
1069 !END DEBUG
1070   ! Create q-mesh
1071   ABI_MALLOC(qmesh,(nqpts))
1072 #if defined DEBUG_VERBOSE
1073   write(msg,'(1x,a)') "[vdW-DF Build] Generating Q-mesh"
1074   call wrtout(std_out,msg,'COLL')
1075 #endif
1076   call vdw_df_create_mesh(qmesh,nqpts,my_vdw_tweaks%qmesh_type, &
1077 &   my_vdw_params%qcut,mesh_ratio=my_vdw_params%qratio, &
1078 &   mesh_tolerance=my_vdw_params%tolerance, mesh_file=qmesh_file)
1079 
1080   ! Build polynomial basis for cubic-spline interpolation
1081   !   * qpoly_basis(:,1) = polynomial basis
1082   !   * qpoly_basis(:,2) = second derivative
1083   ABI_MALLOC(qpoly_basis,(nqpts,nqpts,2))
1084 !  ABI_MALLOC(btmp,(nqpts,nqpts))
1085   ABI_MALLOC(utmp,(nqpts))
1086 #if defined DEBUG_VERBOSE
1087   write(msg,'(1x,a)') "[vdW-DF Build] Generating polynomial basis"
1088   call wrtout(std_out,msg,'COLL')
1089 #endif
1090 !  btmp(:,:) = zero
1091 !  forall(iq1=1:nqpts) btmp(iq1,iq1) = one
1092   utmp(:) = zero
1093   qpoly_basis(:,:,:) = zero
1094   do iq1=1,nqpts
1095     qpoly_basis(iq1,iq1,1) = one
1096     qpoly_basis(1,iq1,2) = zero
1097     utmp(1) = zero
1098     do iq2=2,nqpts-1
1099       ptmp(1) = (qmesh(iq2)-qmesh(iq2-1)) / (qmesh(iq2+1)-qmesh(iq2-1))
1100       ptmp(2) = ptmp(1) * qpoly_basis(iq2-1,iq1,2) + two
1101       qpoly_basis(iq2,iq1,2) = (ptmp(1) - one) / ptmp(2)
1102       utmp(iq2) = ( 6*( (qpoly_basis(iq2+1,iq1,1)-qpoly_basis(iq2,iq1,1))&
1103 &     /(qmesh(iq2+1)-qmesh(iq2)) - (qpoly_basis(iq2,iq1,1)-&
1104 &     qpoly_basis(iq2-1,iq1,1))/(qmesh(iq2)-qmesh(iq2-1)) ) /&
1105 &      (qmesh(iq2+1)-qmesh(iq2-1)) - ptmp(1)*utmp(iq2-1) ) / ptmp(2)
1106 
1107 !      ptmp(1) = (qmesh(iq2) - qmesh(iq2-1)) / &
1108 !&     (qmesh(iq2+1) - qmesh(iq2-1))
1109 !      ptmp(2) = ptmp(1) * qpoly_basis(iq1,iq2-1,2) + two
1110 !      qpoly_basis(iq1,iq2,2) = (ptmp(1) - one) / ptmp(2)
1111 !      utmp(iq2) = ( six * ((btmp(iq1,iq2+1) - btmp(iq1,iq2)) / &
1112 !&     (qmesh(iq2+1) - qmesh(iq2)) - (btmp(iq1,iq2) - &
1113 !&     btmp(iq1,iq2-1)) / (qmesh(iq2) - qmesh(iq2-1))) / &
1114 !&     (qmesh(iq2+1) - qmesh(iq2-1)) - ptmp(1) * utmp(iq2-1)) / &
1115 !&     ptmp(2)
1116     end do
1117 !    utmp(nqpts) = zero
1118     qpoly_basis(nqpts,iq1,2) = zero
1119     do iq2=nqpts-1,1,-1
1120       qpoly_basis(iq2,iq1,2) = qpoly_basis(iq2,iq1,2) * &
1121 &       qpoly_basis(iq2+1,iq1,2) + utmp(iq2)
1122     end do
1123   end do
1124 !  ABI_DEALLOCATE(btmp)
1125   ABI_FREE(utmp)
1126 
1127 !DEBUG-----------------------
1128 ! open(37,file='qmesh-pofq-d2pdq2-check.dat',status='replace')
1129 ! do iq1=1,nqpts
1130 !   do iq2=1,nqpts
1131 !    write(37,'(4e15.6)') qmesh(iq1), qmesh(iq2), qpoly_basis(iq1,iq2,1) &
1132 !    , qpoly_basis(iq1,iq2,2)
1133 !   end do
1134 ! end do
1135 !END DEBUG------------------
1136 
1137   ! Create kernel and its derivatives
1138   ! Note: using 4 neighbours for derivatives
1139   ABI_MALLOC(phi,(ndpts,ndpts,4))
1140 #if defined DEBUG_VERBOSE
1141   write(msg,'(1x,a)') "[vdW-DF Build] Building kernel and its derivatives"
1142   call wrtout(std_out,msg,'COLL')
1143 #endif
1144 
1145   if ( phisoft < zero ) then
1146     phisoft = 0.4_dp  !three * half * exp(-dsoft * two / sqrt(two))
1147     my_vdw_params%phisoft = phisoft
1148   end if
1149 
1150   ! Building of the softened kernel
1151 
1152   sofswt = 1
1153   do id1=1,ndpts
1154     d1 = dmesh(id1)
1155     phi(id1,id1,1) = vdw_df_kernel(d1,d1,dsoft,phisoft,acutmin,aratio,damax)
1156 
1157     ! Delta(d) should be at least 10^-6
1158     ! WARNING: changed tol3 to tol6 and changed max to min
1159     dd = max(0.01_dp*d1,tol3)   !min(my_vdw_params%dratio*d1,tol6)
1160 
1161       phimm = vdw_df_kernel(d1-dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax)
1162       phipm = vdw_df_kernel(d1+dd,d1-dd,dsoft,phisoft,acutmin,aratio,damax)
1163       phipp = vdw_df_kernel(d1+dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax)
1164       phimp = vdw_df_kernel(d1-dd,d1+dd,dsoft,phisoft,acutmin,aratio,damax)
1165 
1166       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1167       phi(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1168       phi(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1169       phi(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1170 
1171     do id2=1,id1-1
1172       d2 = dmesh(id2)
1173       phi(id1,id2,1) = vdw_df_kernel(d1,d2,dsoft,phisoft,acutmin,aratio,damax)
1174       phi(id2,id1,1) = phi(id1,id2,1)
1175 
1176       phimm = vdw_df_kernel(d1-dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax)
1177       phipm = vdw_df_kernel(d1+dd,d2-dd,dsoft,phisoft,acutmin,aratio,damax)
1178       phipp = vdw_df_kernel(d1+dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax)
1179       phimp = vdw_df_kernel(d1-dd,d2+dd,dsoft,phisoft,acutmin,aratio,damax)
1180 
1181       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)
1182       phi(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd)
1183       phi(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd)
1184       phi(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)
1185 
1186       phi(id2,id1,2) = phi(id1,id2,3)
1187       phi(id2,id1,3) = phi(id1,id2,2)
1188       phi(id2,id1,4) = phi(id1,id2,4)
1189     end do
1190 
1191     write(msg,'(1x,a,1x,i3,"% complete")') &
1192 &     '[vdW-DF Build]',int(id1*100.0/ndpts)
1193     call wrtout(std_out,msg,'COLL')
1194 #if defined DEBUG_VERBOSE
1195     call flush_unit(std_out)
1196 #endif
1197   end do
1198   write(msg,'(a)') " "
1199   call wrtout(std_out,msg,'COLL')
1200 
1201   ! Set boundary conditions
1202   ! Note: will have to check that the borders are smooth enough
1203   ! These boundary conditions produce large discontinuities, now testing its removal
1204   phi(ndpts,:,:) = zero
1205   phi(:,ndpts,:) = zero
1206   phi(1,:,2) = zero
1207   phi(:,1,3) = zero
1208   phi(1,:,4) = zero
1209   phi(:,1,4) = zero
1210   
1211 !DEBUG------
1212 ! write(*,*) 'd1 -- d2 -- phi -- dphi/dd1 -- dphi/dd2 -- d2phi/dd1dd2'
1213 ! do id1 = 1,ndpts
1214 !   do id2 = 1,ndpts
1215 !     d1 = dmesh(id1)
1216 !     d2 = dmesh(id2)
1217 !     write(*,'(2f10.6,4e15.6)') d1, d2, phi(id1,id2,1), &
1218 !       phi(id1,id2,2), phi(id1,id2,3), phi(id1,id2,4)
1219 !   end do
1220 ! end do
1221 !END DEBUG-------
1222   ! Ar2 results show better convergence of E_vdW than when boundary conditions were used.
1223   ! Kernel plot show that there are not dicontinuities as well. 
1224 
1225 #if defined DEBUG_VERBOSE
1226   write(msg,'(1x,a)') "[vdW-DF Build] Building filtered kernel"
1227   call wrtout(std_out,msg,'COLL')
1228 #endif
1229 
1230   ! Calculate coefficients for bicubic interpolation
1231   ABI_MALLOC(phi_bicubic,(4,4,ndpts,ndpts))
1232   call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi(:,:,1),phi(:,:,2),&
1233 &   phi(:,:,3),phi(:,:,4),phi_bicubic)
1234 
1235 !DEBUG------------------------------
1236 ! open (unit=41, file='phi-bicubic-check.dat',status='replace')
1237 ! do id1=1,ndpts
1238 !   do id2=1,ndpts
1239 !     write(41,*) id1, id2, ((phi_bicubic(ii,jj,id1,id2),ii=1,4),jj=1,4)
1240 !   end do
1241 ! end do
1242 ! close(unit=41)
1243 !END DEBUG------------------------
1244   ! Build filtered kernel
1245   ABI_MALLOC(phir,(0:nrpts,nqpts,nqpts))
1246   ABI_MALLOC(d2phidr2,(0:nrpts,nqpts,nqpts))
1247   ABI_MALLOC(phig,(0:nrpts,nqpts,nqpts))
1248   ABI_MALLOC(d2phidg2,(0:nrpts,nqpts,nqpts))
1249   call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt)
1250   my_vdw_params%ngpts = ngpts
1251 
1252   ! Find closest indices in dmesh
1253   ! FIXME: something is missing or should be removed here
1254 ! ABI_MALLOC(kern_spline_der,(ndpts,ndpts,3))
1255 ! ABI_MALLOC(deltd,(2))
1256 ! ABI_MALLOC(xde,(2))
1257 
1258 ! kern_spline_der(:,:,:) = zero
1259 
1260 ! do jd1=1,ndpts
1261 !   do jd2=1,ndpts
1262 
1263 !     d1 = dmesh(jd1)
1264 !     d2 = dmesh(jd2)
1265 
1266 !     if ( (d1 < dcut) .or. (d2 < dcut) ) then
1267 !       id1 = vdw_df_indexof(dmesh,ndpts,d1)
1268 !       id2 = vdw_df_indexof(dmesh,ndpts,d2)
1269 
1270 !       deltd(1) = dmesh(id1+1) - dmesh(id1)
1271 !       deltd(2) = dmesh(id2+1) - dmesh(id2)
1272 
1273 !       xde(1) = (d1 - dmesh(id1)) / deltd(1)
1274 !       xde(2) = (d2 - dmesh(id2)) / deltd(2)
1275 
1276 !       kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1)
1277 !       kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2)
1278 !       kern_spline_der(jd1,jd2,3) = &
1279 !         phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) )
1280 !     end if
1281 
1282 !   end do
1283 ! end do
1284 
1285 ! ABI_FREE(kern_spline_der)
1286 ! ABI_FREE(deltd)
1287 ! ABI_FREE(xde)
1288 
1289   ! Building of the unsoftened kernel
1290 
1291   ! Create unsoftened kernel and its derivatives
1292   ! Note: using 4 neighbours for derivatives
1293   ABI_MALLOC(phi_u,(ndpts,ndpts,4))
1294 #if defined DEBUG_VERBOSE
1295   write(msg,'(1x,a)') "[vdW-DF Build] Building unsoftened kernel and its derivatives"
1296   call wrtout(std_out,msg,'COLL')
1297 #endif
1298 !------my debug--------------------------------------------------------
1299   phi_u(:,:,:) = zero
1300 
1301 !   do id1=1,ndpts                                                                          
1302 !     d1 = dmesh(id1)                                                                       
1303 !     phi_u(id1,id1,1) = vdw_df_kernel_value(d1,d1,acutmin,aratio,damax)                    
1304 !                                                                                           
1305 !     ! Delta(d) should be at least 10^-6                                                   
1306 !     ! WARNING: changed tol3 to tol6 and changed max to min                                
1307 !     dd = max(0.01_dp*d1,tol3) !min(my_vdw_params%dratio*d1,tol6)                                                
1308 !                                                                                           
1309 !       phimm = vdw_df_kernel_value(d1-dd,d1-dd,acutmin,aratio,damax)                       
1310 !       phipm = vdw_df_kernel_value(d1+dd,d1-dd,acutmin,aratio,damax)                       
1311 !       phipp = vdw_df_kernel_value(d1+dd,d1+dd,acutmin,aratio,damax)                       
1312 !       phimp = vdw_df_kernel_value(d1-dd,d1+dd,acutmin,aratio,damax)                       
1313 !                                                                                           
1314 !       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)            
1315 !       phi_u(id1,id1,2) = (phipp + phipm - phimp - phimm) / (four * dd)                    
1316 !       phi_u(id1,id1,3) = (phipp - phipm + phimp - phimm) / (four * dd)                    
1317 !       phi_u(id1,id1,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)                
1318 !                                                                                           
1319 !     do id2=1,id1-1                                                                        
1320 !       d2 = dmesh(id2)                                                                     
1321 !       phi_u(id1,id2,1) = vdw_df_kernel_value(d1,d2,acutmin,aratio,damax)                  
1322 !       phi_u(id2,id1,1) = phi_u(id1,id2,1)                                                 
1323 !                                                                                           
1324 !       phimm = vdw_df_kernel_value(d1-dd,d2-dd,acutmin,aratio,damax)                       
1325 !       phipm = vdw_df_kernel_value(d1+dd,d2-dd,acutmin,aratio,damax)                       
1326 !       phipp = vdw_df_kernel_value(d1+dd,d2+dd,acutmin,aratio,damax)                       
1327 !       phimp = vdw_df_kernel_value(d1-dd,d2+dd,acutmin,aratio,damax)                       
1328 !                                                                                           
1329 !       ! Using five point stencil formula for crossed derivative phi(id1,id2,4)            
1330 !       phi_u(id1,id2,2) = (phipp + phipm - phimp - phimm) / (four * dd)                    
1331 !       phi_u(id1,id2,3) = (phipp - phipm + phimp - phimm) / (four * dd)                    
1332 !       phi_u(id1,id2,4) = (phipp - phipm - phimp + phimm) / ((two * dd)**2)                
1333 !                                                                                           
1334 !       phi_u(id2,id1,2) = phi_u(id1,id2,2)                                                 
1335 !       phi_u(id2,id1,3) = phi_u(id1,id2,3)                                                 
1336 !       phi_u(id2,id1,4) = phi_u(id1,id2,4)                                                 
1337 !     end do                                                                                
1338 !                                                                                           
1339 !     write(msg,'(1x,a,1x,i3,"% complete")') &                                              
1340 ! &     '[vdW-DF-Uns Build]',int(id1*100.0/ndpts)                                           
1341 !     call wrtout(std_out,msg,'COLL')                                                       
1342 ! #if defined DEBUG_VERBOSE                                                                 
1343 !     call flush_unit(std_out)                                                              
1344 ! #endif                                                                                    
1345 !   end do                                                                                  
1346 !   write(msg,'(a)') " "                                                                    
1347 !   call wrtout(std_out,msg,'COLL')                                                         
1348 !                                                                                           
1349 !   ! Set boundary conditions                                                               
1350 !   ! Note: will have to check that the borders are smooth enough                           
1351 !   ! These boundary conditions produce large discontinuities, now testing its removal      
1352 !   phi_u(ndpts,:,:) = zero
1353 !   phi_u(:,ndpts,:) = zero
1354 !   phi_u(1,:,2) = zero
1355 !   phi_u(:,1,3) = zero
1356 !   phi_u(1,:,4) = zero
1357 !   phi_u(:,1,4) = zero
1358 ! 
1359 ! 
1360 !   ! Ar2 results show better convergence of E_vdW than when boundary conditions were used. 
1361 !   ! Kernel plot show that there are not dicontinuities as well.                           
1362 !-------end my debug-------------------------------------------------------
1363 #if defined DEBUG_VERBOSE
1364   write(msg,'(1x,a)') "[vdW-DF Build] Building filtered unsoftened kernel"
1365   call wrtout(std_out,msg,'COLL')
1366 #endif
1367 
1368   ! Calculate coefficients for bicubic interpolation
1369   ABI_MALLOC(phi_u_bicubic,(4,4,ndpts,ndpts))
1370   call spline_bicubic(ndpts,ndpts,dmesh,dmesh,phi_u(:,:,1),phi_u(:,:,2),&
1371 &   phi_u(:,:,3),phi_u(:,:,4),phi_u_bicubic)
1372 
1373   ! Build filtered kernel
1374   ABI_MALLOC(phir_u,(0:nrpts,nqpts,nqpts))
1375   ! For the local correction we just need the kernel not its
1376   ! derivatives.
1377   sofswt = 0
1378   !ABI_MALLOC(d2phidr2,(nrpts,nqpts,nqpts))
1379   !ABI_MALLOC(phig,(nrpts,nqpts,nqpts))
1380   !ABI_MALLOC(d2phidg2,(nrpts,nqpts,nqpts))
1381   call vdw_df_filter(nqpts,nrpts,my_vdw_params%rcut,my_vdw_params%gcut,ngpts,sofswt)
1382   !  my_vdw_params%ngpts = ngpts already defined above
1383 
1384   ! The following is not needed for the local correction:
1385   ! Find closest indices in dmesh
1386   ! FIXME: something is missing or should be removed here
1387   ! ABI_MALLOC(kern_spline_der,(ndpts,ndpts,3))
1388   ! ABI_MALLOC(deltd,(2))
1389   ! ABI_MALLOC(xde,(2))
1390 
1391   ! kern_spline_der(:,:,:) = zero
1392 
1393   ! do jd1=1,ndpts
1394   !   do jd2=1,ndpts
1395 
1396   !     d1 = dmesh(jd1)
1397   !     d2 = dmesh(jd2)
1398 
1399   !     if ( (d1 < dcut) .or. (d2 < dcut) ) then
1400   !       id1 = vdw_df_indexof(dmesh,ndpts,d1)
1401   !       id2 = vdw_df_indexof(dmesh,ndpts,d2)
1402 
1403   !       deltd(1) = dmesh(id1+1) - dmesh(id1)
1404   !       deltd(2) = dmesh(id2+1) - dmesh(id2)
1405 
1406   !       xde(1) = (d1 - dmesh(id1)) / deltd(1)
1407   !       xde(2) = (d2 - dmesh(id2)) / deltd(2)
1408 
1409   !       kern_spline_der(jd1,jd2,1) = phi_bicubic(2,1,id1,id2) / deltd(1)
1410   !       kern_spline_der(jd1,jd2,2) = phi_bicubic(1,2,id1,id2) / deltd(2)
1411   !       kern_spline_der(jd1,jd2,3) = &
1412   !&         phi_bicubic(2,2,id1,id2) / ( deltd(1)*deltd(2) )
1413   !     end if
1414 
1415   !   end do
1416   ! end do
1417 
1418   ! ABI_FREE(kern_spline_der)
1419   ! ABI_FREE(deltd)
1420   ! ABI_FREE(xde)
1421 
1422   call vdw_df_internal_checks(1)
1423 
1424   vdw_switch = .true.
1425 
1426   DBG_EXIT("COLL")
1427 
1428 end subroutine xc_vdw_init

ABINIT/m_xc_vdw/xc_vdw_libxc_init [ Functions ]

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

SOURCE

1451 subroutine xc_vdw_libxc_init(ixc_vdw)
1452 
1453 !Arguments ------------------------------------
1454   integer,intent(in) :: ixc_vdw
1455 
1456 !Local variables ------------------------------
1457   character(len=*),parameter :: c_vdw1 = "LDA_C_PW"
1458   character(len=*),parameter :: x_vdw1 = "LDA_X" !"GGA_X_PBE_R"
1459   character(len=*),parameter :: x_vdw2 = "GGA_X_PW86"
1460   character(len=*),parameter :: x_vdw3 = "GGA_X_C09X"
1461   integer :: i,ii,ic,ix,ixc
1462   character(len=500) :: msg
1463 
1464 ! *************************************************************************
1465 
1466   DBG_ENTER("COLL")
1467 
1468   ! Select LibXC functionals
1469   select case (ixc_vdw)
1470     case (1)
1471       ix = libxc_functionals_getid(x_vdw1)
1472       ic = libxc_functionals_getid(c_vdw1)
1473     case (2)
1474       ix = libxc_functionals_getid(x_vdw2)
1475       ic = libxc_functionals_getid(c_vdw1)
1476     case (3)
1477       ix = libxc_functionals_getid(x_vdw3)
1478       ic = libxc_functionals_getid(c_vdw1)
1479       ABI_ERROR('[vdW-DF] C09 not available for now')
1480     case default
1481       ABI_ERROR('[vdW-DF] invalid setting of vdw_xc')
1482   end select
1483   if ( (ix == -1) .or. (ic == -1) ) then
1484     ABI_ERROR('[vdW-DF] unable to set LibXC parameters')
1485   end if
1486   ixc = -(ix * 1000 + ic)
1487 
1488   ! Propagate to internal parameters
1489   my_vdw_params%functional = ixc_vdw !explore what happens with this variable
1490 
1491   ! XC functional init
1492   call libxc_functionals_init(ixc,1,xc_functionals=vdw_funcs)
1493 
1494 !------------my debug------------
1495 ! write(*,*) 'vdw_funcs=', vdw_funcs(1)%id, vdw_funcs(2)%id
1496 ! call libxc_functionals_end(xc_functionals=vdw_funcs)
1497 !-----------end my debug-------
1498 
1499   DBG_EXIT("COLL")
1500 
1501 end subroutine xc_vdw_libxc_init

ABINIT/m_xc_vdw/xc_vdw_memcheck [ Functions ]

[ 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

SOURCE

1517 subroutine xc_vdw_memcheck(unt,vp)
1518 
1519 !Arguments ------------------------------------
1520   integer,intent(in) :: unt
1521   type(xc_vdw_type), intent(in), optional :: vp
1522 
1523 !Local variables ------------------------------
1524   character(len=1536) :: msg
1525   integer :: napts,ndpts,nqpts,nrpts,id1,id2
1526   real(dp) :: dmax,dtmp,mem_perm,mem_temp
1527   type(xc_vdw_type) :: my_vp
1528 
1529 ! *************************************************************************
1530 
1531   if ( present(vp) ) then
1532     my_vp = vp
1533   else
1534     my_vp = my_vdw_params
1535   end if
1536 
1537   ndpts = my_vp%ndpts
1538   nqpts = my_vp%nqpts
1539   nrpts = my_vp%nrpts
1540 
1541   if ( .not. allocated(dmesh) ) return
1542 
1543   dmax = zero
1544   do id1=1,ndpts
1545     do id2=1,id1
1546       dtmp = sqrt(dmesh(id1)**2+dmesh(id2)**2)
1547       if ( dtmp > dmax ) dmax = dtmp
1548     end do
1549   end do
1550   napts = nint(max(my_vp%acutmin,my_vp%aratio*dmax))
1551 
1552   mem_perm = ( &
1553 &   4 * ndpts * (1 + 5 * ndpts ) + &
1554 &   2 * nqpts + (2 + 4 * nrpts) * nqpts * nqpts &
1555 &   ) * sizeof(one) / 1048576.0_dp
1556 
1557   mem_temp = ( &
1558 &   7 * napts + &
1559 &   4 * nqpts + 3* nqpts * nqpts + &
1560 &   5* nrpts + 2 &
1561 &   ) * sizeof(one) / 1048576.0_dp
1562 
1563   write(msg,'(a,1x,3(a),12(3x,a,1x,i8,1x,"elts",a), &
1564 &   a,10x,a,1x,f10.3,1x,"Mb",a,a,13(3x,a,1x,i8,1x,"elts",a), &
1565 &   a,10x,a,1x,f10.3,1x,"Mb",a,a,14x,a,1x,f10.3,1x,"Mb",2(a),1x,2(a))') &
1566 &   ch10, &
1567 &   '==== vdW-DF memory estimation ====',ch10,ch10, &
1568 &   'd2phidg2 .........',nrpts*nqpts*nqpts,ch10, &
1569 &   'd2phidr2 .........',nrpts*nqpts*nqpts,ch10, &
1570 &   'dmesh ............',ndpts * 4,ch10, &
1571 &   'phi ..............',ndpts*ndpts*4,ch10, &
1572 &   'phi_u ............',ndpts*ndpts*4,ch10, &
1573 &   'phi_bicubic ......',4*4*ndpts*ndpts,ch10, &
1574 &   'phi_u_bicubic ....',4*4*ndpts*ndpts,ch10, &
1575 &   'phig .............',nrpts*nqpts*nqpts,ch10, &
1576 &   'phir .............',nrpts*nqpts*nqpts,ch10, &
1577 &   'phir_u ...........',nrpts*nqpts*nqpts,ch10, &
1578 &   'qmesh ............',nqpts * 4,ch10, &
1579 &   'qpoly_basis.......',nqpts * nqpts * 2,ch10,ch10, &
1580 &   'Permanent =',mem_perm,ch10,ch10, &
1581 &   'amesh ............',napts,ch10, &
1582 &   'amesh_cos ........',napts,ch10, &
1583 &   'amesh_sin ........',napts,ch10, &
1584 &   'btmp .............',nqpts*nqpts,ch10, &
1585 &   'dphida ...........',napts,ch10, &
1586 &   'dphidb ...........',napts,ch10, &
1587 &   'ftmp .............',2*(2*nrpts+1),ch10, &
1588 &   'nu1 ..............',napts,ch10, &
1589 &   'nu2 ..............',napts,ch10, &
1590 &   'qpoly ............',nqpts*3,ch10, &
1591 &   'utmp(q) ..........',nqpts,ch10, &
1592 &   'utmp(r) ..........',nrpts,ch10, &
1593 &   'wtmp .............',nqpts*nqpts*2,ch10,ch10, &
1594 &   'Temporary =',mem_temp,ch10,ch10, &
1595 &   'TOTAL =',mem_perm+mem_temp,ch10, &
1596 &   ch10, &
1597 &   '==================================',ch10
1598 
1599   call wrtout(unt,msg,'COLL')
1600 
1601 end subroutine xc_vdw_memcheck

ABINIT/m_xc_vdw/xc_vdw_read [ Functions ]

[ 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

SOURCE

1619 subroutine xc_vdw_read(filename)
1620 
1621 !Arguments ------------------------------------
1622   character(len=*), intent(in)  :: filename
1623 
1624 !Local variables ------------------------------
1625   character(len=512) :: msg
1626   character(len=64) :: format_string
1627   integer :: ncid,dimids(4),varids(13)
1628   integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id
1629   integer :: nderivs,npoly,ndpts,ngpts,nqpts,nrpts
1630 
1631 ! *************************************************************************
1632 
1633   DBG_ENTER("COLL")
1634 
1635 #if defined HAVE_NETCDF
1636   write(msg,'("Reading vdW-DF data from",1x,a)') trim(filename)
1637   call wrtout(std_out,msg,'COLL')
1638 
1639   ! Open file for reading
1640   NETCDF_VDWXC_CHECK(nf90_open(trim(filename),NF90_NOWRITE,ncid=ncid))
1641 
1642   ! Get file format and generator
1643   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'file_format',format_string))
1644   write(msg,'(3x,"File format:",1x,a)') trim(format_string)
1645   call wrtout(std_out,msg,'COLL')
1646   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'generator',msg))
1647   write(msg,'(3x,"Generator:",1x,a)') trim(msg)
1648   call wrtout(std_out,msg,'COLL')
1649 
1650   ! Get attributes
1651   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin))
1652   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio))
1653   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax))
1654   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin))
1655   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut))
1656   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio))
1657   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft))
1658   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional))
1659   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut))
1660   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth))
1661   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft))
1662   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut))
1663   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio))
1664   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut))
1665   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft))
1666   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance))
1667   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks))
1668   NETCDF_VDWXC_CHECK(nf90_get_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab))
1669 
1670   ! Get dimensions
1671   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nderivs',nderivs_id))
1672   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'npoly',npoly_id))
1673   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ndpts',ndpts_id))
1674   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'ngpts',ngpts_id))
1675   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nqpts',nqpts_id))
1676   NETCDF_VDWXC_CHECK(nf90_inq_dimid(ncid,'nrpts',nrpts_id))
1677 
1678   ! Get varids
1679   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qmesh',varids(2)))
1680   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'dmesh',varids(3)))
1681   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi',varids(4)))
1682   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u',varids(5)))
1683   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_bicubic',varids(6)))
1684   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phi_u_bicubic',varids(7)))
1685   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir',varids(8)))
1686   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phir_u',varids(9)))
1687   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidr2',varids(10)))
1688   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'phig',varids(11)))
1689   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'d2phidg2',varids(12)))
1690   NETCDF_VDWXC_CHECK(nf90_inq_varid(ncid,'qpoly_basis',varids(13)))
1691 
1692 ! DEBUG
1693 ! write(*,*) ch10,'Inside xc_vdw_read routine',ch10
1694 ! write(*,*) 'nrpts before =',nrpts,ch10
1695 ! write(*,*) 'varids(8)=',varids(8)
1696 ! write(*,*) 'ndpts before =',ndpts,ch10
1697 !  nrpts = nrpts + 1 
1698 ! ENDDEBUG
1699   ! Get dimensions
1700   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nderivs_id,len=nderivs))
1701   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,npoly_id,len=npoly))
1702   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ndpts_id,len=ndpts))
1703   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,ngpts_id,len=ngpts))
1704   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nqpts_id,len=nqpts))
1705   NETCDF_VDWXC_CHECK(nf90_inquire_dimension(ncid,nrpts_id,len=nrpts))
1706 
1707 ! DEBUG
1708 ! write(*,*) 'nrpts after =',nrpts,ch10
1709 ! write(*,*) 'ndpts after =',ndpts,ch10
1710 !  nrpts = nrpts + 1 
1711 ! ENDDEBUG
1712 
1713   ! Get qmesh
1714   ABI_MALLOC(qmesh,(nqpts))
1715   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(2),qmesh))
1716 
1717   ! Get dmesh
1718   ABI_MALLOC(dmesh,(ndpts))
1719   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(3),dmesh))
1720 
1721   ! Get phi
1722   ABI_MALLOC(phi,(ndpts,ndpts,nderivs))
1723   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(4),phi))
1724 
1725   ! Get phi_u
1726   ABI_MALLOC(phi_u,(ndpts,ndpts,nderivs))
1727   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(5),phi_u))
1728 
1729   ! Get phi_bicubic
1730   ABI_MALLOC(phi_bicubic,(nderivs,nderivs,ndpts,ndpts))
1731   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(6),phi_bicubic))
1732 
1733   ! Get phi_u_bicubic
1734   ABI_MALLOC(phi_u_bicubic,(nderivs,nderivs,ndpts,ndpts))
1735   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(7),phi_u_bicubic))
1736 
1737   ! Get phir
1738   ABI_MALLOC(phir,(0:nrpts-1,nqpts,nqpts))
1739   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(8),phir))
1740 
1741   ! Get phir_u
1742   ABI_MALLOC(phir_u,(0:nrpts-1,nqpts,nqpts))
1743   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(9),phir_u))
1744 
1745   ! Get d2phidr2
1746   ABI_MALLOC(d2phidr2,(0:nrpts-1,nqpts,nqpts))
1747   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(10),d2phidr2))
1748 
1749   ! Get phig
1750   ABI_MALLOC(phig,(0:nrpts-1,nqpts,nqpts))
1751   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(11),phig))
1752 
1753   ! Get d2phidg2
1754   ABI_MALLOC(d2phidg2,(0:nrpts-1,nqpts,nqpts))
1755   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(12),d2phidg2))
1756 
1757   ! Get qpoly_basis
1758   ABI_MALLOC(qpoly_basis,(nqpts,nqpts,npoly))
1759   NETCDF_VDWXC_CHECK(nf90_get_var(ncid,varids(13),qpoly_basis))
1760 
1761   ! Close file
1762   NETCDF_VDWXC_CHECK(nf90_close(ncid))
1763 
1764   ! Update my_vdw_params
1765   my_vdw_params%ndpts = ndpts
1766   my_vdw_params%ngpts = ngpts
1767   my_vdw_params%nqpts = nqpts
1768   my_vdw_params%nrpts = nrpts-1
1769 
1770 ! DEBUG 
1771 ! write(*,*) ch10,'At the end of xc_vdw_read',ch10
1772 ! write(*,*) 'nrpts= ',nrpts
1773 ! END DEBUG
1774 
1775 #else
1776   ABI_ERROR('reading vdW-DF variables requires NetCDF')
1777 #endif
1778 
1779   DBG_EXIT("COLL")
1780 
1781 end subroutine xc_vdw_read

ABINIT/m_xc_vdw/xc_vdw_set_functional [ Functions ]

[ 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

SOURCE

1797 subroutine xc_vdw_set_functional(vdw_func,vdw_zab)
1798 
1799 !Arguments ------------------------------------
1800   integer, intent(in)  :: vdw_func
1801   real(dp), optional, intent(in) :: vdw_zab
1802 
1803 !Local variables ------------------------------
1804   character(len=512) :: msg
1805 
1806 ! *************************************************************************
1807 
1808   my_vdw_params%functional = vdw_func
1809   if ( present(vdw_zab) ) then
1810     my_vdw_params%zab = vdw_zab
1811   end if
1812 
1813 end subroutine xc_vdw_set_functional

ABINIT/m_xc_vdw/xc_vdw_set_params [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_set_params

FUNCTION

  Imports external vdW-DF parameters.

 INPUT
  vdw_params= van der Waals parameters

SOURCE

1828 subroutine xc_vdw_set_params(vdw_params)
1829 
1830 !Arguments ------------------------------------
1831   type(xc_vdw_type), intent(in)  :: vdw_params
1832 
1833 !Local variables ------------------------------
1834   character(len=512) :: msg
1835 
1836 ! *************************************************************************
1837 
1838   my_vdw_params%functional = vdw_params%functional
1839   my_vdw_params%zab        = vdw_params%zab
1840   my_vdw_params%ndpts      = vdw_params%ndpts
1841   my_vdw_params%dcut       = vdw_params%dcut
1842   my_vdw_params%dratio     = vdw_params%dratio
1843   my_vdw_params%dsoft      = vdw_params%dsoft
1844   my_vdw_params%phisoft    = vdw_params%phisoft
1845   my_vdw_params%nqpts      = vdw_params%nqpts
1846   my_vdw_params%qcut       = vdw_params%qcut
1847   my_vdw_params%qratio     = vdw_params%qratio
1848   my_vdw_params%nrpts      = vdw_params%nrpts
1849   my_vdw_params%rcut       = vdw_params%rcut
1850   my_vdw_params%rsoft      = vdw_params%rsoft
1851   my_vdw_params%ngpts      = vdw_params%ngpts
1852   my_vdw_params%gcut       = vdw_params%gcut
1853   my_vdw_params%acutmin    = vdw_params%acutmin
1854   my_vdw_params%aratio     = vdw_params%aratio
1855   my_vdw_params%damax      = vdw_params%damax
1856   my_vdw_params%damin      = vdw_params%damin
1857   my_vdw_params%nsmooth    = vdw_params%nsmooth
1858   my_vdw_params%tolerance  = vdw_params%tolerance
1859   my_vdw_params%tweaks     = vdw_params%tweaks
1860 
1861 end subroutine xc_vdw_set_params

ABINIT/m_xc_vdw/xc_vdw_show [ Functions ]

[ 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

SOURCE

1877 subroutine xc_vdw_show(unt,vp)
1878 
1879 !Arguments ------------------------------------
1880   integer,intent(in) :: unt
1881   type(xc_vdw_type), intent(in), optional :: vp
1882 
1883 !Local variables ------------------------------
1884   character(len=1536) :: msg
1885   type(xc_vdw_type) :: my_vp
1886 
1887 ! *************************************************************************
1888 
1889   if ( present(vp) ) then
1890     my_vp = vp
1891   else
1892     my_vp = my_vdw_params
1893   end if
1894 
1895   write(msg,'(a,1x,3(a),3x,a,1x,i9,a,3x,a,1x,f9.3,a, &
1896 &   5(3x,a,1x,i9,a),14(3x,a,1x,f9.3,a),1(3x,a,1x,i9,a),a,1x,a,a)') &
1897 &   ch10, &
1898 &   '==== vdW-DF internal parameters ====',ch10,ch10, &
1899 &   'functional .............',my_vp%functional,ch10, &
1900 &   'zab ....................',my_vp%zab,ch10, &
1901 &   'd-mesh points ..........',my_vp%ndpts,ch10, &
1902 &   'q-mesh points ..........',my_vp%nqpts,ch10, &
1903 &   'r-mesh points ..........',my_vp%nrpts,ch10, &
1904 &   'g-mesh points ..........',my_vp%ngpts,ch10, &
1905 &   'smoothing iterations ...',my_vp%nsmooth,ch10, &
1906 &   'd-mesh cut-off .........',my_vp%dcut,ch10, &
1907 &   'd-mesh ratio ...........',my_vp%dratio,ch10, &
1908 &   'd-mesh softening .......',my_vp%dsoft,ch10, &
1909 &   'phi softening ..........',my_vp%phisoft,ch10, &
1910 &   'q-mesh cut-off .........',my_vp%qcut,ch10, &
1911 &   'q-mesh ratio ...........',my_vp%qratio,ch10, &
1912 &   'r-mesh cut-off .........',my_vp%rcut,ch10, &
1913 &   'r-mesh softening .......',my_vp%rsoft,ch10, &
1914 &   'g-mesh cut-off .........',my_vp%gcut,ch10, &
1915 &   'a-mesh min cut-off .....',my_vp%acutmin,ch10, &
1916 &   'a-mesh ratio ...........',my_vp%aratio,ch10, &
1917 &   'a-mesh delta max .......',my_vp%damax,ch10, &
1918 &   'a-mesh delta min .......',my_vp%damin,ch10, &
1919 &   'global tolerance .......',my_vp%tolerance,ch10, &
1920 &   'tweaks .................',my_vp%tweaks,ch10, &
1921 &   ch10, &
1922 &   '====================================',ch10
1923 
1924   call wrtout(unt,msg,'COLL')
1925 
1926 end subroutine xc_vdw_show

ABINIT/m_xc_vdw/xc_vdw_status [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_status

FUNCTION

  Returns the status of the main vdW-DF switch.

INPUTS

  None

SOURCE

1941 function xc_vdw_status()
1942 
1943 !Arguments ------------------------------------
1944 
1945 !Local variables ------------------------------
1946   logical :: xc_vdw_status
1947 
1948 ! *************************************************************************
1949 
1950   xc_vdw_status = vdw_switch
1951 
1952 end function xc_vdw_status

ABINIT/m_xc_vdw/xc_vdw_trigger [ Functions ]

[ Top ] [ Functions ]

NAME

  xc_vdw_trigger

FUNCTION

  Switches on and off the calculation of vdW interactions.

INPUTS

  condition= boolean condition to trigger the calculations

SOURCE

1967 subroutine xc_vdw_trigger(condition)
1968 
1969 !Arguments ------------------------------------
1970   logical, intent(in) :: condition
1971 
1972 !Local variables ------------------------------
1973   character(len=512) :: msg
1974 
1975 ! *************************************************************************
1976 
1977   if ( my_vdw_params%functional /= 0 ) then
1978 
1979     if ( vdw_switch ) then
1980       if (.not. condition) then
1981         write (msg,'(1x,a,a)') &
1982 &         "[vdW-DF] xc_vdw_trigger: disabling xc_vdw_aggregate",ch10
1983         vdw_switch = condition
1984       else
1985       write (msg,'(1x,a,a)') &
1986 &       "[vdW-DF] xc_vdw_trigger: keeping xc_vdw_aggregate enabled",ch10
1987       end if
1988     else
1989       if ( condition ) then
1990         write (msg,'(1x,a,a)') &
1991 &         "[vdW-DF] xc_vdw_trigger: enabling xc_vdw_aggregate",ch10
1992       else
1993         write (msg,'(1x,a,a)') &
1994 &         "[vdW-DF] xc_vdw_trigger: disabling xc_vdw_aggregate",ch10
1995       end if
1996 
1997       vdw_switch = condition
1998     end if
1999 
2000     call wrtout(std_out,msg,'COLL')     
2001   end if
2002 
2003 end subroutine xc_vdw_trigger

ABINIT/m_xc_vdw/xc_vdw_type [ Types ]

[ Top ] [ Types ]

NAME

  xc_vdw_type

FUNCTION

  Tunable parameters for calculations relying upon van der Waals XC.

SOURCE

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

ABINIT/m_xc_vdw/xc_vdw_write [ Functions ]

[ 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

SOURCE

2021 subroutine xc_vdw_write(filename)
2022 
2023 !Arguments ------------------------------------
2024   character(len=*), intent(in)  :: filename
2025 
2026 !Local variables ------------------------------
2027   character(len=512) :: msg
2028   integer :: ncid,dimids(4),varids(13)
2029   integer :: nderivs_id,npoly_id,ndpts_id,ngpts_id,nqpts_id,nrpts_id
2030 
2031 ! *************************************************************************
2032 
2033   DBG_ENTER("COLL")
2034 
2035 #if defined HAVE_NETCDF
2036   write(msg,'(a,1x,"Writing vdW-DF data to",1x,a)') ch10,trim(filename)
2037   call wrtout(std_out,msg,'COLL')
2038 
2039   ! Create file (overwriting any existing data)
2040   NETCDF_VDWXC_CHECK(nf90_create(trim(filename),NF90_CLOBBER,ncid))
2041 
2042   ! Write attributes
2043   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'file_format',vdw_format_string))
2044   write(msg,'(a,1x,a,1x,"for",1x,a,1x,"(build",1x,a,")")') &
2045 &   PACKAGE_NAME,ABINIT_VERSION,ABINIT_TARGET,ABINIT_VERSION_BUILD
2046   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'generator',msg))
2047   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'acutmin',my_vdw_params%acutmin))
2048   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'aratio',my_vdw_params%aratio))
2049   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damax',my_vdw_params%damax))
2050   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'damin',my_vdw_params%damin))
2051   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dcut',my_vdw_params%dcut))
2052   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dratio',my_vdw_params%dratio))
2053   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'dsoft',my_vdw_params%dsoft))
2054   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'functional',my_vdw_params%functional))
2055   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'gcut',my_vdw_params%gcut))
2056   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'nsmooth',my_vdw_params%nsmooth))
2057   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'phisoft',my_vdw_params%phisoft))
2058   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qcut',my_vdw_params%qcut))
2059   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'qratio',my_vdw_params%qratio))
2060   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rcut',my_vdw_params%rcut))
2061   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'rsoft',my_vdw_params%rsoft))
2062   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tolerance',my_vdw_params%tolerance))
2063   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'tweaks',my_vdw_params%tweaks))
2064   NETCDF_VDWXC_CHECK(nf90_put_att(ncid,NF90_GLOBAL,'zab',my_vdw_params%zab))
2065 
2066   ! Define dimensions
2067   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nderivs',4,nderivs_id))
2068   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'npoly',2,npoly_id))
2069   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ndpts',my_vdw_params%ndpts,ndpts_id))
2070   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'ngpts',my_vdw_params%ngpts,ngpts_id))
2071   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nqpts',my_vdw_params%nqpts,nqpts_id))
2072   NETCDF_VDWXC_CHECK(nf90_def_dim(ncid,'nrpts',my_vdw_params%nrpts+1,nrpts_id))
2073 
2074   ! Define qmesh
2075   dimids(1) = nqpts_id
2076   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qmesh',NF90_DOUBLE,dimids(1),varids(2)))
2077 
2078   ! Define dmesh
2079   dimids(1) = ndpts_id
2080   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'dmesh',NF90_DOUBLE,dimids(1),varids(3)))
2081 
2082   ! Define phi
2083   dimids(1) = ndpts_id
2084   dimids(2) = ndpts_id
2085   dimids(3) = nderivs_id
2086   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi',NF90_DOUBLE,dimids(1:3),varids(4)))
2087 
2088   ! Define phi_u
2089   dimids(1) = ndpts_id
2090   dimids(2) = ndpts_id
2091   dimids(3) = nderivs_id
2092   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u',NF90_DOUBLE,dimids(1:3),varids(5)))
2093 
2094   ! Define phi_bicubic
2095   dimids(1) = nderivs_id
2096   dimids(2) = nderivs_id
2097   dimids(3) = ndpts_id
2098   dimids(4) = ndpts_id
2099   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_bicubic',NF90_DOUBLE,dimids(1:4),varids(6)))
2100 
2101   ! Define phi_u_bicubic
2102   dimids(1) = nderivs_id
2103   dimids(2) = nderivs_id
2104   dimids(3) = ndpts_id
2105   dimids(4) = ndpts_id
2106   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phi_u_bicubic',NF90_DOUBLE,dimids(1:4),varids(7)))
2107 
2108   ! Define phir
2109   dimids(1) = nrpts_id
2110   dimids(2) = nqpts_id
2111   dimids(3) = nqpts_id
2112   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir',NF90_DOUBLE,dimids(1:3),varids(8)))
2113 
2114   ! Define phir_u
2115   dimids(1) = nrpts_id
2116   dimids(2) = nqpts_id
2117   dimids(3) = nqpts_id
2118   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phir_u',NF90_DOUBLE,dimids(1:3),varids(9)))
2119 
2120   ! Define d2phidr2
2121   dimids(1) = nrpts_id
2122   dimids(2) = nqpts_id
2123   dimids(3) = nqpts_id
2124   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidr2',NF90_DOUBLE,dimids(1:3),varids(10)))
2125 
2126   ! Define phig
2127   dimids(1) = nrpts_id
2128   dimids(2) = nqpts_id
2129   dimids(3) = nqpts_id
2130   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'phig',NF90_DOUBLE,dimids(1:3),varids(11)))
2131 
2132   ! Define d2phidg2
2133   dimids(1) = nrpts_id
2134   dimids(2) = nqpts_id
2135   dimids(3) = nqpts_id
2136   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'d2phidg2',NF90_DOUBLE,dimids(1:3),varids(12)))
2137 
2138   ! Define qpoly_basis
2139   dimids(1) = nqpts_id
2140   dimids(2) = nqpts_id
2141   dimids(3) = npoly_id
2142   NETCDF_VDWXC_CHECK(nf90_def_var(ncid,'qpoly_basis',NF90_DOUBLE,dimids(1:3),varids(13)))
2143 
2144   ! Stop definitions
2145   NETCDF_VDWXC_CHECK(nf90_enddef(ncid))
2146 
2147   ! Write variables
2148   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(2),qmesh))
2149   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(3),dmesh))
2150   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(4),phi))
2151   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(5),phi_u))
2152   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(6),phi_bicubic))
2153   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(7),phi_u_bicubic))
2154   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(8),phir))
2155   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(9),phir_u))
2156   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(10),d2phidr2))
2157   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(11),phig))
2158   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(12),d2phidg2))
2159   NETCDF_VDWXC_CHECK(nf90_put_var(ncid,varids(13),qpoly_basis))
2160 
2161   ! Close file
2162   NETCDF_VDWXC_CHECK(nf90_close(ncid))
2163 #else
2164   ABI_WARNING('writing vdW-DF variables requires NetCDF - skipping')
2165 #endif
2166 
2167   DBG_EXIT("COLL")
2168 
2169 end subroutine xc_vdw_write

m_xc_vdw/vdw_df_netcdf_ioerr [ Functions ]

[ Top ] [ m_xc_vdw ] [ Functions ]

NAME

  vdw_df_netcdf_ioerr

FUNCTION

  Reports errors occurring when allocating memory.

INPUTS

  array_name= name of the array not properly allocated.
  array_size= size of the array.
  file_name= name of the file where the allocation was performed.
  file_line= line number in the file where the allocation was performed.

NOTES

  This routine is usually interfaced with the macros defined in abi_common.h
  and uses this information to define a line offset.

SOURCE

3847 subroutine vdw_df_netcdf_ioerr(ncerr,file_name,file_line)
3848 
3849 !Arguments ------------------------------------
3850  integer,intent(in) :: ncerr
3851  character(len=*),optional,intent(in) :: file_name
3852  integer,optional,intent(in) :: file_line
3853 
3854 !Local variables-------------------------------
3855  character(len=1024) :: msg
3856  character(len=fnlen) :: my_file_name = "N/D"
3857  integer :: my_file_line = -1
3858 
3859 ! *************************************************************************
3860 
3861   if ( present(file_name) ) then
3862     my_file_name = trim(file_name)
3863   end if
3864   if ( present(file_line) ) then
3865     my_file_line = file_line
3866   end if
3867 
3868 #if defined HAVE_NETCDF
3869   write(msg,'(a,a,3x,a)') &
3870 &  'NetCDF returned the error:',ch10,trim(nf90_strerror(ncerr))
3871   if ( ncerr /= NF90_NOERR ) then
3872     call die(msg,trim(my_file_name),my_file_line)
3873   end if
3874 #endif
3875 
3876 end subroutine vdw_df_netcdf_ioerr

m_xc_vdw/vdw_df_saturation [ Functions ]

[ Top ] [ m_xc_vdw ] [ Functions ]

NAME

  vdw_df_saturation

FUNCTION

  Obtain a saturated value of q0 which is always inside
  the interpolation range (eq. 5 RS09.)
  * q0s = q0c * (1 - exp(-Sum_i=1,ns 1/i (q0/q0c)**i))
  * dq0sdq0 = dq0s / dq |q=q0

INPUTS

  q0 = q0 value to be saturated  
  q0c = Saturation value, max q0 value  

 OUTPUTS 
  q0s(1) = saturated q0 value
  q0s(2) = derivative of q0s(1) wrt q at q0 

NOTES

  This routine is usually interfaced with the macros defined in abi_common.h
  and uses this information to define a line offset.

SOURCE

3800 subroutine vdw_df_saturation(q0,q0c,q0s)
3801 
3802 !Arguments ------------------------------------
3803  real(dp),intent(in) :: q0, q0c
3804  real(dp),intent(out) :: q0s(2)  
3805 !integer,intent(in) :: ncerr
3806 !character(len=*),optional,intent(in) :: file_name
3807 !integer,optional,intent(in) :: file_line
3808 
3809 !Local variables-------------------------------
3810  integer  :: is,ns
3811  real(dp) :: ptm, dptmdq
3812 ! *************************************************************************
3813 
3814  ns = my_vdw_params%nsmooth
3815 
3816  ptm = (q0 / q0c) / ns
3817  dptmdq = 1._dp / q0c
3818  do is=ns-1,1,-1
3819    ptm = (ptm + 1._dp / is) * q0 / q0c
3820    dptmdq = (dptmdq * q0 + 1) / q0c ! (dptmdq * q0 + one) / q0c
3821  end do 
3822  q0s(1) = q0c * (1 - exp(-ptm)) !q0c * (one - exp(-ptm))
3823  q0s(2) = q0c * dptmdq * exp(-ptm)
3824 
3825 end subroutine vdw_df_saturation