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

Several Gamma-specific routines merged into PW/ 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@356 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2003-10-29 18:53:40 +00:00
parent ac3144b3d4
commit e34c8c50e3
23 changed files with 184 additions and 2157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
D3/d3toten.f90
View File

@ -29,6 +29,7 @@ program d3toten
call init_clocks (.true.)
call start_clock ('D3TOTEN')
gamma_only = .false.
call startup (nd_nmbr, code, version_number)
!
! Initialization routines

30
Gamma/Makefile
View File

@ -188,18 +188,15 @@ PWOBJS=../PW/pwcom.o \
../PW/gather.o \
../PW/scatter.o \
../PW/init_pool.o \
../PW/bp_bess.o \
../PW/bp_calc_btq.o \
../PW/bp_c_phase.o \
../PW/bp_dbess.o \
../PW/bp_qvan3.o \
../PW/bp_radin.o \
../PW/bp_strings.o \
../PW/bp_ylm_q.o \
../PW/bp_zgedi.o \
../PW/bp_zgefa.o \
../PW/ccalbec.o \
../PW/becmod.o
../PW/addusdens.o \
../PW/allocate_fft.o \
../PW/atomic_rho.o \
../PW/gradcorr.o \
../PW/interpolate.o \
../PW/mix_rho.o \
../PW/set_rhoc.o \
../PW/setlocal.o \
../PW/v_of_rho.o
MODULES = ../Modules/*.o ../PW/error_handler.o \
@ -208,26 +205,17 @@ MODULES = ../Modules/*.o ../PW/error_handler.o \
GAMMA=rbecmod.o \
fake.o \
add_vuspsi.o \
addusdens.o \
allocate_fft.o \
allocate_wfc.o \
atomic_rho.o \
c_bands.o \
rdiaghg.o \
regterg.o \
force_us.o \
gradcorr.o \
h_psi.o \
interpolate.o \
mix_rho.o \
pw_gemm.o \
rotate_wfc.o \
s_psi.o \
set_rhoc.o \
setlocal.o \
stres_us.o \
sum_band.o \
v_of_rho.o \
vloc_psi.o \
wfcinit.o

105
Gamma/addusdens.f90
View File

@ -1,105 +0,0 @@
!
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
subroutine addusdens
!----------------------------------------------------------------------
!
! This routine adds to the charge density the part which is due to
! the US augmentation.
!
#include "machine.h"
use pwcom
USE wavefunctions, ONLY : psic
implicit none
!
! here the local variables
!
integer :: ig, na, nt, ih, jh, ijh, is
! counters
real(kind=DP), allocatable :: qmod (:), ylmk0 (:,:)
! the modulus of G
! the spherical harmonics
complex(kind=DP) :: skk
complex(kind=DP), allocatable :: aux (:,:)
! work space for FFT
! work space for rho(G,nspin)
if (.not.okvan) return
call start_clock ('addusdens')
allocate (aux ( ngm, nspin))
allocate (qmod( ngm))
allocate (ylmk0( ngm, lqx * lqx))
aux (:,:) = (0.d0, 0.d0)
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
do nt = 1, ntyp
if (tvanp (nt) ) then
ijh = 0
do ih = 1, nh (nt)
do jh = ih, nh (nt)
#ifdef DEBUG_ADDUSDENS
call start_clock ('addus:qvan2')
#endif
call qvan2 (ngm, ih, jh, nt, qmod, qgm, ylmk0)
#ifdef DEBUG_ADDUSDENS
call stop_clock ('addus:qvan2')
#endif
ijh = ijh + 1
do na = 1, nat
if (ityp (na) .eq.nt) then
!
! Multiply becsum and qg with the correct structure factor
!
#ifdef DEBUG_ADDUSDENS
call start_clock ('addus:aux')
#endif
do is = 1, nspin
do ig = 1, ngm
skk = eigts1 (ig1 (ig), na) * &
eigts2 (ig2 (ig), na) * &
eigts3 (ig3 (ig), na)
aux(ig,is)=aux(ig,is) + qgm(ig)*skk*becsum(ijh,na,is)
enddo
enddo
#ifdef DEBUG_ADDUSDENS
call stop_clock ('addus:aux')
#endif
endif
enddo
enddo
enddo
endif
enddo
!
deallocate (ylmk0)
deallocate (qmod)
!
! convert aux to real space and add to the charge density
!
do is = 1, nspin
psic(:) = (0.d0, 0.d0)
psic( nl(:) ) = aux(:,is)
if (gamma_only) psic( nlm(:) ) = conjg(aux(:,is))
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
rho (:, is) = rho (:, is) + DREAL (psic (:) )
enddo
deallocate (aux)
call stop_clock ('addusdens')
return
end subroutine addusdens

70
Gamma/allocate_fft.f90
View File

@ -1,70 +0,0 @@
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------
subroutine allocate_fft
!-----------------------------------------------------------------------
! This routine computes the data structure associated to the FFT
! grid and allocate memory for all the arrays which depend upon
! these dimensions
!
#include "machine.h"
use pwcom
USE wavefunctions, ONLY : psic
implicit none
!
! determines the data structure for fft arrays
!
call data_structure( gamma_only )
!
if (nrxx.lt.ngm) then
write (6, '(/,4x," nr1=",i4," nr2= ", i4, " nr3=",i4, &
&" nrxx = ",i8," ngm=",i8)') nr1, nr2, nr3, nrxx, ngm
call errore ('allocate_fft', 'the nr"s are too small!', 1)
endif
if (nrxxs.lt.ngms) then
write (6, '(/,4x," nr1s=",i4," nr2s= ", i4, " nr3s=",i4, &
&" nrxxs = ",i8," ngms=",i8)') nr1s, nr2s, nr3s, nrxxs, ngms
call errore ('allocate_fft', 'the nrs"s are too small!', 1)
endif
if (ngm <= 0) call errore ('allocate_fft', 'wrong ngm', 1)
if (ngms <= 0) call errore ('allocate_fft', 'wrong ngms', 1)
if (nrxx <= 0) call errore ('allocate_fft', 'wrong nrxx', 1)
if (nrxxs<= 0) call errore ('allocate_fft', 'wrong nrxxs', 1)
if (nspin<= 0) call errore ('allocate_fft', 'wrong nspin', 1)
!
! Allocate memory for all kind of stuff.
!
allocate (g( 3, ngm))
allocate (gg( ngm))
allocate (nl( ngm))
if (gamma_only) allocate (nlm(ngm))
allocate (igtongl( ngm))
allocate (ig1( ngm))
allocate (ig2( ngm))
allocate (ig3( ngm))
allocate (rho( nrxx, nspin))
allocate (rho_save( nrxx, nspin))
allocate (vr( nrxx,nspin))
allocate (vltot( nrxx))
allocate (vnew ( nrxx, nspin))
allocate (rho_core( nrxx))
allocate (psic( nrxx))
allocate (vrs( nrxx, nspin))
if (doublegrid) then
allocate (nls( ngms))
if (gamma_only) allocate (nlsm(ngm))
else
nls => nl
if (gamma_only) nlsm=> nlm
endif
return
end subroutine allocate_fft

131
Gamma/atomic_rho.f90
View File

@ -1,131 +0,0 @@
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine atomic_rho (rhoa, nspina)
!-----------------------------------------------------------------------
! This routine calculates rhoa as the superposition of atomic charges.
!
! nspina is the number of spin components to be calculated
!
! if nspina = 1 the total atomic charge density is calculated
! if nspina = 2 the spin up and spin down atomic charge densities are
! calculated assuming an uniform atomic spin-polarization
! equal to starting_magnetization(nt)
!
! NB: nspina may not be equal to nspin because in some cases (as in update)
! the total charge only could be needed, even in a LSDA calculation.
!
!
#include "machine.h"
use pwcom
USE wavefunctions, ONLY : psic
implicit none
integer :: nspina
! the number of spin polarizations
real(kind=DP) :: rhoa (nrxx, nspina)
! the output atomic charge
!
! local variables
!
real(kind=DP) :: rhoneg, rhorea, rhoima, gx
real(kind=DP), allocatable :: rhocgnt (:), aux (:)
complex(kind=DP), allocatable :: rhocg (:,:)
integer :: ir, is, ig, igl, nt
!
! superposition of atomic charges contained in the array rho_at
! (read from pseudopotential files)
!
! allocate work space (psic must already be allocated)
!
allocate (rhocg( ngm, nspina))
allocate (aux( ndm))
allocate (rhocgnt( ngl))
rhoa(:,:) = 0.d0
rhocg(:,:) = (0.d0,0.d0)
do nt = 1, ntyp
!
! Here we compute the G=0 term
!
if (gstart == 2) then
do ir = 1, msh (nt)
aux (ir) = rho_at (ir, nt)
enddo
call simpson (msh (nt), aux, rab (1, nt), rhocgnt (1) )
endif
!
! Here we compute the G<>0 term
!
do igl = gstart, ngl
gx = sqrt (gl (igl) ) * tpiba
do ir = 1, msh (nt)
if (r (ir, nt) .lt.1.0d-8) then
aux(ir) = rho_at(ir,nt)
else
aux(ir) = rho_at(ir,nt) * sin(gx*r(ir,nt)) / (r(ir,nt)*gx)
endif
enddo
call simpson (msh (nt), aux, rab (1, nt), rhocgnt (igl) )
enddo
!
! we compute the 3D atomic charge in reciprocal space
!
if (nspina == 1) then
do ig = 1, ngm
rhocg(ig,1) = rhocg(ig,1) + &
strf(ig,nt) * rhocgnt(igtongl(ig)) / omega
enddo
else
do ig = 1, ngm
rhocg(ig,1) = rhocg(ig,1) + &
0.5d0 * ( 1.d0 + starting_magnetization(nt) ) * &
strf(ig,nt) * rhocgnt(igtongl(ig)) / omega
rhocg(ig,2) = rhocg(ig,2) + &
0.5d0 * ( 1.d0 - starting_magnetization(nt) ) * &
strf(ig,nt) * rhocgnt(igtongl(ig)) / omega
enddo
endif
enddo
deallocate (rhocgnt)
deallocate (aux)
do is = 1, nspina
!
! and we return to real space
!
psic(:) = (0.d0,0.d0)
psic (nl (:) ) = rhocg (:, is)
if (gamma_only) psic ( nlm(:) ) = conjg( rhocg (:, is) )
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! we check that everything is correct
!
rhoneg = 0.d0
rhoima = 0.d0
do ir = 1, nrxx
rhorea = DREAL (psic (ir) )
rhoneg = rhoneg + min (0.d0, rhorea)
rhoima = rhoima + abs (DIMAG (psic (ir) ) )
rhoa (ir, is) = rhorea
enddo
rhoneg = rhoneg / (nr1 * nr2 * nr3)
rhoima = rhoima / (nr1 * nr2 * nr3)
#ifdef __PARA
call reduce (1, rhoneg)
call reduce (1, rhoima)
#endif
if ( rhoneg < -1.0d-4 .or. rhoima > 1.0d-4 ) &
write (6,'(/" Warning: negative or imaginary starting charge ",&
&2f12.6,i3)') rhoneg, rhoima, is
enddo
deallocate (rhocg)
return
end subroutine atomic_rho

8
Gamma/fake.f90
View File

@ -26,3 +26,11 @@ subroutine stres_hub
call errore('Gamma','stres_hub not implemented',1)
end subroutine stres_hub
subroutine c_phase
call errore('Gamma','Berry Phase not implemented',1)
end subroutine c_phase
subroutine kp_strings
call errore('Gamma','Berry Phase not implemented',2)
end subroutine kp_strings

264
Gamma/gradcorr.f90
View File

@ -1,264 +0,0 @@
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------
subroutine gradcorr (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
nrx3, nrxx, nl, ngm, g, alat, omega, nspin, etxc, vtxc, v)
! ===================
!--------------------------------------------------------------------
#include "machine.h"
use parameters
use funct
implicit none
!
integer :: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, nl (ngm), &
nspin
real(kind=DP) :: rho (nrxx, nspin), rho_core (nrxx), v (nrxx, nspin), &
g (3, ngm), vtxc, etxc, alat, omega, zeta, rh, grh2
integer :: k, ipol, is
real(kind=DP), allocatable :: grho (:,:,:), h (:,:,:), dh (:)
real(kind=DP) :: grho2 (2), sx, sc, v1x, v2x, v1c, v2c, v1xup, v1xdw, &
v2xup, v2xdw, v1cup, v1cdw , etxcgc, vtxcgc, segno, arho, fac
real(kind=DP), parameter :: e2 = 2.d0, epsr = 1.0d-6, epsg = 1.0d-10
if (igcx == 0 .and. igcc == 0) return
etxcgc = 0.d0
vtxcgc = 0.d0
allocate (h( 3, nrxx, nspin))
allocate (grho( 3, nrxx, nspin))
! calculate the gradient of rho+rho_core in real space
fac = 1.d0 / float (nspin)
do is = 1, nspin
call DAXPY (nrxx, fac, rho_core, 1, rho (1, is), 1)
call gradient (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, rho (1, is), &
ngm, g, nl, alat, grho (1, 1, is) )
enddo
do k = 1, nrxx
do is = 1, nspin
grho2 (is) = grho(1, k, is)**2 + grho(2, k, is)**2 + grho(3, k, is)**2
enddo
if (nspin == 1) then
!
! This is the spin-unpolarised case
!
arho = abs (rho (k, 1) )
segno = sign (1.d0, rho (k, 1) )
if (arho.gt.epsr.and.grho2 (1) .gt.epsg) then
call gcxc (arho, grho2, sx, sc, v1x, v2x, v1c, v2c)
!
! first term of the gradient correction : D(rho*Exc)/D(rho)
v (k, 1) = v (k, 1) + e2 * (v1x + v1c)
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
do ipol = 1, 3
h (ipol, k, 1) = e2 * (v2x + v2c) * grho (ipol, k, 1)
enddo
vtxcgc = vtxcgc + e2 * (v1x + v1c) * (rho (k, 1) - rho_core(k) )
etxcgc = etxcgc + e2 * (sx + sc) * segno
else
do ipol = 1, 3
h (ipol, k, 1) = 0.d0
enddo
endif
else
!
! spin-polarised case
!
call gcx_spin (rho (k, 1), rho (k, 2), grho2 (1), grho2 (2), &
sx, v1xup, v1xdw, v2xup, v2xdw)
rh = rho (k, 1) + rho (k, 2)
if (rh.gt.epsr) then
zeta = (rho (k, 1) - rho (k, 2) ) / rh
grh2 = (grho (1, k, 1) + grho (1, k, 2) ) **2 + &
(grho (2, k, 1) + grho (2, k, 2) ) **2 + &
(grho (3, k, 1) + grho (3, k, 2) ) **2
call gcc_spin (rh, zeta, grh2, sc, v1cup, v1cdw, v2c)
else
sc = 0.d0
v1cup = 0.d0
v1cdw = 0.d0
v2c = 0.d0
endif
!
! first term of the gradient correction : D(rho*Exc)/D(rho)
!
v (k, 1) = v (k, 1) + e2 * (v1xup + v1cup)
v (k, 2) = v (k, 2) + e2 * (v1xdw + v1cdw)
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
do ipol = 1, 3
h (ipol, k, 1) = e2 * ( (v2xup + v2c) * grho (ipol, k, 1) &
+ v2c * grho (ipol, k, 2) )
h (ipol, k, 2) = e2 * ( (v2xdw + v2c) * grho (ipol, k, 2) &
+ v2c * grho (ipol, k, 1) )
enddo
vtxcgc = vtxcgc + e2 * (v1xup + v1cup) * (rho (k, 1) - &
rho_core (k) * fac)
vtxcgc = vtxcgc + e2 * (v1xdw + v1cdw) * (rho (k, 2) - &
rho_core (k) * fac)
etxcgc = etxcgc + e2 * (sx + sc)
endif
enddo
do is = 1, nspin
call DAXPY (nrxx, - fac, rho_core, 1, rho (1, is), 1)
enddo
deallocate(grho)
allocate (dh( nrxx))
!
! second term of the gradient correction :
! \sum_alpha (D / D r_alpha) ( D(rho*Exc)/D(grad_alpha rho) )
!
do is = 1, nspin
call grad_dot (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, h (1, 1, is), &
ngm, g, nl, alat, dh)
do k = 1, nrxx
v (k, is) = v (k, is) - dh (k)
vtxcgc = vtxcgc - dh (k) * rho (k, is)
enddo
enddo
vtxc = vtxc + omega * vtxcgc / (nr1 * nr2 * nr3)
etxc = etxc + omega * etxcgc / (nr1 * nr2 * nr3)
deallocate (dh)
deallocate (h)
return
end subroutine gradcorr
!--------------------------------------------------------------------
subroutine gradient (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
ngm, g, nl, alat, ga)
!--------------------------------------------------------------------
!
! Calculates ga = \grad a in R-space (a is also in R-space)
use parameters
use gvect, only: nlm
use wvfct, only: gamma_only
implicit none
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, ngm, nl (ngm)
real(kind=DP) :: a (nrxx), g (3, ngm), ga (3, nrxx), alat
integer :: n, ipol
real(kind=DP), allocatable :: aux (:,:), gaux (:,:)
real(kind=DP) :: tpi, tpiba
parameter (tpi = 2.d0 * 3.14159265358979d0)
allocate (aux( 2,nrxx))
allocate (gaux(2,nrxx))
tpiba = tpi / alat
!
! copy a(r) to complex array...
!
aux(2,:) = 0.d0
call DCOPY (nrxx, a, 1, aux, 2)
!
! bring a(r) to G-space, a(G) ...
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
!
! multiply by (iG) to get (\grad_ipol a)(G) ...
!
ga(:,:) = 0.d0
do ipol = 1, 3
gaux(:,:) = 0.d0
do n = 1, ngm
gaux (1, nl (n) ) = - g (ipol, n) * aux (2, nl (n) )
gaux (2, nl (n) ) = g (ipol, n) * aux (1, nl (n) )
enddo
if (gamma_only) then
do n = 1, ngm
gaux (1, nlm(n) ) = gaux (1, nl(n) )
gaux (2, nlm(n) ) = - gaux (2, nl(n) )
enddo
end if
!
! bring back to R-space, (\grad_ipol a)(r) ...
!
call cft3 (gaux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! ...and add the factor 2\pi/a missing in the definition of G
!
call DAXPY (nrxx, tpiba, gaux, 2, ga (ipol, 1), 3)
enddo
deallocate (gaux)
deallocate (aux)
return
end subroutine gradient
!--------------------------------------------------------------------
subroutine grad_dot (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
ngm, g, nl, alat, da)
!--------------------------------------------------------------------
!
! Calculates da = \sum_i \grad_i a_i in R-space
use parameters
use gvect, only: nlm
use wvfct, only: gamma_only
implicit none
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, ngm, nl (ngm)
real(kind=DP) :: a (3, nrxx), g (3, ngm), da (nrxx), alat
integer :: n, ipol
real(kind=DP), allocatable :: aux (:,:), gaux (:,:)
real(kind=DP) :: tpi, tpiba
parameter (tpi = 2.d0 * 3.14159265358979d0)
allocate (aux( 2,nrxx))
allocate (gaux(2,nrxx))
gaux(:,:) = 0.d0
do ipol = 1, 3
!
! copy a(ipol,r) to a complex array...
!
aux(2,:) = 0.d0
call DCOPY (nrxx, a (ipol, 1), 3, aux, 2)
!
! bring a(ipol,r) to G-space, a(G) ...
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
!
! multiply by (iG) to get (\grad_ipol a)(G) ...
!
do n = 1, ngm
gaux (1, nl (n) ) = gaux (1, nl (n) ) - g (ipol, n) * aux (2,nl(n))
gaux (2, nl (n) ) = gaux (2, nl (n) ) + g (ipol, n) * aux (1,nl(n))
enddo
enddo
if (gamma_only) then
do n = 1, ngm
gaux (1, nlm(n) ) = gaux (1, nl (n) )
gaux (2, nlm(n) ) = - gaux (2, nl (n) )
enddo
end if
!
! bring back to R-space, (\grad_ipol a)(r) ...
!
call cft3 (gaux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! ...add the factor 2\pi/a missing in the definition of G and sum
!
tpiba = tpi / alat
do n=1,nrxx
da(n) = gaux(1,n)*tpiba
end do
!
deallocate (gaux)
deallocate (aux)
return
end subroutine grad_dot

165
Gamma/interpolate.f90
View File

@ -1,165 +0,0 @@
!
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
subroutine interpolate (v, vs, iflag)
!
! This subroutine interpolates :
! from the smooth mesh (vs) to a thicker mesh (v) (iflag>0)
! vs is unchanged on output
! from the thick mesh (v ) to a smoother mesh (vs) (iflag<=0)
! v is unchanged on output
! V and Vs are real and in real space . V and Vs may coincide
!
#include"machine.h"
USE parameters, ONLY: DP
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm
USE gsmooth,ONLY: nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,nrxxs,ngms, &
nls, nlsm, doublegrid
implicit none
real(kind=DP) :: v (nrxx), vs (nrxxs)
! function on thick mesh
! function on smooth mesh
complex(kind=DP), allocatable :: aux (:), auxs (:)
! work array on thick mesh
! work array on smooth mesh
integer :: iflag
! gives the direction of the interpolation
integer :: ig, ir
call start_clock ('interpolate')
if (iflag <= 0) then
!
! from thick to smooth
!
if (doublegrid) then
allocate (aux( nrxx))
allocate (auxs(nrxxs))
aux (:) = v (:)
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
auxs (:) = (0.d0, 0.d0)
do ig = 1, ngms
auxs (nls (ig) ) = aux (nl (ig) )
enddo
if (gamma_only) then
do ig = 1, ngms
auxs (nlsm(ig) ) = aux (nlm(ig) )
enddo
end if
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 1)
vs (:) = auxs (:)
deallocate (auxs)
deallocate (aux)
else
do ir = 1, nrxx
vs (ir) = v (ir)
enddo
endif
else
!
! from smooth to thick
!
if (doublegrid) then
allocate (aux( nrxx))
allocate (auxs(nrxxs))
auxs (:) = vs (:)
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 1)
aux (:) = (0.d0, 0.d0)
do ig = 1, ngms
aux (nl (ig) ) = auxs (nls (ig) )
enddo
if (gamma_only) then
do ig = 1, ngms
aux (nlm(ig) ) = aux (nlsm(ig) )
enddo
end if
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
v (:) = aux (:)
deallocate (auxs)
deallocate (aux)
else
do ir = 1, nrxx
v (ir) = vs (ir)
enddo
endif
endif
call stop_clock ('interpolate')
return
end subroutine interpolate
!
subroutine cinterpolate (v, vs, iflag)
!
! This subroutine interpolates :
! from the smooth mesh (vs) to a thicker mesh (v) (iflag>0)
! vs is unchanged on output
! from the thick mesh (v ) to a smoother mesh (vs) (iflag<=0)
! v is unchanged on output
! V and Vs are complex and in real space . V and Vs may coincide
!
#include"machine.h"
USE parameters, ONLY: DP
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm
USE gsmooth,ONLY: nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,nrxxs,ngms, &
nls, nlsm, doublegrid
complex(kind=DP) :: v (nrxx), vs (nrxxs)
! function on thick mesh
! function on smooth mesh
integer :: iflag
! gives the direction of the interpolation
complex(kind=DP), allocatable :: aux (:), auxs (:)
! work array on thick mesh
! work array on smooth mesh
integer :: ig
if (gamma_only) call errore ('cinterpolate','not implemented', 1)
call start_clock ('cinterpolate')
if (iflag <= 0) then
!
! from thick to smooth
!
if (doublegrid) then
allocate (aux ( nrxx))
aux (:) = v(:)
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
vs (:) = (0.d0, 0.d0)
do ig = 1, ngms
vs (nls (ig) ) = aux (nl (ig) )
enddo
call cft3s (vs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 1)
deallocate (aux)
else
call ZCOPY (nrxx, v, 1, vs, 1)
endif
else
!
! from smooth to thick
!
if (doublegrid) then
allocate (auxs ( nrxxs))
auxs (:) = vs(:)
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, -1)
v (:) = (0.d0, 0.d0)
do ig = 1, ngms
v (nl (ig) ) = auxs (nls (ig) )
enddo
call cft3 (v, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
deallocate (auxs)
else
call ZCOPY (nrxx, vs, 1, v, 1)
endif
endif
call stop_clock ('cinterpolate')
return
end subroutine cinterpolate

831
Gamma/mix_rho.f90
View File

@ -1,831 +0,0 @@
!
! Copyright (C) 2002-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
#include "machine.h"
#undef DEBUG
!-----------------------------------------------------------------------
subroutine mix_rho (rhout, rhoin, nsout, nsin, alphamix, dr2, iter, &
n_iter, filename, conv)
!-----------------------------------------------------------------------
!
! Modified Broyden's method for charge density mixing
! d.d. johnson prb 38, 12807 (1988)
! On output: the mixed density is in rhoin, rhout is UNCHANGED
!
use parameters, only : DP
use pwcom
USE wavefunctions, ONLY : psic
!
! First the I/O variable
!
character (len=42) :: &
filename ! (in) I/O filename for mixing history
! if absent everything is kept in memory
integer :: &
iter, &! (in) counter of the number of iterations
n_iter ! (in) numb. of iterations used in mixing
real (kind=DP) :: &
rhout(nrxx,nspin), &! (in) the "out" density; (out) rhout-rhoin
rhoin(nrxx,nspin), &! (in) the "in" density; (out) the new dens.
nsout(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat), &!
nsin(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat), &!
alphamix, &! (in) mixing factor
dr2 ! (out) the estimated errr on the energy
logical :: &
conv ! (out) if true the convergence has been reached
!
integer, parameter :: &
maxmix =25 ! max number of iterations for charge mixing
!
! Here the local variables
!
integer :: &
iunmix, &! I/O unit number of charge density file
iunmix2, &! I/O unit number of ns file
iunit, &! counter on I/O unit numbers
iter_used, &! actual number of iterations used
ipos, &! index of the present iteration
inext, &! index of the next iteration
i, j, &! counters on number of iterations
is, &! counter on spin component
ig, &! counter on G-vectors
iwork(maxmix),&! dummy array used as output by libr. routines
info ! flag saying if the exec. of libr. routines was ok
complex (kind=DP), allocatable :: rhocin(:,:), rhocout(:,:), &
rhoinsave(:), rhoutsave(:), &
nsinsave(:,:,:,:), nsoutsave(:,:,:,:)
complex (kind=DP), allocatable, save :: df(:,:), dv(:,:), &
df_ns(:,:,:,:,:), dv_ns(:,:,:,:,:)
! rhocin(ngm0,nspin)
! rhocout(ngm0,nspin)
! rhoinsave(ngm0*nspin): work space
! rhoutsave(ngm0*nspin): work space
! df(ngm0*nspin,n_iter): information from preceding iterations
! dv(ngm0*nspin,n_iter): " " " " " "
! df_ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat,n_iter):idem
! dv_ns(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat,n_iter):idem
integer :: ldim
real (kind=DP) :: betamix(maxmix,maxmix), gamma0, work(maxmix), dehar
logical :: &
saveonfile, &! save intermediate steps on file "filename"
opnd, &! if true the file is already opened
exst ! if true the file exists
real (kind=DP), external :: rho_dot_product, ns_dot_product, fn_dehar
call start_clock('mix_rho')
if (iter < 1) call errore('mix_rho','iter is wrong',1)
if (n_iter > maxmix) call errore('mix_rho','n_iter too big',1)
if (lda_plus_u) ldim = 2 * Hubbard_lmax + 1
saveonfile=filename.ne.' '
! call DAXPY(nrxx*nspin,-1.d0,rhoin,1,rhout,1)
allocate(rhocin(ngm0,nspin), rhocout(ngm0,nspin))
!
! psic is used as work space - must be already allocated !
!
do is=1,nspin
psic(:) = DCMPLX (rhoin(:,is), 0.d0)
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
do ig=1,ngm0
rhocin(ig,is) = psic(nl(ig))
end do
psic(:) = DCMPLX (rhout(:,is), 0.d0)
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
do ig=1,ngm0
rhocout(ig,is) = psic(nl(ig)) - rhocin(ig,is)
end do
end do
if (lda_plus_u) nsout(:,:,:,:) = nsout(:,:,:,:) - nsin(:,:,:,:)
dr2=rho_dot_product(rhocout,rhocout) + ns_dot_product(nsout,nsout)
conv = (dr2 < tr2)
dehar = fn_dehar(rhocout)
#ifdef DEBUG
! if (lda_plus_u) write (6,*) ' ns_dr2 =', ns_dot_product(nsout,nsout)
if (conv) then
write (6,100) dr2, rho_dot_product(rhocout,rhocout) + &
ns_dot_product(nsout,nsout)
write (6,'(" dehar =",f15.8)') dehar
end if
#endif
if (saveonfile) then
do iunit=99,1,-1
inquire(unit=iunit,opened=opnd)
iunmix=iunit
if (.not.opnd) go to 10
end do
call errore('mix_rho','free unit not found?!?',1)
10 continue
if (lda_plus_u) then
do iunit=iunmix-1,1,-1
inquire(unit=iunit,opened=opnd)
iunmix2=iunit
if (.not.opnd) go to 20
end do
call errore('mix_rho','second free unit not found?!?',1)
20 continue
end if
if (conv) then
call diropn (iunmix, filename, 2*ngm0*nspin, exst)
close (unit=iunmix, status='delete')
if (lda_plus_u) then
call diropn (iunmix2, trim(filename)//'.ns',ldim*ldim*nspin*nat, exst)
close (unit=iunmix2, status='delete')
end if
deallocate (rhocin, rhocout)
call stop_clock('mix_rho')
return
end if
call diropn(iunmix,filename,2*ngm0*nspin,exst)
if (lda_plus_u) call diropn (iunmix2, trim(filename)//'.ns',ldim*ldim*nspin*nat, exst)
if (iter > 1 .and. .not.exst) then
call errore('mix_rho','file not found, restarting',-1)
iter=1
end if
allocate (df(ngm0*nspin,n_iter), dv(ngm0*nspin,n_iter))
if (lda_plus_u) &
allocate (df_ns(ldim,ldim,nspin,nat,n_iter), &
dv_ns(ldim,ldim,nspin,nat,n_iter))
else
if (iter == 1) then
allocate (df(ngm0*nspin,n_iter), dv(ngm0*nspin,n_iter))
if (lda_plus_u) &
allocate (df_ns(ldim,ldim,nspin,nat,n_iter),&
dv_ns(ldim,ldim,nspin,nat,n_iter))
end if
if (conv) then
if (lda_plus_u) deallocate(df_ns, dv_ns)
deallocate (df, dv)
deallocate (rhocin, rhocout)
call stop_clock('mix_rho')
return
end if
allocate (rhoinsave(ngm0*nspin), rhoutsave(ngm0*nspin))
if (lda_plus_u) &
allocate(nsinsave (ldim,ldim,nspin,nat), &
nsoutsave(ldim,ldim,nspin,nat))
end if
!
! copy only the high frequency Fourier component into rhoin
! (NB: rhout=rhout-rhoin)
!
call DCOPY(nrxx*nspin,rhout,1,rhoin,1)
do is=1,nspin
psic(:) = (0.d0, 0.d0)
do ig=1,ngm0
psic(nl(ig)) = rhocin(ig,is)+rhocout(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg ( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
call DAXPY(nrxx,-1.d0,psic,2,rhoin(1,is),1)
end do
!
! iter_used = iter-1 if iter <= n_iter
! iter_used = n_iter if iter > n_iter
!
iter_used=min(iter-1,n_iter)
!
! ipos is the position in which results from the present iteration
! are stored. ipos=iter-1 until ipos=n_iter, then back to 1,2,...
!
ipos =iter-1-((iter-2)/n_iter)*n_iter
!
if (iter > 1) then
if (saveonfile) then
call davcio(df(1,ipos),2*ngm0*nspin,iunmix,1,-1)
call davcio(dv(1,ipos),2*ngm0*nspin,iunmix,2,-1)
if (lda_plus_u) then
call davcio(df_ns(1,1,1,1,ipos),ldim*ldim*nspin*nat,iunmix2,1,-1)
call davcio(dv_ns(1,1,1,1,ipos),ldim*ldim*nspin*nat,iunmix2,2,-1)
end if
end if
call DAXPY(2*ngm0*nspin,-1.d0,rhocout,1,df(1,ipos),1)
call DAXPY(2*ngm0*nspin,-1.d0,rhocin ,1,dv(1,ipos),1)
! norm = sqrt(rho_dot_product(df(1,ipos),df(1,ipos)) + &
! ns_dot_product(df_ns(1,1,1,1,ipos),df_ns(1,1,1,1,ipos)) )
! call DSCAL (2*ngm0*nspin,-1.d0/norm,df(1,ipos),1)
! call DSCAL (2*ngm0*nspin,-1.d0/norm,dv(1,ipos),1)
if (lda_plus_u) then
call DAXPY(ldim*ldim*nspin*nat,-1.d0,nsout,1,df_ns(1,1,1,1,ipos),1)
call DAXPY(ldim*ldim*nspin*nat,-1.d0,nsin ,1,dv_ns(1,1,1,1,ipos),1)
end if
end if
!
if (saveonfile) then
do i=1,iter_used
if (i.ne.ipos) then
call davcio(df(1,i),2*ngm0*nspin,iunmix,2*i+1,-1)
call davcio(dv(1,i),2*ngm0*nspin,iunmix,2*i+2,-1)
if (lda_plus_u) then
call davcio(df_ns(1,1,1,1,i),ldim*ldim*nspin*nat,iunmix2,2*i+1,-1)
call davcio(dv_ns(1,1,1,1,i),ldim*ldim*nspin*nat,iunmix2,2*i+2,-1)
end if
end if
end do
call davcio(rhocout,2*ngm0*nspin,iunmix,1,1)
call davcio(rhocin ,2*ngm0*nspin,iunmix,2,1)
if (iter > 1) then
call davcio(df(1,ipos),2*ngm0*nspin,iunmix,2*ipos+1,1)
call davcio(dv(1,ipos),2*ngm0*nspin,iunmix,2*ipos+2,1)
end if
if (lda_plus_u) then
call davcio(nsout,ldim*ldim*nspin*nat,iunmix2,1,1)
call davcio(nsin ,ldim*ldim*nspin*nat,iunmix2,2,1)
if (iter > 1) then
call davcio(df_ns(1,1,1,1,ipos),ldim*ldim*nspin*nat,iunmix2,2*ipos+1,1)
call davcio(dv_ns(1,1,1,1,ipos),ldim*ldim*nspin*nat,iunmix2,2*ipos+2,1)
end if
end if
else
call DCOPY(2*ngm0*nspin,rhocin ,1,rhoinsave,1)
call DCOPY(2*ngm0*nspin,rhocout,1,rhoutsave,1)
if (lda_plus_u) then
call DCOPY(ldim*ldim*nspin*nat,nsin ,1,nsinsave ,1)
call DCOPY(ldim*ldim*nspin*nat,nsout,1,nsoutsave,1)
end if
end if
!
do i=1,iter_used
do j=i,iter_used
betamix(i,j) = rho_dot_product(df(1,j),df(1,i)) + &
ns_dot_product(df_ns(1,1,1,1,j),df_ns(1,1,1,1,i))
end do
end do
!
call DSYTRF ('U',iter_used,betamix,maxmix,iwork,work,maxmix,info)
call errore('broyden','factorization',info)
call DSYTRI ('U',iter_used,betamix,maxmix,iwork,work,info)
call errore('broyden','DSYTRI',info)
!
do i=1,iter_used
do j=i+1,iter_used
betamix(j,i)=betamix(i,j)
end do
end do
!
do i=1,iter_used
work(i) = rho_dot_product(df(1,i),rhocout) + &
ns_dot_product(df_ns(1,1,1,1,i),nsout)
end do
!
do i=1,iter_used
gamma0=0.d0
do j=1,iter_used
gamma0 = gamma0 + betamix(j,i)*work(j)
end do
call DAXPY(2*ngm0*nspin,-gamma0,dv(1,i),1,rhocin,1)
call DAXPY(2*ngm0*nspin,-gamma0,df(1,i),1,rhocout,1)
if (lda_plus_u) then
call DAXPY(ldim*ldim*nspin*nat,-gamma0,dv_ns(1,1,1,1,i),1,nsin(1,1,1,1) ,1)
call DAXPY(ldim*ldim*nspin*nat,-gamma0,df_ns(1,1,1,1,i),1,nsout(1,1,1,1),1)
end if
end do
!
#ifdef DEBUG
write (6,100) dr2, rho_dot_product(rhocout,rhocout) + &
ns_dot_product(nsout,nsout)
write (6,'(" dehar =",f15.8)') dehar
#endif
100 format (' dr2 =',1pe15.1, ' internal_best_dr2= ', 1pe15.1)
! - auxiliary vectors dv and df not needed anymore
if (saveonfile) then
if (lda_plus_u) then
close(iunmix2,status='keep')
deallocate (df_ns, dv_ns)
end if
close(iunmix,status='keep')
deallocate (df, dv)
else
inext=iter-((iter-1)/n_iter)*n_iter
if (lda_plus_u) then
call DCOPY(ldim*ldim*nspin*nat,nsoutsave,1,df_ns(1,1,1,1,inext),1)
call DCOPY(ldim*ldim*nspin*nat,nsinsave ,1,dv_ns(1,1,1,1,inext),1)
deallocate (nsinsave, nsoutsave)
end if
call DCOPY(2*ngm0*nspin,rhoutsave,1,df(1,inext),1)
call DCOPY(2*ngm0*nspin,rhoinsave,1,dv(1,inext),1)
deallocate (rhoinsave, rhoutsave)
end if
! - preconditioning the new search direction (if imix.gt.0)
if (imix == 1) then
call approx_screening(rhocout)
else if (imix == 2) then
call approx_screening2(rhocout,rhocin)
end if
! - set new trial density
call DAXPY(2*ngm0*nspin,alphamix,rhocout,1,rhocin,1)
do is=1,nspin
psic(:) = (0.d0,0.d0)
do ig=1,ngm0
psic(nl(ig)) = rhocin(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg ( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
call DAXPY(nrxx,1.d0,psic,2,rhoin(1,is),1)
end do
if (lda_plus_u) call DAXPY(ldim*ldim*nspin*nat,alphamix,nsout,1,nsin,1)
! - clean up
deallocate(rhocout)
deallocate(rhocin)
call stop_clock('mix_rho')
return
end subroutine mix_rho
!
!--------------------------------------------------------------------
function rho_dot_product (rho1,rho2)
!--------------------------------------------------------------------
! this function evaluates the dot product between two input densities
!
use parameters, only : DP
use pwcom
!
! I/O variables
!
real (kind=DP) :: rho_dot_product ! (out) the function value
complex (kind=DP) :: rho1(ngm0,nspin), rho2(ngm0,nspin) ! (in) the two densities
!
! and the local variables
!
real (kind=DP) :: fac ! a multiplicative factors
integer :: is, ig
rho_dot_product = 0.d0
if (nspin == 1) then
is=1
do ig = gstart,ngm0
fac = e2*fpi / (tpiba2*gg(ig))
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(ig,is))*rho2(ig,is))
end do
if (gamma_only) rho_dot_product = 2.d0*rho_dot_product
else
do ig = gstart,ngm0
fac = e2*fpi / (tpiba2*gg(ig))
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(ig,1)+rho1(ig,2))* &
(rho2(ig,1)+rho2(ig,2)))
end do
if (gamma_only) rho_dot_product = 2.d0*rho_dot_product
fac = e2*fpi / (tpi**2) ! lambda=1 a.u.
! G=0 term
if (gstart == 2) then
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(1,1)-rho1(1,2))* &
(rho2(1,1)-rho2(1,2)))
end if
if (gamma_only) fac = 2.d0*fac
do ig = gstart,ngm0
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(ig,1)-rho1(ig,2))* &
(rho2(ig,1)-rho2(ig,2)))
end do
end if
rho_dot_product = rho_dot_product * omega / 2.d0
#ifdef __PARA
call reduce(1,rho_dot_product)
#endif
return
end function rho_dot_product
!
!--------------------------------------------------------------------
function ns_dot_product (ns1,ns2)
!--------------------------------------------------------------------
! this function evaluates the dot product between two input densities
!
use parameters, only : DP
use pwcom
!
! I/O variables
!
real (kind=DP) :: ns_dot_product ! (out) the function value
real (kind=DP) :: ns1(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat), &
ns2(2*Hubbard_lmax+1,2*Hubbard_lmax+1,nspin,nat)
! (in) the two ns
!
! and the local variables
!
real (kind=DP) :: sum
integer :: na, nt, is, m1, m2
ns_dot_product = 0.d0
if (.not. lda_plus_u ) return
do na = 1, nat
nt = ityp (na)
if (Hubbard_U(nt).ne.0.d0 .or. Hubbard_alpha(nt).ne.0.d0) then
sum =0.d0
do is = 1,nspin
do m1 = 1, 2 * Hubbard_l(nt) + 1
do m2 = m1, 2 * Hubbard_l(nt) + 1
sum = sum + ns1(m1,m2,is,na)*ns2(m2,m1,is,na)
enddo
enddo
end do
ns_dot_product = ns_dot_product + 0.5d0*Hubbard_U(nt) * sum
endif
end do
if (nspin == 1) ns_dot_product = 2.d0 * ns_dot_product
return
end function ns_dot_product
!--------------------------------------------------------------------
function fn_dehar (drho)
!--------------------------------------------------------------------
! this function evaluates the residual hartree energy of drho
!
use parameters, only : DP
use pwcom
!
! I/O variables
!
real (kind=DP) :: fn_dehar ! (out) the function value
complex (kind=DP) :: drho(ngm0,nspin) ! (in) the density difference
!
! and the local variables
!
real (kind=DP) :: fac ! a multiplicative factors
integer :: is, ig
fn_dehar = 0.d0
if (nspin == 1) then
is=1
do ig = gstart,ngm0
fac = e2*fpi / (tpiba2*gg(ig))
fn_dehar = fn_dehar + fac * abs(drho(ig,is))**2
end do
else
do ig = gstart,ngm0
fac = e2*fpi / (tpiba2*gg(ig))
fn_dehar = fn_dehar + fac * abs(drho(ig,1)+drho(ig,2))**2
end do
end if
if (gamma_only) then
fn_dehar = fn_dehar * omega
else
fn_dehar = fn_dehar * omega / 2.d0
end if
#ifdef __PARA
call reduce(1,fn_dehar)
#endif
return
end function fn_dehar
!--------------------------------------------------------------------
subroutine approx_screening (drho)
!--------------------------------------------------------------------
! apply an average TF preconditioning to drho
!
use parameters, only : DP
use pwcom
!
! I/O
!
complex (kind=DP) drho(ngm0,nspin) ! (in/out)
!
! and the local variables
!
real (kind=DP) :: rrho, rmag, rs, agg0
integer :: is, ig
rs = (3.d0*omega/fpi/nelec)**(1.d0/3.d0)
agg0 = (12.d0/pi)**(2.d0/3.d0)/tpiba2/rs
#ifdef DEBUG
write (6,'(a,f12.6,a,f12.6)') ' avg rs =', rs, ' avg rho =', nelec/omega
#endif
if (nspin == 1) then
is = 1
do ig = 1,ngm0
drho(ig,is) = drho(ig,is) * gg(ig)/(gg(ig)+agg0)
end do
else
do ig = 1,ngm0
rrho = (drho(ig,1) + drho(ig,2)) * gg(ig)/(gg(ig)+agg0)
rmag = (drho(ig,1) - drho(ig,2))
drho(ig,1) = 0.5d0 * (rrho + rmag)
drho(ig,2) = 0.5d0 * (rrho - rmag)
end do
end if
return
end subroutine approx_screening
!
!--------------------------------------------------------------------
subroutine approx_screening2 (drho,rhobest)
!--------------------------------------------------------------------
! apply a local-density dependent TF preconditioning to drho
!
use parameters, only : DP
use pwcom
USE wavefunctions, ONLY : psic
!
! I/O
!
!
complex (kind=DP) :: drho(ngm0,nspin), rhobest(ngm0,nspin)
!
! and the local variables
!
integer :: mmx
parameter (mmx=12)
integer :: iwork(mmx),i,j,m,info, nspin_save
real (kind=DP) :: rs, min_rs, max_rs, avg_rsm1, target, &
dr2_best, ccc, cbest, l2smooth
real (kind=DP) :: aa(mmx,mmx), invaa(mmx,mmx), bb(mmx), work(mmx), &
vec(mmx),agg0
complex (kind=DP) :: rrho, rmag
complex (kind=DP), allocatable :: v(:,:), w(:,:), dv(:), &
vbest(:), wbest(:)
! v(ngm0,mmx), w(ngm0,mmx), dv(ngm0), vbest(ngm0), wbest(ngm0)
real (kind=DP), allocatable :: alpha(:)
! alpha(nrxx)
integer :: is, ir, ig
real (kind=DP), external :: rho_dot_product
if (nspin == 2) then
do ig=1,ngm0
rrho = drho(ig,1) + drho(ig,2)
rmag = drho(ig,1) - drho(ig,2)
drho(ig,1) = rrho
drho(ig,2) = rmag
end do
end if
nspin_save = nspin
nspin = 1
is = 1
target = 0.d0
! write (6,*) ' eccoci qua '
if (gg(1) < 1.d-8) drho(1,is) = (0.d0,0.d0)
allocate (alpha(nrxx), v(ngm0,mmx), w(ngm0,mmx), &
dv(ngm0), vbest(ngm0), wbest(ngm0))
v(:,:) = (0.d0,0.d0)
w(:,:) = (0.d0,0.d0)
dv(:) = (0.d0,0.d0)
vbest(:)= (0.d0,0.d0)
wbest(:)= (0.d0,0.d0)
!
! - calculate alpha from density smoothed with a lambda=0 a.u.
!
l2smooth = 0.d0
psic(:) = (0.d0,0.d0)
if (nspin == 1) then
do ig=1,ngm0
psic(nl(ig)) = rhobest(ig,1) * exp(-0.5*l2smooth*tpiba2*gg(ig))
end do
else
do ig=1,ngm0
psic(nl(ig)) =(rhobest(ig,1) + rhobest(ig,2)) &
* exp(-0.5*l2smooth*tpiba2*gg(ig))
end do
end if
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
alpha(:) = real(psic(:))
min_rs = (3.d0*omega/fpi/nelec)**(1.d0/3.d0)
max_rs = min_rs
avg_rsm1 = 0.d0
do ir=1,nrxx
alpha(ir)=abs(alpha(ir))
rs = (3.d0/fpi/alpha(ir))**(1.d0/3.d0)
min_rs = min(min_rs,rs)
avg_rsm1 =avg_rsm1 + 1.d0/rs
max_rs = max(max_rs,rs)
alpha(ir) = rs
end do
#ifdef __PARA
call reduce (1, avg_rsm1)
call extreme (min_rs, -1)
call extreme (max_rs, +1)
#endif
call DSCAL(nrxx, 3.d0*(tpi/3.d0)**(5.d0/3.d0), alpha, 1)
avg_rsm1 = (nr1*nr2*nr3)/avg_rsm1
rs = (3.d0*omega/fpi/nelec)**(1.d0/3.d0)
agg0 = (12.d0/pi)**(2.d0/3.d0)/tpiba2/avg_rsm1
#ifdef DEBUG
write (6,'(a,5f12.6)') ' min/avgm1/max rs =', min_rs,avg_rsm1,max_rs,rs
#endif
!
! - calculate deltaV and the first correction vector
!
psic(:) = (0.d0,0.d0)
do ig=1,ngm0
psic(nl(ig)) = drho(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir) * alpha(ir)
end do
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
do ig=1,ngm0
dv(ig) = psic(nl(ig))*gg(ig)*tpiba2
v(ig,1)= psic(nl(ig))*gg(ig)/(gg(ig)+agg0)
end do
m=1
ccc = rho_dot_product(dv,dv)
aa(:,:) = 0.d0
bb(:) = 0.d0
3 continue
!
! - generate the vector w
!
do ig=1,ngm0
w(ig,m) = gg(ig)*tpiba2*v(ig,m)
end do
psic(:) = (0.d0,0.d0)
do ig=1,ngm0
psic(nl(ig)) = v(ig,m)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir)*fpi*e2/alpha(ir)
end do
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
do ig=1,ngm0
w(ig,m) = w(ig,m) + psic(nl(ig))
end do
!
! - build the linear system
!
do i=1,m
aa(i,m) = rho_dot_product(w(1,i),w(1,m))
aa(m,i) = aa(i,m)
end do
bb(m) = rho_dot_product(w(1,m),dv)
!
! - solve it -> vec
!
call DCOPY (mmx*mmx,aa,1,invaa,1)
call DSYTRF ('U',m,invaa,mmx,iwork,work,mmx,info)
call errore('BROYDEN','factorization',info)
call DSYTRI ('U',m,invaa,mmx,iwork,work,info)
call errore('broyden','DSYTRI',info)
!
do i=1,m
do j=i+1,m
invaa(j,i)=invaa(i,j)
end do
end do
do i=1,m
vec(i) = 0.d0
do j=1,m
vec(i) = vec(i) + invaa(i,j)*bb(j)
end do
end do
! -
vbest(:) = (0.d0,0.d0)
wbest(:) = dv(:)
do i=1,m
call DAXPY(2*ngm0, vec(i), v(1,i),1, vbest,1)
call DAXPY(2*ngm0,-vec(i), w(1,i),1, wbest,1)
end do
cbest = ccc
do i=1,m
cbest = cbest - bb(i)*vec(i)
end do
dr2_best= rho_dot_product(wbest,wbest)
if (target == 0.d0) target = 1.d-6 * dr2_best
! write (6,*) m, dr2_best, cbest
if (dr2_best < target) then
! write(6,*) ' last', dr2_best/target * 1.d-6
psic(:) = (0.d0,0.d0)
do ig=1,ngm0
psic(nl(ig)) = vbest(ig)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir)/alpha(ir)
end do
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,-1)
do ig=1,ngm0
drho(ig,is) = psic(nl(ig))
end do
nspin = nspin_save
if (nspin == 2) then
do ig=1,ngm0
rrho = drho(ig,1)
rmag = drho(ig,2)
drho(ig,1) = 0.5d0 * ( rrho + rmag )
drho(ig,2) = 0.5d0 * ( rrho - rmag )
end do
end if
deallocate (alpha, v, w, dv, vbest, wbest)
return
else if (m >= mmx) then
! write (6,*) m, dr2_best, cbest
m=1
do ig=1,ngm0
v(ig,m)=vbest(ig)
end do
aa(:,:) = 0.d0
bb(:) = 0.d0
go to 3
end if
m = m + 1
do ig=1,ngm0
v(ig,m)=wbest(ig)/(gg(ig)+agg0)
end do
go to 3
end subroutine approx_screening2

134
Gamma/set_rhoc.f90
View File

@ -1,134 +0,0 @@
!
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------
subroutine set_rhoc
!-----------------------------------------------------------------------
!
! This routine computes the core charge on the real space 3D mesh
!
!
#include "machine.h"
use pwcom
implicit none
!
real(kind=DP), parameter :: eps = 1.d-10
complex(kind=DP) , allocatable :: aux (:)
! used for the fft of the core charge
real(kind=DP) , allocatable :: rhocg(:)
! the radial fourier trasform
real(kind=DP) :: rhoima, rhoneg, rhorea
! used to check the core charge
real(kind=DP) :: vtxcc
! dummy xc energy term
real(kind=DP) , allocatable :: dum(:,:)
! dummy array containing rho=0
integer :: ir, nt, ng
! counter on mesh points
! counter on atomic types
! counter on g vectors
etxcc = 0.d0
do nt = 1, ntyp
if (nlcc (nt) ) goto 10
enddo
rho_core(:) = 0.d0
return
10 continue
allocate (aux( nrxx))
allocate (rhocg( ngl))
aux (:) = 0.d0
!
! the sum is on atom types
!
do nt = 1, ntyp
if (nlcc (nt) ) then
!
! drhoc compute the radial fourier transform for each shell of g vec
!
call drhoc (ngl, gl, omega, tpiba2, numeric (nt), a_nlcc (nt), &
b_nlcc (nt), alpha_nlcc (nt), msh (nt), r (1, nt), rab (1, nt), &
rho_atc (1, nt), rhocg)
!
! multiply by the structure factor and sum
!
do ng = 1, ngm
aux(nl(ng)) = aux(nl(ng)) + strf(ng,nt) * rhocg(igtongl(ng))
enddo
endif
enddo
if (gamma_only) then
do ng = 1, ngm
aux(nlm(ng)) = conjg(aux(nl (ng)))
end do
end if
!
! the core charge in real space
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! test on the charge and computation of the core energy
!
rhoneg = 0.d0
rhoima = 0.d0
do ir = 1, nrxx
rhoneg = rhoneg + min (0.d0, DREAL (aux (ir) ) )
rhoima = rhoima + abs (DIMAG (aux (ir) ) )
rho_core(ir) = DREAL (aux(ir))
!
! NOTE: Core charge is computed in reciprocal space and brought to real
! space by FFT. For non smooth core charges (or insufficient cut-off)
! this may result in negative values in some grid points.
! Up to October 1999 the core charge was forced to be positive definite.
! This induces an error in the force, and probably stress, calculation if
! the number of grid points where the core charge would be otherwise neg
! is large. The error disappears for sufficiently high cut-off, but may be
! rather large and it is better to leave the core charge as it is.
! If you insist to have it positive definite (with the possible problems
! mentioned above) uncomment the following lines. SdG, Oct 15 1999
!
! rhorea = max (DREAL (aux (ir) ), eps)
! rho_core(ir) = rhorea
!
enddo
rhoneg = rhoneg / (nr1 * nr2 * nr3)
rhoima = rhoima / (nr1 * nr2 * nr3)
#ifdef __PARA
call reduce (1, rhoneg)
call reduce (1, rhoima)
#endif
if (rhoneg.lt. - 1.0d-6.or.rhoima.gt.1.0d-6) &
write (6, '(" warning: negative or imaginary core charge ",2f12.6)')&
rhoneg, rhoima
!
! calculate core_only exch-corr energy etxcc=E_xc[rho_core] if required
! The term was present in previous versions of the code but it shouldn't
!
! allocate (dum(nrxx , nspin))
! dum(:,:) = 0.d0
! call v_xc (dum, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
! nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, &
! ehart, etxcc, vtxcc, aux)
! deallocate(dum)
! write (6, 9000) etxcc
! write (6, * ) 'BEWARE it will be subtracted from total energy !'
!
deallocate (rhocg)
deallocate (aux)
!
return
9000 format (5x,'core-only xc energy = ',f15.8,' ryd')
end subroutine set_rhoc

64
Gamma/setlocal.f90
View File

@ -1,64 +0,0 @@
!
! Copyright (C) 2001 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------
subroutine setlocal
!----------------------------------------------------------------------
!
! This routine computes the local potential in real space vltot(ir)
!
#include "machine.h"
USE parameters, ONLY: DP
USE basis, ONLY: ntyp
USE extfield, ONLY: tefield, dipfield
USE gvect, ONLY : nl, nlm, igtongl
USE scf, ONLY: vltot
USE vlocal, ONLY : strf, vloc
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm, ngm
implicit none
complex(kind=DP), allocatable :: aux (:)
! auxiliary variable
integer :: nt, ng, ir
! counter on atom types
! counter on g vectors
! counter on r vectors
!
allocate (aux( nrxx))
!
aux(:)=(0.d0,0.d0)
do nt = 1, ntyp
do ng = 1, ngm
aux (nl(ng))=aux(nl(ng)) + vloc (igtongl (ng), nt) * strf (ng, nt)
enddo
enddo
if (gamma_only) then
do ng = 1, ngm
aux (nlm(ng)) = conjg (aux (nl(ng)))
enddo
end if
!
! aux = potential in G-space . FFT to real space
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
do ir = 1, nrxx
vltot (ir) = DREAL (aux (ir) )
enddo
!
! If required add an electric field to the local potential
!
if (tefield.and.(.not.dipfield)) then
call add_efield(vltot)
endif
deallocate(aux)
return
end subroutine setlocal

306
Gamma/v_of_rho.f90
View File

@ -1,306 +0,0 @@
!
! Copyright (C) 2001-2003 PWSCF group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!--------------------------------------------------------------------
subroutine v_of_rho (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
nrx3, nrxx, nl, ngm, gstart, nspin, g, gg, alat, omega, &
ehart, etxc, vtxc, charge, v)
!--------------------------------------------------------------------
!
! This routine computes the Hartree and Exchange and Correlation
! potential and energies which corresponds to a given charge density
! The XC potential is computed in real space, while the
! Hartree potential is computed in reciprocal space.
!
!
use parameters, only: DP
implicit none
!
! first the dummy variables
!
integer :: nspin, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, nl (ngm),&
gstart
! input: 1=lda, 2=lsda
! input: the FFT indices
! input: the true array dimensions
! input: the total dimension
! input: the number of G vectors
! input: correspondence G <-> FFT
! input: first nonzero G-vector
real(kind=DP) :: rho (nrxx, nspin), rho_core (nrxx), g (3, ngm), &
gg (ngm), alat, omega, vtxc, etxc, ehart, charge, v (nrxx, nspin)
! input: the valence charge
! input: the core charge
! input: the G vectors
! input: the norm of G vector
! input: the length of the cell
! input: the volume of the cell
! output: the integral V_xc * rho
! output: the E_xc energy
! output: the hartree energy
! output: the integral of the charge
! output: the H+xc_up potential
!
integer :: is
!
call start_clock ('v_of_rho')
!
! calculate exchange-correlation potential
!
call v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
nl, ngm, g, nspin, alat, omega, etxc, vtxc, v)
!
! calculate hartree potential
!
call v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
nl, ngm, gg, gstart, nspin, alat, omega, ehart, charge, v)
!
do is=1,nspin
call add_efield(v(1,is))
enddo
call stop_clock ('v_of_rho')
return
end subroutine v_of_rho
!
!--------------------------------------------------------------------
subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, nl, ngm, g, nspin, alat, omega, etxc, vtxc, v)
!--------------------------------------------------------------------
!
! Exchange-Correlation potential Vxc(r) from n(r)
!
use parameters, only : DP
implicit none
!
! input
!
integer :: nspin, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, &
nl (ngm)
! nspin=1 :unpolarized, =2 :spin-polarized
! the FFT indices
! the true FFT array dimensions
! the total dimension
! the number of G vectors
! correspondence G <-> FFT
real(kind=DP) :: rho (nrxx, nspin), rho_core (nrxx), g (3, ngm), &
alat, omega
! the valence charge
! the core charge
! the G vectors
! the length of the cell
! the volume of the cell
!
! output
!
real(kind=DP) :: v (nrxx, nspin), vtxc, etxc
! V_xc potential
! integral V_xc * rho
! E_xc energy
!
! local variables
!
! the square of the e charge
real(kind=DP), parameter :: e2 = 2.d0
real(kind=DP) :: rhox, arhox, zeta, ex, ec, vx (2), vc (2)
! the total charge in each point
! the absolute value of the charge
! the absolute value of the charge
! local exchange energy
! local correlation energy
! local exchange potential
! local correlation potential
integer :: ir, is, ig, neg (3)
! counter on mesh points
! counter on spin polarizations
! counter on G vectors
! number of points with wrong zeta/charge
!
!
! call start_clock('vxc')
!
! initialization
!
etxc = 0.d0
vtxc = 0.d0
v(:,:) = 0.d0
if (nspin == 1) then
!
! spin-unpolarized case
!
do ir = 1, nrxx
rhox = rho (ir, nspin) + rho_core (ir)
arhox = abs (rhox)
if (arhox.gt.1.d-30) then
call xc (arhox, ex, ec, vx, vc)
v(ir,nspin) = e2 * (vx(1) + vc(1) )
etxc = etxc + e2 * (ex + ec) * rhox
vtxc = vtxc + v(ir,nspin) * rho(ir,nspin)
endif
enddo
else
!
! spin-polarized case
!
neg (1) = 0
neg (2) = 0
neg (3) = 0
do ir = 1, nrxx
rhox = rho(ir,1) + rho(ir,2) + rho_core(ir)
arhox = abs(rhox)
if (arhox.gt.1.d-30) then
zeta = ( rho(ir,1) - rho(ir,2) ) / arhox
if (abs(zeta) .gt.1.d0) then
neg(3) = neg(3) + 1
zeta = sign(1.d0,zeta)
endif
if (rho(ir,1) < 0.d0) neg(1) = neg(1) + 1
if (rho(ir,2) < 0.d0) neg(2) = neg(2) + 1
call xc_spin (arhox, zeta, ex, ec, vx(1), vx(2), vc(1), vc(2) )
do is = 1, nspin
v(ir,is) = e2 * (vx(is) + vc(is) )
enddo
etxc = etxc + e2 * (ex + ec) * rhox
vtxc = vtxc + v(ir,1) * rho(ir,1) + v(ir,2) * rho(ir,2)
endif
enddo
#ifdef __PARA
call ireduce (3, neg)
#endif
if (neg(3).gt.0) write (6,'(/,4x," npt with |zeta| > 1: ",i8, &
& ", npt tot ",i8, ",",f10.2, " %" )') neg(3), &
& nr1*nr2*nr3, float(neg(3)*100) / real(nr1*nr2*nr3)
if (neg(1).gt.0) write (6,'(/,4x," npt with rhoup < 0: ",i8, &
& ", npt tot ",i8, ",",f10.2, " %" )') neg(1), &
& nr1*nr2*nr3, float(neg(1)*100) / real(nr1*nr2*nr3)
if (neg(2).gt.0) write (6,'(/,4x," npt with rhodw < 0: ",i8, &
& ", npt tot ",i8, ",",f10.2, " %" )') neg(2), &
& nr1*nr2*nr3, float(neg(2)*100) / real(nr1 * nr2 * nr3)
endif
!
! energy terms, local-density contribution
!
vtxc = omega * vtxc / (nr1 * nr2 * nr3)
etxc = omega * etxc / (nr1 * nr2 * nr3)
!
! add gradient corrections (if any)
!
call gradcorr (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, nl, ngm, g, alat, omega, nspin, etxc, vtxc, v)
#ifdef __PARA
call reduce (1, vtxc)
call reduce (1, etxc)
#endif
! call stop_clock('vxc')
!
return
end subroutine v_xc
!
!--------------------------------------------------------------------
subroutine v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
nrxx, nl, ngm, gg, gstart, nspin, alat, omega, ehart, charge, v)
!--------------------------------------------------------------------
!
! Hartree potential VH(r) from n(r)
!
USE parameters, ONLY: DP
USE gvect, ONLY: nlm
USE wvfct, ONLY: gamma_only
implicit none
!
! input
!
integer :: nspin, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, ngm, &
gstart, nl (ngm)
real(kind=DP) :: rho (nrxx, nspin), gg (ngm), alat, omega
!
! output
!
real(kind=DP) :: v (nrxx, nspin), ehart, charge
!
! local variables
!
real(kind=DP), parameter :: fpi = 4.d0 * 3.14159265358979d0, &
e2 = 2.d0
real(kind=DP) :: tpiba2, fac
real(kind=DP), allocatable :: aux (:,:), aux1 (:,:)
integer :: ir, is, ig
!
! call start_clock('vh')
!
allocate (aux(2,nrxx), aux1(2,ngm) )
tpiba2 = (fpi / 2.d0 / alat) **2
!
! copy total rho in aux
!
aux(2,:) = 0.d0
aux(1,:) = rho(:,1)
if (nspin == 2) aux(1,:) = aux(1,:) + rho(:,2)
!
! bring rho (aux) to G space
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
charge = 0.d0
if (gstart == 2) charge = omega * aux(1,nl(1))
#ifdef __PARA
call reduce (1, charge)
#endif
!
! calculate hartree potential in G-space (NB: V(G=0)=0 )
!
ehart = 0.d0
aux1(:,:) = 0.d0
do ig = gstart, ngm
fac = e2 * fpi / (tpiba2 * gg (ig) )
ehart = ehart + (aux(1,nl(ig))**2 + aux(2,nl(ig))**2) * fac
aux1(1,ig) = fac * aux(1,nl(ig))
aux1(2,ig) = fac * aux(2,nl(ig))
enddo
if (gamma_only) then
ehart = ehart * omega
else
ehart = ehart * omega / 2.d0
end if
#ifdef __PARA
call reduce (1, ehart)
#endif
aux(:,:) = 0.d0
do ig = 1, ngm
aux(1,nl(ig)) = aux1(1,ig)
aux(2,nl(ig)) = aux1(2,ig)
enddo
if (gamma_only) then
do ig = 1, ngm
aux(1,nlm(ig)) = aux1(1,ig)
aux(2,nlm(ig)) = - aux1(2,ig)
enddo
end if
!
! transform hartree potential to real space
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! add hartree potential to the xc potential
!
do is = 1, nspin
do ir = 1, nrxx
v(ir,is) = v(ir,is) + aux(1,ir)
enddo
enddo
deallocate (aux,aux1)
!
! call stop_clock('vh')
!
return
end subroutine v_h

1
PH/phonon.f90
View File

@ -34,6 +34,7 @@ program phonon
call init_clocks (.true.)
call start_clock ('PHONON')
gamma_only = .false.
call startup (nd_nmbr, code, version_number)
write (6, '(/5x,"Ultrasoft (Vanderbilt) Pseudopotentials")')
!

8
PP/start_postproc.f90
View File

@ -11,11 +11,17 @@ subroutine start_postproc (nodenumber)
!
! Wrapper routine for postprocessing initialization
!
use global_version
USE global_version, ONLY: version_number
USE wvfct, ONLY: gamma_only
implicit none
character(len=3) :: nodenumber
character(len=9) :: code = 'POST-PROC'
!
! presently no postprocessing code is expected to work with
! half G-vector sphere
!
gamma_only = .FALSE.
!
call startup (nodenumber, code, version_number)
#ifdef __PARA
call init_pool

3
PW/addusdens.f90
View File

@ -93,8 +93,9 @@ subroutine addusdens
do is = 1, nspin
psic(:) = (0.d0, 0.d0)
psic( nl(:) ) = aux(:,is)
if (gamma_only) psic( nlm(:) ) = conjg(aux(:,is))
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
call DAXPY (nrxx, 1.d0, psic, 2, rho(1,is), 1)
rho (:, is) = rho (:, is) + DREAL (psic (:) )
enddo
deallocate (aux)

17
PW/allocate_fft.f90
View File

@ -15,12 +15,12 @@ subroutine allocate_fft
!
#include "machine.h"
use pwcom
USE wavefunctions, ONLY : evc, psic
USE wavefunctions, ONLY : psic
implicit none
!
! determines the data structure for fft arrays
!
call data_structure( .FALSE. )
call data_structure( gamma_only )
!
if (nrxx.lt.ngm) then
write (6, '(/,4x," nr1=",i4," nr2= ", i4, " nr3=",i4, &
@ -34,17 +34,18 @@ subroutine allocate_fft
call errore ('allocate_fft', 'the nrs"s are too small!', 1)
endif
if (ngm.le.0) call errore ('allocate_fft', 'wrong ngm', 1)
if (ngms.le.0) call errore ('allocate_fft', 'wrong ngms', 1)
if (nrxx.le.0) call errore ('allocate_fft', 'wrong nrxx', 1)
if (nrxxs.le.0) call errore ('allocate_fft', 'wrong nrxxs', 1)
if (nspin.le.0) call errore ('allocate_fft', 'wrong nspin', 1)
if (ngm <= 0) call errore ('allocate_fft', 'wrong ngm', 1)
if (ngms <= 0) call errore ('allocate_fft', 'wrong ngms', 1)
if (nrxx <= 0) call errore ('allocate_fft', 'wrong nrxx', 1)
if (nrxxs<= 0) call errore ('allocate_fft', 'wrong nrxxs', 1)
if (nspin<= 0) call errore ('allocate_fft', 'wrong nspin', 1)
!
! Allocate memory for all kind of stuff.
!
allocate (g( 3, ngm))
allocate (gg( ngm))
allocate (nl( ngm))
if (gamma_only) allocate (nlm(ngm))
allocate (igtongl( ngm))
allocate (ig1( ngm))
allocate (ig2( ngm))
@ -59,8 +60,10 @@ subroutine allocate_fft
allocate (vrs( nrxx, nspin))
if (doublegrid) then
allocate (nls( ngms))
if (gamma_only) allocate (nlsm(ngm))
else
nls => nl
if (gamma_only) nlsm=> nlm
endif
return
end subroutine allocate_fft

57
PW/atomic_rho.f90
View File

@ -17,7 +17,7 @@ subroutine atomic_rho (rhoa, nspina)
! calculated assuming an uniform atomic spin-polarization
! equal to starting_magnetization(nt)
!
! NB: nspina may not be equal to nspin because in some cases (as in upda
! NB: nspina may not be equal to nspin because in some cases (as in update)
! the total charge only could be needed, even in a LSDA calculation.
!
!
@ -28,37 +28,24 @@ subroutine atomic_rho (rhoa, nspina)
implicit none
integer :: nspina
! the number of spin polarizations
real(kind=DP) :: rhoa (nrxx, nspina), rhoneg, rhorea, rhoima, gx
real(kind=DP), allocatable :: rhocgnt (:), aux (:)
real(kind=DP) :: rhoa (nrxx, nspina)
! the output atomic charge
! negative charge
! real charge
! imaginary charge
! the modulus of G
! the value of the integral
! the integrand function
complex(kind=DP), allocatable :: rhocg (:,:)
! auxiliary var: charge dens. in G space
integer :: ir, is, ig, igl, igl0, nt
! counter on mesh points
! counter on spin polarizations
! counter on G vectors
! counter on G vectors shells
! index of first shell with G != 0
! counter on atom types
!
! superposition of atomic charges contained in the array rho_at and
! already set in readin-readvan
! local variables
!
real(kind=DP) :: rhoneg, rhorea, rhoima, gx
real(kind=DP), allocatable :: rhocgnt (:), aux (:)
complex(kind=DP), allocatable :: rhocg (:,:)
integer :: ir, is, ig, igl, nt
!
! superposition of atomic charges contained in the array rho_at
! (read from pseudopotential files)
!
! allocate work space (psic must already be allocated)
!
allocate (rhocg( ngm, nspina))
allocate (aux( ndm))
allocate (rhocgnt( ngl))
! psic is the generic work space
rhoa(:,:) = 0.d0
rhocg(:,:) = (0.d0,0.d0)
@ -66,19 +53,16 @@ subroutine atomic_rho (rhoa, nspina)
!
! Here we compute the G=0 term
!
if (gl (1) .lt.1.0d-8) then
if (gstart == 2) then
do ir = 1, msh (nt)
aux (ir) = rho_at (ir, nt)
enddo
call simpson (msh (nt), aux, rab (1, nt), rhocgnt (1) )
igl0 = 2
else
igl0 = 1
endif
!
! Here we compute the G<>0 term
!
do igl = igl0, ngl
do igl = gstart, ngl
gx = sqrt (gl (igl) ) * tpiba
do ir = 1, msh (nt)
if (r (ir, nt) .lt.1.0d-8) then
@ -92,7 +76,7 @@ subroutine atomic_rho (rhoa, nspina)
!
! we compute the 3D atomic charge in reciprocal space
!
if (nspina.eq.1) then
if (nspina == 1) then
do ig = 1, ngm
rhocg(ig,1) = rhocg(ig,1) + &
strf(ig,nt) * rhocgnt(igtongl(ig)) / omega
@ -109,16 +93,15 @@ subroutine atomic_rho (rhoa, nspina)
endif
enddo
deallocate (rhocgnt)
deallocate (aux)
do is = 1, nspina
!
! and we return to real space
!
psic(:) = (0.d0,0.d0)
do ig = 1, ngm
psic (nl (ig) ) = rhocg (ig, is)
enddo
psic (nl (:) ) = rhocg (:, is)
if (gamma_only) psic ( nlm(:) ) = conjg( rhocg (:, is) )
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
!
! we check that everything is correct
@ -137,7 +120,7 @@ subroutine atomic_rho (rhoa, nspina)
call reduce (1, rhoneg)
call reduce (1, rhoima)
#endif
if ( rhoneg.lt.-1.0d-4 .or. rhoima.gt.1.0d-4 ) &
if ( rhoneg < -1.0d-4 .or. rhoima > 1.0d-4 ) &
write (6,'(/" Warning: negative or imaginary starting charge ",&
&2f12.6,i3)') rhoneg, rhoima, is
enddo

20
PW/gradcorr.f90
View File

@ -26,7 +26,7 @@ subroutine gradcorr (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
v2xup, v2xdw, v1cup, v1cdw , etxcgc, vtxcgc, segno, arho, fac
real(kind=DP), parameter :: e2 = 2.d0, epsr = 1.0d-6, epsg = 1.0d-10
if (igcx.eq.0.and.igcc.eq.0) return
if (igcx == 0 .and. igcc == 0) return
etxcgc = 0.d0
vtxcgc = 0.d0
@ -44,7 +44,7 @@ subroutine gradcorr (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
do is = 1, nspin
grho2 (is) = grho(1, k, is)**2 + grho(2, k, is)**2 + grho(3, k, is)**2
enddo
if (nspin.eq.1) then
if (nspin == 1) then
!
! This is the spin-unpolarised case
!
@ -142,6 +142,8 @@ subroutine gradient (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
!
! Calculates ga = \grad a in R-space (a is also in R-space)
use parameters
use gvect, only: nlm
use wvfct, only: gamma_only
implicit none
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, ngm, nl (ngm)
@ -175,6 +177,12 @@ subroutine gradient (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
gaux (1, nl (n) ) = - g (ipol, n) * aux (2, nl (n) )
gaux (2, nl (n) ) = g (ipol, n) * aux (1, nl (n) )
enddo
if (gamma_only) then
do n = 1, ngm
gaux (1, nlm(n) ) = gaux (1, nl(n) )
gaux (2, nlm(n) ) = - gaux (2, nl(n) )
enddo
end if
!
! bring back to R-space, (\grad_ipol a)(r) ...
!
@ -199,6 +207,8 @@ subroutine grad_dot (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
!
! Calculates da = \sum_i \grad_i a_i in R-space
use parameters
use gvect, only: nlm
use wvfct, only: gamma_only
implicit none
integer :: nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, ngm, nl (ngm)
real(kind=DP) :: a (3, nrxx), g (3, ngm), da (nrxx), alat
@ -229,6 +239,12 @@ subroutine grad_dot (nrx1, nrx2, nrx3, nr1, nr2, nr3, nrxx, a, &
gaux (2, nl (n) ) = gaux (2, nl (n) ) + g (ipol, n) * aux (1,nl(n))
enddo
enddo
if (gamma_only) then
do n = 1, ngm
gaux (1, nlm(n) ) = gaux (1, nl (n) )
gaux (2, nlm(n) ) = - gaux (2, nl (n) )
enddo
end if
!
! bring back to R-space, (\grad_ipol a)(r) ...
!

29
PW/interpolate.f90
View File

@ -16,7 +16,11 @@ subroutine interpolate (v, vs, iflag)
! V and Vs are real and in real space . V and Vs may coincide
!
#include"machine.h"
use pwcom
USE parameters, ONLY: DP
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm
USE gsmooth,ONLY: nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,nrxxs,ngms, &
nls, nlsm, doublegrid
implicit none
real(kind=DP) :: v (nrxx), vs (nrxxs)
! function on thick mesh
@ -27,12 +31,12 @@ subroutine interpolate (v, vs, iflag)
! work array on smooth mesh
integer :: iflag
! gives the direction of the interpolat
! gives the direction of the interpolation
integer :: ig, ir
call start_clock ('interpolate')
if (iflag.le.0) then
if (iflag <= 0) then
!
! from thick to smooth
!
@ -45,6 +49,11 @@ subroutine interpolate (v, vs, iflag)
do ig = 1, ngms
auxs (nls (ig) ) = aux (nl (ig) )
enddo
if (gamma_only) then
do ig = 1, ngms
auxs (nlsm(ig) ) = aux (nlm(ig) )
enddo
end if
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 1)
vs (:) = auxs (:)
deallocate (auxs)
@ -67,6 +76,11 @@ subroutine interpolate (v, vs, iflag)
do ig = 1, ngms
aux (nl (ig) ) = auxs (nls (ig) )
enddo
if (gamma_only) then
do ig = 1, ngms
aux (nlm(ig) ) = aux (nlsm(ig) )
enddo
end if
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
v (:) = aux (:)
deallocate (auxs)
@ -91,7 +105,11 @@ subroutine cinterpolate (v, vs, iflag)
! V and Vs are complex and in real space . V and Vs may coincide
!
#include"machine.h"
use pwcom
USE parameters, ONLY: DP
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm
USE gsmooth,ONLY: nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,nrxxs,ngms, &
nls, nlsm, doublegrid
complex(kind=DP) :: v (nrxx), vs (nrxxs)
! function on thick mesh
! function on smooth mesh
@ -105,8 +123,9 @@ subroutine cinterpolate (v, vs, iflag)
integer :: ig
if (gamma_only) call errore ('cinterpolate','not implemented', 1)
call start_clock ('cinterpolate')
if (iflag.le.0) then
if (iflag <= 0) then
!
! from thick to smooth
!

52
PW/mix_rho.f90
View File

@ -194,6 +194,11 @@ subroutine mix_rho (rhout, rhoin, nsout, nsin, alphamix, dr2, iter, &
do ig=1,ngm0
psic(nl(ig)) = rhocin(ig,is)+rhocout(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg ( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
call DAXPY(nrxx,-1.d0,psic,2,rhoin(1,is),1)
end do
@ -344,6 +349,11 @@ subroutine mix_rho (rhout, rhoin, nsout, nsin, alphamix, dr2, iter, &
do ig=1,ngm0
psic(nl(ig)) = rhocin(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg ( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
call DAXPY(nrxx,1.d0,psic,2,rhoin(1,is),1)
end do
@ -388,6 +398,7 @@ function rho_dot_product (rho1,rho2)
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(ig,is))*rho2(ig,is))
end do
if (gamma_only) rho_dot_product = 2.d0*rho_dot_product
else
do ig = gstart,ngm0
fac = e2*fpi / (tpiba2*gg(ig))
@ -395,9 +406,17 @@ function rho_dot_product (rho1,rho2)
DREAL(conjg(rho1(ig,1)+rho1(ig,2))* &
(rho2(ig,1)+rho2(ig,2)))
end do
if (gamma_only) rho_dot_product = 2.d0*rho_dot_product
fac = e2*fpi / (tpi**2) ! lambda=1 a.u.
do ig = 1,ngm0
rho_dot_product = rho_dot_product + fac * &
! G=0 term
if (gstart == 2) then
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(1,1)-rho1(1,2))* &
(rho2(1,1)-rho2(1,2)))
end if
if (gamma_only) fac = 2.d0*fac
do ig = gstart,ngm0
rho_dot_product = rho_dot_product + fac * &
DREAL(conjg(rho1(ig,1)-rho1(ig,2))* &
(rho2(ig,1)-rho2(ig,2)))
end do
@ -490,8 +509,11 @@ function fn_dehar (drho)
end do
end if
fn_dehar = fn_dehar * omega / 2.d0
if (gamma_only) then
fn_dehar = fn_dehar * omega
else
fn_dehar = fn_dehar * omega / 2.d0
end if
#ifdef __PARA
call reduce(1,fn_dehar)
#endif
@ -619,6 +641,12 @@ end subroutine approx_screening
* exp(-0.5*l2smooth*tpiba2*gg(ig))
end do
end if
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
alpha(:) = real(psic(:))
@ -661,6 +689,12 @@ end subroutine approx_screening
do ig=1,ngm0
psic(nl(ig)) = drho(ig,is)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir) * alpha(ir)
@ -686,6 +720,11 @@ end subroutine approx_screening
do ig=1,ngm0
psic(nl(ig)) = v(ig,m)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir)*fpi*e2/alpha(ir)
@ -747,6 +786,11 @@ end subroutine approx_screening
do ig=1,ngm0
psic(nl(ig)) = vbest(ig)
end do
if (gamma_only) then
do ig=1,ngm0
psic(nlm(ig)) = conjg( psic(nl(ig)) )
end do
end if
call cft3(psic,nr1,nr2,nr3,nrx1,nrx2,nrx3,+1)
do ir=1,nrxx
psic(ir) = psic(ir)/alpha(ir)

5
PW/set_rhoc.f90
View File

@ -67,6 +67,11 @@ subroutine set_rhoc
enddo
endif
enddo
if (gamma_only) then
do ng = 1, ngm
aux(nlm(ng)) = conjg(aux(nl (ng)))
end do
end if
!
! the core charge in real space
!

18
PW/setlocal.f90
View File

@ -13,8 +13,14 @@ subroutine setlocal
! This routine computes the local potential in real space vltot(ir)
!
#include "machine.h"
use pwcom
USE parameters, ONLY: DP
USE basis, ONLY: ntyp
USE extfield, ONLY: tefield, dipfield
USE gvect, ONLY : nl, nlm, igtongl
USE scf, ONLY: vltot
USE vlocal, ONLY : strf, vloc
USE wvfct, ONLY: gamma_only
USE gvect, ONLY: nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, nlm, ngm
implicit none
complex(kind=DP), allocatable :: aux (:)
! auxiliary variable
@ -23,8 +29,7 @@ subroutine setlocal
! counter on g vectors
! counter on r vectors
!
allocate (aux( nrxx))
allocate (aux( nrxx))
!
aux(:)=(0.d0,0.d0)
do nt = 1, ntyp
@ -32,6 +37,11 @@ subroutine setlocal
aux (nl(ng))=aux(nl(ng)) + vloc (igtongl (ng), nt) * strf (ng, nt)
enddo
enddo
if (gamma_only) then
do ng = 1, ngm
aux (nlm(ng)) = conjg (aux (nl(ng)))
enddo
end if
!
! aux = potential in G-space . FFT to real space
!

22
PW/v_of_rho.f90
View File

@ -132,7 +132,7 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
vtxc = 0.d0
v(:,:) = 0.d0
if (nspin.eq.1) then
if (nspin == 1) then
!
! spin-unpolarized case
!
@ -162,8 +162,8 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
neg(3) = neg(3) + 1
zeta = sign(1.d0,zeta)
endif
if (rho(ir,1) .lt.0.d0) neg(1) = neg(1) + 1
if (rho(ir,2) .lt.0.d0) neg(2) = neg(2) + 1
if (rho(ir,1) < 0.d0) neg(1) = neg(1) + 1
if (rho(ir,2) < 0.d0) neg(2) = neg(2) + 1
call xc_spin (arhox, zeta, ex, ec, vx(1), vx(2), vc(1), vc(2) )
do is = 1, nspin
v(ir,is) = e2 * (vx(is) + vc(is) )
@ -212,7 +212,9 @@ subroutine v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
!
! Hartree potential VH(r) from n(r)
!
use parameters, only: DP
USE parameters, ONLY: DP
USE gvect, ONLY: nlm
USE wvfct, ONLY: gamma_only
implicit none
!
! input
@ -263,7 +265,11 @@ subroutine v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
aux1(1,ig) = fac * aux(1,nl(ig))
aux1(2,ig) = fac * aux(2,nl(ig))
enddo
ehart = ehart * omega / 2.d0
if (gamma_only) then
ehart = ehart * omega
else
ehart = ehart * omega / 2.d0
end if
#ifdef __PARA
call reduce (1, ehart)
#endif
@ -272,6 +278,12 @@ subroutine v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
aux(1,nl(ig)) = aux1(1,ig)
aux(2,nl(ig)) = aux1(2,ig)
enddo
if (gamma_only) then
do ig = 1, ngm
aux(1,nlm(ig)) = aux1(1,ig)
aux(2,nlm(ig)) = - aux1(2,ig)
enddo
end if
!
! transform hartree potential to real space
!