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

20 Jan 2003 Added dielectric tensor calculation with USPP 

(experimental) - AdC


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@9 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2003-01-21 15:25:13 +00:00
parent 78ddef1f56
commit d89b07a22a
15 changed files with 236 additions and 154 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
PH/Makefile
View File

@ -9,9 +9,11 @@ include ../make.sys
PHOBJS = phcom.o \
addcore.o \
adddvscf.o \
adddvepsi_us.o \
addnlcc.o \
addusdbec.o \
addusddens.o \
addusddense.o \
addusdynmat.o \
addusldos.o \
add_zstar_ue.o \
@ -29,6 +31,7 @@ compute_becsum.o \
compute_drhous.o \
compute_dvloc.o \
compute_nldyn.o \
compute_qdipol.o \
compute_weight.o \
d2ionq.o \
davcio_drho.o \

14
PH/addusdynmat.f90
View File

@ -59,8 +59,8 @@ subroutine addusdynmat (dynwrk)
nu_j = 3 * (na - 1) + jpol
do is = 1, nspin
do ijh = 1, dim
dynwrk (nu_i, nu_j) = dynwrk (nu_i, nu_j) + int4 (ijh, ipol, &
jpol, na, is) * becsum (ijh, na, is)
dynwrk(nu_i, nu_j)=dynwrk(nu_i, nu_j)+ &
int4(ijh,ipol,jpol,na,is)*becsum (ijh, na, is)
enddo
enddo
enddo
@ -76,8 +76,8 @@ subroutine addusdynmat (dynwrk)
do jh = ih, nh (np)
ijh = ijh + 1
do is = 1, nspin
term (ipol, jpol) = term (ipol, jpol) + conjg (int1 (ih, jh, &
ipol, na, is) ) * alphasum (ijh, jpol, na, is)
term(ipol,jpol) = term(ipol,jpol)+ &
conjg(int1(ih,jh,ipol,na,is))*alphasum(ijh,jpol,na,is)
enddo
enddo
enddo
@ -107,9 +107,9 @@ subroutine addusdynmat (dynwrk)
do jh = ih, nh (np)
ijh = ijh + 1
do is = 1, nspin
dyn1 (nu_i, nu_j) = dyn1 (nu_i, nu_j) + conjg (int2 (ih, jh, &
ipol, nb, na) ) * alphasum (ijh, jpol, na, is) + int5 (ijh, &
ipol, jpol, nb, na) * becsum (ijh, na, is)
dyn1(nu_i,nu_j)=dyn1(nu_i,nu_j)+conjg(int2 (ih, jh, &
ipol,nb,na))*alphasum(ijh,jpol,na,is) &
+ int5 (ijh,ipol,jpol,nb,na)*becsum(ijh,na,is)
enddo
enddo
enddo

1
PH/allocate_phq.f90
View File

@ -64,6 +64,7 @@ subroutine allocate_phq
call mallocate(int3 , nhm , nhm , 3 , nat , nspin)
call mallocate(int4 , nhm * (nhm + 1)/2, 3 , 3 , nat, nspin)
call mallocate(int5 , nhm * (nhm + 1)/2 , 3 , 3 , nat , nat)
call mallocate(dpqq, nhm, nhm, 3, ntyp)
call mallocate(alphasum , nhm * (nhm + 1)/2 , 3 , nat , nspin)
endif
call mallocate(alphap , nkb , nbnd , 3 , nksq)

35
PH/dielec.f90
View File

