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

more cleanup, mainly in PH dir 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@194 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
degironc 2003-04-29 11:20:28 +00:00
parent 236a0a42f3
commit 92070bf269
21 changed files with 395 additions and 516 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
PH/addusddens.f90
View File

@ -29,8 +29,8 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
! input: if zero does not compute drho
! input: the number of perturbations
complex(kind=DP) :: drhoscf (nrxx, nspin, npe), dbecsum (nhm * &
(nhm + 1) / 2, nat, nspin, npe)
complex(kind=DP) :: drhoscf (nrxx, nspin, npe), &
dbecsum (nhm*(nhm+1)/2, nat, nspin, npe)
! inp/out: change of the charge density
!input: sum over kv of bec
integer :: irr, mode0
@ -53,22 +53,14 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
! counter on spin
! counter on combined beta functions
real(kind=DP), allocatable :: qmod (:),&
qpg (:,:),&
ylmk0 (:,:)
real(kind=DP), allocatable :: qmod (:), qpg (:,:), ylmk0 (:,:)
! the modulus of q+G
! the values of q+G
! the spherical harmonics
complex(kind=DP) :: fact, zsum, bb, alpha, alpha_0, u1, u2, u3
complex(kind=DP), allocatable :: sk (:),&
qg (:),&
drhous (:,:),&
aux (:,:,:)
! auxiliary variables
! auxiliary variable
! auxiliary variable
! auxiliary variables
complex(kind=DP), allocatable :: sk (:), drhous (:,:), aux (:,:,:)
! the structure factor
! auxiliary variable for FFT
! contain the charge of drho
@ -78,11 +70,9 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
call start_clock ('addusddens')
allocate (aux( ngm , nspin , 3))
allocate (sk ( ngm))
allocate (qg ( nrxx))
allocate (ylmk0( ngm , lqx * lqx))
allocate (qmod ( ngm))
if (.not.lgamma) allocate (qpg( 3 , ngm))
if (iflag.eq.0) allocate (drhous( nrxx, nspin))
! write(6,*) aux, ylmk0, qmod
!
! And then we compute the additional charge in reciprocal space
@ -98,7 +88,6 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
endif
fact = 0.5d0 * DCMPLX (0.d0, - tpiba)
call setv (6 * ngm * nspin, 0.d0, aux, 1)
@ -116,8 +105,10 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
! calculate the structure factor
!
do ig = 1, ngm
sk (ig) = eigts1 (ig1 (ig), na) * eigts2 (ig2 (ig), na) &
* eigts3 (ig3 (ig), na) * eigqts (na) * qgm (ig)
sk (ig) = eigts1 (ig1 (ig), na) * &
eigts2 (ig2 (ig), na) * &
eigts3 (ig3 (ig), na) * &
eigqts (na) * qgm (ig)
enddo
!
! And qgmq and becp and dbecq
@ -129,27 +120,25 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
u1 = u (mu + 1, mode)
u2 = u (mu + 2, mode)
u3 = u (mu + 3, mode)
if (abs (u1) + abs (u2) + abs (u3) .gt.1d-12.and.iflag.eq.1) &
then
if (abs(u1) + abs(u2) + abs(u3) .gt.1d-12 .and. &
iflag.eq.1) then
bb = becsum (ijh, na, is)
zsum = zsum + 0.5d0 * (alphasum (ijh, 1, na, is) * u1 + &
alphasum (ijh, 2, na, is) * u2 + alphasum (ijh, 3, na, &
is) * u3)
zsum = zsum + 0.5d0 * &
( alphasum (ijh, 1, na, is) * u1 &
+ alphasum (ijh, 2, na, is) * u2 &
+ alphasum (ijh, 3, na, is) * u3)
u1 = u1 * fact
u2 = u2 * fact
u3 = u3 * fact
alpha_0 = xq (1) * u1 + xq (2) * u2 + xq (3) * u3
alpha_0 = xq(1)*u1 + xq(2)*u2 + xq(3)*u3
do ig = 1, ngm
alpha = alpha_0 + g (1, ig) * u1 + g (2, ig) * u2 + g (3, &
ig) * u3
aux (ig, is, ipert) = aux (ig, is, ipert) + (zsum + &
alpha * bb) * sk (ig)
alpha = alpha_0 + &
g(1,ig)*u1 + g(2,ig)*u2 + g(3,ig)*u3
aux(ig,is,ipert) = aux(ig,is,ipert) + &
(zsum + alpha*bb) * sk(ig)
enddo
else
call ZAXPY (ngm, zsum, sk, 1, aux (1, is, ipert), &
1)
call ZAXPY (ngm, zsum, sk, 1, aux(1,is,ipert), 1)
endif
enddo
enddo
@ -164,36 +153,31 @@ subroutine addusddens (drhoscf, dbecsum, irr, mode0, npe, iflag)
!
do ipert = 1, npert (irr)
mu = mode0 + ipert
if (iflag.eq.0) call davcio (drhous, lrdrhous, iudrhous, mu, &
- 1)
do is = 1, nspin
call setv (2 * nrxx, 0.d0, qg, 1)
call setv (2 * nrxx, 0.d0, psic, 1)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is, ipert)
psic (nl (ig) ) = aux (ig, is, ipert)
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
if (iflag.eq.0) then
do ir = 1, nrxx
drhoscf (ir, is, ipert) = drhoscf (ir, is, ipert) + 2.d0 * qg ( &
ir) + drhous (ir, is)
enddo
else
do ir = 1, nrxx
drhoscf (ir, is, ipert) = drhoscf (ir, is, ipert) + 2.d0 * qg ( &
ir)
enddo
endif
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
call DAXPY (2*nrxx, 2.d0, psic, 1, drhoscf(1,is,ipert), 1)
enddo
enddo
if (iflag.eq.0) deallocate (drhous)
if (.not.lgamma) deallocate (qpg)
deallocate (ylmk0)
deallocate (qmod)
deallocate (qg)
deallocate (sk)
deallocate (aux)
if (iflag.eq.0) then
allocate (drhous( nrxx, nspin))
do ipert = 1, npert (irr)
mu = mode0 + ipert
call davcio (drhous, lrdrhous, iudrhous, mu, -1)
call DAXPY (2*nrxx*nspin, 1.d0, drhous, 1, drhoscf(1,1,ipert), 1)
end do
deallocate (drhous)
end if
call stop_clock ('addusddens')
return
end subroutine addusddens

29
PH/addusldos.f90
View File

@ -39,20 +39,18 @@ subroutine addusldos (ldos, becsum1)
! the spherical harmonics
! the modulus of G
complex(kind=DP), allocatable :: qg (:), aux (:,:)
complex(kind=DP), allocatable :: aux (:,:)
! auxiliary variable for FFT
! auxiliary variable for rho(G)
allocate (aux ( ngm , nspin))
allocate (qg ( nrxx))
allocate (qmod( ngm))
allocate (ylmk0 ( ngm , lqx * lqx))
call setv (2 * ngm * nspin, 0.d0, aux, 1)
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
@ -68,9 +66,11 @@ subroutine addusldos (ldos, becsum1)
!
do is = 1, nspin
do ig = 1, ngm
aux (ig, is) = aux (ig, is) + qgm (ig) * becsum1 (ijh, na, &
is) * (eigts1 (ig1 (ig), na) * eigts2 (ig2 (ig), na) &
* eigts3 (ig3 (ig), na) )
aux (ig, is) = aux (ig, is) + &
qgm (ig) * becsum1 (ijh, na, is) * &
( eigts1 (ig1 (ig), na) * &
eigts2 (ig2 (ig), na) * &
eigts3 (ig3 (ig), na) )
enddo
enddo
endif
@ -84,23 +84,16 @@ subroutine addusldos (ldos, becsum1)
!
if (okvan) then
do is = 1, nspin
call setv (2 * nrxx, 0.d0, qg, 1)
psic (:) = (0.d0,0.d0)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
psic (nl (ig) ) = aux (ig, is)
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do ir = 1, nrxx
ldos (ir, is) = ldos (ir, is) + DREAL (qg (ir) )
enddo
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
call DAXPY (nrxx, 1.d0, psic, 2, ldos(1,is), 2 )
enddo
endif
deallocate (ylmk0)
deallocate (qmod)
deallocate (qg)
deallocate (aux)
return
end subroutine addusldos

106
PH/compute_alphasum.f90
View File

