scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
quantum-espresso/Modules/xc_vdW_DF.f90

3232 lines
112 KiB
Fortran
Raw Permalink Blame History

! Copyright (C) 2001-2009 Quantum ESPRESSO group
! Copyright (C) 2019 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 articles 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, Phys. Rev. Lett. 103, 096102 (2009).
!
!! henceforth referred to as SOLER.
!
!
! xc_vdW_DF and xc_vdW_DF_spin are the driver routines for vdW-DF
! calculations and are called from Modules/funct.f90. The routines here
! set up the parallel run (if any) and carry out the calls necessary to
! calculate the non-local correlation contributions to the energy and
! potential.
USE kinds, ONLY : dp
USE constants, ONLY : pi, fpi, e2
USE mp, ONLY : mp_sum, mp_barrier, mp_get, mp_size, mp_rank, mp_bcast
USE mp_images, ONLY : intra_image_comm
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
! ----------------------------------------------------------------------
! No implicit variables
IMPLICIT NONE
! ----------------------------------------------------------------------
! By default everything is private
PRIVATE
! ----------------------------------------------------------------------
! Save all objects in this module
SAVE
! ----------------------------------------------------------------------
! Public functions
PUBLIC :: xc_vdW_DF, xc_vdW_DF_spin, vdW_DF_stress, &
vdW_DF_energy, vdW_DF_potential, &
generate_kernel, interpolate_kernel, &
initialize_spline_interpolation, spline_interpolation, pw,pw_spin
! ----------------------------------------------------------------------
! Public variables
PUBLIC :: inlc, vdW_DF_analysis, Nr_points, r_max, q_min, q_cut, Nqs, q_mesh
! ----------------------------------------------------------------------
! General variables
INTEGER :: inlc = 1
! The non-local correlation
INTEGER :: vdW_DF_analysis = 0
! vdW-DF analysis tool as described in PRB 97, 085115 (2018)
REAL(DP), PARAMETER :: epsr = 1.0D-12
! A small number to cut off densities
INTEGER :: idx
! Indexing variable
! ----------------------------------------------------------------------
! Kernel specific parameters and variables
INTEGER, PARAMETER :: Nr_points = 1024
! The number of radial points (also the number of k points) used in the
! formation of the kernel functions for each pair of q values.
! Increasing this value will help in case you get a run-time error
! saying that you are trying to use a k value that is larger than the
! largest tabulated k point since the largest k point will be 2*pi/r_max
! * Nr_points. Memory usage of the vdW_DF piece of PWSCF will increase
! roughly linearly with this variable.
REAL(DP), PARAMETER :: r_max = 100.0D0
! The value of the maximum radius to use for the real-space kernel
! functions for each pair of q values. The larger this value is the
! smaller the smallest k value will be since the smallest k point value
! is 2*pi/r_max. Be careful though, since this will also decrease the
! maximum k point value and the vdW_DF code will crash if it encounters
! a g-vector with a magnitude greater than 2*pi/r_max *Nr_points.
REAL(DP), PARAMETER :: dr = r_max/Nr_points, dk = 2.0D0*pi/r_max
! Real space and k-space spacing of grid points.
REAL(DP), PARAMETER :: q_min = 1.0D-5, q_cut = 5.0D0
! The maximum and minimum values of q. During a vdW run, values of q0
! found larger than q_cut will be saturated (SOLER equation 5) to q_cut.
INTEGER, PARAMETER :: Nqs = 20
REAL(DP), PARAMETER, DIMENSION(Nqs) :: q_mesh = (/ &
q_min , 0.0449420825586261D0, 0.0975593700991365D0, 0.159162633466142D0, &
0.231286496836006D0, 0.315727667369529D0 , 0.414589693721418D0 , 0.530335368404141D0, &
0.665848079422965D0, 0.824503639537924D0 , 1.010254382520950D0 , 1.227727621364570D0, &
1.482340921174910D0, 1.780437058359530D0 , 2.129442028133640D0 , 2.538050036534580D0, &
3.016440085356680D0, 3.576529545442460D0 , 4.232271035198720D0 , q_cut /)
! The above two parameters define the q mesh to be used in the vdW_DF
! code. These are perhaps the most important to have set correctly.
! Increasing the number of q points will DRAMATICALLY increase the
! memory usage of the vdW_DF code because the memory consumption depends
! quadratically on the number of q points in the mesh. Increasing the
! number of q points may increase accuracy of the vdW_DF code, although,
! in testing it was found to have little effect. The largest value of
! the q mesh is q_cut. All values of q0 (DION equation 11) larger than
! this value during a run will be saturated to this value using equation
! 5 of SOLER. In testing, increasing the value of q_cut was found to
! have little impact on the results, although it is possible that in
! some systems it may be more important. Always make sure that the
! variable Nqs is consistent with the number of q points that are
! actually in the variable q_mesh. Also, do not set any q value to 0.
! This will cause an infinity in the Fourier transform.
INTEGER, PARAMETER :: Nintegration_points = 256
! Number of integration points for real-space kernel generation (see
! DION equation 14). This is how many a's and b's there will be.
REAL(DP), PARAMETER :: a_min = 0.0D0, a_max = 64.0D0
! Min/max values for the a and b integration in DION equation 14.
REAL(DP) :: kernel( 0:Nr_points, Nqs, Nqs ), d2phi_dk2( 0:Nr_points, Nqs, Nqs )
! Matrices holding the Fourier transformed kernel function and its
! second derivative for each pair of q values. The ordering is
! kernel(k_point, q1_value, q2_value).
REAL(DP) :: W_ab( Nintegration_points, Nintegration_points )
! Defined in DION equation 16.
REAL(DP) :: a_points( Nintegration_points ), a_points2( Nintegration_points )
! The values of the "a" points (DION equation 14) and their squares.
CONTAINS
! ####################################################################
! | |
! | functions |
! |_____________|
!
! Functions to be used in get_q0_on_grid, get_q0_on_grid_spin, and
! phi_value.
FUNCTION Fs(s)
IMPLICIT NONE
REAL(DP) :: s, Fs, Z_ab = 0.0D0
IF ( inlc == 1 .OR. inlc == 3 ) THEN
Z_ab = -0.8491D0
ELSE IF ( inlc == 2 .OR. inlc == 4 .OR. inlc == 5 ) THEN
Z_ab = -1.887D0
END IF
Fs = 1.0D0 - Z_ab * s * s / 9.0D0
END FUNCTION Fs
FUNCTION dFs_ds(s)
IMPLICIT NONE
REAL(DP) :: s, dFs_ds, Z_ab = 0.0D0
REAL(DP), PARAMETER :: prefac = -2.0D0/9.0D0
IF ( inlc == 1 .OR. inlc == 3 ) THEN
Z_ab = -0.8491D0
ELSE IF ( inlc == 2 .OR. inlc == 4 .OR. inlc == 5 ) THEN
Z_ab = -1.887D0
END IF
dFs_ds = prefac * s * Z_ab
END FUNCTION dFs_ds
FUNCTION kF(rho)
IMPLICIT NONE
REAL(DP) :: rho, kF
REAL(DP), PARAMETER :: ex = 1.0D0/3.0D0
kF = ( 3.0D0 * pi * pi * rho )**ex
END FUNCTION kF
FUNCTION dkF_drho(rho)
IMPLICIT NONE
REAL(DP) :: rho, dkF_drho
REAL(DP), PARAMETER :: prefac = 1.0D0/3.0D0
dkF_drho = prefac * 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 = 0.5D0 / (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
FUNCTION h_function(y)
IMPLICIT NONE
REAL(DP) :: y, y2, y4, h_function
REAL(DP), PARAMETER :: g1 = fpi/9.0D0 ! vdW-DF1/2
REAL(DP), PARAMETER :: a3 = 0.94950D0, g3 = 1.12D0, g32 = g3*g3 ! vdW-DF3-opt1
REAL(DP), PARAMETER :: a4 = 0.28248D0, g4 = 1.29D0, g42 = g4*g4 ! vdW-DF3-opt2
REAL(DP), PARAMETER :: a5 = 2.01059D0, b5 = 8.17471D0, g5 = 1.84981D0, & ! vdW-DF-C6
AA = ( b5 + a5*(a5/2.0D0-g5) ) / ( 1.0D0+g5-a5 ) !
y2 = y*y
IF ( inlc == 1 .OR. inlc == 2 ) THEN
h_function = 1.0D0 - EXP( -g1*y2 )
ELSE IF ( inlc == 3 ) THEN
y4 = y2*y2
h_function = 1.0D0 - 1.0D0 / ( 1.0D0 + g3*y2 + g32*y4 + a3*y4*y4 )
ELSE IF ( inlc == 4 ) THEN
y4 = y2*y2
h_function = 1.0D0 - 1.0D0 / ( 1.0D0 + g4*y2 + g42*y4 + a4*y4*y4 )
ELSE IF ( inlc == 5 ) THEN
y4 = y2*y2
h_function = 1.0D0 - ( 1.0D0 + ( (a5-g5)*y2 + AA*y4 ) / ( 1.0D0+AA*y2 ) ) * EXP( -a5*y2 )
END IF
END FUNCTION
! ####################################################################
! | |
! | XC_VDW_DF |
! |_____________|
SUBROUTINE xc_vdW_DF (rho_valence, rho_core, etxc, vtxc, v)
!
!! Driver routine for vdW-DF calculations, called from Modules/funct.f90.
!! The routine here sets up the parallel run (if any) and carry out the
!! calls necessary to calculate the non-local correlation contributions
!! to the energy and potential.
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.
! --------------------------------------------------------------------
! Write out the vdW-DF information and initialize the calculation.
IF ( first_iteration ) THEN
IF ( inlc > 5 ) CALL errore( 'xc_vdW_DF', 'inlc not implemented', 1 )
CALL generate_kernel
IF ( ionode ) CALL vdW_info(1)
first_iteration = .FALSE.
END IF
! --------------------------------------------------------------------
! Allocate arrays. nnr is a PWSCF variable that holds the number of
! points assigned to a given processor.
allocate( total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr), &
potential(dfftp%nnr), thetas(dfftp%nnr, Nqs), &
q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
! --------------------------------------------------------------------
! 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 theta arrays 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_DF_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 vdW_DF_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 ( total_rho, grad_rho, potential, thetas, q0, dq0_drho, dq0_dgradrho )
END SUBROUTINE xc_vdW_DF
! ####################################################################
! | |
! | XC_VDW_DF_spin |
! |__________________|
!
!
SUBROUTINE xc_vdW_DF_spin (rho_valence, rho_core, etxc, vtxc, v)
!! This subroutine is as similar to \(\texttt{xc_vdW_DF}\) as possible,
!! but handles the collinear \(\text{nspin}=2\) case.
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.
! --------------------------------------------------------------------
! Write out the vdW-DF information and initialize the calculation.
IF ( first_iteration ) THEN
IF ( inlc > 5 ) CALL errore( 'xc_vdW_DF_spin', 'inlc not implemented', 1 )
CALL generate_kernel
IF ( ionode ) CALL vdW_info(2)
first_iteration = .FALSE.
END IF
! --------------------------------------------------------------------
! Allocate arrays. nnr is a PWSCF variable that holds the number of
! points assigned to a given processor.
ALLOCATE( total_rho(dfftp%nnr), rho_up(dfftp%nnr), rho_down(dfftp%nnr), &
grad_rho(3,dfftp%nnr), grad_rho_up(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr), &
potential_up(dfftp%nnr), potential_down(dfftp%nnr), thetas(dfftp%nnr, Nqs), &
q0(dfftp%nnr), dq0_drho_up(dfftp%nnr), dq0_dgradrho_up(dfftp%nnr), &
dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )
! --------------------------------------------------------------------
! 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) + rho_valence(:,2) + rho_core(:) )*0.5D0
rho_down = ( rho_valence(:,1) - rho_valence(:,2) + rho_core(:) )*0.5D0
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_DF_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 vdW_DF_potential (q0, dq0_drho_up , dq0_dgradrho_up , grad_rho_up , thetas, potential_up )
CALL vdW_DF_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) + &
rho_valence(i_grid,2)) * 0.5_dp * potential_up (i_grid) &
+ e2 * grid_cell_volume * (rho_valence(i_grid,1) - &
rho_valence(i_grid,2)) * 0.5_dp * potential_down(i_grid)
END DO
DEALLOCATE( total_rho, rho_up, rho_down, grad_rho, grad_rho_up, grad_rho_down, &
potential_up, potential_down, thetas, &
q0, dq0_drho_up, dq0_dgradrho_up, dq0_drho_down, dq0_dgradrho_down )
END SUBROUTINE xc_vdW_DF_spin
! ####################################################################
! | |
! | GET_Q0_ON_GRID |
! |__________________|
!
!
SUBROUTINE get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
!! 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:
!
!! * \(\text{q0(ir)}\) = saturated value of q;
!! * \(\text{dq0_drho(ir)} = \text{total_rho} * d\text{q0} /d\text{rho}\)
!! * \(\text{dq0_dgradrho} = \text{total_rho} / |\text{grad_rho}| * d\text{q0}
!! / d |\text{grad_rho}|\)
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
REAL(DP) :: ec
REAL(DP) :: dq0_dq ! The derivative of the saturated
! q0 with respect to q.
INTEGER :: i_grid ! Indexing variable.
! --------------------------------------------------------------------
! 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
!-----------------------------------------------------------------------
SUBROUTINE pw( rs, iflag, ec, vc )
!-----------------------------------------------------------------------
! --A provisional copy of the pw routine in XC lib to avoid external calls--
!! * iflag=1: J.P. Perdew and Y. Wang, PRB 45, 13244 (1992)
!! * iflag=2: G. Ortiz and P. Ballone, PRB 50, 1391 (1994)
!
IMPLICIT NONE
REAL(DP), INTENT(IN) :: rs
INTEGER, INTENT(IN) :: iflag
REAL(DP), INTENT(OUT) :: ec, vc
!
REAL(DP), PARAMETER :: a=0.031091d0, b1=7.5957d0, b2=3.5876d0, c0=a, &
c1=0.046644d0, c2=0.00664d0, c3=0.01043d0, d0=0.4335d0, &
d1=1.4408d0
REAL(DP) :: lnrs, rs12, rs32, rs2, om, dom, olog
REAL(DP) :: a1(2), b3(2), b4(2)
DATA a1 / 0.21370d0, 0.026481d0 /, b3 / 1.6382d0, -0.46647d0 /, &
b4 / 0.49294d0, 0.13354d0 /
! high- and low-density formulae implemented but not used in PW case
! (reason: inconsistencies in PBE/PW91 functionals).
IF ( rs < 1d0 .AND. iflag == 2 ) THEN
! high density formula
lnrs = LOG(rs)
ec = c0 * lnrs - c1 + c2 * rs * lnrs - c3 * rs
vc = c0 * lnrs - (c1 + c0 / 3.d0) + 2.d0 / 3.d0 * c2 * rs * &
lnrs - (2.d0 * c3 + c2) / 3.d0 * rs
ELSEIF ( rs > 100.d0 .AND. iflag == 2 ) THEN
! low density formula
ec = - d0 / rs + d1 / rs**1.5d0
vc = - 4.d0 / 3.d0 * d0 / rs + 1.5d0 * d1 / rs**1.5d0
ELSE
! interpolation formula
rs12 = SQRT(rs)
rs32 = rs * rs12
rs2 = rs**2
om = 2.d0*a*( b1*rs12 + b2*rs + b3(iflag) * rs32 + b4(iflag)*rs2 )
dom = 2.d0*a*( 0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3(iflag) * &
rs32 + 2.d0 * b4(iflag) * rs2 )
olog = LOG( 1.d0 + 1.0d0 / om )
!
ec = - 2.d0 * a * (1.d0 + a1(iflag) * rs) * olog
vc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1(iflag) * rs) &
* olog - 2.d0 / 3.d0 * a * (1.d0 + a1(iflag) * rs) * dom / &
(om * (om + 1.d0) )
ENDIF
!
RETURN
!
END SUBROUTINE pw
!---------------------
! ####################################################################
! | |
! | 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)
!! Find the value of \(\text{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 \(\text{q0}\)s with respect to the
!! charge-density and the gradient of the charge-density. These are
!! needed for the potential.
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) :: fac
REAL(DP) :: ec, vc(2)
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 ! Indexing variable.
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, vc(1), vc(2) )
dqc_drho_up = vc(1)
dqc_drho_down = vc(2)
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
!-----------------------------------------------------------------------
SUBROUTINE pw_spin( rs, zeta, ec, vc_up, vc_dw )
!-----------------------------------------------------------------------
!--A provisional copy of the pw routine in XC lib to avoid external calls--
!! J.P. Perdew and Y. Wang, PRB 45, 13244 (1992).
IMPLICIT NONE
REAL(DP), INTENT(IN) :: rs
!! Wigner-Seitz radius
REAL(DP), INTENT(IN) :: zeta
!! zeta = (rho_up - rho_dw)/rho_tot
REAL(DP), INTENT(OUT) :: ec, vc_up, vc_dw
!
REAL(DP) :: rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz
REAL(DP) :: om, dom, olog, epwc, vpwc
REAL(DP) :: omp, domp, ologp, epwcp, vpwcp
REAL(DP) :: oma, doma, ologa, alpha, vpwca
!
! xc parameters, unpolarised
REAL(DP), PARAMETER :: a = 0.031091d0, a1 = 0.21370d0, b1 = 7.5957d0, b2 = &
3.5876d0, b3 = 1.6382d0, b4 = 0.49294d0, c0 = a, c1 = 0.046644d0, &
c2 = 0.00664d0, c3 = 0.01043d0, d0 = 0.4335d0, d1 = 1.4408d0
! xc parameters, polarised
REAL(DP), PARAMETER :: ap = 0.015545d0, a1p = 0.20548d0, b1p = 14.1189d0, b2p &
= 6.1977d0, b3p = 3.3662d0, b4p = 0.62517d0, c0p = ap, c1p = &
0.025599d0, c2p = 0.00319d0, c3p = 0.00384d0, d0p = 0.3287d0, d1p &
= 1.7697d0
! xc PARAMETERs, antiferro
REAL(DP), PARAMETER :: aa = 0.016887d0, a1a = 0.11125d0, b1a = 10.357d0, b2a = &
3.6231d0, b3a = 0.88026d0, b4a = 0.49671d0, c0a = aa, c1a = &
0.035475d0, c2a = 0.00188d0, c3a = 0.00521d0, d0a = 0.2240d0, d1a &
= 0.3969d0
REAL(DP), PARAMETER :: fz0 = 1.709921d0
!
! if (rs < 0.5d0) then
! high density formula (not implemented)
!
! elseif (rs > 100.d0) then
! low density formula (not implemented)
!
! else
! interpolation formula
!
zeta2 = zeta * zeta
zeta3 = zeta2 * zeta
zeta4 = zeta3 * zeta
rs12 = SQRT(rs)
rs32 = rs * rs12
rs2 = rs**2
! unpolarised
om = 2.d0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2)
dom = 2.d0 * a * (0.5d0 * b1 * rs12 + b2 * rs + 1.5d0 * b3 * rs32 &
+ 2.d0 * b4 * rs2)
olog = LOG(1.d0 + 1.0d0 / om)
epwc = - 2.d0 * a * (1.d0 + a1 * rs) * olog
vpwc = - 2.d0 * a * (1.d0 + 2.d0 / 3.d0 * a1 * rs) * olog - 2.d0 / &
3.d0 * a * (1.d0 + a1 * rs) * dom / (om * (om + 1.d0) )
! polarized
omp = 2.d0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2)
domp = 2.d0 * ap * (0.5d0 * b1p * rs12 + b2p * rs + 1.5d0 * b3p * &
rs32 + 2.d0 * b4p * rs2)
ologp = LOG(1.d0 + 1.0d0 / omp)
epwcp = - 2.d0 * ap * (1.d0 + a1p * rs) * ologp
vpwcp = - 2.d0 * ap * (1.d0 + 2.d0 / 3.d0 * a1p * rs) * ologp - &
2.d0 / 3.d0 * ap * (1.d0 + a1p * rs) * domp / (omp * (omp + 1.d0))
! antiferro
oma = 2.d0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2)
doma = 2.d0 * aa * ( 0.5d0 * b1a * rs12 + b2a * rs + 1.5d0 * b3a * &
rs32 + 2.d0 * b4a * rs2 )
ologa = LOG( 1.d0 + 1.0d0/oma )
alpha = 2.d0 * aa * (1.d0 + a1a*rs) * ologa
vpwca = + 2.d0 * aa * (1.d0 + 2.d0/3.d0 * a1a * rs) * ologa + &
2.d0 / 3.d0 * aa * (1.d0 + a1a*rs) * doma / (oma * (oma + 1.d0))
!
fz = ( (1.d0 + zeta)**(4.d0 / 3.d0) + (1.d0 - zeta)**(4.d0 / &
3.d0) - 2.d0) / (2.d0** (4.d0 / 3.d0) - 2.d0)
dfz = ( (1.d0 + zeta)**(1.d0 / 3.d0) - (1.d0 - zeta)**(1.d0 / &
3.d0) ) * 4.d0 / (3.d0 * (2.d0** (4.d0 / 3.d0) - 2.d0) )
!
ec = epwc + alpha * fz * (1.d0 - zeta4) / fz0 + (epwcp - epwc) &
* fz * zeta4
vc_up = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
* fz * zeta4 + (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
* (1.d0 - zeta)
vc_dw = vpwc + vpwca * fz * (1.d0 - zeta4) / fz0 + (vpwcp - vpwc) &
* fz * zeta4 - (alpha / fz0 * (dfz * (1.d0 - zeta4) - 4.d0 * fz * &
zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.d0 * fz * zeta3) ) &
* (1.d0 + zeta)
RETURN
!
END SUBROUTINE pw_spin
! ####################################################################
! | |
! | saturate_q |
! |______________|
SUBROUTINE saturate_q (q, q_cutoff, q0, dq0_dq)
!! Here, we calculate \(\text{q0}\) by saturating \(\text{q}\) according to
!! equation 5 of SOLER. Also, we find the derivative \(\text{dq0/dq}\) needed
!! for the derivatives \(\text{dq0/drho}\) and \(\text{dq0/dgradrh0}\).
IMPLICIT NONE
REAL(DP), INTENT(IN) :: q ! Input q.
REAL(DP), INTENT(IN) :: q_cutoff ! Cutoff q.
REAL(DP), INTENT(OUT) :: q0 ! Output saturated q.
REAL(DP), INTENT(OUT) :: dq0_dq ! Derivative of dq0/dq.
REAL(DP) :: e_exp ! Exponent.
INTEGER, PARAMETER :: m_cut = 12 ! How many terms to include in
! the sum of SOLER equation 5.
! --------------------------------------------------------------------
e_exp = 0.0D0
dq0_dq = 0.0D0
DO idx = 1, m_cut
e_exp = e_exp + (q/q_cutoff)**idx/idx
dq0_dq = dq0_dq + (q/q_cutoff)**(idx-1)
END Do
q0 = q_cutoff*(1.0D0 - EXP(-e_exp))
dq0_dq = dq0_dq * EXP(-e_exp)
END SUBROUTINE saturate_q
! ####################################################################
! | |
! | vdW_DF_energy |
! |_______________|
!
!
SUBROUTINE vdW_DF_energy (thetas, vdW_xc_energy)
!! This routine carries out the integration of equation 8 of SOLER. It
!! returns the non-local exchange-correlation energy and the \(\text{u_alpha(k)}\)
!! arrays used to find the \(\text{u_alpha(r)}\) arrays via equations 11 and 12 in
!! SOLER.
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).
ALLOCATE ( u_vdW(dfftp%nnr,Nqs), kernel_of_k(Nqs, Nqs) )
vdW_xc_energy = 0.0D0
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_DF_energy
! ####################################################################
! | |
! | vdW_DF_potential |
! |___________________|
!
!
SUBROUTINE vdW_DF_potential (q0, dq0_drho, dq0_dgradrho, grad_rho, u_vdW, 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 \(\text{u_alpha(k)}\) functions were found
!! while calculating the energy. They are passed in as the matrix \(\text{u_vdW}\).
!! Most of the required derivatives were calculated in the \(\texttt{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.
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('vdW_DF_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 vdW_DF_potential
! ####################################################################
! | |
! | SPLINE_INTERPOLATION |
! |________________________|
!
!
SUBROUTINE spline_interpolation (x, evaluation_points, values)
!! This routine is modeled after an algorithm from NUMERICAL_RECIPES 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.
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, 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 |
! |___________________________________|
!
SUBROUTINE initialize_spline_interpolation (x, d2y_dx2)
!! This routine is modeled after an algorithm from NUMERICAL RECIPES It
!! was adapted for Fortran and for the problem at hand.
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 ! 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 |
! |____________________|
!
SUBROUTINE interpolate_kernel (k, kernel_of_k)
!! This routine is modeled after an algorithm from NUMERICAL RECIPES
!! 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.
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
! ####################################################################
! | |
! | VDW_DF_STRESS |
! |_________________|
SUBROUTINE vdW_DF_stress (rho_valence, rho_core, nspin, sigma) ! PH adjusted for wrapper spin/nospin
!! vdW-DF stress calculation.
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 :: rho_up(:) !
real(dp), allocatable :: rho_down(:) !
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 :: q0(:) !
REAL(DP), ALLOCATABLE :: dq0_drho(:) ! q-values.
REAL(DP), ALLOCATABLE :: dq0_dgradrho(:) !
real(dp), allocatable :: dq0_drho_up(:) ! The derivative of the saturated q0
real(dp), allocatable :: dq0_drho_down(:) ! with respect to the spin charge density
real(dp), allocatable :: dq0_dgradrho_up(:) ! The derivative of the saturated q0
real(dp), allocatable :: dq0_dgradrho_down(:) ! with respect to the gradient of the spin charge sensity
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 ( inlc > 5 ) CALL errore( 'xc_vdW_DF', 'inlc not implemented', 1 )
#if defined (__SPIN_BALANCED)
IF ( nspin==2 ) THEN
WRITE(stdout,'(/,/ " Performing spin-balanced Ecnl stress calculation!")')
ELSE IF ( nspin > 2 ) THEN
CALL errore ('vdW_DF_stress', 'noncollinear vdW stress not implemented', 1)
END IF
#else
IF ( nspin>2 ) THEN
CALL errore ('vdW_DF_stress', 'vdW stress not implemented for nspin > 2', 1)
END IF
#endif
sigma(:,:) = 0.0_DP
sigma_grad(:,:) = 0.0_DP
sigma_ker(:,:) = 0.0_DP
! --------------------------------------------------------------------
! Allocations
ALLOCATE( total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr), thetas(dfftp%nnr, Nqs), q0(dfftp%nnr) )
ALLOCATE ( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
#if defined (__SPIN_BALANCED)
#else
IF (nspin==2) THEN
allocate( rho_up(dfftp%nnr), rho_down(dfftp%nnr) )
allocate( grad_rho_up(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr) )
allocate( dq0_drho_up (dfftp%nnr), dq0_dgradrho_up (dfftp%nnr) )
allocate( dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )
ENDIF
#endif
! --------------------------------------------------------------------
! Charge
total_rho = rho_valence(:,1) + rho_core(:)
#if defined (__SPIN_BALANCED)
#else
IF (nspin==2) THEN
rho_up = ( rho_valence(:,1) + rho_valence(:,2) + rho_core(:) )*0.5D0
rho_down = ( rho_valence(:,1) - rho_valence(:,2) + rho_core(:) )*0.5D0
ENDIF
#endif
! --------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT.
CALL fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
#if defined (__SPIN_BALANCED)
#else
IF (nspin==2) THEN
call fft_gradient_r2r (dfftp, rho_up, g, grad_rho_up)
call fft_gradient_r2r (dfftp, rho_down, g, grad_rho_down)
ENDIF
#endif
! --------------------------------------------------------------------
! Get q0.
#if defined (__SPIN_BALANCED)
CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
#else
IF (nspin == 1) THEN
CALL get_q0_on_grid (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, thetas)
ELSEIF (nspin==2) THEN
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)
ENDIF
#endif
! --------------------------------------------------------------------
! Stress
#if defined (__SPIN_BALANCED)
CALL vdW_DF_stress_gradient (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, &
thetas, sigma_grad)
#else
IF (nspin == 1) THEN
CALL vdW_DF_stress_gradient (total_rho, grad_rho, q0, dq0_drho, dq0_dgradrho, &
thetas, sigma_grad)
ELSEIF (nspin == 2) THEN
CALL vdW_DF_stress_gradient_spin (total_rho, grad_rho_up, grad_rho_down, q0, &
dq0_dgradrho_up, dq0_dgradrho_down, &
thetas, sigma_grad)
ENDIF
#endif
CALL vdW_DF_stress_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)
END DO
END DO
DEALLOCATE( total_rho, grad_rho, thetas, q0 )
DEALLOCATE( dq0_drho, dq0_dgradrho )
#if defined (__SPIN_BALANCED)
#else
IF (nspin == 2) THEN
deallocate( rho_up, rho_down )
deallocate( grad_rho_up, grad_rho_down )
deallocate( dq0_drho_up, dq0_drho_down, dq0_dgradrho_up, dq0_dgradrho_down )
ENDIF
#endif
END SUBROUTINE vdW_DF_stress ! PH adjusted for wrapper spin/nospin
! -------------------------------------------------------------------------
! Begin Spin vdW-DF_strees_gradient implemented Per Hyldgaard 2019, GPL. No Waranties
! Adapted from the original nspin = 1 code (subroutine below) by Thonhauser and coauthors
! -------------------------------------------------------------------------
! ####################################################################
! | |
! | VDW_DF_STRESS_GRADIENT_SPIN |
! |_______________________________|
SUBROUTINE vdW_DF_stress_gradient_spin (total_rho, grad_rho_up, grad_rho_down, q0, &
dq0_dgradrho_up, dq0_dgradrho_down, &
thetas, sigma)
USE gvect, ONLY : ngm, g, gg, gstart
USE cell_base, ONLY : omega, tpiba, alat, at, tpiba2
implicit none
real(dp), intent(IN) :: total_rho(:) !
real(dp), intent(IN) :: grad_rho_up (:, :) ! Input variables.
real(dp), intent(IN) :: grad_rho_down(:, :) !
real(dp), intent(inout) :: sigma(:,:) !
real(dp), intent(IN) :: q0(:) !
real(dp), intent(IN) :: dq0_dgradrho_up(:) !
real(dp), intent(IN) :: dq0_dgradrho_down(:) !
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_up, prefactor_down ! Final summation of sigma.
real(dp) :: grad2_up, grad2_down ! 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_up = 0.0_DP
prefactor_down = 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_up = sqrt( grad_rho_up(1,i_grid)**2 &
+ grad_rho_up(2,i_grid)**2 + grad_rho_up(3,i_grid)**2 )
grad2_down = sqrt( grad_rho_down(1,i_grid)**2 &
+ grad_rho_down(2,i_grid)**2 + grad_rho_down(3,i_grid)**2 )
if ( grad2_up == 0.0_dp .or. grad2_down == 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('vdW_DF_stress_gradient_spin','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_up = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho_up(i_grid) / grad2_up
prefactor_down = u_vdW(i_grid,q_i) * dP_dq0 * dq0_dgradrho_down(i_grid) / grad2_down
do l = 1, 3
do m = 1, l
sigma (l, m) = sigma (l, m) - e2 * prefactor_up * &
(grad_rho_up(l,i_grid) * grad_rho_up(m,i_grid))
sigma (l, m) = sigma (l, m) - e2 * prefactor_down * &
(grad_rho_down(l,i_grid) * grad_rho_down(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 vdW_DF_stress_gradient_spin
! ####################################################################
! | |
! | VDW_DF_STRESS_GRADIENT |
! |__________________________|
SUBROUTINE vdW_DF_stress_gradient (total_rho, grad_rho, q0, &
dq0_drho, dq0_dgradrho, thetas, sigma)
USE gvect, ONLY : ngm, g, gstart
USE cell_base, ONLY : omega, tpiba, alat, at, tpiba2
IMPLICIT NONE
REAL(DP), INTENT(IN) :: total_rho(:) !
REAL(DP), INTENT(IN) :: grad_rho(:, :) !
REAL(DP), INTENT(INOUT) :: sigma(:,:) !
REAL(DP), INTENT(IN) :: q0(:) ! Input variables.
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, & ! Iterators.
ix, iy, iz !
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.0D0 / (dfftp%nr1 * dfftp%nr2 * dfftp%nr3), sigma, 1)
DEALLOCATE( d2y_dx2, u_vdW )
END SUBROUTINE vdW_DF_stress_gradient
! ####################################################################
! | |
! | VDW_DF_STRESS_KERNEL |
! |________________________|
SUBROUTINE vdW_DF_stress_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
INTEGER :: last_g, theta_i
REAL(DP) :: g2, ngmod2, g_kernel, G_multiplier
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
END DO
END DO
IF ( g_i < gstart ) sigma(:,:) = sigma(:,:) / G_multiplier
END DO
CALL mp_sum( sigma, intra_bgrp_comm )
DEALLOCATE( dkernel_of_dk )
END SUBROUTINE vdW_DF_stress_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
! ####################################################################
! | |
! | GENERATE_KERNEL |
! |_________________|
!
! The original definition of the kernel function is given in DION
! equations 14-16. The Soler method makes the kernel function a
! function of only 1 variable (r) by first putting it in the form
! phi(q1*r, q2*r). Then, the q-dependence is removed by expanding the
! function in a special way (see SOLER equation 3). This yields a
! separate function for each pair of q points that is a function of r
! alone. There are (Nqs^2+Nqs)/2 unique functions, where Nqs is the
! number of q points used. In the Soler method, the kernel is first
! made in the form phi(d1, d2) but this is not done here. It was found
! that, with q's chosen judiciously ahead of time, the kernel and the
! second derivatives required for interpolation could be tabulated
! ahead of time for faster use of the vdW-DF functional. Through
! testing we found no need to soften the kernel and correct for this
! later (see SOLER eqations 6-7).
!
! The algorithm employed here is "embarrassingly parallel," meaning
! that it parallelizes very well up to (Nqs^2+Nqs)/2 processors,
! where, again, Nqs is the number of q points chosen. However,
! parallelization on this scale is unnecessary. In testing the code
! runs in under a minute on 16 Intel Xeon processors.
!
! IMPORTANT NOTICE: Results are very sensitive to compilation details.
! In particular, the usage of FMA (Fused Multiply-and-Add)
! instructions used by modern CPUs such as AMD Interlagos (Bulldozer)
! and Intel Ivy Bridge may affect quite heavily some components of the
! kernel (communication by Ake Sandberg, Umea University). In practice
! this should not be a problem, since most affected elements are the
! less relevant ones.
!
! Some of the algorithms here are somewhat modified versions of those
! found in:
!
! Numerical Recipes in C; William H. Press, Brian P. Flannery, Saul
! A. Teukolsky, and William T. Vetterling. Cambridge University
! Press (1988).
!
! hereafter referred to as NUMERICAL_RECIPES. The routines were
! translated to Fortran, of course and variable names are generally
! different.
!
! For the calculation of the kernel we have benefited from access to
! earlier vdW-DF implementation into PWscf and ABINIT, written by Timo
! Thonhauser, Valentino Cooper, and David Langreth. These codes, in
! turn, benefited from earlier codes written by Maxime Dion and Henrik
! Rydberg.
SUBROUTINE generate_kernel
IMPLICIT NONE
INTEGER :: a_i, b_i, q1_i, q2_i, r_i
! Indexing variables.
REAL(DP) :: weights( Nintegration_points )
! Array to hold dx values for the Gaussian-Legendre integration of the kernel.
REAL(DP) :: sin_a( Nintegration_points ), cos_a( Nintegration_points )
! Sine and cosine values of the aforementioned points a.
REAL(DP) :: d1, d2, d, integral
! Intermediate values.
! --------------------------------------------------------------------
! The following variables control the parallel environment.
INTEGER :: my_start_q, my_end_q, Ntotal
! Starting and ending q value for each processor, also the total
! number of calculations to do, i.e. (Nqs^2 + Nqs)/2.
REAL(DP), ALLOCATABLE :: phi(:,:), phi_deriv(:,:)
! Arrays to store the kernel functions and their second derivatives.
! They are stored as phi(radial_point, idx).
INTEGER, ALLOCATABLE :: indices(:,:), proc_indices(:,:)
! Indices holds the values of q1 and q2 as partitioned out to the
! processors. It is an Ntotal x 2 array stored as indices(index of
! point number, q1:q2). Proc_indices holds the section of the indices
! array that is assigned to each processor. This is a Nproc x 2
! array, stored as proc_indices(processor_number,
! starting_index:ending_index)
INTEGER :: Nper, Nextra, start_q, end_q
! Baseline number of jobs per processor, number of processors that
! get an extra job in case the number of jobs doesn't split evenly
! over the number of processors, starting index into the indices
! array, ending index into the indices array.
INTEGER :: nproc, mpime
! Number or procs, rank of current processor.
INTEGER :: proc_i, my_Nqs
! --------------------------------------------------------------------
! Start the timer.
CALL start_clock ( 'vdW_kernel' )
! --------------------------------------------------------------------
! The total number of phi_alpha_beta functions that have to be
! calculated.
Ntotal = (Nqs**2 + Nqs)/2
ALLOCATE ( indices(Ntotal, 2) )
! --------------------------------------------------------------------
! This part fills in the indices array. It just loops through the q1
! and q2 values and stores them. Sections of this array will be
! assigned to each of the processors later.
idx = 1
DO q1_i = 1, Nqs
DO q2_i = 1, q1_i
indices(idx, 1) = q1_i
indices(idx, 2) = q2_i
idx = idx + 1
END DO
END DO
! --------------------------------------------------------------------
! Figure out the baseline number of functions to be calculated by each
! processor and how many processors get one extra job.
nproc = mp_size( intra_image_comm )
mpime = mp_rank( intra_image_comm )
Nper = Ntotal/nproc
Nextra = MOD(Ntotal, nproc)
ALLOCATE( proc_indices(nproc, 2) )
start_q = 0
end_q = 0
! --------------------------------------------------------------------
! Loop over all the processors and figure out which section of the
! indices array each processor should do. All processors figure this
! out for every processor so there is no need to communicate results.
DO proc_i = 1, nproc
start_q = end_q + 1
end_q = start_q + (Nper - 1)
IF (proc_i <= Nextra) end_q = end_q + 1
! This is to prevent trouble if number of processors exceeds Ntotal.
IF ( proc_i > Ntotal ) THEN
start_q = Ntotal
end_q = Ntotal
END IF
IF ( proc_i == (mpime+1) ) THEN
my_start_q = start_q
my_end_q = end_q
END IF
proc_indices(proc_i, 1) = start_q
proc_indices(proc_i, 2) = end_q
END DO
! --------------------------------------------------------------------
! Store how many jobs are assigned to me.
my_Nqs = my_end_q - my_start_q + 1
ALLOCATE( phi( 0:Nr_points, my_Nqs ), phi_deriv( 0:Nr_points, my_Nqs ) )
phi = 0.0D0
phi_deriv = 0.0D0
kernel = 0.0D0
d2phi_dk2 = 0.0D0
! --------------------------------------------------------------------
! Find the integration points we are going to use in the
! Gaussian-Legendre integration.
CALL prep_gaussian_quadrature( weights )
! --------------------------------------------------------------------
! Get a, a^2, sin(a), cos(a) and the weights for the Gaussian-Legendre
! integration.
DO a_i=1, Nintegration_points
a_points (a_i) = TAN( a_points(a_i) )
a_points2(a_i) = a_points(a_i)**2
weights(a_i) = weights(a_i)*( 1 + a_points2(a_i) )
cos_a(a_i) = COS( a_points(a_i) )
sin_a(a_i) = SIN( a_points(a_i) )
END DO
! --------------------------------------------------------------------
! Calculate the value of the W function defined in DION equation 16
! for each value of a and b.
DO a_i = 1, Nintegration_points
DO b_i = 1, Nintegration_points
W_ab(a_i, b_i) = 2.0D0 * weights(a_i)*weights(b_i) * ( &
(3.0D0-a_points2(a_i))*a_points(b_i) *sin_a(a_i)*cos_a(b_i) + &
(3.0D0-a_points2(b_i))*a_points(a_i) *cos_a(a_i)*sin_a(b_i) + &
(a_points2(a_i)+a_points2(b_i)-3.0D0)*sin_a(a_i)*sin_a(b_i) - &
3.0D0*a_points(a_i)*a_points(b_i)*cos_a(a_i)*cos_a(b_i) ) / &
(a_points(a_i)*a_points(b_i) )
END DO
END DO
! --------------------------------------------------------------------
! vdW-DF analysis tool as described in PRB 97, 085115 (2018).
IF ( vdW_DF_analysis == 1 ) THEN
DO a_i = 1, Nintegration_points
DO b_i = 1, Nintegration_points
W_ab(a_i, b_i) = weights(a_i)*weights(b_i) * &
a_points(a_i)*a_points(b_i)*sin_a(a_i)*sin_a(b_i)
END DO
END DO
ELSE IF ( vdW_DF_analysis == 2 ) THEN
DO a_i = 1, Nintegration_points
DO b_i = 1, Nintegration_points
W_ab(a_i, b_i) = W_ab(a_i, b_i) - weights(a_i)*weights(b_i) * &
a_points(a_i)*a_points(b_i)*sin_a(a_i)*sin_a(b_i)
END DO
END DO
END IF
! --------------------------------------------------------------------
! Now, we loop over all the pairs (q1,q2) that are assigned to us and
! perform our calculations.
DO idx = 1, my_Nqs
! -----------------------------------------------------------------
! First, get the value of phi(q1*r, q2*r) for each r and the
! particular values of q1 and q2 we are using.
DO r_i = 1, Nr_points
d1 = q_mesh( indices(idx+my_start_q-1, 1) ) * dr * r_i
d2 = q_mesh( indices(idx+my_start_q-1, 2) ) * dr * r_i
phi(r_i, idx) = phi_value(d1, d2)
END DO
! -----------------------------------------------------------------
! Now, perform a radial FFT to turn our phi_alpha_beta(r) into
! phi_alpha_beta(k) needed for SOLER equation 8.
CALL radial_fft( phi(:,idx) )
! -----------------------------------------------------------------
! Determine the spline interpolation coefficients for the Fourier
! transformed kernel function.
CALL set_up_splines( phi(:, idx), phi_deriv(:, idx) )
END DO
! --------------------------------------------------------------------
! Finally, we collect the results after letting everybody catch up.
CALL mp_barrier( intra_image_comm )
DO proc_i = 0, nproc-1
IF ( proc_i >= Ntotal ) EXIT
CALL mp_get ( phi , phi , mpime, 0, proc_i, 0, intra_image_comm )
CALL mp_get ( phi_deriv, phi_deriv, mpime, 0, proc_i, 0, intra_image_comm )
IF ( mpime == 0 ) THEN
DO idx = proc_indices(proc_i+1,1), proc_indices(proc_i+1,2)
q1_i = indices(idx, 1)
q2_i = indices(idx, 2)
kernel (:, q1_i, q2_i) = phi (:, idx - proc_indices(proc_i+1,1) + 1)
d2phi_dk2 (:, q1_i, q2_i) = phi_deriv (:, idx - proc_indices(proc_i+1,1) + 1)
kernel (:, q2_i, q1_i) = kernel (:, q1_i, q2_i)
d2phi_dk2 (:, q2_i, q1_i) = d2phi_dk2 (:, q1_i, q2_i)
END DO
END IF
END DO
CALL mp_bcast ( kernel , 0, intra_image_comm )
CALL mp_bcast ( d2phi_dk2, 0, intra_image_comm )
! --------------------------------------------------------------------
! Keep the lines below for testing and combatibility with the old
! kernel file reading/writing method.
!
! Writing the calculated kernel.
!
! IF ( ionode ) THEN
! WRITE(stdout,'(/ / A)') " vdW-DF kernel table calculated and written to file."
! OPEN(UNIT=21, FILE='kernel_table', STATUS='replace', FORM='formatted', ACTION='write')
! WRITE(21, '(2i5,f13.8)') Nqs, Nr_points
! WRITE(21, '(1p4e23.14)') r_max
! WRITE(21, '(1p4e23.14)') q_mesh
! DO q1_i = 1, Nqs
! DO q2_i = 1, q1_i
! WRITE(21, '(1p4e23.14)') kernel(:, q1_i, q2_i)
! END DO
! END DO
! DO q1_i = 1, Nqs
! DO q2_i = 1, q1_i
! WRITE(21, '(1p4e23.14)') d2phi_dk2(:, q1_i, q2_i)
! END DO
! END DO
! CLOSE (21)
! END IF
!
!
! Reading the kernel from an old kernel file.
!
! IF (ionode) WRITE(stdout,'(/ / A)') " vdW-DF kernel read from file."
! OPEN(UNIT=21, FILE='vdW_kernel_table', STATUS='old', FORM='formatted', ACTION='read')
! read(21, '(/ / / / / /)')
! DO q1_i = 1, Nqs
! DO q2_i = 1, q1_i
! READ(21, '(1p4e23.14)') kernel(:, q1_i, q2_i)
! kernel(:, q2_i, q1_i) = kernel(:, q1_i, q2_i)
! END DO
! END DO
! DO q1_i = 1, Nqs
! DO q2_i = 1, q1_i
! READ(21, '(1p4e23.14)') d2phi_dk2(:, q1_i, q2_i)
! d2phi_dk2(:, q2_i, q1_i) = d2phi_dk2(:, q1_i, q2_i)
! END DO
! END DO
! CLOSE (21)
DEALLOCATE( indices, proc_indices, phi, phi_deriv )
! --------------------------------------------------------------------
! Stop the timer.
CALL stop_clock ( 'vdW_kernel' )
END SUBROUTINE generate_kernel
! ####################################################################
! | |
! | PREP_GAUSSIAN_QUADRATURE |
! |____________________________|
!
SUBROUTINE prep_gaussian_quadrature( weights )
!! Routine to calculate the points and weights for the
!! Gaussian-Legendre integration. This routine is modeled after the
!! routine GAULEG from NUMERICAL RECIPES.
REAL(DP), INTENT(INOUT) :: weights(:)
! The points and weights for the Gaussian-Legendre integration.
INTEGER :: Npoints
! The number of points we actually have to calculate. The rest will
! be obtained from symmetry.
REAL(DP) :: poly_1, poly_2, poly_3
! Temporary storage for Legendre polynomials.
INTEGER :: i_point, i_poly
! Indexing variables.
REAL(DP) :: root, dp_dx, last_root
! The value of the root of a given Legendre polynomial, the derivative
! of the polynomial at that root and the value of the root in the last
! iteration (to check for convergence of Newton's method).
real(dp) :: midpoint, length
! The middle of the x-range and the length to that point.
Npoints = (Nintegration_points + 1)/2
midpoint = 0.5D0 * ( ATAN(a_min) + ATAN(a_max) )
length = 0.5D0 * ( ATAN(a_max) - ATAN(a_min) )
DO i_point = 1, Npoints
! -----------------------------------------------------------------
! Make an initial guess for the root.
root = COS(DBLE(pi*(i_point - 0.25D0)/(Nintegration_points + 0.5D0)))
DO
! --------------------------------------------------------------
! Use the recurrence relations to find the desired polynomial,
! evaluated at the approximate root. See NUMERICAL_RECIPES.
poly_1 = 1.0D0
poly_2 = 0.0D0
DO i_poly = 1, Nintegration_points
poly_3 = poly_2
poly_2 = poly_1
poly_1 = ((2.0D0 * i_poly - 1.0D0)*root*poly_2 - (i_poly-1.0D0)*poly_3)/i_poly
END DO
! --------------------------------------------------------------
! Use the recurrence relations to find the desired polynomial.
! Find the derivative of the polynomial and use it in Newton's
! method to refine our guess for the root.
dp_dx = Nintegration_points * (root*poly_1 - poly_2)/(root**2 - 1.0D0)
last_root = root
root = last_root - poly_1/dp_dx
! --------------------------------------------------------------
! Check for convergence.
IF (abs(root - last_root) <= 1.0D-14) EXIT
END DO
! -----------------------------------------------------------------
! Fill in the array of evaluation points.
a_points(i_point) = midpoint - length*root
a_points(Nintegration_points + 1 - i_point) = midpoint + length*root
! -----------------------------------------------------------------
! Fill in the array of weights.
weights(i_point) = 2.0D0 * length/((1.0D0 - root**2)*dp_dx**2)
weights(Nintegration_points + 1 - i_point) = weights(i_point)
END DO
END SUBROUTINE prep_gaussian_quadrature
! ####################################################################
! | |
! | PHI_VALUE |
! |_____________|
!
REAL(DP) FUNCTION phi_value(d1, d2)
!! This function returns the value of the kernel calculated via DION
!! equation 14.
REAL(DP), INTENT(IN) :: d1, d2
! The point at which to evaluate the kernel. d1 = q1*r and d2 = q2*r.
REAL(DP) :: w, x, y, z, T
! Intermediate values.
REAL(DP) :: nu(Nintegration_points), nu1(Nintegration_points)
! Defined in the discussio below equation 16 of DION.
INTEGER :: a_i, b_i
! Indexing variables.
! --------------------------------------------------------------------
! Loop over all integration points and calculate the value of the nu
! functions defined in the discussion below equation 16 in DION.
DO a_i = 1, Nintegration_points
nu(a_i) = a_points2(a_i)/( 2.0D0 * h_function( a_points(a_i)/d1 ))
nu1(a_i) = a_points2(a_i)/( 2.0D0 * h_function( a_points(a_i)/d2 ))
END DO
! --------------------------------------------------------------------
! Carry out the integration of DION equation 13.
phi_value = 0.0D0
DO a_i = 1, Nintegration_points
w = nu(a_i)
y = nu1(a_i)
DO b_i = 1, Nintegration_points
x = nu(b_i)
z = nu1(b_i)
T = (1.0D0/(w+x) + 1.0D0/(y+z))*(1.0D0/((w+y)*(x+z)) + 1.0D0/((w+z)*(y+x)))
phi_value = phi_value + T * W_ab(a_i, b_i)
END DO
END DO
phi_value = 1.0D0/pi**2*phi_value
END FUNCTION phi_value
! ####################################################################
! | |
! | RADIAL_FFT |
! |______________|
!
SUBROUTINE radial_fft(phi)
!! This subroutine performs a radial Fourier transform on the
!! real-space kernel functions.
! Basically, this is just:
! int(4*pi*r^2*phi*sin(k*r)/(k*r))dr
! integrated from 0 to r_max.
! That is, it is the kernel function phi integrated with the 0^th spherical
! Bessel function radially, with a 4*pi assumed from angular
! integration since we have spherical symmetry. The spherical symmetry
! comes in because the kernel function depends only on the magnitude
! of the vector between two points. The integration is done using the
! trapezoid rule.
REAL(DP), INTENT(INOUT) :: phi(0:Nr_points)
! On input holds the real-space function phi_q1_q2(r).
! On output hold the reciprocal-space function phi_q1_q2(k).
REAL(DP) :: phi_k(0:Nr_points)
! Temporary storage for phi_q1_q2(k).
INTEGER :: k_i, r_i
! Indexing variables.
REAL(DP) :: r, k
! The real and reciprocal space points.
phi_k = 0.0D0
! --------------------------------------------------------------------
! Handle the k=0 point separately.
DO r_i = 1, Nr_points
r = r_i * dr
phi_k(0) = phi_k(0) + phi(r_i)*r**2
END DO
! --------------------------------------------------------------------
! Subtract half of the last value off because of the trapezoid rule.
phi_k(0) = phi_k(0) - 0.5D0 * (Nr_points*dr)**2 * phi(Nr_points)
! --------------------------------------------------------------------
! Integration for the rest of the k-points.
DO k_i = 1, Nr_points
k = k_i * dk
DO r_i = 1, Nr_points
r = r_i * dr
phi_k(k_i) = phi_k(k_i) + phi(r_i) * r * SIN(k*r) / k
END DO
phi_k(k_i) = phi_k(k_i) - 0.5D0 * phi(Nr_points) * r * SIN(k*r) / k
END DO
! --------------------------------------------------------------------
! Add in the 4*pi and the dr factor for the integration.
phi = 4.0D0 * pi * phi_k * dr
END SUBROUTINE radial_fft
! ####################################################################
! | |
! | SET UP SPLINES |
! |__________________|
!
SUBROUTINE set_up_splines(phi, D2)
!! This subroutine accepts a function (phi) and finds at each point the
!! second derivative (D2) for use with spline interpolation. This
!! function assumes we are using the expansion described in SOLER
!! equation 3.
! That is, the derivatives are those needed to interpolate
! Kronecker delta functions at each of the q values. Other than some
! special modification to speed up the algorithm in our particular
! case, this algorithm is taken directly from NUMERICAL_RECIPES.
REAL(DP), INTENT(IN) :: phi(0:Nr_points)
! The k-space kernel function for a particular q1 and q2.
REAL(DP), INTENT(INOUT) :: D2(0:Nr_points)
! The second derivatives to be used in the interpolation expansion
! (SOLER equation 3).
REAL(DP), ALLOCATABLE :: temp_array(:) ! Temporary storage.
REAL(DP) :: temp_1, temp_2
INTEGER :: r_i
! Indexing variable.
ALLOCATE( temp_array(0:Nr_points) )
D2 = 0
temp_array = 0
DO r_i = 1, Nr_points - 1
temp_1 = DBLE(r_i - (r_i - 1))/DBLE( (r_i + 1) - (r_i - 1) )
temp_2 = temp_1 * D2(r_i-1) + 2.0D0
D2(r_i) = (temp_1 - 1.0D0)/temp_2
temp_array(r_i) = ( phi(r_i+1) - phi(r_i))/DBLE( dk*((r_i+1) - r_i) ) - &
( phi(r_i) - phi(r_i-1))/DBLE( dk*(r_i - (r_i-1)) )
temp_array(r_i) = (6.0D0*temp_array(r_i)/DBLE( dk*((r_i+1) - (r_i-1)) )-&
temp_1*temp_array(r_i-1))/temp_2
END DO
D2(Nr_points) = 0.0D0
DO r_i = Nr_points-1, 0, -1
D2(r_i) = D2(r_i)*D2(r_i+1) + temp_array(r_i)
END DO
DEALLOCATE( temp_array )
END SUBROUTINE set_up_splines
! ####################################################################
! | |
! | VDW_INFO |
! |____________|
SUBROUTINE vdW_info( nspin )
USE xc_lib, ONLY : xclib_dft_is
IMPLICIT NONE
INTEGER, INTENT (IN) :: nspin
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,"% If you are calculating stress with vdW-DF, please also cite: %")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"% R. Sabatini et al., J. Phys.: Condens. Matter 24, 424209 (2012). %")')
IF (nspin > 1) WRITE(stdout,'(5x,"% C.M. Frostenson et al., Electr. Struct. 4, 014001 (2022). %")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%",/)')
WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"% vdW-DF NEWS: %")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"% * vdW-DF3 is now available. DOI: 10.1021/acs.jctc.0c00471 %")')
WRITE(stdout,'(5x,"% use with input_dft = ''vdW-DF3-opt1'' or ''vdW-DF3-opt2'' %")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"% * Unscreened and range-separated hybrid vdW-DF-cx functionals %")')
WRITE(stdout,'(5x,"% DOI: 10.1063/1.4986522 and 10.1088/1361-648X/ac2ad2 %")')
WRITE(stdout,'(5x,"% use with input_dft = ''vdW-DF-cx0'' and ''vdW-DF-ahcx'' %")')
WRITE(stdout,'(5x,"% * Unscreened and range-separated hybrid vdW-DF2-b86r functionals %")')
WRITE(stdout,'(5x,"% DOI: 10.1063/1.4986522 and DOI: 10.1103/PhysRevX.12.041003 %")')
WRITE(stdout,'(5x,"% use with input_dft = ''vdW-DF2-br0'' and ''vdW-DF2-ahbr'' %")')
WRITE(stdout,'(5x,"% %")')
WRITE(stdout,'(5x,"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%",//)')
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(idx), idx=1, 4)
WRITE(stdout,'(13X,4F12.8)') (q_mesh(idx), idx=5, Nqs)
END IF
END SUBROUTINE
END MODULE vdW_DF
Reference in New Issue View Git Blame Copy Permalink