@ -31,24 +31,24 @@ complex(kind=DP) :: ZDOTC
call setv (9, 0.0d0, epsilon, 1)
if (nksq.gt.1) rewind (unit = iunigk)
do ik = 1, nksq
if (nksq.gt.1) read (iunigk) npw, igk
weight = wk (ik)
w = fpi * weight / omega
do ipol = 1, 3
nrec = (ipol - 1) * nksq + ik
call davcio (dvpsi, lrbar, iubar, nrec, - 1)
do jpol = 1, 3
nrec = (jpol - 1) * nksq + ik
call davcio (dpsi, lrdwf, iudwf, nrec, - 1)
do ibnd = 1, nbnd_occ (ik)
if (nksq.gt.1) read (iunigk) npw, igk
weight = wk (ik)
w = fpi * weight / omega
do ipol = 1, 3
nrec = (ipol - 1) * nksq + ik
call davcio (dvpsi, lrbar, iubar, nrec, - 1)
do jpol = 1, 3
nrec = (jpol - 1) * nksq + ik
call davcio (dpsi, lrdwf, iudwf, nrec, - 1)
do ibnd = 1, nbnd_occ (ik)
!
! this is the real part of <DeltaV*psi(E)|DeltaPsi(E)>
!
epsilon (ipol, jpol) = epsilon (ipol, jpol) + 4.0 * w * DREAL ( &
ZDOTC (npw, dvpsi (1, ibnd), 1, dpsi (1, ibnd), 1) )
enddo
enddo
enddo
epsilon(ipol,jpol)=epsilon(ipol,jpol)-4.d0*w*DREAL( &
ZDOTC (npw, dvpsi (1, ibnd), 1, dpsi (1, ibnd), 1) )
enddo
enddo
enddo
enddo
#ifdef PARA
call reduce (9, epsilon)
@ -75,14 +75,13 @@ call trntns (epsilon, at, bg, 1)
! add the diagonal part
!
do ipol = 1, 3
epsilon (ipol, ipol) = epsilon (ipol, ipol) + 1.d0
epsilon (ipol, ipol) = epsilon (ipol, ipol) + 1.d0
enddo
!
! and print the result
!
write (6, '(/,10x,"Dielectric constant in cartesian axis ",/)')
write (6, '(10x,"(",3f18.9," )")') ( (epsilon (ipol, jpol) , &
ipol = 1, 3) , jpol = 1, 3)
write (6, '(10x,"(",3f18.9," )")') ((epsilon(ipol,jpol), ipol=1,3), jpol=1,3)
return
end subroutine dielec

85
PH/dvpsi_e.f90
View File

