More harmonization of gradients and second derivatives. BEWARE:

cannot guarantee that what has been modified work as before
This commit is contained in:
Paolo Giannozzi 2018-02-14 13:48:37 +01:00
parent 8faa1a355e
commit 983b85cbb3
6 changed files with 83 additions and 116 deletions

View File

@ -45,7 +45,6 @@ exx_vofr.o \
fft.o \
forces.o \
fromscra.o \
gradrho.o \
gram.o \
gtable.o \
init.o \

View File

@ -1,86 +0,0 @@
!
! Copyright (C) 2002-2020 Quantum ESPRESSO grouo
! 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 .
!
!----------------------------------------------------------------------
SUBROUTINE gradrho( nspin, rhog, d2rho, dxdyrho, dxdzrho, dydzrho )
!----------------------------------------------------------------------
!
! calculates gradient of charge density for gradient corrections
! in: charge density on G-space out: gradient in R-space
!
use cell_base
use gvect, only: g
USE fft_interfaces, ONLY: invfft
USE fft_base, ONLY: dfftp
USE fft_helper_subroutines, ONLY: fftx_oned2threed
!
implicit none
! input
integer nspin
complex(kind=8) rhog(dfftp%ngm,nspin)
! output
real(kind=8):: d2rho(3,dfftp%nnr), dxdyrho(dfftp%nnr), &
dxdzrho(dfftp%nnr), dydzrho(dfftp%nnr)
! local
complex(kind=8), allocatable:: v(:), drhog(:,:)
complex(kind=8) ci
integer iss, ig, ir, j
!
!
allocate(v(dfftp%nnr))
allocate(drhog(dfftp%ngm,3))
ci=(0.0d0,1.0d0)
d2rho = 0.d0
dxdyrho = 0.d0
dxdzrho = 0.d0
dydzrho = 0.d0
do iss=1,nspin
do ig=1,dfftp%ngm
drhog(ig,1) = -1.d0*tpiba**2*g(1,ig)**2*rhog(ig,iss)
drhog(ig,2) = -1.d0*tpiba**2*g(2,ig)**2*rhog(ig,iss)
drhog(ig,3) = -1.d0*tpiba**2*g(3,ig)**2*rhog(ig,iss)
enddo
CALL fftx_oned2threed(dfftp, v, drhog(:,1) )
call invfft('Rho',v, dfftp )
do ir=1,dfftp%nnr
d2rho(1,ir)=d2rho(1,ir)+real(v(ir))
end do
CALL fftx_oned2threed(dfftp, v, drhog(:,2), drhog(:,3) )
call invfft('Rho',v, dfftp )
do ir=1,dfftp%nnr
d2rho(2,ir)=d2rho(2,ir)+real(v(ir))
d2rho(3,ir)=d2rho(3,ir)+aimag(v(ir))
end do
do ig=1,dfftp%ngm
drhog(ig,1) = -1.d0*tpiba**2*g(1,ig)*g(2,ig)*rhog(ig,iss)
drhog(ig,2) = -1.d0*tpiba**2*g(1,ig)*g(3,ig)*rhog(ig,iss)
drhog(ig,3) = -1.d0*tpiba**2*g(2,ig)*g(3,ig)*rhog(ig,iss)
enddo
CALL fftx_oned2threed(dfftp, v, drhog(:,1) )
CALL invfft('Rho',v, dfftp )
do ir=1,dfftp%nnr
dxdyrho(ir)=dxdyrho(ir)+real(v(ir))
end do
CALL fftx_oned2threed(dfftp, v, drhog(:,2), drhog(:,3) )
call invfft('Rho',v, dfftp )
do ir=1,dfftp%nnr
dxdzrho(ir)=dxdzrho(ir)+real(v(ir))
dydzrho(ir)=dydzrho(ir)+aimag(v(ir))
end do
end do
deallocate(v)
deallocate(drhog)
END SUBROUTINE gradrho

View File

