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

New PAW version, now completely working on radial grid. Code on FFT grid can still be compiled using __GRID_PAW flag (but may not work a 

nymore). Other things working: parallel, gamma-only.

Files vxc_t, exc_t and vxcgc moved from atomic to Modules (to prevent cyclic dependencies).

Other random fixes: a kind in init_vloc, a call to infomsg in ./PH/add_for_charges.f90.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4358 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
paulatto 2007-10-22 14:54:29 +00:00
parent 163557457a
commit d858b4a51b
40 changed files with 2140 additions and 1039 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CPV/make.depend
View File

@ -904,7 +904,6 @@ problem_size.o : ../Modules/fft_base.o
problem_size.o : ../Modules/io_global.o
problem_size.o : ../Modules/ions_base.o
problem_size.o : ../Modules/kind.o
problem_size.o : ../Modules/radial_grids.o
problem_size.o : ../Modules/recvec.o
problem_size.o : ../Modules/uspp.o
problem_size.o : electrons.o

3
D3/Makefile
View File

@ -69,6 +69,7 @@ MODULES = \
../Modules/control_flags.o \
../Modules/descriptors.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -103,6 +104,8 @@ MODULES = \
../Modules/uspp.o \
../Modules/timestep.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/xml_io_base.o

3
Gamma/Makefile
View File

@ -18,6 +18,7 @@ MODULES = \
../Modules/descriptors.o \
../Modules/dspev_drv.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -61,6 +62,8 @@ MODULES = \
../Modules/upf_to_internal.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/xml_io_base.o \
../Modules/zhpev_drv.o

3
Modules/Makefile
View File

@ -17,6 +17,7 @@ descriptors.o \
dspev_drv.o \
electrons_base.o \
energies.o \
exc_t.o \
fft_base.o \
fft_scalar.o \
fft_types.o \
@ -68,6 +69,8 @@ timestep.o \
version.o \
upf_to_internal.o \
uspp.o \
vxc_t.o \
vxcgc.o \
wave_base.o \
wavefunctions.o \
wannier.o \

9
Modules/clocks.f90
View File

@ -5,7 +5,6 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!#define __TRACE
!----------------------------------------------------------------------------
MODULE mytime
!----------------------------------------------------------------------------
@ -95,8 +94,8 @@ SUBROUTINE start_clock( label )
! ... store in t0cpu the starting time
!
IF ( t0cpu(n) /= notrunning ) THEN
WRITE( stdout, '("start_clock: clock # ",I2," for ",A12, &
& " already started")' ) n, label
! WRITE( stdout, '("start_clock: clock # ",I2," for ",A12, &
! & " already started")' ) n, label
ELSE
t0cpu(n) = scnds()
IF ( n == 1 ) t0wall = cclock()
@ -159,8 +158,8 @@ SUBROUTINE stop_clock( label )
!
IF ( t0cpu(n) == notrunning ) THEN
!
WRITE( stdout, '("stop_clock: clock # ",I2," for ",A12, &
& " not running")' ) n, label
! WRITE( stdout, '("stop_clock: clock # ",I2," for ",A12, &
! & " not running")' ) n, label
!
ELSE
!

1
Modules/constants.f90
View File

@ -88,6 +88,7 @@ MODULE constants
REAL(DP), PARAMETER :: eps12 = 1.0E-12_DP
REAL(DP), PARAMETER :: eps14 = 1.0E-14_DP
REAL(DP), PARAMETER :: eps16 = 1.0E-16_DP
REAL(DP), PARAMETER :: eps24 = 1.0E-24_DP
REAL(DP), PARAMETER :: eps32 = 1.0E-32_DP
!
REAL(DP), PARAMETER :: gsmall = 1.0E-12_DP

3
atomic/exc_t.f90 → Modules/exc_t.f90
View File

@ -39,6 +39,9 @@ function exc_t(rho,rhoc,lsd)
arho = abs(rhot)
if (arho.gt.1.e-30_DP) then
zeta = (rho(1)-rho(2)) / arho
! In atomic this cannot happen, but in PAW zeta can become
! a little larger than 1, or smaller than -1:
zeta = MAX(MIN(zeta, 1._dp),-1._dp)
call xc_spin(arho,zeta,ex,ec,vx(1),vx(2),vc(1),vc(2))
exc_t=e2*(ex+ec)
endif

60
Modules/grid_paw_variables.f90
View File

@ -4,7 +4,6 @@ module grid_paw_variables
!
! NO spin-orbit
! NO EXX
! NO Parallelism
! NO rinner > 0
!
USE kinds, ONLY : DP
@ -25,22 +24,30 @@ module grid_paw_variables
! Analogous to tvanp in "uspp_param" (Modules/uspp.f90)
LOGICAL :: &
tpawp(npsx) ! if .TRUE. the atom is of PAW type
tpawp(npsx) = .false. ! if .TRUE. the atom is of PAW type
! Analogous to qfunc in "uspp_param" (Modules/uspp.f90)
REAL(DP), TARGET :: &
pfunc(ndmx,nbrx,nbrx,npsx), &! AE: \phi_{mu}(|r|)-\phi_{nu}(|r|)
ptfunc(ndmx,nbrx,nbrx,npsx) ! PS: \tilde{\phi}_{mu}(|r|)-\tilde{\phi}_{nu}(|r|)
REAL(DP), TARGET :: &
!augfun(ndmx,nbrx,nbrx,0:lqmax,npsx), & ! changed to aug of type augfun_t (see below)
pmultipole(nbrx,nbrx,0:lqmax,npsx), &! AE multipoles
ptmultipole(nbrx,nbrx,0:lqmax,npsx) ! PS multipoles
! Augmentation on radial grid:
TYPE augfun_t
REAL(DP), ALLOCATABLE :: fun(:,:,:,:)
END TYPE
TYPE(augfun_t) :: aug(npsx)
! Moments of the augmentation functions
REAL (DP) :: &
augmom(nbrx,nbrx,0:6,npsx) ! moments of PAW augm. functions
INTEGER :: &
nraug(npsx) ! augm. functions cutoff parameter
#ifdef __GRID_PAW
REAL(DP), TARGET :: &
!augfun(ndmx,nbrx,nbrx,0:lqmax,npsx), & ! changed to aug of type augfun_t (see below)
pmultipole(nbrx,nbrx,0:lqmax,npsx), &! AE multipoles
ptmultipole(nbrx,nbrx,0:lqmax,npsx) ! PS multipoles
! Analogous to qq in "uspp_param" (Modules/uspp.f90)
REAL(DP), ALLOCATABLE, TARGET :: &
@ -59,22 +66,10 @@ module grid_paw_variables
prod0p(:,:,:), &! k=0 AE product in reciprocal space
prod0pt(:,:,:) ! k=0 PS product in reciprocal space
!! NEW-AUG !!
! Moments of the augmentation functions
! REAL (DP) :: &
! r2(ndmx,npsx) ! r**2 logarithmic mesh
REAL (DP) :: &
augmom(nbrx,nbrx,0:6,npsx) ! moments of PAW augm. functions
INTEGER :: &
nraug(npsx) ! augm. functions cutoff parameter
!! NEW-AUG !!
! Analogous to rho in "scf" (PW/pwcom.f90) + index scanning atoms
REAL(DP), ALLOCATABLE, TARGET :: &
rho1(:,:,:), &! 1center AE charge density in real space
rho1t(:,:,:) ! 1center PS charge density in real space
!!! No more needed since ptfunc already contains the augmentation charge qfunc
!!! rho1h(:,:,:) ! 1center compensation charge in real space
! Analogous to vr in "scf" (PW/pwcom.f90) + index scanning atoms
REAL(DP), ALLOCATABLE, TARGET :: &
@ -85,26 +80,32 @@ module grid_paw_variables
REAL(DP), ALLOCATABLE, TARGET :: &
int_r2pfunc(:,:,:), &! Integrals of r^2 * pfunc(r) (AE)
int_r2ptfunc(:,:,:) ! Integrals of r^2 * pfunc(r) (PS)
#endif
! Analogous to rho_atc in "atom" (Modules/atom.f90)
REAL(DP), TARGET :: &
aerho_atc(ndmx,npsx), &! radial AE core charge density
psrho_atc(ndmx,npsx) ! radial PS core charge density
! Analogous to rho_core in "scf" (PW/pwcom.f90)
REAL(DP), ALLOCATABLE, TARGET :: &
aerho_core(:,:), &! AE core charge density in real space
psrho_core(:,:) ! PS core charge density in real space
! Analogous to vloc_at in "uspp_param" (Modules/uspp.f90)
! actually pseudopotential (AE and PS) on radial grid.
REAL(DP), TARGET :: &
aevloc_at(ndmx,npsx), &! AE descreened potential
psvloc_at(ndmx,npsx) ! PS descreened potential
! Analogous to vloc in "vlocal" (PW/pwcom.f90)
REAL(DP), ALLOCATABLE, TARGET :: &
#ifdef __GRID_PAW
aevloc(:,:), &! AE local 1-c potential for each atom type
#endif
psvloc(:,:) ! PS local 1-c potential for each atom type
#ifdef __GRID_PAW
! Analogous to rho_core in "scf" (PW/pwcom.f90)
REAL(DP), ALLOCATABLE, TARGET :: &
aerho_core(:,:), &! AE core charge density in real space
psrho_core(:,:) ! PS core charge density in real space
!
REAL(DP), ALLOCATABLE :: &
radial_distance(:), & ! radial distance from origin (minimum image conv)
@ -123,6 +124,7 @@ module grid_paw_variables
ehart1t(:), & ! Hartree energy (PS)
etxc1t (:), & ! XC: energy (PS)
vtxc1t (:) ! XC: Int V*rho (PS)
#endif
! Analogous to dion in "uspp_param" (Modules/uspp.f90)
REAL(DP) :: &
@ -133,19 +135,21 @@ module grid_paw_variables
dpaw_ae(:,:,:,:), &! AE D: D^1_{ij} (except for K.E.)
dpaw_ps(:,:,:,:) ! PS D: \tilde{D}^1_{ij} (except for K.E.)
#ifdef __GRID_PAW
! TMP analogous to rhonew in PW/electrons.f90
REAL(DP), ALLOCATABLE, TARGET :: &
rho1new(:,:,:), &! new 1center AE charge density in real space
rho1tnew(:,:,:) ! new 1center PS charge density in real space
! analogous to deband and descf in PW/electrons.f90
REAL(DP) :: deband_1ae, deband_1ps, descf_1ae, descf_1ps
#endif
! new vectors needed for mixing of augm. channel occupations
REAL(DP), ALLOCATABLE :: &
becnew(:,:,:) ! new augmentation channel occupations
! analogous to deband and descf in PW/electrons.f90
REAL(DP) :: deband_1ae, deband_1ps, descf_1ae, descf_1ps
CONTAINS
CONTAINS
! From PW/init_paw_1.f90
SUBROUTINE step_f(f2,f,r,nrs,nrc,pow,mesh)
USE kinds , ONLY : dp

10
Modules/make.depend
View File

@ -51,6 +51,8 @@ electrons_base.o : kind.o
electrons_base.o : timestep.o
energies.o : io_global.o
energies.o : kind.o
exc_t.o : functionals.o
exc_t.o : kind.o
fft_base.o : fft_types.o
fft_base.o : kind.o
fft_base.o : mp_global.o
@ -254,8 +256,14 @@ uspp.o : constants.o
uspp.o : kind.o
uspp.o : parameters.o
uspp.o : pseudo_types.o
uspp.o : radial_grids.o
uspp.o : random_numbers.o
vxc_t.o : functionals.o
vxc_t.o : io_global.o
vxc_t.o : kind.o
vxcgc.o : constants.o
vxcgc.o : functionals.o
vxcgc.o : kind.o
vxcgc.o : radial_grids.o
wannier.o : kind.o
wave_base.o : kind.o
wave_base.o : mp.o

6
Modules/pseudo_types.f90
View File

@ -67,7 +67,7 @@ END TYPE paw_t
CHARACTER(LEN=80):: comment !
CHARACTER(LEN=2) :: psd ! Element label
CHARACTER(LEN=20) :: typ ! Pseudo type ( NC or US )
LOGICAL :: tvanp ! .true. if Ultrasoft
LOGICAL :: tvanp ! .true. if Ultrasoft
LOGICAL :: nlcc ! Non linear core corrections
CHARACTER(LEN=20) :: dft ! Exch-Corr type
REAL(DP) :: zp ! z valence
@ -105,8 +105,8 @@ END TYPE paw_t
REAL(DP), POINTER :: rab(:) ! rab(mesh) dr(x)/dx (x=linear grid)
REAL(DP), POINTER :: rho_atc(:) ! rho_atc(mesh) atomic core charge
REAL(DP), POINTER :: vloc(:) ! vloc(mesh) local atomic potential
INTEGER, POINTER :: lll(:) ! lll(nbeta) l of each projector
INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) see above kkbeta
INTEGER, POINTER :: lll(:) ! lll(nbeta) l of each projector
INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) see above kkbeta
REAL(DP), POINTER :: beta(:,:) ! beta(mesh,nbeta) projectors
INTEGER :: nd
REAL(DP), POINTER :: dion(:,:) ! dion(nbeta,nbeta) atomic D_{mu,nu}

16
Modules/radial_grids.f90
View File

@ -48,7 +48,10 @@ TYPE radial_grid_type
r(ndmx), & ! the radial mesh
r2(ndmx), & ! the square of the radial mesh
rab(ndmx), & ! d r(x) / d x where x is the linear grid
sqr(ndmx) ! the square root of the radial mesh
sqr(ndmx), & ! the square root of the radial mesh
rm1(ndmx), & ! 1 / r
rm2(ndmx), & ! 1 / r
rm3(ndmx) ! 1 / r
REAL(DP) :: &
xmin, & ! the minimum x
rmax, & ! the maximum radial point
@ -63,7 +66,7 @@ END TYPE radial_grid_type
!============================================================================
!
CONTAINS
CONTAINS
!
! Build the radial (logarithmic) grid
!
@ -106,6 +109,9 @@ CONTAINS
grid%r2(i) = grid%r(i)*grid%r(i)
grid%rab(i) = grid%r(i)*dx
grid%sqr(i) = sqrt(grid%r(i))
grid%rm1(i) = 1._dp/grid%r(i)
grid%rm2(i) = 1._dp/grid%r(i)**2
grid%rm3(i) = 1._dp/grid%r(i)**3
end do
!
grid%mesh = mesh
@ -359,7 +365,7 @@ subroutine read_grid_from_file(iunit,grid)
integer, intent(in) :: iunit
integer :: n
!
READ(iunit,'(i8)') grid%mesh
READ(iunit,'(i8)') grid%mesh
READ(iunit,'(e20.10)') grid%dx
READ(iunit,'(e20.10)') grid%xmin
READ(iunit,'(e20.10)') grid%zmesh
@ -368,6 +374,10 @@ subroutine read_grid_from_file(iunit,grid)
READ(iunit,'(e20.10)') (grid%sqr(n), n=1,grid%mesh)
! READ(iunit,'(e20.10)') (grid%rab(n), n=1,grid%mesh)
grid%rab(1:grid%mesh) = grid%r(1:grid%mesh) * grid%dx
grid%rm1(1:grid%mesh) = 1._dp/grid%r(1:grid%mesh)
grid%rm2(1:grid%mesh) = 1._dp/grid%r2(1:grid%mesh)
grid%rm3(1:grid%mesh) = 1._dp/grid%r(1:grid%mesh)**3
return
end subroutine read_grid_from_file