@ -20,11 +20,12 @@ subroutine dvpsi_e (kpoint, ipol)
use pwcom
use allocate
use parameters, only : DP
use becmod
use phcom
implicit none
!
integer :: ipol, ig, na, ibnd, jbnd, ikb, jkb, nt, lter, kpoint
integer :: ipol, ig, na, ibnd, jbnd, ikb, jkb, nt, lter, kpoint, &
ih, jh, ijkb0
! input: the polarization
! counter on G vectors
! counter on atoms
@ -44,18 +45,21 @@ subroutine dvpsi_e (kpoint, ipol)
logical :: conv_root
! true if convergence has been achieved
complex(kind=DP), pointer :: ps (:,:), dvkb (:,:), dvkb1 (:,:), work (:)
complex(kind=DP), pointer :: ps (:,:), dvkb (:,:), dvkb1 (:,:), work (:,:), &
becp2(:,:), spsi(:,:)
complex(kind=DP) :: ZDOTC
! the scalar products
external ch_psi_all, cg_psi
!
!
call start_clock ('dvpsi_e')
call mallocate(work , npwx)
call mallocate(work , npwx, nkb)
call mallocate(gk , 3, npwx)
call mallocate(h_diag, npwx , nbnd)
call mallocate(ps , 2 , nbnd)
call mallocate(spsi , npwx, nbnd)
call mallocate(eprec , nbnd)
if (nkb.gt.0) call mallocate(becp2 , nkb, nbnd)
if (nkb.gt.0) call mallocate(dvkb , npwx , nkb)
if (nkb.gt.0) call mallocate(dvkb1, npwx , nkb)
call setv (2 * npwx * nkb, 0.d0, dvkb, 1)
@ -65,7 +69,7 @@ subroutine dvpsi_e (kpoint, ipol)
gk (1, ig) = (xk (1, kpoint) + g (1, igk (ig) ) ) * tpiba
gk (2, ig) = (xk (2, kpoint) + g (2, igk (ig) ) ) * tpiba
gk (3, ig) = (xk (3, kpoint) + g (3, igk (ig) ) ) * tpiba
g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2
! g2kin (ig) = gk (1, ig) **2 + gk (2, ig) **2 + gk (3, ig) **2
enddo
!
! this is the kinetic contribution to [H,x]: -2i (k+G)_ipol * psi
@ -74,8 +78,8 @@ subroutine dvpsi_e (kpoint, ipol)
do ig = 1, npw
dpsi (ig, ibnd) = (at(1, ipol) * gk(1, ig) + &
at(2, ipol) * gk(2, ig) + &
at(3, ipol) * gk(3, ig) ) * (0.d0,-2.d0) &
* evc (ig, ibnd)
at(3, ipol) * gk(3, ig) ) &
*(0.d0,-2.d0)*evc (ig, ibnd)
enddo
enddo
!
@ -96,31 +100,47 @@ subroutine dvpsi_e (kpoint, ipol)
enddo
jkb = 0
do nt = 1, nat
do nt = 1, ntyp
do na = 1, nat
if (nt.eq.ityp (na)) then
do ikb = 1, nh (nt)
jkb = jkb + 1
do ig = 1, npw
work (ig) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * &
work (ig,jkb) = dvkb1 (ig, jkb) + dvkb (ig, jkb) * &
(at (1, ipol) * gk (1, ig) + &
at (2, ipol) * gk (2, ig) + &
at (3, ipol) * gk (3, ig) )
enddo
do ibnd = 1, nbnd_occ (kpoint)
ps (1, ibnd) = ZDOTC(npw,work,1,evc(1,ibnd),1) * &
(0.d0,-1.d0) * dvan(ikb,ikb,nt)
ps (2, ibnd) = ZDOTC(npw,vkb(1,jkb),1,evc(1,ibnd),1) * &
(0.d0,-1.d0) * dvan(ikb,ikb,nt)
enddo
endif
enddo
enddo
call ccalbec (nkb, npwx, npw, nbnd, becp2, work, evc)
ijkb0 = 0
do nt = 1, ntyp
do na = 1, nat
if (nt.eq.ityp (na)) then
do ih = 1, nh (nt)
ikb = ijkb0 + ih
ps(:,:)=(0.d0,0.d0)
do jh = 1, nh (nt)
jkb = ijkb0 + jh
do ibnd = 1, nbnd_occ (kpoint)
ps (1, ibnd) = ps(1,ibnd)+ becp2(jkb,ibnd)* &
(0.d0,-1.d0)*(deeq(ih,jh,na,current_spin) &
-et(ibnd,kpoint)*qq(ih,jh,nt))
ps (2, ibnd) = ps(2,ibnd) +becp1(jkb,ibnd,kpoint) * &
(0.d0,-1.d0)*(deeq(ih,jh,na,current_spin)&
-et(ibnd,kpoint)*qq(ih,jh,nt))
enddo
enddo
#ifdef PARA
call reduce (4 * nbnd, ps)
#endif
do ibnd = 1, nbnd_occ (kpoint)
call ZAXPY(npw,ps(1,ibnd),vkb(1,jkb),1,dpsi(1,ibnd),1)
call ZAXPY(npw,ps(2,ibnd),work,1,dpsi(1,ibnd),1)
call ZAXPY(npw,ps(1,ibnd),vkb(1,ikb),1,dpsi(1,ibnd),1)
call ZAXPY(npw,ps(2,ibnd),work(1,ikb),1,dpsi(1,ibnd),1)
enddo
enddo
ijkb0=ijkb0+nh(nt)
end if
end do
end do
@ -139,17 +159,20 @@ subroutine dvpsi_e (kpoint, ipol)
do jbnd = 1, nbnd_occ (kpoint)
call ZAXPY (npw, ps (1, jbnd), evc (1, jbnd), 1, work, 1)
enddo
call DAXPY (2 * npw, 1.0d0, work, 1, dpsi (1, ibnd), 1)
call ccalbec (nkb, npwx, npw, 1, becp, vkb, work)
call s_psi (npwx, npw, 1, work, spsi)
call DAXPY (2 * npw, 1.0d0, spsi, 1, dpsi (1, ibnd), 1)
enddo
!
! dpsi contains now P_c [H,x] psi_v for the three crystal polarizations
! Now solve the linear systems (H-e_v)*(x*psi_v) = P_c [H,x]*psi_v
! dpsi contains now P^+_c [H-eS,x] psi_v for the three crystal
! polarizations
! Now solve the linear systems (H-e_vS)*P_c(x*psi_v)=P_c^+ [H-e_vS,x]*psi_v
!
thresh = 1.d-5
do ibnd = 1, nbnd_occ (kpoint)
conv_root = .true.
do ig = 1, npwq
work (ig) = g2kin (ig) * evc (ig, ibnd)
work (ig,1) = g2kin (ig) * evc (ig, ibnd)
enddo
eprec (ibnd) = 1.35d0 * ZDOTC (npwq, evc (1, ibnd), 1, work, 1)
enddo
@ -165,16 +188,30 @@ subroutine dvpsi_e (kpoint, ipol)
call cgsolve_all (ch_psi_all, cg_psi, et (1, kpoint), dpsi, dvpsi, &
h_diag, npwx, npw, thresh, kpoint, lter, conv_root, anorm, &
nbnd_occ (kpoint) )
!
if (.not.conv_root) write (6, '(5x,"kpoint",i4," ibnd",i4, &
& " linter: root not converged ",e10.3)') &
& kpoint, ibnd, anorm
#ifdef FLUSH
call flush (6)
#endif
!
! In the US case we obtain P_c x |psi>, but we need P_c^+ x | psi>,
! therefore we apply S again, and then subtract the additional term
! furthermore we add the term due to dipole of the augmentation charges.
!
if (okvan) then
call ccalbec(nkb,npwx,npw,nbnd,becp,vkb,dvpsi)
call s_psi(npwx,npw,nbnd,dvpsi,spsi)
call DCOPY(2*npwx*nbnd,spsi,1,dvpsi,1)
call adddvepsi_us(becp2,ipol,kpoint)
endif
if (nkb.gt.0) call mfree (dvkb1)
if (nkb.gt.0) call mfree (dvkb)
if (nkb.gt.0) call mfree (becp2)
call mfree (eprec)
call mfree (spsi)
call mfree (ps)
call mfree (h_diag)
call mfree (gk)