@ -108,7 +108,6 @@ cp_interfaces.o : ../../LAXlib/la_types.o
cp_interfaces.o : ../../Modules/fft_base.o
cp_interfaces.o : ../../Modules/ions_base.o
cp_interfaces.o : ../../Modules/kind.o
cp_interfaces.o : ../../Modules/recvec.o
cp_restart.o : ../../LAXlib/mp_diag.o
cp_restart.o : ../../Modules/cell_base.o
cp_restart.o : ../../Modules/constants.o
@ -517,11 +516,6 @@ fromscra.o : mainvar.o
fromscra.o : modules.o
fromscra.o : ortho_base.o
fromscra.o : printout_base.o
gradrho.o : ../../FFTXlib/fft_helper_subroutines.o
gradrho.o : ../../FFTXlib/fft_interfaces.o
gradrho.o : ../../Modules/cell_base.o
gradrho.o : ../../Modules/fft_base.o
gradrho.o : ../../Modules/recvec.o
gram.o : ../../Modules/electrons_base.o
gram.o : ../../Modules/gvecw.o
gram.o : ../../Modules/ions_base.o

View File

@ -54,8 +54,6 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
real(kind=8) gxl, xyr, xzr, yzr
real(kind=8), allocatable:: vec(:,:,:), aiuto(:,:,:)
real(kind=8), allocatable:: drho(:,:), d2rho(:,:)
real(kind=8), allocatable:: dxdyrho(:), dxdzrho(:)
real(kind=8), allocatable:: dydzrho(:)
real(kind=8), allocatable:: tauv(:,:,:)
complex(kind=8) ci
@ -71,10 +69,7 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
integer displs(nproc), ip, me
#endif
if (abisur) allocate(drho(3,dfftp%nnr))
if (abisur) allocate(d2rho(3,dfftp%nnr))
if (abisur) allocate(dxdyrho(dfftp%nnr))
if (abisur) allocate(dxdzrho(dfftp%nnr))
if (abisur) allocate(dydzrho(dfftp%nnr))
if (abisur) allocate(d2rho(6,dfftp%nnr))
call start_clock( 'vol_clu' )
@ -243,9 +238,9 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
IF (abisur) THEN
DO iss = 1, nspin
CALL fft_gradient_g2r ( dfftp, rhotmp(1,iss), g, drho(1,iss) )
CALL fft_gradient_g2r( dfftp, rhotmp(1,iss), g, drho(1,iss) )
CALL fft_hessian_g2r ( dfftp, rhotmp, g, d2rho(1,iss) )
END DO
CALL gradrho( nspin, rhotmp, d2rho, dxdyrho, dxdzrho, dydzrho )
END IF
CALL rho_g2r( dfftp, rhotmp, rho_gaus )
@ -349,17 +344,15 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
end if
if (abisur) then
modr = 0.d0
lap = 0.d0
gxl = 0.d0
do j = 1,3
modr = modr + drho(j,ir)**2
lap = lap + d2rho(j,ir)
gxl = gxl + drho(j,ir)**2*d2rho(j,ir)
end do
xyr = 2.d0*dxdyrho(ir)*drho(1,ir)*drho(2,ir)
xzr = 2.d0*dxdzrho(ir)*drho(1,ir)*drho(3,ir)
yzr = 2.d0*dydzrho(ir)*drho(2,ir)*drho(3,ir)
! for d2rho: 1=xx, 2=xy, 3=yy, 4=xz, 5=yz, 6=zz
modr = drho(1,ir)**2 + drho(2,ir)**2 + drho(3,ir)**2
lap = d2rho(1,ir) + d2rho(3,ir) + d2rho(6,ir)
gxl = drho(1,ir)**2*d2rho(1,ir) + &
drho(2,ir)**2*d2rho(3,ir) + &
drho(3,ir)**2*d2rho(6,ir)
xyr = 2.d0*d2rho(2,ir)*drho(1,ir)*drho(2,ir)
xzr = 2.d0*d2rho(4,ir)*drho(1,ir)*drho(3,ir)
yzr = 2.d0*d2rho(5,ir)*drho(2,ir)*drho(3,ir)
modr = dsqrt(modr)
surfclu = surfclu + (wpiu-wmeno)*modr
v_vol(ir) = v_vol(ir) -1.d0*Surf_t/dthr * (wpiu-wmeno) * &
@ -387,9 +380,6 @@ SUBROUTINE vol_clu(rho_real,rho_g,flag)
deallocate( tauv )
if ( abisur ) deallocate( drho )
if ( abisur ) deallocate( d2rho )
if ( abisur ) deallocate( dxdyrho )
if ( abisur ) deallocate( dxdzrho )
if ( abisur ) deallocate( dydzrho )
call stop_clock( 'vol_clu' )

