mirror of https://gitlab.com/QEF/q-e.git
1472 lines
52 KiB
Fortran
1472 lines
52 KiB
Fortran
!
|
|
! 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
|
|
|
|
USE kinds, ONLY : dp
|
|
USE constants, ONLY : pi, e2
|
|
USE kernel_table, ONLY : q_mesh, Nr_points, Nqs, r_max
|
|
USE mp, ONLY : mp_bcast, mp_sum, mp_barrier
|
|
USE mp_bands, ONLY : intra_bgrp_comm
|
|
USE io_global, ONLY : ionode
|
|
USE fft_base, ONLY : dfftp
|
|
USE fft_interfaces, ONLY : fwfft, invfft
|
|
USE control_flags, ONLY : gamma_only, iverbosity
|
|
USE io_global, ONLY : stdout
|
|
|
|
IMPLICIT NONE
|
|
|
|
real(dp), parameter :: epsr = 1.d-12, epsg = 1.D-10
|
|
|
|
real(dp) :: b_value = 6.3_DP
|
|
real(dp) :: C_value = 0.0093
|
|
|
|
private
|
|
public :: xc_rVV10, &
|
|
interpolate_kernel, &
|
|
initialize_spline_interpolation, &
|
|
stress_rVV10, b_value
|
|
|
|
CONTAINS
|
|
|
|
!! #################################################################################################
|
|
!! | |
|
|
!! | xc_rVV10 |
|
|
!! |_____________|
|
|
|
|
SUBROUTINE xc_rVV10(rho_valence, rho_core, nspin, etxc, vtxc, v, b_value_)
|
|
|
|
!! Modules to include
|
|
!! -------------------------------------------------------------------------
|
|
|
|
use gvect, ONLY : ngm, g
|
|
USE fft_base, ONLY : dfftp
|
|
USE cell_base, ONLY : omega, tpiba
|
|
!! -------------------------------------------------------------------------
|
|
|
|
!! Local variables
|
|
!! ----------------------------------------------------------------------------------
|
|
! _
|
|
real(dp), intent(IN) :: rho_valence(:) !
|
|
real(dp), intent(IN) :: rho_core(:) ! PWSCF input variables
|
|
INTEGER, INTENT(IN) :: nspin !
|
|
real(dp), intent(inout) :: etxc, vtxc, v(:,:) !_
|
|
real(DP),optional,intent(in) :: b_value_
|
|
|
|
|
|
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.
|
|
|
|
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 > 1) 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
|
|
|
|
|
|
!! #################################################################################################
|
|
!! | |
|
|
!! | STRESS_rVV10 |
|
|
!! |_________________|
|
|
|
|
SUBROUTINE stress_rVV10(rho_valence, rho_core, nspin, sigma)
|
|
|
|
USE fft_base, ONLY : dfftp
|
|
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 :: 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
|
|
!! --------------------------------------------------------------------------------------------------------
|
|
|
|
!call errore('stress_rVV10','vdW 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)
|
|
|
|
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 stress_rVV10_gradient(total_rho, gradient_rho, q0, dq0_drho, &
|
|
dq0_dgradrho, thetas, sigma_grad)
|
|
|
|
CALL stress_rVV10_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 stress_rVV10
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | stress_rVV10_gradient |
|
|
|
|
SUBROUTINE stress_rVV10_gradient (total_rho, gradient_rho, q0, dq0_drho, &
|
|
dq0_dgradrho, 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, 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 stress_rVV10_gradient
|
|
|
|
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | stress_rVV10_kernel |
|
|
!! | |
|
|
|
|
SUBROUTINE stress_rVV10_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 stress_rVV10_kernel
|
|
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | GET_Q0_ON_GRID |
|
|
!! |__________________|
|
|
|
|
SUBROUTINE get_q0_on_grid (total_rho, gradient_rho, q0, dq0_drho, dq0_dgradrho)
|
|
|
|
USE fft_base, ONLY : dfftp
|
|
USE kernel_table, ONLY : q_cut, q_min
|
|
|
|
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(:), 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, 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 |
|
|
!! |___________________________________|
|
|
|
|
|
|
!! This routine is modeled after an algorithm from "Numerical Recipes in C" by Cambridge
|
|
!! University Press, pages 96-97. It was adapted for Fortran and for the problem at hand.
|
|
|
|
SUBROUTINE initialize_spline_interpolation (x, d2y_dx2)
|
|
|
|
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 |
|
|
!! |____________________|
|
|
|
|
|
|
!! 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.
|
|
|
|
|
|
subroutine interpolate_kernel(k, kernel_of_k)
|
|
|
|
USE kernel_table, ONLY : r_max, Nr_points, kernel, d2phi_dk2, dk
|
|
|
|
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)
|
|
|
|
USE kernel_table, ONLY : r_max, Nr_points, kernel, d2phi_dk2, 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
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | GRADIENT_COEFFICIENTS |
|
|
!! |_________________________|
|
|
|
|
|
|
!! This routine returns a pointer to an array holding the coefficients for a derivative expansion to some order.
|
|
!! The derivative is found by multiplying the value of the function at a point + or - n away from the sample point by
|
|
!! the coefficient gradient_coefficients(+ or - n) and dividing by the appropriate dx for that direction.
|
|
|
|
|
|
function gradient_coefficients(N)
|
|
|
|
real(dp), allocatable, target, save:: coefficients(:) !! The local array that will hold the coefficients. A pointer to this
|
|
! !! array will be returned by the function
|
|
|
|
integer, intent(in), optional :: N !! The number of neighbors to use on each side for the gradient
|
|
! !! calculation. Can be between 1 (i.e. 3 point derivative formula)
|
|
! !! and 6 (i.e. 13 point derivative formula).
|
|
|
|
real(dp), pointer :: gradient_coefficients(:) !! Pointer to the coefficients array that will be returned
|
|
|
|
|
|
if (.not. allocated(coefficients) ) then
|
|
|
|
if (.not. present(N) ) call errore('gradient_coefficients', 'Number of neighbors for gradient must be specified',2)
|
|
|
|
allocate( coefficients(-N:N) )
|
|
|
|
select case (N)
|
|
|
|
case (1)
|
|
coefficients(-1:1) = (/-0.5D0, 0.0D0, 0.5D0/)
|
|
case (2)
|
|
coefficients(-2:2) = (/0.0833333333333333D0, -0.6666666666666666D0, 0.0D0, &
|
|
0.6666666666666666D0, -0.0833333333333333D0/)
|
|
case (3)
|
|
coefficients(-3:3) = (/-0.0166666666666666D0, 0.15D0, -0.75D0, 0.0D0, 0.75D0, &
|
|
-0.15D0, 0.016666666666666666D0/)
|
|
case (4)
|
|
coefficients(-4:4) = (/0.00357142857143D0, -0.03809523809524D0, 0.2D0, -0.8D0, 0.0D0, &
|
|
0.8D0, -0.2D0, 0.03809523809524D0, -0.00357142857143D0/)
|
|
case (5)
|
|
coefficients(-5:5) = (/-0.00079365079365D0, 0.00992063492063D0, -0.05952380952381D0, &
|
|
0.23809523809524D0, -0.8333333333333333D0, 0.0D0, 0.8333333333333333D0, &
|
|
-0.23809523809524D0, 0.05952380952381D0, -0.00992063492063D0, 0.00079365079365D0/)
|
|
case (6)
|
|
coefficients(-6:6) = (/0.00018037518038D0, -0.00259740259740D0, 0.01785714285714D0, &
|
|
-0.07936507936508D0, 0.26785714285714D0, -0.85714285714286D0, 0.0D0, &
|
|
0.85714285714286D0, -0.26785714285714D0, 0.07936507936508D0, &
|
|
-0.01785714285714D0, 0.00259740259740D0, -0.00018037518038D0/)
|
|
case default
|
|
|
|
call errore('xc_vdW_DF', 'Order of numerical gradient not implemented', 2)
|
|
|
|
end select
|
|
|
|
end if
|
|
|
|
gradient_coefficients => coefficients
|
|
|
|
|
|
end function gradient_coefficients
|
|
|
|
|
|
!! ###############################################################################################################
|
|
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | GET_3D_INDICES |
|
|
!! |__________________|
|
|
|
|
!! This routine builds a rank 3 array that holds the indices into the FFT grid for a point with a given
|
|
!! set of x, y, and z indices. The array holds an extra 2N points in each dimension (N to the left and N
|
|
!! to the right) so the code can find the neighbors of edge points easily. This is done by just copying the
|
|
!! first N points in each dimension to the end of that dimension and the end N points to the beginning.
|
|
|
|
|
|
function get_3d_indices(N)
|
|
|
|
USE fft_base, ONLY : dfftp
|
|
|
|
|
|
integer, intent(in), optional :: N !! The number of neighbors in each direction that will
|
|
! !! be used for the gradient formula. If not supplied,
|
|
! !! the code just returns the pointer to the already
|
|
! !! allocated rho_3d array.
|
|
|
|
real(dp) :: dx, dy, dz !!
|
|
integer :: ix1, ix2, ix3, i_grid !! Index variables
|
|
|
|
integer, allocatable, target, save :: rho_3d(:,:,:) !! The local array that will store the indices. Only a pointer
|
|
! !! to this array will be returned.
|
|
|
|
integer, pointer :: get_3d_indices(:,:,:) !! The returned pointer to the rho_3d array of indices.
|
|
|
|
|
|
|
|
!! If the routine has not already been run we set up the rho_3d array by looping over it
|
|
!! and assigning indices to its elements. If this routine has already been run we simply
|
|
!! return a pointer to the existing array.
|
|
!! --------------------------------------------------------------------------------
|
|
|
|
if (.not. allocated(rho_3d)) then
|
|
|
|
! Check to make sure we have been given the number of neighbors since the routine has
|
|
! not been run yet.
|
|
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
if (.not. present(N)) then
|
|
|
|
call errore('get_3d_rho','Number of neighbors for numerical derivatives &
|
|
& must be specified',2)
|
|
|
|
end if
|
|
|
|
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
allocate( rho_3d(-N+1:dfftp%nr1x+N, -N+1:dfftp%nr2x+N, -N+1:dfftp%nr3x+N) )
|
|
|
|
i_grid = 0
|
|
|
|
do ix3 = 1, dfftp%nr3x
|
|
do ix2 = 1, dfftp%nr2x
|
|
do ix1 = 1, dfftp%nr1x
|
|
|
|
i_grid = i_grid + 1
|
|
|
|
rho_3d(ix1, ix2, ix3) = i_grid
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
|
|
! Apply periodic boundary conditions to extend the array by N places in each
|
|
! direction
|
|
! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
rho_3d(-N+1:0,:,:) = rho_3d(dfftp%nr1x-N+1:dfftp%nr1x, :, :)
|
|
rho_3d(:,-N+1:0,:) = rho_3d(:, dfftp%nr2x-N+1:dfftp%nr2x, :)
|
|
rho_3d(:,:,-N+1:0) = rho_3d(:, :, dfftp%nr3x-N+1:dfftp%nr3x)
|
|
|
|
rho_3d(dfftp%nr1x+1:dfftp%nr1x+N, :, :) = rho_3d(1:N, :, :)
|
|
rho_3d(:, dfftp%nr2x+1:dfftp%nr2x+N, :) = rho_3d(:, 1:N, :)
|
|
rho_3d(:, :, dfftp%nr3x+1:dfftp%nr3x+N) = rho_3d(:, :, 1:N)
|
|
|
|
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
|
|
end if
|
|
|
|
!! ------------------------------------------------------------------------------------------
|
|
|
|
|
|
!! Return the point to rho_3d
|
|
get_3d_indices => rho_3d
|
|
|
|
|
|
end function get_3d_indices
|
|
|
|
|
|
!! ###############################################################################################################
|
|
!! | |
|
|
!! | INVERT_3X3_MATRIX |
|
|
!! |_____________________|
|
|
|
|
!! This routine is just a hard-wired subroutine to invert a 3x3 matrix. It is used to invert the matrix of
|
|
!! unit cell basis vectors to find the gradient and the derivative of the gradient with respect to the
|
|
!! density.
|
|
|
|
subroutine invert_3x3_matrix(M)
|
|
|
|
real(dp), intent(inout) :: M(3,3) !! On input, the 3x3 matrix to be inverted
|
|
! !! On output, the inverse of the 3x3 matrix given
|
|
|
|
real(dp) :: temp(3,3) !! Temporary storage
|
|
|
|
real(dp) :: determinant_M !! The determinant of the input 3x3 matrix
|
|
|
|
|
|
temp = 0.0D0
|
|
|
|
temp(1,1) = M(2,2)*M(3,3) - M(2,3)*M(3,2)
|
|
temp(1,2) = M(1,3)*M(3,2) - M(1,2)*M(3,3)
|
|
temp(1,3) = M(1,2)*M(2,3) - M(1,3)*M(2,2)
|
|
temp(2,1) = M(2,3)*M(3,1) - M(2,1)*M(3,3)
|
|
temp(2,2) = M(1,1)*M(3,3) - M(1,3)*M(3,1)
|
|
temp(2,3) = M(1,3)*M(2,1) - M(1,1)*M(2,3)
|
|
temp(3,1) = M(2,1)*M(3,2) - M(2,2)*M(3,1)
|
|
temp(3,2) = M(1,2)*M(3,1) - M(1,1)*M(3,2)
|
|
temp(3,3) = M(1,1)*M(2,2) - M(1,2)*M(2,1)
|
|
|
|
determinant_M = M(1,1) * (M(2,2)*M(3,3) - M(2,3)*M(3,2)) &
|
|
- M(1,2) * (M(2,1)*M(3,3) - M(2,3)*M(3,1)) &
|
|
+ M(1,3) * (M(2,1)*M(3,2) - M(2,2)*M(3,1))
|
|
|
|
if (abs(determinant_M) > 1e-6) then
|
|
|
|
M = 1.0D0/determinant_M*temp
|
|
|
|
else
|
|
|
|
call errore('invert_3x3_matrix','Matrix is close to singular',1)
|
|
|
|
end if
|
|
|
|
end subroutine invert_3x3_matrix
|
|
|
|
END MODULE rVV10
|