xc_vdW_DF.f90 71.5 KB
Newer Older
1
! Copyright (C) 2001-2009 Quantum ESPRESSO group
thonhauser's avatar
thonhauser committed
2
! Copyright (C) 2015 Brian Kolb, Timo Thonhauser
3 4 5 6
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
thonhauser's avatar
thonhauser committed
7
! ----------------------------------------------------------------------
8 9 10 11 12


MODULE vdW_DF


thonhauser's avatar
thonhauser committed
13 14 15 16 17 18 19 20 21 22 23 24 25 26 27
! This module calculates the non-local correlation contribution to the
! energy and potential according to
!
!    M. Dion, H. Rydberg, E. Schroeder, D. C. Langreth, and
!    B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
!
! henceforth referred to as DION. Further information about the
! functional and its corresponding potential can be found in:
!
!    T. Thonhauser, V.R. Cooper, S. Li, A. Puzder, P. Hyldgaard,
!    and D.C. Langreth, Phys. Rev. B 76, 125112 (2007).
!
! The proper spin extension of vdW-DF, i.e. svdW-DF, is derived in
!
!    T. Thonhauser, S. Zuluaga, C.A. Arter, K. Berland, E. Schroder,
thonhauser's avatar
thonhauser committed
28
!    and P. Hyldgaard, Phys. Rev. Lett. 115, 136402 (2015).
thonhauser's avatar
thonhauser committed
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
!
! henceforth referred to as THONHAUSER.
!
!
! Two review article show many of the vdW-DF applications:
!
!    D. C. Langreth et al., J. Phys.: Condens. Matter 21, 084203 (2009).
!
!    K. Berland et al, Rep. Prog. Phys. 78, 066501 (2015).
!
!
! The method implemented is based on the method of G. Roman-Perez and
! J. M. Soler described in:
!
!    G. Roman-Perez and J. M. Soler, PRL 103, 096102 (2009).
!
! henceforth referred to as SOLER.
!
!
! There are a number of subroutines in this file. All are used only by
! other subroutines here except for the xc_vdW_DF subroutine, which is
! the driver routine for the vdW-DF calculations and is called from
! v_of_rho. This routine handles setting up the parallel run (if any)
! and carries out the calls necessary to calculate the non-local
! correlation contributions to the energy and potential.
54 55


thonhauser's avatar
thonhauser committed
56 57 58 59
USE kinds,             ONLY : dp
USE constants,         ONLY : pi, e2
USE kernel_table,      ONLY : q_mesh, Nr_points, Nqs, r_max, q_cut, q_min, kernel, d2phi_dk2, dk
USE mp,                ONLY : mp_bcast, mp_sum, mp_barrier
60
USE mp_bands,          ONLY : intra_bgrp_comm
thonhauser's avatar
thonhauser committed
61 62 63 64
USE io_global,         ONLY : stdout, ionode
USE fft_base,          ONLY : dfftp
USE fft_interfaces,    ONLY : fwfft, invfft
USE control_flags,     ONLY : iverbosity, gamma_only
65

thonhauser's avatar
thonhauser committed
66
implicit none
67

thonhauser's avatar
thonhauser committed
68 69
REAL(DP), PARAMETER :: epsr =1.d-12  ! a small number to cut off densities
integer :: vdw_type = 1
70

thonhauser's avatar
thonhauser committed
71 72 73
private
public  :: xc_vdW_DF, xc_vdW_DF_spin, stress_vdW_DF, interpolate_kernel, &
           vdw_type, numerical_gradient, initialize_spline_interpolation
74

thonhauser's avatar
thonhauser committed
75
CONTAINS
76 77 78 79 80 81 82 83








thonhauser's avatar
thonhauser committed
84 85 86 87 88 89
  ! ####################################################################
  !                           |             |
  !                           |  functions  |
  !                           |_____________|
  !
  ! Functions to be used in get_q0_on_grid and get_q0_on_grid_spin().
90

thonhauser's avatar
thonhauser committed
91
  function Fs(s)
92

thonhauser's avatar
thonhauser committed
93 94 95 96
     implicit none
     real(dp) :: s, Fs, Z_ab
     real(dp) :: fa=0.1234D0, fb=17.33D0, fc=0.163D0  ! Reparameterized values
                                                      ! from JCTC 5, 2745 (2009).
97

thonhauser's avatar
thonhauser committed
98 99 100 101 102 103 104 105 106
     if(vdw_type == 4) then
         Fs = ( 1 + 15.0D0*fa*s**2 + fb*s**4 + fc*s**6 )**(1.0D0/15.0D0)
     else
         ! ------------------------------------------------------------
         ! Original functional choice for Fs, as definded in DION
         if (vdw_type == 1) Z_ab = -0.8491D0
         if (vdw_type == 2) Z_ab = -1.887D0
         Fs = 1.0D0 - Z_ab * s**2 / 9.0D0
     end if
107

thonhauser's avatar
thonhauser committed
108
  end function Fs
109 110 111 112




thonhauser's avatar
thonhauser committed
113
  function dFs_ds(s)
114

thonhauser's avatar
thonhauser committed
115 116 117 118
     implicit none
     real(dp) :: s, dFs_ds, Z_ab
     real(dp) :: fa=0.1234D0, fb=17.33D0, fc=0.163D0  ! Reparameterized values
                                                      ! from JCTC 5, 2745 (2009).
119

thonhauser's avatar
thonhauser committed
120 121 122 123 124 125 126 127 128 129
     if(vdw_type == 4) then
         dFs_ds = ( 30.0D0*fa*s + 4.0D0*fb*s**3 + 6.0D0*fc*s**5 ) &
                / ( 15.0D0*( 1.0D0 + 15.0D0*fa*s**2 + fb*s**4 + fc*s**6 )**(14.0D0/15.0D0) )
     else
         ! ------------------------------------------------------------
         ! Original functional choice for Fs, as definded in DION
         if (vdw_type == 1) Z_ab = -0.8491D0
         if (vdw_type == 2) Z_ab = -1.887D0
         dFs_ds =  -2.0D0 * s * Z_ab / 9.0D0
     end if
130

thonhauser's avatar
thonhauser committed
131
  end function dFs_ds
132 133 134 135




thonhauser's avatar
thonhauser committed
136
  function kF(rho)
137

thonhauser's avatar
thonhauser committed
138 139
     implicit none
     real(dp) :: rho, kF
140

thonhauser's avatar
thonhauser committed
141
     kF = ( 3.0D0 * pi**2 * rho )**(1.0D0/3.0D0)
142

thonhauser's avatar
thonhauser committed
143
  end function kF
144 145 146 147




thonhauser's avatar
thonhauser committed
148
  function dkF_drho(rho)
149

thonhauser's avatar
thonhauser committed
150 151
     implicit none
     real(dp) :: rho, dkF_drho
152

thonhauser's avatar
thonhauser committed
153
     dkF_drho = (1.0D0/3.0D0) * kF(rho) / rho
154

thonhauser's avatar
thonhauser committed
155
  end function dkF_drho
156 157 158 159




thonhauser's avatar
thonhauser committed
160
  function ds_drho(rho, s)
161

thonhauser's avatar
thonhauser committed
162 163
     implicit none
     real(dp) :: rho, s, ds_drho
164

thonhauser's avatar
thonhauser committed
165
     ds_drho = -s * ( dkF_drho(rho) / kF(rho) + 1.0D0 / rho )
166

