Commit c2bdfb00 authored by dalcorso's avatar dalcorso

Further cleanup of the phonon code. The indeces of k and k+q are no more

recalculated in every routine.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@5371 c92efa57-630b-4861-b058-cf58834340f0
parent 27c6cdf1
......@@ -25,8 +25,7 @@ subroutine adddvscf (ipert, ik)
USE wvfct, ONLY : nbnd, npwx
USE noncollin_module, ONLY : noncolin, npol
! modules from phcom
USE control_ph, ONLY : lgamma
USE qpoint, ONLY : npwq
USE qpoint, ONLY : npwq, ikks
USE phus, ONLY : int3, int3_nc, becp1, becp1_nc
USE eqv, ONLY : dvpsi
implicit none
......@@ -54,11 +53,7 @@ subroutine adddvscf (ipert, ik)
if (.not.okvan) return
call start_clock ('adddvscf')
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
......
......@@ -23,8 +23,8 @@ subroutine addusdbec (ik, wgt, psi, dbecsum)
USE uspp, only: nkb, vkb, okvan
USE uspp_param, only: upf, nh, nhm
USE phus, ONLY : becp1
USE qpoint, ONLY : npwq
USE control_ph, ONLY : nbnd_occ, lgamma
USE qpoint, ONLY : npwq, ikks
USE control_ph, ONLY : nbnd_occ
implicit none
!
! the dummy variables
......@@ -61,11 +61,7 @@ subroutine addusdbec (ik, wgt, psi, dbecsum)
call start_clock ('addusdbec')
allocate (dbecq( nkb, nbnd))
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
!
! First compute the product of psi and vkb
!
......
......@@ -24,9 +24,9 @@ subroutine addusdbec_nc (ik, wgt, psi, dbecsum_nc)
USE uspp, only: nkb, vkb, okvan
USE noncollin_module, ONLY : noncolin, npol
USE uspp_param, only: upf, nh, nhm
USE qpoint, ONLY : npwq
USE qpoint, ONLY : npwq, ikks
USE phus, ONLY : becp1, becp1_nc
USE control_ph, ONLY : nbnd_occ, lgamma
USE control_ph, ONLY : nbnd_occ
implicit none
!
......@@ -63,11 +63,7 @@ subroutine addusdbec_nc (ik, wgt, psi, dbecsum_nc)
call start_clock ('addusdbec_nc')
allocate (dbecq_nc( nkb,npol, nbnd))
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
!
! First compute the product of psi and vkb
!
......
......@@ -20,8 +20,9 @@ subroutine cch_psi_all (n, h, ah, e, ik, m)
USE uspp, ONLY: nkb, vkb
USE wvfct, ONLY : npwx, nbnd
USE control_ph, ONLY : alpha_pv, nbnd_occ, lgamma
USE control_ph, ONLY : alpha_pv, nbnd_occ
USE eqv, ONLY : evq
USE qpoint, ONLY : ikqs
USE mp_global, ONLY: intra_pool_comm
USE mp, ONLY: mp_sum
......@@ -75,11 +76,7 @@ subroutine cch_psi_all (n, h, ah, e, ik, m)
!
! Here we compute the projector in the valence band
!
if (lgamma) then
ikq = ik
else
ikq = 2 * ik
endif
ikq = ikqs(ik)
ps (:,:) = (0.d0, 0.d0)
call ZGEMM ('C', 'N', nbnd_occ (ikq) , m, n, (1.d0, 0.d0) , evq, &
......
......@@ -21,8 +21,9 @@ subroutine ch_psi_all (n, h, ah, e, ik, m)
USE uspp, ONLY: nkb, vkb
USE noncollin_module, ONLY : noncolin, npol
USE control_ph, ONLY : alpha_pv, nbnd_occ, lgamma
USE control_ph, ONLY : alpha_pv, nbnd_occ
USE eqv, ONLY : evq
USE qpoint, ONLY : ikqs
USE mp_global, ONLY: intra_pool_comm
USE mp, ONLY: mp_sum
......@@ -84,11 +85,7 @@ subroutine ch_psi_all (n, h, ah, e, ik, m)
!
! Here we compute the projector in the valence band
!
if (lgamma) then
ikq = ik
else
ikq = 2 * ik
endif
ikq = ikqs(ik)
ps (:,:) = (0.d0, 0.d0)
IF (noncolin) THEN
......
......@@ -28,8 +28,8 @@ subroutine compute_alphasum
USE phus, ONLY : alphasum, alphasum_nc, becp1, becp1_nc, alphap, &
alphap_nc
USE qpoint, ONLY : nksq
USE control_ph, ONLY : nbnd_occ, lgamma
USE qpoint, ONLY : nksq, ikks, ikqs
USE control_ph, ONLY : nbnd_occ
implicit none
......@@ -51,13 +51,8 @@ subroutine compute_alphasum
alphasum = 0.d0
IF (noncolin) alphasum_nc=(0.d0,0.d0)
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
ikk = ikks(ik)
ikq = ikqs(ik)
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
......
......@@ -29,7 +29,7 @@ subroutine compute_becalp (becq, alpq)
USE control_ph, ONLY : lgamma
USE eqv, ONLY : evq
USE units_ph, ONLY : lrwfc, iuwfc
USE qpoint, ONLY : nksq, npwq, igkq
USE qpoint, ONLY : nksq, npwq, igkq, ikqs
implicit none
......@@ -51,7 +51,7 @@ subroutine compute_becalp (becq, alpq)
if (nksq.gt.1) rewind (iunigk)
do ik = 1, nksq
ikq = 2 * ik
ikq = ikqs(ik)
if (nksq.gt.1) then
read (iunigk, err = 100, iostat = ios) npw, igk
......
......@@ -28,8 +28,8 @@ subroutine compute_becsum_ph
USE phus, ONLY : alphasum, alphasum_nc, becp1, becp1_nc, alphap, &
alphap_nc, becsum_nc
USE qpoint, ONLY : nksq
USE control_ph, ONLY : nbnd_occ, lgamma
USE qpoint, ONLY : nksq, ikks, ikqs
USE control_ph, ONLY : nbnd_occ
implicit none
......@@ -42,13 +42,8 @@ subroutine compute_becsum_ph
IF (noncolin) becsum_nc = (0.d0,0.d0)
becsum = 0.d0
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
ikk = ikks(ik)
ikq = ikqs(ik)
if (lsda) current_spin = isk (ikk)
ijkb0 = 0
do nt = 1, ntyp
......
......@@ -27,7 +27,7 @@ subroutine compute_drhous (drhous, dbecsum, wgg, becq, alpq)
USE gsmooth, ONLY : nrxxs, nr1s,nr2s,nr3s, nrx1s,nrx2s,nrx3s, nls
USE wvfct, ONLY : npw, nbnd, igk
USE qpoint, ONLY : nksq, igkq, npwq
USE qpoint, ONLY : nksq, igkq, npwq, ikks, ikqs
USE eqv, ONLY : evq
USE units_ph, ONLY : iuwfc, lrwfc
USE control_ph, ONLY : lgamma
......@@ -80,14 +80,9 @@ subroutine compute_drhous (drhous, dbecsum, wgg, becq, alpq)
read (iunigk, err = 110, iostat = ios) npw, igk
110 call errore ('compute_drhous', 'reading igk', abs (ios) )
endif
if (lgamma) then
ikk = ik
ikq = ik
npwq = npw
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lgamma) npwq = npw
ikk = ikks(ik)
ikq = ikqs(ik)
weight = wk (ikk)
if (lsda) current_spin = isk (ikk)
......
......@@ -29,7 +29,7 @@ subroutine compute_drhous_nc (drhous, dbecsum, wgg, becq, alpq)
USE uspp, ONLY: okvan, nkb, vkb
USE uspp_param, ONLY: nhm
USE qpoint, ONLY : nksq, igkq, npwq
USE qpoint, ONLY : nksq, igkq, npwq, ikks, ikqs
USE eqv, ONLY : evq
USE units_ph, ONLY : lrwfc, iuwfc
USE control_ph, ONLY : lgamma
......@@ -83,14 +83,9 @@ subroutine compute_drhous_nc (drhous, dbecsum, wgg, becq, alpq)
read (iunigk, err = 110, iostat = ios) npw, igk
110 call errore ('compute_drhous', 'reading igk', abs (ios) )
endif
if (lgamma) then
ikk = ik
ikq = ik
npwq = npw
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lgamma) npwq = npw
ikk = ikks(ik)
ikq = ikqs(ik)
weight = wk (ikk)
if (lsda) current_spin = isk (ikk)
......
......@@ -27,11 +27,11 @@ subroutine compute_nldyn (wdyn, wgg, becq, alpq)
USE spin_orb, ONLY : lspinorb
USE wvfct, ONLY : nbnd, et
USE qpoint, ONLY : nksq
USE qpoint, ONLY : nksq, ikks, ikqs
USE modes, ONLY : u
USE phus, ONLY : becp1, becp1_nc, alphap, alphap_nc, int1, int2, &
int2_so, int1_nc
USE control_ph, ONLY : nbnd_occ, lgamma
USE control_ph, ONLY : nbnd_occ
USE mp_global, ONLY: intra_pool_comm
USE mp, ONLY: mp_sum
......@@ -75,13 +75,8 @@ subroutine compute_nldyn (wdyn, wgg, becq, alpq)
dynwrk (:,:) = (0.d0, 0.d0)
call divide (nbnd, startb, lastb)
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
ikk = ikks(ik)
ikq = ikqs(ik)
if (lsda) current_spin = isk (ikk)
IF (noncolin) THEN
......
......@@ -18,8 +18,7 @@ subroutine compute_weight (wgg)
USE klist, ONLY : wk, lgauss, degauss, ngauss
USE ener, ONLY : ef
USE wvfct, ONLY : nbnd, wg, et
USE qpoint, ONLY : nksq
USE control_ph, ONLY : lgamma
USE qpoint, ONLY : nksq, ikks, ikqs
implicit none
real(DP) :: wgg (nbnd, nbnd, nksq)
......@@ -35,13 +34,8 @@ subroutine compute_weight (wgg)
! the weights are computed for each k point ...
!
do ik = 1, nksq
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
ikk = ikks(ik)
ikq = ikqs(ik)
!
! each band v ...
!
......
......@@ -17,7 +17,7 @@ subroutine deallocate_phq
USE ramanm, ONLY: ramtns
USE modes, ONLY : tmq, t, npert, u, ubar, rtau
USE qpoint, ONLY : eigqts, igkq
USE qpoint, ONLY : eigqts, igkq, ikks, ikqs
USE efield_mod, ONLY : zstareu, zstarue, zstarue0, zstareu0, &
zstarue0_rec
USE phus, ONLY : int1, int1_nc, int2, int2_so, int3, int3_nc, int3_paw, &
......@@ -52,6 +52,8 @@ subroutine deallocate_phq
if(allocated(vlocq)) deallocate (vlocq)
if(allocated(dmuxc)) deallocate (dmuxc)
!
if(allocated(ikks)) deallocate (ikks)
if(allocated(ikqs)) deallocate (ikqs)
if(allocated(eigqts)) deallocate (eigqts)
if(allocated(rtau)) deallocate (rtau)
if(associated(u)) deallocate (u)
......
......@@ -29,7 +29,7 @@ subroutine drhodv (nu_i0, nper, drhoscf)
USE dynmat, ONLY : dyn, dyn_rec
USE modes, ONLY : u, npertx
USE qpoint, ONLY : nksq, npwq, igkq
USE qpoint, ONLY : nksq, npwq, igkq, ikks, ikqs
USE eqv, ONLY : dpsi
USE units_ph, ONLY : lrdwf, iuwfc, iudwf
USE control_ph, ONLY : lgamma
......@@ -76,13 +76,11 @@ subroutine drhodv (nu_i0, nper, drhoscf)
if (nksq > 1) rewind (unit = iunigk)
do ik = 1, nksq
if (nksq > 1) read (iunigk) npw, igk
ikk = ikks(ik)
ikq = ikqs(ik)
if (lgamma) then
ikk = ik
ikq = ik
npwq = npw
else
ikk = 2 * ik - 1
ikq = ikk + 1
if (nksq > 1) read (iunigk) npwq, igkq
endif
if (lsda) current_spin = isk (ikk)
......
......@@ -31,8 +31,7 @@ subroutine dvqpsi_us (ik, mode, uact, addnlcc)
USE wavefunctions_module, ONLY: evc
USE nlcc_ph, ONLY : nlcc_any, drc
USE eqv, ONLY : dvpsi, dmuxc, vlocq
USE qpoint, ONLY : npwq, igkq, xq, eigqts
USE control_ph, ONLY : lgamma
USE qpoint, ONLY : npwq, igkq, xq, eigqts, ikks
implicit none
!
......@@ -81,11 +80,7 @@ subroutine dvqpsi_us (ik, mode, uact, addnlcc)
! reciprocal space while the product with the wavefunction is done in
! real space
!
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
dvpsi(:,:) = (0.d0, 0.d0)
aux1(:) = (0.d0, 0.d0)
do na = 1, nat
......
......@@ -28,7 +28,7 @@ subroutine dvqpsi_us_only (ik, mode, uact)
USE noncollin_module, ONLY : noncolin, npol
USE uspp, ONLY: okvan, nkb, vkb, qq, qq_so, deeq, deeq_nc
USE uspp_param, ONLY: nh
USE qpoint, ONLY : igkq, npwq
USE qpoint, ONLY : igkq, npwq, ikks, ikqs
USE phus, ONLY : int1, int1_nc, int2, int2_so, alphap, alphap_nc, &
becp1, becp1_nc
USE eqv, ONLY : dvpsi
......@@ -82,13 +82,8 @@ subroutine dvqpsi_us_only (ik, mode, uact)
allocate (ps2 ( nkb , nbnd , 3))
end if
allocate (aux ( npwx))
if (lgamma) then
ikk = ik
ikq = ik
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
ikk = ikks(ik)
ikq = ikqs(ik)
if (lsda) current_spin = isk (ikk)
!
! we first compute the coefficients of the vectors
......
......@@ -34,7 +34,7 @@ SUBROUTINE dynmat_us()
USE becmod, ONLY : calbec
USE io_global, ONLY : stdout
USE qpoint, ONLY : npwq, nksq, igkq
USE qpoint, ONLY : npwq, nksq, igkq, ikks
USE modes, ONLY : u
USE dynmat, ONLY : dyn
USE phus, ONLY : becp1, becp1_nc, alphap, alphap_nc
......@@ -125,11 +125,7 @@ SUBROUTINE dynmat_us()
!
IF (nksq > 1) REWIND (unit = iunigk)
DO ik = 1, nksq
IF (lgamma) THEN
ikk = ik
ELSE
ikk = 2 * ik - 1
ENDIF
ikk = ikks(ik)
IF (lsda) current_spin = isk (ikk)
IF (nksq > 1) READ (iunigk) npw, igk
......
......@@ -23,8 +23,8 @@ subroutine incdrhoscf (drhoscf, weight, ik, dbecsum)
USE uspp_param,ONLY: nhm
USE wavefunctions_module, ONLY: evc
USE eqv, ONLY : dpsi
USE qpoint, ONLY : npwq, igkq
USE control_ph, ONLY : nbnd_occ, lgamma
USE qpoint, ONLY : npwq, igkq, ikks
USE control_ph, ONLY : nbnd_occ
implicit none
......@@ -54,11 +54,7 @@ subroutine incdrhoscf (drhoscf, weight, ik, dbecsum)
allocate (dpsic( nrxxs))
allocate (psi ( nrxxs))
wgt = 2.d0 * weight / omega
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
!
! dpsi contains the perturbed wavefunctions of this k point
! evc contains the unperturbed wavefunctions of this k point
......
......@@ -26,9 +26,9 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum)
USE uspp_param,ONLY : nhm
USE wvfct, ONLY : npw, npwx, igk
USE wavefunctions_module, ONLY: evc
USE qpoint, ONLY : npwq, igkq
USE qpoint, ONLY : npwq, igkq, ikks
USE eqv, ONLY : dpsi
USE control_ph, ONLY : nbnd_occ, lgamma
USE control_ph, ONLY : nbnd_occ
implicit none
......@@ -58,11 +58,7 @@ subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum)
allocate (dpsic( nrxxs, npol))
allocate (psi ( nrxxs, npol))
wgt = 2.d0 * weight / omega
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
!
! dpsi contains the perturbed wavefunctions of this k point
! evc contains the unperturbed wavefunctions of this k point
......
......@@ -24,9 +24,9 @@ subroutine incdrhous (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
USE uspp, ONLY : nkb, qq
USE uspp_param,ONLY : nhm, nh
USE wvfct, ONLY : nbnd, npwx
USE qpoint, ONLY : nksq, igkq, npwq
USE qpoint, ONLY : nksq, igkq, npwq, ikks
USE phus, ONLY : becp1, alphap
USE control_ph, ONLY: nbnd_occ, lgamma
USE control_ph, ONLY: nbnd_occ
USE eqv, ONLY : evq, dpsi
USE modes, ONLY : u
USE mp_global, ONLY: intra_pool_comm
......@@ -70,11 +70,7 @@ subroutine incdrhous (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
call divide (nbnd, startb, lastb)
ps1 (:,:) = (0.d0, 0.d0)
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk=ikks(ik)
!
! Here we prepare the two terms
!
......
......@@ -28,9 +28,9 @@ subroutine incdrhous_nc (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
USE uspp_param,ONLY : nhm, nh
USE wvfct, ONLY : nbnd, npwx
USE modes, ONLY : u
USE qpoint, ONLY : npwq, nksq, igkq
USE qpoint, ONLY : npwq, nksq, igkq, ikks
USE eqv, ONLY : dpsi, evq
USE control_ph, ONLY : nbnd_occ, lgamma
USE control_ph, ONLY : nbnd_occ
USE phus, ONLY : becp1_nc, alphap_nc
USE mp_global, ONLY: intra_pool_comm
USE mp, ONLY: mp_sum
......@@ -74,11 +74,7 @@ subroutine incdrhous_nc (drhoscf, weight, ik, dbecsum, evcr, wgg, becq, &
call divide (nbnd, startb, lastb)
ps1 (:,:) = (0.d0, 0.d0)
if (lgamma) then
ikk = ik
else
ikk = 2 * ik - 1
endif
ikk = ikks(ik)
!
! Here we prepare the two terms
!
......
......@@ -74,6 +74,7 @@ END MODULE dynmat
!
MODULE qpoint
USE kinds, ONLY : DP
USE parameters, ONLY : npk
!
! ... The q point
!
......@@ -84,6 +85,9 @@ MODULE qpoint
INTEGER :: nksq, npwq
! the real number of k points
! the number of plane waves for q
INTEGER, ALLOCATABLE :: ikks(:), ikqs(:)
! the index of k point in the list of k
! the index of k+q point in the list of k
REAL (DP) :: xq(3)
! the coordinates of the q point
COMPLEX (DP), ALLOCATABLE :: eigqts(:) ! nat)
......@@ -192,7 +196,9 @@ MODULE phus
REAL (DP), ALLOCATABLE :: &
alphasum(:,:,:,:), &! nhm*(nhm+1)/2,3,nat,nspin)
! used to compute modes
dpqq(:,:,:,:) ! dipole moment of each Q
dpqq(:,:,:,:) ! (nhm, nhm, 3, ntyp)
! alphasum contains \sum_i <psi_i| d/du (|\beta_n><beta_m|) | psi_i> + (m-n)
! dipole moment of each Q
COMPLEX (DP), ALLOCATABLE :: &
int1(:,:,:,:,:), &! nhm, nhm, 3, nat, nspin),&
int2(:,:,:,:,:), &! nhm, nhm, 3,nat, nat),&
......@@ -201,27 +207,38 @@ MODULE phus
int4(:,:,:,:,:), &! nhm*(nhm+1)/2, 3, 3, nat, nspin),&
int5(:,:,:,:,:), &! nhm*(nhm+1)/2, 3, 3, nat, nat),&
int1_nc(:,:,:,:,:), &! nhm, nhm, 3, nat, nspin),&
int2_so(:,:,:,:,:,:), &! nhm, nhm, 3,nat,nat,nspin),&
int2_so(:,:,:,:,:,:), &! nhm, nhm, 3, nat,nat,nspin),&
int3_nc(:,:,:,:,:), &! nhm, nhm, max_irr_dim, nat, nspin),&
int4_nc(:,:,:,:,:,:), &! nhm, nhm, 3, 3, nat, nspin),&
int5_so(:,:,:,:,:,:,:), &! nhm*(nhm+1)/2, 3, 3, nat, nat, nspin),&
!
! These variables contains the five integrals defined in PRB 64, 35118 (2001)
! int1 -> \int V_eff d/du (Q) d^3r
! int2 -> \int d/du (V_loc) Q d^3r
! int3 -> \int d\du (V_Hxc) Q d^3r
! int4 -> \int V_eff d^2/dudu (Q) d^3r
! int5 -> \int d/du (V_loc) d/du (Q) d^3r
!
! int3_paw contains d/du (D^1-\tilde D^1)
!
!
becsum_nc(:,:,:,:), &! nhm*(nhm+1)/2,nat,npol,npol)
becsumort(:,:,:,:), &! nhm*(nhm+1)/2,nat,nspin,3*nat)
alphasum_nc(:,:,:,:,:), &! nhm*(nhm+1)/2,3,nat,npol,npol)
dpqq_so(:,:,:,:,:) ! dipole moment of each Q and the fcoef factors
dpqq_so(:,:,:,:,:) ! nhm, nhm, nspin, 3, ntyp
!
! becsum contains \sum_i <\psi_i | \beta_n><\beta_m| \psi_i > + (m-n)
! besumort contains alphasum+\sum_i <\psi_i | \beta_n><\beta_m| \delta \psi_i >
! dpqq_so dipole moment of each Q multiplied by the fcoef factors
!
COMPLEX (DP), ALLOCATABLE, TARGET :: &
becp1(:,:,:), &! nkbtot, nbnd, nksq),&
becp1_nc(:,:,:,:), &! nkbtot, npol, nbnd, nksq),&
alphap(:,:,:,:), &! nkbtot, nbnd, 3, nksq)
alphap_nc(:,:,:,:,:) ! nkbtot, npol, nbnd, 3, nksq)
! integrals of dQ and V_eff
! integrals of dQ and V_loc
! integrals of Q and dV_Hxc
! integrals of d^2Q and V
! integrals of dQ and dV_lo
! the becq used in ch_psi
! the derivative of the bec
!
! becp1 contains < beta_n | \psi_i >
! alphap contains < d\du (\beta_n) | psi_i>
!
END MODULE phus
!
......
......@@ -28,7 +28,7 @@ PROGRAM phonon
USE start_k, ONLY : xk_start, wk_start, nks_start
USE noncollin_module,ONLY : noncolin
USE control_flags, ONLY : restart
USE qpoint, ONLY : xq, nksq
USE qpoint, ONLY : xq, nksq, ikks, ikqs
USE modes, ONLY : nirr
USE partial, ONLY : done_irr, comp_irr
USE disp, ONLY : nqs, x_q, done_iq, rep_iq, done_rep_iq
......@@ -49,7 +49,7 @@ PROGRAM phonon
!
IMPLICIT NONE
!
INTEGER :: iq, iq_start, ierr, iu
INTEGER :: iq, iq_start, ierr, iu, ik
INTEGER :: irr
LOGICAL :: exst, do_band
CHARACTER (LEN=9) :: code = 'PHONON'
......@@ -273,10 +273,20 @@ PROGRAM phonon
IF ( lgamma ) THEN
!
nksq = nks
ALLOCATE(ikks(nksq), ikqs(nksq))
DO ik=1,nksq
ikks(ik) = ik
ikqs(ik) = ik
ENDDO
!
ELSE
!
nksq = nks / 2
ALLOCATE(ikks(nksq), ikqs(nksq))
DO ik=1,nksq
ikks(ik) = 2 * ik - 1
ikqs(ik) = 2 * ik
ENDDO
!
END IF
!
......
......@@ -57,7 +57,7 @@ SUBROUTINE phq_init()
USE nlcc_ph, ONLY : nlcc_any
USE control_ph, ONLY : zue, epsil, lgamma, all_done
USE units_ph, ONLY : lrwfc, iuwfc
USE qpoint, ONLY : xq, igkq, npwq, nksq, eigqts
USE qpoint, ONLY : xq, igkq, npwq, nksq, eigqts, ikks, ikqs
USE paw_onecenter, ONLY : PAW_potential, PAW_symmetrize
USE paw_variables, ONLY : okpaw, ddd_paw
!
......@@ -125,13 +125,8 @@ SUBROUTINE phq_init()
!
DO ik = 1, nksq
!
IF ( lgamma ) THEN
ikk = ik
ikq = ik
ELSE
ikk = 2 * ik - 1
ikq = ikk + 1
END IF
ikk = ikks(ik)
ikq = ikqs(ik)
!
IF ( lsda ) current_spin = isk( ikk )
!
......
......@@ -56,11 +56,10 @@ subroutine solve_linter (irr, imode0, npe, drhoscf)
USE output, ONLY : fildrho, fildvscf
USE phus, ONLY : int1, int2, int3, int3_paw, becsumort
USE eqv, ONLY : dvpsi, dpsi, evq
USE qpoint, ONLY : npwq, igkq, nksq
USE qpoint, ONLY : xq, npwq, igkq, nksq, ikks, ikqs
USE modes, ONLY : npert, u, t, max_irr_dim, irotmq, tmq, &
minus_q, irgq, nsymq, rtau
USE efield_mod, ONLY : zstareu0, zstarue0
USE qpoint, ONLY : xq
! used oly to write the restart file
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
USE mp, ONLY : mp_sum
......@@ -239,14 +238,9 @@ subroutine solve_linter (irr, imode0, npe, drhoscf)
read (iunigk, err = 100, iostat = ios) npw, igk
100 call errore ('solve_linter', 'reading igk', abs (ios) )
endif
if (lgamma) then
ikk = ik
ikq = ik
npwq = npw
else
ikk = 2 * ik - 1
ikq = ikk + 1
endif
if (lgamma) npwq = npw
ikk = ikks(ik)
ikq = ikqs(ik)
if (lsda) current_spin = isk (ikk)
if (.not.lgamma.and.nksq.gt.1) then
read (iunigk, err = 200, iostat = ios) npwq, igkq
......