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

sp3 fix (allocation with zero length). Manual updates. 

Gamma: fixed occupations, electric fields, cleanup.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@289 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2003-08-21 15:14:22 +00:00
parent 847cae092d
commit 319ae12049
12 changed files with 162 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
Gamma/Makefile
View File

@ -14,6 +14,7 @@ PWOBJS=../PW/pwcom.o \
../PW/aainit.o \
../PW/addusforce.o \
../PW/addusstress.o \
../PW/add_efield.o \
../PW/allocate_locpot.o \
../PW/allocate_nlpot.o \
../PW/atomic_wfc.o \
@ -29,6 +30,7 @@ PWOBJS=../PW/pwcom.o \
../PW/checkallsym.o \
../PW/checksym.o \
../PW/clean_pw.o \
../PW/compute_dip.o \
../PW/constrain.o \
../PW/conv_to_num.o \
../PW/coset.o \
@ -64,6 +66,7 @@ PWOBJS=../PW/pwcom.o \
../PW/g_psi.o \
../PW/gen_us_dj.o \
../PW/gen_us_dy.o \
../PW/ggen.o \
../PW/gk_sort.o \
../PW/gweights.o \
../PW/hexsym.o \
@ -164,6 +167,7 @@ PWOBJS=../PW/pwcom.o \
../PW/w0gauss.o \
../PW/w1gauss.o \
../PW/wgauss.o \
../PW/wsweight.o \
../PW/which_dft.o \
../PW/write_ns.o \
../PW/ylmr2.o \
@ -194,9 +198,8 @@ PWOBJS=../PW/pwcom.o \
../PW/bp_ylm_q.o \
../PW/bp_zgedi.o \
../PW/bp_zgefa.o \
../PW/ggen.o \
../PW/ccalbec.o \
../PW/becmod.o
../PW/becmod.o
MODULES = ../Modules/*.o ../PW/error_handler.o \

9
Gamma/atomic_rho.f90
View File

@ -40,7 +40,7 @@ subroutine atomic_rho (rhoa, nspina)
! the integrand function
complex(kind=DP), allocatable :: rhocg (:,:)
! auxiliary var: charge dens. in G spac
! auxiliary var: charge dens. in G space
integer :: ir, is, ig, igl, igl0, nt
! counter on mesh points
@ -58,10 +58,9 @@ subroutine atomic_rho (rhoa, nspina)
allocate (aux( ndm))
allocate (rhocgnt( ngl))
! psic is the generic work space
call setv (nrxx, 0.d0, rhoa, 1)
rhoa(:,:) = 0.d0
rhocg(:,:) = (0.d0,0.d0)
call setv (2 * nspina * ngm, 0.d0, rhocg, 1)
do nt = 1, ntyp
!
! Here we compute the G=0 term
@ -118,7 +117,7 @@ subroutine atomic_rho (rhoa, nspina)
!
! and we return to real space
!
call setv (2 * nrxx, 0.d0, psic, 1)
psic(:) = (0.d0, 0.d0)
do ig = 1, ngm
psic (nl (ig) ) = rhocg (ig, is)
psic (nlm(ig) ) = conjg( rhocg (ig, is) )

32
Gamma/interpolate.f90
View File

@ -23,12 +23,12 @@ subroutine interpolate (v, vs, iflag)
! function on thick mesh
! function on smooth mesh
complex(kind=DP),pointer :: aux (:), auxs (:)
complex(kind=DP), allocatable :: aux (:), auxs (:)
! work array on thick mesh
! work array on smooth mesh
integer :: iflag
! gives the direction of the interpolat
! gives the direction of the interpolation
integer :: ig, ir
@ -40,23 +40,21 @@ subroutine interpolate (v, vs, iflag)
if (doublegrid) then
allocate (aux( nrxx))
allocate (auxs(nrxxs))
do ir = 1, nrxx
aux (ir) = v (ir)
enddo
aux (:) = v (:)
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
call setv (2 * nrxxs, 0.d0, auxs, 1)
auxs (:) = (0.d0, 0.d0)
do ig = 1, ngms
auxs (nls (ig) ) = aux (nl (ig) )
auxs (nlsm(ig) ) = aux (nlm(ig) )
enddo
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 1)
do ir = 1, nrxxs
vs (ir) = auxs (ir)
enddo
vs (:) = auxs (:)
deallocate (auxs)
deallocate (aux)
else
call DCOPY (nrxx, v, 1, vs, 1)
do ir = 1, nrxx
vs (ir) = v (ir)
enddo
endif
else
!
@ -65,23 +63,21 @@ subroutine interpolate (v, vs, iflag)
if (doublegrid) then
allocate (aux( nrxx))
allocate (auxs(nrxxs))
do ir = 1, nrxxs
auxs (ir) = vs (ir)
enddo
auxs (:) = vs (:)
call cft3s (auxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 1)
call setv (2 * nrxx, 0.d0, aux, 1)
aux (:) = (0.d0, 0.d0)
do ig = 1, ngms
aux (nl (ig) ) = auxs (nls (ig) )
aux (nlm(ig) ) = auxs (nlsm(ig) )
enddo
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do ir = 1, nrxx
v (ir) = aux (ir)
enddo
v (:) = aux (:)
deallocate (auxs)
deallocate (aux)
else
call DCOPY (nrxx, vs, 1, v, 1)
do ir = 1, nrxx
v (ir) = vs (ir)
enddo
endif
endif
call stop_clock ('interpolate')

6
Gamma/set_rhoc.f90
View File

@ -41,7 +41,7 @@ subroutine set_rhoc
do nt = 1, ntyp
if (nlcc (nt) ) goto 10
enddo
call setv (nrxx, 0.d0, rho_core, 1)
rho_core(:) = 0.d0
return
10 continue
@ -88,10 +88,10 @@ subroutine set_rhoc
! 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 i
! 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
! rather large and it is better to leave the core charge as it is.
! 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
!

9
Gamma/setlocal.f90
View File

@ -32,7 +32,7 @@ subroutine setlocal
!
! here we compute the local potential in real space
!
call setv (2 * nrxx, 0.d0, aux, 1)
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)
@ -49,6 +49,13 @@ subroutine setlocal
vltot (ir) = DREAL (aux (ir) )
enddo
deallocate(aux)
!
! If required add an electric field to the local potential
!
if (tefield.and.(.not.dipfield)) then
call add_efield(vltot)
endif
!
return
end subroutine setlocal

14
Gamma/sum_band.f90
View File

@ -36,14 +36,14 @@ subroutine sum_band
! weight
!
call start_clock ('sum_band')
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum, 1)
call setv (nrxx * nspin, 0.0d0, rho, 1)
becsum(:,:,:) = 0.d0
rho(:,:) = 0.d0
eband = 0.d0
demet = 0.d0
!
! calculate weights for the insulator case
!
if (.not.lgauss.and..not.ltetra) then
if (.not.lgauss.and..not.ltetra.and..not.tfixed_occ) then
call iweights (nks, wk, nbnd, nelec, wg)
!
! calculate weights for the metallic case
@ -64,6 +64,12 @@ subroutine sum_band
elseif (lgauss) then
call gweights (nks, wk, nbnd, nelec, degauss, ngauss, &
et, ef, demet, wg)
elseif (tfixed_occ) then
do is=1,nspin
do ibnd=1,nbnd
wg(ibnd,is)=f_inp(ibnd,is)
enddo
enddo
endif
!
! Needed for LDA+U
@ -94,7 +100,7 @@ subroutine sum_band
!
end do
do ibnd = 1, nbnd, 2
call setv (2 * nrxx, 0.d0, psic, 1)
psic(:) = (0.d0,0.d0)
if (ibnd.lt.nbnd) then
! two ffts at the same time
do ig = 1, npw

182
Gamma/v_of_rho.f90
View File

@ -18,9 +18,7 @@ subroutine v_of_rho (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
! Hartree potential is computed in reciprocal space.
!
!
#include "machine.h"
use parameters
use gamma, only: nlm
use parameters, only: DP
implicit none
!
! first the dummy variables
@ -50,24 +48,8 @@ subroutine v_of_rho (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
! output: the integral of the charge
! output: the H+xc_up potential
!
real(kind=DP), parameter :: fpi = 4.d0 * 3.14159265358979d0, &
e2 = 2.d0
integer :: is
!
! and the local variables
!
real(kind=DP) :: tpiba2, fac
! the measure unit in reciprocal space
! a multiplicative factors
real(kind=DP), allocatable :: aux (:,:), aux1 (:,:)
! used to do the fft
! auxiliary variable for the potential
integer :: ir, is, ig
! counter on mesh points
! counter on spin polarizations
! counter on G vectors
call start_clock ('v_of_rho')
!
! calculate exchange-correlation potential
@ -75,61 +57,15 @@ subroutine v_of_rho (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
call v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
nl, ngm, g, nspin, alat, omega, etxc, vtxc, v)
!
allocate (aux(2,nrxx),aux1(2,ngm) )
tpiba2 = (fpi / 2.d0 / alat) **2
! calculate hartree potential
!
! copy total rho in aux
call v_h (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, &
nl, ngm, gg, gstart, nspin, alat, omega, ehart, charge, v)
!
call setv (nrxx, 0.d0, aux (2, 1), 2)
call DCOPY (nrxx, rho (1, 1), 1, aux (1, 1), 2)
if (nspin.eq.2) call DAXPY (nrxx, 1.0d0, rho (1, 2), 1, aux (1, 1), 2)
!
! bring rho (aux) to G space
!
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
charge = 0.d0
if (gg (1) .lt.1.0d-8) 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
ehart = ehart * omega
#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)
aux (1, nlm(ig) ) = aux1 (1, ig)
aux (2, nlm(ig) ) =-aux1 (2, ig)
do is=1,nspin
call add_efield(v(1,is))
enddo
!
! 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 ('v_of_rho')
return
end subroutine v_of_rho
@ -176,8 +112,7 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
! local variables
!
! the square of the e charge
real(kind=DP) :: e2
parameter (e2 = 2.d0)
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
@ -191,8 +126,7 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
! counter on mesh points
! counter on spin polarizations
! counter on G vectors
! number of points with wrong zeta/cha
!
! number of points with wrong zeta/charge
!
! call start_clock('vxc')
!
@ -200,9 +134,9 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
!
etxc = 0.d0
vtxc = 0.d0
v(:,:) = 0.d0
call setv (nspin * nrxx, 0.d0, v, 1)
if (nspin.eq.1) then
if (nspin == 1) then
!
! spin-unpolarized case
!
@ -232,8 +166,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) )
@ -274,3 +208,93 @@ subroutine v_xc (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, nrx3, &
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 gamma, only: nlm
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
ehart = ehart * omega
#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)
aux (1, nlm(ig) ) = aux1 (1, ig)
aux (2, nlm(ig) ) =-aux1 (2, ig)
enddo
!
! 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

6
Gamma/vloc_psi.f90
View File

@ -27,8 +27,8 @@ subroutine vloc_psi(lda, n, m, psi, v, hpsi)
! the local potential V_Loc psi. First bring psi to real space
!
do ibnd = 1, m, 2
call setv (2 * nrxxs, 0.d0, psic, 1)
if (ibnd.lt.m) then
psic(:) = (0.d0, 0.d0)
if (ibnd < m) then
! two ffts at the same time
do j = 1, n
psic (nls (igk(j))) = psi(j, ibnd) + (0.0,1.d0)*psi(j, ibnd+1)
@ -54,7 +54,7 @@ subroutine vloc_psi(lda, n, m, psi, v, hpsi)
!
! addition to the total product
!
if (ibnd.lt.m) then
if (ibnd < m) then
! two ffts at the same time
do j = 1, n
fp = (psic (nls(igk(j))) + psic (nlsm(igk(j))))*0.5d0

8
PW/add_efield.f90
View File

@ -72,18 +72,18 @@
elseif (edir.eq.3) then
npoints=nr3
else
call errore('setlocal',' wrong edir',1)
call errore('add_efield',' wrong edir',1)
endif
length=alat/bmod
deltal=length/npoints
nmax =int(real(npoints,dp)*(emaxpos-eps))+1
if (nmax.lt.1.or.nmax.gt.npoints) &
call errore('setlocal','nmax out of range',1)
call errore('add_efield','nmax out of range',1)
ndesc=int(real(npoints,dp)*(eopreg-eps))+1
if (ndesc.lt.1.or.ndesc.gt.npoints) &
call errore('setlocal','ndesc out of range',1)
call errore('add_efield','ndesc out of range',1)
dip=0.d0
dipion=0.d0
@ -234,7 +234,7 @@
end do
end do
else
call errore('setlocal', 'wrong edir', 1)
call errore('add_efield', 'wrong edir', 1)
endif
first=.false.
return

4
PW/allocate_nlpot.f90
View File

@ -63,7 +63,6 @@ subroutine allocate_nlpot
endif
endif
enddo
lqx = 2*lmaxkb+1
!
! calculate the maximum number of beta functions
!
@ -88,8 +87,9 @@ subroutine allocate_nlpot
!
nqxq = ( (sqrt(gcutm) + sqrt(xqq(1)**2 + xqq(2)**2 + xqq(3)**2) ) &
/ dq + 4) * cell_factor
lqx = 2*lmaxkb+1
!
allocate (qrad( nqxq, nbrx*(nbrx+1)/2, lqx, ntyp))
if (lqx > 0) allocate (qrad( nqxq, nbrx*(nbrx+1)/2, lqx, ntyp))
allocate (vkb( npwx, nkb))
allocate (qgm( ngm))
allocate (becsum( nhm * (nhm + 1)/2, nat, nspin))

12
PW/input.f90
View File

@ -415,7 +415,7 @@ subroutine iosys
ntcheck=nstep+1
CASE DEFAULT
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': ion_dymanics='//trim(ion_dynamics)//' not supported', 1 )
& ': ion_dynamics='//trim(ion_dynamics)//' not supported', 1 )
END SELECT
endif
if ( TRIM(calculation) == 'md' ) then
@ -430,7 +430,7 @@ subroutine iosys
ntcheck=nstep+1
CASE DEFAULT
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': ion_dymanics='//trim(ion_dynamics)//' not supported', 1 )
& ': ion_dynamics='//trim(ion_dynamics)//' not supported', 1 )
END SELECT
endif
if ( TRIM(calculation) == 'vc-relax' ) then
@ -455,11 +455,11 @@ subroutine iosys
ntcheck=nstep+1
CASE DEFAULT
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': cell_dymanics='//trim(cell_dynamics)//' not supported', 1 )
& ': cell_dynamics='//trim(cell_dynamics)//' not supported', 1 )
END SELECT
if ( TRIM(ion_dynamics) .ne. 'damp' ) then
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': ion_dymanics='//trim(ion_dynamics)//' not supported', 1 )
& ': ion_dynamics='//trim(ion_dynamics)//' not supported', 1 )
end if
endif
if ( TRIM(calculation) == 'vc-md' ) then
@ -478,11 +478,11 @@ subroutine iosys
ntcheck=nstep+1
CASE DEFAULT
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': ion_dymanics='//trim(ion_dynamics)//' not supported', 1 )
& ': ion_dynamics='//trim(ion_dynamics)//' not supported', 1 )
END SELECT
if ( TRIM(ion_dynamics) .ne. 'beeman' ) then
CALL errore(' iosys ','calculation='//trim(calculation)// &
& ': ion_dymanics='//trim(ion_dynamics)//' not supported', 1 )
& ': ion_dynamics='//trim(ion_dynamics)//' not supported', 1 )
end if
endif

1
PW/v_of_rho.f90
View File

@ -18,7 +18,6 @@ subroutine v_of_rho (rho, rho_core, nr1, nr2, nr3, nrx1, nrx2, &
! Hartree potential is computed in reciprocal space.
!
!
#include "machine.h"
use parameters, only: DP
implicit none
!