mirror of https://gitlab.com/QEF/q-e.git
2043 lines
70 KiB
Fortran
2043 lines
70 KiB
Fortran
! Copyright (C) 2001-2009 Quantum ESPRESSO group
|
||
! Copyright (C) 2015 Brian Kolb, Timo Thonhauser
|
||
! 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 .
|
||
! ----------------------------------------------------------------------
|
||
|
||
|
||
MODULE vdW_DF
|
||
|
||
|
||
! 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,
|
||
! and P. Hyldgaard, Phys. Rev. Lett. 115, 136402 (2015).
|
||
!
|
||
! 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.
|
||
|
||
|
||
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
|
||
USE mp_bands, ONLY : intra_bgrp_comm
|
||
USE io_global, ONLY : stdout, ionode
|
||
USE fft_base, ONLY : dfftp
|
||
USE fft_interfaces, ONLY : fwfft, invfft
|
||
USE control_flags, ONLY : iverbosity, gamma_only
|
||
|
||
implicit none
|
||
|
||
REAL(DP), PARAMETER :: epsr =1.d-12 ! a small number to cut off densities
|
||
integer :: vdw_type = 1
|
||
|
||
private
|
||
public :: xc_vdW_DF, xc_vdW_DF_spin, stress_vdW_DF, interpolate_kernel, &
|
||
vdw_type, initialize_spline_interpolation
|
||
|
||
CONTAINS
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | functions |
|
||
! |_____________|
|
||
!
|
||
! Functions to be used in get_q0_on_grid and get_q0_on_grid_spin().
|
||
|
||
function Fs(s)
|
||
|
||
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).
|
||
|
||
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
|
||
|
||
end function Fs
|
||
|
||
|
||
|
||
|
||
function dFs_ds(s)
|
||
|
||
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).
|
||
|
||
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
|
||
|
||
end function dFs_ds
|
||
|
||
|
||
|
||
|
||
function kF(rho)
|
||
|
||
implicit none
|
||
real(dp) :: rho, kF
|
||
|
||
kF = ( 3.0D0 * pi**2 * rho )**(1.0D0/3.0D0)
|
||
|
||
end function kF
|
||
|
||
|
||
|
||
|
||
function dkF_drho(rho)
|
||
|
||
implicit none
|
||
real(dp) :: rho, dkF_drho
|
||
|
||
dkF_drho = (1.0D0/3.0D0) * kF(rho) / rho
|
||
|
||
end function dkF_drho
|
||
|
||
|
||
|
||
|
||
function ds_drho(rho, s)
|
||
|
||
implicit none
|
||
real(dp) :: rho, s, ds_drho
|
||
|
||
ds_drho = -s * ( dkF_drho(rho) / kF(rho) + 1.0D0 / rho )
|
||
|
||
end function ds_drho
|
||
|
||
|
||
|
||
|
||
function ds_dgradrho(rho)
|
||
|
||
implicit none
|
||
real(dp) :: rho, ds_dgradrho
|
||
|
||
ds_dgradrho = 1.0D0 / (2.0D0 * kF(rho) * rho)
|
||
|
||
end function ds_dgradrho
|
||
|
||
|
||
|
||
|
||
function dqx_drho(rho, s)
|
||
|
||
implicit none
|
||
real(dp) :: rho, s, dqx_drho
|
||
|
||
dqx_drho = dkF_drho(rho) * Fs(s) + kF(rho) * dFs_ds(s) * ds_drho(rho, s)
|
||
|
||
end function dqx_drho
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | XC_VDW_DF |
|
||
! |_____________|
|
||
|
||
SUBROUTINE xc_vdW_DF (rho_valence, rho_core, etxc, vtxc, v)
|
||
|
||
USE gvect, ONLY : ngm, g
|
||
USE cell_base, ONLY : omega, tpiba
|
||
|
||
implicit none
|
||
|
||
! --------------------------------------------------------------------
|
||
! Local variables
|
||
! _
|
||
real(dp), intent(IN) :: rho_valence(:,:) !
|
||
real(dp), intent(IN) :: rho_core(:) ! PWSCF input variables
|
||
real(dp), intent(inout) :: etxc, vtxc, v(:,:) !_
|
||
|
||
|
||
integer :: i_grid, theta_i, i_proc ! Indexing variables over grid points,
|
||
! theta functions, and processors.
|
||
|
||
real(dp) :: grid_cell_volume ! The volume of the unit cell per G-grid point.
|
||
|
||
real(dp), allocatable :: q0(:) ! The saturated value of q (equations 11 and 12
|
||
! of DION). This saturation is that of
|
||
! equation 5 in SOLER.
|
||
|
||
real(dp), allocatable :: grad_rho(:,:) ! The gradient of the charge density. The
|
||
! format is as follows:
|
||
! grad_rho(cartesian_component,grid_point).
|
||
|
||
real(dp), allocatable :: potential(:) ! The vdW contribution to the potential
|
||
|
||
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).
|
||
|
||
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).
|
||
|
||
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).
|
||
|
||
real(dp) :: Ec_nl ! The non-local vdW contribution to the energy.
|
||
|
||
real(dp), allocatable :: total_rho(:) ! This is the sum of the valence and core
|
||
! charge. This just holds the piece assigned
|
||
! to this processor.
|
||
|
||
logical, save :: first_iteration = .TRUE. ! Whether this is the first time this
|
||
! routine has been called.
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Check that the requested non-local functional is implemented.
|
||
|
||
if ( vdW_type /= 1 .AND. vdW_type /= 2) call errore('xc_vdW_DF','E^nl_c not implemented',1)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Write out the vdW-DF imformation.
|
||
|
||
if ( ionode .AND. first_iteration ) call vdW_info
|
||
first_iteration = .FALSE.
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Allocate arrays. nnr is a PWSCF variable that holds the number of
|
||
! points assigned to a given processor.
|
||
|
||
allocate( q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr), &
|
||
grad_rho(3,dfftp%nnr) )
|
||
allocate( total_rho(dfftp%nnr), potential(dfftp%nnr), thetas(dfftp%nnr, Nqs) )
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
total_rho = rho_valence(:,1) + rho_core(:)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Here we calculate the gradient in reciprocal space using FFT.
|
||
|
||
call fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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
|
||
call mp_sum(Ec_nl, intra_bgrp_comm)
|
||
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
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
do theta_i = 1, Nqs
|
||
CALL invfft('Rho', thetas(:,theta_i), dfftp)
|
||
end do
|
||
|
||
call get_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, thetas, potential)
|
||
v(:,1) = v(:,1) + e2 * potential(:)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! The integral of rho(r)*potential(r) for the vtxc output variable.
|
||
|
||
grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
vtxc = vtxc + e2 * grid_cell_volume * rho_valence(i_grid,1) * potential(i_grid)
|
||
end do
|
||
|
||
deallocate ( potential, q0, grad_rho, dq0_drho, dq0_dgradrho, total_rho, thetas )
|
||
|
||
END SUBROUTINE xc_vdW_DF
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | XC_VDW_DF_spin |
|
||
! |__________________|
|
||
!
|
||
! This subroutine is as similar to xc_vdW_DF() as possible, but
|
||
! handles the collinear nspin=2 case.
|
||
|
||
SUBROUTINE xc_vdW_DF_spin (rho_valence, rho_core, etxc, vtxc, v)
|
||
|
||
USE gvect, ONLY : ngm, g
|
||
USE cell_base, ONLY : omega, tpiba
|
||
|
||
implicit none
|
||
|
||
! --------------------------------------------------------------------
|
||
! Local variables
|
||
! _
|
||
real(dp), intent(IN) :: rho_valence(:,:) !
|
||
real(dp), intent(IN) :: rho_core(:) ! PWSCF input variables
|
||
real(dp), intent(inout) :: etxc, vtxc, v(:,:) !_
|
||
|
||
|
||
integer :: i_grid, theta_i, i_proc ! Indexing variables over grid points,
|
||
! theta functions, and processors, and a
|
||
! generic index.
|
||
|
||
real(dp) :: grid_cell_volume ! The volume of the unit cell per G-grid point.
|
||
|
||
real(dp), allocatable :: q0(:) ! The saturated value of q (equations 11 and 12
|
||
! of DION). This saturation is that of
|
||
! equation 5 in SOLER.
|
||
|
||
real(dp), allocatable :: grad_rho(:,:) ! The gradient of the charge density. The
|
||
! format is as follows:
|
||
! grad_rho(cartesian_component,grid_point)
|
||
real(dp), allocatable :: grad_rho_up(:,:) ! The gradient of the up charge density.
|
||
! Same format as grad_rho
|
||
real(dp), allocatable :: grad_rho_down(:,:) ! The gradient of the down charge density.
|
||
! Same format as grad_rho
|
||
|
||
real(dp), allocatable :: potential_up(:) ! The vdW contribution to the potential.
|
||
real(dp), allocatable :: potential_down(:) ! The vdW contribution to the potential.
|
||
|
||
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).
|
||
|
||
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).
|
||
|
||
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).
|
||
|
||
real(dp) :: Ec_nl ! The non-local vdW contribution to the energy.
|
||
|
||
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.
|
||
|
||
logical, save :: first_iteration = .TRUE. ! Whether this is the first time this
|
||
! routine has been called.
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Check that the requested non-local functional is implemented.
|
||
|
||
if ( vdW_type /= 1 .AND. vdW_type /= 2) call errore('xc_vdW_DF','E^nl_c not implemented',1)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Write out the vdW-DF imformation.
|
||
|
||
if ( ionode .AND. first_iteration ) call vdW_info
|
||
first_iteration = .FALSE.
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Allocate arrays. nnr is a PWSCF variable that holds the number of
|
||
! points assigned to a given processor.
|
||
|
||
allocate( q0(dfftp%nnr), total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr) )
|
||
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(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr) )
|
||
allocate( potential_up(dfftp%nnr), potential_down(dfftp%nnr) )
|
||
allocate( thetas(dfftp%nnr, Nqs) )
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
rho_up = rho_valence(:,1) + 0.5D0*rho_core(:)
|
||
rho_down = rho_valence(:,2) + 0.5D0*rho_core(:)
|
||
total_rho = rho_up + rho_down
|
||
|
||
#if defined (__SPIN_BALANCED)
|
||
rho_up = total_rho*0.5D0
|
||
rho_down = rho_up
|
||
write(stdout,'(/,/," Performing spin-balanced Ecnl calculation!")')
|
||
#endif
|
||
|
||
! --------------------------------------------------------------------
|
||
! Here we calculate the gradient in reciprocal space using FFT.
|
||
|
||
call fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
|
||
call fft_gradient_r2r (dfftp, rho_up, g, grad_rho_up)
|
||
call fft_gradient_r2r (dfftp, rho_down, g, grad_rho_down)
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
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)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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
|
||
call mp_sum(Ec_nl, intra_bgrp_comm)
|
||
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
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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.
|
||
|
||
do theta_i = 1, Nqs
|
||
CALL invfft('Rho', thetas(:,theta_i), dfftp)
|
||
end do
|
||
|
||
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(:)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! The integral of rho(r)*potential(r) for the vtxc output variable
|
||
|
||
grid_cell_volume = omega/(dfftp%nr1*dfftp%nr2*dfftp%nr3)
|
||
|
||
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
|
||
|
||
|
||
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 )
|
||
|
||
END SUBROUTINE xc_vdW_DF_spin
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | 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|
|
||
|
||
SUBROUTINE get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(IN) :: total_rho(:), grad_rho(:,:) ! Input variables needed
|
||
|
||
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.
|
||
|
||
integer, parameter :: m_cut = 12 ! How many terms to include in the sum
|
||
! of SOLER equation 5.
|
||
|
||
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.
|
||
|
||
integer :: i_grid, idx ! Indexing variables.
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Initialize q0-related arrays.
|
||
|
||
q0(:) = q_cut
|
||
dq0_drho(:) = 0.0D0
|
||
dq0_dgradrho(:) = 0.0D0
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
|
||
rho = total_rho(i_grid)
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
if ( rho < epsr ) cycle
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! Calculate some intermediate values needed to find q.
|
||
|
||
r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
|
||
|
||
s = sqrt( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 ) &
|
||
/ (2.0D0 * kF(rho) * rho )
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
call pw(r_s, 1, ec, dq0_drho(i_grid))
|
||
q = -4.0D0*pi/3.0D0 * ec + kF(rho) * Fs(s)
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! Bring q into its proper bounds.
|
||
|
||
CALL saturate_q ( q, q_cut, q0(i_grid), dq0_dq )
|
||
if (q0(i_grid) < q_min) q0(i_grid) = q_min
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
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)
|
||
|
||
end do
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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)
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
|
||
end do
|
||
|
||
do idx = 1, Nqs
|
||
CALL fwfft ('Rho', thetas(:,idx), dfftp)
|
||
end do
|
||
|
||
END SUBROUTINE get_q0_on_grid
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | GET_Q0_ON_GRID_spin |
|
||
! |_______________________|
|
||
|
||
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)
|
||
|
||
implicit none
|
||
|
||
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.
|
||
|
||
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.
|
||
|
||
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
|
||
|
||
|
||
|
||
|
||
fac = 2.0D0**(-1.0D0/3.0D0)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Initialize q0-related arrays.
|
||
|
||
q0(:) = q_cut
|
||
dq0_drho_up(:) = 0.0D0
|
||
dq0_drho_down(:) = 0.0D0
|
||
dq0_dgradrho_up(:) = 0.0D0
|
||
dq0_dgradrho_down(:) = 0.0D0
|
||
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
|
||
rho = total_rho(i_grid)
|
||
up = rho_up(i_grid)
|
||
down = rho_down(i_grid)
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
if ( rho < epsr ) cycle
|
||
|
||
calc_qx_up = .TRUE.
|
||
calc_qx_down = .TRUE.
|
||
|
||
if ( up < epsr/2.0D0 ) calc_qx_up = .FALSE.
|
||
if ( down < epsr/2.0D0 ) calc_qx_down = .FALSE.
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
q0x_up = 0.0D0
|
||
q0x_down = 0.0D0
|
||
dqx_drho_up = 0.0D0
|
||
dqx_drho_down = 0.0D0
|
||
|
||
|
||
if (calc_qx_up) then
|
||
s_up = sqrt( grad_rho_up(1,i_grid)**2 + grad_rho_up(2,i_grid)**2 + &
|
||
grad_rho_up(3,i_grid)**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
|
||
|
||
if (calc_qx_down) then
|
||
s_down = sqrt( grad_rho_down(1,i_grid)**2 + grad_rho_down(2,i_grid)**2 + &
|
||
grad_rho_down(3,i_grid)**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
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! This is the q value defined in equations 11 and 12 of DION and
|
||
! equation 8 of THONHAUSER (also see text above that equation).
|
||
|
||
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)
|
||
|
||
qx = ( up * q0x_up + down * q0x_down ) / rho
|
||
qc = -4.0D0*pi/3.0D0 * ec
|
||
q = qx + qc
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! Bring q into its proper bounds.
|
||
|
||
CALL saturate_q (q, q_cut, q0(i_grid), dq0_dq)
|
||
if (q0(i_grid) < q_min) q0(i_grid) = q_min
|
||
|
||
|
||
! -----------------------------------------------------------------
|
||
! 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.
|
||
|
||
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
|
||
|
||
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
|
||
|
||
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
|
||
|
||
dqc_drho_up = -4.0D0*pi/3.0D0 * (dqc_drho_up - ec)
|
||
dqc_drho_down = -4.0D0*pi/3.0D0 * (dqc_drho_down - ec)
|
||
|
||
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)
|
||
|
||
end do
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! 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)
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
thetas(i_grid,:) = thetas(i_grid,:) * total_rho(i_grid)
|
||
end do
|
||
|
||
do idx = 1, Nqs
|
||
CALL fwfft ('Rho', thetas(:,idx), dfftp)
|
||
end do
|
||
|
||
END SUBROUTINE get_q0_on_grid_spin
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | saturate_q |
|
||
! |______________|
|
||
|
||
SUBROUTINE saturate_q (q, q_cut, q0, dq0_dq)
|
||
|
||
implicit none
|
||
|
||
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 : 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,kind=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(dfftp%nl(g_i),:)
|
||
|
||
do q2_i = 1, Nqs
|
||
do q1_i = 1, Nqs
|
||
u_vdW(dfftp%nl(g_i),q2_i) = u_vdW(dfftp%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(dfftp%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(dfftp%nlm(:),:) = CONJG(u_vdW(dfftp%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 : g
|
||
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.
|
||
|
||
real(dp) :: gradient2 ! Squared gradient.
|
||
|
||
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.
|
||
|
||
if (.not. allocated( d2y_dx2) ) then
|
||
|
||
allocate( d2y_dx2(Nqs, Nqs) )
|
||
call initialize_spline_interpolation (q_mesh, d2y_dx2(:,:))
|
||
|
||
end if
|
||
|
||
|
||
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, kind=dp )
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
gradient2 = grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2
|
||
if ( gradient2 > 0.0D0 ) h(i_grid) = h(i_grid) / SQRT( gradient2 )
|
||
end do
|
||
|
||
CALL fwfft ('Rho', h, dfftp)
|
||
h(dfftp%nl(:)) = CMPLX(0.0_DP,1.0_DP, kind=dp) * tpiba * g(icar,:) * h(dfftp%nl(:))
|
||
if (gamma_only) h(dfftp%nlm(:)) = CONJG(h(dfftp%nl(:)))
|
||
CALL invfft ('Rho', h, dfftp)
|
||
potential(:) = potential(:) - REAL(h(:))
|
||
|
||
end do
|
||
|
||
deallocate (h_prefactor, h)
|
||
|
||
END SUBROUTINE get_potential
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | SPLINE_INTERPOLATION |
|
||
! |________________________|
|
||
!
|
||
! This routine is modeled after an algorithm from "Numerical Recipes
|
||
! in C" by Cambridge University press, page 97. It was adapted for
|
||
! Fortran, of course and for the problem at hand, in that it finds the
|
||
! bin a particular x value is in and then loops over all the P_i
|
||
! functions so we only have to find the bin once.
|
||
|
||
SUBROUTINE spline_interpolation (x, evaluation_points, values)
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(in) :: x(:), evaluation_points(:) ! Input variables. The x values used to
|
||
! form the interpolation (q_mesh in this
|
||
! case) and the values of q0 for which we
|
||
! are interpolating the function.
|
||
|
||
complex(dp), intent(inout) :: values(:,:) ! An output array (allocated outside this
|
||
! routine) that stores the interpolated
|
||
! values of the P_i (SOLER equation 3)
|
||
! polynomials. The format is
|
||
! values(grid_point, P_i).
|
||
|
||
integer :: Ngrid_points, Nx ! Total number of grid points to evaluate
|
||
! and input x points.
|
||
|
||
real(dp), allocatable, save :: d2y_dx2(:,:) ! The second derivatives required to do
|
||
! the interpolation.
|
||
|
||
integer :: i_grid, lower_bound, upper_bound, idx, P_i ! Some indexing variables.
|
||
|
||
real(dp), allocatable :: y(:) ! Temporary variables needed for the
|
||
real(dp) :: a, b, c, d, dx ! interpolation.
|
||
|
||
|
||
|
||
|
||
Nx = size(x)
|
||
Ngrid_points = size(evaluation_points)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Allocate the temporary array.
|
||
|
||
allocate( y(Nx) )
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! If this is the first time this routine has been called we need to
|
||
! get the second derivatives (d2y_dx2) required to perform the
|
||
! interpolations. So we allocate the array and call
|
||
! initialize_spline_interpolation to get d2y_dx2.
|
||
|
||
if (.not. allocated(d2y_dx2) ) then
|
||
|
||
allocate( d2y_dx2(Nx,Nx) )
|
||
call initialize_spline_interpolation(x, d2y_dx2)
|
||
|
||
end if
|
||
|
||
do i_grid=1, Ngrid_points
|
||
|
||
lower_bound = 1
|
||
upper_bound = Nx
|
||
|
||
do while ( (upper_bound - lower_bound) > 1 )
|
||
|
||
idx = (upper_bound+lower_bound) / 2
|
||
|
||
if ( evaluation_points(i_grid) > x(idx) ) then
|
||
lower_bound = idx
|
||
else
|
||
upper_bound = idx
|
||
end if
|
||
|
||
end do
|
||
|
||
dx = x(upper_bound)-x(lower_bound)
|
||
|
||
a = (x(upper_bound) - evaluation_points(i_grid))/dx
|
||
b = (evaluation_points(i_grid) - x(lower_bound))/dx
|
||
c = ((a**3-a)*dx**2)/6.0D0
|
||
d = ((b**3-b)*dx**2)/6.0D0
|
||
|
||
do P_i = 1, Nx
|
||
|
||
y = 0
|
||
y(P_i) = 1
|
||
|
||
values(i_grid, P_i) = a*y(lower_bound) + b*y(upper_bound) &
|
||
+ (c*d2y_dx2(P_i,lower_bound) + d*d2y_dx2(P_i, upper_bound))
|
||
end do
|
||
|
||
end do
|
||
|
||
deallocate( y )
|
||
|
||
END SUBROUTINE spline_interpolation
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | INITIALIZE_SPLINE_INTERPOLATION |
|
||
! |___________________________________|
|
||
!
|
||
! This routine is modeled after an algorithm from "Numerical Recipes
|
||
! in C" by Cambridge University Press, pages 96-97. It was adapted
|
||
! for Fortran and for the problem at hand.
|
||
|
||
SUBROUTINE initialize_spline_interpolation (x, d2y_dx2)
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(in) :: x(:) ! The input abscissa values.
|
||
real(dp), intent(inout) :: d2y_dx2(:,:) ! The output array (allocated outside this routine)
|
||
! that holds the second derivatives required for
|
||
! interpolating the function.
|
||
|
||
integer :: Nx, P_i, idx ! The total number of x points and some indexing
|
||
! variables.
|
||
|
||
real(dp), allocatable :: temp_array(:), y(:) ! Some temporary arrays required. y is the array
|
||
! that holds the funcion values (all either 0 or
|
||
! 1 here).
|
||
|
||
real(dp) :: temp1, temp2 ! Some temporary variables required.
|
||
|
||
|
||
|
||
|
||
Nx = size(x)
|
||
|
||
allocate( temp_array(Nx), y(Nx) )
|
||
|
||
do P_i=1, Nx
|
||
|
||
! -----------------------------------------------------------------
|
||
! In the Soler method, the polynomicals that are interpolated are Kroneker
|
||
! delta funcions at a particular q point. So, we set all y values to 0
|
||
! except the one corresponding to the particular function P_i.
|
||
|
||
y = 0.0D0
|
||
y(P_i) = 1.0D0
|
||
|
||
d2y_dx2(P_i,1) = 0.0D0
|
||
temp_array(1) = 0.0D0
|
||
|
||
do idx = 2, Nx-1
|
||
|
||
temp1 = (x(idx)-x(idx-1))/(x(idx+1)-x(idx-1))
|
||
temp2 = temp1 * d2y_dx2(P_i,idx-1) + 2.0D0
|
||
d2y_dx2(P_i,idx) = (temp1-1.0D0)/temp2
|
||
|
||
temp_array(idx) = (y(idx+1)-y(idx))/(x(idx+1)-x(idx)) &
|
||
- (y(idx)-y(idx-1))/(x(idx)-x(idx-1))
|
||
temp_array(idx) = (6.0D0*temp_array(idx)/(x(idx+1)-x(idx-1)) &
|
||
- temp1*temp_array(idx-1))/temp2
|
||
|
||
end do
|
||
|
||
d2y_dx2(P_i,Nx) = 0.0D0
|
||
|
||
do idx=Nx-1, 1, -1
|
||
|
||
d2y_dx2(P_i,idx) = d2y_dx2(P_i,idx) * d2y_dx2(P_i,idx+1) + temp_array(idx)
|
||
|
||
end do
|
||
|
||
end do
|
||
|
||
deallocate( temp_array, y)
|
||
|
||
END SUBROUTINE initialize_spline_interpolation
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | INTERPOLATE_KERNEL |
|
||
! |____________________|
|
||
!
|
||
! This routine is modeled after an algorithm from "Numerical Recipes in C" by
|
||
! Cambridge University Press, page 97. Adapted for Fortran and the problem at
|
||
! hand. This function is used to find the Phi_alpha_beta needed for equations
|
||
! 8 and 11 of SOLER.
|
||
|
||
SUBROUTINE interpolate_kernel (k, kernel_of_k)
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(in) :: k ! Input value, the magnitude of the g-vector
|
||
! for the current point.
|
||
|
||
real(dp), intent(inout) :: kernel_of_k(:,:) ! An output array (allocated outside this routine)
|
||
! that holds the interpolated value of the kernel
|
||
! for each pair of q points (i.e. the phi_alpha_beta
|
||
! of the Soler method.
|
||
|
||
integer :: q1_i, q2_i, k_i ! Indexing variables.
|
||
|
||
real(dp) :: A, B, C, D ! Intermediate values for the interpolation.
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Check to make sure that the kernel table we have is capable of
|
||
! dealing with this value of k. If k is larger than
|
||
! Nr_points*2*pi/r_max then we can't perform the interpolation. In
|
||
! that case, a kernel file should be generated with a larger number of
|
||
! radial points.
|
||
|
||
if ( k >= Nr_points*dk ) then
|
||
|
||
write(*,'(A,F10.5,A,F10.5)') "k = ", k, " k_max = ",Nr_points*dk
|
||
call errore('interpolate kernel', 'k value requested is out of range',1)
|
||
|
||
end if
|
||
|
||
kernel_of_k = 0.0D0
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! This integer division figures out which bin k is in since the kernel
|
||
! is set on a uniform grid.
|
||
|
||
k_i = int(k/dk)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Test to see if we are trying to interpolate a k that is one of the
|
||
! actual function points we have. The value is just the value of the
|
||
! function in that case.
|
||
|
||
if (mod(k,dk) == 0) then
|
||
|
||
do q1_i = 1, Nqs
|
||
do q2_i = 1, q1_i
|
||
|
||
kernel_of_k(q1_i, q2_i) = kernel(k_i,q1_i, q2_i)
|
||
kernel_of_k(q2_i, q1_i) = kernel(k_i,q2_i, q1_i)
|
||
|
||
end do
|
||
end do
|
||
|
||
return
|
||
|
||
end if
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! If we are not on a function point then we carry out the
|
||
! interpolation.
|
||
|
||
A = (dk*(k_i+1.0D0) - k)/dk
|
||
B = (k - dk*k_i)/dk
|
||
C = (A**3-A)*dk**2/6.0D0
|
||
D = (B**3-B)*dk**2/6.0D0
|
||
|
||
do q1_i = 1, Nqs
|
||
do q2_i = 1, q1_i
|
||
|
||
kernel_of_k(q1_i, q2_i) = A*kernel(k_i, q1_i, q2_i) + B*kernel(k_i+1, q1_i, q2_i) &
|
||
+(C*d2phi_dk2(k_i, q1_i, q2_i) + D*d2phi_dk2(k_i+1, q1_i, q2_i))
|
||
|
||
kernel_of_k(q2_i, q1_i) = kernel_of_k(q1_i, q2_i)
|
||
|
||
end do
|
||
end do
|
||
|
||
END SUBROUTINE interpolate_kernel
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | STRESS_VDW_DF |
|
||
! |_________________|
|
||
|
||
SUBROUTINE stress_vdW_DF (rho_valence, rho_core, nspin, sigma)
|
||
|
||
use gvect, ONLY : ngm, g
|
||
USE cell_base, ONLY : tpiba
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(IN) :: rho_valence(:,:) !
|
||
real(dp), intent(IN) :: rho_core(:) ! Input variables
|
||
integer, intent(IN) :: nspin !
|
||
real(dp), intent(inout) :: sigma(3,3) !
|
||
|
||
real(dp), allocatable :: grad_rho(:,:) !
|
||
real(dp), allocatable :: total_rho(:) ! Rho values
|
||
|
||
real(dp), allocatable :: q0(:) !
|
||
real(dp), allocatable :: dq0_drho(:) ! Q-values
|
||
real(dp), allocatable :: dq0_dgradrho(:) !
|
||
|
||
complex(dp), allocatable :: thetas(:,:) ! Thetas
|
||
integer :: i_proc, theta_i, l, m
|
||
|
||
real(dp) :: sigma_grad(3,3)
|
||
real(dp) :: sigma_ker(3,3)
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Tests
|
||
|
||
#if defined (__SPIN_BALANCED)
|
||
if ( nspin==2 ) then
|
||
write(stdout,'(/,/ " Performing spin-balanced Ecnl stress calculation!")')
|
||
else if ( nspin > 2 ) then
|
||
call errore ('stres_vdW_DF','noncollinear vdW stress not implemented',1)
|
||
end if
|
||
#else
|
||
if ( nspin>=2 ) then
|
||
call errore ('stres_vdW_DF', 'vdW stress not implemented for nspin > 1',1)
|
||
end if
|
||
#endif
|
||
|
||
sigma(:,:) = 0.0_DP
|
||
sigma_grad(:,:) = 0.0_DP
|
||
sigma_ker(:,:) = 0.0_DP
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Allocations
|
||
|
||
allocate( grad_rho(3,dfftp%nnr) )
|
||
allocate( total_rho(dfftp%nnr) )
|
||
allocate( q0(dfftp%nnr) )
|
||
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
|
||
allocate( thetas(dfftp%nnr, Nqs) )
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Charge
|
||
|
||
total_rho = rho_valence(:,1) + rho_core(:)
|
||
#if defined (__SPIN_BALANCED)
|
||
if ( nspin == 2 ) then
|
||
total_rho = rho_valence(:,2) + total_rho(:)
|
||
end if
|
||
#endif
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Here we calculate the gradient in reciprocal space using FFT.
|
||
|
||
call fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
|
||
|
||
! --------------------------------------------------------------------
|
||
! Get q0.
|
||
|
||
CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Stress
|
||
|
||
CALL stress_vdW_DF_gradient (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas, sigma_grad)
|
||
CALL stress_vdW_DF_kernel (total_rho, q0, thetas, sigma_ker)
|
||
|
||
sigma = - (sigma_grad + sigma_ker)
|
||
|
||
do l = 1, 3
|
||
do m = 1, l - 1
|
||
sigma (m, l) = sigma (l, m)
|
||
enddo
|
||
enddo
|
||
|
||
deallocate( grad_rho, total_rho, q0, dq0_drho, dq0_dgradrho, thetas )
|
||
|
||
END SUBROUTINE stress_vdW_DF
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | STRESS_VDW_DF_GRADIENT |
|
||
! |__________________________|
|
||
|
||
SUBROUTINE stress_vdW_DF_gradient (total_rho, grad_rho, q0, dq0_drho, &
|
||
dq0_dgradrho, thetas, sigma)
|
||
|
||
USE gvect, ONLY : ngm, g, gg, igtongl, &
|
||
gl, ngl, gstart
|
||
USE cell_base, ONLY : omega, tpiba, alat, at, tpiba2
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(IN) :: total_rho(:) !
|
||
real(dp), intent(IN) :: grad_rho(:, :) ! Input variables.
|
||
real(dp), intent(inout) :: sigma(:,:) !
|
||
real(dp), intent(IN) :: q0(:) !
|
||
real(dp), intent(IN) :: dq0_drho(:) !
|
||
real(dp), intent(IN) :: dq0_dgradrho(:) !
|
||
complex(dp), intent(IN) :: thetas(:,:) !
|
||
|
||
complex(dp), allocatable :: u_vdW(:,:) !
|
||
|
||
real(dp), allocatable :: d2y_dx2(:,:) !
|
||
real(dp) :: y(Nqs), dP_dq0, P, a, b, c, d, e, f ! Interpolation.
|
||
real(dp) :: dq !
|
||
|
||
integer :: q_low, q_hi, q, q1_i, q2_i , g_i ! Loop and q-points.
|
||
|
||
integer :: l, m
|
||
real(dp) :: prefactor ! Final summation of sigma.
|
||
real(dp) :: grad2 ! Magnitude of density gradient.
|
||
|
||
integer :: i_proc, theta_i, i_grid, q_i, & !
|
||
ix, iy, iz ! Iterators.
|
||
|
||
character(LEN=1) :: intvar
|
||
|
||
|
||
|
||
|
||
allocate( d2y_dx2(Nqs, Nqs) )
|
||
allocate( u_vdW(dfftp%nnr, Nqs) )
|
||
|
||
sigma(:,:) = 0.0_DP
|
||
prefactor = 0.0_DP
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Get u in k-space.
|
||
|
||
call thetas_to_uk(thetas, u_vdW)
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Get u in real space.
|
||
|
||
do theta_i = 1, Nqs
|
||
CALL invfft('Rho', u_vdW(:,theta_i), dfftp)
|
||
end do
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Get the second derivatives for interpolating the P_i.
|
||
|
||
call initialize_spline_interpolation(q_mesh, d2y_dx2(:,:))
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Do the real space integration to obtain the stress component.
|
||
|
||
do i_grid = 1, dfftp%nnr
|
||
|
||
if ( total_rho(i_grid) < epsr ) cycle
|
||
|
||
q_low = 1
|
||
q_hi = Nqs
|
||
grad2 = sqrt( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 )
|
||
|
||
if ( grad2 == 0.0_dp ) cycle
|
||
! -----------------------------------------------------------------
|
||
! Figure out which bin our value of q0 is 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('stress_vdW_gradient','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 q_i = 1, Nqs
|
||
|
||
y(:) = 0.0D0
|
||
y(q_i) = 1.0D0
|
||
|
||
dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(q_i,q_low) + f*d2y_dx2(q_i,q_hi)
|
||
|
||
prefactor = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho(i_grid) / grad2
|
||
|
||
do l = 1, 3
|
||
do m = 1, l
|
||
|
||
sigma (l, m) = sigma (l, m) - e2 * prefactor * &
|
||
(grad_rho(l,i_grid) * grad_rho(m,i_grid))
|
||
end do
|
||
end do
|
||
|
||
end do
|
||
|
||
end do
|
||
call mp_sum( sigma, intra_bgrp_comm )
|
||
|
||
call dscal (9, 1.d0 / (dfftp%nr1 * dfftp%nr2 * dfftp%nr3), sigma, 1)
|
||
|
||
deallocate( d2y_dx2, u_vdW )
|
||
|
||
END SUBROUTINE stress_vdW_DF_gradient
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | STRESS_VDW_DF_KERNEL |
|
||
! |__________________________|
|
||
|
||
SUBROUTINE stress_vdW_DF_kernel (total_rho, q0, thetas, sigma)
|
||
|
||
USE gvect, ONLY : ngm, g, gg, igtongl, gl, ngl, gstart
|
||
USE cell_base, ONLY : omega, tpiba, tpiba2
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(IN) :: q0(:)
|
||
real(dp), intent(IN) :: total_rho(:)
|
||
real(dp), intent(inout) :: sigma(3,3)
|
||
complex(dp), intent(IN) :: thetas(:,:)
|
||
|
||
real(dp), allocatable :: dkernel_of_dk(:,:)
|
||
|
||
integer :: l, m, q1_i, q2_i , g_i
|
||
real(dp) :: g2, ngmod2, g_kernel, G_multiplier
|
||
integer :: last_g, theta_i
|
||
|
||
|
||
|
||
|
||
allocate( dkernel_of_dk(Nqs, Nqs) )
|
||
|
||
sigma(:,:) = 0.0_DP
|
||
|
||
! --------------------------------------------------------------------
|
||
! Integration in g-space.
|
||
|
||
last_g = -1
|
||
|
||
G_multiplier = 1.0D0
|
||
|
||
if (gamma_only) G_multiplier = 2.0D0
|
||
|
||
do g_i = gstart, ngm
|
||
|
||
g2 = gg (g_i) * tpiba2
|
||
g_kernel = sqrt(g2)
|
||
|
||
if ( igtongl(g_i) .ne. last_g) then
|
||
|
||
call interpolate_Dkernel_Dk(g_kernel, dkernel_of_dk) ! Gets the derivatives.
|
||
last_g = igtongl(g_i)
|
||
|
||
end if
|
||
|
||
do q2_i = 1, Nqs
|
||
do q1_i = 1, Nqs
|
||
do l = 1, 3
|
||
do m = 1, l
|
||
|
||
sigma (l, m) = sigma (l, m) - G_multiplier * 0.5 * e2 * thetas(dfftp%nl(g_i),q1_i) * &
|
||
dkernel_of_dk(q1_i,q2_i)*conjg(thetas(dfftp%nl(g_i),q2_i))* &
|
||
(g (l, g_i) * g (m, g_i) * tpiba2) / g_kernel
|
||
end do
|
||
end do
|
||
enddo
|
||
end do
|
||
|
||
if (g_i < gstart ) sigma(:,:) = sigma(:,:) / G_multiplier
|
||
|
||
enddo
|
||
call mp_sum( sigma, intra_bgrp_comm )
|
||
|
||
deallocate( dkernel_of_dk )
|
||
|
||
END SUBROUTINE stress_vdW_DF_kernel
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | INTERPOLATE_DKERNEL_DK |
|
||
! |________________________|
|
||
|
||
|
||
SUBROUTINE interpolate_Dkernel_Dk (k, dkernel_of_dk)
|
||
|
||
implicit none
|
||
|
||
real(dp), intent(in) :: k ! Input value, the magnitude of the g-vector
|
||
! for the current point.
|
||
|
||
real(dp), intent(inout) :: dkernel_of_dk(Nqs,Nqs) ! An output array (allocated outside this
|
||
! routine) that holds the interpolated value of
|
||
! the kernel for each pair of q points (i.e. the
|
||
! phi_alpha_beta of the Soler method.
|
||
|
||
integer :: q1_i, q2_i, k_i ! Indexing variables.
|
||
|
||
real(dp) :: A, B, dAdk, dBdk, dCdk, dDdk ! Intermediate values for the interpolation.
|
||
|
||
|
||
|
||
|
||
if ( k >= Nr_points*dk ) then
|
||
|
||
write(*,'(A,F10.5,A,F10.5)') "k = ", k, " k_max = ",Nr_points*dk
|
||
call errore('interpolate kernel', 'k value requested is out of range',1)
|
||
|
||
end if
|
||
|
||
dkernel_of_dk = 0.0D0
|
||
|
||
k_i = int(k/dk)
|
||
|
||
A = (dk*(k_i+1.0D0) - k)/dk
|
||
B = (k - dk*k_i)/dk
|
||
|
||
dAdk = -1.0D0/dk
|
||
dBdk = 1.0D0/dk
|
||
dCdk = -((3*A**2 -1.0D0)/6.0D0)*dk
|
||
dDdk = ((3*B**2 -1.0D0)/6.0D0)*dk
|
||
|
||
do q1_i = 1, Nqs
|
||
do q2_i = 1, q1_i
|
||
|
||
dkernel_of_dk(q1_i, q2_i) = dAdk*kernel(k_i, q1_i, q2_i) + dBdk*kernel(k_i+1, q1_i, q2_i) &
|
||
+ dCdk*d2phi_dk2(k_i, q1_i, q2_i) + dDdk*d2phi_dk2(k_i+1, q1_i, q2_i)
|
||
|
||
dkernel_of_dk(q2_i, q1_i) = dkernel_of_dk(q1_i, q2_i)
|
||
|
||
end do
|
||
end do
|
||
|
||
END SUBROUTINE interpolate_Dkernel_Dk
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | thetas_to_uk |
|
||
! |______________|
|
||
|
||
|
||
SUBROUTINE thetas_to_uk (thetas, u_vdW)
|
||
|
||
USE gvect, ONLY : gg, ngm, igtongl, gl, ngl, gstart
|
||
USE cell_base, ONLY : tpiba, omega
|
||
|
||
implicit none
|
||
|
||
complex(dp), intent(in) :: thetas(:,:) ! On input this variable holds the theta functions
|
||
! (equation 8, SOLER) in the format
|
||
! thetas(grid_point, theta_i).
|
||
complex(dp), intent(out) :: u_vdW(:,:) ! On output this array holds u_alpha(k) =
|
||
! Sum_j[theta_beta(k)phi_alpha_beta(k)].
|
||
|
||
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
|
||
integer :: last_g, g_i, q1_i, q2_i, i_grid ! Index variables.
|
||
|
||
complex(dp) :: theta(Nqs) ! Temporary storage vector used since we are
|
||
! overwriting the thetas array here.
|
||
|
||
|
||
|
||
|
||
allocate( kernel_of_k(Nqs, Nqs) )
|
||
|
||
u_vdW(:,:) = CMPLX(0.0_DP,0.0_DP,kind=dp)
|
||
|
||
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(dfftp%nl(g_i),:)
|
||
|
||
do q2_i = 1, Nqs
|
||
do q1_i = 1, Nqs
|
||
u_vdW(dfftp%nl(g_i),q2_i) = u_vdW(dfftp%nl(g_i),q2_i) + kernel_of_k(q2_i,q1_i)*theta(q1_i)
|
||
end do
|
||
end do
|
||
|
||
end do
|
||
|
||
if (gamma_only) u_vdW(dfftp%nlm(:),:) = CONJG(u_vdW(dfftp%nl(:),:))
|
||
|
||
deallocate( kernel_of_k )
|
||
|
||
END SUBROUTINE thetas_to_uk
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
! ####################################################################
|
||
! | |
|
||
! | VDW_INFO |
|
||
! |____________|
|
||
|
||
SUBROUTINE vdW_info
|
||
|
||
implicit none
|
||
|
||
integer :: I
|
||
|
||
|
||
|
||
|
||
! --------------------------------------------------------------------
|
||
! Here we output some of the parameters being used in the run. This is
|
||
! important because these parameters are read from the
|
||
! vdW_kernel_table file. The user should ensure that these are the
|
||
! parameters they were intending to use on each run.
|
||
|
||
WRITE(stdout,'(/)')
|
||
WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"% You are using vdW-DF, which was implemented by the Thonhauser group. %")')
|
||
WRITE(stdout,'(5x,"% Please cite the following two papers that made this development %")')
|
||
WRITE(stdout,'(5x,"% possible and the two reviews that describe the various versions: %")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"% T. Thonhauser et al., PRL 115, 136402 (2015). %")')
|
||
WRITE(stdout,'(5x,"% T. Thonhauser et al., PRB 76, 125112 (2007). %")')
|
||
WRITE(stdout,'(5x,"% K. Berland et al., Rep. Prog. Phys. 78, 066501 (2015). %")')
|
||
WRITE(stdout,'(5x,"% D.C. Langreth et al., J. Phys.: Condens. Matter 21, 084203 (2009). %")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"% If you are calculating the stress with vdW-DF, please also cite: %")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"% R. Sabatini et al., J. Phys.: Condens. Matter 24, 424209 (2012). %")')
|
||
WRITE(stdout,'(5x,"% %")')
|
||
WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")')
|
||
WRITE(stdout,'(/)')
|
||
|
||
if (iverbosity > 0) then
|
||
WRITE(stdout,'(5x,"Carrying out vdW-DF run using the following parameters:")')
|
||
WRITE(stdout,'(5X,A,I3,A,I5,A,F8.3)' ) "Nqs = ",Nqs," Npoints = ",Nr_points," r_max = ",r_max
|
||
WRITE(stdout,'(5X,"q_mesh =",4F12.8)') (q_mesh(I), I=1, 4)
|
||
WRITE(stdout,'(13X,4F12.8)') (q_mesh(I), I=5, Nqs)
|
||
end if
|
||
|
||
END SUBROUTINE
|
||
|
||
|
||
END MODULE vdW_DF
|