@ -34,69 +34,67 @@ integer :: ik, ikk, ikq, ijkb0, ijh, ikb, jkb, ih, jh, na, nt, &
! counter on polarizations
! counter on bands
real(kind=DP) :: wgg1
! auxiliary weight
if (.not.okvan) return
call setv ( (nhm * (nhm + 1) ) / 2 * nat * 3 * nspin, 0.d0, &
alphasum, 1)
call setv ( (nhm*(nhm+1))/2*nat*3*nspin, 0.d0, alphasum, 1)
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
if (tvanp (nt) ) then
do na = 1, nat
if (ityp (na) .eq.nt) then
ijh = 0
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ijh = ijh + 1
do ibnd = 1, nbnd_occ (ikk)
wgg1 = wg (ibnd, ikk)
do ipol = 1, 3
alphasum (ijh, ipol, na, current_spin) = alphasum (ijh, &
ipol, na, current_spin) + 2.d0 * wgg1 * DREAL (conjg ( &
alphap (ikb, ibnd, ipol, ik) ) * becp1 (ikb, ibnd, ik) )
enddo
enddo
do jh = 1, nh (nt)
jkb = ijkb0 + jh
if (jh.gt.ih) ijh = ijh + 1
do ibnd = 1, nbnd
if (jh.gt.ih) then
wgg1 = wg (ibnd, ikk)
do ipol = 1, 3
alphasum (ijh, ipol, na, current_spin) = alphasum (ijh, &
ipol, na, current_spin) + 2.d0 * wgg1 * DREAL (conjg ( &
alphap (ikb, ibnd, ipol, ik) ) * becp1 (jkb, ibnd, ik) &
+ conjg (becp1 (ikb, ibnd, ik) ) * alphap (jkb, ibnd, &
ipol, ik) )
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
if (tvanp (nt) ) then
do na = 1, nat
if (ityp (na) .eq.nt) then
ijh = 0
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ijh = ijh + 1
do ibnd = 1, nbnd_occ (ikk)
wgg1 = wg (ibnd, ikk)
do ipol = 1, 3
alphasum(ijh,ipol,na,current_spin) = &
alphasum(ijh,ipol,na,current_spin) + 2.d0 * wgg1 * &
DREAL (conjg (alphap (ikb,ibnd,ipol,ik) ) * &
becp1 (ikb,ibnd,ik) )
enddo
enddo
do jh = 1, nh (nt)
jkb = ijkb0 + jh
if (jh.gt.ih) ijh = ijh + 1
do ibnd = 1, nbnd
if (jh.gt.ih) then
wgg1 = wg (ibnd, ikk)
do ipol = 1, 3
alphasum(ijh,ipol,na,current_spin) = &
alphasum(ijh,ipol,na,current_spin) + &
2.d0 * wgg1 * &
DREAL (conjg ( alphap (ikb,ibnd,ipol,ik) ) * &
becp1 (jkb,ibnd,ik) + &
conjg ( becp1 (ikb,ibnd,ik) ) * &
alphap (jkb,ibnd,ipol,ik) )
enddo
endif
enddo
enddo
enddo
ijkb0 = ijkb0 + nh (nt)
endif
enddo
else
do na = 1, nat
if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
enddo
enddo
ijkb0 = ijkb0 + nh (nt)
endif
enddo
else
do na = 1, nat
if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
enddo
! do na=1,nat
! nt=ityp(na)

104
PH/compute_becsum.f90
View File