thonhauser's avatar
thonhauser committed
167
  end function ds_drho
168 169 170 171




thonhauser's avatar
thonhauser committed
172
  function ds_dgradrho(rho)
173

thonhauser's avatar
thonhauser committed
174 175
     implicit none
     real(dp) :: rho, ds_dgradrho
176

thonhauser's avatar
thonhauser committed
177
     ds_dgradrho = 1.0D0 / (2.0D0 * kF(rho) * rho)
178

thonhauser's avatar
thonhauser committed
179
  end function ds_dgradrho
180 181 182 183




thonhauser's avatar
thonhauser committed
184
  function dqx_drho(rho, s)
185

thonhauser's avatar
thonhauser committed
186 187
     implicit none
     real(dp) :: rho, s, dqx_drho
188

thonhauser's avatar
thonhauser committed
189
     dqx_drho = dkF_drho(rho) * Fs(s) + kF(rho) * dFs_ds(s) * ds_drho(rho, s)
190

thonhauser's avatar
thonhauser committed
191
  end function dqx_drho
192 193 194 195 196 197 198 199








thonhauser's avatar
thonhauser committed
200 201 202 203
  ! ####################################################################
  !                           |             |
  !                           |  XC_VDW_DF  |
  !                           |_____________|
204

thonhauser's avatar
thonhauser committed
205
  SUBROUTINE xc_vdW_DF (rho_valence, rho_core, etxc, vtxc, v)
206

thonhauser's avatar
thonhauser committed
207 208
  USE gvect,                 ONLY : ngm, nl, g, nlm
  USE cell_base,             ONLY : omega, tpiba
209

thonhauser's avatar
thonhauser committed
210
  implicit none
211

thonhauser's avatar
thonhauser committed
212 213 214 215 216 217
  ! --------------------------------------------------------------------
  ! Local variables
  !                                               _
  real(dp), intent(IN) :: rho_valence(:,:)       !
  real(dp), intent(IN) :: rho_core(:)            !  PWSCF input variables
  real(dp), intent(inout) :: etxc, vtxc, v(:,:)  !_
218 219


thonhauser's avatar
thonhauser committed
220 221
  integer :: i_grid, theta_i, i_proc             ! Indexing variables over grid points,
                                                 ! theta functions, and processors.
222

thonhauser's avatar
thonhauser committed
223
  real(dp) :: grid_cell_volume                   ! The volume of the unit cell per G-grid point.
224

thonhauser's avatar
thonhauser committed
225 226 227
  real(dp), allocatable ::  q0(:)                ! The saturated value of q (equations 11 and 12
                                                 ! of DION). This saturation is that of
                                                 ! equation 5 in SOLER.
228

thonhauser's avatar
thonhauser committed
229 230 231
  real(dp), allocatable :: grad_rho(:,:)         ! The gradient of the charge density. The
                                                 ! format is as follows:
                                                 ! grad_rho(grid_point, cartesian_component).
232

thonhauser's avatar
thonhauser committed
233
  real(dp), allocatable :: potential(:)          ! The vdW contribution to the potential
234

thonhauser's avatar
thonhauser committed
235 236 237 238
  real(dp), allocatable :: dq0_drho(:)           ! The derivative of the saturated q0
                                                 ! (equation 5 of SOLER) with respect
                                                 ! to the charge density (see
                                                 ! get_q0_on_grid subroutine for details).
239

thonhauser's avatar
thonhauser committed
240 241 242 243
  real(dp), allocatable :: dq0_dgradrho(:)       ! The derivative of the saturated q0
                                                 ! (equation 5 of SOLER) with respect
                                                 ! to the gradient of the charge density
                                                 ! (again, see get_q0_on_grid subroutine).
244

thonhauser's avatar
thonhauser committed
245 246 247 248 249 250
  complex(dp), allocatable :: thetas(:,:)        ! These are the functions of equation 8 of
                                                 ! SOLER. They will be forward Fourier transformed
                                                 ! in place to get theta(k) and worked on in
                                                 ! place to get the u_alpha(r) of equation 11
                                                 ! in SOLER. They are formatted as follows:
                                                 ! thetas(grid_point, theta_i).
251

thonhauser's avatar
thonhauser committed
252
  real(dp) :: Ec_nl                              ! The non-local vdW contribution to the energy.
253

thonhauser's avatar
thonhauser committed
254 255 256
  real(dp), allocatable :: total_rho(:)          ! This is the sum of the valence and core
                                                 ! charge. This just holds the piece assigned
                                                 ! to this processor.
257

thonhauser's avatar
thonhauser committed
258 259
  logical, save :: first_iteration = .TRUE.      ! Whether this is the first time this
                                                 ! routine has been called.
260 261 262 263




thonhauser's avatar
thonhauser committed
264 265
  ! --------------------------------------------------------------------
  ! Check that the requested non-local functional is implemented.
266

thonhauser's avatar
thonhauser committed
267
  if ( vdW_type /= 1 .AND. vdW_type /= 2) call errore('xc_vdW_DF','E^nl_c not implemented',1)
268 269


thonhauser's avatar
thonhauser committed
270 271
  ! --------------------------------------------------------------------
  ! Write out the vdW-DF imformation.
272

thonhauser's avatar
thonhauser committed
273 274
  if ( ionode .AND. first_iteration ) call vdW_info
  first_iteration = .FALSE.
275 276


thonhauser's avatar
thonhauser committed
277 278 279
  ! --------------------------------------------------------------------
  ! Allocate arrays. nnr is a PWSCF variable that holds the number of
  ! points assigned to a given processor.
280

thonhauser's avatar
thonhauser committed
281 282
  allocate( q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr), grad_rho(dfftp%nnr, 3) )
  allocate( total_rho(dfftp%nnr), potential(dfftp%nnr), thetas(dfftp%nnr, Nqs) )
283 284


thonhauser's avatar
thonhauser committed
285 286 287 288 289
  ! --------------------------------------------------------------------
  ! Add together the valence and core charge densities to get the total
  ! charge density. Note that rho_core is not the true core density and
  ! it is only non-zero for pseudopotentials with non-local core
  ! corrections.
290

thonhauser's avatar
thonhauser committed
291
  total_rho = rho_valence(:,1) + rho_core(:)
292 293


thonhauser's avatar
thonhauser committed
294 295
  ! --------------------------------------------------------------------
  ! Here we calculate the gradient in reciprocal space using FFT.
296

thonhauser's avatar
thonhauser committed
297
  call numerical_gradient (total_rho, grad_rho)
298 299


thonhauser's avatar
thonhauser committed
300 301 302 303 304 305 306
  ! --------------------------------------------------------------------
  ! Find the value of q0 for all assigned grid points. q is defined in
  ! equations 11 and 12 of DION and q0 is the saturated version of q
  ! defined in equation 5 of SOLER. This routine also returns the
  ! derivatives of the q0s with respect to the charge-density and the
  ! gradient of the charge-density. These are needed for the potential
  ! calculated below. This routine also calculates the thetas.
307

thonhauser's avatar
thonhauser committed
308
  CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
309 310


thonhauser's avatar
thonhauser committed
311 312 313 314 315
  ! --------------------------------------------------------------------
  ! Carry out the integration in equation 8 of SOLER. This also turns
  ! the thetas array into the precursor to the u_i(k) array which is
  ! inverse fourier transformed to get the u_i(r) functions of SOLER
  ! equation 11. Add the energy we find to the output variable etxc.
