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

"numerical_gradient" in vdW-DF and rVV10 replaced by "fft_gradent_r2r". 

NOTA BENE: the indices of the array containing the gradient for vdW-DF and
rVV10 functionals have been reverted: now the index of components is the first,
the index  of grid points is the last last. Not sure this is the best choice
but it is the choice done (almost) everywhere else in QE
This commit is contained in:
Paolo Giannozzi 2018-01-18 09:50:35 +01:00
parent 6e01664329
commit a88e82dd2a
4 changed files with 64 additions and 184 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
LR_Modules/dv_rVV10.f90
View File

@ -20,7 +20,7 @@ MODULE ph_rVV10
USE control_flags, ONLY : iverbosity, gamma_only
USE io_global, ONLY : stdout
USE rVV10, ONLY : b_value, initialize_spline_interpolation, &
numerical_gradient, interpolate_kernel
interpolate_kernel
USE gc_lr, ONLY : grho
@ -145,7 +145,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
!! Global variables
allocate(total_rho(dfftp%nnr) )
allocate(gradient_rho(dfftp%nnr, 3))
allocate(gradient_rho(3,dfftp%nnr))
allocate(gradient_drho(dfftp%nnr, 3))
allocate(q0(dfftp%nnr), q(dfftp%nnr))
allocate(dq0_dq(dfftp%nnr), d2q0_dq2(dfftp%nnr))
@ -180,7 +180,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
!! -------------------------------------------------------------------------
total_rho(:) = rho(:,1)
call numerical_gradient(total_rho,gradient_rho)
call fft_gradient_r2r( dfftp, total_rho, g, gradient_rho )
CALL qgradient (q_point, dfftp%nnr, drho(:,1), ngm, g, dfftp%nl, alat, gradient_drho)
@ -203,7 +203,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
do i_grid = 1,dfftp%nnr
gmod2 = gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2
gmod2 = gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2
if (total_rho(i_grid) <= epsr ) cycle
CALL get_abcdef (q0, i_grid, q_hi, q_low, dq, a,b,c,d,e,f )
@ -271,7 +271,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
do i_grid = 1,dfftp%nnr
gmod2 = gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2
gmod2 = gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2
if (total_rho(i_grid) <= epsr) cycle
CALL get_abcdef (q0, i_grid, q_hi, q_low, dq, a,b,c,d,e,f )
@ -299,7 +299,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
allocate (delta_h2_aux(dfftp%nnr))
do icar = 1,3
delta_h(:) = (h1t(:) * gradient_rho(:,icar)+ h2t(:) * gradient_drho(:,icar))
delta_h(:) = (h1t(:) * gradient_rho(icar,:)+ h2t(:) * gradient_drho(:,icar))
CALL fwfft ('Rho', delta_h, dfftp)
@ -442,10 +442,10 @@ end subroutine get_delta_v
!!
!! Fractions
!!
gmod = sqrt(gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2)
gradn_graddeltan = gradient_rho(i_grid,1)*gradient_drho(i_grid,1) + &
gradient_rho(i_grid,2)*gradient_drho(i_grid,2) + &
gradient_rho(i_grid,3)*gradient_drho(i_grid,3)
gmod = sqrt(gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2)
gradn_graddeltan = gradient_rho(1,i_grid)*gradient_drho(i_grid,1) + &
gradient_rho(2,i_grid)*gradient_drho(i_grid,2) + &
gradient_rho(3,i_grid)*gradient_drho(i_grid,3)
END SUBROUTINE get_thetas_exentended
@ -516,7 +516,7 @@ end subroutine get_delta_v
if (total_rho(i_grid) < epsr)cycle
gmod2 = gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2
gmod2 = gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2
gmod = sqrt(gmod2)
@ -750,6 +750,8 @@ subroutine qgradient (xq, nrxx, a, ngm, g, nl, alat, ga)
end subroutine qgradient
! ####################################################################
! | |
! | thetas_to_uk |
END MODULE ph_rVV10

26
LR_Modules/dv_vdW_DF.f90
View File

