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_rVV10.f90

1405 lines
48 KiB
Fortran
Raw Permalink Blame History

!
! Copyright (C) 2001-2009 Quantum ESPRESSO group
! 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 rVV10
!--------------------------------------------------------------------------
!! This module is modeled after the vdW-DF implementation in
!! 'Modules/xc_vdW_DF.f90'. See that file for references, explanations, and
!! many useful comments.
USE kinds, ONLY : dp
USE constants, ONLY : pi
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
USE io_global, ONLY : ionode, stdout
USE fft_base, ONLY : dfftp
USE fft_interfaces, ONLY : fwfft, invfft
USE control_flags, ONLY : gamma_only, iverbosity
IMPLICIT NONE
SAVE
real(dp), parameter :: epsr = 1.d-12
real(dp), parameter :: epsg = 1.D-10
integer, parameter :: Nr_points = 1024
real(dp), parameter :: r_max = 100.0D0
real(dp), parameter :: dr = r_max/Nr_points
real(dp), parameter :: dk = 2.0D0*pi/r_max
real(dp), parameter :: q_min = 1.0D-4
real(dp), parameter :: q_cut = 0.5D0
integer, parameter :: Nqs = 20
real(dp), parameter, dimension(Nqs):: q_mesh= (/ q_min, 3.0D-4, 5.893850845618885D-4, 1.008103720396345D-3, &
1.613958359589310D-3, 2.490584839564653D-3, 3.758997979748929D-3, 5.594297198907115D-3, &
8.249838297569416D-3, 1.209220822453922D-2, 1.765183095571029D-2, 2.569619042667097D-2, &
3.733577865542191D-2, 5.417739477463518D-2, 7.854595729872216D-2, 0.113805449932145D0, &
0.164823306218807D0 , 0.238642339497217D0 , 0.345452975434964D0 , q_cut /)
real(dp) :: kernel( 0:Nr_points, Nqs, Nqs ), d2phi_dk2( 0:Nr_points, Nqs, Nqs )
real(dp) :: b_value = 6.3_DP
real(dp) :: C_value = 0.0093
private
public :: xc_rVV10, &
interpolate_kernel, &
initialize_spline_interpolation, &
rVV10_stress, b_value, &
q_mesh, Nr_points, r_max, q_min, q_cut, Nqs
CONTAINS
! #################################################################################################
! | |
! | xc_rVV10 |
! |_____________|
SUBROUTINE xc_rVV10(rho_valence, rho_core, nspin, etxc, vtxc, v, b_value_)
!! Calculate exchange-correlation energy and potential for rVV10.
! Modules to include
! -------------------------------------------------------------------------
use gvect, ONLY : ngm, g
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : omega, tpiba
! -------------------------------------------------------------------------
real(dp), intent(IN) :: rho_valence(:)
!! valence charge density
real(dp), intent(IN) :: rho_core(:)
!! core charge density
INTEGER, INTENT(IN) :: nspin
!! number of spin components
real(dp), intent(inout) :: etxc
!! total XC energy
real(dp), intent(inout) :: vtxc
!! total XC potential
real(dp), intent(inout) :: v(:,:)
!! XC potential on rho grid
real(DP),optional,intent(in) :: b_value_
!
! Local variables
! ----------------------------------------------------------------------------------
!
integer :: i_grid, theta_i, i_proc, I
real(dp) :: grid_cell_volume
real(dp), allocatable :: total_rho(:)
real(dp), allocatable :: gradient_rho(:,:)
real(dp), allocatable :: q0(:)
real(dp), allocatable :: dq0_drho(:)
real(dp), allocatable :: dq0_dgradrho(:)
complex(dp), allocatable :: thetas(:,:)
real(dp) :: Ec_nl
real(dp), allocatable :: potential(:)
logical, save :: first_iteration = .true.
real(dp) :: beta
! ---------------------------------------------------------------------------------------------
! Begin calculations
!call errore('xc_rVV10','rVV10 functional not implemented for spin polarized runs', size(rho_valence,2)-1)
if (nspin>2) call errore('xc_vdW_DF','vdW functional not implemented for nspin > 2', nspin)
if(present(b_value_)) b_value = b_value_
! --------------------------------------------------------------------------------------------------------
call start_clock( 'rVV10' )
beta = 0.0625d0 * (3.0d0 / (b_value**2.0D0) )**(0.75d0)
! Write parameters during the first iteratio
!
if (first_iteration) then
first_iteration = .false.
CALL generate_kernel
if (ionode .and. iverbosity > -1 ) then
WRITE(stdout,'(/ /A )') "---------------------------------------------------------------------------------"
WRITE(stdout,'(A)') "Carrying out rVV10 run using the following parameters:"
WRITE(stdout,'(A,I6,A,I6,A,F8.3)') "Nqs = ",Nqs, " Nr_points = ", Nr_points," r_max = ",r_max
WRITE(stdout, '(A, F8.5, A, F8.5 )') "b_value = ", b_value, " beta = ", beta
WRITE(stdout,'(5X,"q_mesh =",4F12.8)') (q_mesh(I), I=1, 4)
WRITE(stdout,'(13X,4F12.8)') (q_mesh(I), I=5, Nqs)
WRITE(stdout,'(/ A )') "Gradients computed in Reciprocal space"
WRITE(stdout,'(/ A / /)') "---------------------------------------------------------------------------------"
end if
end if
! --------------------------------------------------------------------------------------------------
! Allocate arrays.
! ---------------------------------------------------------------------------------------
allocate( q0(dfftp%nnr) )
allocate( gradient_rho(3,dfftp%nnr) )
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
allocate( total_rho(dfftp%nnr) )
! ---------------------------------------------------------------------------------------
! Add together the valence and core charge densities to get the total charge density
!
total_rho = rho_valence(:) + rho_core(:)
! -------------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT
! -------------------------------------------------------------------------
call fft_gradient_r2r( dfftp, total_rho, g, gradient_rho)
! -------------------------------------------------------------------------
! Get Q and all the derivatives
! -------------------------------------------------------------------------
CALL get_q0_on_grid(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho)
! ---------------------------------------------------------------------------------
allocate( thetas(dfftp%nnr, Nqs) )
CALL get_thetas_on_grid(total_rho, q0, thetas)
call start_clock( 'rVV10_energy')
call vdW_energy(thetas, Ec_nl)
Ec_nl = Ec_nl + beta * SUM(total_rho) * (omega/(dfftp%nr1x*dfftp%nr2x*dfftp%nr3x))
etxc = etxc + Ec_nl
call stop_clock( 'rVV10_energy')
! Print stuff if verbose run
!
if (iverbosity > 0) then
call mp_sum(Ec_nl,intra_bgrp_comm)
if (ionode) write(*,'(/ / A /)') " ----------------------------------------------------------------"
if (ionode) write(*,'(A, F22.15 /)') " Non-local correlation energy = ", Ec_nl
if (ionode) write(*,'(A /)') " ----------------------------------------------------------------"
end if
! ----------------------------------------------------------------------------------------
! Inverse Fourier transform the u_i(k) to get the u_i(r)
!---------------------------------------------------------------------------------------
call start_clock( 'rVV10_ffts')
do theta_i = 1, Nqs
CALL invfft('Rho', thetas(:,theta_i), dfftp)
end do
call stop_clock( 'rVV10_ffts')
! -------------------------------------------------------------------------
call start_clock( 'rVV10_v' )
allocate( potential(dfftp%nnr) )
call get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, gradient_rho, thetas, potential)
! -------------------------------------------------------------------------
! Add beta
! -------------------------------------------------------------------------
potential = potential + beta
v(:,1) = v(:,1) + potential(:)
if (nspin==2) v(:,2) = v(:,2) + potential(:)
call stop_clock( 'rVV10_v' )
! -----------------------------------------------------------------------
! 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 + grid_cell_volume*rho_valence(i_grid)*potential(i_grid)
end do
deallocate(potential)
! ----------------------------------------------------------------------
! Deallocate all arrays.
deallocate(q0, gradient_rho, dq0_drho, dq0_dgradrho, total_rho, thetas)
call stop_clock('rVV10')
END SUBROUTINE xc_rVV10
! #################################################################################################
! | |
! | rVV10_STRESS |
! |_________________|
SUBROUTINE rVV10_stress (rho_valence, rho_core, nspin, sigma)
!! Calculate the stress tensor for rVV10.
USE fft_base, ONLY : dfftp
use gvect, ONLY : ngm, g
USE cell_base, ONLY : tpiba
implicit none
real(dp), intent(IN) :: rho_valence(:)
!! valence charge density
real(dp), intent(IN) :: rho_core(:)
!! core charge density
INTEGER, INTENT(IN) :: nspin
!! number of spin components
real(dp), intent(inout) :: sigma(3,3)
!! stress tensor
real(dp), allocatable :: gradient_rho(:,:) !
real(dp), allocatable :: total_rho(:) ! Rho values
real(dp), allocatable :: q0(:) !
real(dp), allocatable :: dq0_drho(:) ! Q-values
real(dp), allocatable :: dq0_dgradrho(:) !
complex(dp), allocatable :: thetas(:,:) ! Thetas
integer :: i_proc, theta_i, l, m
real(dp) :: sigma_grad(3,3)
real(dp) :: sigma_ker(3,3)
! ---------------------------------------------------------------------------------------------
! Tests
! --------------------------------------------------------------------------------------------------------
if (nspin>2) call errore('rV10_stress',' rVV10 stress not implemented for nspin > 2', nspin)
sigma(:,:) = 0.0_DP
sigma_grad(:,:) = 0.0_DP
sigma_ker(:,:) = 0.0_DP
! ---------------------------------------------------------------------------------------
! Allocations
! ---------------------------------------------------------------------------------------
allocate( gradient_rho(3,dfftp%nnr) )
allocate( total_rho(dfftp%nnr) )
allocate( q0(dfftp%nnr) )
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
allocate( thetas(dfftp%nnr, Nqs) )
! ---------------------------------------------------------------------------------------
! Charge
! ---------------------------------------------------------------------------------------
total_rho = rho_valence(:) + rho_core(:)
! -------------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT
! -------------------------------------------------------------------------
call fft_gradient_r2r( dfftp, total_rho, g, gradient_rho)
! -------------------------------------------------------------------------------------------------------------
! Get q0.
! ---------------------------------------------------------------------------------
CALL get_q0_on_grid(total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho)
! ---------------------------------------------------------------------------------
! Get thetas in reciprocal space.
! ---------------------------------------------------------------------------------
CALL get_thetas_on_grid(total_rho, q0, thetas)
! ---------------------------------------------------------------------------------------
! Stress
! ---------------------------------------------------------------------------------------
CALL rVV10_stress_gradient(total_rho, gradient_rho, q0, dq0_drho, &
dq0_dgradrho, thetas, sigma_grad)
CALL rVV10_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)
enddo
enddo
deallocate( gradient_rho, total_rho, q0, dq0_drho, dq0_dgradrho, thetas )
END SUBROUTINE rVV10_stress
! ###############################################################################################################
! | |
! | rVV10_stress_gradient |
SUBROUTINE rVV10_stress_gradient (total_rho, gradient_rho, q0, dq0_drho, &
dq0_dgradrho, thetas, sigma)
!! Calculate rVV10 stress with gradient correction.
!-----------------------------------------------------------------------------------
! Modules to include
! ----------------------------------------------------------------------------------
use gvect, ONLY : ngm, g, gg, igtongl, &
gl, ngl, gstart
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : omega, tpiba, alat, at, tpiba2
! ----------------------------------------------------------------------------------
implicit none
real(dp), intent(IN) :: total_rho(:) !
real(dp), intent(IN) :: gradient_rho(:, :) ! Input variables
real(dp), intent(inout) :: sigma(:,:) !
real(dp), intent(IN) :: q0(:) !
real(dp), intent(IN) :: dq0_drho(:) !
real(dp), intent(IN) :: dq0_dgradrho(:) !
complex(dp), intent(IN) :: thetas(:,:) !
complex(dp), allocatable :: u_vdW(:,:) !
real(dp), allocatable :: d2y_dx2(:,:) !
real(dp) :: y(Nqs), dP_dq0, P, a, b, c, d, e, f ! Interpolation
real(dp) :: dq !
integer :: q_low, q_hi, q, q1_i, q2_i , g_i ! Loop and q-points
integer :: l, m
real(dp) :: prefactor ! Final summation of sigma
integer :: i_proc, theta_i, i_grid, q_i, & !
ix, iy, iz ! Iterators
character(LEN=1) :: intvar
real(dp) :: const
!real(dp) :: at_inverse(3,3)
allocate( d2y_dx2(Nqs, Nqs) )
allocate( u_vdW(dfftp%nnr, Nqs) )
const = 1.0D0 / (3.0D0 * b_value**(3.0D0/2.0D0) * pi**(5.0D0/4.0D0) )
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.
! ---------------------------------------------------------------------------------------------------
call start_clock( 'rVV10_ffts')
do theta_i = 1, Nqs
CALL invfft('Rho', u_vdW(:,theta_i), dfftp)
end do
call stop_clock( 'rVV10_ffts')
! --------------------------------------------------------------------------------------------------
! Get the second derivatives for interpolating the P_i
! ---------------------------------------------------------------------------------------------------
call initialize_spline_interpolation(q_mesh, d2y_dx2(:,:))
! ---------------------------------------------------------------------------------------------
i_grid = 0
! ----------------------------------------------------------------------------------------------------
! Do the real space integration to obtain the stress component
! ----------------------------------------------------------------------------------------------------
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('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)
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (total_rho(i_grid) > epsr) then
prefactor = u_vdW(i_grid,q_i) * const * total_rho(i_grid)**(3.0D0/4.0D0) * dP_dq0 * dq0_dgradrho(i_grid)
do l = 1, 3
do m = 1, l
sigma (l, m) = sigma (l, m) - prefactor * &
(gradient_rho(l,i_grid) * gradient_rho(m,i_grid))
enddo
enddo
endif
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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 rVV10_stress_gradient
! ###############################################################################################################
! | |
! | rVV10_stress_kernel |
! | |
SUBROUTINE rVV10_stress_kernel (total_rho, q0, thetas, sigma)
! Modules to include
! ----------------------------------------------------------------------------------
use gvect, ONLY : ngm, g, gg, igtongl, gl, ngl, gstart
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : omega, tpiba, tpiba2
USE constants, ONLY: pi
implicit none
real(dp), intent(IN) :: q0(:)
real(dp), intent(IN) :: total_rho(:)
real(dp), intent(inout) :: sigma(3,3) !
complex(dp), intent(IN) :: thetas(:,:)
real(dp), allocatable :: dkernel_of_dk(:,:) !
integer :: l, m, q1_i, q2_i , g_i !
real(dp) :: g2, ngmod2, g_kernel, G_multiplier !
integer :: last_g, theta_i
allocate( dkernel_of_dk(Nqs, Nqs) )
sigma(:,:) = 0.0_DP
! --------------------------------------------------------------------------------------------------
! Integration in g-space
! ---------------------------------------------------------------------------------------------------
last_g = -1
G_multiplier = 1.0D0
if (gamma_only) G_multiplier = 2.0D0
do g_i = gstart, ngm
g2 = gg (g_i) * tpiba2
g_kernel = sqrt(g2)
if ( igtongl(g_i) .ne. last_g) then
call interpolate_Dkernel_Dk(g_kernel, dkernel_of_dk) ! Gets the derivatives
last_g = igtongl(g_i)
end if
do q2_i = 1, Nqs
do q1_i = 1, Nqs
do l = 1, 3
do m = 1, l
sigma (l, m) = sigma (l, m) - G_multiplier * 0.5 * &
thetas(dfftp%nl(g_i),q1_i)*dkernel_of_dk(q1_i,q2_i)*conjg(thetas(dfftp%nl(g_i),q2_i))* &
(g (l, g_i) * g (m, g_i) * tpiba2) / g_kernel
end do
end do
enddo
end do
if (g_i < gstart ) sigma(:,:) = sigma(:,:) / G_multiplier
enddo
call mp_sum( sigma, intra_bgrp_comm )
deallocate( dkernel_of_dk )
END SUBROUTINE rVV10_stress_kernel
! ###############################################################################################################
! | |
! | GET_Q0_ON_GRID |
! |__________________|
SUBROUTINE get_q0_on_grid (total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho)
USE fft_base, ONLY : dfftp
real(dp), intent(IN) :: total_rho(:), gradient_rho(:,:)
real(dp), intent(OUT) :: q0(:), dq0_drho(:), dq0_dgradrho(:)
integer, parameter :: m_cut = 12
real(dp) :: dw0_dn, dk_dn, gmod2
real(dp) :: mod_grad, wp2, wg2, w0, k
real(dp) :: q, exponent, dq0_dq
integer :: i_grid, index, count=0
! initialize q0-related arrays ...
q0(:) = q_cut
dq0_drho(:) = 0.0_DP
dq0_dgradrho(:) = 0.0_DP
do i_grid = 1, dfftp%nnr
if (total_rho(i_grid) > epsr) then
gmod2 = gradient_rho(1,i_grid)**2 + &
gradient_rho(2,i_grid)**2 + &
gradient_rho(3,i_grid)**2
! Calculate some intermediate values needed to find q
! ------------------------------------------------------------------------------------
mod_grad = sqrt(gmod2)
wp2= 16.0_dp*pi*total_rho(i_grid)
wg2 = 4.0_dp*C_value * (mod_grad/total_rho(i_grid))**4
k = b_value*3.0_dp*pi* ((total_rho(i_grid)/(9.0_dp*pi))**(1.0_dp/6.0_dp))
w0 = sqrt( wg2 + wp2/3.0_dp )
q = w0 / k
! Here, we calculate q0 by saturating q according
! ---------------------------------------------------------------------------------------
exponent = 0.0_dp
dq0_dq = 0.0_dp
do index = 1, m_cut
exponent = exponent + ( (q/q_cut)**index)/index
dq0_dq = dq0_dq + ( (q/q_cut)**(index-1))
end do
q0(i_grid) = q_cut*(1.0_dp - exp(-exponent))
dq0_dq = dq0_dq * exp(-exponent)
! ---------------------------------------------------------------------------------------
if (q0(i_grid) < q_min) then
q0(i_grid) = q_min
end if
!---------------------------------Final values---------------------------------
dw0_dn = 1.0_dp/(2.0_dp*w0) * (16.0_dp/3.0_dp*pi - 4.0_dp*wg2 / total_rho(i_grid) )
dk_dn = k / ( 6.0_dp * total_rho(i_grid) )
dq0_drho(i_grid) = dq0_dq / (k**2) * (dw0_dn * k - dk_dn * w0 )
IF ( gmod2 > epsr) THEN
dq0_dgradrho(i_grid) = dq0_dq / ( 2.0_dp*k*w0 ) * 4.0_dp*wg2 / (mod_grad**2)
ELSE
dq0_dgradrho(i_grid) = 0.0_dp
ENDIF
endif
end do
end SUBROUTINE get_q0_on_grid
! ###############################################################################################################
! | |
! | GET_THETAS_ON_GRID |
SUBROUTINE get_thetas_on_grid (total_rho, q0_on_grid, thetas)
real(dp), intent(in) :: total_rho(:), q0_on_grid(:)
complex(dp), intent(inout):: thetas(:,:)
integer :: i_grid, Ngrid_points
integer :: theta_i
Ngrid_points = size(q0_on_grid)
! Interpolate the P_i polynomials
CALL spline_interpolation(q_mesh, q0_on_grid, thetas)
! Form the thetas where theta is defined as rho*p_i(q0)
! ------------------------------------------------------------------------------------
do i_grid = 1, Ngrid_points
if (total_rho(i_grid) > epsr ) then
thetas(i_grid,:) = thetas(i_grid,:) * (1.0 / (3.0 * sqrt(pi) &
* ( b_value**(3.0/2.0) ) ) ) * (total_rho(i_grid) / pi)**(3.0/4.0)
else
thetas(i_grid,:) = 0.0d0
endif
end do
! ------------------------------------------------------------------------------------
call start_clock( 'rVV10_ffts')
do theta_i = 1, Nqs
CALL fwfft ('Rho', thetas(:,theta_i), dfftp)
end do
call stop_clock( 'rVV10_ffts')
END SUBROUTINE get_thetas_on_grid
! ###############################################################################################################
! | |
! | SPLINE_INTERPOLATION |
! |________________________|
SUBROUTINE spline_interpolation (x, evaluation_points, values)
real(dp), intent(in) :: x(:)
!! The x values used to form the interpolation
real(dp), intent(in) :: evaluation_points(:)
!! (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, index, P_i ! Some indexing variables
real(dp), allocatable :: y(:) ! Temporary variables needed for the interpolation
real(dp) :: a, b, c, d, dx !
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 )
index = (upper_bound+lower_bound)/2
if ( evaluation_points(i_grid) > x(index) ) then
lower_bound = index
else
upper_bound = index
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 in C" by Cambridge
!! University Press, pages 96-97. It was adapted for Fortran and for the problem at hand.
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, index ! 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 index = 2, Nx-1
temp1 = (x(index)-x(index-1))/(x(index+1)-x(index-1))
temp2 = temp1 * d2y_dx2(P_i,index-1) + 2.0D0
d2y_dx2(P_i,index) = (temp1-1.0D0)/temp2
temp_array(index) = (y(index+1)-y(index))/(x(index+1)-x(index)) &
- (y(index)-y(index-1))/(x(index)-x(index-1))
temp_array(index) = (6.0D0*temp_array(index)/(x(index+1)-x(index-1)) &
- temp1*temp_array(index-1))/temp2
end do
d2y_dx2(P_i,Nx) = 0.0D0
do index=Nx-1, 1, -1
d2y_dx2(P_i,index) = d2y_dx2(P_i,index) * d2y_dx2(P_i,index+1) + temp_array(index)
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 in C" by Cambridge
!! University Press, page 97. Adapted for Fortran and the problem at hand. This function is used to
!! find the Phi_alpha_beta needed for equations 11 and 14 of SOLER.
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
! ###############################################################################################################
! | |
! | INTERPOLATE_DKERNEL_DK |
! |________________________|
subroutine interpolate_Dkernel_Dk(k, dkernel_of_dk)
implicit none
real(dp), intent(in) :: k
real(dp), intent(inout) :: dkernel_of_dk(Nqs,Nqs)
integer :: q1_i, q2_i, k_i
real(dp) :: A, B, dAdk, dBdk, dCdk, dDdk
! -------------------------------------------------------------------------------------
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 fft_base, ONLY : dfftp
USE cell_base, ONLY : tpiba, omega
complex(dp), intent(in) :: thetas(:,:)
complex(dp), intent(out) :: u_vdW(:,:)
real(dp), allocatable :: kernel_of_k(:,:)
real(dp) :: g
integer :: last_g, g_i, q1_i, q2_i, count, i_grid
complex(dp) :: theta(Nqs)
! -------------------------------------------------------------------------------------------------
allocate( kernel_of_k(Nqs, Nqs) )
u_vdW(:,:) = CMPLX(0.0_DP,0.0_DP, kind=dp)
last_g = -1
do g_i = 1, ngm
if ( igtongl(g_i) .ne. last_g) then
g = sqrt(gl(igtongl(g_i))) * tpiba
call interpolate_kernel(g, kernel_of_k)
last_g = igtongl(g_i)
end if
theta = thetas(dfftp%nl(g_i),:)
do q2_i = 1, Nqs
do q1_i = 1, Nqs
u_vdW(dfftp%nl(g_i),q2_i) = u_vdW(dfftp%nl(g_i),q2_i) + kernel_of_k(q2_i,q1_i)*theta(q1_i)
end do
end do
end do
if (gamma_only) u_vdW(dfftp%nlm(:),:) = CONJG(u_vdW(dfftp%nl(:),:))
deallocate( kernel_of_k )
! -----------------------------------------------------------------------------------------------
end subroutine thetas_to_uk
! #################################################################################################
! | |
! | VDW_ENERGY |
! |_____________|
subroutine vdW_energy(thetas, vdW_xc_energy)
USE gvect, ONLY : gg, ngm, igtongl, gl, ngl, gstart
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : tpiba, omega
complex(dp), intent(inout) :: thetas(:,:)
real(dp), intent(out) :: vdW_xc_energy
real(dp), allocatable :: kernel_of_k(:,:)
real(dp) :: g
integer :: last_g
integer :: g_i, q1_i, q2_i, count, i_grid
complex(dp) :: theta(Nqs), thetam(Nqs), theta_g(Nqs)
real(dp) :: G0_term, G_multiplier
complex(dp), allocatable :: u_vdw(:,:)
vdW_xc_energy = 0.0D0
allocate (u_vdW(dfftp%nnr,Nqs))
u_vdW(:,:) = CMPLX(0.0_DP,0.0_DP, kind=dp)
allocate( kernel_of_k(Nqs, Nqs) )
!
! Here we should use gstart,ngm but all the cases are handeld 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(:),:))
! Final value
vdW_xc_energy = 0.5D0 * omega * vdW_xc_energy
deallocate( kernel_of_k )
thetas(:,:) = u_vdW(:,:)
deallocate (u_vdW)
! ---------------------------------------------------------------------------------------------------
end subroutine vdW_energy
! ###############################################################################################################
! | |
! | GET_POTENTIAL |
! |_________________|
subroutine get_potential(q0, dq0_drho, dq0_dgradrho, total_rho, gradient_rho, u_vdW, potential)
use gvect, ONLY : g
USE fft_base, ONLY : dfftp
USE cell_base, ONLY : alat, tpiba
real(dp), intent(in) :: q0(:), gradient_rho(:,:)
real(dp), intent(in) :: dq0_drho(:), dq0_dgradrho(:)
real(dp), intent(in) :: total_rho(:)
complex(dp), intent(in) :: u_vdW(:,:)
real(dp), intent(inout) :: potential(:)
real(dp), allocatable, save :: d2y_dx2(:,:)
integer :: i_grid, P_i,icar
integer :: q_low, q_hi, q
real(dp) :: dq, a, b, c, d, e, f
real(dp) :: y(Nqs), dP_dq0, P
!
real(dp), allocatable ::h_prefactor(:)
complex(dp), allocatable ::h(:)
real(dp) :: dtheta_dn, dtheta_dgradn
real(dp) :: const
allocate (h_prefactor(dfftp%nnr),h(dfftp%nnr))
const = 1.0D0 / (3.0D0 * b_value**(3.0D0/2.0D0) * pi**(5.0D0/4.0D0) )
potential = 0.0D0
h_prefactor = 0.0D0
! -------------------------------------------------------------------------------------------
! Get the second derivatives of the P_i functions for interpolation
! ---------------------------------------------------------------------------------------------
if (.not. allocated( d2y_dx2) ) then
allocate( d2y_dx2(Nqs, Nqs) )
call initialize_spline_interpolation(q_mesh, d2y_dx2(:,:))
end if
! ---------------------------------------------------------------------------------------------
do i_grid = 1,dfftp%nnr
q_low = 1
q_hi = Nqs
! Figure out which bin our value of q0 is in in the q_mesh
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
do while ( (q_hi - q_low) > 1)
q = int((q_hi + q_low)/2)
if (q_mesh(q) > q0(i_grid)) then
q_hi = q
else
q_low = q
end if
end do
if (q_hi == q_low) call errore('get_potential','qhi == qlow',1)
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
dq = q_mesh(q_hi) - q_mesh(q_low)
a = (q_mesh(q_hi) - q0(i_grid))/dq
b = (q0(i_grid) - q_mesh(q_low))/dq
c = (a**3 - a)*dq**2/6.0D0
d = (b**3 - b)*dq**2/6.0D0
e = (3.0D0*a**2 - 1.0D0)*dq/6.0D0
f = (3.0D0*b**2 - 1.0D0)*dq/6.0D0
do P_i = 1, Nqs
y = 0.0D0
y(P_i) = 1.0D0
dP_dq0 = (y(q_hi) - y(q_low))/dq - e*d2y_dx2(P_i,q_low) + f*d2y_dx2(P_i,q_hi)
P = a*y(q_low) + b*y(q_hi) + c*d2y_dx2(P_i,q_low) + d*d2y_dx2(P_i,q_hi)
! IF THE CHARGE DENSITY IS NEGATIVE WE PUT POTENTIAL = 0, OUTSIDE THE SUBROUTINE WE ADD BETA.
if (total_rho(i_grid) > epsr) then
dtheta_dn = const * (3.0D0/4.0D0) / (total_rho(i_grid)**(1.0D0/4.0D0)) * P + &
const * total_rho(i_grid)**(3.0D0/4.0D0) * dP_dq0 * dq0_drho(i_grid)
dtheta_dgradn = const * total_rho(i_grid)**(3.0D0/4.0D0) * dP_dq0 * dq0_dgradrho(i_grid)
potential(i_grid) = potential(i_grid) + u_vdW(i_grid,P_i)* dtheta_dn
if (q0(i_grid) .ne. q_mesh(Nqs)) then
h_prefactor(i_grid) = h_prefactor(i_grid) + u_vdW(i_grid,P_i)* dtheta_dgradn
end if
end if
end do
end do
do icar = 1,3
h(:) = CMPLX( h_prefactor(:)*gradient_rho(icar,:), 0.0_DP, kind=dp)
CALL fwfft ('Rho', h, dfftp)
h(dfftp%nl(:)) = CMPLX(0.0_DP,1.0_DP,kind=dp)*tpiba*g(icar,:)*h(dfftp%nl(:))
if (gamma_only) h(dfftp%nlm(:)) = CONJG(h(dfftp%nl(:)))
CALL invfft ('Rho', h, dfftp)
potential(:) = potential(:) - REAL(h(:))
end do
! ------------------------------------------------------------------------------------------------------------------------
deallocate (h_prefactor,h)
end subroutine get_potential
! ###############################################################################################################
! | |
! | generate_kernel |
! |___________________|
SUBROUTINE generate_kernel
implicit none
integer :: q1_i, q2_i, r_i ! Indexing variables
real(dp) :: d1, d2 ! Intermediate values
kernel = 0.0D0
d2phi_dk2 = 0.0D0
do q1_i = 1, Nqs
do q2_i = 1, q1_i
do r_i = 1, Nr_points
d1 = q_mesh(q1_i) * (dr * r_i)**2 ! Different definition of d1 and d2 for vv10
d2 = q_mesh(q2_i) * (dr * r_i)**2 ! Different definition of d1 and d2 for vv10
kernel(r_i, q1_i, q2_i) = -24.0D0 / ( ( d1+1.0 ) * ( d2+1.0 ) * ( d1+d2+2.0 ) )
end do
call radial_fft( kernel(:, q1_i, q2_i) )
call set_up_splines( kernel(:, q1_i, q2_i), d2phi_dk2(:, q1_i, q2_i) )
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 do
END SUBROUTINE generate_kernel
! ###############################################################################################################
! | |
! | radial_fft |
! |______________|
SUBROUTINE radial_fft(phi)
REAL(DP), INTENT(INOUT) :: phi(0:Nr_points)
REAL(DP) :: phi_k(0:Nr_points)
INTEGER :: k_i, r_i
REAL(DP) :: r, k
phi_k = 0.0D0
DO r_i = 1, Nr_points
r = r_i * dr
phi_k(0) = phi_k(0) + phi(r_i)*r**2
END DO
phi_k(0) = phi_k(0) - 0.5D0 * (Nr_points*dr)**2 * phi(Nr_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
phi = 4.0D0 * pi * phi_k * dr
END SUBROUTINE radial_fft
! ###############################################################################################################
! | |
! | set_up_splines |
! |__________________|
SUBROUTINE set_up_splines(phi, D2)
REAL(DP), INTENT(IN) :: phi(0:Nr_points)
REAL(DP), INTENT(INOUT) :: D2(0:Nr_points)
REAL(DP), ALLOCATABLE :: temp_array(:)
REAL(DP) :: temp_1, temp_2
INTEGER :: r_i
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
END MODULE rVV10
Reference in New Issue View Git Blame Copy Permalink