316

thonhauser's avatar
thonhauser committed
317 318
  call vdW_energy (thetas, Ec_nl)
  etxc = etxc + Ec_nl
319

thonhauser's avatar
thonhauser committed
320
  if (iverbosity > 0) then
321
     call mp_sum(Ec_nl, intra_bgrp_comm)
thonhauser's avatar
thonhauser committed
322 323 324 325 326 327
     if (ionode) then
        write(stdout,'(/ / A)')       "     -----------------------------------------------"
        write(stdout,'(A, F15.8, A)') "     Non-local corr. energy    =  ", Ec_nl, " Ry"
        write(stdout,'(A /)')         "     -----------------------------------------------"
     end if
  end if
328 329


thonhauser's avatar
thonhauser committed
330 331 332 333 334 335 336
  ! --------------------------------------------------------------------
  ! Here we calculate the potential. This is calculated via equation 10
  ! of SOLER, using the u_i(r) calculated from quations 11 and 12 of
  ! SOLER. Each processor allocates the array to be the size of the full
  ! grid  because, as can be seen in SOLER equation 10, processors need
  ! to access grid points outside their allocated regions. Begin by
  ! FFTing the u_i(k) to get the u_i(r) of SOLER equation 11.
337

thonhauser's avatar
thonhauser committed
338 339 340
  do theta_i = 1, Nqs
     CALL invfft('Dense', thetas(:,theta_i), dfftp)
  end do
341

thonhauser's avatar
thonhauser committed
342 343
  call get_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, thetas, potential)
  v(:,1) = v(:,1) + e2 * potential(:)
344 345


thonhauser's avatar
thonhauser committed
346 347
  ! --------------------------------------------------------------------
  ! The integral of rho(r)*potential(r) for the vtxc output variable.
348

thonhauser's avatar
thonhauser committed
349
  grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
350

thonhauser's avatar
thonhauser committed
351 352 353
  do i_grid = 1, dfftp%nnr
     vtxc = vtxc + e2 * grid_cell_volume * rho_valence(i_grid,1) * potential(i_grid)
  end do
354

thonhauser's avatar
thonhauser committed
355
  deallocate ( potential, q0, grad_rho, dq0_drho, dq0_dgradrho, total_rho, thetas )
356

thonhauser's avatar
thonhauser committed
357
  END SUBROUTINE xc_vdW_DF
358 359 360 361 362 363 364 365








thonhauser's avatar
thonhauser committed
366 367 368 369 370 371 372
  ! ####################################################################
  !                          |                  |
  !                          |  XC_VDW_DF_spin  |
  !                          |__________________|
  !
  ! This subroutine is as similar to xc_vdW_DF() as possible, but
  ! handles the collinear nspin=2 case.
373

thonhauser's avatar
thonhauser committed
374
  SUBROUTINE xc_vdW_DF_spin (rho_valence, rho_core, etxc, vtxc, v)
375

thonhauser's avatar
thonhauser committed
376 377
  USE gvect,                 ONLY : ngm, nl, g, nlm
  USE cell_base,             ONLY : omega, tpiba
378

thonhauser's avatar
thonhauser committed
379
  implicit none
380

thonhauser's avatar
thonhauser committed
381 382 383 384 385 386
  ! --------------------------------------------------------------------
  ! Local variables
  !                                              _
  real(dp), intent(IN) :: rho_valence(:,:)      !
  real(dp), intent(IN) :: rho_core(:)           !  PWSCF input variables
  real(dp), intent(inout) :: etxc, vtxc, v(:,:) !_
387 388


thonhauser's avatar
thonhauser committed
389 390 391
  integer :: i_grid, theta_i, i_proc            ! Indexing variables over grid points,
                                                ! theta functions, and processors, and a
                                                ! generic index.
392

thonhauser's avatar
thonhauser committed
393
  real(dp) :: grid_cell_volume                  ! The volume of the unit cell per G-grid point.
394

thonhauser's avatar
thonhauser committed
395 396 397
  real(dp), allocatable :: q0(:)                ! The saturated value of q (equations 11 and 12
                                                ! of DION). This saturation is that of
                                                ! equation 5 in SOLER.
398

thonhauser's avatar
thonhauser committed
399 400 401 402 403 404 405 406 407
  real(dp), allocatable :: grad_rho(:,:)        ! The gradient of the charge density. The
                                                ! format is as follows:
                                                ! grad_rho(grid_point, cartesian_component)
  real(dp), allocatable :: grad_rho_up(:,:)     ! The gradient of the up charge density. The
                                                ! format is as follows:
                                                ! grad_rho(grid_point, cartesian_component)
  real(dp), allocatable :: grad_rho_down(:,:)   ! The gradient of the down charge density. The
                                                ! format is as follows:
                                                ! grad_rho(grid_point, cartesian_component)
408

thonhauser's avatar
thonhauser committed
409 410
  real(dp), allocatable :: potential_up(:)      ! The vdW contribution to the potential.
  real(dp), allocatable :: potential_down(:)    ! The vdW contribution to the potential.
411

thonhauser's avatar
thonhauser committed
412 413 414 415
  real(dp), allocatable :: dq0_drho_up(:)       ! The derivative of the saturated q0
  real(dp), allocatable :: dq0_drho_down(:)     ! (equation 5 of SOLER) with respect
                                                ! to the charge density (see
                                                ! get_q0_on_grid subroutine for details).
416

thonhauser's avatar
thonhauser committed
417 418 419 420
  real(dp), allocatable :: dq0_dgradrho_up(:)   ! The derivative of the saturated q0
  real(dp), allocatable :: dq0_dgradrho_down(:) ! (equation 5 of SOLER) with respect
                                                ! to the gradient of the charge density
                                                ! (again, see get_q0_on_grid subroutine).
421

thonhauser's avatar
thonhauser committed
422 423 424 425 426 427
  complex(dp), allocatable :: thetas(:,:)       ! These are the functions of equation 8 of
                                                ! SOLER. They will be forward Fourier transformed
                                                ! in place to get theta(k) and worked on in
                                                ! place to get the u_alpha(r) of equation 11
                                                ! in SOLER. They are formatted as follows:
                                                ! thetas(grid_point, theta_i).
428

thonhauser's avatar
thonhauser committed
429
  real(dp) :: Ec_nl                             ! The non-local vdW contribution to the energy.
430

thonhauser's avatar
thonhauser committed
431 432 433 434 435 436 437 438 439
  real(dp), allocatable :: total_rho(:)         ! This is the sum of the valence (up and down)
                                                ! and core charge. This just holds the piece
                                                ! assigned to this processor.
  real(dp), allocatable :: rho_up(:)            ! This is the just the up valence charge.
                                                ! This just holds the piece assigned
                                                ! to this processor.
  real(dp), allocatable :: rho_down(:)          ! This is the just the down valence charge.
                                                ! This just holds the piece assigned
                                                ! to this processor.
440

thonhauser's avatar
thonhauser committed
441 442
  logical, save :: first_iteration = .TRUE.     ! Whether this is the first time this
                                                ! routine has been called.
443 444 445 446




thonhauser's avatar
thonhauser committed
447 448
  ! --------------------------------------------------------------------
  ! Check that the requested non-local functional is implemented.
449

thonhauser's avatar
thonhauser committed
450
  if ( vdW_type /= 1 .AND. vdW_type /= 2) call errore('xc_vdW_DF','E^nl_c not implemented',1)
451 452