80
Modules/read_paw.f90
View File

@ -49,7 +49,7 @@
write(un,'(a)') "scalars:"
WRITE(un,'(e20.10)') pawset_%zval
WRITE(un,'(e20.10)') pawset_%z
WRITE(un,'(L1)') pawset_%nlcc
WRITE(un,'(L)') pawset_%nlcc
WRITE(un,'(i8)') pawset_%nwfc
WRITE(un,'(i8)') pawset_%irc
WRITE(un,'(i8)') pawset_%lmax
@ -135,7 +135,7 @@
read(un, '(a)') dummy
READ (un,'(e20.10)') pawset_%zval
READ (un,'(e20.10)') pawset_%z
READ (un,'(L1)') pawset_%nlcc
READ (un,'(L)') pawset_%nlcc
READ (un,'(i8)') pawset_%nwfc
READ (un,'(i8)') pawset_%irc
READ (un,'(i8)') pawset_%lmax
@ -323,6 +323,7 @@ END SUBROUTINE deallocate_pseudo_paw
!---------------------------------------------------------------------
!#define __DO_NOT_CUTOFF_PAW_FUNC
#define __SHARP_CUTOFF_PAW_FUNC
subroutine set_pseudo_paw (is, pawset)
!---------------------------------------------------------------------
!
@ -333,18 +334,15 @@ subroutine set_pseudo_paw (is, pawset)
!
USE atom, ONLY: rgrid, msh, chi, oc, nchi, lchi, jchi, rho_at, &
rho_atc, nlcc
! USE pseud, ONLY: lloc, lmax
USE uspp_param, ONLY: upf
USE funct, ONLY: set_dft_from_name, dft_is_meta, dft_is_hybrid
!
! USE spin_orb, ONLY: lspinorb
USE pseudo_types
USE constants, ONLY: FPI
!
USE grid_paw_variables, ONLY : tpawp, pfunc, ptfunc, aevloc_at, psvloc_at, &
aerho_atc, psrho_atc, kdiff, &
augmom, nraug, step_f,aug !!NEW-AUG
!USE grid_paw_routines, ONLY : step_f
augmom, nraug, step_f, aug
!
implicit none
!
@ -361,6 +359,10 @@ subroutine set_pseudo_paw (is, pawset)
!
#if defined __DO_NOT_CUTOFF_PAW_FUNC
PRINT '(A)', 'WARNING __DO_NOT_CUTOFF_PAW_FUNC'
#else
#ifdef __SHARP_CUTOFF_PAW_FUNC
PRINT '(A)', 'WARNING __SHARP_CUTOFF_PAW_FUNC'
#endif
#endif
!
! Cutoffing: WARNING: arbitrary right now, for grid calculation
@ -378,7 +380,6 @@ subroutine set_pseudo_paw (is, pawset)
tpawp(is)=.true.
nlcc(is) = pawset%nlcc
call set_dft_from_name( pawset%dft )
! call which_dft( pawset%dft )
!
IF ( dft_is_meta() ) &
CALL errore( 'upf_to_internal', 'META-GGA not implemented in PWscf', 1 )
@ -397,17 +398,10 @@ subroutine set_pseudo_paw (is, pawset)
!
nchi(is)=0
do i=1, pawset%nwfc
#if defined __DEBUG_UPF_TO_INTERNAL
! take only occupied wfcs (to have exactly the same as for US)
if (pawset%oc(i)>0._dp) then
#endif
nchi(is)=nchi(is)+1
lchi(nchi(is),is)=pawset%l(i)
oc(nchi(is),is)=MAX(pawset%oc(i),0._DP)
chi(1:pawset%grid%mesh, nchi(is), is) = pawset%pswfc(1:pawset%grid%mesh, i)
#if defined __DEBUG_UPF_TO_INTERNAL
end if
#endif
end do
!
upf(is)%nbeta= pawset%nwfc
@ -415,25 +409,27 @@ subroutine set_pseudo_paw (is, pawset)
do nb=1,pawset%nwfc
upf(is)%kbeta(nb)=pawset%ikk(nb)
end do
! kkbeta is the maximum distance from nucleus where AE
! and PS wavefunction may not match:
upf(is)%kkbeta=MAXVAL (upf(is)%kbeta(:))
allocate (upf(is)%beta(1:pawset%grid%mesh, 1:pawset%nwfc))
upf(is)%beta(1:pawset%grid%mesh, 1:pawset%nwfc) = &
pawset%proj(1:pawset%grid%mesh, 1:pawset%nwfc)
pawset%proj(1:pawset%grid%mesh, 1:pawset%nwfc)
allocate(upf(is)%dion(1:pawset%nwfc, 1:pawset%nwfc))
upf(is)%dion(1:pawset%nwfc, 1:pawset%nwfc) = pawset%dion(1:pawset%nwfc, 1:pawset%nwfc)
kdiff(1:pawset%nwfc, 1:pawset%nwfc, is) = pawset%kdiff(1:pawset%nwfc, 1:pawset%nwfc)
! HOPE!
! lmax(is) = pawset%lmax
upf(is)%nqlc = 2*pawset%lmax+1
upf(is)%nqf = 0 !! no rinner, all numeric
allocate (upf(is)%lll(pawset%nwfc) )
upf(is)%lll(1:pawset%nwfc) = pawset%l(1:pawset%nwfc)
allocate (upf(is)%rinner(upf(is)%nqlc))
upf(is)%rinner(1:upf(is)%nqlc) = 0._dp !! no rinner, all numeric
!
! integral of augmentation charges vanishes for different values of l
!
allocate ( upf(is)%qqq(pawset%nwfc,pawset%nwfc))
do i = 1, pawset%nwfc
do j = 1, pawset%nwfc
@ -444,7 +440,8 @@ subroutine set_pseudo_paw (is, pawset)
end if
end do
end do
!! NEW-AUG !!
! import augmentation charge:
nraug(is) = pawset%irc
rgrid(is)%dx = pawset%grid%dx
rgrid(is)%r2 (1:pawset%grid%mesh) = pawset%grid%r2 (1:pawset%grid%mesh)
@ -465,13 +462,14 @@ subroutine set_pseudo_paw (is, pawset)
pawset%augfun(1:pawset%grid%mesh,nb,mb,0)
enddo
enddo
! augfun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax,is) = &
! pawset%augfun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax)
! new sparse allocation for augmentation functions
if (allocated(aug(is)%fun)) deallocate(aug(is)%fun)
allocate(aug(is)%fun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax))
aug(is)%fun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax) = &
pawset%augfun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax)
!
aug(is)%fun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax) &
= &
pawset%augfun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax)
!
@ -485,6 +483,14 @@ subroutine set_pseudo_paw (is, pawset)
pawset%aewfc(1:pawset%grid%mesh, i) * pawset%aewfc(1:pawset%grid%mesh, j)
ptfunc (1:pawset%grid%mesh, i, j, is) = &
pawset%pswfc(1:pawset%grid%mesh, i) * pawset%pswfc(1:pawset%grid%mesh, j)
#else
#if defined __SHARP_CUTOFF_PAW_FUNC
pfunc(:,i,j,is) = 0._dp
ptfunc(:,i,j,is) = 0._dp
pfunc (1:pawset%irc, i, j, is) = &
pawset%aewfc(1:pawset%irc, i) * pawset%aewfc(1:pawset%irc, j)
ptfunc (1:pawset%irc, i, j, is) = &
pawset%pswfc(1:pawset%irc, i) * pawset%pswfc(1:pawset%irc, j)
#else
aux(1:pawset%grid%mesh) = pawset%aewfc(1:pawset%grid%mesh, i) * &
pawset%aewfc(1:pawset%grid%mesh, j)
@ -494,6 +500,7 @@ subroutine set_pseudo_paw (is, pawset)
pawset%pswfc(1:pawset%grid%mesh, j)
CALL step_f( ptfunc(1:pawset%grid%mesh,i,j,is), aux(1:pawset%grid%mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh)
#endif
#endif
end do
end do
@ -519,10 +526,14 @@ subroutine set_pseudo_paw (is, pawset)
rgrid(is)%r2(1:pawset%grid%mesh) = pawset%grid%r2(1:pawset%grid%mesh)
rgrid(is)%rab(1:pawset%grid%mesh) = pawset%grid%rab(1:pawset%grid%mesh)
rgrid(is)%sqr(1:pawset%grid%mesh) = sqrt(pawset%grid%r(1:pawset%grid%mesh))
! these speed up a lot a few calculations:
rgrid(is)%rm1(1:pawset%grid%mesh) = pawset%grid%rm1(1:pawset%grid%mesh)
rgrid(is)%rm2(1:pawset%grid%mesh) = pawset%grid%rm2(1:pawset%grid%mesh)
rgrid(is)%rm3(1:pawset%grid%mesh) = pawset%grid%rm3(1:pawset%grid%mesh)
! NO spin orbit PAW implemented right now (oct 2005)
!!$ if (lspinorb.and..not.pawset%has_so) &
!!$ call infomsg ('pawset_to_internal','At least one non s.o. pseudo', -1)
!!$ call infomsg ('pawset_to_internal','At least one non s.o. pseudo')
!!$
!!$ lspinorb=lspinorb.and.pawset%has_so
!!$ if (pawset%has_so) then
@ -558,13 +569,18 @@ subroutine set_pseudo_paw (is, pawset)
aevloc_at(1:pawset%grid%mesh,is) = pawset%aeloc(1:pawset%grid%mesh)
psvloc_at(1:pawset%grid%mesh,is) = pawset%psloc(1:pawset%grid%mesh)
#else
aux(1:pawset%grid%mesh) = pawset%aeloc(1:pawset%grid%mesh)
CALL step_f( aevloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh)
aux(1:pawset%grid%mesh) = pawset%psloc(1:pawset%grid%mesh)
CALL step_f( psvloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh)
DEALLOCATE (aux)
#if defined __SHARP_CUTOFF_PAW_FUNC
aevloc_at(1:pawset%grid%mesh,is) = pawset%aeloc(1:pawset%grid%mesh)
psvloc_at(1:pawset%grid%mesh,is) = pawset%psloc(1:pawset%grid%mesh)
aux(1:pawset%grid%mesh) = pawset%aeloc(1:pawset%grid%mesh)
#else
CALL step_f( aevloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh)
aux(1:pawset%grid%mesh) = pawset%psloc(1:pawset%grid%mesh)
CALL step_f( psvloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh)
DEALLOCATE (aux)
#endif
#endif
do ir = 1, rgrid(is)%mesh

0
atomic/vxc_t.f90 → Modules/vxc_t.f90
View File

10
atomic/vxcgc.f90 → Modules/vxcgc.f90
View File

@ -21,7 +21,6 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc,iflag)
use kinds, only : DP
use constants, only : fpi
use funct, only : gcxc, gcx_spin, gcc_spin
use ld1inc, only : grid
implicit none
integer, intent(in) :: ndm,mesh,nspin,iflag
real(DP), intent(in) :: r(mesh), r2(mesh), rho(ndm,2), rhoc(ndm)
@ -52,7 +51,7 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc,iflag)
do i=1, mesh
rhoaux(i,is)=(rho(i,is)+rhoc(i)/nspin)/fpi/r2(i)
enddo
call gradient(rhoaux(1,is),grho(1,is),r,mesh,iflag)
call radial_gradient(rhoaux(1,is),grho(1,is),r,mesh,iflag)
enddo
if (nspin.eq.1) then
@ -121,7 +120,7 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc,iflag)
! and correlation potential.
!
do is=1,nspin
call gradient(h(1,is),dh,r,mesh,iflag)
call radial_gradient(h(1,is),dh,r,mesh,iflag)
!
! Finally we compute the total exchange and correlation energy and
! potential. We put the original values on the charge and multiply
@ -144,7 +143,7 @@ subroutine vxcgc(ndm,mesh,nspin,r,r2,rho,rhoc,vgc,egc,iflag)
return
end subroutine vxcgc
subroutine gradient(f,gf,r,mesh,iflag)
subroutine radial_gradient(f,gf,r,mesh,iflag)
!
! This subroutine calculates the derivative with respect to r of a
! radial function defined on the mesh r. If iflag=0 it uses all mesh
@ -153,7 +152,6 @@ subroutine gradient(f,gf,r,mesh,iflag)
! is too smooth.
!
use kinds, only : DP
use ld1inc, only: zed
use radial_grids, only : series
implicit none
integer, intent(in) :: mesh, iflag
@ -248,7 +246,7 @@ do i=1,imin
gf(i)=b(1)+r(i)*(b(2)+r(i)*(b(3)+r(i)*b(4)))
enddo
return
end subroutine gradient
end subroutine radial_gradient
subroutine fit_pol(xdata,ydata,n,degree,b)
!

3
PH/Makefile
View File

@ -152,6 +152,7 @@ MODULES = \
../Modules/descriptors.o \
../Modules/dspev_drv.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -198,6 +199,8 @@ MODULES = \
../Modules/upf_to_internal.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/wave_base.o \
../Modules/xml_io_base.o \

2
PH/add_for_charges.f90
View File

@ -80,7 +80,7 @@ subroutine add_for_charges (ik, uact)
ikk = ik
ikq = ik
else
call infomsg ('add_for_charges', 'called for lgamma .eq. false', -1)
call infomsg ('add_for_charges', 'called for lgamma .eq. false')
endif
if (lsda) current_spin = isk (ikk)
!

1
PH/make.depend
View File