@ -20,7 +20,7 @@ MODULE ph_vdW_DF
USE control_flags, ONLY : iverbosity, gamma_only
USE io_global, ONLY : stdout
USE vdW_DF, ONLY : vdw_type, initialize_spline_interpolation, &
numerical_gradient, interpolate_kernel
interpolate_kernel
USE gc_lr, ONLY : grho
@ -145,8 +145,8 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
!! Global variables
allocate(total_rho(dfftp%nnr) )
allocate(gradient_rho(dfftp%nnr, 3))
allocate(gradient_drho(dfftp%nnr, 3))
allocate(gradient_rho(3,dfftp%nnr))
allocate(gradient_drho(dfftp%nnr,3))
allocate(q0(dfftp%nnr), q(dfftp%nnr))
allocate(dq0_dq(dfftp%nnr), d2q0_dq2(dfftp%nnr))
allocate(dq_dn_n(dfftp%nnr), dn_dq_dn_n_n(dfftp%nnr), dq_dgradn_n_gmod(dfftp%nnr))
@ -179,7 +179,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
!! -------------------------------------------------------------------------
total_rho(:) = rho(:,1)
call numerical_gradient(total_rho,gradient_rho)
call fft_gradient_r2r(dfftp,total_rho,g,gradient_rho)
CALL qgradient (q_point, dfftp%nnr, drho(:,1), ngm, g, dfftp%nl, alat, gradient_drho)
!! -------------------------------------------------------------------------
@ -297,7 +297,7 @@ subroutine get_delta_v(rho, drho, nspin, q_point, delta_v)
allocate(delta_h_aux(dfftp%nnr))
do icar = 1,3
delta_h(:) = (h1t(:) * gradient_rho(:,icar)+ h2t(:) * gradient_drho(:,icar))
delta_h(:) = (h1t(:) * gradient_rho(icar,:)+ h2t(:) * gradient_drho(:,icar))
CALL fwfft ('Rho', delta_h, dfftp)
@ -423,10 +423,10 @@ end subroutine get_delta_v
!!
!! Fractions
!!
gmod = sqrt(gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2)
gradn_graddeltan = gradient_rho(i_grid,1)*gradient_drho(i_grid,1) + &
gradient_rho(i_grid,2)*gradient_drho(i_grid,2) + &
gradient_rho(i_grid,3)*gradient_drho(i_grid,3)
gmod = sqrt(gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2)
gradn_graddeltan = gradient_rho(1,i_grid)*gradient_drho(i_grid,1) + &
gradient_rho(2,i_grid)*gradient_drho(i_grid,2) + &
gradient_rho(3,i_grid)*gradient_drho(i_grid,3)
END SUBROUTINE get_thetas_exentended
@ -503,9 +503,9 @@ end subroutine get_delta_v
sqrt_r_s = sqrt(r_s)
gc = -Z_ab/(36.0D0*kF*total_rho(i_grid)**2) &
* (gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2)
* (gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2)
gmod = sqrt(gradient_rho(i_grid,1)**2+gradient_rho(i_grid,2)**2+gradient_rho(i_grid,3)**2)
gmod = sqrt(gradient_rho(1,i_grid)**2+gradient_rho(2,i_grid)**2+gradient_rho(3,i_grid)**2)
LDA_1 = 8.0D0*pi/3.0D0*(LDA_A*(1.0D0+LDA_a1*r_s))
LDA_2 = 2.0D0*LDA_A * (LDA_b1*sqrt_r_s + LDA_b2*r_s + LDA_b3*r_s*sqrt_r_s + LDA_b4*r_s*r_s)
@ -725,5 +725,9 @@ subroutine qgradient (xq, nrxx, a, ngm, g, nl, alat, ga)
end subroutine qgradient
! ####################################################################
! | |
! | thetas_to_uk |
END MODULE ph_vdW_DF

72
Modules/xc_rVV10.f90
View File