thonhauser's avatar
thonhauser committed
453 454
  ! --------------------------------------------------------------------
  ! Write out the vdW-DF imformation.
455

thonhauser's avatar
thonhauser committed
456 457
  if ( ionode .AND. first_iteration ) call vdW_info
  first_iteration = .FALSE.
458 459


thonhauser's avatar
thonhauser committed
460 461 462
  ! --------------------------------------------------------------------
  ! Allocate arrays. nnr is a PWSCF variable that holds the number of
  ! points assigned to a given processor.
463

thonhauser's avatar
thonhauser committed
464 465 466 467 468 469 470
  allocate( q0(dfftp%nnr), total_rho(dfftp%nnr), grad_rho(dfftp%nnr, 3) )
  allocate( rho_up(dfftp%nnr), rho_down(dfftp%nnr) )
  allocate( dq0_drho_up (dfftp%nnr), dq0_dgradrho_up  (dfftp%nnr) )
  allocate( dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )
  allocate( grad_rho_up(dfftp%nnr, 3), grad_rho_down(dfftp%nnr, 3) )
  allocate( potential_up(dfftp%nnr), potential_down(dfftp%nnr) )
  allocate( thetas(dfftp%nnr, Nqs) )
471 472


thonhauser's avatar
thonhauser committed
473 474 475 476 477
  ! --------------------------------------------------------------------
  ! Add together the valence and core charge densities to get the total
  ! charge density. Note that rho_core is not the true core density and
  ! it is only non-zero for pseudopotentials with non-local core
  ! corrections.
478

thonhauser's avatar
thonhauser committed
479 480 481
  rho_up    = rho_valence(:,1) + 0.5D0*rho_core(:)
  rho_down  = rho_valence(:,2) + 0.5D0*rho_core(:)
  total_rho = rho_up + rho_down
482 483


thonhauser's avatar
thonhauser committed
484 485
  ! --------------------------------------------------------------------
  ! Here we calculate the gradient in reciprocal space using FFT.
486

thonhauser's avatar
thonhauser committed
487 488 489
  call numerical_gradient ( total_rho, grad_rho      )
  call numerical_gradient ( rho_up,    grad_rho_up   )
  call numerical_gradient ( rho_down,  grad_rho_down )
490 491


thonhauser's avatar
thonhauser committed
492 493 494 495 496 497 498 499
  ! --------------------------------------------------------------------
  ! Find the value of q0 for all assigned grid points. q is defined in
  ! equations 11 and 12 of DION and q0 is the saturated version of q
  ! defined in equation 5 of SOLER. In the spin case, q0 is defined by
  ! equation 8 (and text above that equation) of THONHAUSER. This
  ! routine also returns the derivatives of the q0s with respect to the
  ! charge-density and the gradient of the charge-density. These are
  ! needed for the potential calculated below.
500

thonhauser's avatar
thonhauser committed
501 502
  CALL get_q0_on_grid_spin (total_rho, rho_up, rho_down, grad_rho, grad_rho_up, grad_rho_down, &
       q0, dq0_drho_up, dq0_drho_down, dq0_dgradrho_up, dq0_dgradrho_down, thetas)
503 504


thonhauser's avatar
thonhauser committed
505 506 507 508 509 510 511 512 513 514
  ! --------------------------------------------------------------------
  ! Carry out the integration in equation 8 of SOLER. This also turns
  ! the thetas array into the precursor to the u_i(k) array which is
  ! inverse fourier transformed to get the u_i(r) functions of SOLER
  ! equation 11. Add the energy we find to the output variable etxc.

  call vdW_energy(thetas, Ec_nl)
  etxc = etxc + Ec_nl

  if (iverbosity > 0) then
515
     call mp_sum(Ec_nl, intra_bgrp_comm)
thonhauser's avatar
thonhauser committed
516 517 518 519 520
     if (ionode) then
        write(stdout,'(/ / A)')       "     -----------------------------------------------"
        write(stdout,'(A, F15.8, A)') "     Non-local corr. energy    =  ", Ec_nl, " Ry"
        write(stdout,'(A /)')         "     -----------------------------------------------"
     end if
521 522 523
  end if


thonhauser's avatar
thonhauser committed
524 525 526 527 528 529 530
  ! --------------------------------------------------------------------
  ! Here we calculate the potential. This is calculated via equation 10
  ! of SOLER, using the u_i(r) calculated from quations 11 and 12 of
  ! SOLER. Each processor allocates the array to be the size of the full
  ! grid because, as can be seen in SOLER equation 10, processors need
  ! to access grid points outside their allocated regions. Begin by
  ! FFTing the u_i(k) to get the u_i(r) of SOLER equation 11.
531

thonhauser's avatar
thonhauser committed
532 533
  do theta_i = 1, Nqs
     CALL invfft('Dense', thetas(:,theta_i), dfftp)
534 535
  end do

thonhauser's avatar
thonhauser committed
536 537 538 539
  call get_potential (q0, dq0_drho_up  , dq0_dgradrho_up  , grad_rho_up  , thetas, potential_up  )
  call get_potential (q0, dq0_drho_down, dq0_dgradrho_down, grad_rho_down, thetas, potential_down)
  v(:,1) = v(:,1) + e2 * potential_up  (:)
  v(:,2) = v(:,2) + e2 * potential_down(:)
540 541


thonhauser's avatar
thonhauser committed
542 543
  ! --------------------------------------------------------------------
  ! The integral of rho(r)*potential(r) for the vtxc output variable
544

thonhauser's avatar
thonhauser committed
545
  grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
546

thonhauser's avatar
thonhauser committed
547 548 549 550
  do i_grid = 1, dfftp%nnr
     vtxc = vtxc + e2 * grid_cell_volume * rho_valence(i_grid,1) * potential_up  (i_grid) &
                 + e2 * grid_cell_volume * rho_valence(i_grid,2) * potential_down(i_grid)
  end do
551 552


thonhauser's avatar
thonhauser committed
553 554 555
  deallocate( potential_up, potential_down, q0, grad_rho, grad_rho_up, &
              grad_rho_down, dq0_drho_up, dq0_drho_down, thetas, &
              dq0_dgradrho_up, dq0_dgradrho_down, total_rho, rho_up, rho_down )
556

thonhauser's avatar
thonhauser committed
557
  END SUBROUTINE xc_vdW_DF_spin
558 559 560 561 562 563 564 565








thonhauser's avatar
thonhauser committed
566 567 568 569 570 571 572 573 574 575 576 577
  ! ####################################################################
  !                       |                  |
  !                       |  GET_Q0_ON_GRID  |
  !                       |__________________|
  !
  ! This routine first calculates the q value defined in (DION equations
  ! 11 and 12), then saturates it according to (SOLER equation 5). More
  ! specifically it calculates the following:
  !
  !     q0(ir) = q0 as defined above
  !     dq0_drho(ir) = total_rho * d q0 /d rho
  !     dq0_dgradrho = total_rho / |grad_rho| * d q0 / d |grad_rho|
578

thonhauser's avatar
thonhauser committed
579
  SUBROUTINE get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
580

thonhauser's avatar
thonhauser committed
581
  implicit none
582

thonhauser's avatar
thonhauser committed
583
  real(dp),  intent(IN)     :: total_rho(:), grad_rho(:,:)         ! Input variables needed
584