View File

@ -516,6 +516,76 @@ END SUBROUTINE fft_laplacian
! Routines computing hessian via FFT
!--------------------------------------------------------------------
!
!----------------------------------------------------------------------
SUBROUTINE fft_hessian_g2r ( dfft, a, g, ha )
!----------------------------------------------------------------------
!
! ... Calculates ha = hessian(a)
! ... input : dfft FFT descriptor
! ... a(:) a real function on the real-space FFT grid
! ... g(3,:) G-vectors, in (2\pi/a)^2 units
! ... output: ha(6,:) hessian(a), real, on the real-space FFT grid
! ... lower-packed matrix indeces 1-6 correspond to:
! ... 1 = xx, 2 = yx, 3 = yy, 4 = zx, 5 = zy, 6 = zz
!
USE kinds, ONLY : DP
USE cell_base, ONLY : tpiba
USE fft_types, ONLY : fft_type_descriptor
USE fft_interfaces,ONLY : fwfft, invfft
USE fft_helper_subroutines, ONLY: fftx_oned2threed
!
IMPLICIT NONE
!
TYPE(fft_type_descriptor),INTENT(IN) :: dfft
REAL(DP), INTENT(IN) :: g(3,dfft%ngm)
COMPLEX(DP), INTENT(IN) :: a(dfft%ngm)
REAL(DP), INTENT(OUT) :: ha( 6, dfft%nnr )
!
INTEGER :: ig, ir
COMPLEX(DP), ALLOCATABLE :: aux(:), haux(:,:)
!
IF ( .NOT. dfft%lgamma ) CALL errore ('fft_hessian_g2r',&
'only gamma case is implemented',1)
ALLOCATE ( aux(dfft%nnr))
ALLOCATE (haux(dfft%ngm,2))
! xx, yx
DO ig=1,dfft%ngm
haux(ig,1) = -tpiba**2*g(1,ig)**2 *a(ig)
haux(ig,2) = -tpiba**2*g(1,ig)*g(2,ig)*a(ig)
END DO
CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) )
CALL invfft('Rho', aux, dfft)
DO ir=1,dfft%nnr
ha(1,ir) = DBLE(aux(ir))
ha(2,ir) =AIMAG(aux(ir))
END DO
! yy, zx
DO ig=1,dfft%ngm
haux(ig,1) = -tpiba**2*g(2,ig)**2 *a(ig)
haux(ig,2) = -tpiba**2*g(1,ig)*g(3,ig)*a(ig)
END DO
CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) )
CALL invfft('Rho', aux, dfft)
DO ir=1,dfft%nnr
ha(3,ir) = DBLE(aux(ir))
ha(4,ir) =AIMAG(aux(ir))
END DO
! zy, zz
DO ig=1,dfft%ngm
haux(ig,1) = -tpiba**2*g(2,ig)*g(3,ig)*a(ig)
haux(ig,2) = -tpiba**2*g(3,ig)**2 *a(ig)
END DO
CALL fftx_oned2threed( dfft, aux, haux(:,1), haux(:,2) )
CALL invfft('Rho', aux, dfft)
DO ir=1,dfft%nnr
ha(5,ir) = DBLE(aux(ir))
ha(6,ir) =AIMAG(aux(ir))
END DO
!
DEALLOCATE(aux)
DEALLOCATE(haux)
END SUBROUTINE fft_hessian_g2r
!--------------------------------------------------------------------
SUBROUTINE fft_hessian( dfft, a, g, ga, ha )
!--------------------------------------------------------------------

View File

@ -91,7 +91,6 @@ fft_rho.o : ../FFTXlib/fft_interfaces.o
fft_rho.o : ../FFTXlib/fft_types.o
fft_rho.o : control_flags.o
fft_rho.o : kind.o
fox_as_iotk.o : kind.o
funct.o : io_global.o
funct.o : kind.o
funct.o : xc_rVV10.o
@ -107,6 +106,7 @@ generate_function.o : io_global.o
generate_function.o : kind.o
generate_function.o : mp_bands.o
generate_k_along_lines.o : kind.o
gradutils.o : ../FFTXlib/fft_helper_subroutines.o
gradutils.o : ../FFTXlib/fft_interfaces.o
gradutils.o : ../FFTXlib/fft_types.o
gradutils.o : cell_base.o