Commit 509459fc authored by giannozz's avatar giannozz

Cleanup in PH/ and D3/, setv removed (please verify collateral damages!)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@698 c92efa57-630b-4861-b058-cf58834340f0
parent b8b20875
......@@ -146,7 +146,6 @@ PHOBJS = \
../PH/setlocq.o \
../PH/setqmod.o \
../PH/setup_dgc.o \
../PH/setv.o \
../PH/smallgq.o \
../PH/solve_e.o \
../PH/solve_linter.o \
......
......@@ -46,15 +46,15 @@ subroutine allocate_d3
allocate (d3dyn_aux7 ( 3 * nat, 3 * nat, 3 * nat))
allocate (d3dyn_aux8 ( 3 * nat, 3 * nat, 3 * nat))
allocate (d3dyn_aux9 ( 3 * nat, 3 * nat, 3 * nat))
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux1, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux2, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux3, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux4, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux5, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux6, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux7, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux8, 1)
call setv (54 * nat * nat * nat, 0.d0, d3dyn_aux9, 1)
d3dyn_aux1 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux2 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux3 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux4 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux5 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux6 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux7 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux8 (:,:,:) = (0.d0, 0.d0)
d3dyn_aux9 (:,:,:) = (0.d0, 0.d0)
return
end subroutine allocate_d3
......@@ -5,9 +5,7 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine ch_psi_all2 (n, h, ah, e, ik, m)
!-----------------------------------------------------------------------
!
......@@ -52,8 +50,8 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m)
allocate (hpsi( npwx, m))
allocate (spsi( npwx, m))
call setv (2 * npwx * m, 0.d0, hpsi, 1)
call setv (2 * npwx * m, 0.d0, spsi, 1)
hpsi = (0.d0, 0.d0)
spsi = (0.d0, 0.d0)
!
! compute the product of the hamiltonian with the h vector
!
......@@ -71,24 +69,24 @@ subroutine ch_psi_all2 (n, h, ah, e, ik, m)
!
! Here we compute the projector in the valence band
!
call setv (2 * npwx * m, 0.d0, hpsi, 1)
hpsi = (0.d0, 0.d0)
if (lgamma) then
ikq = ik
else
ikq = 2 * ik
endif
call setv (2 * nbnd * m, 0.d0, ps, 1)
ps = (0.d0, 0.d0)
call ZGEMM ('C', 'N', nbnd, m, n, (1.d0, 0.d0) , evq, npwx, spsi, &
npwx, (0.d0, 0.d0) , ps, nbnd)
call DSCAL (2 * nbnd * m, alpha_pv, ps, 1)
ps = ps * alpha_pv
#ifdef __PARA
call reduce (2 * nbnd * m, ps)
#endif
call ZGEMM ('N', 'N', n, m, nbnd, (1.d0, 0.d0) , evq, npwx, ps, &
nbnd, (1.d0, 0.d0) , hpsi, npwx)
call ZCOPY (npwx * m, hpsi, 1, spsi, 1)
spsi = hpsi
!
! And apply S again
!
......
......@@ -62,16 +62,16 @@ subroutine d0rhod2v (ipert, drhoscf)
allocate (work5(npwx))
allocate (work6(npwx))
call setv (2*9*nat*nat,0.0d0,d3dywrk,1)
d3dywrk (:,:) = (0.d0, 0.d0)
!
! Here the contribution deriving from the local part of the potential
#ifdef __PARA
! ... computed only by the first pool (no sum over k needed)
!
if (mypool.ne.1) goto 100
if (mypool /= 1) goto 100
#endif
!
call ZCOPY (nrxx, drhoscf, 1, work0, 1)
work0 (:) = drhoscf (:)
call cft3 (work0, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
do na = 1, nat
do icart = 1,3
......@@ -97,7 +97,6 @@ subroutine d0rhod2v (ipert, drhoscf)
WRITE( stdout,'(3(2f10.6,2x))') &
((d3dywrk(3*(na-1)+icart,3*(na-1)+jcart), &
jcart=1,3),icart=1,3)
enddo
#ifdef __PARA
call reduce(2*9*nat*nat,d3dywrk)
......@@ -136,7 +135,7 @@ subroutine d0rhod2v (ipert, drhoscf)
! In the metallic case corrects dpsi so as that the density matrix
! will be: Sum_{k,nu} 2 * | dpsi > < psi |
!
if (degauss.ne.0.d0) then
if (degauss /= 0.d0) then
nrec = ipert + (ik - 1) * 3 * nat
call davcio (psidqvpsi, lrpdqvp, iupd0vp, nrec, - 1)
call dpsi_corr (evc, psidqvpsi, ikk, ikk, ipert)
......@@ -155,7 +154,7 @@ subroutine d0rhod2v (ipert, drhoscf)
jkb=0
do nt = 1, ntyp
do na = 1, nat
if (ityp (na).eq.nt) then
if (ityp (na) == nt) then
na_icart = 3 * (na - 1) + icart
na_jcart = 3 * (na - 1) + jcart
do ikb = 1, nh (nt)
......@@ -211,8 +210,8 @@ subroutine d0rhod2v (ipert, drhoscf)
d3dyn(nu_i,nu_j,nu_k) = d3dyn(nu_i,nu_j,nu_k) + work
endif
enddo
enddo
deallocate (work6)
deallocate (work5)
deallocate (work4)
......
......@@ -39,16 +39,16 @@ subroutine d3_exc
!
! Calculates third derivative of Exc
!
call setv (nrxx, 0.d0, d2muxc, 1)
d2muxc(:) = 0.d0
do ir = 1, nrxx
rhotot = rho (ir, 1) + rho_core (ir)
if (rhotot.gt.1.d-30) d2muxc (ir) = d2mxc (rhotot)
if (rhotot.lt. - 1.d-30) d2muxc (ir) = - d2mxc ( - rhotot)
if (rhotot > 1.d-30) d2muxc (ir) = d2mxc (rhotot)
if (rhotot < - 1.d-30) d2muxc (ir) = - d2mxc ( - rhotot)
enddo
!
! Calculates the contribution to d3dyn
!
call setv (2 * 27 * nat * nat * nat, 0.d0, d3dyn1, 1)
d3dyn1 (:,:,:) = (0.d0, 0.d0)
do ipert = 1, 3 * nat
if (q0mode (ipert) ) then
call davcio_drho (work1, lrdrho, iud0rho, ipert, - 1)
......@@ -75,8 +75,8 @@ subroutine d3_exc
call poolbcast (2 * 27 * nat * nat * nat, d3dyn1)
#endif
call DAXPY (2 * 27 * nat * nat * nat, 1.d0, d3dyn1, 1, d3dyn, 1)
call ZCOPY (27 * nat * nat * nat, d3dyn1, 1, d3dyn_aux9, 1)
d3dyn = d3dyn + d3dyn1
d3dyn_aux9 = d3dyn1
deallocate (d2muxc)
deallocate (work1)
......
......@@ -43,9 +43,9 @@ subroutine d3_init
if (.not.lgamma) then
allocate (d0rc( ngm, ntyp))
call setv (3, 0.d0, work, 1)
work = 0.d0
call set_drhoc (work)
call ZCOPY (ngm * ntyp, drc, 1, d0rc, 1)
d0rc (:,:) = drc (:,:)
else
d0rc => drc
endif
......@@ -62,8 +62,8 @@ subroutine d3_init
! the fourier components of the local potential at q+G for q=0
!
if (.not.lgamma) then
call setv (ngm * ntyp, 0.d0, vlocg0, 1)
call setv (3, 0.d0, work, 1)
vlocg0 (:,:) = 0.d0
work = 0.d0
do nt = 1, ntyp
call setlocq (work, lloc(nt), lmax(nt), numeric(nt), &
mesh(nt), msh(nt), rab(1,nt), r(1,nt), vnl(1,lloc(nt),nt), &
......@@ -75,7 +75,7 @@ subroutine d3_init
! Reads the q=0 variation of the charge --d0rho-- and symmetrizes it
!
#ifdef __PARA
! if (mypool.ne.1) goto 100
! if (mypool /= 1) goto 100
#endif
do irr = 1, nirrg0
imode0 = 0
......
......@@ -65,23 +65,11 @@ subroutine d3_setup
! computes derivative of xc potential
! working array
integer :: ir, table (48, 48), isym, jsym, iinv, irot, jrot, ik, &
integer :: table (48, 48)
! the multiplication table of the point group
integer :: ir, isym, jsym, iinv, irot, jrot, ik, &
ibnd, ipol, mu, nu, imode0, irr, ipert, nt, ii, nu_i
! counter on mesh points
! the multiplication table of the point g
! counter on symmetries
! counter on symmetries
! the index of the inverse
! counter on rotations
! counter on rotations
! counter on k points
! counter on bands
! counter on polarizations
! counter on modes
! the starting mode
! counter on representation and perturbat
! counter on atomic type
! counters
logical :: sym (48)
! the symmetry operations
......@@ -96,25 +84,25 @@ subroutine d3_setup
!
! 2) Computes the derivative of the xc potential
!
call setv (nrxx * nspin * nspin, 0.d0, dmuxc, 1)
dmuxc (:,:,:) = 0.d0
if (lsda) then
do ir = 1, nrxx
rhoup = rho (ir, 1) + 0.5d0 * rho_core (ir)
rhodw = rho (ir, 2) + 0.5d0 * rho_core (ir)
call dmxc_spin (rhoup, rhodw, dmuxc (ir, 1, 1), dmuxc (ir, 2, &
1), dmuxc (ir, 1, 2), dmuxc (ir, 2, 2) )
call dmxc_spin (rhoup, rhodw, dmuxc (ir, 1, 1), &
dmuxc (ir, 2, 1), dmuxc (ir, 1, 2), dmuxc (ir, 2, 2) )
enddo
else
do ir = 1, nrxx
rhotot = rho (ir, nspin) + rho_core (ir)
if (rhotot.gt.1.d-30) dmuxc (ir, 1, 1) = dmxc (rhotot)
if (rhotot.lt. - 1.d-30) dmuxc (ir, 1, 1) = - dmxc ( - rhotot)
if (rhotot > 1.d-30) dmuxc (ir, 1, 1) = dmxc (rhotot)
if (rhotot < - 1.d-30) dmuxc (ir, 1, 1) = - dmxc ( - rhotot)
enddo
endif
!
! 3) Computes the number of occupated bands for each k point
!
if (degauss.ne.0.d0) then
if (degauss /= 0.d0) then
!
! discard conduction bands such that w0gauss(x,n) < small
!
......@@ -131,23 +119,21 @@ subroutine d3_setup
!
! - limit appropriated for Fermi-Dirac
!
if (ngauss.eq. - 99) then
if (ngauss == - 99) then
fac = 1.d0 / sqrt (small)
xmax = 2.d0 * log (0.5 * (fac + sqrt (fac * fac - 4.0) ) )
endif
target = ef + xmax * degauss
do ik = 1, nks
do ibnd = 1, nbnd
if (et (ibnd, ik) .lt.target) nbnd_occ (ik) = ibnd
if (et (ibnd, ik) < target) nbnd_occ (ik) = ibnd
enddo
if (nbnd_occ (ik) .eq.nbnd) &
if (nbnd_occ (ik) == nbnd) &
WRITE( stdout, '(5x,/,"Possibly too few bands at point ", &
& i4,3f10.5)') ik, (xk (ipol, ik) , ipol = 1, 3)
enddo
else
if (lsda) call errore ('d3_setup', 'occupation numbers probably wro &
&ng', - 1)
if (lsda) call errore ('d3_setup', 'occupation numbers probably wrong',-1)
do ik = 1, nks
nbnd_occ (ik) = nint (nelec) / degspin
enddo
......@@ -199,15 +185,15 @@ subroutine d3_setup
call sgam_ph (at, bg, nsym, s, irt, tau, rtau, nat, sym)
nmodes = 3 * nat
! if minus_q=.t. set_irr will search for
minus_q = (iswitch.gt. - 3)
minus_q = (iswitch > - 3)
! Sq=-q+G symmetry. On output minus_q=.t.
! if such a symmetry has been found
if (iswitch.eq. - 4) then
if (iswitch == - 4) then
call set_irr_mode (nat, at, bg, xq, s, invs, nsym, rtau, irt, &
irgq, nsymq, minus_q, irotmq, t, tmq, max_irr_dim, u, &
npert, nirr, gi, gimq, iverbosity, modenum)
else
if (nsym.gt.1) then
if (nsym > 1) then
call io_pattern(fildrho,nirr,npert,u,-1)
call set_sym_irr (nat, at, bg, xq, s, invs, nsym, rtau, irt, &
irgq, nsymq, minus_q, irotmq, t, tmq, max_irr_dim, u, &
......@@ -241,7 +227,7 @@ subroutine d3_setup
call multable (nsymg0, s, table)
do irot = 1, nsymg0
do jrot = 1, nsymg0
if (table (irot, jrot) .eq.1) invs (irot) = jrot
if (table (irot, jrot) == 1) invs (irot) = jrot
enddo
enddo
!
......@@ -263,7 +249,6 @@ subroutine d3_setup
nlcc_any = .false.
do nt = 1, ntyp
nlcc_any = nlcc_any.or.nlcc (nt)
enddo
if (nlcc_any) allocate (drc( ngm, ntyp))
......@@ -276,7 +261,7 @@ subroutine d3_setup
nlength_w = (3 * nat) / npool
nresto = 3 * nat - nlength_w * npool
do ii = 1, npool
if (ii.le.nresto) then
if (ii <= nresto) then
nlength (ii) = nlength_w + 1
else
nlength (ii) = nlength_w
......@@ -293,7 +278,7 @@ subroutine d3_setup
! 8) Sets up variables needed to calculate only selected
! modes at q=0 --the first index of the third order matrix--
!
if (q0mode_todo (1) .le.0) then
if (q0mode_todo (1) <= 0) then
do ii = 1, 3 * nat
q0mode (ii) = .true.
enddo
......@@ -302,7 +287,7 @@ subroutine d3_setup
q0mode (ii) = .false.
enddo
ii = 1
do while (q0mode_todo (ii) .gt.0)
do while (q0mode_todo (ii) > 0)
q0mode (q0mode_todo (ii) ) = .true.
ii = ii + 1
enddo
......@@ -312,7 +297,7 @@ subroutine d3_setup
! the calculation can be simplyfied, in this case allmodes
! is set .true.
!
allmodes = lgamma.and.q0mode_todo (1) .le.0
allmodes = lgamma.and.q0mode_todo (1) <= 0
!
! Sets up variables needed to write only selected
! modes at q=0 --the first index of the third order matrix--
......@@ -325,8 +310,7 @@ subroutine d3_setup
endif
enddo
wrmode (ii) = .false.
if (wrk.gt.1.d-8) wrmode (ii) = .true.
if (wrk > 1.d-8) wrmode (ii) = .true.
enddo
call stop_clock ('d3_setup')
return
......
......@@ -43,18 +43,7 @@ subroutine d3_symdyn (d3dyn, u, ug0, xq, s, invs, rtau, irt, irgq, &
! input: the patterns
integer :: i, j, i1, icart, jcart, kcart, na, nb, nc, mu, nu, om
! counter on modes
! counter on modes
! counter on modes
! counter on cartesian coordinates
! counter on cartesian coordinates
! counter on cartesian coordinates
! counter on atoms
! counter on atoms
! counter on atoms
! counter on modes
! counter on modes
! counter on modes
! counters
complex (kind = dp) :: work, wrk (3, 3)
! auxiliary variables
......@@ -65,7 +54,7 @@ subroutine d3_symdyn (d3dyn, u, ug0, xq, s, invs, rtau, irt, irgq, &
!
! First we transform in the cartesian coordinates
!
call setv (2 * 27 * nat * nat * nat, 0.d0, phi, 1)
phi = (0.d0, 0.d0)
do i1 = npert_i, npert_f
nc = (i1 - 1) / 3 + 1
kcart = i1 - 3 * (nc - 1)
......@@ -79,8 +68,8 @@ subroutine d3_symdyn (d3dyn, u, ug0, xq, s, invs, rtau, irt, irgq, &
do om = 1, 3 * nat
do mu = 1, 3 * nat
do nu = 1, 3 * nat
work = work + conjg (ug0 (i1, om) ) * u (i, mu) * d3dyn (om, mu, &
nu) * conjg (u (j, nu) )
work = work + conjg (ug0 (i1, om) ) * u (i, mu) * &
d3dyn (om, mu, nu) * conjg (u (j, nu) )
enddo
enddo
enddo
......
......@@ -18,6 +18,7 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
!
#include "machine.h"
USE kinds, only : DP
USE constants, only : tpi
implicit none
!
! The dummy variables
......@@ -38,32 +39,13 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
logical :: minus_q
! input: true if a symmetry q->-q+G
complex (kind = dp) :: phi (3, 3, 3, nat, nat, nat)
! inp/out: the matrix to symmetr
! inp/out: the matrix to symmetrize
!
! One parameter
! local variables
!
real (kind = dp) :: tpi
parameter (tpi = 2.0d0 * 3.14159265358979d0)
!
! and the local variables
!
integer :: isymq, sna, snb, snc, irot, na, nb, nc, ipol, jpol, &
lpol, kpol, mpol, npol
! counter on symmetries
! the rotated of the a atom
! the rotated of the b atom
! the rotated of the b atom
! counter on rotations
! counter on atoms
! counter on atoms
! counter on atoms
! counter on polarizations
! counter on polarizations
! counter on polarizations
! counter on polarizations
! counter on polarizations
! counter on polarizations
! counters
integer, allocatable:: iflb (:,:,:)
! used to account for symmetrized elements
......@@ -71,8 +53,8 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
! the argument of the phase
complex (kind = dp), allocatable :: phip (:,:,:,:,:,:)
! work space
complex (kind = dp) :: work (3, 3, 3), fase, faseq (48)
! working space
! the phase factor
! the phases for each symmetry
......@@ -85,10 +67,11 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
do kpol = 1, 3
do ipol = 1, 3
do jpol = 1, 3
phi (kpol, ipol, jpol, nc, na, nb) = 0.5d0 * (phi (kpol, ipol, &
jpol, nc, na, nb) + conjg (phi (kpol, jpol, ipol, nc, nb, na) ) )
phi (kpol, jpol, ipol, nc, nb, na) = conjg (phi (kpol, ipol, jpol, &
nc, na, nb) )
phi (kpol, ipol, jpol, nc, na, nb) = 0.5d0 * &
(phi (kpol, ipol, jpol, nc, na, nb) + &
conjg (phi (kpol, jpol, ipol, nc, nb, na) ) )
phi (kpol, jpol, ipol, nc, nb, na) = &
conjg (phi (kpol, ipol, jpol, nc, na, nb) )
enddo
enddo
enddo
......@@ -99,7 +82,7 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
! If no other symmetry is present we quit here
!
if ( (nsymq.eq.1) .and. (.not.minus_q) ) return
if ( (nsymq == 1) .and. (.not.minus_q) ) return
allocate (phip( 3, 3, 3, nat, nat, nat))
!
! Then we impose the symmetry q -> -q+G if present
......@@ -111,42 +94,47 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
do mpol = 1, 3
do ipol = 1, 3
do jpol = 1, 3
call setv (54, 0.d0, work, 1)
work = (0.d0, 0.d0)
snc = irt (irotmq, nc)
sna = irt (irotmq, na)
snb = irt (irotmq, nb)
arg = 0.d0
do kpol = 1, 3
arg = arg + (xq (kpol) * (rtau (kpol, irotmq, na) - rtau (kpol, &
irotmq, nb) ) )
arg = arg + (xq (kpol) * (rtau (kpol, irotmq, na) - &
rtau (kpol, irotmq, nb) ) )
enddo
arg = arg * tpi
fase = DCMPLX (cos (arg), sin (arg) )
do npol = 1, 3
do kpol = 1, 3
do lpol = 1, 3
work (mpol, ipol, jpol) = work (mpol, ipol, jpol) + fase * s ( &
ipol, kpol, irotmq) * s (jpol, lpol, irotmq) * s (mpol, npol, &
irotmq) * phi (npol, kpol, lpol, snc, sna, snb)
work (mpol, ipol, jpol) = work (mpol, ipol, jpol) + &
fase * s (ipol, kpol, irotmq) * &
s (jpol, lpol, irotmq) * &
s (mpol, npol, irotmq) * &
phi (npol, kpol, lpol, snc, sna, snb)
enddo
enddo
enddo
phip (mpol, ipol, jpol, nc, na, nb) = (phi (mpol, ipol, jpol, &
nc, na, nb) + conjg (work (mpol, ipol, jpol) ) ) * 0.5d0
phip (mpol, ipol, jpol, nc, na, nb) = &
(phi (mpol, ipol, jpol, nc, na, nb) + &
conjg (work (mpol, ipol, jpol) ) ) * 0.5d0
enddo
enddo
enddo
enddo
enddo
enddo
call ZCOPY (27 * nat * nat * nat, phip, 1, phi, 1)
phi = phip
endif
deallocate (phip)
!
! Here we symmetrize with respect to the small group of q
!
if (nsymq.eq.1) return
if (nsymq == 1) return
allocate (iflb( nat, nat, nat))
do na = 1, nat
......@@ -155,13 +143,13 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
iflb (nc, na, nb) = 0
enddo
enddo
enddo
do nc = 1, nat
do na = 1, nat
do nb = 1, nat
if (iflb (nc, na, nb) .eq.0) then
call setv (54, 0.d0, work, 1)
work = (0.d0, 0.d0)
do isymq = 1, nsymq
irot = irgq (isymq)
snc = irt (irot, nc)
......@@ -169,8 +157,8 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
snb = irt (irot, nb)
arg = 0.d0
do ipol = 1, 3
arg = arg + (xq (ipol) * (rtau (ipol, irot, na) - rtau (ipol, &
irot, nb) ) )
arg = arg + (xq (ipol) * (rtau (ipol, irot, na) - &
rtau (ipol, irot, nb) ) )
enddo
arg = arg * tpi
faseq (isymq) = DCMPLX (cos (arg), sin (arg) )
......@@ -180,9 +168,12 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
do npol = 1, 3
do kpol = 1, 3
do lpol = 1, 3
work (mpol, ipol, jpol) = work (mpol, ipol, jpol) + s (ipol, &
kpol, irot) * s (jpol, lpol, irot) * s (mpol, npol, irot) &
* phi (npol, kpol, lpol, snc, sna, snb) * faseq (isymq)
work (mpol, ipol, jpol) = work (mpol, ipol, jpol) + &
s (ipol, kpol, irot) * &
s (jpol, lpol, irot) * &
s (mpol, npol, irot) * &
phi (npol, kpol, lpol, snc, sna, snb) &
* faseq (isymq)
enddo
enddo
enddo
......@@ -202,10 +193,13 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
do npol = 1, 3
do kpol = 1, 3
do lpol = 1, 3
phi (mpol, ipol, jpol, snc, sna, snb) = phi (mpol, ipol, jpol, &
snc, sna, snb) + s (mpol, npol, invs (irot) ) * s (ipol, kpol, &
invs (irot) ) * s (jpol, lpol, invs (irot) ) * work (npol, &
kpol, lpol) * conjg (faseq (isymq) )
phi (mpol, ipol, jpol, snc, sna, snb) = &
phi (mpol, ipol, jpol, snc, sna, snb) +&
s (mpol, npol, invs (irot) ) * &
s (ipol, kpol, invs (irot) ) * &
s (jpol, lpol, invs (irot) ) * &
work (npol, kpol, lpol) * &
conjg (faseq (isymq) )
enddo
enddo
enddo
......@@ -218,7 +212,7 @@ subroutine d3_symdynph (xq, phi, s, invs, rtau, irt, irgq, nsymq, &
enddo
enddo
enddo
call DSCAL (54 * nat * nat * nat, 1.d0 / nsymq, phi, 1)
phi = phi / float(nsymq)
deallocate (iflb)
return
end subroutine d3_symdynph
......@@ -19,14 +19,14 @@ subroutine d3_valence
integer :: ik, ikk, ikq, nu_i, nu_j, nu_k, ibnd, jbnd, kbnd, nrec
real (kind = dp) :: de1, de2, de3, wg1, wg2, wg3, wwg1, wwg2, d_dos, wrk, &
wgauss, wga (nbnd), wgq (nbnd), w0gauss, w0g (nbnd), w1g (nbnd), &
w_1gauss
wga (nbnd), wgq (nbnd), w0g (nbnd), w1g (nbnd)
real (kind = dp), external :: wgauss, w0gauss, w_1gauss
complex (kind = dp) :: wrk1, aux (3 * nat)
complex (kind = dp), allocatable :: pdvp_i (:,:), pdvp_j (:,:), dpsidvpsi (:,:), &
pdvp_k (:,:), aux1 (:,:,:), aux2 (:,:,:), aux3 (:,:,:), aux4 (:,:,:)
if (degauss.eq.0.d0) return
if (degauss == 0.d0) return