thonhauser's avatar
thonhauser committed
585 586 587
  real(dp),  intent(OUT)    :: q0(:), dq0_drho(:), dq0_dgradrho(:) ! Output variables that have been allocated
                                                                   ! outside this routine but will be set here.
  complex(dp), intent(inout):: thetas(:,:)                         ! The thetas from SOLER.
588

thonhauser's avatar
thonhauser committed
589 590
  integer,   parameter      :: m_cut = 12                          ! How many terms to include in the sum
                                                                   ! of SOLER equation 5.
591

thonhauser's avatar
thonhauser committed
592 593 594 595 596 597
  real(dp)                  :: rho                                 ! Local variable for the density.
  real(dp)                  :: r_s                                 ! Wigner–Seitz radius.
  real(dp)                  :: s                                   ! Reduced gradient.
  real(dp)                  :: q, ec
  real(dp)                  :: dq0_dq                              ! The derivative of the saturated
                                                                   ! q0 with respect to q.
598

thonhauser's avatar
thonhauser committed
599
  integer                   :: i_grid, idx                         ! Indexing variables.
600 601 602 603




thonhauser's avatar
thonhauser committed
604 605
  ! --------------------------------------------------------------------
  ! Initialize q0-related arrays.
606

thonhauser's avatar
thonhauser committed
607 608 609
  q0(:)           = q_cut
  dq0_drho(:)     = 0.0D0
  dq0_dgradrho(:) = 0.0D0
610

thonhauser's avatar
thonhauser committed
611
  do i_grid = 1, dfftp%nnr
612

thonhauser's avatar
thonhauser committed
613
     rho = total_rho(i_grid)
614 615


thonhauser's avatar
thonhauser committed
616 617 618 619 620 621 622
     ! -----------------------------------------------------------------
     ! This prevents numerical problems. If the charge density is
     ! negative (an unphysical situation), we simply treat it as very
     ! small. In that case, q0 will be very large and will be saturated.
     ! For a saturated q0 the derivative dq0_dq will be 0 so we set q0 =
     ! q_cut and dq0_drho = dq0_dgradrho = 0 and go on to the next
     ! point.
623

thonhauser's avatar
thonhauser committed
624
     if ( rho < epsr ) cycle
625 626


thonhauser's avatar
thonhauser committed
627 628
     ! -----------------------------------------------------------------
     ! Calculate some intermediate values needed to find q.
629

thonhauser's avatar
thonhauser committed
630
     r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
631

thonhauser's avatar
thonhauser committed
632 633
     s   = sqrt( grad_rho(i_grid,1)**2 + grad_rho(i_grid,2)**2 + grad_rho(i_grid,3)**2 ) &
         / (2.0D0 * kF(rho) * rho )
634 635


thonhauser's avatar
thonhauser committed
636 637 638
     ! -----------------------------------------------------------------
     ! This is the q value defined in equations 11 and 12 of DION.
     ! Use pw() from flib/functionals.f90 to get qc = kf/eps_x * eps_c.
639

thonhauser's avatar
thonhauser committed
640 641
     call pw(r_s, 1, ec, dq0_drho(i_grid))
     q = -4.0D0*pi/3.0D0 * ec + kF(rho) * Fs(s)
642 643


thonhauser's avatar
thonhauser committed
644 645
     ! -----------------------------------------------------------------
     ! Bring q into its proper bounds.
646

thonhauser's avatar
thonhauser committed
647 648
     CALL saturate_q ( q, q_cut, q0(i_grid), dq0_dq )
     if (q0(i_grid) < q_min) q0(i_grid) = q_min
649 650


thonhauser's avatar
thonhauser committed
651 652 653 654 655 656 657 658 659 660 661 662 663 664 665
     ! -----------------------------------------------------------------
     ! Here we find derivatives. These are actually the density times
     ! the derivative of q0 with respect to rho and grad_rho. The
     ! density factor comes in since we are really differentiating
     ! theta = (rho)*P(q0) with respect to density (or its gradient)
     ! which will be
     !
     !    dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho]
     !
     ! and
     !
     !    dtheta_dgrad_rho = dP_dq0 * [rho * dq0_dq * dq_dgrad_rho]
     !
     ! The parts in square brackets are what is calculated here. The
     ! dP_dq0 term will be interpolated later.
666

thonhauser's avatar
thonhauser committed
667 668 669
     dq0_drho(i_grid)     = dq0_dq * rho * ( -4.0D0*pi/3.0D0 * &
                            (dq0_drho(i_grid) - ec)/rho + dqx_drho(rho, s) )
     dq0_dgradrho(i_grid) = dq0_dq * rho * kF(rho) * dFs_ds(s) * ds_dgradrho(rho)
670

thonhauser's avatar
thonhauser committed
671
  end do
672 673