@ -611,7 +611,6 @@ set_asr_c.o : ../Modules/kind.o
set_drhoc.o : ../Modules/atom.o
set_drhoc.o : ../Modules/ions_base.o
set_drhoc.o : ../Modules/kind.o
set_drhoc.o : ../Modules/radial_grids.o
set_drhoc.o : ../PW/pwcom.o
set_drhoc.o : phcom.o
set_dvscf.o : ../Modules/kind.o

3
PP/Makefile
View File

@ -50,6 +50,7 @@ MODULES = \
../Modules/control_flags.o \
../Modules/descriptors.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -84,6 +85,8 @@ MODULES = \
../Modules/upf_to_internal.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/wannier.o \
../Modules/xml_io_base.o \

3
PW/Makefile
View File

@ -274,6 +274,7 @@ MODULES = \
../Modules/descriptors.o \
../Modules/dspev_drv.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -320,6 +321,8 @@ MODULES = \
../Modules/upf_to_internal.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/wave_base.o \
../Modules/xml_io_base.o \

2
PW/c_bands.f90
View File

@ -39,7 +39,7 @@ SUBROUTINE c_bands( iter, ik_, dr2 )
INTEGER :: ik_, iter
! k-point already done
! current iterations
REAL(DP) :: dr2
REAL(DP),INTENT(IN) :: dr2
! current accuracy of self-consistency
!
! ... local variables

158
PW/electrons.f90
View File