@ -25,69 +25,59 @@ use parameters, only : DP
use phcom
implicit none
integer :: ik, ikk, ikq, ijkb0, ijh, ikb, jkb, ih, jh, na, nt, &
ibnd
! counter on k points
! counters on beta functions
! counters on beta functions
! counters for atoms
! counter on bands
integer :: ik, ikk, ikq, ijkb0, ijh, ikb, jkb, ih, jh, na, nt, ibnd
! counter on k points, beta functions, atoms and bands
real(kind=DP) :: wgg1 ! auxiliary weight
real(kind=DP) :: wgg1
! auxiliary weight
if (.not.okvan) return
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum, 1)
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
if (tvanp (nt) ) then
do na = 1, nat
if (ityp (na) .eq.nt) then
ijh = 0
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ijh = ijh + 1
do ibnd = 1, nbnd_occ (ikk)
wgg1 = wg (ibnd, ikk)
becsum (ijh, na, current_spin) = becsum (ijh, na, &
current_spin) + wgg1 * DREAL (conjg (becp1 (ikb, ibnd, ik) ) &
* becp1 (ikb, ibnd, ik) )
enddo
do jh = 1, nh (nt)
jkb = ijkb0 + jh
if (jh.gt.ih) ijh = ijh + 1
do ibnd = 1, nbnd
if (jh.gt.ih) then
wgg1 = wg (ibnd, ikk)
becsum (ijh, na, current_spin) = becsum (ijh, na, &
current_spin) + wgg1 * 2.d0 * DREAL (conjg (becp1 (ikb, &
ibnd, ik) ) * becp1 (jkb, ibnd, ik) )
endif
enddo
enddo
enddo
ijkb0 = ijkb0 + nh (nt)
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
if (tvanp (nt) ) then
do na = 1, nat
if (ityp (na) .eq.nt) then
ijh = 0
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ijh = ijh + 1
do ibnd = 1, nbnd_occ (ikk)
wgg1 = wg (ibnd, ikk)
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + wgg1 * &
DREAL ( conjg(becp1(ikb,ibnd,ik)) * becp1(ikb,ibnd,ik) )
enddo
do jh = 1, nh (nt)
jkb = ijkb0 + jh
if (jh.gt.ih) ijh = ijh + 1
do ibnd = 1, nbnd
if (jh.gt.ih) then
wgg1 = wg (ibnd, ikk)
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + wgg1 * 2.d0 * &
DREAL ( conjg(becp1(ikb,ibnd,ik)) * &
becp1(jkb,ibnd,ik) )
endif
enddo
enddo
enddo
ijkb0 = ijkb0 + nh (nt)
endif
enddo
else
do na = 1, nat
if (ityp(na).eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
else
do na = 1, nat
if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
enddo
! do na=1,nat
! nt=ityp(na)

3
PH/drho.f90
View File

@ -23,8 +23,7 @@ subroutine drho
use phcom
implicit none
integer :: nt, mode, mu, na, is, ir, irr, iper, npe, nrstot, nu_i, &
nu_j
integer :: nt, mode, mu, na, is, ir, irr, iper, npe, nrstot, nu_i, nu_j
! counter on atomic types
! counter on modes
! counter on atoms and polarizations

30
PH/dv_of_drho.f90
View File

@ -55,16 +55,14 @@ subroutine dv_of_drho (mode, dvscf, flag)
if (nlcc_any.and.flag) then
call addcore (mode, drhoc)
do is = 1, nspin
call DAXPY (nrxx, fac, rho_core, 1, rho (1, is), 1)
call DAXPY (2 * nrxx, fac, drhoc, 1, dvscf (1, is), 1)
call DAXPY (nrxx, fac, rho_core, 1, rho (1, is), 1)
call DAXPY (2*nrxx, fac, drhoc, 1, dvscf (1, is), 1)
enddo
endif
do is = 1, nspin
do is1 = 1, nspin
do ir = 1, nrxx
dvaux (ir, is) = dvaux (ir, is) + dmuxc (ir, is, is1) * dvscf (ir, &
is1)
dvaux(ir,is) = dvaux(ir,is) + dmuxc(ir,is,is1) * dvscf(ir,is1)
enddo
enddo
enddo
@ -78,37 +76,33 @@ subroutine dv_of_drho (mode, dvscf, flag)
alat, omega, dvaux)
if (nlcc_any.and.flag) then
do is = 1, nspin
call DAXPY (nrxx, - fac, rho_core, 1, rho (1, is), 1)
call DAXPY (2 * nrxx, - fac, drhoc, 1, dvscf (1, is), 1)
call DAXPY (nrxx, -fac, rho_core, 1, rho (1, is), 1)
call DAXPY (2*nrxx, -fac, drhoc, 1, dvscf (1, is), 1)
enddo
endif
!
! copy the total (up+down) delta rho in dvscf(*,1) and go to G-space
!
do is = 2, nspin
call DAXPY (2 * nrxx, 1.d0, dvscf (1, is), 1, dvscf (1, 1), &
1)
call DAXPY (2 * nrxx, 1.d0, dvscf(1,is), 1, dvscf(1,1), 1)
enddo
call cft3 (dvscf, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
call cft3 (dvscf, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
!
! hartree contribution is computed in reciprocal space
!
do is = 1, nspin
call cft3 (dvaux (1, is), nr1, nr2, nr3, nrx1, nrx2, nrx3, &
- 1)
call cft3 (dvaux (1, is), nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
do ig = 1, ngm
qg2 = (g (1, ig) + xq (1) ) **2 + (g (2, ig) + xq (2) ) **2 + &
(g (3, ig) + xq (3) ) **2
qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
if (qg2.gt.1.d-8) then
dvaux (nl (ig), is) = dvaux (nl (ig), is) + e2 * fpi * dvscf ( &
nl (ig), 1) / (tpiba2 * qg2)
dvaux(nl(ig),is) = dvaux(nl(ig),is) + &
e2 * fpi * dvscf(nl(ig),1) / (tpiba2 * qg2)
endif
enddo
!
! and transformed back to real space
!
call cft3 (dvaux (1, is), nr1, nr2, nr3, nrx1, nrx2, nrx3, &
+ 1)
call cft3 (dvaux (1, is), nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)
enddo
!
! at the end the two contributes are added

4
PH/dvanqq.f90
View File

@ -90,7 +90,6 @@ subroutine dvanqq
call ylmr2 (lqx * lqx, ngm, g, gg, ylmk0)
do ig = 1, ngm
qmodg (ig) = sqrt (gg (ig) )
enddo
if (.not.lgamma) then
call setqmod (ngm, xq, g, qmod, qpg)
@ -125,7 +124,6 @@ subroutine dvanqq
aux5 (ig, na, ipol) = sk (ig) * (g (ipol, ig) + xq (ipol) )
enddo
enddo
enddo
do ntb = 1, ntyp
if (tvanp (ntb) ) then
@ -150,7 +148,6 @@ subroutine dvanqq
aux1 (ig) = qgmq (ig) * eigts1 (ig1 (ig), nb) &
* eigts2 (ig2 (ig), nb) &
* eigts3 (ig3 (ig), nb)
enddo
do na = 1, nat
fact = eigqts (na) * conjg (eigqts (nb) )
@ -158,7 +155,6 @@ subroutine dvanqq
! nb is the atom of the augmentation function
!
do ipol = 1, 3
int2 (ih, jh, ipol, na, nb) = fact * fact1 * &
ZDOTC (ngm, aux1, 1, aux5(1,na,ipol), 1)
do jpol = 1, 3

104
PH/ef_shift.f90
View File

@ -25,19 +25,18 @@ implicit none
integer :: npe
! input: the number of perturbation
complex(kind=DP) :: drhoscf (nrxx, nspin, npe), ldos (nrxx, nspin), &
ldoss (nrxxs, nspin)
! inp/out:the change of the charge
! inp: local DOS at Ef
! inp: local DOS at Ef without augme
complex(kind=DP) :: drhoscf(nrxx,nspin,npe), &
ldos(nrxx,nspin), ldoss(nrxxs,nspin)
! inp/out:the change of the charge
! inp: local DOS at Ef
! inp: local DOS at Ef without augme
real(kind=DP) :: dos_ef
! inp: density of states at Ef
! inp: density of states at Ef
integer :: irr
! inp: index of the current irr. rep.
! inp: index of the current irr. rep.
logical :: flag
! inp: if true the eigenfunctions are updated
! inp: if true the eigenfunctions are updated
!
! local variables
!
@ -68,86 +67,61 @@ call start_clock ('ef_shift')
if (.not.flag) then
write (6, * )
do ipert = 1, npert (irr)
delta_n = (0.d0, 0.d0)
do is = 1, nspin
call cft3 (drhoscf (1, is, ipert), nr1, nr2, nr3, nrx1, nrx2, &
nrx3, - 1)
if (gg (1) .lt.1.0d-8) delta_n = delta_n + omega * drhoscf (nl &
(1), is, ipert)
call cft3 (drhoscf (1, is, ipert), nr1, nr2, nr3, nrx1, nrx2, &
nrx3, + 1)
enddo
call reduce (2, delta_n)
def (ipert) = - delta_n / dos_ef
!check
! write (6,*) 'delta_n , dos_ef, def(ipert)'
!scal write (6,*) 'delta_n , dos_ef, def(ipert)'
! write (6,*) delta_n , dos_ef, def(ipert)
!check
delta_n = (0.d0, 0.d0)
do is = 1, nspin
call cft3 (drhoscf(1,is,ipert), nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
if (gg(1).lt.1.0d-8) delta_n = delta_n + omega*drhoscf(nl(1),is,ipert)
call cft3 (drhoscf(1,is,ipert), nr1, nr2, nr3, nrx1, nrx2, nrx3, +1)
enddo
call reduce (2, delta_n)
def (ipert) = - delta_n / dos_ef
enddo
!
! symmetrizes the Fermi energy shift
!
call sym_def (def, irr)
write (6, '(5x,"Pert. #",i3, &
& ": Fermi energy shift (Ryd) =", &
& 2f10.4)') (ipert, def (ipert) , ipert = 1, npert (irr) )
write (6, '(5x,"Pert. #",i3,": Fermi energy shift (Ryd) =", 2f10.4)') &
(ipert, def (ipert) , ipert = 1, npert (irr) )
!
! corrects the density response accordingly...
!
do ipert = 1, npert (irr)
call ZAXPY (nrxx * nspin, def (ipert), ldos, 1, drhoscf (1, 1, &
ipert), 1)
!check
! delta_n= (0.d0,0.d0)
! do is=1,nspin
! call cft3(drhoscf(1,is,ipert),nr1,nr2,nr3,nrx1,nrx2,nrx3,-1
! if (gg(1).lt.1.0d-8)
! + delta_n = delta_n + omega*drhoscf(nl(1),is,ipert)
! call cft3(drhoscf(1,is,ipert),nr1,nr2,nr3,nrx1,nrx2,nrx3,+1
! end do
! call reduce(2,delta_n)
! write (6,*) 'new delta_n , dos_ef, def(ipert)',nd_nmbr
! write (6,*) delta_n , dos_ef, def(ipert)
!check
call ZAXPY (nrxx*nspin, def(ipert), ldos, 1, drhoscf(1,1,ipert), 1)
enddo
else
!
! does the same for perturbed wfc
!
do ik = 1, nksq
npw = ngk (ik)
npw = ngk (ik)
!
! reads unperturbed wavefuctions psi_k in G_space, for all bands
!
ikrec = ik
if (nksq.gt.1) call davcio (evc, lrwfc, iuwfc, ikrec, - 1)
ikrec = ik
if (nksq.gt.1) call davcio (evc, lrwfc, iuwfc, ikrec, - 1)
!
! reads delta_psi from iunit iudwf, k=kpoint
!
do ipert = 1, npert (irr)
nrec = (ipert - 1) * nksq + ik
if (nksq.gt.1.or.npert (irr) .gt.1) call davcio (dpsi, lrdwf, &
iudwf, nrec, - 1)
do ibnd = 1, nbnd_occ (ik)
wfshift = 0.5d0 * def (ipert) * w0gauss ( (ef - et (ibnd, ik) ) &
/ degauss, ngauss) / degauss
call ZAXPY (npw, wfshift, evc (1, ibnd), 1, dpsi (1, ibnd), &
1)
enddo
do ipert = 1, npert (irr)
nrec = (ipert - 1) * nksq + ik
if (nksq.gt.1.or.npert(irr).gt.1) &
call davcio (dpsi, lrdwf, iudwf, nrec, -1)
do ibnd = 1, nbnd_occ (ik)
wfshift = 0.5d0 * def(ipert) * &
w0gauss( (ef-et(ibnd,ik))/degauss, ngauss) / degauss
call ZAXPY (npw, wfshift, evc(1,ibnd), 1, dpsi(1,ibnd), 1)
enddo
!
! writes corrected delta_psi to iunit iudwf, k=kpoint,
!
if (nksq.gt.1.or.npert (irr) .gt.1) call davcio (dpsi, lrdwf, &
iudwf, nrec, + 1)
enddo
enddo
do ipert = 1, npert (irr)
do is = 1, nspin
call ZAXPY (nrxxs, def (ipert), ldoss (1, is), 1, drhoscf (1, &
is, ipert), 1)
enddo
if (nksq.gt.1.or.npert(irr).gt.1) &
call davcio (dpsi, lrdwf, iudwf, nrec, +1)
enddo
enddo
do ipert = 1, npert (irr)
do is = 1, nspin
call ZAXPY (nrxxs, def(ipert), ldoss(1,is), 1, drhoscf(1,is,ipert), 1)
enddo
enddo
endif
call stop_clock ('ef_shift')

16
PH/localdos.f90
View File

@ -75,8 +75,7 @@ subroutine localdos (ldos, ldoss, dos_ef)
call start_clock ('localdos')
allocate (becsum1( (nhm * (nhm + 1)) / 2, nat, nspin))
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum1, &
1)
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum1, 1)
call setv (2 * nrxx * nspin, 0.d0, ldos, 1)
call setv (2 * nrxxs * nspin, 0.d0, ldoss, 1)
dos_ef = 0.d0
@ -100,29 +99,23 @@ subroutine localdos (ldos, ldoss, dos_ef)
call ccalbec (nkb, npwx, npw, nbnd, becp, vkb, evc)
do ibnd = 1, nbnd_occ (ik)
wdelta = w0gauss ( (ef - et (ibnd, ik) ) / degauss, ngauss) &
/ degauss
wdelta = w0gauss ( (ef-et(ibnd,ik)) / degauss, ngauss) / degauss
!
! unperturbed wf from reciprocal to real space
!
call setv (2 * nrxxs, 0.d0, psic, 1)
do ig = 1, npw
psic (nls (igk (ig) ) ) = evc (ig, ibnd)
enddo
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 1)
w1 = weight * wdelta / omega
do j = 1, nrxxs
ldoss (j, current_spin) = ldoss (j, current_spin) + w1 * (DREAL ( &
psic (j) ) **2 + DIMAG (psic (j) ) **2)
ldoss (j, current_spin) = ldoss (j, current_spin) + &
w1 * (DREAL ( psic (j) ) **2 + DIMAG (psic (j) ) **2)
enddo
!
! If we have a US pseudopotential we compute here the sumbec term
!
w1 = weight * wdelta
ijkb0 = 0
do nt = 1, ntyp
@ -152,7 +145,6 @@ subroutine localdos (ldos, ldoss, dos_ef)
if (ityp (na) .eq.nt) ijkb0 = ijkb0 + nh (nt)
enddo
endif
enddo
dos_ef = dos_ef + weight * wdelta
enddo

30
PH/newdq.f90
View File

@ -82,7 +82,6 @@ subroutine newdq (dvscf, npe)
do ig = 1, ngm
qmod (ig) = sqrt (gg (ig) )
enddo
endif
!
! and for each perturbation of this irreducible representation
@ -90,6 +89,7 @@ subroutine newdq (dvscf, npe)
! the Q functions
!
do ipert = 1, npe
do is = 1, nspin
do ir = 1, nrxx
veff (ir) = dvscf (ir, is, ipert)
@ -98,8 +98,8 @@ subroutine newdq (dvscf, npe)
do ig = 1, ngm
aux2 (ig, is) = veff (nl (ig) )
enddo
enddo
do nt = 1, ntyp
if (tvanp (nt) ) then
do ih = 1, nh (nt)
@ -108,48 +108,40 @@ subroutine newdq (dvscf, npe)
do na = 1, nat
if (ityp (na) .eq.nt) then
do ig = 1, ngm
aux1(ig)=qgm(ig)*eigts1(ig1(ig),na)*eigts2(ig2( &
ig),na)*eigts3(ig3(ig),na)*eigqts(na)
aux1(ig) = qgm(ig) * eigts1(ig1(ig),na) * &
eigts2(ig2(ig),na) * &
eigts3(ig3(ig),na) * &
eigqts(na)
enddo
do is = 1, nspin
int3(ih,jh,ipert,na,is)=omega*ZDOTC(ngm,aux1,1, &
aux2 (1, is), 1)
int3(ih,jh,ipert,na,is) = omega * &
ZDOTC(ngm,aux1,1,aux2(1,is),1)
enddo
!
! ps contain the integral of V_loc and Q_nm
!
endif
enddo
enddo
enddo
do na = 1, nat
if (ityp (na) .eq.nt) then
if (ityp(na) .eq.nt) then
!
! We use the symmetry properties of the ps factor
!
do ih = 1, nh (nt)
do jh = ih, nh (nt)
do is = 1, nspin
int3 (jh,ih,ipert,na,is)=int3(ih,jh,ipert,na,is)
int3(jh,ih,ipert,na,is) = int3(ih,jh,ipert,na,is)
enddo
enddo
enddo
endif
enddo
endif
enddo
enddo
#ifdef __PARA
call reduce (2 * nhm * nhm * 3 * nat * nspin, int3)
#endif
! do ih = 1,nh(1)
! do jh=1,nh(1)
! write(6,*) int3(jh,ih,1,1,1)
! enddo
! enddo
! call stop_ph(.true.)
if (.not.lgamma) deallocate (qg)
deallocate (qmod)
deallocate (ylmk0)

22
PH/phq_init.f90
View File

@ -124,10 +124,9 @@ subroutine phq_init
if (.not.lgamma) then
call gk_sort (xk (1, ikq), ngm, g, ecutwfc / tpiba2, npwq, igkq, g2kin)
if (nksq.gt.1) write (iunigk) npwq, igkq
if (abs (xq (1) - (xk (1, ikq) - xk (1, ikk) ) ) &
.gt.1.d-8.or.abs (xq (2) - (xk (2, ikq) - xk (2, ikk) ) ) &
.gt.1.d-8.or.abs (xq (3) - (xk (3, ikq) - xk (3, ikk) ) ) &
.gt.1.d-8) then
if (abs (xq (1) - (xk (1, ikq) - xk (1, ikk) ) ) .gt.1.d-8 .or. &
abs (xq (2) - (xk (2, ikq) - xk (2, ikk) ) ) .gt.1.d-8 .or. &
abs (xq (3) - (xk (3, ikq) - xk (3, ikk) ) ) .gt.1.d-8) then
write (6, * ) ikk, ikq, nksq
write (6, * ) (xq (ipol), ipol = 1, 3)
write (6, * ) (xk (ipol, ikq), ipol = 1, 3)
@ -140,17 +139,12 @@ subroutine phq_init
!
call init_us_2 (npw, igk, xk (1, ikk), vkb)
!
! e) we compute also the becp terms which are used in the rest of
! the code
!
! read the wavefunctions at k
!
call davcio (evc, lrwfc, iuwfc, ikk, - 1)
!
! if there is only one k-point the wavefunctions are read once here
!
if (nksq.eq.1.and..not.lgamma) call davcio (evq, lrwfc, iuwfc, ikq, -1)
! e) we compute the becp terms which are used in the rest of
! the code
call ccalbec (nkb, npwx, npw, nbnd, becp1 (1, 1, ik), vkb, evc)
!
! e') we compute the derivative of the becp term with respect to an
@ -164,9 +158,13 @@ subroutine phq_init
enddo
enddo
call ccalbec (nkb, npwx, npw, nbnd, alphap(1,1,ipol,ik), vkb, aux1)
enddo
!
! if there is only one k-point the k+q wavefunctions are read once here
!
if (nksq.eq.1.and..not.lgamma) call davcio (evq, lrwfc, iuwfc, ikq, -1)
enddo
deallocate (aux1)

201
PP/chdens.f90
View File

@ -225,10 +225,9 @@ subroutine do_chdens
! reading the rest of input (spanning vectors, origin, number-of points)
if (iflag.lt.4) then
read (inunit, *, err = 1100, iostat = ios) (e (ipol, 1), &
ipol = 1, 3)
if (e (1, 1) **2 + e (2, 1) **2 + e (3, 1) **2.lt.1d-3) call &
errore ('chdens', 'zero vector', 1)
read (inunit, *, err = 1100, iostat = ios) (e (ipol,1), ipol = 1, 3)
if (e(1,1)**2 + e(2,1)**2 + e(3,1)**2 .lt. 1d-3) &
call errore ('chdens', 'zero vector', 1)
endif
if (iflag.eq.1) then
!
@ -242,44 +241,39 @@ subroutine do_chdens
!
! reading for the 2D and 3D plots
!
read (inunit, *, err = 1100, iostat = ios) (e (ipol, 2), &
ipol = 1, 3)
read (inunit, *, err = 1100, iostat = ios) (e(ipol,2), ipol=1,3)
!
! here we control that the vectors are not on the same line
!
if ( (abs (e (1, 1) * e (2, 2) - e (2, 1) * e (1, 2) ) .lt.1e-7 ) &
.and. (abs (e (3, 1) * e (1, 2) - e (1, 1) * e (3, 2) ) .lt.1e-7) &
.and. (abs (e (3, 1) * e (2, 2) - e (2, 1) * e (3, 2) ) .lt.1e-7) ) &
call errore ('chdens', 'vectors on the same line', 1)
if ( (abs(e(1,1)*e(2,2) - e(2,1)*e(1,2) ) .lt. 1e-7) .and. &
(abs(e(3,1)*e(1,2) - e(1,1)*e(3,2) ) .lt. 1e-7) .and. &
(abs(e(3,1)*e(2,2) - e(2,1)*e(3,2) ) .lt. 1e-7) ) &
call errore ('chdens', 'vectors on the same line', 1)
!
! and here that they are orthogonal
!
if (abs (e (1, 1) * e (1, 2) + e (2, 1) * e (2, 2) + e (3, 1) &
& * e (3, 2) ) .gt.1e-4) call errore ('chdens', &
'vectors are not orthogonal', 1)
if (abs(e(1,1)*e(1,2) + e(2,1)*e(2,2) + e(3,1)*e(3,2)) .gt. 1e-4) &
call errore ('chdens', 'vectors are not orthogonal', 1)
!
if (iflag.eq.3) then
!
! reading for the 3D plot
!
read (inunit, *, err = 1100, iostat = ios) (e (ipol, 3), &
ipol = 1, 3)
read (inunit, *, err=1100, iostat=ios) (e(ipol,3), ipol=1,3)
!
! here we control that the vectors are not on the same line
!
if ( (abs (e (1, 1) * e (2, 3) - e (2, 1) * e (1, 3) ) &
.lt.1e-7) .and. (abs (e (3, 1) * e (1, 3) - e (1, 1) &
* e (3, 3) ) .lt.1e-7) .and. (abs (e (3, 1) * e (2, 3) &
- e (2, 1) * e (3, 3) ) .lt.1e-7) ) call errore ('chdens', &
'vectors on the same line', 2)
if ( (abs(e(1,1)*e(2,3) - e(2,1)*e(1,3)) .lt. 1e-7) .and. &
(abs(e(3,1)*e(1,3) - e(1,1)*e(3,3)) .lt. 1e-7) .and. &
(abs(e(3,1)*e(2,3) - e(2,1)*e(3,3)) .lt. 1e-7) ) &
call errore ('chdens', 'vectors on the same line', 2)
!
! and here that they are orthogonal
!
if (abs (e (1, 1) * e (1, 3) + e (2, 1) * e (2, 3) + e (3, &
1) * e (3, 3) ) .gt.1e-4.or.abs (e (1, 2) * e (1, 3) &
+ e (2, 2) * e (2, 3) + e (3, 2) * e (3, 3) ) .gt.1e-4) &
call errore ('chdens', 'vectors are not orthogonal', 2)
if (abs(e(1,1)*e(1,3) + e(2,1)*e(2,3) + e(3,1)*e(3,3)) .gt. 1e-4 .or. &
abs(e(1,2)*e(1,3) + e(2,2)*e(2,3) + e(3,2)*e(3,3)) .gt. 1e-4) &
call errore ('chdens', 'vectors are not orthogonal', 2)
endif
read (inunit, *, err = 1100, iostat = ios) (x0 (ipol), ipol = &
1, 3)
@ -296,15 +290,13 @@ subroutine do_chdens
endif
! check for plot_out
if (plot_out.lt.0.or.plot_out.gt.4) call errore ('chdens', &
'plot_out wrong', 1)
if (plot_out.lt.0.or.plot_out.gt.4) call errore ('chdens','plot_out wrong',1)
!
! Read the header and allocate objects
!
call plot_io (filepp (1), title, nrx1, nrx2, nrx3, nr1, nr2, &
nr3, nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rhodum, 0)
call plot_io (filepp (1), title, nrx1, nrx2, nrx3, nr1, nr2, nr3, &
nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rhodum, 0)
!
allocate(tau (3, nat))
allocate(ityp(nat))
@ -320,11 +312,9 @@ subroutine do_chdens
endif
nspin = 1
if (ibrav.gt.0) call latgen (ibrav, celldm, at (1, 1), &
& at (1, 2), at (1, 3) )
call recips (at (1, 1), at (1, 2), at (1, 3), bg (1, 1), bg (1, 2) &
, bg (1, 3) )
call volume (alat, at (1, 1), at (1, 2), at (1, 3), omega)
if (ibrav.gt.0) call latgen (ibrav, celldm, at(1,1), at(1,2), at(1,3) )
call recips (at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
call volume (alat, at(1,1), at(1,2), at(1,3), omega)
call set_fft_dim
@ -334,9 +324,9 @@ subroutine do_chdens
!
! Read first file
!
call plot_io (filepp (1), title, nrx1, nrx2, nrx3, nr1, nr2, &
nr3, nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rho, -1)
call plot_io (filepp (1), title, nrx1, nrx2, nrx3, nr1, nr2, nr3, &
nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rho, -1)
!
do ir = 1, nrxx
psic (ir) = weight (1) * cmplx (rho (ir, 1),0.d0)
@ -370,8 +360,7 @@ subroutine do_chdens
enddo
!
do ir = 1, nrxx
psic (ir) = psic (ir) + weight (ifile) * cmplx (rho (ir, 1), &
0.d0)
psic(ir) = psic(ir) + weight(ifile) * cmplx(rho(ir,1),0.d0)
enddo
enddo
@ -380,8 +369,7 @@ subroutine do_chdens
!
if (fileout.ne.' ') then
ounit = 1
open (unit = ounit, file = fileout, form = 'formatted', status &
= 'unknown')
open (unit=ounit, file=fileout, form='formatted', status='unknown')
write (6, '(5x,"Writing data on file ",a)') fileout
else
ounit = 6
@ -410,11 +398,8 @@ subroutine do_chdens
!
if (iflag.lt.4) then
m1 = sqrt (e (1, 1) **2 + e (2, 1) **2 + e (3, 1) **2)
if (iflag.ge.2) m2 = sqrt (e (1, 2) **2 + e (2, 2) **2 + e (3, &
2) **2)
if (iflag.eq.3) m3 = sqrt (e (1, 3) **2 + e (2, 3) **2 + e (3, &
3) **2)
if (iflag.ge.2) m2 = sqrt(e(1,2)**2 + e(2,2)**2 + e(3,2)**2)
if (iflag.eq.3) m3 = sqrt(e(1,3)**2 + e(2,3)**2 + e(3,3)**2)
do ipol = 1, 3
e (ipol, 1) = e (ipol, 1) / m1
@ -450,7 +435,7 @@ subroutine do_chdens
!
! bring the quantity in real space and write the output file
!
call setv (2 * nrxx, 0.d0, psic, 1)
psic(:) = (0.d0,0.d0)
do ig = 1, ngm
psic (nl (ig) ) = vgc (ig)
enddo
@ -459,19 +444,17 @@ subroutine do_chdens
do ir = 1, nrxx
rho (ir, 1) = real (psic (ir) )
enddo
call plot_io (filepol, title, nrx1, nrx2, nrx3, nr1, nr2, &
nr3, nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rho, + 1)
call plot_io (filepol, title, nrx1, nrx2, nrx3, nr1, nr2, nr3, &
nat, ntyp, ibrav, celldm, at, gcutm, dual, ecutwfc, &
plot_num, atm, ityp, zv, tau, rho, + 1)
endif
endif
!
! And now the plot
!
if (iflag.eq.1) then
call plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
ounit)
call plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, ounit)
elseif (iflag.eq.2) then
@ -479,8 +462,7 @@ subroutine do_chdens
at, nat, tau, atm, ityp, output_format, ounit)
if (output_format.eq.2) then
write (ounit, '(i4)') nat
write (ounit, '(3f8.4,i3)') ( (tau (ipol, na) , ipol = 1, 3) &
, 1, na = 1, nat)
write (ounit, '(3f8.4,i3)') ( (tau(ipol,na), ipol=1,3), 1, na=1,nat)
write (ounit, '(f10.6)') celldm (1)
write (ounit, '(3(3f12.6/))') at
endif
@ -517,9 +499,7 @@ subroutine do_chdens
elseif (iflag.eq.4) then
radius = radius / alat
call plot_2ds (nx, ny, radius, ngm, g, vgc, output_format, &
ounit)
call plot_2ds (nx, ny, radius, ngm, g, vgc, output_format, ounit)
else
call errore ('chdens', 'wrong iflag', 1)
@ -531,8 +511,7 @@ subroutine do_chdens
end subroutine do_chdens
!
!-----------------------------------------------------------------------
subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
ounit)
subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, ounit)
!-----------------------------------------------------------------------
!
use parameters, only : DP
@ -567,7 +546,7 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
complex(kind=DP) :: rho0g, carica (nx)
deltax = m1 / (nx - 1)
call setv (2 * nx, 0.d0, carica, 1)
carica(:) = (0.d0,0.d0)
if (plot_out.eq.1) then
do i = 1, nx
xi = x0 (1) + (i - 1) * deltax * e (1)
@ -580,10 +559,8 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
!
! NB: G are in 2pi/alat units, r are in alat units
!
arg = 2.d0 * pi * (xi * g (1, ig) + yi * g (2, ig) + zi * g (3, &
ig) )
carica (i) = carica (i) + vgc (ig) * cmplx (cos (arg), sin ( &
arg) )
arg = 2.d0 * pi * ( xi*g(1,ig) + yi*g(2,ig) + zi*g(3,ig) )
carica(i) = carica(i) + vgc (ig) * cmplx(cos(arg),sin(arg))
enddo
enddo
else
@ -597,16 +574,15 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
enddo
! G!=0 terms
do ig = 2, ngm
arg = 2.d0 * pi * (x0 (1) * g (1, ig) + x0 (2) * g (2, ig) &
+ x0 (3) * g (3, ig) )
arg = 2.d0 * pi * ( x0(1)*g(1,ig) + x0(2)*g(2,ig) + x0(3)*g(3,ig) )
! This displaces the origin into x0
rho0g = vgc (ig) * cmplx (cos (arg), sin (arg) )
rho0g = vgc (ig) * cmplx(cos(arg),sin(arg))
! r =0 term
carica (1) = carica (1) + 4.d0 * pi * rho0g
! r!=0 terms
do i = 2, nx
gr = 2.d0 * pi * sqrt (g (1, ig) **2 + g (2, ig) **2 + g (3, &
ig) **2) * (i - 1) * deltax
gr = 2.d0 * pi * sqrt(g(1,ig)**2 + g(2,ig)**2 + g(3,ig)**2) * &
(i-1) * deltax
carica (i) = carica (i) + 4.d0 * pi * rho0g * sin (gr) / gr
enddo
@ -615,27 +591,24 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
!
! Here we check the value of the resulting charge
!
rhomin = 1.0d10
rhomax = - 1.0d10
rhomin = 1.0d10
rhomax = -1.0d10
rhoim = 0.d0
do i = 1, nx
rhomin = min (rhomin, DREAL (carica (i) ) )
rhomax = max (rhomax, DREAL (carica (i) ) )
rhoim = rhoim + abs (DIMAG (carica (i) ) )
enddo
rhoim = rhoim / nx
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, &
rhomax, rhoim
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, rhomax, rhoim
!
! we print the charge on output
!
if (plot_out.eq.1) then
do i = 1, nx
write (ounit, '(2f20.10)') deltax * float (i - 1) , real ( &
carica (i) )
write (ounit, '(2f20.10)') deltax*float(i-1), real(carica(i))
enddo
else
rhoint = 0.d0
@ -643,10 +616,8 @@ subroutine plot_1d (nx, m1, x0, e, ngm, g, vgc, alat, plot_out, &
!
! simple trapezoidal rule: rhoint=int carica(i) r^2(i) dr
!
rhoint = rhoint + real (carica (i) ) * ( (i - 1) * deltax) **2 &
* deltax * alat**3
write (ounit, '(3f20.10)') deltax * float (i - 1) , real ( &
carica (i) ) , rhoint
rhoint = rhoint + real(carica(i)) * (i-1)**2 * (deltax*alat)**3
write (ounit, '(3f20.10)') deltax*float(i-1), real(carica(i)), rhoint
enddo
endif
@ -671,8 +642,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
! output unit
! output format
character(len=3) :: atm(*) ! atomic symbols
real(kind=DP) :: e (3, 2), x0 (3), m1, m2, g (3, ngm), &
alat, tau (3, nat), at (3, 3)
real(kind=DP) :: e(3,2), x0(3), m1, m2, g(3,ngm), alat, tau(3,nat), at(3,3)
! vectors e1, e2 defining the plane
! origin
! modulus of e1
@ -701,7 +671,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
deltax = m1 / (nx - 1)
deltay = m2 / (ny - 1)
call setv (2 * nx * ny, 0.d0, carica, 1)
carica(:,:) = (0.d0,0.d0)
do ig = 1, ngm
!
! eigx=exp(iG*e1+iGx0), eigy=(iG*e2)
@ -725,8 +695,8 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
!
! Here we check the value of the resulting charge
!
rhomin = 1.0d10
rhomax = - 1.0d10
rhomin = 1.0d10
rhomax = -1.0d10
rhoim = 0.d0
do i = 1, nx
@ -739,8 +709,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
enddo
rhoim = rhoim / nx / ny
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, &
rhomax, rhoim
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, rhomax, rhoim
print '(5x,"Output format: ",i3)', output_format
!
@ -752,7 +721,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
!
! write(ounit,'(2i6)') nx,ny
do i = 1, nx
write (ounit, '(e25.14)') (DREAL (carica (i, j) ) , j = 1, ny)
write (ounit, '(e25.14)') ( DREAL(carica(i,j)), j = 1, ny )
write (ounit, * )
enddo
elseif (output_format.eq.1) then
@ -760,8 +729,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
! contour.x format
!
write (ounit, '(3i5,2e25.14)') nx, ny, 1, deltax, deltay
write (ounit, '(4e25.14)') ( (DREAL (carica (i, j) ) , j = 1, &
ny) , i = 1, nx)
write (ounit, '(4e25.14)') ( ( DREAL(carica(i,j)), j = 1, ny ), i = 1, nx )
elseif (output_format.eq.2) then
!
! plotrho format
@ -769,8 +737,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
write (ounit, '(2i4)') nx - 1, ny - 1
write (ounit, '(8f8.4)') (deltax * (i - 1) , i = 1, nx)
write (ounit, '(8f8.4)') (deltay * (j - 1) , j = 1, ny)
write (ounit, '(6e12.4)') ( (DREAL (carica (i, j) ) , i = 1, &
nx) , j = 1, ny)
write (ounit, '(6e12.4)') ( ( DREAL(carica(i,j)), i = 1, nx ), j = 1, ny )
write (ounit, '(3f8.4)') x0
write (ounit, '(3f8.4)') (m1 * e (i, 1) , i = 1, 3)
write (ounit, '(3f8.4)') (m2 * e (i, 2) , i = 1, 3)
@ -792,8 +759,7 @@ subroutine plot_2d (nx, ny, m1, m2, x0, e, ngm, g, vgc, alat, &
end subroutine plot_2d
!
!-----------------------------------------------------------------------
subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
ounit)
subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, ounit)
!-----------------------------------------------------------------------
use parameters, only : DP
!
@ -813,8 +779,7 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
integer :: i, j, ig
real(kind=DP), allocatable :: r (:,:,:)
real(kind=DP) :: theta, phi, rhomin, rhomax, rhoim, &
deltax, deltay
real(kind=DP) :: theta, phi, rhomin, rhomax, rhoim, deltax, deltay
! the point in space
! the position on the sphere
! minimum value of the charge
@ -835,7 +800,7 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
deltay = pi / (ny - 1)
call setv (2 * nx * ny, 0.d0, carica, 1)
carica(:,:) = (0.d0,0.d0)
do j = 1, ny
do i = 1, nx
phi = (i - 1) * deltax
@ -852,8 +817,8 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
!
do j = 1, ny
do i = 1, nx
eig = exp ( (0.d0,1.d0) * 2.d0 * pi * (r (1, i, j) * g (1, ig) &
+ r (2, i, j) * g (2, ig) + r (3, i, j) * g (3, ig) ) )
eig = exp ( (0.d0,1.d0) * 2.d0 * pi * &
( r(1,i,j)*g(1,ig) + r(2,i,j)*g(2,ig) + r(3,i,j)*g(3,ig) ) )
carica (i, j) = carica (i, j) + vgc (ig) * eig
enddo
enddo
@ -861,8 +826,8 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
!
! Here we check the value of the resulting charge
!
rhomin = 1.0d10
rhomax = - 1.0d10
rhomin = 1.0d10
rhomax = -1.0d10
rhoim = 0.d0
do i = 1, nx
@ -875,8 +840,7 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
enddo
rhoim = rhoim / nx / ny
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, &
rhomax, rhoim
print '(5x,"Min, Max, imaginary charge: ",3f12.6)', rhomin, rhomax, rhoim
!
! and we print the charge on output
!
@ -886,15 +850,14 @@ subroutine plot_2ds (nx, ny, x0, ngm, g, vgc, output_format, &
!
write (ounit, '(2i8)') nx, ny
do i = 1, nx
write (ounit, '(e25.14)') (DREAL (carica (i, j) ) , j = 1, ny)
write (ounit, '(e25.14)') ( DREAL(carica(i,j)), j = 1, ny )
enddo
elseif (output_format.eq.1) then
!
! contour.x format
!
write (ounit, '(3i5,2e25.14)') nx, ny, 1, deltax, deltay
write (ounit, '(4e25.14)') ( (DREAL (carica (i, j) ) , j = 1, &
ny) , i = 1, nx)
write (ounit, '(4e25.14)') ( ( DREAL(carica(i,j)), j = 1, ny ), i = 1, nx )
else
call errore ('plot_2ds', 'not implemented plot', 1)
@ -921,8 +884,8 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
! output unit
character(len=3) :: atm(*)
real(kind=DP) :: alat, tau (3, nat), at (3, 3), g (3, ngm), e (3, 3), &
x0 (3), m1, m2, m3, dipol(0:3)
real(kind=DP) :: alat, tau(3,nat), at(3,3), g(3,ngm), e(3,3), x0(3), &
m1, m2, m3, dipol(0:3)
! lattice parameter
! atomic positions
! lattice vectors
@ -961,16 +924,16 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
! These factors are calculated and stored in order to save CPU time
!
do i = 1, nx
eigx (i) = exp((0.d0, 1.d0)*2.d0*pi*((i - 1) * deltax * &
(e(1,1)*g(1,ig)+e(2,1)*g(2,ig)+e(3,1)*g(3,ig)) &
+(x0(1)*g(1,ig)+x0(2)*g(2,ig)+x0(3)*g(3,ig) ) ) )
eigx (i) = exp( (0.d0,1.d0) * 2.d0 * pi * ( (i-1) * deltax * &
(e(1,1)*g(1,ig)+e(2,1)*g(2,ig)+e(3,1)*g(3,ig)) + &
( x0(1)*g(1,ig)+ x0(2)*g(2,ig)+ x0(3)*g(3,ig)) ) )
enddo
do j = 1, ny
eigy (j) = exp ( (0.d0,1.d0)*2.d0*pi*(j-1)*deltay* &
(e(1,2)*g(1,ig)+e(2,2)*g(2,ig)+e(3,2)*g(3,ig) ) )
eigy (j) = exp( (0.d0,1.d0) * 2.d0 * pi * (j-1) * deltay * &
(e(1,2)*g(1,ig)+e(2,2)*g(2,ig)+e(3,2)*g(3,ig)) )
enddo
do k = 1, nz
eigz (k) = exp ( (0.d0,1.d0)*2.d0*pi*(k-1)*deltaz* &
eigz (k) = exp( (0.d0,1.d0) * 2.d0 * pi * (k-1) * deltaz * &
(e(1,3)*g(1,ig)+e(2,3)*g(2,ig)+e(3,3)*g(3,ig)) )
enddo
do k = 1, nz
@ -1037,8 +1000,8 @@ subroutine plot_3d (alat, at, nat, tau, atm, ityp, ngm, g, vgc, &
do ipol=1,3
dipol(ipol)=dipol(ipol) / suma * omega * alat
enddo
print '(/5x,"Min, Max, Total, Abs charge: ",4f10.6)', rhomin, &
rhomax, rhotot, rhoabs
print '(/5x,"Min, Max, Total, Abs charge: ",4f10.6)', rhomin, rhomax, &
rhotot, rhoabs
if (output_format.eq.4) then
!
@ -1075,7 +1038,7 @@ subroutine plot_fast (alat, at, nat, tau, atm, ityp,&
!
use parameters, only : DP
implicit none
integer :: nat, ityp (nat), nrx1, nrx2, nrx3, nr1, nr2, nr3, &
integer :: nat, ityp(nat), nrx1, nrx2, nrx3, nr1, nr2, nr3, &
output_format, ounit
character(len=3) :: atm(*)

34
PP/elf.f90
View File

@ -55,9 +55,8 @@ subroutine do_elf (elf)
do is = 2, nspin
call DAXPY (nrxx, 1.d0, rho (1, is), 1, rho (1, 1), 1)
enddo
call setv (2 * nrxx, 0d0, aux, 1)
call setv (nrxx, 0d0, kkin, 1)
aux(:) = (0.d0,0.d0)
kkin(:) = 0.d0
!
! Calculates local kinetic energy, stored in kkin
!
@ -75,19 +74,16 @@ subroutine do_elf (elf)
do ibnd = 1, nbnd
do j = 1, 3
call setv (2 * nrxx, 0d0, aux, 1)
aux(:) = (0.d0,0.d0)
w1 = wg (ibnd, ik) / omega
do i = 1, npw
gv (j) = (xk (j, ik) + g (j, igk (i) ) ) * tpiba
aux (nl (igk (i) ) ) = cmplx (0d0, gv (j) ) * evc (i, ibnd)
enddo
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do i = 1, nrxx
kkin (i) = kkin (i) + w1 * (real (aux (i) ) **2 + DIMAG (aux (i) ) &
**2)
kkin(i) = kkin(i) + w1 * (real(aux(i))**2 + DIMAG(aux(i))**2)
enddo
! j
@ -103,47 +99,43 @@ subroutine do_elf (elf)
! reduce local kinetic energy across pools
!
call poolreduce (nrxx, kkin)
call psymrho (kkin, nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s, &
ftau)
call psymrho (kkin, nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s, ftau)
#else
call symrho (kkin, nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s, ftau)
call symrho (kkin, nrx1, nrx2, nrx3, nr1, nr2, nr3, nsym, s, ftau)
#endif
!
! Calculates the bosonic kinetic density, stored in tbos
! aux --> charge density in Fourier space
! aux2 --> iG * rho(G)
!
call setv (nrxx, 0d0, tbos, 1)
call setv (2 * nrxx, 0d0, aux, 1)
tbos(:) = 0.d0
aux(:) = (0.d0,0.d0)
call DCOPY (nrxx, rho, 1, aux, 2)
call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
do j = 1, 3
call setv (2 * nrxx, 0d0, aux2, 1)
aux2(:) = (0.d0,0.d0)
do i = 1, ngm
aux2 (nl (i) ) = aux (nl (i) ) * cmplx (0.0d0, g (j, i) * tpiba)
aux2(nl(i)) = aux(nl(i)) * cmplx (0.0d0, g(j,i)*tpiba)
enddo
call cft3 (aux2, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
do i = 1, nrxx
tbos (i) = tbos (i) + real (aux2 (i) ) **2
tbos (i) = tbos (i) + real(aux2(i))**2
enddo
enddo
!
! Calculates ELF
!
fac = 5.d0 / (3.d0 * (3.d0 * pi**2) ** (2.d0 / 3.d0) )
call setv (nrxx, 0d0, elf, 1)
elf(:) = 0.d0
do i = 1, nrxx
arho = abs (rho (i, 1) )
if (arho.gt.1.d-30) then
d = fac / (rho (i, 1) ** (5d0 / 3d0) ) * (kkin (i) - 0.25d0 * &
tbos (i) / rho (i, 1) )
d = fac / rho(i,1)**(5d0/3d0) * (kkin(i)-0.25d0*tbos(i)/rho(i,1))
elf (i) = 1.0d0 / (1.0d0 + d**2)
endif
enddo
deallocate (aux)
deallocate (aux2)

7
PP/ggen1d.f90
View File

@ -37,13 +37,12 @@ subroutine ggen1d (ngm1d, g1d, gg1d, ig1dto3d, nl1d, igtongl1d)
parameter (eps = 1.d-12)
call setv (nr3 * 3, 0.d0, g1d, 1)
call setv (nr3, 0.d0, gg1d, 1)
g1d(:,:) = 0.d0
gg1d(:) = 0.d0
ig1d = 0
do ig = 1, ngm
if ( (abs (g (1, ig) ) .lt.eps) .and. (abs (g (2, ig) ) .lt.eps) ) &
then
if ( (abs(g(1,ig)).lt.eps) .and. (abs(g(2,ig)) .lt.eps) ) then
!
! a vector of the 1D grid has been found
!

8
PP/local_dos.f90
View File

@ -51,9 +51,9 @@ subroutine local_dos (iflag, lsign, kpoint, kband, emin, emax, dos)
logical :: lgamma
external w0gauss, w1gauss
!
call setv (nrxx * nspin, 0.d0, rho, 1)
call setv (nrxx, 0.d0, dos, 1)
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum, 1)
rho(:,:) = 0.d0
dos(:) = 0.d0
becsum(:,:,:) = 0.d0
!
! calculate the correct weights
!
@ -105,7 +105,7 @@ subroutine local_dos (iflag, lsign, kpoint, kband, emin, emax, dos)
!
do ibnd = 1, nbnd
if (ibnd.eq.kband.or.iflag.ne.0) then
call setv (2 * nrxxs, 0.d0, psic, 1)
psic(1:nrxxs) = (0.d0,0.d0)
do ig = 1, npw
psic (nls (igk (ig) ) ) = evc (ig, ibnd)
enddo

22
PP/local_dos1d.f90
View File

@ -61,9 +61,9 @@ subroutine local_dos1d (ik, kband, plan)
allocate (prho(nrxx))
allocate (aux(nrxx))
call setv (nrxx, 0.d0, aux, 1)
aux(:) = 0.d0
becsum(:,:,:) = 0.d0
call setv ( (nhm * (nhm + 1) ) / 2 * nat * nspin, 0.d0, becsum, 1)
wg (kband, ik) = 1.d0
!
!
@ -71,17 +71,15 @@ subroutine local_dos1d (ik, kband, plan)
! mesh
!
call setv (2 * nrxxs, 0.d0, psic, 1)
psic(1:nrxxs) = (0.d0,0.d0)
do ig = 1, npw
psic (nls (igk (ig) ) ) = evc (ig, kband)
enddo
call cft3s (psic, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, 2)
w1 = wg (kband, ik) / omega
do ir = 1, nrxxs
aux (ir) = aux (ir) + w1 * (real (psic (ir) ) **2 + DIMAG (psic ( &
ir) ) **2)
aux(ir) = aux(ir) + w1 * (real(psic(ir))**2 + DIMAG(psic(ir))**2)
enddo
!
! If we have a US pseudopotential we compute here the sumbec term
@ -97,15 +95,15 @@ subroutine local_dos1d (ik, kband, plan)
ijh = 1
do ih = 1, nh (np)
ikb = ijkb0 + ih
becsum (ijh, na, current_spin) = becsum (ijh, na, &
current_spin) + w1 * real (conjg (becp (ikb, ibnd) ) &
* becp (ikb, ibnd) )
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + w1 * &
real ( conjg(becp(ikb,ibnd)) * becp(ikb,ibnd) )
ijh = ijh + 1
do jh = ih + 1, nh (np)
jkb = ijkb0 + jh
becsum (ijh, na, current_spin) = becsum (ijh, na, &
current_spin) + w1 * 2.d0 * real (conjg (becp (ikb, ibnd) ) &
* becp (jkb, ibnd) )
becsum(ijh,na,current_spin) = &
becsum(ijh,na,current_spin) + w1 * 2.d0 * &
real( conjg(becp(ikb,ibnd)) * becp(jkb,ibnd) )
ijh = ijh + 1
enddo
enddo

13
PP/plan_avg.f90
View File

@ -59,15 +59,14 @@ subroutine plan_avg (averag, plan, ninter)
! Compute the number of planes and the coordinates on the mesh of th
! points which define each plane
!
call setv (nat, 0.d0, avg, 1)
avg(:) = 0.d0
ninter = 1
z1 (ninter) = tau (3, 1)
avg (ninter) = tau (3, 1)
ntau (ninter) = 1
do na = 2, nat
do iin = 1, ninter
if (abs (mod (z1 (iin) - tau (3, na), celldm (3) ) ) .lt.sp_min) &
then
if (abs (mod (z1(iin)-tau(3,na), celldm(3)) ) .lt. sp_min) then
avg (iin) = avg (iin) + tau (3, na)
ntau (iin) = ntau (iin) + 1
goto 100
@ -112,8 +111,8 @@ subroutine plan_avg (averag, plan, ninter)
!
! for each state compute the planar average
!
call setv (nat * nbnd * nkstot, 0.d0, averag, 1)
call setv (nr3 * nbnd * nkstot, 0.d0, plan, 1)
averag(:,:,:) = 0.d0
plan(:,:,:) = 0.d0
do ik = 1, nks
if (lsda) current_spin = isk (ik)
call gk_sort (xk (1, ik), ngm, g, ecutwfc / tpiba2, npw, igk, g2kin)
@ -136,8 +135,7 @@ subroutine plan_avg (averag, plan, ninter)
sum = averag (1, ibnd, ik)
do iin = 2, ninter
do ir = i1 (iin - 1), i1 (iin) - 1
averag (iin, ibnd, ik) = averag (iin, ibnd, ik) + plan (ir, ibnd, &
ik)
averag(iin,ibnd,ik) = averag(iin,ibnd,ik) + plan(ir,ibnd,ik)
enddo
averag (iin, ibnd, ik) = averag (iin, ibnd, ik) * zdim / nr3
sum = sum + averag (iin, ibnd, ik)
@ -147,7 +145,6 @@ subroutine plan_avg (averag, plan, ninter)
#ifdef __PARA
call poolrecover (plan, nr3 * nbnd, nkstot, nks)
call poolrecover (averag, nat * nbnd, nkstot, nks)
call poolrecover (xk, 3, nkstot, nks)
#endif
return

26
PW/addusdens.f90
View File

@ -28,7 +28,7 @@ subroutine addusdens
! the spherical harmonics
complex(kind=DP) :: skk
complex(kind=DP), allocatable :: qg (:), aux (:,:)
complex(kind=DP), allocatable :: aux (:,:)
! work space for FFT
! work space for rho(G,nspin)
@ -50,13 +50,22 @@ subroutine addusdens
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) * &
@ -65,6 +74,9 @@ subroutine addusdens
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
@ -77,16 +89,12 @@ subroutine addusdens
!
! convert aux to real space and add to the charge density
!
allocate (qg( nrxx))
do is = 1, nspin
qg(:) = (0.d0, 0.d0)
do ig = 1, ngm
qg (nl (ig) ) = aux (ig, is)
enddo
call cft3 (qg, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
rho (:, is) = rho (:, is) + DREAL (qg (:) )
psic(:) = (0.d0, 0.d0)
psic( nl(:) ) = aux(:,is)
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, 1)
call DAXPY (nrxx, 1.d0, psic, 2, rho(1,is), 1)
enddo
deallocate (qg)
deallocate (aux)
call stop_clock ('addusdens')

48
PW/newd.f90
View File

@ -20,7 +20,7 @@ subroutine newd
implicit none
integer :: ig, nt, ih, jh, na, is
! counters on g vectors, atom type, beta functions x 2, atoms, spin
complex(kind=DP), allocatable :: aux (:,:), vg (:)
complex(kind=DP), allocatable :: aux (:,:), qgm_na (:)
! work space
real(kind=DP), allocatable :: ylmk0 (:,:), qmod (:)
! spherical harmonics, modulus of G
@ -49,7 +49,7 @@ subroutine newd
fact = 1.d0
end if
call start_clock ('newd')
allocate ( aux(ngm,nspin), vg(nrxx), qmod(ngm), ylmk0(ngm, lqx*lqx) )
allocate ( aux(ngm,nspin), qgm_na(ngm), qmod(ngm), ylmk0(ngm, lqx*lqx) )
!
deeq(:,:,:,:) = 0.d0
!
@ -64,10 +64,10 @@ subroutine newd
call start_clock ('newd:fftvg')
#endif
do is = 1, nspin
vg (:) = vltot (:) + vr (:, is)
call cft3 (vg, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
psic (:) = vltot (:) + vr (:, is)
call cft3 (psic, nr1, nr2, nr3, nrx1, nrx2, nrx3, - 1)
do ig = 1, ngm
aux (ig, is) = vg (nl (ig) )
aux (ig, is) = psic (nl (ig) )
enddo
enddo
#ifdef DEBUG_NEWD
@ -84,48 +84,50 @@ subroutine newd
#ifdef DEBUG_NEWD
call start_clock ('newd:qvan2')
#endif
!
! The Q(r) for this atomic species without structure factor
!
call qvan2 (ngm, ih, jh, nt, qmod, qgm, ylmk0)
#ifdef DEBUG_NEWD
call stop_clock ('newd:qvan2')
#endif
do na = 1, nat
if (ityp (na) .eq.nt) then
!
! The product of the potential and the structure factor
!
do is = 1, nspin
#ifdef DEBUG_NEWD
call start_clock ('newd:int1')
#endif
do ig = 1, ngm
vg (ig) = aux(ig, is) * conjg(eigts1 (ig1(ig), na) &
* eigts2 (ig2(ig), na) &
* eigts3 (ig3(ig), na) )
enddo
!
! The Q(r) for this specific atom
!
do ig = 1, ngm
qgm_na (ig) = qgm(ig) * eigts1 (ig1(ig), na) &
* eigts2 (ig2(ig), na) &
* eigts3 (ig3(ig), na)
enddo
#ifdef DEBUG_NEWD
call stop_clock ('newd:int1')
#endif
!
! and the product with the Q functions
!
#ifdef DEBUG_NEWD
call start_clock ('newd:int2')
#endif
!
! and the product with the Q functions
!
do is = 1, nspin
deeq (ih, jh, na, is) = fact * omega * &
DDOT (2 * ngm, vg, 1, qgm, 1)
DDOT (2 * ngm, aux(1,is), 1, qgm_na, 1)
if (gamma_only .and. gstart==2) &
deeq (ih, jh, na, is) = &
deeq (ih, jh, na, is) - omega*real(vg(1)*qgm(1))
deeq (ih, jh, na, is) = deeq (ih, jh, na, is) - &
omega*real(aux(1,is)*qgm_na(1))
enddo
#ifdef DEBUG_NEWD
call stop_clock ('newd:int2')
#endif
enddo
endif
enddo
enddo
enddo
endif
enddo
#ifdef __PARA
call reduce (nhm * nhm * nat * nspin, deeq)
@ -150,7 +152,7 @@ subroutine newd
! end do
enddo
deallocate ( aux, vg, qmod, ylmk0 )
deallocate ( aux, qgm_na, qmod, ylmk0 )
call stop_clock ('newd')
return

9
PW/print_clock_pw.f90
View File

@ -29,7 +29,8 @@ subroutine print_clock_pw
call print_clock ('v_of_rho')
call print_clock ('newd')
#ifdef DEBUG_NEWD
write (*,*) nhm*(nhm+1)/2, nbrx*(nbrx+1)/2*lqx
write (*,*) "nhm*(nhm+1)/2 = ", nhm*(nhm+1)/2, nhm
write (*,*) "nbrx*(nbrx+1)/2*lqx = ", nbrx*(nbrx+1)/2*lqx, nbrx,lqx
call print_clock ('newd:fftvg')
call print_clock ('newd:qvan2')
call print_clock ('newd:int1')
@ -55,6 +56,12 @@ subroutine print_clock_pw
call print_clock ('sumbec')
call print_clock ('addusdens')
#ifdef DEBUG_ADDUSDENS
call print_clock ('addus:qvan2')
call print_clock ('addus:strf')
call print_clock ('addus:aux2')
call print_clock ('addus:aux')
#endif
write (6, * )
if (isolve.eq.0) then
call print_clock ('cegterg')

11
PW/struct_fact.f90
View File

@ -35,8 +35,10 @@ subroutine struc_fact (nat, tau, ntyp, ityp, ngm, g, bg, nr1, nr2, &
! input: the positions of the atoms in the c
! input: the coordinates of the g vectors
complex(kind=DP) :: strf (ngm, ntyp), eigts1 ( - nr1:nr1, nat), &
eigts2 ( - nr2:nr2, nat), eigts3 ( - nr3:nr3, nat)
complex(kind=DP) :: strf (ngm, ntyp), &
eigts1 ( -nr1:nr1, nat), &
eigts2 ( -nr2:nr2, nat), &
eigts3 ( -nr3:nr3, nat)
! output: the structure factor
!
! output: the phases e^{-iG\tau_s}
@ -74,8 +76,9 @@ subroutine struc_fact (nat, tau, ntyp, ityp, ngm, g, bg, nr1, nr2, &
do na = 1, nat
do ipol = 1, 3
bgtau (ipol) = bg (1, ipol) * tau (1, na) + bg (2, ipol) * tau (2, &
na) + bg (3, ipol) * tau (3, na)
bgtau (ipol) = bg (1, ipol) * tau (1, na) + &
bg (2, ipol) * tau (2, na) + &
bg (3, ipol) * tau (3, na)
enddo
do n1 = - nr1, nr1
arg = tpi * n1 * bgtau (1)