@ -31,7 +31,6 @@ MODULE rVV10
public :: xc_rVV10, &
interpolate_kernel, &
initialize_spline_interpolation, &
numerical_gradient, &
stress_rVV10, b_value
CONTAINS
@ -120,7 +119,7 @@ CONTAINS
!! ---------------------------------------------------------------------------------------
allocate( q0(dfftp%nnr) )
allocate( gradient_rho(dfftp%nnr, 3) )
allocate( gradient_rho(3,dfftp%nnr) )
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
allocate( total_rho(dfftp%nnr) )
@ -137,7 +136,7 @@ CONTAINS
!! -------------------------------------------------------------------------
!! Here we calculate the gradient in reciprocal space using FFT
!! -------------------------------------------------------------------------
call numerical_gradient(total_rho,gradient_rho)
call fft_gradient_r2r( dfftp, total_rho, g, gradient_rho)
!! -------------------------------------------------------------------------
!! Get Q and all the derivatives
@ -271,7 +270,7 @@ CONTAINS
!! Allocations
!! ---------------------------------------------------------------------------------------
allocate( gradient_rho(dfftp%nnr, 3) )
allocate( gradient_rho(3,dfftp%nnr) )
allocate( total_rho(dfftp%nnr) )
allocate( q0(dfftp%nnr) )
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
@ -291,7 +290,7 @@ CONTAINS
!! -------------------------------------------------------------------------
!! Here we calculate the gradient in reciprocal space using FFT
!! -------------------------------------------------------------------------
call numerical_gradient(total_rho,gradient_rho)
call fft_gradient_r2r( dfftp, total_rho, g, gradient_rho)
!! -------------------------------------------------------------------------------------------------------------
!! Get q0.
@ -460,7 +459,7 @@ CONTAINS
do m = 1, l
sigma (l, m) = sigma (l, m) - prefactor * &
(gradient_rho(i_grid,l) * gradient_rho(i_grid,m))
(gradient_rho(l,i_grid) * gradient_rho(m,i_grid))
enddo
enddo
endif
@ -586,9 +585,9 @@ CONTAINS
if (total_rho(i_grid) > epsr) then
gmod2 = gradient_rho(i_grid,1)**2 + &
gradient_rho(i_grid,2)**2 + &
gradient_rho(i_grid,3)**2
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
!! ------------------------------------------------------------------------------------
@ -1011,59 +1010,6 @@ subroutine interpolate_Dkernel_Dk(k, dkernel_of_dk)
end subroutine interpolate_Dkernel_Dk
!! ###############################################################################################################
!! | |
!! | NUMERICAL_GRADIENT |
!! |_______________________|
!! Calculates the gradient of the charge density numerically on the grid. We use
!! the PWSCF gradient style.
subroutine numerical_gradient(total_rho, gradient_rho)
use gvect, ONLY : ngm, g
USE cell_base, ONLY : tpiba
USE fft_base, ONLY : dfftp
USE fft_interfaces, ONLY : fwfft, invfft
!
! I/O variables
!
real(dp), intent(in) :: total_rho(:) !! Input array holding total charge density.
real(dp), intent(out) :: gradient_rho(:,:) !! Output array that will holds the gradient
! !! of the charge density.
! local variables
!
integer :: icar !! counter on cartesian components
complex(dp), allocatable :: c_rho(:) !! auxiliary complex array for rho
complex(dp), allocatable :: c_grho(:) !! auxiliary complex array for grad rho
! rho in G space
allocate ( c_rho(dfftp%nnr), c_grho(dfftp%nnr) )
c_rho(1:dfftp%nnr) = CMPLX(total_rho(1:dfftp%nnr),0.0_DP)
CALL fwfft ('Rho', c_rho, dfftp)
do icar=1,3
! compute gradient in G space
c_grho(:) =CMPLX(0.0_DP,0.0_DP)
c_grho(dfftp%nl(:)) = CMPLX (0.0_DP,1.0_DP) * tpiba * g(icar,:) * c_rho(dfftp%nl(:))
if (gamma_only) c_grho( dfftp%nlm(:) ) = CONJG( c_grho( dfftp%nl(:) ) )
! back in real space
CALL invfft ('Rho', c_grho, dfftp)
gradient_rho(:,icar) = REAL( c_grho(:) )
end do
deallocate ( c_rho, c_grho )
!gradient_rho = 0.0D0
return
end subroutine numerical_gradient
!! #################################################################################################
!! | |
!! | thetas_to_uk |
@ -1311,7 +1257,7 @@ end subroutine vdW_energy
end do
do icar = 1,3
h(:) = CMPLX(h_prefactor(:) * gradient_rho(:,icar),0.0_DP)
h(:) = CMPLX(h_prefactor(:) * gradient_rho(icar,:),0.0_DP)
CALL fwfft ('Rho', h, dfftp)
h(dfftp%nl(:)) = CMPLX(0.0_DP,1.0_DP) * tpiba * g(icar,:) * h(dfftp%nl(:))
if (gamma_only) h(dfftp%nlm(:)) = CONJG(h(dfftp%nl(:)))