5
PH/dynmatrix.f90
View File

@ -109,9 +109,8 @@ call q2qstar_ph (dyn, at, bg, nat, nsym, s, invs, irt, rtau, nq, &
!
! Writes (if the case) results for quantities involving electric field
!
if (epsil) call write_epsilon_and_zeu (zstareu, epsilon, nat, &
iudyn)
if (zue) call sym_and_write_zue
if (epsil) call write_epsilon_and_zeu (zstareu, epsilon, nat, iudyn)
if (zue.and..not.okvan) call sym_and_write_zue
!
! Diagonalizes the dynamical matrix at q
!

17
PH/newdq.f90
View File

@ -7,7 +7,7 @@
!
!
!----------------------------------------------------------------------
subroutine newdq (dvscf, irr, npe)
subroutine newdq (dvscf, npe)
!----------------------------------------------------------------------
!
! This routine computes the contribution of the selfconsistent
@ -30,9 +30,7 @@ subroutine newdq (dvscf, irr, npe)
complex(kind=DP) :: dvscf (nrxx, nspin, npe)
! input: the change of the self
!consistent pot.
integer :: irr
! input: the irreducible representation
! consistent pot.
!
! And the local variables
!
@ -92,7 +90,7 @@ subroutine newdq (dvscf, irr, npe)
! integrate the change of the self consistent potential and
! the Q functions
!
do ipert = 1, npert (irr)
do ipert = 1, npe
do is = 1, nspin
do ir = 1, nrxx
veff (ir) = dvscf (ir, is, ipert)
@ -111,12 +109,11 @@ subroutine newdq (dvscf, irr, 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, &
int3(ih,jh,ipert,na,is)=omega*ZDOTC(ngm,aux1,1, &
aux2 (1, is), 1)
enddo
!
@ -135,7 +132,7 @@ subroutine newdq (dvscf, irr, npe)
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

8
PH/phcom.f90
View File

@ -147,14 +147,15 @@ end module gc_ph
!
module phus
use parameters, only : DP
real (kind=DP), pointer :: alphasum (:,:,:,:) ! nhm*(nhm+1)/2,3,nat,nspin)
! used to compute modes
real (kind=DP), pointer :: alphasum (:,:,:,:),&! nhm*(nhm+1)/2,3,nat,nspin)
! used to compute modes
dpqq(:,:,:,:) ! dipole moment of each Q
complex (kind=DP), pointer :: int1 (:,:,:,:,:),&! nhm, nhm, 3, nat, nspin),&
int2 (:,:,:,:,:),& ! nhm, nhm, 3,nat, nat),&
int3 (:,:,:,:,:),& ! nhm, nhm, 3, nat, nspin),&
int4 (:,:,:,:,:),& ! nhm*(nhm+1)/2, 3, 3, nat, nspin),&
int5 (:,:,:,:,:),& ! nhm*(nhm+1)/2, 3, 3, nat, nat),&
becp1 (:,:,:),& ! nkbtot, nbnd, nksq),&
becp1 (:,:,:), & ! nkbtot, nbnd, nksq),&
alphap (:,:,:,:) ! nkbtot, nbnd, 3, nksq)
! integrals of dQ and V_eff
! integrals of dQ and V_loc
@ -164,7 +165,6 @@ module phus
! the becq used in ch_psi
! the derivative of the bec
end module phus
!
! the variables needed for partial computation of dynamical matrix
!

3
PH/phq_init.f90
View File

@ -181,6 +181,9 @@ subroutine phq_init
call drho
endif
if (epsil.and.okvan) then
call compute_qdipol
endif
call stop_clock ('phq_init')
return
end subroutine phq_init

2
PH/phq_readin.f90
View File

@ -161,7 +161,7 @@ subroutine phq_readin
nksq = nks / 2
endif
!
if ( (epsil.or.zue) .and.okvan) then
if ( (zue) .and.okvan) then
call error ('phq_readin', 'No electric field in US case', &
- 1)
epsil = .false.