thonhauser's avatar
thonhauser committed
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693
  ! --------------------------------------------------------------------
  ! Here we calculate the theta functions of SOLER equation 8. These are
  ! defined as
  !
  !    rho * P_i(q0(rho, grad_rho))
  !
  ! where P_i is a polynomial that interpolates a Kroneker delta
  ! function at the point q_i (taken from the q_mesh) and q0 is the
  ! saturated version of q. q is defined in equations 11 and 12 of DION
  ! and the saturation proceedure is defined in equation 5 of SOLER.
  ! This is the biggest memory consumer in the method since the thetas
  ! array is (total # of FFT points)*Nqs complex numbers. In a parallel
  ! run, each processor will hold the values of all the theta functions
  ! on just the points assigned to it. thetas are stored in reciprocal
  ! space as theta_i(k) because this is the way they are used later for
  ! the convolution (equation 8 of SOLER). Start by interpolating the
  ! P_i polynomials defined in equation 3 in SOLER for the particular q0
  ! values we have.

  CALL spline_interpolation (q_mesh, q0, thetas)
694

thonhauser's avatar
thonhauser committed
695 696 697
  do i_grid = 1, dfftp%nnr
     thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
  end do
698

thonhauser's avatar
thonhauser committed
699 700 701
  do idx = 1, Nqs
     CALL fwfft ('Dense', thetas(:,idx), dfftp)
  end do
702

thonhauser's avatar
thonhauser committed
703
  END SUBROUTINE get_q0_on_grid
704 705 706 707 708 709 710 711








thonhauser's avatar
thonhauser committed
712 713 714 715
  ! ####################################################################
  !                       |                       |
  !                       |  GET_Q0_ON_GRID_spin  |
  !                       |_______________________|
716

thonhauser's avatar
thonhauser committed
717 718 719
  SUBROUTINE get_q0_on_grid_spin (total_rho, rho_up, rho_down, grad_rho, &
             grad_rho_up, grad_rho_down, q0, dq0_drho_up, dq0_drho_down, &
             dq0_dgradrho_up, dq0_dgradrho_down, thetas)
720

thonhauser's avatar
thonhauser committed
721
  implicit none
722

thonhauser's avatar
thonhauser committed
723 724 725
  real(dp),  intent(IN)      :: total_rho(:), grad_rho(:,:)              ! Input variables.
  real(dp),  intent(IN)      :: rho_up(:), grad_rho_up(:,:)              ! Input variables.
  real(dp),  intent(IN)      :: rho_down(:), grad_rho_down(:,:)          ! Input variables.
726

thonhauser's avatar
thonhauser committed
727 728 729
  real(dp),  intent(OUT)     :: q0(:), dq0_drho_up(:), dq0_drho_down(:)  ! Output variables.
  real(dp),  intent(OUT)     :: dq0_dgradrho_up(:), dq0_dgradrho_down(:) ! Output variables.
  complex(dp), intent(inout) :: thetas(:,:)                              ! The thetas from SOLER.
730

thonhauser's avatar
thonhauser committed
731 732 733 734 735 736 737 738 739 740 741 742
  real(dp)                   :: rho, up, down                            ! Local copy of densities.
  real(dp)                   :: zeta                                     ! Spin polarization.
  real(dp)                   :: r_s                                      ! Wigner-Seitz radius.
  real(dp)                   :: q, qc, qx, qx_up, qx_down                ! q for exchange and correlation.
  real(dp)                   :: q0x_up, q0x_down                         ! Saturated q values.
  real(dp)                   :: ec, fac
  real(dp)                   :: dq0_dq, dq0x_up_dq, dq0x_down_dq         ! Derivative of q0 w.r.t q.
  real(dp)                   :: dqc_drho_up, dqc_drho_down               ! Intermediate values.
  real(dp)                   :: dqx_drho_up, dqx_drho_down               ! Intermediate values.
  real(dp)                   :: s_up, s_down                             ! Reduced gradients.
  integer                    :: i_grid, idx                              ! Indexing variables
  logical                    :: calc_qx_up, calc_qx_down
743 744 745 746




thonhauser's avatar
thonhauser committed
747
  fac = 2.0D0**(-1.0D0/3.0D0)
748 749


thonhauser's avatar
thonhauser committed
750 751
  ! --------------------------------------------------------------------
  ! Initialize q0-related arrays.
752

thonhauser's avatar
thonhauser committed
753 754 755 756 757
  q0(:)                = q_cut
  dq0_drho_up(:)       = 0.0D0
  dq0_drho_down(:)     = 0.0D0
  dq0_dgradrho_up(:)   = 0.0D0
  dq0_dgradrho_down(:) = 0.0D0
758 759


thonhauser's avatar
thonhauser committed
760
  do i_grid = 1, dfftp%nnr
761

thonhauser's avatar
thonhauser committed
762 763 764
     rho  = total_rho(i_grid)
     up   = rho_up(i_grid)
     down = rho_down(i_grid)
765 766


thonhauser's avatar
thonhauser committed
767 768 769 770 771 772 773
     ! -----------------------------------------------------------------
     ! This prevents numerical problems. If the charge density is
     ! negative (an unphysical situation), we simply treat it as very
     ! small. In that case, q0 will be very large and will be saturated.
     ! For a saturated q0 the derivative dq0_dq will be 0 so we set q0 =
     ! q_cut and dq0_drho = dq0_dgradrho = 0 and go on to the next
     ! point.
774

thonhauser's avatar
thonhauser committed
775
     if ( rho < epsr ) cycle
776

thonhauser's avatar
thonhauser committed
777 778
     calc_qx_up   = .TRUE.
     calc_qx_down = .TRUE.
779

thonhauser's avatar
thonhauser committed
780 781
     if ( up   < epsr/2.0D0 ) calc_qx_up   = .FALSE.
     if ( down < epsr/2.0D0 ) calc_qx_down = .FALSE.
782 783


thonhauser's avatar
thonhauser committed
784 785 786 787 788
     ! -----------------------------------------------------------------
     ! The spin case is numerically even more tricky and we have to
     ! saturate each spin channel separately. Note that we are
     ! saturating at a higher value here, so that very large q values
     ! get saturated to exactly q_cut in the second, overall saturation.
789

thonhauser's avatar
thonhauser committed
790 791 792 793
     q0x_up        = 0.0D0
     q0x_down      = 0.0D0
     dqx_drho_up   = 0.0D0
     dqx_drho_down = 0.0D0
794 795


thonhauser's avatar
thonhauser committed
796 797 798 799 800 801
     if (calc_qx_up) then
        s_up    = sqrt( grad_rho_up(i_grid,1)**2 + grad_rho_up(i_grid,2)**2 + &
                  grad_rho_up(i_grid,3)**2 ) / (2.0D0 * kF(up) * up)
        qx_up   = kF(2.0D0*up) * Fs(fac*s_up)
        CALL saturate_q (qx_up, 4.0D0*q_cut, q0x_up, dq0x_up_dq)
     end if
802

thonhauser's avatar
thonhauser committed
803 804 805 806 807 808
     if (calc_qx_down) then
        s_down  = sqrt( grad_rho_down(i_grid,1)**2 + grad_rho_down(i_grid,2)**2 + &
                  grad_rho_down(i_grid,3)**2) / (2.0D0 * kF(down) * down)
        qx_down = kF(2.0D0*down) * Fs(fac*s_down)
        CALL saturate_q (qx_down, 4.0D0*q_cut, q0x_down, dq0x_down_dq)
     end if
809 810


thonhauser's avatar
thonhauser committed
811 812 813
     ! -----------------------------------------------------------------
     ! This is the q value defined in equations 11 and 12 of DION and
     ! equation 8 of THONHAUSER (also see text above that equation).
814

thonhauser's avatar
thonhauser committed
815 816 817 818
     r_s  = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
     zeta = (up - down) / rho
     IF ( ABS(zeta) > 1.0D0 ) zeta = SIGN(1.0D0, zeta)
     call pw_spin(r_s, zeta, ec, dqc_drho_up, dqc_drho_down)
819

thonhauser's avatar
thonhauser committed
820 821 822
     qx = ( up * q0x_up + down * q0x_down ) / rho
     qc = -4.0D0*pi/3.0D0 * ec
     q  = qx + qc
823 824


thonhauser's avatar
thonhauser committed
825 826
     ! -----------------------------------------------------------------
     ! Bring q into its proper bounds.
827

thonhauser's avatar
thonhauser committed
828 829
     CALL saturate_q (q, q_cut, q0(i_grid), dq0_dq)
     if (q0(i_grid) < q_min) q0(i_grid) = q_min
830 831


thonhauser's avatar
thonhauser committed
832 833 834 835 836 837 838 839 840 841 842 843 844 845 846
     ! -----------------------------------------------------------------
     ! Here we find derivatives. These are actually the density times
     ! the derivative of q0 with respect to rho and grad_rho. The
     ! density factor comes in since we are really differentiating
     ! theta = (rho)*P(q0) with respect to density (or its gradient)
     ! which will be
     !
     !    dtheta_drho = P(q0) + dP_dq0 * [rho * dq0_dq * dq_drho]
     !
     ! and
     !
     !    dtheta_dgrad_rho = dP_dq0 * [rho * dq0_dq * dq_dgrad_rho]
     !
     ! The parts in square brackets are what is calculated here. The
     ! dP_dq0 term will be interpolated later.
847

thonhauser's avatar
thonhauser committed
848 849 850 851 852
     if (calc_qx_up) then
        dqx_drho_up   = 2.0D0*dq0x_up_dq*up*dqx_drho(2.0D0*up, fac*s_up) + q0x_up*down/rho
        dq0_dgradrho_up (i_grid) = 2.0D0 * dq0_dq * dq0x_up_dq * up * kF(2.0D0*up) * &
                        dFs_ds(fac*s_up) * ds_dgradrho(2.0D0*up)
     end if
853

thonhauser's avatar
thonhauser committed
854 855 856 857 858
     if (calc_qx_down) then
        dqx_drho_down = 2.0D0*dq0x_down_dq*down*dqx_drho(2.0D0*down, fac*s_down) + q0x_down*up/rho
        dq0_dgradrho_down(i_grid) = 2.0D0 * dq0_dq * dq0x_down_dq * down * kF(2.0D0*down) * &
                        dFs_ds(fac*s_down) * ds_dgradrho(2.0D0*down)
     end if
859

thonhauser's avatar
thonhauser committed
860 861
     if (calc_qx_down) dqx_drho_up   = dqx_drho_up   - q0x_down*down/rho
     if (calc_qx_up)   dqx_drho_down = dqx_drho_down - q0x_up  *up  /rho
862

thonhauser's avatar
thonhauser committed
863 864
     dqc_drho_up   = -4.0D0*pi/3.0D0 * (dqc_drho_up   - ec)
     dqc_drho_down = -4.0D0*pi/3.0D0 * (dqc_drho_down - ec)
865

thonhauser's avatar
thonhauser committed
866 867
     dq0_drho_up  (i_grid) = dq0_dq * (dqc_drho_up   + dqx_drho_up  )
     dq0_drho_down(i_grid) = dq0_dq * (dqc_drho_down + dqx_drho_down)
868

thonhauser's avatar
thonhauser committed
869
  end do
870 871


thonhauser's avatar
thonhauser committed
872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891
  ! --------------------------------------------------------------------
  ! Here we calculate the theta functions of SOLER equation 8. These are
  ! defined as
  !
  !    rho * P_i(q0(rho, grad_rho))
  !
  ! where P_i is a polynomial that interpolates a Kroneker delta
  ! function at the point q_i (taken from the q_mesh) and q0 is the
  ! saturated version of q. q is defined in equations 11 and 12 of DION
  ! and the saturation proceedure is defined in equation 5 of SOLER.
  ! This is the biggest memory consumer in the method since the thetas
  ! array is (total # of FFT points)*Nqs complex numbers. In a parallel
  ! run, each processor will hold the values of all the theta functions
  ! on just the points assigned to it. thetas are stored in reciprocal
  ! space as theta_i(k) because this is the way they are used later for
  ! the convolution (equation 8 of SOLER). Start by interpolating the
  ! P_i polynomials defined in equation 3 in SOLER for the particular q0
  ! values we have.

  CALL spline_interpolation (q_mesh, q0, thetas)
892

thonhauser's avatar
thonhauser committed
893 894 895
  do i_grid = 1, dfftp%nnr
     thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
  end do
896

thonhauser's avatar
thonhauser committed
897 898 899
  do idx = 1, Nqs
     CALL fwfft ('Dense', thetas(:,idx), dfftp)
  end do
900

thonhauser's avatar
thonhauser committed
901
  END SUBROUTINE get_q0_on_grid_spin
902 903 904 905 906 907 908 909








thonhauser's avatar
thonhauser committed
910 911 912 913
  ! ####################################################################
  !                            |              |
  !                            |  saturate_q  |
  !                            |______________|
914

thonhauser's avatar
thonhauser committed
915
  SUBROUTINE saturate_q (q, q_cut, q0, dq0_dq)
916

thonhauser's avatar
thonhauser committed
917
  implicit none
918

thonhauser's avatar
thonhauser committed
919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130
  real(dp),  intent(IN)      :: q             ! Input q.
  real(dp),  intent(IN)      :: q_cut         ! Cutoff q.
  real(dp),  intent(OUT)     :: q0            ! Output saturated q.
  real(dp),  intent(OUT)     :: dq0_dq        ! Derivative of dq0/dq.

  integer,    parameter      :: m_cut = 12    ! How many terms to include
                                              ! in the sum of SOLER equation 5.

  real(dp)                   :: e             ! Exponent.
  integer                    :: idx           ! Indexing variable.




  ! --------------------------------------------------------------------
  ! Here, we calculate q0 by saturating q according to equation 5 of
  ! SOLER. Also, we find the derivative dq0_dq needed for the
  ! derivatives dq0_drho and dq0_dgradrh0 discussed below.

  e      = 0.0D0
  dq0_dq = 0.0D0

  do idx = 1, m_cut
     e      = e + (q/q_cut)**idx/idx
     dq0_dq = dq0_dq + (q/q_cut)**(idx-1)
  end do

  q0     = q_cut*(1.0D0 - exp(-e))
  dq0_dq = dq0_dq * exp(-e)

  END SUBROUTINE saturate_q








  ! ####################################################################
  !                            |             |
  !                            | VDW_ENERGY  |
  !                            |_____________|
  !
  ! This routine carries out the integration of equation 8 of SOLER.  It
  ! returns the non-local exchange-correlation energy and the u_alpha(k)
  ! arrays used to find the u_alpha(r) arrays via equations 11 and 12 in
  ! SOLER.

  SUBROUTINE vdW_energy (thetas, vdW_xc_energy)

  USE gvect,           ONLY : nl, nlm, gg, ngm, igtongl, gl, ngl, gstart
  USE cell_base,       ONLY : tpiba, omega

  implicit none

  complex(dp), intent(inout) :: thetas(:,:)            ! On input this variable holds the theta
                                                       ! functions (equation 8, SOLER) in the format
                                                       ! thetas(grid_point, theta_i). On output
                                                       ! this array holds u_alpha(k) =
                                                       ! Sum_j[theta_beta(k)phi_alpha_beta(k)]

  real(dp), intent(out) :: vdW_xc_energy               ! The non-local correlation energy.

  real(dp), allocatable :: kernel_of_k(:,:)            ! This array will hold the interpolated kernel
                                                       ! values for each pair of q values in the q_mesh.

  real(dp)    :: g                                     ! The magnitude of the current g vector.
  integer     :: last_g                                ! The shell number of the last g vector.


  integer     :: g_i, q1_i, q2_i, i_grid               ! Index variables.

  complex(dp) :: theta(Nqs), thetam(Nqs), theta_g(Nqs) ! Temporary storage arrays used since we
                                                       ! are overwriting the thetas array here.
  real(dp)    :: G0_term, G_multiplier

  complex(dp), allocatable :: u_vdw(:,:)               ! Temporary array holding u_alpha(k).




  vdW_xc_energy = 0.0D0
  allocate (u_vdW(dfftp%nnr,Nqs), kernel_of_k(Nqs, Nqs))
  u_vdW(:,:) = CMPLX(0.0_DP,0.0_DP)


  ! --------------------------------------------------------------------
  ! Loop over PWSCF's array of magnitude-sorted g-vector shells. For
  ! each shell, interpolate the kernel at this magnitude of g, then find
  ! all points on the shell and carry out the integration over those
  ! points. The PWSCF variables used here are ngm = number of g-vectors
  ! on this processor, nl = an array that gives the indices into the FFT
  ! grid for a particular g vector, igtongl = an array that gives the
  ! index of which shell a particular g vector is in, gl = an array that
  ! gives the magnitude of the g vectors for each shell. In essence, we
  ! are forming the reciprocal-space u(k) functions of SOLER equation
  ! 11. These are kept in thetas array. Here we should use gstart,ngm
  ! but all the cases are handled by conditionals inside the loop

  G_multiplier = 1.0D0
  if (gamma_only) G_multiplier = 2.0D0

  last_g = -1

  do g_i = 1, ngm

     if ( igtongl(g_i) .ne. last_g) then

        g = sqrt(gl(igtongl(g_i))) * tpiba
        call interpolate_kernel(g, kernel_of_k)
        last_g = igtongl(g_i)

     end if

     theta = thetas(nl(g_i),:)

     do q2_i = 1, Nqs
        do q1_i = 1, Nqs
           u_vdW(nl(g_i),q2_i)  = u_vdW(nl(g_i),q2_i) + kernel_of_k(q2_i,q1_i)*theta(q1_i)
        end do
        vdW_xc_energy = vdW_xc_energy + G_multiplier * (u_vdW(nl(g_i),q2_i)*conjg(theta(q2_i)))
     end do

     if (g_i < gstart ) vdW_xc_energy = vdW_xc_energy / G_multiplier

  end do

  if (gamma_only) u_vdW(nlm(:),:) = CONJG(u_vdW(nl(:),:))


  ! --------------------------------------------------------------------
  ! Apply scaling factors. The e2 comes from PWSCF's choice of units.
  ! This should be 0.5 * e2 * vdW_xc_energy * (2pi)^3/omega * (omega)^2,
  ! with the (2pi)^3/omega being the volume element for the integral
  ! (the volume of the reciprocal unit cell) and the 2 factors of omega
  ! being used to cancel the factor of 1/omega PWSCF puts on forward
  ! FFTs of the 2 theta factors. 1 omega cancels and the (2pi)^3
  ! cancels because there should be a factor of 1/(2pi)^3 on the radial
  ! Fourier transform of phi that was left out to cancel with this
  ! factor.

  vdW_xc_energy = 0.5D0 * e2 * omega * vdW_xc_energy

  thetas(:,:) = u_vdW(:,:)
  deallocate (u_vdW, kernel_of_k)

  END SUBROUTINE vdW_energy








  ! ####################################################################
  !                          |                 |
  !                          |  GET_POTENTIAL  |
  !                          |_________________|
  !
  ! This routine finds the non-local correlation contribution to the
  ! potential (i.e. the derivative of the non-local piece of the energy
  ! with respect to density) given in SOLER equation 10. The u_alpha(k)
  ! functions were found while calculating the energy. They are passed
  ! in as the matrix u_vdW. Most of the required derivatives were
  ! calculated in the "get_q0_on_grid" routine, but the derivative of
  ! the interpolation polynomials, P_alpha(q), (SOLER equation 3) with
  ! respect to q is interpolated here, along with the polynomials
  ! themselves.

  SUBROUTINE get_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, u_vdW, potential)

  USE gvect,               ONLY : nl, g, nlm
  USE cell_base,           ONLY : alat, tpiba

  implicit none

  real(dp), intent(in) ::  q0(:), grad_rho(:,:)       ! Input arrays holding the value of q0 for
                                                      ! all points assigned to this processor and
                                                      ! the gradient of the charge density for
                                                      ! points assigned to this processor.

  real(dp), intent(in) :: dq0_drho(:), dq0_dgradrho(:)! The derivative of q0 with respect to the
                                                      ! charge density and gradient of the charge
                                                      ! density (almost). See comments in the
                                                      ! get_q0_on_grid subroutine above.

  complex(dp), intent(in) :: u_vdW(:,:)               ! The functions u_alpha(r) obtained by
                                                      ! inverse transforming the functions
                                                      ! u_alph(k). See equations 11 and 12 in SOLER

  real(dp), intent(inout) :: potential(:)             ! The non-local correlation potential for
                                                      ! points on the grid over the whole cell (not
                                                      ! just those assigned to this processor).

  real(dp), allocatable, save :: d2y_dx2(:,:)         ! Second derivatives of P_alpha polynomials
                                                      ! for interpolation.

  integer :: i_grid, P_i,icar                         ! Index variables.

  integer :: q_low, q_hi, q                           ! Variables to find the bin in the q_mesh that
                                                      ! a particular q0 belongs to (for interpolation).
  real(dp) :: dq, a, b, c, d, e, f                    ! Intermediate variables used in the
                                                      ! interpolation of the polynomials.

  real(dp) :: y(Nqs), dP_dq0, P                       ! The y values for a given polynomial (all 0
                                                      ! exept for element i of P_i) The derivative
                                                      ! of P at a given q0 and the value of P at a
                                                      ! given q0. Both of these are interpolated
                                                      ! below.

1131 1132
  real(dp) :: gradient2                               ! Squared gradient.

thonhauser's avatar
thonhauser committed
1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148
  real(dp)   , allocatable ::h_prefactor(:)
  complex(dp), allocatable ::h(:)




  allocate (h_prefactor(dfftp%nnr), h(dfftp%nnr))

  potential     = 0.0D0
  h_prefactor   = 0.0D0


  ! --------------------------------------------------------------------
  ! Get the second derivatives of the P_i functions for interpolation.
  ! We have already calculated this once but it is very fast and it's
  ! just as easy to calculate it again.
1149 1150

  if (.not. allocated( d2y_dx2) ) then
thonhauser's avatar
thonhauser committed
1151

1152
     allocate( d2y_dx2(Nqs, Nqs) )
thonhauser's avatar
thonhauser committed
1153 1154
     call initialize_spline_interpolation (q_mesh, d2y_dx2(:,:))

1155 1156
  end if

thonhauser's avatar
thonhauser committed
1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210

  do i_grid = 1, dfftp%nnr

     q_low = 1
     q_hi = Nqs


     ! -----------------------------------------------------------------
     ! Figure out which bin our value of q0 is in in the q_mesh.

     do while ( (q_hi - q_low) > 1)

        q = int((q_hi + q_low)/2)

        if (q_mesh(q) > q0(i_grid)) then
           q_hi = q
        else
           q_low = q
        end if

     end do

     if (q_hi == q_low) call errore('get_potential','qhi == qlow',1)

     dq = q_mesh(q_hi) - q_mesh(q_low)

     a = (q_mesh(q_hi) - q0(i_grid))/dq
     b = (q0(i_grid) - q_mesh(q_low))/dq
     c = (a**3 - a)*dq**2/6.0D0
     d = (b**3 - b)*dq**2/6.0D0
     e = (3.0D0*a**2 - 1.0D0)*dq/6.0D0
     f = (3.0D0*b**2 - 1.0D0)*dq/6.0D0

     do P_i = 1, Nqs
        y = 0.0D0
        y(P_i) = 1.0D0

        P      = a*y(q_low) + b*y(q_hi)  + c*d2y_dx2(P_i,q_low) + d*d2y_dx2(P_i,q_hi)
        dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(P_i,q_low) + f*d2y_dx2(P_i,q_hi)


        ! --------------------------------------------------------------
        ! The first term in equation 10 of SOLER.

        potential(i_grid) = potential(i_grid) + u_vdW(i_grid,P_i)* (P + dP_dq0 * dq0_drho(i_grid))
        if (q0(i_grid) .ne. q_mesh(Nqs)) then
           h_prefactor(i_grid) = h_prefactor(i_grid) + u_vdW(i_grid,P_i)*dP_dq0*dq0_dgradrho(i_grid)
        end if
     end do
  end do

  do icar = 1,3

     h(:) = CMPLX( h_prefactor(:) * grad_rho(:,icar), 0.0_DP )
1211 1212 1213 1214

     do i_grid = 1, dfftp%nnr
        gradient2 = grad_rho(i_grid,1)**2 + grad_rho(i_grid,2)**2 + grad_rho(i_grid,3)**2
        if ( gradient2 > 0.0D0 ) h(i_grid) = h(i_grid) / SQRT( gradient2 )
1215
     end do
1216

thonhauser's avatar
thonhauser committed
1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255