@ -62,17 +62,19 @@ SUBROUTINE electrons()
USE mp_global, ONLY : intra_pool_comm, npool
USE mp, ONLY : mp_sum
!
!!PAW]
#ifdef __GRID_PAW
USE grid_paw_variables, ONLY : really_do_paw, okpaw, tpawp, &
ehart1, ehart1t, etxc1, etxc1t, deband_1ae, deband_1ps, &
descf_1ae, descf_1ps, rho1, rho1t, rho1new, rho1tnew, &
vr1, vr1t, becnew
USE grid_paw_routines, ONLY : compute_onecenter_potentials, &
compute_onecenter_charges, delta_e_1, delta_e_1scf
USE rad_paw_routines, ONLY : PAW_energy
USE uspp, ONLY : becsum
USE uspp_param, ONLY : nhm
!!PAW]
#else
USE grid_paw_variables, ONLY : really_do_paw, okpaw, tpawp, becnew
#endif
USE rad_paw_routines, ONLY : PAW_potential,PAW_integrate
USE uspp, ONLY : becsum ! used for PAW
USE uspp_param, ONLY : nhm ! used for PAW
!
IMPLICIT NONE
!
@ -111,16 +113,17 @@ SUBROUTINE electrons()
REAL(DP), EXTERNAL :: ewald, get_clock
!
! ... additional variables for PAW
REAL(DP), ALLOCATABLE :: becstep(:,:,:)
REAL(DP), ALLOCATABLE :: becstep(:,:,:) ! cross-band occupations, used for mixing
REAL (DP) :: deband_PAW, descf_PAW, e_PAW ! deband, descf and E corrections from PAW
REAL (DP) :: correction1c ! total PAW correction
#ifdef __GRID_PAW
INTEGER :: na
INTEGER,PARAMETER :: AE=1,PS=2,XC=1,H=2
REAL (DP) :: correction1c
REAL (DP), ALLOCATABLE :: deband_1ps_na(:), deband_1ae_na(:), & ! auxiliary info on
descf_1ps_na(:), descf_1ae_na(:) ! one-center corrections
REAL (DP) :: paw_correction(nat,2,2) ! {# of atoms}, {XC|H}, {AE|PS}
!
IF (okpaw) ALLOCATE ( deband_1ae_na(nat), deband_1ps_na(nat), &
descf_1ae_na(nat), descf_1ps_na(nat) )
#endif
!
iter = 0
ik_ = 0
@ -135,7 +138,9 @@ SUBROUTINE electrons()
!
IF ( output_drho /= ' ' ) CALL remove_atomic_rho ()
!
#ifdef __GRID_PAW
IF (okpaw) DEALLOCATE(deband_1ae_na, deband_1ps_na, descf_1ae_na, descf_1ps_na)
#endif
RETURN
!
END IF
@ -156,7 +161,9 @@ SUBROUTINE electrons()
!
conv_elec = .TRUE.
!
#ifdef __GRID_PAW
IF (okpaw) DEALLOCATE(deband_1ae_na, deband_1ps_na, descf_1ae_na, descf_1ps_na)
#endif
!
RETURN
!
@ -190,10 +197,11 @@ SUBROUTINE electrons()
!
IF (okpaw) THEN
IF ( .not. ALLOCATED(becstep) ) ALLOCATE (becstep(nhm*(nhm+1)/2,nat,nspin))
becstep (:,:,:) = 0.d0
DO na = 1, nat
IF (tpawp(ityp(na))) becstep(:,na,:) = becsum(:,na,:)
END DO
becstep (:,:,:) = 0._dp
! DO na = 1, nat
! IF (tpawp(ityp(na))) becstep(:,na,:) = becsum(:,na,:)
! END DO
becstep(:,:,:) = becsum(:,:,:)
END IF
!
DO idum = 1, niter
@ -354,19 +362,24 @@ SUBROUTINE electrons()
!
! ... update core occupations for PAW
IF (okpaw) THEN
#ifdef __GRID_PAW
CALL compute_onecenter_charges (becsum,rho1,rho1t)
!
deband_1ae = delta_e_1(rho1, vr1, deband_1ae_na) !AE
deband_1ps = delta_e_1(rho1t,vr1t,deband_1ps_na) !PS
!
CALL infomsg ('electrons','mixing several times ns if lda_plus_U')
deband_1ae = delta_e_1(rho1, vr1, deband_1ae_na) !AE
deband_1ps = delta_e_1(rho1t,vr1t,deband_1ps_na) !PS
write(6,"(a,2f12.6)") " == deband (grid PAW): ", deband_1ae-deband_1ps
#endif
deband_PAW = - PAW_integrate(becsum)
!write(6,"(a,2f13.8)") " == deband (radial PAW): ", deband_PAW
!CALL infomsg ('electrons','mixing several times ns if lda_plus_U')
IF (lda_plus_U) STOP 'electrons - not implemented'
!
ALLOCATE (becnew(nhm*(nhm+1)/2, nat, nspin) )
becnew(:,:,:) = 0.d0
DO na = 1, nat
IF (tpawp(ityp(na))) becnew(:,na,:) = becsum(:,na,:)
END DO
becnew(:,:,:) = 0._DP
! DO na = 1, nat
! IF (tpawp(ityp(na))) becnew(:,na,:) = becsum(:,na,:)
! END DO
becnew(:,:,:) = becsum(:,:,:)
END IF
!
CALL mix_rho( rhognew, rhog, taukgnew, taukg, becnew, becstep, &
@ -417,6 +430,7 @@ SUBROUTINE electrons()
!
END IF
!
not_converged_electrons : &
IF ( .NOT. conv_elec ) THEN
!
! ... bring mixed charge density (rhog) from G- to R-space (rhonew)
@ -484,58 +498,33 @@ SUBROUTINE electrons()
!
! ... compute PAW corrections to descf
IF (okpaw) THEN
! ... radial PAW:
CALL PAW_energy(becstep,paw_correction)
!
! ... grid PAW:
#ifdef __GRID_PAW
ALLOCATE (rho1new (nrxx,nspin,nat), rho1tnew(nrxx,nspin,nat) )
!
CALL compute_onecenter_charges (becstep,rho1new,rho1tnew)
WRITE(*,*) "== rho max diff: ", MAXVAL(ABS(rho1(:,:,:)-rho1t(:,:,:)))
CALL compute_onecenter_potentials(becstep,rho1new,rho1tnew)
!
descf_1ae = delta_e_1scf(rho1, rho1new, vr1, descf_1ae_na) ! AE
descf_1ps = delta_e_1scf(rho1t,rho1tnew,vr1t,descf_1ps_na) ! PS
!#define __PAW_REPORT
#ifdef __PAW_REPORT
WRITE(6,"(a,f15.10)") "==PAW (HARTREE): GRID RAD"
DO i = 1,nat
WRITE(6,"(a,i2,a,2f15.10)") "==AE",i," :", ehart1(i),paw_correction(i,H,AE)
WRITE(6,"(a,i2,a,2f15.10)") "==PS",i," :", ehart1t(i),paw_correction(i,H,PS)
WRITE(6,"(a,i2,a,2f15.10)") "==AE-PS",i," :", ehart1(i)-ehart1t(i),&
paw_correction(i,H,AE)-paw_correction(i,H,PS)
ENDDO
WRITE(6,"(a,2f15.10)") "==AE tot :", SUM(ehart1(:)),SUM(paw_correction(:,H,AE))
WRITE(6,"(a,2f15.10)") "==PS tot :", SUM(ehart1t(:)),SUM(paw_correction(:,H,PS))
WRITE(6,"(a,2f15.10)") "==AE-PS tot:", SUM(ehart1(:))-SUM(ehart1t(:)),&
SUM(paw_correction(:,H,AE))-SUM(paw_correction(:,H,PS))
WRITE(6,"(a,2f15.10)") "=="
!
WRITE(6,"(a,2f15.10)") "==PAW (XC): GRID RAD"
DO i = 1,nat
WRITE(6,"(a,i2,a,2f15.10)") "==AE",i," :", etxc1(i),paw_correction(i,XC,AE)
WRITE(6,"(a,i2,a,2f15.10)") "==PS",i," :", etxc1t(i),paw_correction(i,XC,PS)
WRITE(6,"(a,i2,a,2f15.10)") "==AE-PS",i," :", etxc1(i)-etxc1t(i),&
paw_correction(i,XC,AE)-paw_correction(i,XC,PS)
ENDDO
WRITE(6,"(a,2f15.10)") "==AE tot :", SUM(etxc1(:)),SUM(paw_correction(:,XC,AE))
WRITE(6,"(a,2f15.10)") "==PS tot :", SUM(etxc1t(:)),SUM(paw_correction(:,XC,PS))
WRITE(6,"(a,2f15.10)") "==AE-PS tot:", SUM(etxc1(:))-SUM(etxc1t(:)),&
SUM(paw_correction(:,XC,AE))-SUM(paw_correction(:,XC,PS))
WRITE(6,"(a,2f15.10)") "=="
WRITE(6,"(a,2f15.10)") "==TOTAL:", SUM(etxc1(:))-SUM(etxc1t(:))+SUM(ehart1(:))-SUM(ehart1t(:)),&
SUM(paw_correction(:,XC,AE))-SUM(paw_correction(:,XC,PS))+&
SUM(paw_correction(:,H,AE))-SUM(paw_correction(:,H,PS))
#endif
!
DEALLOCATE (rho1new, rho1tnew)
DO na = 1,nat
write(*,"(a,2i2,2f10.5,a)") " == ",1,na, ehart1(na), etxc1(na),"(gu)"
write(*,"(a,2i2,2f10.5,a)") " == ",2,na, ehart1t(na), etxc1t(na),"(gu)"
ENDDO
#endif
CALL PAW_potential(becstep, e_PAW)
descf_PAW = - PAW_integrate(becstep-becsum)
!write(6,"(a,2f12.6)") " == descf: ", descf_1ae-descf_1ps, descf_PAW
END IF
!
! ... write the charge density to file
!
CALL write_rho( rho, nspin )
!
ELSE
ELSE not_converged_electrons
!
! ... convergence reached:
! ... 1) the output HXC-potential is saved in vr
@ -550,25 +539,35 @@ SUBROUTINE electrons()
!
! CHECKME: is it becsum or becstep??
IF (okpaw) THEN
CALL PAW_energy(becstep, paw_correction)
#ifdef __GRID_PAW
CALL compute_onecenter_potentials(becsum,rho1,rho1t)
CALL infomsg ('electrons','PAW forces missing')
DO na = 1,nat
write(*,"(a,2i2,2f10.5,a)") " == ",1,na, ehart1(na), etxc1(na),"(g)"
write(*,"(a,2i2,2f10.5,a)") " == ",2,na, ehart1t(na), etxc1t(na),"(g)"
ENDDO
#endif
CALL PAW_potential(becsum, e_PAW)
!write(*,"(a,f13.8)") " == energy (radial PAW)",e_PAW
!CALL infomsg ('electrons','PAW forces missing')
ENDIF
!
! ... note that rho is the output, not mixed, charge density
! ... so correction for variational energy is no longer needed
!
descf = 0.D0
descf = 0._dp
!
#ifdef __GRID_PAW
IF (okpaw) THEN
descf_1ae=0._DP
descf_1ps=0._DP
descf_1ae_na(:)=0.d0
descf_1ps_na(:)=0.d0
descf_1ae=0._dp
descf_1ps=0._dp
descf_1ae_na(:)=0._dp
descf_1ps_na(:)=0._dp
ENDIF
#endif
IF (okpaw) descf_PAW = 0._dp
!
!
END IF
END IF not_converged_electrons
!
#if defined (EXX)
IF ( exx_is_active() ) THEN
@ -652,18 +651,29 @@ SUBROUTINE electrons()
etot = eband + ( etxc - etxcc ) + ewld + ehart + deband + demet + descf
!
IF (okpaw) THEN
WRITE(stdout,'(A,F15.8)') 'US energy before PAW additions', etot
#ifdef __GRID_PAW
DO na = 1, nat
IF (tpawp(ityp(na))) THEN
correction1c = ehart1(na) -ehart1t(na) +etxc1(na) -etxc1t(na) + &
deband_1ae_na(na) - deband_1ps_na(na) + &
descf_1ae_na(na) - descf_1ps_na(na)
PRINT '(4x,A,I3,A,F15.8)', 'atom #: ', na, ' grid correction: ', correction1c
PRINT '(4x,A,I3,A,F15.8)', 'atom #: ', na, ' rad correction: ', &
SUM(paw_correction(na,:,AE)-paw_correction(na,:,PS))
IF (really_do_paw) etot = etot + correction1c
deband_1ae_na(na) - deband_1ps_na(na) + &
descf_1ae_na(na) - descf_1ps_na(na)
PRINT *, "== "
PRINT '(A,I1,A,F13.8)', ' == atom #: ', na, ' grid correction: ', correction1c
PRINT '(A,4f13.8)', ' == ',deband_1ae_na(na) - deband_1ps_na(na),&
descf_1ae_na(na) - descf_1ps_na(na),&
ehart1(na) -ehart1t(na) +etxc1(na) -etxc1t(na)
PRINT *, "== == == == == == =="
END IF
END DO
#endif
correction1c = (deband_PAW + descf_PAW + e_PAW)
PRINT '(5x,A,f10.6,A)', 'PAW correction: ',correction1c, ', composed of: '
PRINT '(5x,A,f10.6,A,f10.6,A,f10.6)',&
'de_band: ',deband_PAW,', de_scf: ',descf_PAW,', 1-center E: ', e_PAW
IF (really_do_paw) THEN
etot = etot + correction1c
hwf_energy = hwf_energy + correction1c
ENDIF
END IF
!
#if defined (EXX)
@ -736,6 +746,7 @@ SUBROUTINE electrons()
IF ( tefield ) WRITE( stdout, 9061 ) etotefield
IF ( lda_plus_u ) WRITE( stdout, 9065 ) eth
IF ( ABS (descf) > eps8 ) WRITE( stdout, 9069 ) descf
IF( okpaw) WRITE( stdout, 9067 ) correction1c
!
! ... With Fermi-Dirac population factor, etot is the electronic
! ... free energy F = E - TS , demet is the -TS contribution
@ -804,7 +815,9 @@ SUBROUTINE electrons()
!DEALLOCATE (rhog) |should be elsewhere
IF (okpaw) THEN
DEALLOCATE (becstep)
#ifdef __GRID_PAW
DEALLOCATE(deband_1ae_na, deband_1ps_na, descf_1ae_na, descf_1ps_na)
#endif
END IF
!
RETURN
@ -849,8 +862,9 @@ SUBROUTINE electrons()
9062 FORMAT( ' Fock energy 1 =',F15.8,' Ry' )
9063 FORMAT( ' Fock energy 2 =',F15.8,' Ry' )
9064 FORMAT( ' Half Fock energy 2 =',F15.8,' Ry' )
9066 FORMAT( ' dexx =',F15.8,' Ry' )
9065 FORMAT( ' Hubbard energy =',F15.8,' Ry' )
9066 FORMAT( ' dexx =',F15.8,' Ry' )
9067 FORMAT( ' one-center paw contrib. =',F15.8,' Ry' )
9069 FORMAT( ' scf correction =',F15.8,' Ry' )
9070 FORMAT( ' smearing contrib. (-TS) =',F15.8,' Ry' )
9071 FORMAT( ' Magnetic field =',3F12.7,' Ry' )

336
PW/grid_paw_routines.f90
View File

@ -5,6 +5,12 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!!! *****************************************************
!!! WARNING! *
!!! This module contains only two functions unless the *
!!! flag __GRID_PAW is specified at compilation time! *
!!! *****************************************************
#include "f_defs.h"
!
MODULE grid_paw_routines
@ -13,16 +19,73 @@ MODULE grid_paw_routines
!
! NO spin-orbit
! NO EXX
! NO Parallelism
! NO rinner > 0
! NO Gamma (?)
!
IMPLICIT NONE
PUBLIC! <===
SAVE
!!!=========================================================================
CONTAINS
CONTAINS
! Initialize becsum with atomic occupations (for PAW atoms only)
! Notice: requires exact correspondence chi <--> beta in the atom,
! that is that all wavefunctions considered for PAW generation are
! counted in chi (otherwise the array "oc" does not correspond to beta)
SUBROUTINE atomic_becsum()
USE kinds, ONLY : DP
USE uspp, ONLY : becsum, nhtol, indv
USE uspp_param, ONLY : nh
USE ions_base, ONLY : nat, ityp
USE lsda_mod, ONLY : nspin, starting_magnetization
USE atom, ONLY : oc
USE grid_paw_variables, ONLY : tpawp, okpaw
IMPLICIT NONE
INTEGER :: ispin, na, nt, ijh, ih, jh, nb, mb
!
IF (.NOT. okpaw) RETURN
!
if (nspin.GT.2) STOP 'atomic_becsum not implemented'
!
na_loop: DO na = 1, nat
nt = ityp(na)
is_paw: IF (tpawp(nt)) THEN
!
ijh = 1
ih_loop: DO ih = 1, nh(nt)
nb = indv(ih,nt)
!
IF (nspin==1) THEN
!
becsum(ijh,na,1) = oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
!
ELSE IF (nspin==2) THEN
!
becsum(ijh,na,1) = 0.5d0 * (1.d0 + starting_magnetization(nt))* &
oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
becsum(ijh,na,2) = 0.5d0 * (1.d0 - starting_magnetization(nt))* &
oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
!
END IF
ijh = ijh + 1
!
jh_loop: DO jh = ( ih + 1 ), nh(nt)
!mb = indv(jh,nt)
DO ispin = 1, nspin
becsum(ijh,na,ispin) = 0._DP
END DO
ijh = ijh + 1
!
END DO jh_loop
END DO ih_loop
END IF is_paw
END DO na_loop
#if defined __DEBUG_ATOMIC_BECSUM
PRINT '(1f20.10)', becsum(:,1,1)
#endif
END SUBROUTINE atomic_becsum
! Analogous to PW/allocate_nlpot.f90
SUBROUTINE allocate_paw_internals
@ -33,35 +96,41 @@ CONTAINS
USE uspp_param, ONLY : lmaxq, nhm, nbetam
USE gvect, ONLY : ngl
!
USE grid_paw_variables, ONLY : pp, ppt, prad, ptrad, rho1, rho1t, &
vr1, vr1t, int_r2pfunc, int_r2ptfunc, ehart1, etxc1, vtxc1, &
ehart1t, etxc1t, vtxc1t, aerho_core, psrho_core, radial_distance, radial_r,&
dpaw_ae, dpaw_ps, aevloc, psvloc, aevloc_r, psvloc_r, &
prodp, prodpt, prod0p, prod0pt
! USE grid_paw_variables, ONLY : pp, ppt, prad, ptrad, rho1, rho1t, &
! vr1, vr1t, int_r2pfunc, int_r2ptfunc, ehart1, etxc1, vtxc1, &
! ehart1t, etxc1t, vtxc1t, aerho_core, psrho_core, radial_distance, radial_r,&
! dpaw_ae, dpaw_ps, aevloc, psvloc, aevloc_r, psvloc_r, &
! prodp, prodpt, prod0p, prod0pt
USE grid_paw_variables
!
IMPLICIT NONE
call start_clock('allocate_paw')
!
#ifdef __GRID_PAW
ALLOCATE (pp( nhm, nhm, nsp))
ALLOCATE (ppt( nhm, nhm, nsp))
!
ALLOCATE (int_r2pfunc( nhm, nhm, nsp))
ALLOCATE (int_r2ptfunc( nhm, nhm, nsp))
ALLOCATE (int_r2pfunc( nhm, nhm, nsp))
ALLOCATE (int_r2ptfunc( nhm, nhm, nsp))
!
IF (lmaxq > 0) ALLOCATE (prad( nqxq, nbetam*(nbetam+1)/2, lmaxq, nsp))
IF (lmaxq > 0) ALLOCATE (ptrad( nqxq, nbetam*(nbetam+1)/2, lmaxq, nsp))
!
ALLOCATE (aevloc( ngl, ntyp))
#endif
ALLOCATE (psvloc( ngl, ntyp))
#ifdef __GRID_PAW
!
ALLOCATE (aevloc_r(nrxx,ntyp))
ALLOCATE (psvloc_r(nrxx,ntyp))
ALLOCATE (radial_distance(nrxx))
ALLOCATE (radial_r(3,nrxx))
!
ALLOCATE (prodp(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
ALLOCATE (prodpt(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
ALLOCATE (prod0p(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
ALLOCATE (prod0pt(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
! ALLOCATE (prodp(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
! ALLOCATE (prodpt(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
! ALLOCATE (prod0p(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
! ALLOCATE (prod0pt(nhm*(nhm+1)/2,nhm*(nhm+1)/2,ntyp))
!
ALLOCATE (rho1(nrxx, nspin, nat))
ALLOCATE (rho1t(nrxx, nspin, nat))
@ -77,10 +146,12 @@ CONTAINS
ALLOCATE (vtxc1t (nat))
ALLOCATE (aerho_core(nrxx,ntyp))
ALLOCATE (psrho_core(nrxx,ntyp))
#endif
!
ALLOCATE(dpaw_ae( nhm, nhm, nat, nspin))
ALLOCATE(dpaw_ps( nhm, nhm, nat, nspin))
!
!
END SUBROUTINE allocate_paw_internals
! Called from clean_pw
@ -92,15 +163,17 @@ CONTAINS
USE uspp_param, ONLY : lmaxq, nhm, nbetam
USE gvect, ONLY : ngl
!
USE grid_paw_variables, ONLY : pp, ppt, prad, ptrad, rho1, rho1t, &
vr1, vr1t, int_r2pfunc, int_r2ptfunc, ehart1, etxc1, vtxc1, &
ehart1t, etxc1t, vtxc1t, aerho_core, psrho_core, radial_distance, radial_r,&
dpaw_ae, dpaw_ps, aevloc, psvloc, aevloc_r, psvloc_r, &
prodp, prodpt, prod0p, prod0pt, aug
! USE grid_paw_variables, ONLY : pp, ppt, prad, ptrad, rho1, rho1t, &
! vr1, vr1t, int_r2pfunc, int_r2ptfunc, ehart1, etxc1, vtxc1, &
! ehart1t, etxc1t, vtxc1t, aerho_core, psrho_core, radial_distance, radial_r,&
! dpaw_ae, dpaw_ps, aevloc, psvloc, aevloc_r, psvloc_r, &
! prodp, prodpt, prod0p, prod0pt, aug
USE grid_paw_variables
!
IMPLICIT NONE
INTEGER :: nt
!
#ifdef __GRID_PAW
IF(allocated(pp)) DEALLOCATE (pp)
IF(allocated(ppt)) DEALLOCATE (ppt)
!
@ -111,8 +184,10 @@ CONTAINS
IF(allocated(ptrad)) DEALLOCATE (ptrad)
!
IF(allocated(aevloc)) DEALLOCATE (aevloc)
#endif
IF(allocated(psvloc)) DEALLOCATE (psvloc)
!
#ifdef __GRID_PAW
IF(allocated(aevloc_r)) DEALLOCATE (aevloc_r)
IF(allocated(psvloc_r)) DEALLOCATE (psvloc_r)
IF(allocated(radial_distance)) &
@ -138,12 +213,67 @@ CONTAINS
IF(allocated(vtxc1t )) DEALLOCATE (vtxc1t )
IF(allocated(aerho_core)) DEALLOCATE (aerho_core)
IF(allocated(psrho_core)) DEALLOCATE (psrho_core)
#endif
!
IF(allocated(dpaw_ae)) DEALLOCATE (dpaw_ae)
IF(allocated(dpaw_ps)) DEALLOCATE (dpaw_ps)
!
DO nt = 1,ntyp
IF(allocated(aug(nt)%fun)) DEALLOCATE (aug(nt)%fun)
ENDDO
!
END SUBROUTINE deallocate_paw_internals
#ifdef __GRID_PAW
! Analogous to PW/init_vloc.f90
SUBROUTINE init_paw_vloc
!
USE kinds, ONLY: DP
USE atom, ONLY : rgrid, msh
USE ions_base, ONLY : ntyp => nsp
USE cell_base, ONLY : omega, tpiba2
USE gvect, ONLY : ngl, gl
USE uspp_param, ONLY : upf
USE grid_paw_variables, ONLY: aevloc_at, psvloc_at, aevloc, psvloc
!
IMPLICIT NONE
!
INTEGER :: nt
! counter on atomic types
!
REAL(DP), POINTER :: vloc_at_(:,:)
REAL(DP), POINTER :: vloc_(:,:)
INTEGER :: i_what
!
call start_clock ('init_pawvloc')
whattodo: DO i_what=1, 2
NULLIFY(vloc_at_,vloc_)
IF (i_what==1) THEN
vloc_at_ => aevloc_at
vloc_ => aevloc
ELSE IF (i_what==2) THEN
vloc_at_ => psvloc_at
vloc_ => psvloc
END IF
! associate a pointer to the AE or PS part
vloc_(:,:) = 0.d0
DO nt = 1, ntyp
!
! compute V_loc(G) for a given type of atom
!
CALL vloc_of_g_noerf (rgrid(nt)%mesh, &
msh (nt), rgrid(nt)%rab(1), rgrid(nt)%r(1), vloc_at_ (:, nt), &
upf(nt)%zp, tpiba2, ngl, gl, omega, vloc_ (:, nt) )
END DO
END DO whattodo
call stop_clock ('init_pawvloc')
!
RETURN
!
END SUBROUTINE init_paw_vloc
!!$=========================================================================
! Analogous to part of PW/init_us_1.f90
! Notice integration performed not only up to r(kkbeta(:)) but up to r(msh(:))
@ -202,8 +332,6 @@ CONTAINS
!
whattodo: DO i_what=1, 2
! associate a pointer to the AE or PS part
write (*,*) "entering init_prad"
write (*,*) lmaxq, lqmax
NULLIFY(pfunc_,pp_,prad_)
IF (i_what==1) THEN
pfunc_=> pfunc
@ -629,6 +757,8 @@ CONTAINS
nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
nl, ngm, g, nspin, alat, omega, &
etxc1_(na), vtxc1_(na), vr1_(:,:,na) )
#ifdef __NO_HARTREE
#else
!
spheropole=0.d0
DO nt = 1, ntyp
@ -649,7 +779,7 @@ CONTAINS
END DO
END DO
END DO
write (*,*) "calculate tot_multipoles", lqmax
!write (*,*) "calculate tot_multipoles", lqmax
tot_multipole(:)=0.d0
DO nt = 1, ntyp
IF (ityp(na)==nt) THEN
@ -682,13 +812,15 @@ CONTAINS
END DO
END IF
END DO
write (*,*) "done"
!write (*,*) "done"
!
!PRINT *, 'SPHEROPOLE:',na, spheropole
CALL v_h_grid( rho1_(:,:,na), nr1,nr2,nr3, nrx1,nrx2,nrx3, nrxx, &
nl, ngm, gg, gstart, nspin, alat, omega, &
ehart1_(na), charge, vr1_(:,:,na), &
spheropole, tot_multipole, na)
#endif
END IF
END DO
END DO whattodo
@ -878,6 +1010,10 @@ CONTAINS
INTEGER :: ir, & ! counter on mesh points
nt ! counter on atomic types
WRITE(*,*) "**************************************"
WRITE(*,*) "* set_paw_rhoc: I SHOULD NOT BE HERE *"
WRITE(*,*) "**************************************"
!
call start_clock ('set_paw_rhoc')
ALLOCATE (aux(nrxx), rhocg(ngl))
@ -988,6 +1124,11 @@ SUBROUTINE v_h_grid( rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
REAL(DP) :: dummyx, c(3), r2, gaussrho
INTEGER :: ir1,ir2,ir3, i,j,k
real(DP), external :: erf
WRITE(*,*) "**********************************"
WRITE(*,*) "* v_h_grid: I SHOULD NOT BE HERE *"
WRITE(*,*) "**********************************"
!
call start_clock ('v_h_grid')
CALL infomsg ('v_h_grid','alpha set manually')
@ -1015,6 +1156,7 @@ SUBROUTINE v_h_grid( rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
charge = charge * omega / (nr1*nr2*nr3)
dipole(:) = dipole(:) * omega / (nr1*nr2*nr3)
quad(:) = quad(:) * omega / (nr1*nr2*nr3)
#ifdef __PRINT_MULTIPOLES
write (*,'(a,2f10.5)') 'charge=', charge, tot_multipole(1)*sqrt(fpi)
write (*,'(a,2f10.5)') 'dipole=', dipole(1),-tot_multipole(3)*sqrt(fpi/3.0_DP)
write (*,'(a,2f10.5)') 'dipole=', dipole(2),-tot_multipole(4)*sqrt(fpi/3.0_DP)
@ -1024,6 +1166,7 @@ SUBROUTINE v_h_grid( rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
write (*,'(a,2f10.5)') 'quad =', quad(3), tot_multipole(7)*sqrt(fpi/5.0_DP)
write (*,'(a,2f10.5)') 'quad =', quad(4), tot_multipole(8)*sqrt(fpi/5.0_DP)
write (*,'(a,2f10.5)') 'quad =', quad(5), tot_multipole(9)*sqrt(fpi/5.0_DP)
#endif
!
! ... bring rho (aux) to G space
!
@ -1035,7 +1178,7 @@ SUBROUTINE v_h_grid( rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
!
CALL reduce( 1, charge )
!
PRINT *, 'charge=', charge, tot_multipole(1)*sqrt(fpi)
!PRINT *, 'charge=', charge, tot_multipole(1)*sqrt(fpi)
!
! ... calculate hartree potential in G-space (NB: only G/=0 )
!
@ -1098,7 +1241,7 @@ SUBROUTINE v_h_grid( rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
! Add analytic contribution from gaussian charges to Hartree energy
ehart = ehart + e2 * (charge**2 + p2*alpha / 3.0_DP) * SQRT(alpha/TPI)
edip_r = e2 * p2*alpha / 3.0_DP * SQRT(alpha/TPI)
PRINT '(A,4f12.8)','EDIP',edip_g0,edip_g,edip_r, edip_g0+edip_g+edip_r
!PRINT '(A,4f12.8)','EDIP',edip_g0,edip_g,edip_r, edip_g0+edip_g+edip_r
!
#if defined __DEBUG_NEWD_PAW_GRID
! PRINT '(A,3f20.10)', 'SPHEROPOLE,CHARGE: ', charge, &
@ -1255,17 +1398,23 @@ SUBROUTINE newd_paw_grid(na)
INTEGER,INTENT(IN) :: na
INTEGER :: i, ir, is, nir
REAL(DP) :: deltav, average, rms
WRITE(*,*) "***************************************"
WRITE(*,*) "* newd_paw_grid: I SHOULD NOT BE HERE *"
WRITE(*,*) "***************************************"
!
IF ( .NOT. okpaw ) RETURN
call start_clock ('newd_paw')
!
!PRINT '(A)', 'WARNING newd_paw_grid contains only H+xc potentials'
!
! write(6,"(a)") "Now AE part:"
CALL integrate_potential_x_charge (prad, vr1, aevloc_r, dpaw_ae, &
SIZE(prad,1),SIZE(prad,2),SIZE(prad,3),SIZE(prad,4), &
SIZE(vr1,1),SIZE(vr1,2),SIZE(vr1,3), &
SIZE(aevloc_r,1),SIZE(aevloc_r,2), &
SIZE(dpaw_ae,1),SIZE(dpaw_ae,2),SIZE(dpaw_ae,3),SIZE(dpaw_ae,4))
! write(6,"(a)") "Now PS part:"
CALL integrate_potential_x_charge (ptrad, vr1t, psvloc_r, dpaw_ps, &
SIZE(prad,1),SIZE(prad,2),SIZE(prad,3),SIZE(prad,4), &
SIZE(vr1,1),SIZE(vr1,2),SIZE(vr1,3), &
@ -1273,10 +1422,11 @@ SUBROUTINE newd_paw_grid(na)
SIZE(dpaw_ae,1),SIZE(dpaw_ae,2),SIZE(dpaw_ae,3),SIZE(dpaw_ae,4))
call stop_clock ('newd_paw')
!
CONTAINS
CONTAINS
!
SUBROUTINE integrate_potential_x_charge (prad_, vr_, vl_, dpaw_, &
sp1,sp2,sp3,sp4, sv1,sv2,sv3, sl1,sl2, sd1,sd2,sd3,sd4)
#define __DEBUG_NEWD_PAW_GRID
USE kinds, ONLY : DP
USE ions_base, ONLY : nat, ntyp => nsp, ityp
USE cell_base, ONLY : omega
@ -1307,10 +1457,15 @@ CONTAINS
! spherical harmonics, modulus of G
REAL(DP) :: fact, DDOT
!
nt = ityp(na)
! PRINT *, 'BEFORE:'
! PRINT '(8f15.7)', ((dpaw_(jh,ih,na,1),jh=1,nh(nt)),ih=1,nh(nt))
!
IF ( gamma_only ) THEN
fact = 2.D0
fact = 2._DP
ELSE
fact = 1.D0
fact = 1._DP
END IF
!
ALLOCATE( aux(ngm,nspin), qgm(ngm), qmod(ngm), ylmk0(ngm,lmaxq*lmaxq) )
@ -1321,14 +1476,17 @@ CONTAINS
!
qmod(1:ngm) = SQRT( gg(1:ngm) )
!
! LOOP on NA removed to complain with new newd.f0 structure!!
! LOOP on NA removed to complain with new newd.f90 structure!!
!!$ DO na=1, nat
! The type is fixed
nt = ityp(na)
!
! ... fourier transform of the total effective potential
DO is = 1, nspin
psic(:) = vl_(:,nt) + vr_(:,is,na)
#ifdef __NO_LOCAL
psic(:) = vr_(:,is,na)
#else
psic(:) = vr_(:,is,na) + vl_(:,nt)
#endif
CALL cft3( psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1 )
aux(1:ngm,is) = psic( nl(1:ngm) )
END DO
@ -1374,8 +1532,8 @@ CONTAINS
CALL reduce( nhm * nhm * nat * nspin, dpaw_ )
!
#if defined __DEBUG_NEWD_PAW_GRID
! PRINT *, 'D - D1 or D1~'
! PRINT '(8f15.7)', ((dpaw_(jh,ih,1,1),jh=1,nh(nt)),ih=1,nh(nt))
! PRINT *, 'AFTER:'
! PRINT '(8f15.7)', ((dpaw_(jh,ih,na,1),jh=1,nh(nt)),ih=1,nh(nt))
#endif
!
DEALLOCATE( aux, qgm, qmod, ylmk0 )
@ -1384,116 +1542,6 @@ CONTAINS
END SUBROUTINE newd_paw_grid
! Initialize becsum with atomic occupations (for PAW atoms only)
! Notice: requires exact correspondence chi <--> beta in the atom,
! that is that all wavefunctions considered for PAW generation are
! counted in chi (otherwise the array "oc" does not correspond to beta)
!#define __DEBUG_ATOMIC_BECSUM
SUBROUTINE atomic_becsum()
USE kinds, ONLY : DP
USE uspp, ONLY : becsum, nhtol, indv
USE uspp_param, ONLY : nh
USE ions_base, ONLY : nat, ityp
USE lsda_mod, ONLY : nspin, starting_magnetization
USE atom, ONLY : oc
USE grid_paw_variables, ONLY : tpawp, okpaw
IMPLICIT NONE
INTEGER :: ispin, na, nt, ijh, ih, jh, nb, mb
!
IF (.NOT. okpaw) RETURN
!
if (nspin.GT.2) STOP 'atomic_becsum not implemented'
!
na_loop: DO na = 1, nat
nt = ityp(na)
is_paw: IF (tpawp(nt)) THEN
!
ijh = 1
ih_loop: DO ih = 1, nh(nt)
nb = indv(ih,nt)
!
#if defined __DEBUG_ATOMIC_BECSUM
PRINT *, ijh,ih,nb,oc(nb,nt),nhtol(ih,nt)
#endif
IF (nspin==1) THEN
!
becsum(ijh,na,1) = oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
!
ELSE IF (nspin==2) THEN
!
becsum(ijh,na,1) = 0.5d0 * (1.d0 + starting_magnetization(nt))* &
oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
becsum(ijh,na,2) = 0.5d0 * (1.d0 - starting_magnetization(nt))* &
oc(nb,nt) / REAL(2*nhtol(ih,nt)+1,DP)
!
END IF
ijh = ijh + 1
!
jh_loop: DO jh = ( ih + 1 ), nh(nt)
!mb = indv(jh,nt)
DO ispin = 1, nspin
becsum(ijh,na,ispin) = 0._DP
END DO
ijh = ijh + 1
!
END DO jh_loop
END DO ih_loop
END IF is_paw
END DO na_loop
#if defined __DEBUG_ATOMIC_BECSUM
PRINT '(1f20.10)', becsum(:,1,1)
#endif
END SUBROUTINE atomic_becsum
! Analogous to PW/init_vloc.f90
SUBROUTINE init_paw_vloc
!
USE kinds, ONLY: DP
USE atom, ONLY : rgrid, msh
USE ions_base, ONLY : ntyp => nsp
USE cell_base, ONLY : omega, tpiba2
USE gvect, ONLY : ngl, gl
USE uspp_param, ONLY : upf
USE grid_paw_variables, ONLY: aevloc_at, psvloc_at, aevloc, psvloc
!
IMPLICIT NONE
!
INTEGER :: nt
! counter on atomic types
!
REAL(DP), POINTER :: vloc_at_(:,:)
REAL(DP), POINTER :: vloc_(:,:)
INTEGER :: i_what
!
call start_clock ('init_pawvloc')
whattodo: DO i_what=1, 2
! associate a pointer to the AE or PS part
NULLIFY(vloc_at_,vloc_)
IF (i_what==1) THEN
vloc_at_ => aevloc_at
vloc_ => aevloc
ELSE IF (i_what==2) THEN
vloc_at_ => psvloc_at
vloc_ => psvloc
END IF
vloc_(:,:) = 0.d0
DO nt = 1, ntyp
!
! compute V_loc(G) for a given type of atom
!
CALL vloc_of_g_noerf (rgrid(nt)%mesh, &
msh (nt), rgrid(nt)%rab(1), rgrid(nt)%r(1), vloc_at_ (:, nt), &
upf(nt)%zp, tpiba2, ngl, gl, omega, vloc_ (:, nt) )
END DO
END DO whattodo
call stop_clock ('init_pawvloc')
!
RETURN
!
END SUBROUTINE init_paw_vloc
! Adapted from vloc_of_g.f90, for bounded potentials: don't add erf(r)/r
! because we don't assume anymore that v(r)\approx 2*e^2*zp/r at large r
@ -1748,4 +1796,6 @@ END FUNCTION delta_e_1
!
END FUNCTION delta_e_1scf
#endif
END MODULE grid_paw_routines

10
PW/hinit0.f90
View File

@ -22,7 +22,9 @@ SUBROUTINE hinit0()
USE wvfct, ONLY : npw, g2kin, igk
USE io_files, ONLY : iunigk
USE realus, ONLY : tqr, qpointlist
#ifdef __GRID_PAW
USE grid_paw_routines, ONLY : init_prad, set_paw_rhoc, init_paw_vloc, paw_grid_setlocal
#endif
USE grid_paw_variables, ONLY : okpaw
!
IMPLICIT NONE
@ -33,12 +35,16 @@ SUBROUTINE hinit0()
! ... calculate the local part of the pseudopotentials
!
CALL init_vloc()
#ifdef __GRID_PAW
IF (okpaw) CALL init_paw_vloc() !!PAW!!
#endif
!
! ... k-point independent parameters of non-local pseudopotentials
!
CALL init_us_1()
#ifdef __GRID_PAW
IF (okpaw) CALL init_prad() !!PAW!!
#endif
CALL init_at_1()
!
REWIND( iunigk )
@ -83,12 +89,16 @@ SUBROUTINE hinit0()
! ... calculate the total local potential
!
CALL setlocal()
#ifdef __GRID_PAW
IF (okpaw) CALL paw_grid_setlocal() !!PAW!!
#endif
!
! ... calculate the core charge (if any) for the nonlinear core correction
!
CALL set_rhoc()
#ifdef __GRID_PAW
IF (okpaw) CALL set_paw_rhoc() !!PAW!!
#endif
!
IF ( tqr ) CALL qpointlist()
!

2
PW/init_run.f90
View File

@ -15,6 +15,7 @@ SUBROUTINE init_run()
USE dynamics_module, ONLY : allocate_dyn_vars
USE grid_paw_routines, ONLY : allocate_paw_internals
USE grid_paw_variables, ONLY : okpaw
USE rad_paw_routines, ONLY : paw_init
USE bp, ONLY : lberry, lelfield
!
IMPLICIT NONE
@ -38,6 +39,7 @@ SUBROUTINE init_run()
!
CALL allocate_nlpot()
if (okpaw) CALL allocate_paw_internals()
if (okpaw) CALL paw_init()
CALL allocate_locpot()
CALL allocate_wfc()
CALL allocate_bp_efield()

36
PW/init_us_1.f90
View File

@ -41,9 +41,7 @@ subroutine init_us_1
qq_so, dvan_so, okvan
USE uspp_param, ONLY : upf, lmaxq, nbetam, nh, nhm, lmaxkb
USE spin_orb, ONLY : lspinorb, so, rot_ylm, fcoef
!! NEW-AUG !!
USE grid_paw_variables, ONLY : really_do_paw, okpaw, tpawp, aug
!! NEW-AUG !!
!
implicit none
!
@ -243,28 +241,28 @@ subroutine init_us_1
do nb = 1, upf(nt)%nbeta
do mb = nb, upf(nt)%nbeta
ijv = mb * (mb-1) / 2 + nb
paw:& ! in PAW formalism aug. charge is computed elsewhere
paw : & ! in PAW formalism aug. charge is computed elsewhere
if (tpawp(nt)) then
qtot(1:upf(nt)%kkbeta,ijv) = aug(nt)%fun(1:upf(nt)%kkbeta,nb,mb,l)
else
if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. &
(mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then
do ir = 1, upf(nt)%kkbeta
if (rgrid(nt)%r(ir) >=upf(nt)%rinner (l+1) ) then
qtot (ir, ijv) = upf(nt)%qfunc(ir,ijv)
else
ilast = ir
endif
enddo
if ( upf(nt)%rinner (l+1) > 0.0_dp) &
call setqf(upf(nt)%qfcoef (1, l+1, nb, mb), &
qtot(1,ijv), rgrid(nt)%r, upf(nt)%nqf, &
l, ilast)
if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. &
(mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then
do ir = 1, upf(nt)%kkbeta
if (rgrid(nt)%r(ir) >=upf(nt)%rinner (l+1) ) then
qtot (ir, ijv) = upf(nt)%qfunc(ir,ijv)
else
ilast = ir
endif
enddo
if ( upf(nt)%rinner (l+1) > 0.0_dp) &
call setqf(upf(nt)%qfcoef (1, l+1, nb, mb), &
qtot(1,ijv), rgrid(nt)%r, upf(nt)%nqf, &
l, ilast)
endif
endif paw
enddo
enddo
enddo ! mb
enddo ! nb
!
! here we compute the spherical bessel function for each |g|
!

5
PW/init_vloc.f90
View File

@ -15,6 +15,7 @@ subroutine init_vloc()
! potential vloc(ig,it) for each type of atom
!
USE atom, ONLY : msh, rgrid
USE kinds, ONLY : dp
USE uspp_param, ONLY : upf
USE ions_base, ONLY : ntyp => nsp
USE cell_base, ONLY : omega, tpiba2
@ -23,12 +24,12 @@ subroutine init_vloc()
!
implicit none
!
integer :: nt
integer :: nt,k
! counter on atomic types
!
call start_clock ('init_vloc')
vloc(:,:) = 0.d0
vloc(:,:) = 0._dp
do nt = 1, ntyp
!
! compute V_loc(G) for a given type of atom

4
PW/make.depend
View File

@ -602,6 +602,7 @@ init_run.o : ../Modules/grid_paw_variables.o
init_run.o : dynamics_module.o
init_run.o : grid_paw_routines.o
init_run.o : pwcom.o
init_run.o : rad_paw_routines.o
init_us_1.o : ../Modules/atom.o
init_us_1.o : ../Modules/cell_base.o
init_us_1.o : ../Modules/constants.o
@ -704,6 +705,7 @@ mix_rho.o : ../Modules/kind.o
mix_rho.o : ../Modules/uspp.o
mix_rho.o : ../Modules/wavefunctions.o
mix_rho.o : pwcom.o
mix_rho.o : rad_paw_routines.o
mode_group.o : ../Modules/constants.o
mode_group.o : ../Modules/kind.o
move_ions.o : ../Modules/basic_algebra_routines.o
@ -741,6 +743,7 @@ newd.o : ../Modules/wavefunctions.o
newd.o : grid_paw_routines.o
newd.o : noncol.o
newd.o : pwcom.o
newd.o : rad_paw_routines.o
newd.o : realus.o
non_scf.o : ../Modules/control_flags.o
non_scf.o : ../Modules/io_files.o
@ -817,6 +820,7 @@ potinit.o : io_rho_xml.o
potinit.o : noncol.o
potinit.o : pw_restart.o
potinit.o : pwcom.o
potinit.o : rad_paw_routines.o
print_clock_pw.o : ../Modules/control_flags.o
print_clock_pw.o : ../Modules/grid_paw_variables.o
print_clock_pw.o : ../Modules/io_global.o

19
PW/mix_rho.f90
View File

@ -7,7 +7,7 @@
!
#include "f_defs.h"
!
#define ZERO ( 0.D0, 0.D0 )
#define ZERO ( 0._dp, 0._dp )
#define ONLY_SMOOTH_G
!
!----------------------------------------------------------------------------
@ -31,19 +31,18 @@ SUBROUTINE mix_rho( rhocout, rhocin, taukout, taukin, becout, becin, &
USE control_flags, ONLY : imix, ngm0, tr2, io_level
USE io_files, ONLY : find_free_unit
USE cell_base, ONLY : omega
!
!!PAW]
! ... for PAW:
USE uspp_param, ONLY : nhm
USE grid_paw_variables, ONLY : okpaw
!!PAW]
USE rad_paw_routines, ONLY : paw_ddot
!
IMPLICIT NONE
!
! ... First the I/O variable
!
INTEGER :: &
iter, &! (in) counter of the number of iterations
n_iter ! (in) numb. of iterations used in mixing
iter, &! (in) counter of the number of iterations
n_iter ! (in) numb. of iterations used in mixing
COMPLEX(DP) :: &
rhocin(ngm,nspin), rhocout(ngm,nspin)
COMPLEX(DP) :: &
@ -147,7 +146,7 @@ SUBROUTINE mix_rho( rhocout, rhocin, taukout, taukin, becout, becin, &
dr2 = rho_ddot( rhocout, rhocout, ngm, ngm, nspin, ngm )
IF ( tmeta ) dr2 = dr2 + tauk_ddot( taukout, taukout, ngm, ngm, nspin, ngm )
IF ( lda_plus_u ) dr2 = dr2 + ns_ddot( nsout, nsout, nspin )
IF ( okpaw ) dr2 = dr2 + rho1_ddot ( becout, becout ) !PAW
IF ( okpaw ) dr2 = dr2 + PAW_ddot ( becout, becout ) !PAW
!
conv = ( dr2 < tr2 )
!
@ -447,7 +446,7 @@ SUBROUTINE mix_rho( rhocout, rhocin, taukout, taukin, becout, becin, &
!
IF ( okpaw ) &
betamix(i,j) = betamix(i,j) + &
rho1_ddot( df_bec(1,1,1,j), df_bec(1,1,1,i) )
PAW_ddot( df_bec(1,1,1,j), df_bec(1,1,1,i) )
!
betamix(j,i) = betamix(i,j) !symmetrize
!
@ -476,7 +475,7 @@ SUBROUTINE mix_rho( rhocout, rhocin, taukout, taukin, becout, becin, &
work(i) = work(i) + ns_ddot( df_ns(1,1,1,1,i), nsout, nspin )
!
IF ( okpaw ) &
work(i) = work(i) + rho1_ddot( df_bec(1,1,1,i), becout )
work(i) = work(i) + PAW_ddot( df_bec(1,1,1,i), becout )
!
END DO
!
@ -828,6 +827,7 @@ FUNCTION tauk_ddot( tauk1, tauk2, ngm1, ngm2, nspin, gf )
END FUNCTION tauk_ddot
!
!
#ifdef __DONT_DO_THAT_THEN
!----------------------------------------------------------------------------
FUNCTION rho1_ddot( bec1, bec2 )
!----------------------------------------------------------------------------
@ -978,6 +978,7 @@ FUNCTION rho1_ddot( bec1, bec2 )
RETURN
!
END FUNCTION rho1_ddot
#endif
!
!----------------------------------------------------------------------------
FUNCTION ns_ddot( ns1, ns2, nspin )

50
PW/newd.f90
View File

@ -41,7 +41,10 @@ SUBROUTINE newd_g()
!
USE grid_paw_variables, ONLY : really_do_paw, okpaw, tpawp, &
& kdiff, dpaw_ae, dpaw_ps
#ifdef __GRID_PAW
USE grid_paw_routines, ONLY : newd_paw_grid
#endif
USE rad_paw_routines, ONLY : PAW_newd
USE uspp, ONLY : nhtol, nhtolm
!
IMPLICIT NONE
@ -192,6 +195,38 @@ SUBROUTINE newd_g()
!
CALL reduce( nhm * nhm * nat * nspin0, deeq )
!
IF (okpaw) THEN
! prepare non-kinetic paw contribution to D coefficients
! (they are added later in the "atoms" loop)
#ifdef __GRID_PAW
dpaw_ae=0._dp
dpaw_ps=0._dp
CALL newd_paw_grid(na)
IF (na==1) THEN
PRINT '(a)','--------------------------------------------------'
PRINT *, 'GRID AE 1:'
PRINT '(8f15.7)', (((dpaw_ae(jh,ih,na,1)),jh=1,nh(nt)),ih=1,nh(nt))
PRINT '(8f15.7)', 0.,0.,0.,0.,0.,0.,0.,0.
PRINT '(8f15.7)', (((dpaw_ae(jh,ih,na,2)),jh=1,nh(nt)),ih=1,nh(nt))
ENDIF
dpaw_ae=0._dp
dpaw_ps=0._dp
#endif
CALL PAW_newd(dpaw_ae, dpaw_ps)
!#define __VERBOSE_PAW_NEWD
#ifdef __VERBOSE_PAW_NEWD
PRINT *, 'deeq 1:'
PRINT '(8f15.7)', ((deeq(jh,ih,na,1),jh=1,nh(nt)),ih=1,nh(nt))
PRINT *, 'ddd AE 1:'
PRINT '(8f15.7)', ((dpaw_ae(jh,ih,na,1),jh=1,nh(nt)),ih=1,nh(nt))
PRINT *, 'ddd PS 1:'
PRINT '(8f15.7)', ((dpaw_ps(jh,ih,na,1),jh=1,nh(nt)),ih=1,nh(nt))
! STOP
#endif
ENDIF
atoms : &
DO na = 1, nat
!
nt = ityp(na)
@ -209,9 +244,8 @@ SUBROUTINE newd_g()
END IF
!
ELSE if_noncolin
if_paw:&
paw:&
IF (okpaw) THEN
CALL newd_paw_grid(na)
!
DO is = 1, nspin
!
@ -221,10 +255,10 @@ SUBROUTINE newd_g()
DO jh = ih, nh(nt)
mb = indv(jh,nt)
!
IF ( nhtol (ih, nt) == nhtol (jh, nt) .AND. &
nhtolm(ih, nt) == nhtolm(jh, nt) ) THEN
IF ( nhtol (ih, nt) == nhtol (jh, nt) .and. &
nhtolm(ih, nt) == nhtolm(jh, nt) ) THEN
deeq(ih,jh,na,is) = deeq(ih,jh,na,is) + &
kdiff(nb,mb,nt)
kdiff(nb,mb,nt)
END if
!
deeq(ih,jh,na,is) = deeq(ih,jh,na,is) + &
@ -237,7 +271,7 @@ SUBROUTINE newd_g()
END DO
!
END DO
ELSE if_paw
ELSE paw
!
DO is = 1, nspin
!
@ -253,11 +287,11 @@ SUBROUTINE newd_g()
END DO
!
END DO
END IF if_paw
END IF paw
!
END IF if_noncolin
!
END DO
END DO atoms
!
DEALLOCATE( aux, qgm_na, qgm, qmod, ylmk0 )
!

33
PW/paw_xc.f90
View File

@ -6,6 +6,12 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
#if ! defined __GRID_PAW
! empty subroutine because some compilers refuses to compile an empty file
SUBROUTINE unuseful_subroutine
RETURN
END SUBROUTINE unuseful_subroutine
#else
!----------------------------------------------------------------------------
SUBROUTINE paw_v_xc( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, nl, ngm, g, nspin, alat, omega, etxc, vtxc, v )
@ -26,6 +32,7 @@ SUBROUTINE paw_v_xc( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nspin, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, ngm, nl(ngm)
@ -70,11 +77,17 @@ SUBROUTINE paw_v_xc( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
vanishing_mag = 1.D-20
!
!
#ifdef __NO_XC
write(*,*) "##### skipping XC in PAW_xc"
RETURN
#endif
etxc = 0.D0
vtxc = 0.D0
v(:,:) = 0.D0
rhoneg = 0.D0
!
!goto 666
IF ( nspin == 1 .OR. ( nspin == 4 .AND. .NOT. domag ) ) THEN
!
! ... spin-unpolarized case
@ -219,8 +232,14 @@ SUBROUTINE paw_v_xc( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
! ... add gradient corrections (if any)
!
!
#ifdef __ONLY_GCXC
v(:,:) = 0._dp
#endif
#ifdef __NO_GCXC
#else
CALL paw_gradcorr( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, nl, ngm, g, alat, omega, nspin, etxc, vtxc, v )
#endif
!
CALL reduce( 1, vtxc )
CALL reduce( 1, etxc )
@ -391,6 +410,7 @@ SUBROUTINE paw_gradcorr( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
IF ( rh > epsr ) THEN
!
IF ( get_igcc() == 3 ) THEN
write(*,*) "USING MORE SPIN"
!
rup = rhoout(k,1)
rdw = rhoout(k,2)
@ -489,7 +509,17 @@ SUBROUTINE paw_gradcorr( rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
CALL paw_grad_dot( nrx1, nrx2, nrx3, nr1, nr2, nr3, &
nrxx, h(1,1,is), ngm, g, nl, alat, dh )
!
v(:,is) = v(:,is) - dh(:)
#ifdef __PAW_ONLYHAM
write(0,*) "######-- Only using hamiltonian contribution to gcxc"
v(:,is) = - dh(:)
#else
#ifdef __PAW_NONHAM
write(0,*) "######-- Only including NON HAM in gcxc"
#else
v(:,is) = v(:,is) - dh(:)
#endif
#endif
!
vtxcgc = vtxcgc - SUM( dh(:) * rhoout(:,is) )
!
@ -658,3 +688,4 @@ SUBROUTINE paw_grad_dot( nrx1, nrx2, nrx3, nr1, nr2, nr3, &
!
END SUBROUTINE paw_grad_dot
#endif

19
PW/potinit.f90
View File

@ -50,8 +50,14 @@ SUBROUTINE potinit()
USE io_rho_xml, ONLY : read_rho
!
USE uspp, ONLY : becsum
USE grid_paw_variables, ONLY : rho1, rho1t
#ifdef __GRID_PAW
USE grid_paw_variables, ONLY : rho1, rho1t, okpaw
USE grid_paw_routines, ONLY : atomic_becsum, compute_onecenter_charges, compute_onecenter_potentials
#else
USE grid_paw_variables, ONLY : okpaw
USE grid_paw_routines, ONLY : atomic_becsum
#endif
USE rad_paw_routines, ONLY : PAW_potential
!
IMPLICIT NONE
!
@ -239,9 +245,14 @@ SUBROUTINE potinit()
! ... PAW initialization: from atomic augmentation channel occupations
! ... compute corresponding one-center charges and potentials
!
CALL atomic_becsum()
CALL compute_onecenter_charges(becsum,rho1,rho1t)
CALL compute_onecenter_potentials(becsum,rho1,rho1t)
IF ( okpaw ) THEN
CALL atomic_becsum()
CALL PAW_potential(becsum)
#ifdef __GRID_PAW
CALL compute_onecenter_charges(becsum,rho1,rho1t)
CALL compute_onecenter_potentials(becsum,rho1,rho1t)
#endif
ENDIF
!
RETURN
!

52
PW/print_clock_pw.f90
View File

@ -178,25 +178,39 @@ SUBROUTINE print_clock_pw()
CALL print_clock ('cycleig')
#endif
!
IF ( okpaw ) THEN
WRITE( stdout, '(5X,"PAW routines")' )
call print_clock ('init_prad')
call print_clock ('paw_prod_p')
call print_clock ('one-charge')
call print_clock ('one-pot')
call print_clock ('pvan2')
call print_clock ('set_paw_rhoc')
call print_clock ('v_h_grid')
call print_clock ('newd_paw')
call print_clock ('init_pawvloc')
call print_clock ('vloc_of_g_no')
call print_clock ('paw_setlocal')
! radial routines:
CALL print_clock ('PAW_energy')
CALL print_clock ('PAW_rho_lm')
CALL print_clock ('PAW_h_energy')
CALL print_clock ('PAW_sph_int')
END IF
IF ( okpaw ) THEN
WRITE( stdout, '(5X,"PAW routines")' )
#ifdef __GRID_PAW
call print_clock ('init_prad')
call print_clock ('paw_prod_p')
call print_clock ('one-charge')
call print_clock ('one-pot')
call print_clock ('pvan2')
call print_clock ('set_paw_rhoc')
call print_clock ('v_h_grid')
call print_clock ('newd_paw')
call print_clock ('init_pawvloc')
call print_clock ('vloc_of_g_no')
call print_clock ('paw_setlocal')
#endif
! radial routines:
CALL print_clock ('PAW_pot')
CALL print_clock ('PAW_newd')
CALL print_clock ('PAW_int')
CALL print_clock ('PAW_ddot')
CALL print_clock ('PAW_rad_init')
CALL print_clock ('PAW_energy')
! second level routines:
CALL print_clock ('PAW_rho_lm')
CALL print_clock ('PAW_h_pot')
CALL print_clock ('PAW_xc_pot')
CALL print_clock ('PAW_lm2rad')
CALL print_clock ('PAW_rad2lm')
! third level, or deeper:
CALL print_clock ('PAW_gcxc_v')
CALL print_clock ('PAW_div')
CALL print_clock ('PAW_grad')
END IF
!
RETURN
!

1790
PW/rad_paw_routines.f90
View File

File diff suppressed because it is too large Load Diff

393
PW/rad_paw_trash.f90
View File

@ -18,6 +18,399 @@
! WRITE(6,"(a,2f15.7)") "== LM=8:",ecomps(8,1,1,1),ecomps(8,1,1,2)
! WRITE(6,"(a,2f15.7)") "== LM=9:",ecomps(9,1,1,1),ecomps(9,1,1,2)
! WRITE(6,"(a,2f15.7)") "=============================================="
!___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
!!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!!
!!! This is the main driver of PAW routines, it uses directly or indirectly
!!! all the other routines of the module
!!
SUBROUTINE PAW_energy(becsum,e_atom)
USE kinds, ONLY : DP
USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp
USE atom, ONLY : g => rgrid
USE grid_paw_variables, ONLY : pfunc, ptfunc, tpawp, &
aerho_atc, psrho_atc, aug
USE uspp_param, ONLY : nhm, lmaxq
REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupations
REAL(DP), INTENT(OUT),OPTIONAL :: &
e_atom(nat,2,2) ! {# of atoms}, {XC|H}, {AE|PS}
INTEGER, PARAMETER :: AE = 1, PS = 2,& ! All-Electron and Pseudo
XC = 1, H = 2 ! XC and Hartree
INTEGER :: i_what ! counter on AE and PS
INTEGER :: na,first_nat,last_nat ! atoms counters and indexes
! hartree energy scalar fields expanded on Y_lm
REAL(DP), ALLOCATABLE :: rho_lm(:,:,:) ! radial density expanded on Y_lm
! REAL(DP), ALLOCATABLE :: v_h_lm(:,:) ! hartree potential, summed on spins
!
REAL(DP) :: e_h(lmaxq**2,nspin)! hartree energy components
REAL(DP) :: e,e1,e2 ! placeholders
! xc variables:
REAL(DP) :: e_xc(lmaxq**2,nspin)! UNUSED! XC energy components
REAL(DP), POINTER :: rho_core(:,:) ! pointer to AE/PS core charge density
TYPE(paw_info) :: i ! minimal info on atoms
! BEWARE THAT HARTREE ONLY DEPENDS ON THE TOTAL RHO NOT ON RHOUP AND RHODW SEPARATELY...
! TREATMENT OF NSPIN>1 MUST BE CHECKED AND CORRECTED
CALL start_clock ('PAW_energy')
! nullify energy components:
IF( present(e_atom) ) e_atom(:,:,:) = 0._dp
! The following operations will be done, first on AE, then on PS part
! furthermore they will be done for one atom at a time (to reduce memory usage)
! in the future code will be parallelized on atoms:
!
! 1. build rho_lm (PAW_rho_lm)
! 2. compute v_h_lm and use it to compute HARTREE
! energy (PAW_h_energy)
! 3. compute XC energy
! a. build rho_rad(theta, phi) from rho_lm
! + if using gradient correction compute also grad(rho)
! b. compute XC energy from rho_rad
! c. iterate on enough directions to compute the
! spherical integral
! CHECK: maybe we don't need to alloc/dealloc rho_lm every time
first_nat = 1
last_nat = nat
! Operations from here on are (will be) parallelized on atoms, for the
! moment this OpenMP directive gives single-core parallelization:
!$OMP PARALLEL DO default(shared) reduction(+:e_atom) &
!$OMP private(rho_lm, e_xc, e_h, e,e1,e2, rho_core, na, i_what, i)
atoms: DO na = first_nat, last_nat
!
i%a = na ! the index of the atom
i%t = ityp(na) ! the type of atom na
i%m = g(i%t)%mesh! radial mesh size for atom na
!
ifpaw: IF (tpawp(i%t)) THEN
!
! Arrays are allocated inside the cycle to allow reduced
! memory usage as differnt atoms have different meshes
ALLOCATE(rho_lm(i%m,lmaxq**2,nspin))
!
whattodo: DO i_what = AE, PS
i%w = i_what ! spherical_gradient likes to know
! STEP: 1 [ build rho_lm (PAW_rho_lm) ]
NULLIFY(rho_core)
IF (i_what == AE) THEN
! to pass the atom indes is dirtyer but faster and
! uses less memory than to pass only a hyperslice of the array
CALL PAW_rho_lm(i, becsum, pfunc, rho_lm)
! used later for xc energy:
rho_core => aerho_atc
ELSE
CALL PAW_rho_lm(i, becsum, ptfunc, rho_lm, aug)
! optional argument for pseudo part --> ^^^
! used later for xc energy:
rho_core => psrho_atc
ENDIF
! STEP: 2 [ compute Hartree energy ]
! 2a. use rho_lm to compute hartree potential (PAW_v_h)
e = PAW_h_energy(i, rho_lm, e_h)
IF( present(e_atom) ) e_atom(i%a, H, i_what) = e
! STEP: 3 [ compute XC energy ]
e1 = PAW_xc_energy(i, rho_lm, rho_core, e_xc,0)
IF( present(e_atom) ) e_atom(i%a, XC, i_what) = e1
e2 = PAW_xc_energy(i, rho_lm, rho_core, e_xc,1)
write(*,"(a,2i2,2f10.5,a)") " == ", i_what, na, e, e1," (2)"
ENDDO whattodo
! cleanup
DEALLOCATE(rho_lm)
ENDIF ifpaw
ENDDO atoms
!$OMP END PARALLEL DO
CALL stop_clock ('PAW_energy')
END SUBROUTINE PAW_energy
!___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
!!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!!
!!! add gradient correction to energy, code mostly adapted from ../atomic/vxcgc.f90
SUBROUTINE PAW_gcxc_energy(i, rho,core,grho,e)
USE kinds, ONLY : DP
USE ions_base, ONLY : ityp
USE lsda_mod, ONLY : nspin
USE atom, ONLY : g => rgrid
USE parameters, ONLY : ntypx
USE constants, ONLY : fpi,pi,e2
USE funct, ONLY : gcxc, gcx_spin, gcc_spin
!
TYPE(paw_info) :: i ! atom's minimal info
REAL(DP), INTENT(IN) :: rho(i%m,nspin) ! radial density,
REAL(DP), INTENT(IN) :: grho(i%m,nspin)! square gradient of rho
REAL(DP), INTENT(IN) :: core(i%m,ntypx) ! spherical core density
REAL(DP), INTENT(INOUT):: e(i%m) ! radial local xc energy
!
INTEGER :: k ! counter for mesh
! workspaces:
REAL(DP) :: arho, sgn
REAL(DP) :: sx,sc ! x and c energy from gcxc
REAL(DP) :: v1x,v2x,v1c,v2c ! potentials from gcxc (unused)
REAL(DP),PARAMETER :: eps = 1.D-12 ! 1.e-12 may be small enough
! for nspin>1
REAL(DP) :: rh,grh2,zeta
REAL(DP) :: v1cup, v1cdw,v1xup, v1xdw, v2xup, v2xdw !(unused)
OPTIONAL_CALL start_clock ('PAW_gcxc_e')
IF (nspin.eq.1) THEN
DO k=1,i%m
arho = rho(k,1)*g(i%t)%rm2(k) + core(k,i%t)
sgn = sign(1._dp,arho)
arho = abs(arho)
IF (arho.gt.eps.and.abs(grho(k,1)).gt.eps) THEN
CALL gcxc(arho,grho(k,1),sx,sc,v1x,v2x,v1c,v2c)
e(k) = e(k) + sgn *e2 *(sx+sc)*g(i%t)%r2(k) !&
ENDIF
ENDDO
ELSE
! this is the \sigma-GGA case
DO k=1,i%m
CALL gcx_spin (rho(k,1), rho(k,2), grho(k,1), grho(k,2), &
sx, v1xup, v1xdw, v2xup, v2xdw)
rh = rho(k,1) + rho(k,2)
IF (rh.gt.eps) THEN
zeta = (rho (k,1) - rho (k,2) ) / rh
! FIXME: next line is wrong because it looses sign!!
grh2 = (sqrt(grho(k,1)) + sqrt(grho(k,2)) ) **2
CALL gcc_spin (rh, zeta, grh2, sc, v1cup, v1cdw, v2c)
ELSE
sc = 0._dp
! v1cup = 0.0_dp
! v1cdw = 0.0_dp
! v2c = 0.0_dp
ENDIF
e(k) = e(k) + e2*(sx+sc)
! vgc(i,1)= v1xup+v1cup
! vgc(i,2)= v1xdw+v1cdw
! h(i,1) =((v2xup+v2c)*grho(i,1)+v2c*grho(i,2))*r2(i)
! h(i,2) =((v2xdw+v2c)*grho(i,2)+v2c*grho(i,1))*r2(i)
ENDDO
ENDIF
OPTIONAL_CALL stop_clock ('PAW_gcxc_e')
END SUBROUTINE PAW_gcxc_energy
!___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
!!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!!
!!! use the density produced by sum_rad_rho to compute hartree potential
!!! the potential is then directly integrated to compute hartree energy
FUNCTION PAW_h_energy(i, rho_lm, e_lm)
USE kinds, ONLY : DP
USE constants, ONLY : fpi
USE parameters, ONLY : ntypx
USE lsda_mod, ONLY : nspin
USE uspp_param, ONLY : nh, lmaxq
USE atom, ONLY : g => rgrid
REAL(DP) :: PAW_h_energy ! total hartree energy
!
TYPE(paw_info) :: i ! the number of the atom
REAL(DP), INTENT(IN) :: rho_lm(i%m,lmaxq**2,nspin) ! charge density as lm components
REAL(DP), OPTIONAL,INTENT(OUT) :: e_lm(lmaxq**2) ! out: energy components
!
REAL(DP) :: aux(i%m) ! workspace
REAL(DP) :: par_energy ! workspace
REAL(DP) :: v_lm(i%m,lmaxq**2) ! potential as lm components
REAL(DP) :: trho_lm(i%m,lmaxq**2) ! total spin-summed charge density
INTEGER :: ispin, & ! counter on spins
lm,l ! counter on composite angmom lm = l**2 +m
OPTIONAL_CALL start_clock ('PAW_h_energy')
! init total energy and its lm components
PAW_h_energy = 0._dp
IF (present(e_lm)) e_lm(:) = 0._dp
! sum up spin components, respecting column-majorness
trho_lm(:,:) = rho_lm(:,:,1)
DO ispin = 2,nspin ! if nspin < 2 it jumps to ENDDO
trho_lm(:,:) = trho_lm(:,:) + rho_lm(:,:,ispin)
ENDDO
! compute the hartree potential
CALL PAW_h_potential(i, rho_lm, v_lm)
DO lm = 1, lmaxq**2
!
! now energy is computed integrating v_h^{lm} * \sum_{spin} rho^{lm}
! and summing on lm
aux(:) = trho_lm(:,lm) * v_lm(:,lm)
CALL simpson (i%m, aux, g(i%t)%rab, par_energy)
!
PAW_h_energy = PAW_h_energy + par_energy
IF (present(e_lm)) e_lm(lm) = par_energy
ENDDO ! lm
OPTIONAL_CALL stop_clock ('PAW_h_energy')
END FUNCTION PAW_h_energy
!___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___ ___
!!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!! !!!!
!!! use the density produced by sum_rad_rho to compute xc potential and energy, as
!!! xc functional is not diagonal on angular momentum numerical integration is performed
FUNCTION PAW_xc_energy(i, rho_lm, rho_core, e_lm, task)
USE kinds, ONLY : DP
USE constants, ONLY : fpi
USE parameters, ONLY : ntypx
USE lsda_mod, ONLY : nspin
USE uspp_param, ONLY : lmaxq
USE ions_base, ONLY : ityp
USE radial_grids, ONLY : ndmx
USE atom, ONLY : g => rgrid
REAL(DP) :: PAW_xc_energy ! total xc energy
!
TYPE(paw_info) :: i ! atom's minimal info
INTEGER, INTENT(IN) :: task!remove me
REAL(DP), INTENT(IN) :: rho_lm(i%m,lmaxq**2,nspin)! charge density as lm components
REAL(DP), INTENT(IN) :: rho_core(ndmx,ntypx) ! core charge, radial and spherical
! TODO:
REAL(DP), OPTIONAL,INTENT(OUT) :: e_lm(lmaxq**2) ! out: energy components
!
INTEGER :: lsd ! switch to control local spin density
!
REAL(DP) :: rho_loc(2) = (/0._dp, 0._dp/)
! local density (workspace), up and down
REAL(DP) :: e ! workspace
REAL(DP) :: e_rad(i%m) ! radial energy (to be integrated)
REAL(DP) :: rho_rad(i%m,nspin) ! workspace (radial slice of rho)
INTEGER :: ix,k ! counters on directions and radial grid
! for gradient correction only:
REAL(DP),ALLOCATABLE :: grho_rad(:,:)! workspace (radial slice of grad(rho))
! functions from atomic code:
REAL(DP),EXTERNAL :: exc_t
! to prevent warnings (this quantity is not implemented,
! and maybe it never will):
e_lm = 0._dp
OPTIONAL_CALL start_clock ('PAW_xc_nrg')
lsd = nspin-1
! init for gradient correction
IF (do_gcxc) ALLOCATE(grho_rad(i%m,nspin))
PAW_xc_energy = 0._dp
DO ix = 1, nx
! LDA (and LSDA) part (no gradient correction):
CALL PAW_lm2rad(i, ix, rho_lm, rho_rad)
!
DO k = 1,i%m
rho_loc(1:nspin) = rho_rad(k,1:nspin)*g(i%t)%rm2(k)
!
e_rad(k) = exc_t(rho_loc, rho_core(k,i%t), lsd)&
* (SUM(rho_rad(k,1:nspin))+rho_core(k,i%t)*g(i%t)%r2(k))
ENDDO
gradient_correction:& ! add it!
IF (do_gcxc) THEN
IF (task == 1) e_rad(:) = 0._dp ! reset the energy, so that only the correction is displayed
CALL PAW_gradient(i, ix, rho_lm, rho_rad, rho_core, grho_rad)
! v-------------^
CALL PAW_gcxc_energy(i, rho_rad,rho_core,grho_rad, e_rad)
ENDIF gradient_correction
!
! integrate radial slice of xc energy:
CALL simpson (i%m, e_rad, g(i%t)%rab, e)
! integrate on sph. surface v-----^
PAW_xc_energy = PAW_xc_energy + e * ww(ix)
ENDDO
! cleanup for gc
IF (do_gcxc) DEALLOCATE(grho_rad)
OPTIONAL_CALL stop_clock ('PAW_xc_nrg')
END FUNCTION PAW_xc_energy
SUBROUTINE PAW_brute_radial_ddd(becsum)
USE kinds, ONLY : DP
USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp
USE uspp_param, ONLY : nhm, nh
USE grid_paw_variables, ONLY : tpawp
REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupations
REAL(DP) :: ddd(nhm,nhm,nat,nspin,1:2)! cross band occupations
INTEGER, PARAMETER :: AE = 1, PS = 2, AEmPS=0 ! All-Electron and Pseudo
INTEGER :: ih, jh, ijh
INTEGER :: na, nt, ispin
!
REAL(DP), PARAMETER :: eps = 1.D-6
!
REAL(DP) :: b1,b2,delta
REAL(DP) :: bectmp(nhm*(nhm+1)/2,nat,nspin)
REAL(DP) :: e(nat,2,2), ee ! {# of atoms}, {XC|H}, {AE|PS}
CALL start_clock('PAW_brddd')
!
ddd = 0._dp
ispin = 1 ! FIXME
bectmp = becsum
DO na = 1,nat
nt = ityp(na)
IF ( tpawp(nt) ) THEN
!
DO ih = 1, nh(nt)
DO jh = ih, nh(nt)
! write(6,"(a,3i4)") "now doing:",ih,jh,ijh
!
ijh = jh * (jh-1) / 2 + ih
!
b1 = becsum(ijh,na,ispin) + eps
b2 = becsum(ijh,na,ispin) - eps
delta = (b1-b2)
!
bectmp(ijh,na,ispin) = b1
CALL PAW_potential(bectmp, ee, e)
ddd(ih,jh,na,ispin,AE) = - SUM(e(na,:,AE))
ddd(ih,jh,na,ispin,PS) = - SUM(e(na,:,PS))
!
bectmp(ijh,na,ispin) = b2
CALL PAW_potential(bectmp,ee, e)
ddd(ih,jh,na,ispin,AE) = ddd(ih,jh,na,ispin,AE) &
+ SUM(e(na,:,AE))
ddd(ih,jh,na,ispin,PS) = ddd(ih,jh,na,ispin,PS) &
+ SUM(e(na,:,PS))
!
ddd(ih,jh,na,ispin,:) = ddd(ih,jh,na,ispin,:) / delta
! write(6,*) ddd(ih,jh,na,ispin,:)
ddd(jh,ih,na,ispin,AE) = ddd(ih,jh,na,ispin,AE)
ddd(jh,ih,na,ispin,PS) = ddd(ih,jh,na,ispin,PS)
!
bectmp(ijh,na,ispin) = becsum(ijh,na,ispin)
!
END DO ! jh
END DO ! ih
!
ENDIF
ENDDO
!
CALL stop_clock('PAW_brddd')
PRINT *, "RADIAL VERSION"
PRINT *, 'D - D1'
PRINT '(8f15.7)', ((ddd(jh,ih,1,1,AE),jh=1,nh(nt)),ih=1,nh(nt))
PRINT *, 'D - D1~'
PRINT '(8f15.7)', ((ddd(jh,ih,1,1,PS),jh=1,nh(nt)),ih=1,nh(nt))
END SUBROUTINE PAW_brute_radial_ddd
SUBROUTINE integrate_pfunc

2
PW/read_pseudo.f90
View File

@ -57,7 +57,7 @@ subroutine readpp
END IF
ALLOCATE ( upf(ntyp) )
do nt = 1, ntyp
tpawp(nt) = .false.
!tpawp(nt) = .false.
!
! obsolescent variables, not read from UPF format, no longer used
!

3
PWCOND/Makefile
View File

@ -54,6 +54,7 @@ MODULES = \
../Modules/control_flags.o \
../Modules/descriptors.o \
../Modules/electrons_base.o \
../Modules/exc_t.o \
../Modules/fft_base.o \
../Modules/fft_scalar.o \
../Modules/fft_types.o \
@ -88,6 +89,8 @@ MODULES = \
../Modules/upf_to_internal.o \
../Modules/uspp.o \
../Modules/version.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../Modules/wavefunctions.o \
../Modules/xml_io_base.o

6
atomic/Makefile
View File

@ -37,7 +37,6 @@ elsd_highv.o \
elsdps.o \
elsdps_paw.o \
esic.o \
exc_t.o \
find_qi.o \
gener_pseudo.o \
green.o \
@ -78,8 +77,6 @@ trou.o \
vdpack.o \
vext.o \
vpack.o \
vxcgc.o \
vxc_t.o \
v_of_rho_at.o \
write_cpmd.o \
write_paw_recon.o \
@ -98,6 +95,7 @@ MODULES = \
../Modules/cell_base.o \
../Modules/constants.o \
../Modules/control_flags.o \
../Modules/exc_t.o \
../Modules/functionals.o \
../Modules/grid_paw_variables.o \
../Modules/io_global.o \
@ -116,6 +114,8 @@ MODULES = \
../Modules/mp_global.o \
../Modules/mp.o \
../Modules/uspp.o \
../Modules/vxc_t.o \
../Modules/vxcgc.o \
../PW/startup.o
TLDEPS= bindir mods libs

24
atomic/atomic_paw.f90
View File

@ -201,7 +201,7 @@ CONTAINS
REAL(dp) :: twosigma2, rm, gaussian(ndmx) ! needed for "gaussian" augfun
!
mesh = grid%mesh
irc = maxval(ikk(1:nbeta))
irc = maxval(ikk(1:nbeta))+8
CALL nullify_pseudo_paw(pawset_)
CALL allocate_pseudo_paw(pawset_,ndmx,nwfsx,lmaxx)
pawset_%symbol = atom_name(nint(zed))
@ -360,7 +360,15 @@ CONTAINS
pawset_%augfun(1:ircm,ns,ns1,l3) = &
pawset_%augmom(ns,ns1,l3) * grid%r2(1:ircm) * &
(xc(1) * j1(1:ircm,1) + xc(2) * j1(1:ircm,2))
pawset_%augfun(1:mesh,ns1,ns,l3)=pawset_%augfun(1:mesh,ns,ns1,l3)
! change the sign of augfun if they are negative, and symmetrize them:
! IF ( SUM(pawset_%augfun(1:mesh,ns,ns1,l3)) > 0._dp ) THEN
pawset_%augfun(1:mesh,ns1,ns,l3)=pawset_%augfun(1:mesh,ns,ns1,l3)
! ELSE
! CALL infomsg('us2paw', 'WARNING: augmentation functions are negative, changing their sign!')
! pawset_%augfun(1:mesh,ns1,ns,l3)=-pawset_%augfun(1:mesh,ns,ns1,l3)
! pawset_%augfun(1:mesh,ns,ns1,l3)= pawset_%augfun(1:mesh,ns1,ns,l3)
! ENDIF
!
end do
DEALLOCATE (j1)
@ -769,6 +777,8 @@ CONTAINS
INTEGER :: is, ns, ns1, l
REAL(dp) :: aux(ndmx), dd
REAL(DP), EXTERNAL :: int_0_inf_dr
REAL(dp):: dddd(nwfsx,nwfsx,3)
!
! D^ = Int Q*v~
! D1-D1~ = kindiff + Int[ae*v1*ae - ps*v1~*ps - Q*v1~]
@ -782,6 +792,7 @@ CONTAINS
pawset_%augfun(1:pawset_%grid%mesh,ns,ns1,0) * &
veffps_(1:pawset_%grid%mesh,is)
dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
dddd(ns,ns1,1) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
! Int[ae*v1*ae]
aux(1:pawset_%grid%mesh) = &
pawset_%aewfc(1:pawset_%grid%mesh,ns ) * &
@ -789,6 +800,7 @@ CONTAINS
veff1_(1:pawset_%grid%mesh,is)
dd = dd + &
int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
dddd(ns,ns1,2) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
! Int[ps*v1~*ps + aufun*v1~]
aux(1:pawset_%grid%mesh) = &
( pawset_%pswfc(1:pawset_%grid%mesh,ns ) * &
@ -797,6 +809,7 @@ CONTAINS
veff1ps_(1:pawset_%grid%mesh,is)
dd = dd - &
int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
dddd(ns,ns1,3) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
! collect
ddd_(ns,ns1,is) = pawset_%kdiff(ns,ns1) + dd
ddd_(ns1,ns,is) = ddd_(ns,ns1,is)
@ -804,6 +817,13 @@ CONTAINS
END DO
END DO
END DO
write(*,*) 'deeq'
write(*,'(4f13.8)') dddd(1:4,1:4,1)
write(*,*) 'ddd ae'
write(*,'(4f13.8)') dddd(1:4,1:4,2)
write(*,*) 'ddd ps'
write(*,'(4f13.8)') dddd(1:4,1:4,3)
END SUBROUTINE compute_nonlocal_coeff
!
!============================================================================

5
atomic/gener_pseudo.f90
View File

@ -212,13 +212,10 @@ subroutine gener_pseudo
! compute the phi functions
!
if (lpaw.and.lnc2paw) then
! first compute possibly harder NC pseudowfcts to be
! first compute possibly harder NC pseudowfcs to be
! used as AE reference for PAW generation
nnode=0
!call compute_phi(lam,ik, psi_in,phis(1,ns),xc,1,occ,enls(ns),els(ns)) CURRENT
call compute_phi(lam,iknc2paw,psi_in,phis(1,ns),xc,1,occ,enls(ns),els(ns))
!call compute_phi(lam,iknc2paw,ik,nwf0,ns,xc,1,nnode,ocs(ns)) PAW
!call compute_phi(lam,ik, ik,nwf0,ns,xc,1,nnode,ocs(ns)) OLD
psipaw(1:grid%mesh,ns)=phis(1:grid%mesh,ns)
endif
!

13
atomic/make.depend
View File

@ -16,6 +16,9 @@ ascheqps_drv.o : ../Modules/kind.o
ascheqps_drv.o : ../Modules/radial_grids.o
ascheqps_drv.o : ld1inc.o
ascheqps_drv.o : parameters.o
ascheqps_old.o : ../Modules/io_global.o
ascheqps_old.o : ../Modules/kind.o
ascheqps_old.o : ../Modules/radial_grids.o
atomic_paw.o : ../Modules/constants.o
atomic_paw.o : ../Modules/functionals.o
atomic_paw.o : ../Modules/kind.o
@ -125,8 +128,6 @@ elsdps_paw.o : parameters.o
esic.o : ../Modules/kind.o
esic.o : ../Modules/radial_grids.o
esic.o : ld1inc.o
exc_t.o : ../Modules/functionals.o
exc_t.o : ../Modules/kind.o
find_qi.o : ../Modules/kind.o
find_qi.o : ld1inc.o
gener_pseudo.o : ../Modules/io_global.o
@ -279,14 +280,6 @@ v_of_rho_at.o : ld1inc.o
vdpack.o : ../Modules/kind.o
vext.o : ../Modules/kind.o
vpack.o : ../Modules/kind.o
vxc_t.o : ../Modules/functionals.o
vxc_t.o : ../Modules/io_global.o
vxc_t.o : ../Modules/kind.o
vxcgc.o : ../Modules/constants.o
vxcgc.o : ../Modules/functionals.o
vxcgc.o : ../Modules/kind.o
vxcgc.o : ../Modules/radial_grids.o
vxcgc.o : ld1inc.o
write_cpmd.o : ../Modules/constants.o
write_cpmd.o : ../Modules/functionals.o
write_cpmd.o : ../Modules/kind.o