179
PH/solve_e.f90
View File

@ -10,12 +10,12 @@
subroutine solve_e
!-----------------------------------------------------------------------
!
! This routine is a driver for the solution of the linear system whic
! defines the change of the wavefunction due to the electric field.
! This routine is a driver for the solution of the linear system which
! defines the change of the wavefunction due to an electric field.
! It performs the following tasks:
! a) It computes the kinetic energy
! b) It adds the term Delta V_{SCF} | psi >
! c) It applies P_c to the known part
! c) It applies P_c^+ to the known part
! d) It calls linter to solve the linear system
! e) It computes Delta rho, Delta V_{SCF} and symmetrize them
!
@ -23,6 +23,7 @@ subroutine solve_e
use pwcom
use allocate
use parameters, only : DP
use becmod
use phcom
implicit none
@ -36,23 +37,21 @@ subroutine solve_e
! cut-off for preconditioning
! convergence limit
complex(kind=DP) , pointer :: dvscfin (:,:), dvscfins (:,:),&
dvscfout (:,:), auxg (:), aux1 (:), ps (:)
! change of the scf potential (input
! change of the scf potential (smoot
! change of the scf potential (outpu
! auxiliary space
! the psi function
! the scalar product
complex(kind=DP) :: dbecsum, ZDOTC
! dummy
! the scalar product function
complex(kind=DP) , pointer :: &
dvscfin (:,:,:), & ! change of the scf potential (input)
dvscfins (:,:,:), & ! change of the scf potential (smooth)
dvscfout (:,:,:), & ! change of the scf potential (output)
dbecsum(:,:,:,:), & ! the becsum with dpsi
auxg (:), aux1 (:), spsi(:), ps (:)
complex(kind=DP) :: ZDOTC ! the scalar product function
logical :: conv_root, exst
! true if linter is converged
! used to open the recover file
integer :: kter, ipol, ibnd, jbnd, iter, lter, ltaver, lintercall, &
ik, ig, irr, ir, nrec, nrec1, ios
ik, ig, irr, ir, is, nrec, nrec1, ios
! counter on iterations
! counter on perturbations
! counter on bands
@ -79,28 +78,32 @@ subroutine solve_e
external ch_psi_all, cg_psi
if (lsda) call error ('solve_e', ' LSDA not implemented', 1)
if (lsda) call error ('solve_e', ' LSDA not implemented', 1)
call start_clock ('solve_e')
call mallocate(dvscfin, nrxx , 3)
call mallocate(dvscfin, nrxx, nspin, 3)
if (doublegrid) then
call mallocate(dvscfins, nrxxs , 3)
call mallocate(dvscfins, nrxxs, nspin, 3)
else
dvscfins => dvscfin
endif
call mallocate(dvscfout, nrxx , 3)
call mallocate(auxg , npwx)
call mallocate(aux1 , nrxxs)
call mallocate(ps , nbnd)
call mallocate(h_diag , npwx , nbnd)
call mallocate(eprec , nbnd)
call mallocate(dvscfout, nrxx , nspin, 3)
call mallocate(dbecsum, nhm*(nhm+1)/2, nat, nspin, 3)
call mallocate(auxg, npwx)
call mallocate(aux1, nrxxs)
call mallocate(spsi, npwx)
call mallocate(ps, nbnd)
call mallocate(h_diag, npwx , nbnd)
call mallocate(eprec, nbnd)
if (iter0.ne.0) then
if (okvan) read(iunrec) int3
read (iunrec) dr2, dvscfin
close (unit = iunrec, status = 'keep')
if (doublegrid) then
do ipol = 1, 3
call cinterpolate (dvscfin (1, ipol), dvscfins (1, ipol), &
- 1)
do is=1,nspin
do ipol=1,3
call cinterpolate (dvscfin(1,is,ipol), dvscfins(1,is,ipol), -1)
enddo
enddo
endif
@ -126,10 +129,12 @@ subroutine solve_e
ltaver = 0
lintercall = 0
call setv (2 * nrxx * 3, 0.d0, dvscfout, 1)
dvscfout(:,:,:)=(0.d0,0.d0)
dbecsum(:,:,:,:)=(0.d0,0.d0)
if (nksq.gt.1) rewind (unit = iunigk)
do ik = 1, nksq
if (lsda) current_spin = isk (ik)
if (nksq.gt.1) then
read (iunigk, err = 100, iostat = ios) npw, igk
100 call error ('solve_e', 'reading igk', abs (ios) )
@ -159,12 +164,13 @@ subroutine solve_e
! At the first iteration dpsi and dvscfin are set to zero,
! [H,x]*psi_kpoint is calculated and written to file
!
call setv (2 * nbnd * npwx, 0.d0, dpsi, 1)
call setv (2 * nrxx, 0.d0, dvscfin (1, ipol), 1)
dpsi(:,:)=(0.d0,0.d0)
dvscfin(:,:,:)=(0.d0,0.d0)
call dvpsi_e (ik, ipol)
call davcio (dvpsi, lrbar, iubar, nrec, 1)
!
! starting threshold for the iterative solution of the linear sistem (li
! starting threshold for the iterative solution of the linear
! system
!
thresh = 1.d-2
else
@ -179,19 +185,18 @@ subroutine solve_e
do ibnd = 1, nbnd_occ (ik)
call setv (2 * nrxxs, 0.d0, aux1, 1)
do ig = 1, npw
aux1 (nls (igk (ig) ) ) = evc (ig, ibnd)
aux1 (nls(igk(ig)))=evc(ig,ibnd)
enddo
call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, &
+ 2)
call cft3s (aux1,nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,+2)
do ir = 1, nrxxs
aux1 (ir) = aux1 (ir) * dvscfins (ir, ipol)
aux1(ir)=aux1(ir)*dvscfins(ir,current_spin,ipol)
enddo
call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, &
- 2)
call cft3s (aux1,nr1s,nr2s,nr3s,nrx1s,nrx2s,nrx3s,-2)
do ig = 1, npwq
dvpsi (ig, ibnd) = dvpsi (ig, ibnd) - aux1 (nls (igkq (ig) ) )
dvpsi(ig,ibnd)=dvpsi(ig,ibnd)+aux1(nls(igkq(ig)))
enddo
enddo
call adddvscf(ipol,ik)
!
! starting value for delta_psi is read from iudwf
!
@ -208,45 +213,50 @@ subroutine solve_e
do ibnd = 1, nbnd_occ (ik)
call setv (2 * npwx, 0.d0, auxg, 1)
do jbnd = 1, nbnd_occ (ik)
ps (jbnd) = - ZDOTC (npwq, evc (1, jbnd), 1, dvpsi (1, ibnd), &
1)
ps(jbnd)=-ZDOTC(npw,evc(1,jbnd),1,dvpsi(1,ibnd),1)
enddo
#ifdef PARA
call reduce (2 * nbnd, ps)
#endif
do jbnd = 1, nbnd_occ (ik)
call ZAXPY (npwq, ps (jbnd), evc (1, jbnd), 1, auxg, 1)
call ZAXPY (npw, ps (jbnd), evc (1, jbnd), 1, auxg, 1)
enddo
call DAXPY (2 * npwq, 1.0d0, auxg, 1, dvpsi (1, ibnd), 1)
call ccalbec (nkb, npwx, npw, 1, becp, vkb, auxg)
call s_psi (npwx, npw, 1, auxg, spsi)
call DAXPY (2*npw, 1.0d0, spsi, 1, dvpsi (1, ibnd), 1)
enddo
!
! Here we change the sign of the known term
!
call DSCAL (2*npwx*nbnd, -1.d0, dvpsi, 1)
!
! iterative solution of the linear system (H-e)*dpsi=dvpsi
! dvpsi=-P_c+ (dvbare+dvscf)*psi , dvscf fixed.
!
do ibnd = 1, nbnd_occ (ik)
do ig = 1, npwq
do ig = 1, npw
auxg (ig) = g2kin (ig) * evc (ig, ibnd)
enddo
eprec (ibnd) = 1.35d0 * ZDOTC (npwq, evc (1, ibnd), 1, auxg, 1)
eprec (ibnd) = 1.35d0*ZDOTC(npwq,evc(1,ibnd),1,auxg,1)
enddo
#ifdef PARA
call reduce (nbnd_occ (ik), eprec)
#endif
do ibnd = 1, nbnd_occ (ik)
do ig = 1, npwq
h_diag (ig, ibnd) = 1.d0 / max (1.0d0, g2kin (ig) / eprec (ibnd) )
do ig = 1, npw
h_diag(ig,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd))
enddo
enddo
conv_root = .true.
call cgsolve_all (ch_psi_all, cg_psi, et (1, ik), dvpsi, dpsi, &
h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, nbnd_occ ( &
ik) )
call cgsolve_all (ch_psi_all,cg_psi,et(1,ik),dvpsi,dpsi, &
h_diag,npwx,npw,thresh,ik,lter,conv_root,anorm,nbnd_occ(ik) )
ltaver = ltaver + lter
lintercall = lintercall + 1
if (.not.conv_root) write (6, '(5x,"kpoint",i4," ibnd",i4, &
& " linter: root not converged ",e10.3)') ik &
if (.not.conv_root) write (6, '(5x,''kpoint'',i4,'' ibnd'',i4, &
& '' linter: root not converged '',e10.3)') ik &
&, ibnd, anorm
!
! writes delta_psi on iunit iudwf, k=kpoint,
@ -258,24 +268,37 @@ subroutine solve_e
! calculates dvscf, sum over k => dvscf_q_ipert
!
weight = wk (ik)
call incdrhoscf (dvscfout (1, ipol), weight, ik, dbecsum, 1, 1)
call incdrhoscf (dvscfout (1, current_spin, ipol), weight, &
ik, dbecsum(1,1,current_spin,ipol), 1, 1)
enddo ! on perturbation
enddo ! on k points
#ifdef PARA
call reduce (nhm * (nhm + 1) * nat * 3* nspin, dbecsum)
#endif
if (doublegrid) then
do is=1,nspin
do ipol=1,3
call cinterpolate (dvscfout(1,is,ipol),dvscfout(1,is,ipol),1)
enddo
enddo
endif
enddo
do ipol = 1, 3
if (doublegrid) call cinterpolate (dvscfout (1, ipol), dvscfout ( &
1, ipol), 1)
if (fildrho.ne.' ') call davcio_drho (dvscfout (1, ipol) , lrdrho, &
iudrho, ipol, + 1)
call dv_of_drho (0, dvscfout (1, ipol), .false.)
call addusddense (dvscfout, dbecsum)
enddo
!
! After the loop over the perturbations we have the change of the pote
! for all the modes of this representation. We symmetrize this potenti
!
#ifdef PARA
call poolreduce (2 * 3 * nrxx, dvscfout)
call poolreduce (2 * 3 * nrxx *nspin, dvscfout)
#endif
do ipol=1,3
if (fildrho.ne.' ') call davcio_drho(dvscfout(1,1,ipol),lrdrho, &
iudrho,ipol,+1)
call dv_of_drho (0, dvscfout (1, 1, ipol), .false.)
enddo
#ifdef PARA
call psyme (dvscfout)
#else
call syme (dvscfout)
@ -284,21 +307,24 @@ subroutine solve_e
! And we mix with the old potential
!
call mix_potential (2 * 3 * nrxx, dvscfout, dvscfin, alpha_mix ( &
call mix_potential (2 * 3 * nrxx *nspin, dvscfout, dvscfin, alpha_mix ( &
kter), dr2, 3 * tr2_ph, iter, nmix_ph, flmixdpot, convt)
if (doublegrid) then
do ipol = 1, 3
call cinterpolate (dvscfin (1, ipol), dvscfins (1, ipol), &
- 1)
do is=1,nspin
do ipol = 1, 3
call cinterpolate (dvscfin(1,is,ipol),dvscfins(1,is,ipol),-1)
enddo
enddo
endif
call newdq(dvscfin,3)
averlt = dfloat (ltaver) / dfloat (lintercall)
write (6, '(//,5x," iter # ",i3, &
& " av.it.: ",f5.1)') iter, averlt
write (6, '(//,5x,'' iter # '',i3, &
& '' av.it.: '',f5.1)') iter, averlt
dr2 = dr2 / 3
write (6, '(5x," thresh=",e10.3, " alpha_mix = ",f6.3, &
& " |ddv_scf|^2 = ",e10.3 )') thresh, alpha_mix (kter) , dr2
write (6, '(5x,'' thresh='',e10.3, '' alpha_mix = '',f6.3, &
& ''|ddv_scf|^2 = '',e10.3 )') thresh, alpha_mix (kter), dr2
#ifdef FLUSH
call flush (6)
#endif
@ -309,11 +335,11 @@ subroutine solve_e
write (iunrec) dyn, dyn00, epsilon, zstareu, zstarue, zstareu0, &
zstarue0
if (reduce_io) then
write (iunrec) irr, 0, convt, done_irr, comp_irr, ifat
write(iunrec) irr, 0, convt, done_irr, comp_irr, ifat
else
write (iunrec) irr, iter, convt, done_irr, comp_irr, ifat
write (iunrec) dr2, dvscfin
write(iunrec) irr, iter, convt, done_irr, comp_irr, ifat
if (okvan) write(iunrec) int3
write(iunrec) dr2, dvscfin
endif
close (unit = iunrec, status = 'keep')
@ -323,16 +349,17 @@ subroutine solve_e
enddo
155 continue
if (tcpu.gt.time_max) then
write (6, '(/,5x,"Stopping for time limit ",2f10.0)') tcpu, &
write (6, '(/,5x,''Stopping for time limit '',2f10.0)') tcpu, &
time_max
call stop_ph (.false.)
endif
call mfree (eprec)
call mfree (h_diag)
call mfree (ps)
call mfree (spsi)
call mfree (aux1)
call mfree (auxg)
call mfree (dbecsum)
call mfree (dvscfout)
if (doublegrid) call mfree (dvscfins)
call mfree (dvscfin)

10
PH/solve_linter.f90
View File

@ -355,14 +355,14 @@ subroutine solve_linter (irr, imode0, npe, drhoscf)
#endif
do ibnd = 1, nbnd_occ (ikk)
do ig = 1, npwq
h_diag (ig, ibnd) = 1.d0 / max (1.0d0, g2kin (ig) / eprec (ibnd) )
h_diag(ig,ibnd)=1.d0/max(1.0d0,g2kin(ig)/eprec(ibnd))
enddo
enddo
conv_root = .true.
call cgsolve_all (ch_psi_all, cg_psi, et (1, ikk), dvpsi, dpsi, &
h_diag, npwx, npwq, thresh, ik, lter, conv_root, anorm, nbnd_occ ( &
ikk) )
call cgsolve_all (ch_psi_all,cg_psi,et(1,ikk),dvpsi,dpsi, &
h_diag,npwx,npwq,thresh,ik,lter,conv_root,anorm,nbnd_occ(ikk))
ltaver = ltaver + lter
lintercall = lintercall + 1
@ -457,7 +457,7 @@ subroutine solve_linter (irr, imode0, npe, drhoscf)
! of the change of potential and Q
!
call newdq (dvscfin, irr, npe)
call newdq (dvscfin, npe)
#ifdef PARA
aux_avg (1) = dfloat (ltaver)
aux_avg (2) = dfloat (lintercall)

14
PH/zstar_eu.f90
View File

@ -16,7 +16,6 @@ subroutine zstar_eu
use pwcom
use allocate
use parameters, only : DP
use phcom
implicit none
@ -41,6 +40,14 @@ subroutine zstar_eu
complex(kind=DP) :: ZDOTC
! scalar product
call setv (2 * 9 * nat, 0.d0, zstareu0, 1)
call setv (9 * nat, 0.d0, zstareu, 1)
if (okvan) then
call error('zstar_eu', 'Effective charges not implemented',-1)
return
endif
if (nksq.gt.1) rewind (iunigk)
do ik = 1, nksq
if (nksq.gt.1) read (iunigk) npw, igk
@ -64,8 +71,8 @@ subroutine zstar_eu
!
call davcio (dpsi, lrdwf, iudwf, nrec, - 1)
do ibnd = 1, nbnd
zstareu0 (jpol, mode) = zstareu0 (jpol, mode) + 2.d0 * weight * &
ZDOTC (npw, dpsi (1, ibnd), 1, dvpsi (1, ibnd), 1)
zstareu0(jpol,mode)=zstareu0(jpol, mode)-2.d0*weight*&
ZDOTC(npw,dpsi(1,ibnd),1,dvpsi(1,ibnd),1)
enddo
enddo
enddo
@ -77,7 +84,6 @@ subroutine zstar_eu
call poolreduce (18 * nat, zstareu0)
#endif
call setv (9 * nat, 0.d0, zstareu, 1)
!
! bring the mode index to cartesian coordinates
!

7
configure vendored
View File

@ -22,9 +22,14 @@ then
echo OSHOME=$PWD > make.sys
if [ "$1" = "pc_ifc" ]
then
echo creating the library file
echo creating catalog file for Intel compiler
echo work.pc > $PWD/Modules/intel.pcl
echo work.pc > $PWD/upftools/intel.pcl
echo work.pc > $PWD/flib/intel.pcl
echo $PWD/Modules/work.pc >> $PWD/flib/intel.pcl
echo work.pc > $PWD/PW/intel.pcl
echo $PWD/Modules/work.pc >> $PWD/PW/intel.pcl

7
install/configure vendored
View File

@ -22,9 +22,14 @@ then
echo OSHOME=$PWD > make.sys
if [ "$1" = "pc_ifc" ]
then
echo creating the library file
echo creating catalog file for Intel compiler
echo work.pc > $PWD/Modules/intel.pcl
echo work.pc > $PWD/upftools/intel.pcl
echo work.pc > $PWD/flib/intel.pcl
echo $PWD/Modules/work.pc >> $PWD/flib/intel.pcl
echo work.pc > $PWD/PW/intel.pcl
echo $PWD/Modules/work.pc >> $PWD/PW/intel.pcl