124
Modules/xc_vdW_DF.f90
View File

@ -70,7 +70,7 @@ integer :: vdw_type = 1
private
public :: xc_vdW_DF, xc_vdW_DF_spin, stress_vdW_DF, interpolate_kernel, &
vdw_type, numerical_gradient, initialize_spline_interpolation
vdw_type, initialize_spline_interpolation
CONTAINS
@ -228,7 +228,7 @@ CONTAINS
real(dp), allocatable :: grad_rho(:,:) ! The gradient of the charge density. The
! format is as follows:
! grad_rho(grid_point, cartesian_component).
! grad_rho(cartesian_component,grid_point).
real(dp), allocatable :: potential(:) ! The vdW contribution to the potential
@ -278,7 +278,8 @@ CONTAINS
! Allocate arrays. nnr is a PWSCF variable that holds the number of
! points assigned to a given processor.
allocate( q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr), grad_rho(dfftp%nnr, 3) )
allocate( q0(dfftp%nnr), dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr), &
grad_rho(3,dfftp%nnr) )
allocate( total_rho(dfftp%nnr), potential(dfftp%nnr), thetas(dfftp%nnr, Nqs) )
@ -294,7 +295,7 @@ CONTAINS
! --------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT.
call numerical_gradient (total_rho, grad_rho)
call fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
! --------------------------------------------------------------------
@ -398,13 +399,11 @@ CONTAINS
real(dp), allocatable :: grad_rho(:,:) ! The gradient of the charge density. The
! format is as follows:
! grad_rho(grid_point, cartesian_component)
real(dp), allocatable :: grad_rho_up(:,:) ! The gradient of the up charge density. The
! format is as follows:
! grad_rho(grid_point, cartesian_component)
real(dp), allocatable :: grad_rho_down(:,:) ! The gradient of the down charge density. The
! format is as follows:
! grad_rho(grid_point, cartesian_component)
! 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.
@ -461,11 +460,11 @@ CONTAINS
! Allocate arrays. nnr is a PWSCF variable that holds the number of
! points assigned to a given processor.
allocate( q0(dfftp%nnr), total_rho(dfftp%nnr), grad_rho(dfftp%nnr, 3) )
allocate( q0(dfftp%nnr), total_rho(dfftp%nnr), grad_rho(3,dfftp%nnr) )
allocate( rho_up(dfftp%nnr), rho_down(dfftp%nnr) )
allocate( dq0_drho_up (dfftp%nnr), dq0_dgradrho_up (dfftp%nnr) )
allocate( dq0_drho_down(dfftp%nnr), dq0_dgradrho_down(dfftp%nnr) )
allocate( grad_rho_up(dfftp%nnr, 3), grad_rho_down(dfftp%nnr, 3) )
allocate( grad_rho_up(3,dfftp%nnr), grad_rho_down(3,dfftp%nnr) )
allocate( potential_up(dfftp%nnr), potential_down(dfftp%nnr) )
allocate( thetas(dfftp%nnr, Nqs) )
@ -489,10 +488,9 @@ CONTAINS
! --------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT.
call numerical_gradient ( total_rho, grad_rho )
call numerical_gradient ( rho_up, grad_rho_up )
call numerical_gradient ( rho_down, grad_rho_down )
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
@ -634,7 +632,7 @@ CONTAINS
r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
s = sqrt( grad_rho(i_grid,1)**2 + grad_rho(i_grid,2)**2 + grad_rho(i_grid,3)**2 ) &
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 )
@ -799,15 +797,15 @@ CONTAINS
if (calc_qx_up) then
s_up = sqrt( grad_rho_up(i_grid,1)**2 + grad_rho_up(i_grid,2)**2 + &
grad_rho_up(i_grid,3)**2 ) / (2.0D0 * kF(up) * up)
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(i_grid,1)**2 + grad_rho_down(i_grid,2)**2 + &
grad_rho_down(i_grid,3)**2) / (2.0D0 * kF(down) * down)
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
@ -1212,10 +1210,10 @@ CONTAINS
do icar = 1,3
h(:) = CMPLX( h_prefactor(:) * grad_rho(:,icar), 0.0_DP )
h(:) = CMPLX( h_prefactor(:) * grad_rho(icar,:), 0.0_DP )
do i_grid = 1, dfftp%nnr
gradient2 = grad_rho(i_grid,1)**2 + grad_rho(i_grid,2)**2 + grad_rho(i_grid,3)**2
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
@ -1587,7 +1585,7 @@ CONTAINS
! --------------------------------------------------------------------
! Allocations
allocate( grad_rho(dfftp%nnr, 3) )
allocate( grad_rho(3,dfftp%nnr) )
allocate( total_rho(dfftp%nnr) )
allocate( q0(dfftp%nnr) )
allocate( dq0_drho(dfftp%nnr), dq0_dgradrho(dfftp%nnr) )
@ -1608,8 +1606,7 @@ CONTAINS
! --------------------------------------------------------------------
! Here we calculate the gradient in reciprocal space using FFT.
call numerical_gradient (total_rho, grad_rho)
call fft_gradient_r2r (dfftp, total_rho, g, grad_rho)
! --------------------------------------------------------------------
! Get q0.
@ -1720,7 +1717,7 @@ CONTAINS
q_low = 1
q_hi = Nqs
grad2 = sqrt( grad_rho(i_grid,1)**2 + grad_rho(i_grid,2)**2 + grad_rho(i_grid,3)**2 )
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
! -----------------------------------------------------------------
@ -1762,7 +1759,7 @@ CONTAINS
do m = 1, l
sigma (l, m) = sigma (l, m) - e2 * prefactor * &
(grad_rho(i_grid,l) * grad_rho(i_grid,m))
(grad_rho(l,i_grid) * grad_rho(m,i_grid))
end do
end do
@ -1922,75 +1919,6 @@ CONTAINS
END SUBROUTINE interpolate_Dkernel_Dk
! ####################################################################
! | |
! | NUMERICAL_GRADIENT |
! |_______________________|
!
! Calculates the gradient of the charge density numerically on the
! grid. We use the PWSCF gradient style.
SUBROUTINE numerical_gradient (total_rho, grad_rho)
USE gvect, ONLY : ngm, g
USE cell_base, ONLY : tpiba
implicit none
real(dp), intent(in) :: total_rho(:) ! Input array holding total charge density.
real(dp), intent(out) :: grad_rho(:,:) ! Output array that will holds the gradient
! of the charge density.
integer :: icar ! Counter on cartesian components.
complex(dp), allocatable :: c_rho(:) ! auxiliary complex array for rho.
complex(dp), allocatable :: c_grho(:) ! auxiliary complex array for grad rho.
! --------------------------------------------------------------------
! Rho in G space.
allocate ( c_rho(dfftp%nnr), c_grho(dfftp%nnr) )
c_rho(1:dfftp%nnr) = CMPLX(total_rho(1:dfftp%nnr),0.0_DP)
CALL fwfft ('Rho', c_rho, dfftp)
do icar=1,3
! -----------------------------------------------------------------
! Compute gradient in G space.
c_grho(:) =CMPLX(0.0_DP,0.0_DP)
c_grho(dfftp%nl(:)) = CMPLX (0.0_DP,1.0_DP) * tpiba * g(icar,:) * c_rho(dfftp%nl(:))
if (gamma_only) c_grho( dfftp%nlm(:) ) = CONJG( c_grho( dfftp%nl(:) ) )
! -----------------------------------------------------------------
! Back in real space.
CALL invfft ('Rho', c_grho, dfftp)
grad_rho(:,icar) = REAL( c_grho(:) )
end do
deallocate ( c_rho, c_grho )
END SUBROUTINE numerical_gradient
! ####################################################################
! | |
! | thetas_to_uk |