!-----------------------------------------------------------------------
      subroutine aainit(lli,ap,lpx,lpl)
!-----------------------------------------------------------------------
!
      use van_parameters
      use cnst
!
      implicit none
! 
!     25 = combined angular momentum for 2*lli-1 (for s,p and d states)
!
      integer lli
      real(kind=8) ap(25,nlx,nlx)
      integer lpx(nlx,nlx), lpl(nlx,nlx,mx)
! local variables
      integer k, l, ll, li, mi, lj, mj, lp, mp, np, m, n, p, ik, il, ii,&
     &        mpp, ilp
      real(kind=8)    cc(lix,lx,lix,lx,lx), f, y, t, s, fs, ss, fts, ts
      complex(kind=8) aa(lx,mx,lix,lx,lix,lx), u(lx,mx,mx), sum
!
!      cc = < lm | limi | ljmj > = sqrt(2l+1/4pi) (-1)^mj c^l(li-mi,ljmj)
!
!      coeff. c^l from weissbluth 'atoms and molec..'  page 246
!      cc(li,mi,lj,mj,l)  l  is free index, m=mj+mj   (input limi ljmj)
!
      ll=2*lli-1
      write(6,*)
      write(6,*) ' aainit:  lli ll ',lli,ll
      write(6,*)
      if (lli.gt.lix) call error(' aainit ',' lli .gt. lix ',lli)
!
      call zero(lix*lx*lix*lx*lx,cc)
!
!      limi  ljmj  -  lm    corresponding to each cc are given
!
!      s     s     -   s
       cc(1,1,1,1,1)= 1. 
!
      if(lli.gt.1) then
!      s    p_-1   -  p_-1
!      s    p_0    -  p_0
!      s    p_1    -  p_1
       cc(1,1,2,1,2)= 1.
       cc(1,1,2,2,2)= 1.
       cc(1,1,2,3,2)= 1.
! 
!      p_-1  p_-1   -  d_-2
!      p_-1  p_0    -  d_-1
!      p_-1  p_1    -   s
!      p_-1  p_1    -  d_0
!      p_0   p_0    -   s
!      p_0   p_0    -  d_0
!      p_0   p_1    -  d_1
!      p_1   p_1    -  d_2
       cc(2,1,2,1,3)= sqrt(6./5.)
       cc(2,1,2,2,3)= sqrt(3./5.)
       cc(2,1,2,3,1)=-1.         
       cc(2,1,2,3,3)= sqrt(1./5.)
       cc(2,2,2,2,1)= 1.         
       cc(2,2,2,2,3)= sqrt(4./5.)
       cc(2,2,2,3,3)= sqrt(3./5.)
       cc(2,3,2,3,3)= sqrt(6./5.)
      endif
! 
      if(lli.gt.2) then
!      s    d_-2   -  d_-2
!      s    d_-1   -  d_-1
!      s    d_0    -  d_0
!      s    d_1    -  d_1
!      s    d_2    -  d_2
       cc(1,1,3,1,3)= 1.
       cc(1,1,3,2,3)= 1.
       cc(1,1,3,3,3)= 1.
       cc(1,1,3,4,3)= 1.
       cc(1,1,3,5,3)= 1.
! 
!      p_-1  d_-2  -  f_-3
!      p_-1  d_-1  -  f_-2
!      p_-1  d_0   -  p_-1
!      p_-1  d_0   -  f_-1
!      p_-1  d_1   -  p_0
!      p_-1  d_1   -  f_0
!      p_-1  d_2   -  p_1
!      p_-1  d_2   -  f_1
       cc(2,1,3,1,4)= sqrt(45./35.)
       cc(2,1,3,2,4)= sqrt(30./35.)
       cc(2,1,3,3,2)=-sqrt( 1./ 5.)
       cc(2,1,3,3,4)= sqrt(18./35.)
       cc(2,1,3,4,2)=-sqrt( 3./ 5.)
       cc(2,1,3,4,4)= sqrt( 9./35.)
       cc(2,1,3,5,2)=-sqrt( 6./ 5.)
       cc(2,1,3,5,4)= sqrt( 3./35.)
!
!      p_0   d_-2  -  f_-2
!      p_0   d_-1  -  p_-1
!      p_0   d_-1  -  f_-1
!      p_0   d_0   -  p_0
!      p_0   d_0   -  f_0
!      p_0   d_1   -  p_1
!      p_0   d_1   -  f_1
!      p_0   d_2   -  f_2
       cc(2,2,3,1,4)= sqrt(15./35.)
       cc(2,2,3,2,2)= sqrt( 3./ 5.)
       cc(2,2,3,2,4)= sqrt(24./35.)
       cc(2,2,3,3,2)= sqrt( 4./ 5.)
       cc(2,2,3,3,4)= sqrt(27./35.)
       cc(2,2,3,4,2)= sqrt( 3./ 5.)
       cc(2,2,3,4,4)= sqrt(24./35.)
       cc(2,2,3,5,4)= sqrt(15./35.)
!
!      p_1   d_-2  -  p_-1
!      p_1   d_-2  -  f_-1
!      p_1   d_-1  -  p_0
!      p_1   d_-1  -  f_0
!      p_1   d_0   -  p_1
!      p_1   d_0   -  f_1
!      p_1   d_1   -  f_2
!      p_1   d_2   -  f_3
       cc(2,3,3,1,2)=-sqrt( 6./ 5.)
       cc(2,3,3,1,4)= sqrt( 3./35.)
       cc(2,3,3,2,2)=-sqrt( 3./ 5.)
       cc(2,3,3,2,4)= sqrt( 9./35.)
       cc(2,3,3,3,2)=-sqrt( 1./ 5.)
       cc(2,3,3,3,4)= sqrt(18./35.)
       cc(2,3,3,4,4)= sqrt(30./35.)
       cc(2,3,3,5,4)= sqrt(45./35.)
!
!      d_-2  d_-2  -  g_-4
!      d_-2  d_-1  -  g_-3
!      d_-2  d_0   -  d_-2
!      d_-2  d_0   -  g_-2
!      d_-2  d_1   -  d_1
!      d_-2  d_1   -  g_1
!      d_-2  d_2   -  s
!      d_-2  d_2   -  d_0
!      d_-2  d_2   -  g_0
       cc(3,1,3,1,5)= sqrt(70./49.)
       cc(3,1,3,2,5)= sqrt(35./49.)
       cc(3,1,3,3,3)=-sqrt(20./49.)
       cc(3,1,3,3,5)= sqrt(15./49.)
       cc(3,1,3,4,3)=-sqrt(30./49.)
       cc(3,1,3,4,5)= sqrt( 5./49.)
       cc(3,1,3,5,1)= 1.           
       cc(3,1,3,5,3)=-sqrt(20./49.)
       cc(3,1,3,5,5)= sqrt( 1./49.)
! 
!      d_-1  d_-1  -  d_-2
!      d_-1  d_-1  -  g_-2
!      d_-1  d_0   -  d_-1
!      d_-1  d_0   -  g_-1
!      d_-1  d_1   -  s
!      d_-1  d_1   -  d_0
!      d_-1  d_1   -  g_0
!      d_-1  d_2   -  d_1
!      d_-1  d_2   -  g_1
       cc(3,2,3,2,3)= sqrt(30./49.)
       cc(3,2,3,2,5)= sqrt(40./49.)
       cc(3,2,3,3,3)= sqrt( 5./49.)
       cc(3,2,3,3,5)= sqrt(30./49.)
       cc(3,2,3,4,1)=-1.           
       cc(3,2,3,4,3)=-sqrt( 5./49.)
       cc(3,2,3,4,5)= sqrt(16./49.)
       cc(3,2,3,5,3)=-sqrt(30./49.)
       cc(3,2,3,5,5)= sqrt( 5./49.)
!
!      d_0   d_0   -  s
!      d_0   d_0   -  d_0
!      d_0   d_0   -  g_0
!      d_0   d_1   -  d_1
!      d_0   d_1   -  g_1
!      d_0   d_2   -  d_2
!      d_0   d_2   -  g_2
       cc(3,3,3,3,1)= 1.
       cc(3,3,3,3,3)= sqrt(20./49.)
       cc(3,3,3,3,5)= sqrt(36./49.)
       cc(3,3,3,4,3)= sqrt( 5./49.)
       cc(3,3,3,4,5)= sqrt(30./49.)
       cc(3,3,3,5,3)=-sqrt(20./49.)
       cc(3,3,3,5,5)= sqrt(15./49.)
!
!      d_1   d_1   -  d_2
!      d_1   d_1   -  g_2
!      d_1   d_2   -  g_3
       cc(3,4,3,4,3)= sqrt(30./49.)
       cc(3,4,3,4,5)= sqrt(40./49.)
       cc(3,4,3,5,5)= sqrt(35./49.)
!
!      d_2   d_2   -  g_4
       cc(3,5,3,5,5)= sqrt(70./49.)
      endif
!
      do l=1,ll
         do li=1,lli
            do mi=1,2*li-1      
!     symmetry   limi  <-> ljmj
               do lj=1,li-1
                  do mj=1,2*lj-1
                     cc(li,mi,lj,mj,l)=cc(lj,mj,li,mi,l)
                  end do
               end do
!
               do mj=1,mi-1
                  cc(li,mi,li,mj,l)=cc(li,mj,li,mi,l)      
               end do
!     li = lj
            end do
         end do
      end do
!
      call SCAL(lix*lx*lix*lx*lx,sqrt(1./fpi),cc,1)
!
!     transform between real spherical harmonics and the original ones
!                (  y^r_lm = sum_n u(l,m,n) y_ln )
!
!          u(l,m,n)  see  weissbluth  pages 128 - 130
!          errors of weissbluth have been corrected:
!          1.)   i/2 has been changed to i/4 for u(4,6,*)      
!          2.)  -i/4 has been changed to i/4 for u(5,8,*)
!
      call zero(2*lx*mx*mx,u)
!     l=0
      u(1,1,1)=cmplx(1.,0.)
!
      if(lli.gt.1) then
!      l = 1    ordering   x   z   y
!                   m     -1   0   1
!                   m      1   2   3
         y=1./sqrt(2.)
         u(2,1,1)=cmplx( y,0.) 
         u(2,1,3)=cmplx(-y,0.)
         u(2,2,2)=cmplx(1.,0.)
         u(2,3,1)=cmplx(0., y)
         u(2,3,3)=cmplx(0., y)
!
!       l = 2
!       ordering  xy   xz  z^2  yz  x^2-y^2
!          m      -2   -1   0    1     2
!          m       1    2   3    4     5
         u(3,1,1)=cmplx(0., y)
         u(3,1,5)=cmplx(0.,-y)
         u(3,2,2)=cmplx( y,0.)
         u(3,2,4)=cmplx(-y,0.)
         u(3,3,3)=cmplx(1.,0.)
         u(3,4,2)=cmplx(0., y)
         u(3,4,4)=cmplx(0., y)
         u(3,5,1)=cmplx( y,0.)
         u(3,5,5)=cmplx( y,0.)
      endif
!
      if(lli.gt.2) then
!      l = 3
!      order
!      x  y  xyz  z  z(x^2-y^2) y(z^2-x^2) x(y^2-z^2)
!      1  2   3   4       5          6          7
         f=sqrt(5.)/4.
         t=sqrt(3.)/4.
         u(4,1,1)=cmplx(f ,0.)
         u(4,1,3)=cmplx(-t,0.)
         u(4,1,5)=cmplx(t ,0.)
         u(4,1,7)=cmplx(-f,0.)
         u(4,2,1)=cmplx(0.,-f)
         u(4,2,3)=cmplx(0.,-t)
         u(4,2,5)=cmplx(0.,-t)
         u(4,2,7)=cmplx(0.,-f)
         u(4,3,2)=cmplx(0., y)
         u(4,3,6)=cmplx(0.,-y)
         u(4,4,4)=cmplx(1.,0.)
         u(4,5,2)=cmplx( y,0.) 
         u(4,5,6)=cmplx( y,0.) 
!     
!     the following 4 are missprinted in weissbluth (page 129)
         u(4,6,1)=cmplx(0.,-t)
         u(4,6,3)=cmplx(0., f)
         u(4,6,5)=cmplx(0., f)
         u(4,6,7)=cmplx(0.,-t)
!
         u(4,7,1)=cmplx(-t,0.)
         u(4,7,3)=cmplx(-f,0.)
         u(4,7,5)=cmplx( f,0.)
         u(4,7,7)=cmplx( t,0.)
! 
!      l = 4
         s  =sqrt( 7.   )/4.0
         fs =sqrt( 5./6.)/2.0
         ss =sqrt( 7./6.)/2.0
         fts=sqrt(14./6.)/2.0
         ts =sqrt(10./6.)/2.0
         u(5,1,1)=cmplx(fs,0.)       
         u(5,1,5)=cmplx(fts,0.)
         u(5,1,9)=cmplx(fs,0.)
         u(5,2,2)=cmplx(0.,-0.25)
         u(5,2,4)=cmplx(0.,-s)
         u(5,2,6)=cmplx(0.,-s)
         u(5,2,8)=cmplx(0.,-0.25)
         u(5,3,2)=cmplx(-0.25,0.)
         u(5,3,4)=cmplx( s,0.)
         u(5,3,6)=cmplx(-s,0.)
         u(5,3,8)=cmplx( 0.25,0.)
         u(5,4,3)=cmplx(-y,0.)
         u(5,4,7)=cmplx(-y,0.)
         u(5,5,1)=cmplx(0., y)
         u(5,5,9)=cmplx(0.,-y)
         u(5,6,3)=cmplx(0., y)
         u(5,6,7)=cmplx(0.,-y)
         u(5,7,2)=cmplx(-s,0.)
         u(5,7,4)=cmplx(-0.25,0.)
         u(5,7,6)=cmplx( 0.25,0.)
         u(5,7,8)=cmplx( s,0.)
!
!      the following 4 are missprinted in weissbluth (page 129)
         u(5,8,2)=cmplx(0., s)   
         u(5,8,4)=cmplx(0.,-0.25)
         u(5,8,6)=cmplx(0.,-0.25)
         u(5,8,8)=cmplx(0., s)   
!     
         u(5,9,1)=cmplx(-ss,0.)
         u(5,9,5)=cmplx( ts,0.)
         u(5,9,9)=cmplx(-ss,0.)
      endif
!
!
      call zero(2*lix*lx*lix*lx*lx*mx,aa)
      do lp=1,ll
         do l=1,lli
            do m=1,(2*l-1)
               do k=1,lli
                  do n=1,(2*k-1)
                     do np=1,(2*k-1)
                        do p=1,(2*l-1)
                           mp=(p-l)+(np-k)+lp
!         m'= p   + n'
                           if((mp.ge.1).and.(mp.le.(2*lp-1))) then
                              aa(lp,mp,l,m,k,n)=aa(lp,mp,l,m,k,n)+      &
     &                          u(l,m,p)*u(k,n,np)*cc(l,p,k,np,lp)
                           endif
                        end do
                     end do
                  end do
               end do
            end do
         end do
      end do
!
      do ik=1,nlx
         do il=1,nlx
            do ii=1,mx
               lpl(il,ik,ii)=0
            end do
            lpx(il,ik)=0
         end do
      end do
      do lp=1,ll
         do mp=1,(2*lp-1)
            do l=1,lli
               do m=1,(2*l-1)
                  do k=1,lli
                     do n=1,(2*k-1)
                        sum=cmplx(0.,0.)
                        do mpp=1,(2*lp-1)
                           sum = sum +                                  &
     &                          conjg(u(lp,mp,mpp))*aa(lp,mpp,l,m,k,n)
                        end do
                        if(abs(sum).gt.0.001) then
                           il =(l-1)*(l-1)+m
                           ik =(k-1)*(k-1)+n
                           ilp=(lp-1)*(lp-1)+mp
                           if(abs(aimag(sum)).gt.0.001) then
                              write(6,'(a,2f8.5,1x,3i3,2x,i3)')         &
     &                        'ap(ilp,il,ik) =  ', sum,ilp,il,ik
                              call error('aainit','see above',1)
                           endif
                           ap(ilp,il,ik)=real(sum)
                           lpx(il,ik)=lpx(il,ik)+1
                           lpl(il,ik,lpx(il,ik))=ilp
                        endif
                     end do
                  end do
               end do
            end do
         end do
      end do
!
      return
      end
!-----------------------------------------------------------------------
      subroutine bess(xg,l,mmax,r,jl)
!-----------------------------------------------------------------------
!     calculates spherical bessel functions j_l(qr)
!     NOTA BENE: it is assumed that r(1)=0 always
!
      implicit none
      integer l, mmax
      real(kind=8)    xg, jl(mmax), r(mmax)
! local variables
      real(kind=8)    eps, xrg, xrg2
      parameter(eps=1.e-8)
      integer i, ir
!
!    l=-1 (for derivative  calculations)
!
      if(l.eq.0) then
         if(xg.lt.eps) then 
            do i=1,mmax
               jl(i)=0.0
            end do
         else
            jl(1)=0.
            do ir=2,mmax
               xrg=r(ir)*xg
               jl(ir)=cos(xrg)/xrg
            end do
         end if
      end if
!
!    s part
!
      if(l.eq.1) then
         if(xg.lt.eps) then 
            do i=1,mmax
               jl(i)=1.0
            end do
         else
            jl(1)=1.
            do ir=2,mmax
               xrg=r(ir)*xg
               jl(ir)=sin(xrg)/xrg
            end do
         endif
      endif
!
!     p-part
! 
      if(l.eq.2) then
         if(xg.lt.eps) then
            do i=1,mmax
               jl(i)=0.0
            end do
         else
            jl(1)=0.
            do ir=2,mmax
               xrg=r(ir)*xg
               jl(ir)=(sin(xrg)/xrg-cos(xrg))/xrg
            end do
         endif
      endif
!
!     d part
! 
      if(l.eq.3) then
         if(xg.lt.eps) then
            do i=1,mmax
               jl(i)=0.0
            end do
         else
            jl(1)=0.
            do ir=2,mmax
               xrg=r(ir)*xg
               jl(ir)=(sin(xrg)*(3./(xrg*xrg)-1.)                       &
     &              -3.*cos(xrg)/xrg) /xrg
            end do
         endif
      endif
!
!     f part
!
      if(l.eq.4) then
         if(xg.lt.eps) then
            do i=1,mmax
               jl(i)=0.0
            end do
         else
            jl(1)=0.
            do ir=2,mmax
               xrg=r(ir)*xg
               xrg2=xrg*xrg
               jl(ir)=( sin(xrg)*(15./(xrg2*xrg)-6./xrg)                &
     &              +cos(xrg)*(1.-15./xrg2)           )/xrg
            end do
         endif
      endif
!
!     g part
!
      if(l.eq.5) then
         if(xg.lt.eps) then
            do i=1,mmax
               jl(i)=0.0
            end do
         else
            jl(1)=0.
            do ir=2,mmax
               xrg=r(ir)*xg
               xrg2=xrg*xrg
               jl(ir)=( sin(xrg)*(105./(xrg2*xrg2)-45./xrg2+1.)         &
     &              +cos(xrg)*(10./xrg-105./(xrg2*xrg)) )/xrg
            end do
         endif
      endif
!
      return
      end
!-----------------------------------------------------------------------
      subroutine calbec (nspmn,nspmx,eigr,c,bec)
!-----------------------------------------------------------------------
!     this routine calculates array bec
!
!        < psi_n | beta_i,i > = c_n(0) beta_i,i(0) +
!                 2 sum_g> re(c_n*(g) (-i)**l beta_i,i(g) e^-ig.r_i)
!
!     routine makes use of c(-g)=c*(g)  and  beta(-g)=beta*(g)
!
      use ions
      use elct
      use control
      use cvan
!
      implicit none
      integer nspmn, nspmx
      real(kind=8)  bec(nhsa,n)
      complex(kind=8) c(ngw,n), eigr(ngw,nas,nspmx)
! local variables
      integer is, ia, i , iv
!
!
      call tictac(25,0)
      call nlsm1(n,nspmn,nspmx,eigr,c,bec)
!
      if (iprsta.gt.2) then
         write(6,*)
         do is=1,nspmx
            if(nspmx.gt.1) then
               write(6,'(33x,a,i4)') ' calbec: bec (is)',is
               write(6,'(8f9.4)')                                       &
     &              ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
            else
               do ia=1,na(is)
                  write(6,'(33x,a,i4)') ' calbec: bec (ia)',ia
                  write(6,'(8f9.4)')                                    &
     &             ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
               end do
            end if
         end do
      endif
      call tictac(25,1)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine calphi(c0,ema0bg,bec,betae,phi)
!-----------------------------------------------------------------------
!     input: c0 (orthonormal with s(r(t)), bec=<c0|beta>, betae=|beta>
!     computes the matrix phi (with the old positions)
!       where  |phi> = s'|c0> = |c0> + sum q_ij |i><j|c0>
!     where s'=s(r(t))  
!
      use ions
      use cvan
      use elct
      use parm
      use cnst
      use control
!
      implicit none
      complex(kind=8) c0(ngw,n), phi(ngw,n), betae(ngw,nhsa)
      real(kind=8)    ema0bg(ngw), bec(nhsa,n), emtot
! local variables
      integer is, iv, jv, ia, inl, jnl, i, j
      real(kind=8) qtemp(nhsavb,n) ! automatic array
!
      call tictac(9,0)
      call zero(2*ngw*n,phi)
!
      if (nvb.gt.0) then
         call zero(nhsavb*n,qtemp)
         do is=1,nvb
            do iv=1,nh(is)
               do jv=1,nh(is)
                  if(abs(qq(iv,jv,is)).gt.1.e-5) then
                     do ia=1,na(is)
                        inl=ish(is)+(iv-1)*na(is)+ia
                        jnl=ish(is)+(jv-1)*na(is)+ia
                        do i=1,n
                           qtemp(inl,i) = qtemp(inl,i) +                &
     &                                    qq(iv,jv,is)*bec(jnl,i)
                        end do
                     end do
                  endif
               end do
            end do
         end do
!
         call MXMA                                                     &
     &       (betae,1,2*ngw,qtemp,1,nhsavb,phi,1,2*ngw,2*ngw,nhsavb,n)
      end if
!
      do j=1,n
         do i=1,ngw
            phi(i,j)=(phi(i,j)+c0(i,j))*ema0bg(i)
         end do
      end do
!     =================================================================
      if(iprsta.gt.2) then
         emtot=0.
         do j=1,n
            do i=1,ngw
               emtot=emtot                                              &
     &        +2.*real(phi(i,j)*conjg(c0(i,j)))*ema0bg(i)**(-2.)
            end do
         end do
         emtot=emtot/n
#ifdef PARA
         call reduce(1,emtot)
#endif      
         write(6,*) 'in calphi sqrt(emtot)=',sqrt(emtot)
         write(6,*)
         do is=1,nsp
            if(nsp.gt.1) then
               write(6,'(33x,a,i4)') ' calphi: bec (is)',is
               write(6,'(8f9.4)')                                       &
     &            ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
            else
               do ia=1,na(is)
                  write(6,'(33x,a,i4)') ' calphi: bec (ia)',ia
                  write(6,'(8f9.4)')                                    &
     &               ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
               end do
            end if
         end do
      endif
      call tictac(9,1)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine cofmass(tau,cdm)
!-----------------------------------------------------------------------
!
      use ions
!
      implicit none
      real(kind=8) tau(3,nax,nsp), cdm(3)
! local variables
      real(kind=8) tmas
      integer is,i,ia
!
      tmas=0.0
      do is=1,nsp
         tmas=tmas+na(is)*pmass(is)
      end do
!
      do i=1,3
         cdm(i)=0.0
         do is=1,nsp
            do ia=1,na(is)
               cdm(i)=cdm(i)+tau(i,ia,is)*pmass(is)
            end do
         end do
         cdm(i)=cdm(i)/tmas
      end do
!
      return
      end
!-----------------------------------------------------------------------
      real(kind=8) function cscnorm(bec,cp,i)
!-----------------------------------------------------------------------
!     requires in input the updated bec(i)
!
      use ions
      use elct
      use cvan
!
      implicit none
      integer i
      real(kind=8) bec(nhsa,n)
      complex(kind=8) cp(ngw,n)
!
      integer ig, is, iv, jv, ia, inl, jnl
      real(kind=8) sum, SSUM
      real(kind=8), allocatable:: temp(:)
!
!
      allocate(temp(ngw))
      do ig=1,ngw
         temp(ig)=real(conjg(cp(ig,i))*cp(ig,i))
      end do
      sum=2.*SSUM(ngw,temp,1)
      if (ng0.eq.2) sum=sum-temp(1)
#ifdef PARA
      call reduce(1,sum)
#endif
      do is=1,nvb
         do iv=1,nh(is)
            do jv=1,nh(is)
               if(abs(qq(iv,jv,is)).gt.1.e-5) then 
                  do ia=1,na(is)
                     inl=ish(is)+(iv-1)*na(is)+ia
                     jnl=ish(is)+(jv-1)*na(is)+ia
                     sum = sum +                                        &
     &                    qq(iv,jv,is)*bec(inl,i)*bec(jnl,i)
                  end do
               endif
            end do
         end do
      end do
!
      cscnorm=sqrt(sum)
!
      deallocate(temp)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine dforce (bec,deeq,betae,i,c,ca,df,da,v,tbuff)
!-----------------------------------------------------------------------
!omputes: the generalized force df=cmplx(dfr,dfi) acting on the i-th
!         electron state at the gamma point of the brillouin zone
!         represented by the vector c=cmplx(cr,ci)
!
!     d_n(g) = f_n { 0.5 g^2 c_n(g) + [vc_n](g) +
!              sum_i,ij d^q_i,ij (-i)**l beta_i,i(g) 
!                                 e^-ig.r_i < beta_i,j | c_n >}
#ifdef SOUNCORNO
      use pres_mod
#endif
      use gvec
      use gvecs
      use cvan
      use parm
      use parms
      use elct
      use cnst
      use ions
      use work1
!
      implicit none
!
      complex(kind=8) betae(ngw,nhsa), c(ngw), ca(ngw), df(ngw), da(ngw)
      real(kind=8) bec(nhsa,n), deeq(nat,nhx,nhx,nspin), v(nnrs,nspin)
      integer i
      logical tbuff
! local variables
      integer iv, jv, ia, is, isa, ism, ios, iss1, iss2, ir, ig, inl, jnl
      real(kind=8) fi, fip, dd
      complex(kind=8) fp,fm,ci
      real(kind=8) af(nhsa), aa(nhsa) ! automatic arrays
      complex(kind=8)  dtemp(ngw)    !
      complex(kind=8), pointer:: psi(:)
!
!
      psi => wrk1
!
!     important: if n is odd => c(*,n+1)=0.
! 
      if (mod(n,2).ne.0.and.i.eq.n) then
         do ig=1,ngw
            ca(ig)=(0.,0.)
         end do
      endif
!
      call tictac(18,0)
      ci=(0.0,1.0)
!
      if (.not.tbuff) then
!
         call zero(2*nnrs,psi)
         do ig=1,ngw
            psi(nms(ig))=conjg(c(ig)-ci*ca(ig))
            psi(nps(ig))=c(ig)+ci*ca(ig)
         end do
!
         call ivfftw(psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!     
      else
!
!     read psi from buffer 21
!
!!!!#if defined(CRAYY) || defined(NEC)
! A.Bongiorno
#if defined(CRAYY)
         buffer in(21,0) (psi(1),psi(nnrs))
         ios = unit(21)
#else
         read(21,iostat=ios) psi
#endif
         if(ios.ne.0) call error                                        &
     &       (' dforce',' error in reading unit 21',ios)
!
      endif
! 
      iss1=ispin(i)
!
! the following avoids a potential out-of-bounds error
!
      if (i.ne.n) then
         iss2=ispin(i+1)
      else
         iss2=iss1
      end if
!
      do ir=1,nnrs
         psi(ir)=cmplx(v(ir,iss1)* real(psi(ir)),                       &
     &                 v(ir,iss2)*aimag(psi(ir)) )
      end do
!
      call fwfftw(psi,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!
!     note : the factor 0.5 appears 
!       in the kinetic energy because it is defined as 0.5*g**2
!       in the potential part because of the logics
!
      fi =-  f(i)*0.5
      fip=-f(i+1)*0.5
      do ig=1,ngw
         fp= psi(nps(ig)) + psi(nms(ig))
         fm= psi(nps(ig)) - psi(nms(ig))
#ifdef SOUNCORNO
         df(ig)= fi*(tpiba2*ggp(ig)* c(ig)+cmplx(real(fp), aimag(fm)))
         da(ig)=fip*(tpiba2*ggp(ig)*ca(ig)+cmplx(aimag(fp),-real(fm)))
#else
         df(ig)= fi*(tpiba2*g(ig)* c(ig)+cmplx(real(fp), aimag(fm)))
         da(ig)=fip*(tpiba2*g(ig)*ca(ig)+cmplx(aimag(fp),-real(fm)))
#endif
      end do
      call tictac(18,1)
!
      call tictac(23,0)
!
!     aa_i,i,n = sum_j d_i,ij <beta_i,j|c_n>
! 
      if(nhsa.gt.0)then
         do inl=1,nhsa
            af(inl)=0.
            aa(inl)=0.
         end do
!
         do is=1,nsp
            do iv=1,nh(is)
               do jv=1,nh(is)
                  isa=0
                  do ism=1,is-1
                     isa=isa+na(ism)
                  end do
                  do ia=1,na(is)
                     inl=ish(is)+(iv-1)*na(is)+ia
                     jnl=ish(is)+(jv-1)*na(is)+ia
                     isa=isa+1
                     dd = deeq(isa,iv,jv,iss1)+dvan(iv,jv,is)
                     af(inl)=af(inl)-  f(i)*dd*bec(jnl,  i)
                     dd = deeq(isa,iv,jv,iss2)+dvan(iv,jv,is)
                     if (i.ne.n) aa(inl)=aa(inl)-f(i+1)*dd*bec(jnl,i+1)
                  end do
               end do
            end do
         end do
!
         do ig=1,ngw
            dtemp(ig)=(0.,0.)
         end do
         call MXMA                                                      &
     &        (betae,1,2*ngw,af,1,nhsa,dtemp,1,2*ngw,2*ngw,nhsa,1)
         do ig=1,ngw
            df(ig)=df(ig)+dtemp(ig)
         end do
!
         do ig=1,ngw
            dtemp(ig)=(0.,0.)
         end do
         call MXMA                                                      &
     &        (betae,1,2*ngw,aa,1,nhsa,dtemp,1,2*ngw,2*ngw,nhsa,1)
         do ig=1,ngw
            da(ig)=da(ig)+dtemp(ig)
         end do
      endif
!
      call tictac(23,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine dftname(xctype)
!-----------------------------------------------------------------------
      use dft_mod
!
      implicit none
      character*20 xctype
!
                        xctype = ' unknown exch.-corr.'
      if (dft.eq.lda)   xctype = '      ceperley-alder'
      if (dft.eq.-1.)   xctype = '              wigner'
      if (dft.eq.-2.)   xctype = '     hedin-lundqvist'
      if (dft.eq.-3.)   xctype = ' gunnarson-lundqvist'
      if (dft.eq.-4.)   xctype = ' no exch-corr at all'
      if (dft.eq.blyp)  xctype = ' C-A + B88gx + LYPgc'
      if (dft.eq.becke) xctype = ' C-A + B88gx        '
      if (dft.eq.bp88)  xctype = ' C-A + B88gx + P86gc'
      if (dft.eq.pw91)  xctype = ' Perdew Wang 1991   '
      if (dft.eq.pbe)   xctype = ' PBE - GGA          '
!
      return
      end
!-----------------------------------------------------------------------
      subroutine displace(nsp,na,tau,taui,pmass,dis,discm)
!-----------------------------------------------------------------------
!
      use ion_parameters
      implicit none
!
      integer nsp, na(nsp)
      real(kind=8) tau(3,nax,nsp), taui(3,nax,nsp), pmass(nsp), dis, discm
! local variables
      integer i, k, j
      real(kind=8) cdm(3),cdmi(3),r2(3)
!
! 
      dis=0.0
      discm=0.0
      call cofmass(tau,cdm)
      call cofmass(taui,cdmi)
      do i=1,3
         r2(i)=0.0
         do k=1,nsp
            do j=1,na(k)
               r2(i)=r2(i)+(tau(i,j,k)-cdm(i)-(taui(i,j,k)-cdmi(i)))**2
            end do
            r2(i)=r2(i)/float(na(k))
         end do
         dis=dis+r2(i)
         discm=discm+(cdm(i)-cdmi(i))**2
      end do
      dis  =sqrt(dis  )
      discm=sqrt(discm) 
      return
      end
!-----------------------------------------------------------------------
      subroutine dotcsc(eigr,cp)
!-----------------------------------------------------------------------
!
      use ions
      use elct
      use cvan
!
      implicit none
!
      complex(kind=8)  eigr(ngw,nas,nsp), cp(ngw,n)
! local variables
      real(kind=8)  sum, csc(n) ! automatic array
      complex(kind=8) CSUM, temp(ngw) ! automatic array
 
      real(kind=8), allocatable::  becp(:,:)
      integer i,kmax,nnn,k,ig,is,ia,iv,jv,inl,jnl
!
      allocate(becp(nhsa,n))
!
!     < beta | phi > is real. only the i lowest:
!
      nnn=min(12,n)
      do i=nnn,1,-1
         kmax=i
         call nlsm1(i,1,nvb,eigr,cp,becp)
!
         do k=1,kmax
            do ig=1,ngw
               temp(ig)=conjg(cp(ig,k))*cp(ig,i)
            end do
            csc(k)=2.*real(CSUM(ngw,temp,1))
            if (ng0.eq.2) csc(k)=csc(k)-real(temp(1))
         end do
#ifdef PARA
         call reduce(kmax,csc)
#endif
         do k=1,kmax
            sum=0.
            do is=1,nvb
               do iv=1,nh(is)
                  do jv=1,nh(is)
                     do ia=1,na(is)
                        inl=ish(is)+(iv-1)*na(is)+ia
                        jnl=ish(is)+(jv-1)*na(is)+ia
                        sum = sum +                                      &
     &                   qq(iv,jv,is)*becp(inl,i)*becp(jnl,k)
                     end do
                  end do
               end do
            end do
            csc(k)=csc(k)+sum
         end do
!
         write(6,'(a,12f18.15)')' dotcsc = ',(csc(k),k=1,i)
!
      end do
      write(6,*)
!
      deallocate(becp)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine eigs(nspin,nx,nupdwn,iupdwn,f,lambda)
!-----------------------------------------------------------------------
!     computes: the eigenvalues of the real symmetric matrix lambda
!               the eigenvalues are printed out in electron volts.
!
      implicit none
      integer nspin, nx, nupdwn(nspin), iupdwn(nspin)
      real(kind=8) lambda(nx,nx), f(nx)
! local variables
      real(kind=8) lambdar(nx,nx), wr(nx), fv1(nx),fm1(2,nx) ! automatic arrays
      real(kind=8) zr,au
      integer iss,j,i,ierr
!
      au=27.212
!
      do iss=1,nspin
         do i=1,nupdwn(iss)
            do j=1,nupdwn(iss)
               lambdar(j,i)=lambda(iupdwn(iss)-1+j,iupdwn(iss)-1+i)
            end do
         end do
         call rs(nx,nupdwn(iss),lambdar,wr,0,zr,fv1,fm1,ierr)
         do i=1,nupdwn(iss)
            if (f(iupdwn(iss)-1+i).gt.1.e-6) then
               wr(i)=au*wr(i)/f(iupdwn(iss)-1+i)
            else
               wr(i)=0.0
            end if
         end do
!
!     print out eigenvalues
!
         write(6,12) 0., 0., 0.
         write(6,14) (wr(i),i=1,nupdwn(iss))
      end do
   12 format(//' eigenvalues at k-point: ',3f6.3)
!   14 format(10f8.2)
    14 format(8f11.5)
      write(6,*)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine error(a,b,n)
!-----------------------------------------------------------------------
      character*(*) a,b
#ifdef PARA
      include 'mpif.h'
#endif
!
      write(6,1) a,b,n
    1 format(//' program ',a,':',a,'.',8x,i8,8x,'stop')
#ifdef PARA
      call mpi_abort( MPI_COMM_WORLD, ierr, ierr)
      call mpi_finalize(ierr)
#endif
!
      stop 1 !ATTENZIONE   
      end
!----------------------------------------------------------------------
      subroutine expxc(nnr,nspin,rhor,exc)
!----------------------------------------------------------------------
!
!       ceperley & alder's correlation energy
!       after j.p. perdew & a. zunger prb 23, 5048 (1981)
!
!       rhor contains rho(r) on input, vxc(r) on output
!
      use cnst
!
      implicit none
!
      integer nnr, nspin
      real(kind=8) rhor(nnr,nspin), exc
! local variables
      integer ir, iflg, isup, isdw
      real(kind=8) roe, aroe, rs, rsl, rsq, ecca, vcca, eccp, vccp,    &
     &     zeta, onemz, zp, zm, fz, dfzdz, exc1, vxc1, vxc2
! constants
      real(kind=8) x76, x43, x13
      parameter(x76=7.d0/6.d0, x43=4.d0/3.d0, x13=1.d0/3.d0)
      real(kind=8) ax
      parameter (ax = -0.916330586d0)
! Perdew and Zunger parameters
      real(kind=8) ap, bp, cp, dp, af, bf, cf, df,                      &
     &     bp1, cp1, dp1, bf1, cf1, df1
      parameter                                                         &
     &( ap=0.03110*2.0, bp=-0.0480*2.0, cp=0.0020*2.0, dp=-0.0116*2.0   &
     &, af=0.01555*2.0, bf=-0.0269*2.0, cf=0.0007*2.0, df=-0.0048*2.0   &
     &, bp1=bp-ap/3.0, cp1=2.0*cp/3.0, dp1=(2.0*dp-cp)/3.0              &
     &, bf1=bf-af/3.0, cf1=2.0*cf/3.0, df1=(2.0*df-cf)/3.0 )
      real(kind=8) va(2), vb(2), vc(2), vd(2), vbt1(2), vbt2(2)
      real(kind=8)  a(2), b(2), c(2), d(2), g(2), b1(2), b2(2)
      data va/ap ,af /, vb/bp1,bf1/, vc/cp1,cf1/, vd/dp1,df1/,          &
     &     vbt1/1.0529,1.3981/, vbt2/0.3334,0.2611/
      data a/0.0622,0.0311/, b/-0.096,-0.0538/, c/0.0040,0.0014/,       &
     &     d/-0.0232,-0.0096/, b1/1.0529,1.3981/, b2/0.3334,0.2611/,    &
     &     g/-0.2846,-0.1686/
!
      if (nspin.eq.1) then
!
! iflg=1: paramagnetic (unpolarised) results
!
         iflg=1
         do ir=1,nnr
            roe=rhor(ir,1)
            if(roe.lt.1.0d-30) goto 10
            aroe=abs(roe)
            rs= (3.d0/aroe/fpi)**x13
            if(rs.le.1.d0) then
               rsl=log(rs)
               ecca= a(iflg)*rsl+ b(iflg)+ c(iflg)*rs*rsl+ d(iflg)*rs
               vcca=va(iflg)*rsl+vb(iflg)+vc(iflg)*rs*rsl+vd(iflg)*rs
            else
               rsq=sqrt(rs)
               ecca=g(iflg)/(1.d0+b1(iflg)*rsq+b2(iflg)*rs)
               vcca=ecca*(1.d0+x76*vbt1(iflg)*rsq+x43*vbt2(iflg)*rs)/   &
     &                   (1.d0+    vbt1(iflg)*rsq+    vbt2(iflg)*rs)
            end if
            exc1 = ( ax/rs + ecca )/2.
            exc = exc + exc1*roe
            rhor(ir,1)= ( x43*ax/rs + vcca )/2.
 10         continue
         end do
      else
         isup=1
         isdw=2
         do ir=1,nnr
            roe=rhor(ir,isup)+rhor(ir,isdw)
            if(roe.lt.1.0d-30) goto 20
            aroe=abs(roe)
            rs= (3.d0/aroe/fpi)**x13
            zeta=abs(rhor(ir,isup)-rhor(ir,isdw))/aroe
            zp = (1.d0+zeta)**x13
            onemz=max(0.d0,1.d0-zeta)
            zm = onemz**x13
            fz= ((1.d0+zeta)*zp + onemz*zm - 2.d0)/                     &
     &           (2.d0**x43 -2.d0)
            dfzdz= x43*(zp - zm)/(2.d0**x43-2.d0)
!
! iflg=1:  paramagnetic (unpolarised) results
! iflg=2: ferromagnetic (  polarised) results
!
            if(rs.le.1.d0) then
               rsl=log(rs)
               ecca= a(1)*rsl+ b(1)+ c(1)*rs*rsl+ d(1)*rs
               vcca=va(1)*rsl+vb(1)+vc(1)*rs*rsl+vd(1)*rs
               eccp= a(2)*rsl+ b(2)+ c(2)*rs*rsl+ d(2)*rs
               vccp=va(2)*rsl+vb(2)+vc(2)*rs*rsl+vd(2)*rs
            else
               rsq=sqrt(rs)
               ecca=g(1)/(1.d0+b1(1)*rsq+b2(1)*rs)
               vcca=ecca*(1.d0+x76*vbt1(1)*rsq+x43*vbt2(1)*rs)/         &
     &                   (1.d0+    vbt1(1)*rsq+    vbt2(1)*rs)
               eccp=g(2)/(1.d0+b1(2)*rsq+b2(2)*rs)
               vccp=eccp*(1.d0+x76*vbt1(2)*rsq+x43*vbt2(2)*rs)/         &
     &                   (1.d0+    vbt1(2)*rsq+    vbt2(2)*rs)
            end if
! exchange part
            exc1 = ax/rs*((1.d0+zeta)*zp+(1.d0-zeta)*zm)/2.
            vxc1 = x43*ax/rs*zp
            vxc2 = x43*ax/rs*zm
! correlation part
            vxc1 = vxc1 + vcca + fz*(vccp-vcca)                         &
     &           + dfzdz*(eccp-ecca)*( 1.d0-zeta)
            vxc2 = vxc2 + vcca + fz*(vccp-vcca)                         &
     &           + dfzdz*(eccp-ecca)*(-1.d0-zeta)
            exc  = exc + (exc1 + ecca+fz*(eccp-ecca))*roe/2.
            rhor(ir,isup)=vxc1/2.
            rhor(ir,isdw)=vxc2/2.
 20         continue
         end do
      end if
#ifdef PARA
      call reduce(1,exc)
#endif
      return
      end

!-----------------------------------------------------------------------
      subroutine fwfft(f,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!-----------------------------------------------------------------------
! forward fourier transform of potentials and charge density 
! on the dense grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(14,0)
#ifdef PARA
      call cfftp(f,nr1,nr2,nr3,nr1x,nr2x,nr3x,-1)
#else
      call cfft3(f,nr1,nr2,nr3,nr1x,nr2x,nr3x,-1)
#endif
      call tictac(14,1)
      return
      end
!-----------------------------------------------------------------------
      subroutine fwfftb(f,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!-----------------------------------------------------------------------
! forward fourier transform of Q functions (Vanderbilt pseudopotentials)
! on the  box grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(16,0)
      call cfft3b(f,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,-1)
      call tictac(16,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine fwffts(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!-----------------------------------------------------------------------
! forward fourier transform of potentials and charge density 
! on the smooth grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(15,0)
#ifdef PARA
      call cfftps(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,-1)
#else
      call cfft3s(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,-1)
#endif
      call tictac(15,1)
      return
      end
!-----------------------------------------------------------------------
      subroutine fwfftw(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!-----------------------------------------------------------------------
! forward fourier transform of potentials and charge density 
! on the smooth grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(15,0)
#ifdef PARA
      call cfftps(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,-2)
#else
      call cfft3s(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,-1)
#endif
      call tictac(15,1)
      return
      end
!-----------------------------------------------------------------------
      subroutine gausin(eigr,cm)
!-----------------------------------------------------------------------
!
! initialize wavefunctions with gaussians - edit to fit your system
!
      use ions
      use elct, only: n, ngw
      use gvec
!
      implicit none
!
      complex(kind=8) eigr(ngw,nas,nsp), cm(ngw,n)
      real(kind=8)    sigma, auxf
      integer nband, is, ia, ig
!
      sigma=12.0
      nband=0
!!!      do is=1,nsp
      is=1
         do ia=1,na(is)
! s-like gaussians
            nband=nband+1
            do ig=1,ngw
               auxf=exp(-g(ig)/sigma**2)
               cm(ig,nband)=auxf*eigr(ig,ia,is)
            end do
! px-like gaussians
            nband=nband+1
            do ig=1,ngw
               auxf=exp(-g(ig)/sigma**2)
               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,1)
            end do
! py-like gaussians
            nband=nband+1
            do ig=1,ngw
               auxf=exp(-g(ig)/sigma**2)
               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,2)
            end do
! pz-like gaussians
            nband=nband+1
            do ig=1,ngw
               auxf=exp(-g(ig)/sigma**2)
               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,3)
            end do
         end do
      is=2
         do ia=1,na(is)
! s-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)
!            end do
! px-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,1)
!            end do
! py-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,2)
!            end do
! pz-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,3)
!            end do
! dxy-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,1)*gx(ig,2)
!            end do
! dxz-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,1)*gx(ig,3)
!            end do
! dxy-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*gx(ig,2)*gx(ig,3)
!            end do
! dx2-y2-like gaussians
!            nband=nband+1
!            do ig=1,ngw
!               auxf=exp(-g(ig)/sigma**2)
!               cm(ig,nband)=auxf*eigr(ig,ia,is)*                        &
!     &              (gx(ig,1)**2-gx(ig,2)**2)
!            end do
         end do
!!!      end do
      return
      end
!            
!-------------------------------------------------------------------------
      subroutine ggen                                                   &
     &     ( b1, b2, b3, nr1, nr2, nr3, nr1s, nr2s, nr3s,               &
     &      nr1x, nr2x, nr3x, nr1sx, nr2sx, nr3sx, gcut, gcuts, gcutw )
!-----------------------------------------------------------------------
!   generates the reciprocal lattice vectors (g>) with length squared 
!   less than gcut and returns them in order of increasing length.
!      g=i*b1+j*b2+k*b3,
!   where b1,b2,b3 are the vectors defining the reciprocal lattice
!
!   Only half of the g vectors (g>) are stored:
!   if g is present, -g is not (with the exception of g=0)
!   The set g> is defined by
!          g> = line(i=j=0,k>0)+plane(i=0,j>0)+space(i>0)
!
!   n1p,n2p, and n3p are the fast-fourier transform indexes of g> :
!      n1p=i+1      if i.ge.0
!      n1p=i+1+nr1  if i.lt.0
!   and the similar definitions for n2p and n3p.
!
!   n1m,n2m, and n3m are the fft indexes for g<, that is, the set
!   of vectors g=-i*b1-j*b2-k*b3 . These can be shown to be:
!      n1m=1          if i.eq.0 (or n1p.eq.1)
!      n1m=nr1-n1p+2  if i.ne.0
!   and the similar definitions for n2m and n3m.
!
!   the indexes (n1p,n2p,n3p) are collapsed into a one-dimensional
!   index np, and the same applies to negative vectors indexes
!
!   The fft indices are displaced by one unit so that g=0 corresponds
!   to element (1,1,1) (and not (0,0,0)) for purely historical reasons.
!   Negative coefficients are refolded to positive coefficients,
!   introducing a factor  exp(m*2pi*i)=1 in the fourier transform.
!
!   For a transform length n and for a single axis, if n odd:
!   -n-1                       n-1 n+1
!   ----, ..., -1, 0, 1, ...., ---,---,....,n-1  is the "true" index i,
!     2         |  |  |         2   2
!     |         |  |  |         |
!     |         |  V  V         V
!     |         |              n+1 n+3
!     |         |  1, 2, ...., ---,---,....,n    is the fft index n1 of G
!     |         |               2   2  
!     | folding  \_________________ | ______|
!     |_____________________________|
!
! so: if (n1.le.(n+1)/2) i=n1-1 , otherwise, i=n1-n-1
! 
! If n is even:
!     n                         n  n
!    -- , ..., -1, 0, 1, ....,  - ,-+1,....,n-1  is the "real" index i,
!     2         |  |  |         2  2
!     |         |  |  |         |
!     |         |  V  V         V
!     |         |              n   n 
!     |         |  1, 2, ...., -+1,-+2,....,n    is the fft index n1 of G
!     |         |              2   2  
!     | folding  \_____________ | __________|
!     |_________________________|   
!
! so: if (n1.le.n/2+1) i=n1-1 ; if(n1.gt.n/2+1) i=n1-n-1 ;
!     if (n1.eq.n/2+1) i=n1-1 or i=n1-n-1, depending on how
!     the G vectors are refolded
!
!   The indices in1p, in2p, in3p are the i,j,k values.
!   They are used to quickly calculate the structure factors 
!      eigt=exp(-i*g*tau)        (i=imaginary unit!)
!   by decomposing eigt into products of exponentials:
!      eigt=ei1(i)*ei2(j)*ei3(k) (i=index, see above!).
!
!   ng is the total number of vectors with length squared less than gcut.   
!
!   The smooth grid of g with length squared less than gcuts 
!   (gcuts.le.gcut) is calculated in this routine.
!   Smooth grid variables have an "s" appended.
!
!   ngw is the total number of vectors with length squared less than gcutw
!   (gcutw.le.gcut).
!
!   the g's are in units of 2pi/a.
!
      use gvec
      use gvecs
      use elct
#ifdef PARA
      use para_mod
#endif     
      implicit none
      integer nr1,nr2,nr3, nr1s,nr2s,nr3s
      integer nr1x,nr2x,nr3x, nr1sx,nr2sx,nr3sx
      real(kind=8) b1(3), b2(3), b3(3), gcut, gcuts, gcutw
      integer n1p, n2p, n3p, n1m, n2m, n3m
      integer n1ps, n2ps, n3ps, n1ms, n2ms, n3ms
!       cray:
!      integer jwork(257)
      integer, allocatable:: index(:)
      integer it, icurr, nr1m1, nr2m1, nr3m1, nrefold, ir, ig, i,j,k
      real(kind=8) t(3), g2
!
      if(gcut.lt.gcuts) call error('ggen','gcut .lt. gcuts',1)
      ng=0
      ngs=0
      ngw=0
!
! NOTA BENE: these limits are larger than those actually needed
! (-nr/2,..,+nr/2  for nr even; -(nr-1)/2,..,+(nr-1)/2  for nr odd).
! This allows to use a slightly undersized fft grid, with some degree 
! of G-vector refolding, at your own risk
!
      nr1m1=nr1-1
      nr2m1=nr2-1
      nr3m1=nr3-1
!
!     first step : count the number of vectors with g2 < gcut
!
!     exclude space with x<0
!
      do i= 0,nr1m1
         do j=-nr2m1,nr2m1
!
!     exclude plane with x=0, y<0
!
            if(i.eq.0.and.j.lt.0) go to 10
!
            do k=-nr3m1,nr3m1
!
!     exclude line with x=0, y=0, z<0
!
               if(i.eq.0.and.j.eq.0.and.k.lt.0) go to 20
#ifdef PARA
!
!     consider only columns that belong to this node
!
               n1p = i + 1
               if (n1p.lt.1) n1p = n1p + nr1
               n2p = j + 1
               if (n2p.lt.1) n2p = n2p + nr2
               if (ipc(n1p+(n2p-1)*nr1x).eq.0) go to 20
#endif
               g2=0.d0
               do ir=1,3
                  t(ir) = dfloat(i)*b1(ir)                              &
     &                   +dfloat(j)*b2(ir)                              &
     &                   +dfloat(k)*b3(ir)
                  g2=g2+t(ir)*t(ir)
               end do
               if(g2.gt.gcut) go to 20
               ng=ng+1
               if(g2.lt.gcutw) ngw=ngw+1
               if(g2.lt.gcuts) ngs=ngs+1
 20            continue
            end do
 10         continue
         end do
      end do
!
!     second step: allocate space
!
      allocate(gx(ng,3))
      allocate(g(ng))
      allocate(np(ng))
      allocate(nm(ng))
      allocate(igl(ng))
      allocate(in1p(ng))
      allocate(in2p(ng))
      allocate(in3p(ng))
      allocate(index(ng))
!
!     third step : find the vectors with g2 < gcut
!
      ng=0
!
!     exclude space with x<0
!
      do i= 0,nr1m1
         do j=-nr2m1,nr2m1
!
!     exclude plane with x=0, y<0
!
            if(i.eq.0.and.j.lt.0) go to 15
!
            do k=-nr3m1,nr3m1
!
!     exclude line with x=0, y=0, z<0
!
               if(i.eq.0.and.j.eq.0.and.k.lt.0) go to 25
#ifdef PARA
!
!     consider only columns that belong to this node
!
               n1p = i + 1
               if (n1p.lt.1) n1p = n1p + nr1
               n2p = j + 1
               if (n2p.lt.1) n2p = n2p + nr2
               if (ipc(n1p+(n2p-1)*nr1x).eq.0) go to 25
#endif
               g2=0.d0
               do ir=1,3
                  t(ir) = dfloat(i)*b1(ir)                              &
     &                   +dfloat(j)*b2(ir)                              &
     &                   +dfloat(k)*b3(ir)
                  g2=g2+t(ir)*t(ir)
               end do
               if(g2.gt.gcut) go to 25
               ng=ng+1
               g(ng)=g2
               in1p(ng)=i
               in2p(ng)=j
               in3p(ng)=k
 25            continue
            end do
 15         continue
         end do
      end do
!
      write(6,*)
      write(6,150) ng
 150  format(' ggen:  # of g vectors < gcut   ng= ',i6)
      write(6,160) ngs
 160  format(' ggen:  # of g vectors < gcuts ngs= ',i6)
      write(6,170) ngw
 170  format(' ggen:  # of g vectors < gcutw ngw= ',i6)
!
!     reorder the g's in order of increasing magnitude.
!
!     cray:
!     call orders (2,jwork,g,index,ng,1,8,1)
!     generic:
      call kb07ad(g,ng,index)  
!
      do ig=1,ng-1
         icurr=ig
 30      if(index(icurr).ne.ig) then
!     comment if not using cray orders from here
!     g2=g(icurr)
!        g(icurr)=g(index(icurr))
!        g(index(icurr))=g2
!     to here.
            it=in1p(icurr)
            in1p(icurr)=in1p(index(icurr))
            in1p(index(icurr))=it
            it=in2p(icurr)
            in2p(icurr)=in2p(index(icurr))
            in2p(index(icurr))=it
            it=in3p(icurr)
            in3p(icurr)=in3p(index(icurr))
            in3p(index(icurr))=it
!
            it=icurr
            icurr=index(icurr)
            index(it)=it
            if(index(icurr).eq.ig) then
               index(icurr)=icurr
               goto 35
            endif
            goto 30
         endif
 35      continue
      end do
!
! check for the presence of refolded G-vectors (dense grid)
!
      nrefold=0
      if (mod(nr1,2).eq.0) then
         nr1m1=nr1/2-1
      else
         nr1m1=(nr1-1)/2
      end if
      if (mod(nr2,2).eq.0) then
         nr2m1=nr2/2-1
      else
         nr2m1=(nr2-1)/2
      end if
      if (mod(nr3,2).eq.0) then
         nr3m1=nr3/2-1
      else
         nr3m1=(nr3-1)/2
      end if
      do ig=1,ng
         if ( in1p(ig).lt.-nr1m1.or.in1p(ig).gt.nr1m1 .or.              &
     &        in2p(ig).lt.-nr2m1.or.in2p(ig).gt.nr2m1 .or.              &
     &        in3p(ig).lt.-nr3m1.or.in3p(ig).gt.nr3m1      )            &
     &        nrefold=nrefold+1
      end do
      if (nrefold.ne.0) write(6, '('' WARNING: '',i6,                   &
     &     '' G-vectors refolded into dense FFT grid'')') nrefold
!
! costruct fft indexes (n1,n2,n3) for the dense grid
!
      do ig=1,ng
         i=in1p(ig)
         j=in2p(ig)
         k=in3p(ig)
!
! n1p,n2p,n3p: indexes of G
! negative indexes are refolded (note that by construction i.ge.0)
!
         n1p=i+1
         n2p=j+1
         n3p=k+1
         if(i.lt.0) n1p=n1p+nr1
         if(j.lt.0) n2p=n2p+nr2
         if(k.lt.0) n3p=n3p+nr3
!
! n1m,n2m,n3m: indexes of -G
!
         if(i.eq.0) then
            n1m=1
         else
            n1m=nr1-n1p+2
         end if
         if(j.eq.0) then
            n2m=1
         else
            n2m=nr2-n2p+2
         end if
         if(k.eq.0) then
            n3m=1
         else
            n3m=nr3-n3p+2
         end if
#ifdef PARA
!
! conversion from (i,j,k) index to combined 1-d ijk index
! for the parallel case: columns along z are stored contiguously
!
         np(ig) = n3p + (ipc(n1p+(n2p-1)*nr1x)-1)*nr3x
         nm(ig) = n3m + (ipc(n1m+(n2m-1)*nr1x)-1)*nr3x
#else
!
! conversion from (i,j,k) index to combined 1-d ijk index:
! ijk = 1 + (i-1)+(j-1)*ix+(k-1)*ix*jx
! where the (i,j,k) array is assumed to be dimensioned (ix,jx,kx)
!
         np(ig) = n1p + (n2p-1)*nr1x + (n3p-1)*nr1x*nr2x
         nm(ig) = n1m + (n2m-1)*nr1x + (n3m-1)*nr1x*nr2x
#endif
      end do
!
! check for the presence of refolded G-vectors (smooth  grid)
!
      nrefold=0
      if (mod(nr1s,2).eq.0) then
         nr1m1=nr1s/2-1
      else
         nr1m1=(nr1s-1)/2
      end if
      if (mod(nr2s,2).eq.0) then
         nr2m1=nr2s/2-1
      else
         nr2m1=(nr2s-1)/2
      end if
      if (mod(nr3s,2).eq.0) then
         nr3m1=nr3s/2-1
      else
         nr3m1=(nr3s-1)/2
      end if
      do ig=1,ngs
         if ( in1p(ig).lt.-nr1m1.or.in1p(ig).gt.nr1m1 .or.              &
     &        in2p(ig).lt.-nr2m1.or.in2p(ig).gt.nr2m1 .or.              &
     &        in3p(ig).lt.-nr3m1.or.in3p(ig).gt.nr3m1      )            &
     &        nrefold=nrefold+1
      end do
      if (nrefold.ne.0) write(6, '('' WARNING: '',i6,                   &
     &     '' G-vectors refolded into smooth FFT grid'')') nrefold
!
! costruct fft indexes (n1s,n2s,n3s) for the small grid
!
      allocate(nps(ngs))
      allocate(nms(ngs))
!
      do ig=1,ngs
         i=in1p(ig)
         j=in2p(ig)
         k=in3p(ig)
!
! n1ps,n2ps,n3ps: indexes of G
! negative indexes are refolded (note that by construction i.ge.0)
!
         n1ps=i+1
         n2ps=j+1
         n3ps=k+1
         if(i.lt.0) n1ps=n1ps+nr1s
         if(j.lt.0) n2ps=n2ps+nr2s
         if(k.lt.0) n3ps=n3ps+nr3s
!
! n1ms,n2ms,n3ms: indexes of -G
!
         if(i.eq.0) then
            n1ms=1
         else
            n1ms=nr1s-n1ps+2
         end if
         if(j.eq.0) then
            n2ms=1
         else
            n2ms=nr2s-n2ps+2
         end if
         if(k.eq.0) then
            n3ms=1
         else
            n3ms=nr3s-n3ps+2
         end if
#ifdef PARA
!
! conversion from (i,j,k) index to combined 1-d ijk index
! for the parallel case: columns along z are stored contiguously
!
         nps(ig) = n3ps + (ipcs(n1ps+(n2ps-1)*nr1sx)-1)*nr3sx
         nms(ig) = n3ms + (ipcs(n1ms+(n2ms-1)*nr1sx)-1)*nr3sx
#else
!
! conversion from (i,j,k) index to combined 1-d ijk index:
!
         nps(ig) = n1ps+(n2ps-1)*nr1sx+(n3ps-1)*nr1sx*nr2sx
         nms(ig) = n1ms+(n2ms-1)*nr1sx+(n3ms-1)*nr1sx*nr2sx
#endif
      end do
! 
! shells of G - first calculate their number and position
!
      ngl=1
      igl(1)=ngl
      do ig=2,ng
         if(abs(g(ig)-g(ig-1)).gt.1.e-6)then
            ngl=ngl+1
            if (g(ig).lt.gcuts) ngsl=ngl
            if (g(ig).lt.gcutw) ngwl=ngl
         endif
         igl(ig)=ngl
      end do
! 
! then allocate the array gl
!
      allocate(gl(ngl))
! 
! and finally fill gl with the values of the shells
!
      gl(igl(1))=g(1)
      do ig=2,ng
         if(igl(ig).ne.igl(ig-1)) gl(igl(ig))=g(ig)
      end do
!
! ng0 is the index of the first nonzero G-vector
! needed in the parallel case (G=0 is found on one node only!)
!
      if (g(1).lt.1.e-6) then
         ng0=2
      else
         ng0=1
      end if
!
      write(6,180) ngl
 180  format(' ggen:  # of g shells  < gcut  ngl= ',i6)
      write(6,*)
!
! calculation of G-vectors
!
      do ig=1,ng
         i=in1p(ig)
         j=in2p(ig)
         k=in3p(ig)
         gx(ig,1)=i*b1(1)+j*b2(1)+k*b3(1)
         gx(ig,2)=i*b1(2)+j*b2(2)+k*b3(2)
         gx(ig,3)=i*b1(3)+j*b2(3)+k*b3(3)
      end do
!
      return
      end
!-----------------------------------------------------------------------
      subroutine ggenb (b1b, b2b, b3b,                                  &
     &                  nr1b ,nr2b, nr3b, nr1bx ,nr2bx, nr3bx, gcutb )
!-----------------------------------------------------------------------
!
! As ggen, for the box grid. A "b" is appended to box variables.
! The documentation for ggen applies
!
      use gvecb
!
      implicit none
!
      integer nr1b, nr2b, nr3b, nr1bx, nr2bx, nr3bx
      real(kind=8) b1b(3), b2b(3), b3b(3), gcutb
!
      integer, allocatable:: index(:)
      integer n1pb, n2pb, n3pb, n1mb, n2mb, n3mb
      integer it, icurr, nr1m1, nr2m1, nr3m1, ir, ig, i,j,k
!       cray:
!      integer jwork(257)
      real(kind=8) t(3), g2
!
      nr1m1=nr1b-1
      nr2m1=nr2b-1
      nr3m1=nr3b-1
      ngb=0
!
!     first step : count the number of vectors with g2 < gcutb
!
!     exclude space with x<0
!
      do i= 0,nr1m1
         do j=-nr2m1,nr2m1
!
!     exclude plane with x=0, y<0
!
            if(i.eq.0.and.j.lt.0) go to 10
!
            do k=-nr3m1,nr3m1
!
!     exclude line with x=0, y=0, z<0
!
               if(i.eq.0.and.j.eq.0.and.k.lt.0) go to 20
               g2=0.d0
               do ir=1,3
                  t(ir) = dfloat(i)*b1b(ir)                             &
     &                  + dfloat(j)*b2b(ir)                             &
     &                  + dfloat(k)*b3b(ir)
                  g2=g2+t(ir)*t(ir)
               end do
               if(g2.gt.gcutb) go to 20
               ngb=ngb+1
 20            continue
            end do
 10         continue
         end do
      end do
!
!     second step: allocate space
!
      allocate(gxb(ngb,3))
      allocate(gxnb(ngb,3))
      allocate(gb(ngb))
      allocate(npb(ngb))
      allocate(nmb(ngb))
      allocate(iglb(ngb))
      allocate(in1pb(ngb))
      allocate(in2pb(ngb))
      allocate(in3pb(ngb))
      allocate(index(ngb))
!
!     third step : find the vectors with g2 < gcutb
!
      ngb=0
!
!     exclude space with x<0
!
      do i= 0,nr1m1
         do j=-nr2m1,nr2m1
!
!     exclude plane with x=0, y<0
!
            if(i.eq.0.and.j.lt.0) go to 15
!
            do k=-nr3m1,nr3m1
!
!     exclude line with x=0, y=0, z<0
!
               if(i.eq.0.and.j.eq.0.and.k.lt.0) go to 25
               g2=0.d0
               do ir=1,3
                  t(ir) = dfloat(i)*b1b(ir)                             &
     &                  + dfloat(j)*b2b(ir)                             &
     &                  + dfloat(k)*b3b(ir)
                  g2=g2+t(ir)*t(ir)
               end do
               if(g2.gt.gcutb) go to 25
               ngb=ngb+1
               gb(ngb)=g2
               in1pb(ngb)=i
               in2pb(ngb)=j
               in3pb(ngb)=k
 25            continue
            end do
 15         continue
         end do
      end do
!
      write(6,*)
      write(6,170) ngb
 170  format(' ggenb: # of gb vectors < gcutb ngb = ',i6)
!
!   reorder the g's in order of increasing magnitude.
!       cray:
!       call orders (2,jwork,gb,index,ngb,1,8,1)
!       generic:
      call kb07ad (gb,ngb,index)  
      do ig=1,ngb-1
         icurr=ig
 30      if(index(icurr).ne.ig) then
!     comment if not using cray orders from here 
!         g2=gb(icurr)
!         gb(icurr)=gb(index(icurr))
!         gb(index(icurr))=g2
!       to here.
            it=in1pb(icurr)
            in1pb(icurr)=in1pb(index(icurr))
            in1pb(index(icurr))=it
            it=in2pb(icurr)
            in2pb(icurr)=in2pb(index(icurr))
            in2pb(index(icurr))=it
            it=in3pb(icurr)
            in3pb(icurr)=in3pb(index(icurr))
            in3pb(index(icurr))=it
!
            it=icurr
            icurr=index(icurr)
            index(it)=it
            if(index(icurr).eq.ig) then
               index(icurr)=icurr
               goto 35
            endif
            goto 30
         endif
 35      continue
      end do
!
      deallocate(index)
!
! costruct fft indexes (n1b,n2b,n3b) for the box grid
!
      do ig=1,ngb
         i=in1pb(ig)
         j=in2pb(ig)
         k=in3pb(ig)
         n1pb=i+1
         n2pb=j+1
         n3pb=k+1
!
! n1pb,n2pb,n3pb: indexes of G
! negative indexes are refolded (note that by construction i.ge.0)
!
         if(i.lt.0) n1pb=n1pb+nr1b
         if(j.lt.0) n2pb=n2pb+nr2b
         if(k.lt.0) n3pb=n3pb+nr3b
!
! n1mb,n2mb,n3mb: indexes of -G
!
         if(i.eq.0) then
            n1mb=1
         else
            n1mb=nr1b-n1pb+2
         end if
         if(j.eq.0) then
            n2mb=1
         else
            n2mb=nr2b-n2pb+2
         end if
         if(k.eq.0) then
            n3mb=1
         else
            n3mb=nr3b-n3pb+2
         end if
!
! conversion from (i,j,k) index to combined 1-d ijk index:
! ijk = 1 + (i-1)+(j-1)*ix+(k-1)*ix*jx
! where the (i,j,k) array is assumed to be dimensioned (ix,jx,kx)
!
         npb(ig) = n1pb+(n2pb-1)*nr1bx+(n3pb-1)*nr1bx*nr2bx
         nmb(ig) = n1mb+(n2mb-1)*nr1bx+(n3mb-1)*nr1bx*nr2bx
      end do
! 
! shells of G - first calculate their number and position
!
      nglb=1
      iglb(1)=nglb
      do ig=2,ngb
         if(abs(gb(ig)-gb(ig-1)).gt.1.e-6)then
            nglb=nglb+1
         endif
         iglb(ig)=nglb
      end do
      write(6,180) nglb
 180  format(' ggenb: # of gb shells  < gcutb nglb= ',i6)
! 
! then allocate the array glb
!
      allocate(glb(nglb))
! 
! and finally fill glb with the values of the shells
!
      glb(iglb(1))=gb(1)
      do ig=2,ngb
         if(iglb(ig).ne.iglb(ig-1)) glb(iglb(ig))=gb(ig)
      end do
!
! calculation of G-vectors
!
      do ig=1,ngb
         i=in1pb(ig)
         j=in2pb(ig)
         k=in3pb(ig)
         gxb(ig,1)=i*b1b(1)+j*b2b(1)+k*b3b(1)
         gxb(ig,2)=i*b1b(2)+j*b2b(2)+k*b3b(2)
         gxb(ig,3)=i*b1b(3)+j*b2b(3)+k*b3b(3)
      end do
!
      do i=1,3
         do ig=2,ngb
            gxnb(ig,i)=gxb(ig,i)/sqrt(gb(ig))
         end do
      end do
!
      return
      end
!-------------------------------------------------------------------------
      subroutine gracsc(bec,betae,cp,i,csc)
!-----------------------------------------------------------------------
!     requires in input the updated bec(k) for k<i
!     on output: bec(i) is recalculated
!
      use ions
      use cvan
      use elct
!
      implicit none
!
      integer i
      complex(kind=8) betae(ngw,nhsa)
      real(kind=8)  bec(nhsa,n), cp(2,ngw,n)
      real(kind=8)  csc(nx)
      integer k, kmax,ig, is, iv, jv, ia, inl, jnl
      real(kind=8) sum, SSUM
      real(kind=8) temp(ngw) ! automatic array
!
!     calculate csc(k)=<cp(i)|cp(k)>,  k<i
!
      kmax=i-1
      do k=1,kmax
         csc(k)=0.
         if (ispin(i).eq.ispin(k)) then
            do ig=1,ngw
               temp(ig)=cp(1,ig,k)*cp(1,ig,i)+cp(2,ig,k)*cp(2,ig,i)
            end do
            csc(k)=2.*SSUM(ngw,temp,1)
            if (ng0.eq.2) csc(k)=csc(k)-temp(1)
         endif
      end do
#ifdef PARA
      call reduce(kmax,csc)
#endif
!
!     calculate bec(i)=<cp(i)|beta>
!
      do inl=1,nhsavb
         do ig=1,ngw
            temp(ig)=cp(1,ig,i)* real(betae(ig,inl))+             &
     &               cp(2,ig,i)*aimag(betae(ig,inl))
         end do
         bec(inl,i)=2.*SSUM(ngw,temp,1)
         if (ng0.eq.2) bec(inl,i)= bec(inl,i)-temp(1)
      end do
#ifdef PARA
      call reduce(nhsavb,bec(1,i))
#endif
!
!     calculate csc(k)=<cp(i)|S|cp(k)>,  k<i
!
      do k=1,kmax
         if (ispin(i).eq.ispin(k)) then
            sum=0.
            do is=1,nvb
               do iv=1,nh(is)
                  do jv=1,nh(is)
                     if(abs(qq(iv,jv,is)).gt.1.e-5) then 
                        do ia=1,na(is)
                           inl=ish(is)+(iv-1)*na(is)+ia
                           jnl=ish(is)+(jv-1)*na(is)+ia
                           sum = sum + qq(iv,jv,is)*bec(inl,i)*bec(jnl,k)
                        end do
                     endif
                  end do
               end do
            end do
            csc(k)=csc(k)+sum
         endif
      end do
!
!     orthogonalized cp(i) : |cp(i)>=|cp(i)>-\sum_k<i csc(k)|cp(k)>
!
!     corresponing bec:  bec(i)=<cp(i)|beta>-csc(k)<cp(k)|beta>
!
      do k=1,kmax
         do inl=1,nhsavb
            bec(inl,i)=bec(inl,i)-csc(k)*bec(inl,k)
         end do
      end do
!
      return
      end
!-------------------------------------------------------------------------
      subroutine graham(betae,bec,cp)
!-----------------------------------------------------------------------
!     gram-schmidt orthogonalization of the set of wavefunctions cp
!
      use cvan
      use elct
!
      implicit none
!
      real(kind=8)  bec(nhsa,n)
      complex(kind=8)   cp(ngw,n), betae(ngw,nhsa)
!
      real(kind=8) anorm, cscnorm
      real(kind=8) csc(nx) ! automatic array
      integer i,k,is,iv,ia,inl
      external cscnorm
!
      do i=1,n
         call gracsc(bec,betae,cp,i,csc)
!
! calculate orthogonalized cp(i) : |cp(i)>=|cp(i)>-\sum_k<i csc(k)|cp(k)>
!
         do k=1,i-1
            call AXPY(2*ngw,-csc(k),cp(1,k),1,cp(1,i),1)
         end do
         anorm =cscnorm(bec,cp,i)
         call SCAL(2*ngw,1.0/anorm,cp(1,i),1)
!
!         these are the final bec's
!
         call SCAL(nhsavb,1.0/anorm,bec(1,i),1)
      end do
!
      return
      end
!
!-----------------------------------------------------------------------
      subroutine herman_skillman_grid(mesh,z,cmesh,r)
!-----------------------------------------------------------------------
!
      implicit none
!
      integer mesh
      real(kind=8) z, cmesh, r(mesh)
!
      real(kind=8) deltax
      integer nblock,i,j,k
!
      nblock = mesh/40
      i=1
      r(i)=0.0
      cmesh=0.88534138/z**(1.0/3.0)
      deltax=0.0025*cmesh
      do j=1,nblock
         do k=1,40
            i=i+1
            r(i)=r(i-1)+deltax
         end do
         deltax=deltax+deltax
      end do
!
      return
      end
!
!-----------------------------------------------------------------------
      subroutine herman_skillman_int(mesh,cmesh,func,asum)
!-----------------------------------------------------------------------
!     simpsons rule integration for herman skillman mesh
!     mesh - # of mesh points
!     c    - 0.8853418/z**(1/3.)
!
      implicit none
      integer mesh
      real(kind=8) cmesh, func(mesh), asum
!
      integer i, j, k, i1, nblock
      real(kind=8) a1, a2e, a2o, a2es, h
!
      a1=0.0
      a2e=0.0
      asum=0.0
      h=0.0025*cmesh
      nblock=mesh/40
      i=1
      func(1)=0.0
      do j=1,nblock
         do k=1,20
            i=i+2
            i1=i-1
            a2es=a2e
            a2o=func(i1)/12.0
            a2e=func(i)/12.0
            a1=a1+5.0*a2es+8.0*a2o-a2e
            func(i1)=asum+a1*h
            a1=a1-a2es+8.0*a2o+5.0*a2e
            func(i)=asum+a1*h
         end do
         asum=func(i)
         a1=0.0
         h=h+h
      end do
!
      return
      end
!
!-----------------------------------------------------------------------
      subroutine initbox ( nr1, nr2, nr3, nr1b, nr2b, nr3b, nsp, na,    &
     &                     a1,a2,a3, ainv, tau0, taub, irb )
!-----------------------------------------------------------------------
!
!     sets the indexes irb and positions taub for the small boxes 
!     around atoms
!
      use ion_parameters
      implicit none
! input
      integer nr1, nr2, nr3, nr1b, nr2b, nr3b, nsp, na(nsp)
      real(kind=8) a1(3), a2(3), a3(3), ainv(3,3), tau0(3,nax,nsx)
! output
      integer irb(3,nax,nsx)
      real(kind=8) taub(3,nax,nsx)
! local
      real(kind=8) x(3), xmod
      integer nr(3), nrb(3), xint, is, ia, i
!
      call tictac(2,0)
      nr (1)=nr1
      nr (2)=nr2
      nr (3)=nr3
      nrb(1)=nr1b
      nrb(2)=nr2b
      nrb(3)=nr3b
!
      do is=1,nsp
         do ia=1,na(is)
!
            do i=1,3
!
! bring atomic positions in crystal axis
!
               x(i) = ainv(i,1)*tau0(1,ia,is) +                         &
     &                ainv(i,2)*tau0(2,ia,is) +                         &
     &                ainv(i,3)*tau0(3,ia,is)
!
! bring x(i) in the range between 0 and 1
!
               x(i) = mod(x(i),1.d0)
               if (x(i).lt.0.d0) x(i)=x(i)+1.d0
! nr(i) even
               if (mod(nr(i),2).eq.0) then
!
! irb is the index of the grid point at the corner of the small box
! Note that the indices of the small box run from irb to irb+nrb-1
!
                  xint=int(x(i)*nr(i))
                  irb (i,ia,is)=xint+1-nrb(i)/2+1
                  if(irb(i,ia,is).lt.1) irb(i,ia,is)=irb(i,ia,is)+nr(i)
!
! x(i) are the atomic positions in crystal coordinates, where the
! "crystal lattice" is the small box lattice and the origin is at
! the corner of the small box. Used to calculate phases exp(iG*taub)
!
                  xmod=x(i)*nr(i)-xint
                  x(i)=(xmod+nrb(i)/2-1)/nr(i)
               else
!
! nr(i) odd - see above for comments
!
                  xint=nint(x(i)*nr(i))
                  irb (i,ia,is)=xint+1-(nrb(i)+1)/2
                  if(irb(i,ia,is).lt.1) irb(i,ia,is)=irb(i,ia,is)+nr(i)
                  xmod=x(i)*nr(i)-xint
                  x(i)=(xmod+(nrb(i)-1)/2)/nr(i)
               end if
            end do
!
! bring back taub in cartesian coordinates
!
            do i=1,3
               taub(i,ia,is)= x(1)*a1(i) + x(2)*a2(i) + x(3)*a3(i)
            end do
         end do
      end do
      call tictac(2,1)
!
      return
      end
!
!-----------------------------------------------------------------------
      subroutine invfft(f,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!-----------------------------------------------------------------------
! inverse fourier transform of potentials and charge density 
! on the dense grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(14,0)
#ifdef PARA
      call cfftp(f,nr1,nr2,nr3,nr1x,nr2x,nr3x,1)
#else
      call cfft3(f,nr1,nr2,nr3,nr1x,nr2x,nr3x,1)
#endif
      call tictac(14,1)
!
      return
      end 
!-----------------------------------------------------------------------
      subroutine ivfftb(f,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!-----------------------------------------------------------------------
! inverse fourier transform of Q functions (Vanderbilt pseudopotentials)
! on the  box grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(16,0)
      call cfft3b(f,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx,1)
      call tictac(16,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine ivffts(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!-----------------------------------------------------------------------
! inverse fourier transform of  potentials and charge density
! on the smooth grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(15,0)
#ifdef PARA
      call cfftps(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,1)
#else
      call cfft3s(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,1)
#endif
      call tictac(15,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine ivfftw(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx)
!-----------------------------------------------------------------------
! inverse fourier transform of wavefunctions 
! on the smooth grid . On output, f is overwritten
!
      complex(kind=8) f(*)
      call tictac(15,0)
#ifdef PARA
      call cfftps(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,2)
#else
      call cfft3s(f,nr1s,nr2s,nr3s,nr1sx,nr2sx,nr3sx,1)
#endif
      call tictac(15,1)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine kb07ad(count,n,index)
!-------------------------------------------------------------------------
!
!             kb07ad      handles double precision variables
!  standard fortran 66 (a verified pfort subroutine)
!  the work-space 'mark' of length 50 permits up to 2**(50/2) numbers
!  to be sorted. this is more than the ibm virtual memory space
!  will hold .
      implicit real(kind=8)(a-h,o-z)
      integer n, index(n)
      real(kind=8) count(n)
      integer mark(50)
!  set index array to original order .
      do i=1,n
         index(i)=i
      end do
!  check that a trivial case has not been entered .
      if(n.eq.1) go to 10
      if(n.gt.1) go to 30
      write(6,20)
   20 format(///20x,'***kb07ad***no numbers to be sorted ** return to', &
     & ' calling program' )
      goto 10
!  'm' is the length of segment which is short enough to enter
!  the final sorting routine. it may be easily changed.
   30 m=12
!  set up initial values.
      la=2
      is=1
      if=n
      do 190 mloop=1,n
!  if segment is short enough sort with final sorting routine .
      ifka=if-is
      if((ifka+1).gt.m)goto 70
!********* final sorting ***
!  ( a simple bubble sort )
      is1=is+1
      do 60 j=is1,if
      i=j
   40 if(count(i-1).lt.count(i))goto 60
      if(count(i-1).gt.count(i))goto 50
      if(index(i-1).lt.index(i))goto 60
   50 av=count(i-1)
      count(i-1)=count(i)
      count(i)=av
      int=index(i-1)
      index(i-1)=index(i)
      index(i)=int
      i=i-1
      if(i.gt.is)goto  40
   60 continue
      la=la-2
      goto 170
!             *******  quicksort  ********
!  select the number in the central position in the segment as
!  the test number.replace it with the number from the segment's
!  highest address.
   70 iy=(is+if)/2
      x=count(iy)
      intest=index(iy)
      count(iy)=count(if)
      index(iy)=index(if)
!  the markers 'i' and 'ifk' are used for the beginning and end
!  of the section not so far tested against the present value
!  of x .
      k=1
      ifk=if
!  we alternate between the outer loop that increases i and the
!  inner loop that reduces ifk, moving numbers and indices as
!  necessary, until they meet .
      do 110 i=is,if
      if(x.gt.count(i))goto 110
      if(x.lt.count(i))goto 80
      if(intest.gt.index(i))goto 110
   80 if(i.ge.ifk)goto 120
      count(ifk)=count(i)
      index(ifk)=index(i)
      k1=k
      do 100 k=k1,ifka
      ifk=if-k
      if(count(ifk).gt.x)goto 100
      if(count(ifk).lt.x)goto 90
      if(intest.le.index(ifk))goto 100
   90 if(i.ge.ifk)goto 130
      count(i)=count(ifk)
      index(i)=index(ifk)
      go to 110
  100 continue
      goto 120
  110 continue
!  return the test number to the position marked by the marker
!  which did not move last. it divides the initial segment into
!  2 parts. any element in the first part is less than or equal
!  to any element in the second part, and they may now be sorted
!  independently .
  120 count(ifk)=x
      index(ifk)=intest
      ip=ifk
      goto 140
  130 count(i)=x
      index(i)=intest
      ip=i
!  store the longer subdivision in workspace.
  140 if((ip-is).gt.(if-ip))goto 150
      mark(la)=if
      mark(la-1)=ip+1
      if=ip-1
      goto 160
  150 mark(la)=ip-1
      mark(la-1)=is
      is=ip+1
!  find the length of the shorter subdivision.
  160 lngth=if-is
      if(lngth.le.0)goto 180
!  if it contains more than one element supply it with workspace .
      la=la+2
      goto 190
  170 if(la.le.0)goto 10
!  obtain the address of the shortest segment awaiting quicksort
  180 if=mark(la)
      is=mark(la-1)
  190 continue
   10 return
      end
!-------------------------------------------------------------------------
      subroutine latgen(ibrav,celldm,a1,a2,a3,omega)
!-----------------------------------------------------------------------
!     sets up the crystallographic vectors a1,a2, and a3.
!
!     ibrav is the structure index:
!       1  cubic p (sc)                8  orthorhombic p
!       2  cubic f (fcc)               9  one face centered orthorhombic
!       3  cubic i (bcc)              10  all face centered orthorhombic
!       4  hexagonal and trigonal p   11  body centered orthorhombic
!       5  trigonal r                 12  monoclinic p
!       6  tetragonal p (st)          13  one face centered monoclinic
!       7  tetragonal i (bct)         14  triclinic p
!
!     celldm are parameters which fix the shape of the unit cell
!
!     NOTA BENE: the cases ibrav=5, 7, 10 have been modified so as to yield
!     right-handed axis triplets. Boxes for US PPs do not seem to work with
!     the original (left-handed) axis triplets. PG 27-Jan-2001
!
      implicit none
      integer ibrav
      real(kind=8) celldm(6),a1(3),a2(3),a3(3)
!
      integer i,j,k,l,iperm,ir
      real(kind=8) sr2,sr3,term,cbya,s,term1,omega,term2,singam,sen
!
      sr2=sqrt(2.d0)
      sr3=sqrt(3.d0)
!
      do ir=1,3
         a1(ir)=0.d0
         a2(ir)=0.d0
         a3(ir)=0.d0
      end do
!
      go to (2,4,6,8,10,12,14,16,18,20,22,24,26,28),ibrav
      call error('latgen',' bravais lattice not programmed',ibrav)
!
    2 a1(1)=celldm(1)
      a2(2)=celldm(1)
      a3(3)=celldm(1)
      go to 100
!
    4 term=celldm(1)/2.d0
      a1(1)=-term
      a1(3)=term
      a2(2)=term
      a2(3)=term
      a3(1)=-term
      a3(2)=term
      go to 100
!
    6 term=celldm(1)/2.d0
      do ir=1,3
         a1(ir)=term
         a2(ir)=term
         a3(ir)=term
      end do
      a2(1)=-term
      a3(1)=-term
      a3(2)=-term
      go to 100
!
    8 cbya=celldm(3)
      a1(1)=celldm(1)
      a2(1)=-celldm(1)/2.d0
      a2(2)=celldm(1)*sr3/2.d0
      a3(3)=celldm(1)*cbya
      go to 100
!
   10 term1=sqrt(1.d0+2.d0*celldm(4))
      term2=sqrt(1.d0-celldm(4))
      a2(2)=sr2*celldm(1)*term2/sr3
      a2(3)=celldm(1)*term1/sr3
      a1(1)=celldm(1)*term2/sr2
      a1(2)=-a1(1)/sr3
      a1(3)= a2(3)
      a3(1)=-a1(1)
      a3(2)= a1(2)
      a3(3)= a2(3)
      go to 100
!
   12 cbya=celldm(3)
      a1(1)=celldm(1)
      a2(2)=celldm(1)
      a3(3)=celldm(1)*cbya
      go to 100
!
   14 cbya=celldm(3)
      a2(1)=celldm(1)/2.d0
      a2(2)=a2(1)
      a2(3)=cbya*celldm(1)/2.d0
      a1(1)= a2(1)
      a1(2)=-a2(1)
      a1(3)= a2(3)
      a3(1)=-a2(1)
      a3(2)=-a2(1)
      a3(3)= a2(3)
      go to 100
!
   16 a1(1)=celldm(1)
      a2(2)=celldm(1)*celldm(2)
      a3(3)=celldm(1)*celldm(3)
      go to 100
!
   18 a1(1) = 0.5 * celldm(1)
      a1(2) = a1(1) * celldm(2)
      a2(1) = - a1(1)
      a2(2) = a1(2)
      a3(3) = celldm(1) * celldm(3)
      go to 100
!
   20 a2(1) = 0.5 * celldm(1)
      a2(2) = a2(1) * celldm(2)
      a1(1) = a2(1)
      a1(3) = a2(1) * celldm(3)
      a3(2) = a2(1) * celldm(2)
      a3(3) = a1(3)
      go to 100
!
   22 a1(1) = 0.5 * celldm(1)
      a1(2) = a1(1) * celldm(2)
      a1(3) = a1(1) * celldm(3)
      a2(1) = - a1(1)
      a2(2) = a1(2)
      a2(3) = a1(3)
      a3(1) = - a1(1)
      a3(2) = - a1(2)
      a3(3) = a1(3)
      go to 100
!
   24 sen=sqrt(1.d0-celldm(4)**2)
      a1(1)=celldm(1)
      a2(1)=celldm(1)*celldm(2)*celldm(4)
      a2(2)=celldm(1)*celldm(2)*sen
      a3(3)=celldm(1)*celldm(3)
      go to 100
!
   26 sen = sqrt( 1.d0 - celldm(4) ** 2 )
      a1(1) = celldm(1) * celldm(4)
      a1(3) = celldm(1) * sen
      a2(1) = a1(1)
      a2(3) = - a1(3)
      a3(1) = celldm(1) * celldm(2)
      a3(2) = celldm(1) * celldm(3)
      go to 100
!
   28 singam=sqrt(1.d0-celldm(6)**2)
      term=sqrt((1.d0+2.d0*celldm(4)*celldm(5)*celldm(6)               &
     & -celldm(4)**2-celldm(5)**2-celldm(6)**2)/(1.d0-celldm(6)**2))
      a1(1)=celldm(1)
      a2(1)=celldm(1)*celldm(2)*celldm(6)
      a2(2)=celldm(1)*celldm(2)*singam
      a3(1)=celldm(1)*celldm(3)*celldm(5)
      a3(2)=celldm(1)*celldm(3)*(celldm(4)-celldm(5)*celldm(6))/singam
      a3(3)=celldm(1)*celldm(3)*term
!
  100 omega=0.d0
      s=1.d0
      i=1
      j=2
      k=3
!
 101  do iperm=1,3
         omega=omega+s*a1(i)*a2(j)*a3(k)
         l=i
         i=j
         j=k
         k=l
      end do
!
      i=2
      j=1
      k=3
      s=-s
      if(s.lt.0.d0) go to 101
      omega=abs(omega)
      return
!
      end
! CORRECTED BY A.BONGIORNO
!  - F90 COMPILER ON SX5 DOESN'T INTERPRET CORRECTLY SOME
!    STRUCTURES!
!  - IN ORDER TO USE THE HIGHEST VECT-OPT LEVEL THE FOLLOWING
!    ROUTINE HAS BEEN SLIGHTLY MODIFIED
!-------------------------------------------------------------------------
      subroutine newd(vr,tfor,irb,eigrb,rhovan,deeq,fion)
!-----------------------------------------------------------------------
!     this routine calculates array deeq
!
!      deeq_i,ij = omega( v(g=0) q_i,ij(g=0) +
!              2 sum_g> re[ v*(g) q_i,ij(g) e^-ig.r_i ] )
!
!     routine makes use of q(-g) = q*(g)  and  v(-g) = v*(g)
!
      use cvan, only:nh, nhx, nvb
      use ions
      use cnst
      use parm
      use gvecb
      use parmb
      use qgb_mod
      use elct
      use control
      use work_box
#ifdef PARA
      use para_mod
#endif
!
      implicit none
!
      logical tfor
      integer irb(3,nax,nsx)
      real(kind=8)  vr(nnr,nspin), rhovan(nat,nhx,nhx,nspin)
      real(kind=8)  deeq(nat,nhx,nhx,nspin), fion(3,nax,nsp)
      complex(kind=8) eigrb(ngb,nas,nsp)
!
      integer isup, isdw, ir1,ir2,ir3, ir,ibig1,ibig2,ibig3, ib,        &
     &     iv, ijv, iac, ik, jv, ism, ia, nc, isa, iss, is, ig, ic
      real(kind=8)    qsum, c
      complex(kind=8) fp, fm, ci, sum, CSUM
      complex(kind=8), allocatable:: vtotg(:,:), vqg(:)
      real(kind=8), parameter :: Z_0 = 0.d0
!
      call tictac(6,0)
      ci=(0.,1.)
      qsum=0.
!
      allocate(vtotg(ngb,2))
      allocate(vqg(ngb))
!
      if(nspin.eq.1)then
!     =================================================================
!     case nspin=1
!     -----------------------------------------------------------------
         iss=1
         do is=1,nvb
            isa=0
            do ism=1,is-1
               isa=isa+na(ism)
            end do
            do ia=1,na(is),2
               nc=2
               if( ia.eq.na(is).and.na(is)/2.ne.real(na(is))/2.) nc=1
               call zero(2*nnrb,qv)
               do ir3=1,nr3b
                  ibig3=irb(3,ia,is)+ir3-1
                  ibig3=1+mod(ibig3-1,nr3)
#ifdef PARA
                  ibig3=ibig3-n3(me)
                  if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                  do ir2=1,nr2b
                     ibig2=irb(2,ia,is)+ir2-1
                     ibig2=1+mod(ibig2-1,nr2)
                     do ir1=1,nr1b
                        ibig1=irb(1,ia,is)+ir1-1
                        ibig1=1+mod(ibig1-1,nr1)
                        ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                        ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                        qv(ir) = cmplx(vr(ib,iss),Z_0)
                     end do
                  end do
#ifdef PARA
                  endif
#endif
               end do
!
               if(nc.eq.2)then
                  do ir3=1,nr3b
                     ibig3=irb(3,ia+1,is)+ir3-1
                     ibig3=1+mod(ibig3-1,nr3)
#ifdef PARA
                     ibig3=ibig3-n3(me)
                     if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                     do ir2=1,nr2b
                        ibig2=irb(2,ia+1,is)+ir2-1
                        ibig2=1+mod(ibig2-1,nr2)
                        do ir1=1,nr1b
                           ibig1=irb(1,ia+1,is)+ir1-1
                           ibig1=1+mod(ibig1-1,nr1)
                           ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                           ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                           qv(ir)= qv(ir) + cmplx(Z_0,vr(ib,iss))
                        end do
                     end do
#ifdef PARA
                     endif
#endif
                  end do
               endif
#ifdef PARA
!
! this is a simple (but not efficient nor elegant) way to
! ensure that every node has its complete copy of the small box
!
               call reduce(2*nnrb,qv)
#endif
               call fwfftb(qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!
               do ig=1,ngb
                  fp=  qv(npb(ig)) + qv(nmb(ig))
                  fm=  qv(npb(ig)) - qv(nmb(ig))
                  vtotg(ig,1)=cmplx(0.5*real(fp) , 0.5*aimag(fm))
                  vtotg(ig,2)=cmplx(0.5*aimag(fp),-0.5*real(fm))
               end do
!
               do ic=1,nc
                  iac=ia+ic-1
                  isa=isa+1
                  ijv=0
                  do iv= 1,nh(is)
                     do jv=iv,nh(is)
                        ijv=ijv+1
                        if(tfor) then
                           if(iv.ne.jv) then
                              c=2.*omegab*tpibab*rhovan(isa,iv,jv,iss)
                           else
                              c=   omegab*tpibab*rhovan(isa,iv,jv,iss)        
                           endif
                           do ik=1,3
                              do ig=1,ngb
                                 vqg(ig)=conjg(vtotg(ig,ic))            &
     &                                *qgb(ig,ijv,is)*eigrb(ig,iac,is)  &
     &                                *cmplx(0.,-gxb(ig,ik))
                              end do
                              sum=2.0*CSUM(ngb-1,vqg(2),1)+vqg(1)
                              fion(ik,iac,is)=fion(ik,iac,is)           &
     &                             -c*real(sum)
                           end do
                        endif
!
                        do ig=1,ngb
                           vqg(ig)=qgb(ig,ijv,is)*conjg(vtotg(ig,ic))   &
     &                          *eigrb(ig,iac,is)
                        end do
                        sum=2.0*CSUM(ngb-1,vqg(2),1)+vqg(1)
                        deeq(isa,iv,jv,iss)=omegab*real(sum)
                        deeq(isa,jv,iv,iss)=deeq(isa,iv,jv,iss)
                        if(iv.ne.jv) then
                           qsum=qsum+2.*deeq(isa,iv,jv,iss)
                        else
                           qsum=qsum+   deeq(isa,iv,jv,iss)
                        endif
                     end do
                  end do
               end do
            end do
         end do
!
      else
!     =================================================================
!     case nspin=2
!     -----------------------------------------------------------------
         isup=1
         isdw=2
! isup=spin up, isdw=spin down
         do is=1,nvb
            isa=0
            do ism=1,is-1
               isa=isa+na(ism)
            end do
            do ia=1,na(is)
               call zero(2*nnrb,qv)
               do ir3=1,nr3b
                  ibig3=irb(3,ia,is)+ir3-1
                  ibig3=1+mod(ibig3-1,nr3)
#ifdef PARA
                  ibig3=ibig3-n3(me)
                  if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                  do ir2=1,nr2b
                     ibig2=irb(2,ia,is)+ir2-1
                     ibig2=1+mod(ibig2-1,nr2)
                     do ir1=1,nr1b
                        ibig1=irb(1,ia,is)+ir1-1
                        ibig1=1+mod(ibig1-1,nr1)
                        ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                        ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                        qv(ir)= cmplx(vr(ib,isup),vr(ib,isdw))
                     end do
                  end do
#ifdef PARA
                  end if
#endif
               end do
#ifdef PARA
!
! this is a simple (but not efficient nor elegant) way to
! ensure that every node has its complete copy of the small box
!
               call reduce(2*nnrb,qv)
#endif
!     
               call fwfftb(qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!
               do ig=1,ngb
                  fp=  qv(npb(ig)) + qv(nmb(ig))
                  fm=  qv(npb(ig)) - qv(nmb(ig))
                  vtotg(ig,1)=cmplx(0.5*real(fp) , 0.5*aimag(fm))
                  vtotg(ig,2)=cmplx(0.5*aimag(fp),-0.5*real(fm))
               end do
!
               isa=isa+1
               do iss=1,nspin
                  ijv=0
                  do iv= 1,nh(is)
                     do jv=iv,nh(is)
                        ijv=ijv+1
!
                        if(tfor) then
                           if(iv.ne.jv) then
                              c=2.*omegab*tpibab*rhovan(isa,iv,jv,iss)
                           else
                              c=   omegab*tpibab*rhovan(isa,iv,jv,iss)        
                           endif
                           do ik=1,3
                              do ig=1,ngb
                                 vqg(ig) = conjg(vtotg(ig,iss))         &
     &                             *qgb(ig,ijv,is)*eigrb(ig,ia,is)      &
     &                                *cmplx(0.,-gxb(ig,ik))
                              end do
                              sum=2.0*CSUM(ngb-1,vqg(2),1)+vqg(1)
                              fion(ik,ia,is)=fion(ik,ia,is)-c*real(sum)
                           end do
                        endif
!
                        do ig=1,ngb
                           vqg(ig)=qgb(ig,ijv,is)*                      &
     &                          conjg(vtotg(ig,iss))*eigrb(ig,ia,is)
                        end do
                        sum=2.0*CSUM(ngb-1,vqg(2),1)+vqg(1)
                        deeq(isa,iv,jv,iss)=omegab*real(sum)
                        deeq(isa,jv,iv,iss)=deeq(isa,iv,jv,iss)
                        if(iv.ne.jv) then
                           qsum=qsum+2.*deeq(isa,iv,jv,iss)
                        else
                           qsum=qsum+   deeq(isa,iv,jv,iss)
                        endif
                     end do
                  end do
               end do
            end do
         end do
      endif
!     -----------------------------------------------------------------
      if(iprsta.gt.2) then
         write(6,'(a,f12.6)') ' newd:  sum_i,ij deeq(i,ij) ',qsum
         write(6,*)
      endif
      if(iprsta.gt.4) then
         do iss=1,nspin
            isa=0
            do is=1,nvb
               do ia=1,na(is)
                  isa=isa+1
                  write(6,*)
                  write(6,'(33x,a,i4)') '    deeq(isa) iss ',isa,iss
                  do iv=1,nh(is)
                     write(6,'(8f9.4)')                                 &
     &                    (deeq(isa,iv,jv,iss),jv=1,nh(is))
                  end do
               end do
            end do
         end do
      endif
!
      deallocate(vqg)
      deallocate(vtotg)
      call tictac(6,1)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine nlfl(bec,becdr,lambda,fion)
!-----------------------------------------------------------------------
!     contribution to fion due to the orthonormality constraint
! 
!
      use ions
      use gvec
      use cvan
      use elct
      use cnst
      use parm
!
      implicit none
      real(kind=8) bec(nhsa,n), becdr(nhsa,n,3), lambda(nx,nx)
      real(kind=8) fion(3,nax,nsp)
!
      integer k, is, ia, iv, jv, i, j, inl
      real(kind=8) tt, SSUM
      real(kind=8) temp(nx,nx), tmpbec(nhx,nx),tmpdr(nx,nhx) ! automatic arrays
!
      do k=1,3
         do is=1,nvb
            do ia=1,na(is)
!
               call zero(nhx*n,tmpbec)
               call zero(nhx*n,tmpdr)
!
               do iv=1,nh(is)
                  do jv=1,nh(is)
                     inl=ish(is)+(jv-1)*na(is)+ia
                     if(abs(qq(iv,jv,is)).gt.1.e-5) then
                        do i=1,n
                           tmpbec(iv,i)=tmpbec(iv,i)                    &
     &                          + qq(iv,jv,is)*bec(inl,i)
                        end do
                     endif
                  end do
               end do
!
               do iv=1,nh(is)
                  inl=ish(is)+(iv-1)*na(is)+ia
                  do i=1,n
                     tmpdr(i,iv)=becdr(inl,i,k)
                  end do
               end do
!
               if(nh(is).gt.0)then
                  call zero(nx*n,temp)
!
                  call MXMA                                             &
     &                 (tmpdr,1,nx,tmpbec,1,nhx,temp,1,nx,n,nh(is),n)
!
                  do j=1,n
                     do i=1,n
                        temp(i,j)=temp(i,j)*lambda(i,j)
                     end do
                  end do
!
                  tt=SSUM(nx*n,temp,1)
!
                  fion(k,ia,is)=fion(k,ia,is)+2.*tt
               endif
!
            end do
         end do
      end do
!
!     end of x/y/z loop
!
      return
      end
!-----------------------------------------------------------------------
      subroutine nlfq(c,deeq,eigr,bec,becdr,fion)
!-----------------------------------------------------------------------
!     contribution to fion due to nonlocal part
!
      use gvec
      use cvan
      use ions
      use elct
      use cnst
      use parm
! 
      implicit none
      real(kind=8) deeq(nat,nhx,nhx,nspin), bec(nhsa,n), becdr(nhsa,n,3),&
     &       c(2,ngw,n)
      complex(kind=8) eigr(ngw,nas,nsp)
      real(kind=8) fion(3,nax,nsx)
!
      integer k, is, ia, isa, iss, inl, iv, jv, i, j
      real(kind=8) tmpbec(nhx,n), tmpdr(nhx,n) ! automatic arrays
      real(kind=8) temp, tt, SSUM
!
!     nlsm2 fills becdr
!
      call tictac(8,0)
      call nlsm2(eigr,c,becdr)
!
      do k=1,3
!
         isa=0
         do is=1,nsp
            do ia=1,na(is)
               isa=isa+1
!
               call zero(nhx*n,tmpbec)
               call zero(nhx*n,tmpdr )
!
               do iv=1,nh(is)
                  do jv=1,nh(is)
                     inl=ish(is)+(jv-1)*na(is)+ia
                     do i=1,n
                        iss=ispin(i)
                        temp=dvan(iv,jv,is)+deeq(isa,jv,iv,iss)
                        tmpbec(iv,i)=tmpbec(iv,i)+temp*bec(inl,i)
                     end do
                  end do
               end do
!  
               do iv=1,nh(is)
                  inl=ish(is)+(iv-1)*na(is)+ia
                  do i=1,n
                     tmpdr(iv,i)=f(i)*becdr(inl,i,k)
                  end do
               end do
!
               do i=1,n
                  do iv=1,nh(is)
                     tmpdr(iv,i)=tmpdr(iv,i)*tmpbec(iv,i)
                  end do
               end do
!
               tt=SSUM(nhx*n,tmpdr,1)
!     
               fion(k,ia,is)=fion(k,ia,is)-2.*tt
!
            end do
         end do
      end do
!
!     end of x/y/z loop
!
      call tictac(8,1)
!
      return
      end
! CORRECTED BY A.BONGIORNO
!  - F90 COMPILER ON SX5 DOESN'T INTERPRET CORRECTLY SOME
!    STRUCTURES!
!  - IN ORDER TO USE THE HIGHEST VECT-OPT LEVEL THE FOLLOWING
!    ROUTINE HAS BEEN SLIGHTLY MODIFIED
!-----------------------------------------------------------------------
      subroutine nlsm1 (n,nspmn,nspmx,eigr,c,becp)
!-----------------------------------------------------------------------
!     computes: the array becp
!     becp(ia,n,iv,is)=
!         = sum_g [(-i)**l beta(g,iv,is) e^(-ig.r_ia)]^* c(g,n)
!         = delta_l0 beta(g=0,iv,is) c(g=0,n)
!          +sum_g> beta(g,iv,is) 2 re[(i)**l e^(ig.r_ia) c(g,n)]
!
!     routine makes use of c*(g)=c(-g)  (g> see routine ggen)
!     input : beta(ig,l,is), eigr, c
!     output: becp as parameter
!
      use ions, only: na, nas
      use elct, only: ngw, ng0
      use cnst
      use cvan
      use work2
!
      implicit none
      integer n, nspmn, nspmx
      real(kind=8)  eigr(2,ngw,nas,nspmx)
      complex(kind=8) c(ngw,n)
      real(kind=8)  becp(nhsa,n)
!
      integer ig, is, iv, ia, l, ixr(3), ixi(3), inl, ism, i
      real(kind=8) signre(3), signim(3), argre, argim
!
      l=0
      ixr(l+1)=1
      ixi(l+1)=2
      signre(l+1) = 1.d0
      signim(l+1) = 1.d0
!
      l=1
      ixr(l+1)=2
      ixi(l+1)=1
      signre(l+1) = 1.d0
      signim(l+1) = -1.d0
!
      l=2
      ixr(l+1)=1
      ixi(l+1)=2
      signre(l+1) = -1.d0
      signim(l+1) = -1.d0
!
      do is=nspmn,nspmx
         do iv=1,nh(is)
            l=nhtol(iv,is)
!
            do ia=1,na(is)
               if (ng0.eq.2) then
!                   q = 0   component (with weight 1.0)
                    argre = signre(l+1)*beta(1,iv,is)*eigr(ixr(l+1),1,ia,is)
                    argim = signim(l+1)*beta(1,iv,is)*eigr(ixi(l+1),1,ia,is)
                    wrk2(1,ia)= cmplx( argre , argim ) 
               end if
               do ig=ng0,ngw
!                   q > 0   components (with weight 2.0)
                    argre = 2.0*beta(ig,iv,is)*signre(l+1)*eigr(ixr(l+1),ig,ia,is)
                    argim = 2.0*beta(ig,iv,is)*signim(l+1)*eigr(ixi(l+1),ig,ia,is)
                    wrk2(ig,ia) = cmplx( argre , argim )
               end do
            end do
            inl=ish(is)+(iv-1)*na(is)+1
            call MXMA(wrk2,2*ngw,1,c,1,2*ngw,becp(inl,1),1,nhsa,       &
     &           na(is),2*ngw,n)
         end do
#ifdef PARA
         inl=ish(is)+1
         do i=1,n
            call reduce(na(is)*nh(is),becp(inl,i))
         end do
#endif
      end do
      return
      end
! CORRECTED BY A.BONGIORNO
!  - F90 COMPILER ON SX5 DOESN'T INTERPRET CORRECTLY SOME
!    STRUCTURES!
!  - IN ORDER TO USE THE HIGHEST VECT-OPT LEVEL THE FOLLOWING
!    ROUTINE HAS BEEN SLIGHTLY MODIFIED
!-------------------------------------------------------------------------
      subroutine nlsm2(eigr,c,becdr)
!-----------------------------------------------------------------------
!     computes: the array becdr 
!     becdr(ia,n,iv,is,k)
!      =2.0 sum_g> g_k beta(g,iv,is) re[ (i)**(l+1) e^(ig.r_ia) c(g,n)]
! 
!     routine makes use of  c*(g)=c(-g)  (g> see routine ggen)
!     input : eigr, c
!     output: becdr
!
      use ions
      use elct
      use gvec
      use cnst
      use cvan
      use work2
!
      implicit none
!
      real(kind=8)  eigr(2,ngw,nas,nsp),c(2,ngw,n), becdr(nhsa,n,3)
      integer ig, is, iv, ia, k, l, inl
      integer ixr(3), ixi(3)
      real(kind=8) signre(3), signim(3), argre, argim
      real(kind=8), allocatable:: gk(:)
!
      call tictac(22,0)
      allocate(gk(ngw))
      call zero(3*nhsa*n,becdr)
!
      l=0
      ixr(l+1)=2
      ixi(l+1)=1
      signre(l+1) = 1.d0
      signim(l+1) = -1.d0
!
      l=1
      ixr(l+1)=1
      ixi(l+1)=2
      signre(l+1) = -1.d0
      signim(l+1) = -1.d0
!
      l=2
      ixr(l+1)=2
      ixi(l+1)=1
      signre(l+1) = -1.d0
      signim(l+1) = 1.d0
!
      do k=1,3
         do ig=1,ngw
            gk(ig)=gx(ig,k)*tpiba
         end do
!
         do is=1,nsp
            do iv=1,nh(is)
!
!     order of states:  s_1  p_x1  p_z1  p_y1  s_2  p_x2  p_z2  p_y2
!
               l=nhtol(iv,is)
               do ia=1,na(is)
                  if (ng0.eq.2) then
!                             q = 0   component (with weight 1.0)
        argre = signre(l+1)*gk(1)*beta(1,iv,is)*eigr(ixr(l+1),1,ia,is)
        argim = signim(l+1)*gk(1)*beta(1,iv,is)*eigr(ixi(l+1),1,ia,is)
        wrk2(1,ia) = cmplx ( argre , argim )
                  end if
                  do ig=ng0,ngw
!                            q > 0   components (with weight 2.0)
        argre = 2.0*gk(ig)*beta(ig,iv,is)*signre(l+1)*eigr(ixr(l+1),ig,ia,is)
        argim = 2.0*gk(ig)*beta(ig,iv,is)*signim(l+1)*eigr(ixi(l+1),ig,ia,is)
        wrk2(ig,ia) = cmplx ( argre , argim )
                  end do
               end do
               inl=ish(is)+(iv-1)*na(is)+1
               call MXMA(wrk2,2*ngw,1,c,1,2*ngw,becdr(inl,1,k),1,       &
     &                   nhsa,na(is),2*ngw,n)
            end do
         end do
      end do
#ifdef PARA
      call reduce(3*nhsa*n,becdr)
#endif
      deallocate(gk)
      call tictac(22,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine ortho                                                  &
     &      (eigr,cp,phi,x0,diff,iter,ccc,eps,max,delt,bephi,becp)
!-----------------------------------------------------------------------
!     input = cp (non-orthonormal), beta
!     input = phi |phi>=s'|c0>
!     output= cp (orthonormal with s( r(t+dt) ) )
!     output= bephi, becp
!     the method used is similar to the version in les houches 1988
!     'simple molecular systems at..'  p. 462-463  (18-22)
!      xcx + b x + b^t x^t + a = 1
!     where c = <s'c0|s|s'c0>   b = <s'c0|s cp>   a = <cp|s|cp>
!     where s=s(r(t+dt)) and s'=s(r(t))  
!     for vanderbilt pseudo pot - kl & ap
!
      use ions
      use cvan
      use elct
      use control
!
      implicit none
!
      complex(kind=8)   cp(ngw,n), phi(ngw,n), eigr(ngw,nas,nsp)
      real(kind=8) x0(nx,nx), diff, ccc, eps, delt
      integer iter, max
      real(kind=8) bephi(nhsa,n), becp(nhsa,n)
!
      real(kind=8) diag(nx),work1(nx),work2(nx),                        &
     &     xloc(nx,nx),tmp1(nx,nx),tmp2(nx,nx),dd(nx,nx),               &
     &     x1(nx,nx),rhos(nx,nx),rhor(nx,nx),con(nx,nx),  u(nx,nx)
      real(kind=8) sig(nx,nx), rho(nx,nx), tau(nx,nx)
! the above are all automatic arrays
      integer istart, nss, ifail, i, j, iss, iv, jv, ia, is, inl, jnl
      real(kind=8), allocatable:: qbephi(:,:), qbecp(:,:)
!
!     calculation of becp and bephi
!
      call tictac(17,0)
      call nlsm1(n,1,nvb,eigr, cp, becp)
      call nlsm1(n,1,nvb,eigr,phi,bephi)
      call tictac(17,1)
!
      call tictac(24,0)
!
!     calculation of qbephi and qbecp
!
      allocate(qbephi(nhsa,n))
      allocate(qbecp (nhsa,n))
      call zero(nhsa*n,qbephi)
      call zero(nhsa*n,qbecp )
!
      do is=1,nvb
         do iv=1,nh(is)
            do jv=1,nh(is)
               if(abs(qq(iv,jv,is)).gt.1.e-5) then
                  do ia=1,na(is)
                     inl=ish(is)+(iv-1)*na(is)+ia
                     jnl=ish(is)+(jv-1)*na(is)+ia
                     do i=1,n
                        qbephi(inl,i)= qbephi(inl,i)                    &
     &                       +qq(iv,jv,is)*bephi(jnl,i)
                        qbecp (inl,i)=qbecp (inl,i)                     &
     &                       +qq(iv,jv,is)*becp (jnl,i)
                     end do
                  end do
               endif
            end do
         end do
      end do
!
      call tictac(24,1)
!
      do iss=1,nspin
         nss=nupdwn(iss)
         istart=iupdwn(iss)
!
!     rho = <s'c0|s|cp>
!     sig = 1-<cp|s|cp>
!     tau = <s'c0|s|s'c0>
!
         call tictac(19,0)
         call rhoset(cp,phi,bephi,qbecp,nss,istart,rho)
         call sigset(cp,becp,qbecp,nss,istart,sig)
         call tauset(phi,bephi,qbephi,nss,istart,tau)
         call tictac(19,1)
!
         if(iprsta.gt.4) then
            write(6,*)
            write(6,'(26x,a)') '    rho '
            do i=1,nss
               write(6,'(7f11.6)') (rho(i,j),j=1,nss)
            end do
            write(6,*)
            write(6,'(26x,a)') '    sig '
            do i=1,nss
               write(6,'(7f11.6)') (sig(i,j),j=1,nss)
            end do
            write(6,*)
            write(6,'(26x,a)') '    tau '
            do i=1,nss
               write(6,'(7f11.6)') (tau(i,j),j=1,nss)
            end do
         endif
!
!
!----------------------------------------------------------------by ap--
! 
         do j=1,nss
            do i=1,nss
               xloc(i,j) = x0(istart-1+i,istart-1+j)*ccc
               dd(i,j) = 0.d0
               x1(i,j) = 0.d0
               tmp1(i,j)=0.d0
               rhos(i,j)=0.5d0*( rho(i,j)+rho(j,i) )
!
! on some machines (IBM RS/6000 for instance) the following test allows
! to distinguish between Numbers and Sodium Nitride (NaN, Not a Number).
! If a matrix of Not-Numbers is passed to rsg, the most likely outcome is
! that the program goes on forever doing nothing and writing nothing.
!
               if (rhos(i,j).ne.rhos(i,j))                                &
     &             call error('ortho','ortho went bananas',1)
               rhor(i,j)=rho(i,j)-rhos(i,j)
            end do
         end do
!     
         do i=1,nss
            tmp1(i,i)=1.d0
         end do
         ifail=0

         call tictac(20,0)
         call rsg(nx,nss,rhos,tmp1,diag,1,u,work1,work2,ifail)
         call tictac(20,1)
! 
!
!                calculation of lagranges multipliers
!
!
         call tictac(21,0)
         do iter=1,max
!
!       the following 4 MXMA-calls do the following matrix 
!       multiplications:
!                       tmp1 = x0*rhor    (1st call)
!                       dd   = x0*tau*x0  (2nd and 3rd call)
!                       tmp2 = x0*rhos    (4th call)
!
            call MXMA( xloc,1,nx,rhor,1,nx,tmp1,1,nx,nss,nss,nss)
            call MXMA( tau ,1,nx,xloc,1,nx,tmp2,1,nx,nss,nss,nss)
            call MXMA( xloc,1,nx,tmp2,1,nx,  dd,1,nx,nss,nss,nss)
            call MXMA( xloc,1,nx,rhos,1,nx,tmp2,1,nx,nss,nss,nss)
            do i=1,nss
               do j=1,nss
                  x1(i,j) = sig(i,j)-tmp1(i,j)-tmp1(j,i)-dd(i,j)       
                  con(i,j)= x1(i,j)-tmp2(i,j)-tmp2(j,i)
               end do
            end do
!
!         x1      = sig      -x0*rho    -x0*rho^t  -x0*tau*x0
!
            diff=0.d0
            do i=1,nss
               do j=1,nss
                  if(abs(con(i,j)).gt.diff) diff=abs(con(i,j))
               end do
            end do
!
            if( diff.le.eps ) go to 20
!     
!     the following two MXMA-calls do:   
!                       tmp1 = x1*u
!                       tmp2 = ut*x1*u
!
            call MXMA(x1,1,nx,   u,1,nx,tmp1,1,nx,nss,nss,nss)
            call MXMA(u ,nx,1,tmp1,1,nx,tmp2,1,nx,nss,nss,nss)   
!
!       g=ut*x1*u/d  (g is stored in tmp1)
! 
            do i=1,nss
               do j=1,nss
                  tmp1(i,j)=tmp2(i,j)/(diag(i)+diag(j))
               end do
            end do
!     
!       the following two MXMA-calls do:   
!                       tmp2 = g*ut
!                       x0 = u*g*ut
!
            call MXMA(tmp1,1,nx,  u,nx,1,tmp2,1,nx,nss,nss,nss)
            call MXMA(   u,1,nx,tmp2,1,nx,xloc,1,nx,nss,nss,nss)
         end do
         write(6,*) ' diff= ',diff,' iter= ',iter
         call error('ortho','max number of iterations exceeded',iter)
!
 20      continue
!
!-----------------------------------------------------------------------
!
         call tictac(21,1)
!
         if(iprsta.gt.4) then
            write(6,*)
            write(6,'(26x,a)') '    lambda '
            do i=1,nss
               write(6,'(7f11.6)') (xloc(i,j)/f(i+istart-1),j=1,nss)
            end do
         endif
!     
         if(iprsta.gt.2) then
            write(6,*) ' diff= ',diff,' iter= ',iter
         endif
!     
!     lagrange multipliers
!
         do i=1,nss
            do j=1,nss
               x0(istart-1+i,istart-1+j)=xloc(i,j)/ccc
               if (xloc(i,j).ne.xloc(i,j))                                &
     &             call error('ortho','ortho went bananas',2)
            end do
         end do
!
      end do
!
      deallocate(qbecp )
      deallocate(qbephi)

      return
      end
!
!-----------------------------------------------------------------------
      subroutine pbc(rin,a1,a2,a3,ainv,rout)
!-----------------------------------------------------------------------
!
!     brings atoms inside the unit cell
!
      implicit none
! input
      real(kind=8) rin(3), a1(3),a2(3),a3(3), ainv(3,3)
! output
      real(kind=8) rout(3)
! local
      real(kind=8) x,y,z
!
! bring atomic positions to crystal axis
!
      x = ainv(1,1)*rin(1)+ainv(1,2)*rin(2)+ainv(1,3)*rin(3)
      y = ainv(2,1)*rin(1)+ainv(2,2)*rin(2)+ainv(2,3)*rin(3)
      z = ainv(3,1)*rin(1)+ainv(3,2)*rin(2)+ainv(3,3)*rin(3)
!
! bring x,y,z in the range between -0.5 and 0.5
!
      x = x - nint(x)
      y = y - nint(y)
      z = z - nint(z)
!
! bring atomic positions back in cartesian axis
!
      rout(1) = x*a1(1)+y*a2(1)+z*a3(1)
      rout(2) = x*a1(2)+y*a2(2)+z*a3(2)
      rout(3) = x*a1(3)+y*a2(3)+z*a3(3)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine phbox(taub,eigrb)
!-----------------------------------------------------------------------
!     calculates the phase factors for the g's of the little box
!     eigrt=exp(-i*g*tau) .
!     Uses the same logic for fast calculation as in phfac (see below)
!
      use ions
      use gvecb
      use parmb
      use cnst
      use control
!
      implicit none
      real(kind=8) taub(3,nax,nsx)
      complex(kind=8) eigrb(ngb,nas,nsp)
! local
      integer i,j,k, is, ia, ig
      real(kind=8) taup(3),sum,ar1,ar2,ar3
      complex(kind=8), allocatable:: ei1b(:,:,:), ei2b(:,:,:), ei3b(:,:,:)
      complex(kind=8) ctep1,ctep2,ctep3,ctem1,ctem2,ctem3
!
      call tictac(2,0)
!
      if(nr1b.lt.3) call error(' phbox ',' nr1b too small ',nr1b)
      if(nr2b.lt.3) call error(' phbox ',' nr2b too small ',nr2b)
      if(nr3b.lt.3) call error(' phbox ',' nr3b too small ',nr3b)
!
      allocate(ei1b(-nr1b:nr1b,nas,nsp))
      allocate(ei2b(-nr2b:nr2b,nas,nsp))
      allocate(ei3b(-nr3b:nr3b,nas,nsp))
!
      if(iprsta.gt.3) then 
         write(6,*) ' phbox: taub '
         write(6,*)  (((taub(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp)
      endif
      do is=1,nsp
         do ia=1,na(is)
            do i=1,3
               sum=0.d0
               do j=1,3
                  sum=sum+ainvb(i,j)*taub(j,ia,is)
               end do
               taup(i)=sum
            end do
!
            ar1=2.d0*pi*taup(1)
            ctep1=cmplx(cos(ar1),-sin(ar1))
            ctem1=conjg(ctep1)
            ei1b( 0,ia,is)=cmplx(1.d0,0.d0)
            ei1b( 1,ia,is)=ctep1
            ei1b(-1,ia,is)=ctem1
            do i=2,nr1b-1
               ei1b( i,ia,is)=ei1b( i-1,ia,is)*ctep1
               ei1b(-i,ia,is)=ei1b(-i+1,ia,is)*ctem1
            end do
!
            ar2=2.d0*pi*taup(2)
            ctep2=cmplx(cos(ar2),-sin(ar2))
            ctem2=conjg(ctep2)
            ei2b( 0,ia,is)=cmplx(1.d0,0.d0)
            ei2b( 1,ia,is)=ctep2
            ei2b(-1,ia,is)=ctem2
            do j=2,nr2b-1
               ei2b( j,ia,is)=ei2b( j-1,ia,is)*ctep2
               ei2b(-j,ia,is)=ei2b(-j+1,ia,is)*ctem2
            end do
!
            ar3=2.d0*pi*taup(3)
            ctep3=cmplx(cos(ar3),-sin(ar3))
            ctem3=conjg(ctep3)
            ei3b( 0,ia,is)=cmplx(1.d0,0.d0)
            ei3b( 1,ia,is)=ctep3
            ei3b(-1,ia,is)=ctem3
            do k=2,nr3b-1
               ei3b( k,ia,is)=ei3b( k-1,ia,is)*ctep3
               ei3b(-k,ia,is)=ei3b(-k+1,ia,is)*ctem3
            end do
!
         end do
      end do
!
!     calculation of eigrb(g,ia,is)=e^(-ig.r(ia,is))
!
      do is=1,nsp
         do ia=1,na(is)
            do ig=1,ngb
               eigrb(ig,ia,is) = ei1b(in1pb(ig),ia,is) *                 &
     &                           ei2b(in2pb(ig),ia,is) *                 &
     &                           ei3b(in3pb(ig),ia,is)
            end do
         end do
      end do
!
      if(iprsta.gt.4) then
         write(6,*)
         if(nsp.gt.1) then
            do is=1,nsp
               write(6,'(33x,a,i4)') ' ei1b, ei2b, ei3b (is)',is
               do ig=1,4
                  write(6,'(6f9.4)')                                    &
     &                 ei1b(ig,1,is),ei2b(ig,1,is),ei3b(ig,1,is)
               end do
               write(6,*)
            end do
         else
            do ia=1,na(1)
               write(6,'(33x,a,i4)') ' ei1b, ei2b, ei3b (ia)',ia
               do ig=1,4
                  write(6,'(6f9.4)')                                    &
     &                 ei1b(ig,ia,1),ei2b(ig,ia,1),ei3b(ig,ia,1)
               end do
               write(6,*)
            end do
         endif
      endif
!
      deallocate(ei3b)
      deallocate(ei2b)
      deallocate(ei1b)
      call tictac(2,1)
!     
      return
      end
!-----------------------------------------------------------------------
      subroutine phfac(tau0,ei1,ei2,ei3,eigr)
!-----------------------------------------------------------------------
!  this subroutine generates the complex matrices ei1, ei2, and ei3
!  used to compute the structure factor and forces on atoms :
!     ei1(n1,ia,is) = exp(-i*n1*b1*tau(ia,is)) -nr1<n1<nr1
!  and similar definitions for ei2 and ei3 ; and :
!     eigr(n,ia,is) = ei1*ei2*ei3 = exp(-i g*tau(ia,is))
!  The value of n1,n2,n3 for a vector g is supplied by arrays in1p,in2p,in3p
!  calculated in ggen .
!
      use ions
      use elct, only:ngw
      use parm
      use cnst
      use gvec
      use control
!
      implicit none
      real(kind=8) tau0(3,nax,nsx)
!
      complex(kind=8) ei1(-nr1:nr1,nas,nsp), ei2(-nr2:nr2,nas,nsp),      &
     &                ei3(-nr3:nr3,nas,nsp), eigr(ngw,nas,nsp)
!
      integer i,j,k, ia, is, ig
      real(kind=8) taup(3), sum, ar1,ar2,ar3
      complex(kind=8) ctep1,ctep2,ctep3,ctem1,ctem2,ctem3
!
      call tictac(3,0)
!
      if(nr1.lt.3) call error(' phfac ',' nr1 too small ',nr1)
      if(nr2.lt.3) call error(' phfac ',' nr1 too small ',nr2)
      if(nr3.lt.3) call error(' phfac ',' nr1 too small ',nr3)
!
      if(iprsta.gt.3) then
         write(6,*) ' phfac: tau0 '
         write(6,*) (((tau0(i,ia,is),i=1,3),ia=1,na(is)),is=1,nsp)
      endif
      do is=1,nsp
         do ia=1,na(is)
            do i=1,3
               sum=0.d0
               do j=1,3
                  sum=sum+ainv(i,j)*tau0(j,ia,is)
               end do
               taup(i)=sum
!
! tau0=x1*a1+x2*a2+x3*a3 => taup(1)=x1=tau0*b1 and so on
!
            end do
!
            ar1=2.d0*pi*taup(1)
            ctep1=cmplx(cos(ar1),-sin(ar1))
            ctem1=conjg(ctep1)
            ei1( 0,ia,is)=cmplx(1.d0,0.d0)
            ei1( 1,ia,is)=ctep1
            ei1(-1,ia,is)=ctem1
            do i=2,nr1-1
               ei1( i,ia,is)=ei1( i-1,ia,is)*ctep1
               ei1(-i,ia,is)=ei1(-i+1,ia,is)*ctem1
            end do
!
            ar2=2.d0*pi*taup(2)
            ctep2=cmplx(cos(ar2),-sin(ar2))
            ctem2=conjg(ctep2)
            ei2( 0,ia,is)=cmplx(1.d0,0.d0)
            ei2( 1,ia,is)=ctep2
            ei2(-1,ia,is)=ctem2
            do j=2,nr2-1
               ei2( j,ia,is)=ei2( j-1,ia,is)*ctep2
               ei2(-j,ia,is)=ei2(-j+1,ia,is)*ctem2
            end do
!
            ar3=2.d0*pi*taup(3)
            ctep3=cmplx(cos(ar3),-sin(ar3))
            ctem3=conjg(ctep3)
            ei3( 0,ia,is)=cmplx(1.d0,0.d0)
            ei3( 1,ia,is)=ctep3
            ei3(-1,ia,is)=ctem3
            do k=2,nr3-1
               ei3( k,ia,is)=ei3( k-1,ia,is)*ctep3
               ei3(-k,ia,is)=ei3(-k+1,ia,is)*ctem3
            end do
!
         end do
      end do
!
      if(iprsta.gt.4) then
         write(6,*)
         if(nsp.gt.1) then
            do is=1,nsp
               write(6,'(33x,a,i4)') ' ei1, ei2, ei3 (is)',is
               do ig=1,4
                  write(6,'(6f9.4)')                                    &
     &                 ei1(ig,1,is),ei3(ig,1,is),ei3(ig,1,is)
               end do
               write(6,*)
            end do
         else
            do ia=1,na(1)
               write(6,'(33x,a,i4)') ' ei1, ei2, ei3 (ia)',ia
               do ig=1,4
                  write(6,'(6f9.4)')                                    &
     &                 ei1(ig,ia,1),ei3(ig,ia,1),ei3(ig,ia,1)
               end do
               write(6,*)
            end do
         endif
      endif
!
!     calculation of eigr(g,ia,is)=e^(-ig.r(ia,is))
!
      do is=1,nsp
         do ia=1,na(is)
            do ig=1,ngw
               eigr(ig,ia,is) = ei1(in1p(ig),ia,is)                     &
     &                        * ei2(in2p(ig),ia,is)                     &
     &                        * ei3(in3p(ig),ia,is)
            end do
         end do
      end do
      call tictac(3,1)
!
      return
      end
!
!-------------------------------------------------------------------------
      subroutine prefor(eigr,betae)
!-----------------------------------------------------------------------
!
!     input :        eigr =  e^-ig.r_i
!     output:        betae_i,i(g) = (-i)**l beta_i,i(g) e^-ig.r_i 
!
      use ions
      use cvan 
      use elct
!
      implicit none
      complex(kind=8) eigr(ngw,nas,nsp)
      complex(kind=8) betae(ngw,nhsa)
!
      integer is, iv, ia, inl, ig
      complex(kind=8) ci
!
      call tictac(26,0)
      do is=1,nsp
         do iv=1,nh(is)
            ci=(0.,-1.)**nhtol(iv,is)
            do ia=1,na(is)
               inl=ish(is)+(iv-1)*na(is)+ia
               do ig=1,ngw
                  betae(ig,inl)=ci*beta(ig,iv,is)*eigr(ig,ia,is)
               end do
            end do
         end do
      end do
      call tictac(26,1)
!
      return
      end
!
!---------------------------------------------------------------------
      subroutine randin(nmin,nmax,ng0,ngw,ampre,c)
!---------------------------------------------------------------------
!
      implicit none
! input
      integer nmin, nmax, ng0, ngw
      real(kind=8) ampre
! output
      complex(kind=8) c(ngw,nmax)
! local
      integer i,j
      real(kind=8) ranf1, randy, ranf2, ampexp
!
      do i=nmin,nmax
         do j=ng0,ngw
            ranf1=.5-randy()
            ranf2=.5-randy()
            ampexp=ampre*exp(-(4.*j)/ngw)
            c(j,i)=c(j,i)+ampexp*cmplx(ranf1,ranf2)
         end do
      end do
!
      return
      end
!
!---------------------------------------------------------------------
      subroutine readpp
!---------------------------------------------------------------------
!
      use dft_mod!ATTENZIONE questa xe' una zonta
      use ncprm
      use cvan
      use ions
!
      implicit none
!
      integer is
!
!     -----------------------------------------------------------------
!     first vanderbilt species , then norm-conserving are read !!
!
!     ipp=-1 ultrasoft vanderbilt pp , AdC format(read from file 14)
!     ipp= 0 ultrasoft vanderbilt pp , old format(read from file 14)
!     ipp= 1 ultrasoft vanderbilt pp (read from file 14)
!     ipp= 2 norm-conserving hsc pp  (read from file 14)
!     ipp= 3 norm-conserving bhs pp  (read from file 18)
!     -----------------------------------------------------------------
!
      rewind 14
      rewind 15
      rewind 18
      nvb=0
!
      do is=1,nsp
!
!     check on ipp value and input order
!
         if (ipp(is).lt.-1.or.ipp(is).gt.3)                            &
     &       call error('readpp', 'wrong pseudo type',ipp(is))
         if(is.ge.2)then
            if(ipp(is).lt.ipp(is-1))then
               call error('readpp', 'first vdb, then nnc',10)
            endif
         endif
!
!     ------------------------------------------------------
!     counting of u-s vanderbilt species 
!     ------------------------------------------------------
         if(ipp(is).le.1) nvb=nvb+1
!
         write(6,*) 'reading ppot for species is = ',is
!
         if(ipp(is).ne.3) then
!     =================================================================
!     ultrasoft and hsc species  
!     =================================================================
            write(6,*) 'from file 14'
            if (ipp(is).eq.-1) then
               call readAdC(is,14)
            else
               call readvan(is,ipp(is),14)
            end if
         else
!     ==================================================================
!     bhs species  
!     ==================================================================
            write(6,*) 'from file 18'
            call readbhs(is,18)
         endif
      end do
      close(14)
      close(15)
      close(18)
!
!      dft=lda!ATTENZIONE sole per sto debugging 
      return
      end
!     
!---------------------------------------------------------------------
      subroutine readbhs( is, iunps )
!---------------------------------------------------------------------
!
      use ncprm
      use bhs
      use dft_mod
      use ions, only: zv
!
      implicit none
!
      integer is, iunps
!
      integer meshp, ir, ib, il, i, j, jj
      character*20 xctype
      real(kind=8), allocatable:: fint(:)
      real(kind=8) rdum, alpha, z, zval, cmeshp, exfact
!
! ifpcor is unfortunately not read from file
!
      ifpcor(is)=0
      read(iunps,*) z,zval,nbeta(is),lloc(is),exfact
      if (zval.ne.zv(is)) then
         write(6,*) ' readpp: z,zv,nbeta,lloc: ',                      &
              &              z,zval,nbeta(is),lloc(is)
         call error('readpp','wrong potential read',15)
      endif
      if (is.eq.1) then
         dft=exfact
!     else if (dft.ne.exfact) then
!        call error('readpp','inconsistent DFT',is)
! CORRECTION INTRODUCED BY A.BONGIORNO
!   - UPPERCASE LETTERS ARE USED TO DISTINGUISH FROM THE
!     ORIGINAL VERSION!
!
      ELSE IF (DFT.NE.EXFACT) THEN
         WRITE(6,*)'DFT CODE INCONSISTENT !!!!!!!'
         WRITE(6,*)'SPECIES NO. :',IS
         WRITE(6,*)'WARNING ... WARNING ... WARNING ...'
         WRITE(6,*)'WARNING ... WARNING ... WARNING ...'
         WRITE(6,*)'WARNING ... WARNING ... WARNING ...'
         WRITE(6,*)'WARNING ... WARNING ... WARNING ...'
      end if
!
      if(lloc(is).eq.2)then 
         lll(1,is)=0
         lll(2,is)=1
      else if(lloc(is).ne.2) then
         call error('readbhs','kb-ization for lloc=2 only',10)
      endif
!     
!     see eqs. (2.21) and (2.22) of bhs, prb 26, 4199 (1982).
!
!     wrc1  =c_core(1)
!     wrc2  =c_core(2)
!     rc1   =alpha_core(1)
!     rc2   =alpha_core(2)
!     al(i) =a(i)          i=1,3
!     bl(i) =a(i+3)        i=1,3
!     rcl(i)=alpha(i)      i=1,3 
!
!     ------------------------------------------------------------------
!     pp parameters are read from file iunps
!     bhs 's coefficients have been turned into lengths
!     ------------------------------------------------------------------
      read(iunps,*) wrc1(is),rc1(is),wrc2(is),rc2(is)  
      rc1(is)=1.0/sqrt(rc1(is))
      rc2(is)=1.0/sqrt(rc2(is))
      do il=1,3
         do ib=1,3
            read(iunps,*) rcl(ib,is,il),al(ib,is,il),bl(ib,is,il)
            rcl(ib,is,il)=1.0/sqrt(rcl(ib,is,il))
         end do
      end do
!
!     ------------------------------------------------------------------
!     wavefunctions are read from file iunps
!     ------------------------------------------------------------------
      do il=1,nbeta(is)
         read(iunps,*) mesh(is),cmesh(is)
!
! kkbeta is for compatibility with Vanderbilt PP
!
         kkbeta(is)=mesh(is)
         do j=1,mesh(is)
            read(iunps,*) jj,r(j,is),betar(j,il,is)
         end do
      end do
!     
!     ------------------------------------------------------------------
!     core charge is read from unit 15
!     ------------------------------------------------------------------
!
      if(ifpcor(is).eq.1) then
         read(15,*) meshp,cmeshp
         if(meshp.ne.mesh(is).or.cmeshp.ne.cmesh(is))then
            call error('readpp','core charge mesh minmatch',is)
         endif
         do ir=1,mesh(is)
            read(15,*) rdum,rscore(ir,is)
         end do
      endif
!
!  rab(i) is the derivative of the radial mesh
!
      cmesh(is)=log(cmesh(is))
      do ir=1,mesh(is)
         rab(ir,is)=r(ir,is)*cmesh(is)
      end do
!
!     ------------------------------------------------------------------
!     local potential 
!     ------------------------------------------------------------------
      lloc(is)=lloc(is)+1
      do ir=1,mesh(is)
         rucore(ir,lloc(is),is)=0.
         do i=1,3
            rucore(ir,lloc(is),is)=rucore(ir,lloc(is),is)               &
     &            +(al(i,is,lloc(is))+bl(i,is,lloc(is))*r(ir,is)**2)    &
     &            *exp(-(r(ir,is)/rcl(i,is,lloc(is)))**2)
         end do
      end do
!
!     ------------------------------------------------------------------
!     nonlocal potentials: kleinman-bylander form 
!     (1) definition of betar   (2) calculation of dion 
!     ------------------------------------------------------------------
      allocate(fint(mesh(is)))
      do il=1,nbeta(is)
         do ir=1,mesh(is)
            rucore(ir,il,is)=0.
            do i=1,3
               rucore(ir,il,is)=rucore(ir,il,is)                        &
     &                    +(al(i,is,il)+bl(i,is,il)*r(ir,is)**2)        &
     &                    *exp(-(r(ir,is)/rcl(i,is,il))**2)
            end do
            rucore(ir,il,is) = rucore(ir,il,is)                         &
     &                       - rucore(ir,lloc(is),is)
            fint(ir)=betar(ir,il,is)**2*rucore(ir,il,is)
            betar(ir,il,is)=rucore(ir,il,is)*betar(ir,il,is)
         end do
         call simpson(mesh(is),fint,rab(1,is),dion(il,il,is))
         dion(il,il,is)=1.0/dion(il,il,is)
      end do
      deallocate(fint)
!     
!     ------------------------------------------------------------------
!     output: pp info 
!     ------------------------------------------------------------------
      write(6,3000) z,zval 
3000  format(2x,'bhs pp for z=',f3.0,2x,'zv=',f3.0)
      call dftname(xctype)
      write(6,'(2x,a20)') xctype
      write(6,3002) lloc(is)-1 
3002  format(2x,'   local angular momentum: l=',i3)
      write(6,3005) nbeta(is)
3005  format(2x,'number of nl ang. mom. nbeta=',i3)
      do il=1,nbeta(is)
         write(6,3010) lll(il,is)
3010     format(2x,'nonlocal angular momentum: l=',i3)
      end do
      write(6,3030) 
3030  format(2x,'pseudopotential parameters:')
      write(6,3035) wrc1(is),1.0/rc1(is)**2
3035  format(2x,'core:',2x,'c1_c=',f7.4,' alpha1_c=',f7.4)
      write(6,3036) wrc2(is),1.0/rc2(is)**2
3036  format(2x,'     ',2x,'c2_c=',f7.4,' alpha2_c=',f7.4)
      write(6,3038)
3038  format(2x,'other table parameters:')
      do il=1,3
         write(6,3040) il-1
3040     format(2x,'l=',i3)
         do i =1,3
            alpha=1.0/rcl(i,is,il)**2
            write(6,3050) i,alpha,i,al(i,is,il),i+3,bl(i,is,il)
         end do
      end do
3050  format(2x,'alpha',i1,'=',f6.2,'  a',i1,'=',f16.7,                 &
     &           '  a',i1,'=',f16.7)
      write(6,*)
!     
      return
      end
!     
!---------------------------------------------------------------------
      subroutine readAdC( is, iunps )
!---------------------------------------------------------------------
!
!     This routine reads Vanderbilt pseudopotentials produced by the
!     code of Andrea Dal Corso. Hard PPs are first generated
!     according to the Rabe Rappe Kaxiras Johannopoulos recipe.
!     Ultrasoft PP's are subsequently generated from the hard PP's.
!
!     Output parameters in module "ncprm"
!     info on DFT level in module "dft"
!
      use ion_parameters
      use ncprm
      use dft_mod
!
      implicit none
!
!    First the arguments passed to the subroutine
!
      integer                                                           &
     &      is,        &! The number of the pseudopotential
     &      iunps       ! the unit with the pseudopotential
!
!    Local variables
!
      integer                                                           &
     &       nb,mb,     &! counters on beta functions
     &       n,         &! counter on mesh points
     &       ir,        &! counters on mesh points
     &       pseudotype,&! the type of pseudopotential
     &       ios,       &! I/O control
     &       nwfs,      &! the number of pseudo wavefunctions
     &       ndum,      &! dummy integer variable
     &       l,         &! counter on angular momentum
     &       ikk         ! the kkbeta for each beta
      real(kind=8)                                                            &
     &       x,         &! auxiliary variable
     &       etotps,    &! total energy of the pseudoatom
     &       rdum        ! dummy real variable
!
      logical                                                           &
     &       rel      ! if true the atomic calculation is relativistic
!
      character*75                                                      &
     &       titleps    ! the title of the pseudo
!
!   the following variables are not actually used or are translated
!   into corresponding variables used in the Vanderbilt format
!   Short Vanderbilt-AdC dictionary: 
!   nwfs => nchi        nchix => nvalx    rho_at  => rsatom
!   lchi => lll         vnl   => rucore   rho_atc => rscore
!
      real(kind=8)  xmin, zmesh, dx,&! mesh parameters
     &        zp(nsx),        &! ionic charges
     &        oc(nvalx,nsx),  &! occupancies
     &        rsatom(mmaxx),  &! charge density of pseudoatom
     &        chi(mmaxx,nvalx) ! atomic wavefunctions (not used)
      integer mfxcx, mfxcc,   &!
     &        mgcx, mgcc,     &! exch-corr functional indices 
     &        exfact,         &
     &        lmin, lmax,     &! min and max l
     &        nchi(nsx)        ! number of pseudo wavefunctions
      logical nlcc             ! core correction flag (same as ifpcor)
!
!
      if (is.lt.0 .or. is.gt.nsx)                                       &
     &   call error('readAdC','Wrong is number', 1)
!
      read( iunps, '(a75)', err=100, iostat=ios ) titleps
!
      read( iunps, '(i5)',err=100, iostat=ios ) pseudotype
      if (pseudotype.eq.3) then
         write(6,'('' RRKJ3 Ultrasoft PP for '',a2)') titleps(7:8)
      else
         write(6,'('' RRKJ3 norm-conserving PP for '',a2)') titleps(7:8)
      endif
      read( iunps, '(2l5)',err=100, iostat=ios ) rel, nlcc
      if (nlcc) then
         ifpcor(is)=1
      else
         ifpcor(is)=0
      end if
      read( iunps, '(4i5)',err=100, iostat=ios ) mfxcx, mfxcc,          &
     &                                            mgcx, mgcc
!
      if (mfxcx.eq.1.and.mfxcc.eq.1.and.mgcx.eq.0.and.mgcc.eq.0) then
         exfact=lda
      else if (mfxcx.eq.1.and.mfxcc.eq.1.and.mgcx.eq.1.and.mgcc.eq.1) then
         exfact=bp88
      else if (mfxcx.eq.1.and.mfxcc.eq.4.and.mgcx.eq.2.and.mgcc.eq.2) then
         exfact=pw91
      else if(mfxcx.eq.1.and.mfxcc.eq.4.and.mgcx.eq.3.and.mgcc.eq.4) then
         exfact=pbe
      else if(mfxcx.eq.1.and.mfxcc.eq.3.and.mgcx.eq.1.and.mgcc.eq.3) then
         exfact=blyp
      else
         call error('readAdC','which functional is this?',1)
      end if
!
      if (is.eq.1) then
         dft=exfact
      else if (dft.ne.exfact) then
         call error('readAdC','inconsistent DFT',is)
      end if
!
      read( iunps, '(2e17.11,i5)') zp(is), etotps, lmax
!
      read( iunps, '(4e17.11,i5)',err=100, iostat=ios )                 &
     &                       xmin,rdum,zmesh,dx,mesh(is)
!
      if (mesh(is).gt.mmaxx.or. mesh(is).lt.0)                          &
     &   call error('readAdC', 'wrong mesh',is)
!
      read( iunps, '(2i5)', err=100, iostat=ios ) nchi(is), nbeta(is)
!
      if (nbeta(is).gt.nbrx.or. nbeta(is).lt.0)                         &
     &   call error('readAdC', 'wrong nbeta', is)
      if (nchi(is).gt.nvalx.or.nchi(is).lt.1)                           &
     &   call error('readAdC', 'wrong nchi', is)
!
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                              ( rdum, nb=1,nchi(is) )
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                              ( rdum, nb=1,nchi(is) )
!
      do nb=1,nchi(is)
         read(iunps,'(a2,2i3,f6.2)',err=100,iostat=ios)                 &
     &                rdum, ndum, lll(nb,is), oc(nb,is)
      enddo
!
      kkbeta(is)=0
      do nb=1,nbeta(is)
         read ( iunps, '(i6)',err=100, iostat=ios ) ikk
         kkbeta(is)=max(kkbeta(is),ikk)
         read ( iunps, '(1p4e19.11)',err=100, iostat=ios )              &
     &                      ( betar(ir,nb,is), ir=1,ikk)
         do ir=ikk+1,mesh(is)
            betar(ir,nb,is)=0.d0
         enddo
         do mb=1,nb
            read( iunps, '(1p4e19.11)', err=100, iostat=ios )           &
     &            dion(nb,mb,is)
            dion(mb,nb,is)=dion(nb,mb,is)
            if (pseudotype.eq.3) then
               read(iunps,'(1p4e19.11)',err=100,iostat=ios)             &
     &              qqq(nb,mb,is)
               qqq(mb,nb,is)=qqq(nb,mb,is)
               read(iunps,'(1p4e19.11)',err=100,iostat=ios)             &
     &              (qfunc(n,nb,mb,is),n=1,mesh(is))
               do n=1,mesh(is)
                  qfunc(n,mb,nb,is)=qfunc(n,nb,mb,is)
               enddo
            else
               qqq(nb,mb,is)=0.d0
               qqq(mb,nb,is)=0.d0
               do n=1,mesh(is)
                  qfunc(n,nb,mb,is)=0.d0
                  qfunc(n,mb,nb,is)=0.d0
               enddo
            endif
         enddo
      enddo
!
!   reads the local potential 
!
      read( iunps, '(1p4e19.11)',err=100, iostat=ios ) rdum,            &
     &                       ( rucore(ir,1,is), ir=1,mesh(is) )
!
!   reads the atomic charge
!
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                         ( rsatom(ir), ir=1,mesh(is) )
!
!   if present reads the core charge
!
      if ( nlcc ) then 
         read( iunps, '(1p4e19.11)', err=100, iostat=ios )              &
     &                         ( rscore(ir,is), ir=1,mesh(is) )
      endif
!
!   read the pseudo wavefunctions of the atom
!  
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &             ((chi(ir,nb),ir=1,mesh(is)),nb=1,nchi(is))
!
!    set several variables for compatibility with the rest of the
!    code
!
      nqlc(is)=2*lmax+1
      if ( nqlc(is).gt.lqx .or. nqlc(is).lt.0 )                         &
     &     call error(' readAdC', 'Wrong  nqlc', nqlc(is) )
      do l=1,nqlc(is)
         rinner(l,is)=0.d0
      enddo
!
!    fill the q(r)
!
      do nb=1,nbeta(is)
         do mb=nb,nbeta(is)
            lmin=abs(lll(mb,is)-lll(nb,is))+1
            lmax=lmin+2*lll(nb,is)
            do l=lmin,lmax,2
               do ir=1,kkbeta(is)
                  qrl(ir,nb,mb,l,is)=qfunc(ir,nb,mb,is)
               end do
            end do
         end do
      end do
!
!    compute the radial mesh
!
      do ir = 1, mesh(is)
        x = xmin + float(ir-1) * dx
        r(ir,is) = exp(x) / zmesh
        rab(ir,is) = dx * r(ir,is)
      end do
!
!    For compatibility with Vandebilt format
!
      do ir = 1, mesh(is)
         rucore(ir,1,is)= rucore(ir,1,is)*r(ir,is)
      end do
!
!    For compatibility with rho_atc in the non-US case
!
!!!      if ( nlcc(is) ) then
!!!         do ir=1,mesh(is)
!!!            rho_atc(ir,is) = rho_atc(ir,is)/fpi/r(ir,is)**2
!!!         enddo
!!!      end if
!
      return
100   call error('readAdC','Reading pseudo file',abs(ios))
      stop
      end
!     
!---------------------------------------------------------------------
      subroutine readvan( is, ipp, iunps )
!---------------------------------------------------------------------
!
!     Read Vanderbilt pseudopotential of type "ipp" for species "is" 
!     from unit "iunps"
!     Output parameters in module "ncprm"
!     info on DFT level in module "dft"
!
!
!     ------------------------------------------------------
!     Important:
!     ------------------------------------------------------
!     The order of all l-dependent objects is always s,p,d
!     ------------------------------------------------------
!     potentials, e.g. rucore, are really r*v(r)
!     wave funcs, e.g. chi, are really proportional to r*psi(r)
!     and are normalized so int (chi**2) dr = 1
!     thus psi(r-vec)=(1/r)*chi(r)*y_lm(theta,phi)
!     conventions carry over to beta, etc
!     charge dens, e.g. rscore, really 4*pi*r**2*rho
!
!     ------------------------------------------------------
!     Notes on qfunc and qfcoef:
!     ------------------------------------------------------
!     Since Q_ij(r) is the product of two orbitals like
!     psi_{l1,m1}^star * psi_{l2,m2}, it can be decomposed by
!     total angular momentum L, where L runs over | l1-l2 | ,
!     | l1-l2 | +2 , ... , l1+l2.  (L=0 is the only component
!     needed by the atomic program, which assumes spherical
!     charge symmetry.)
!
!     Recall  qfunc(r) = y1(r) * y2(r)  where y1 and y2 are the
!     radial parts of the wave functions defined according to
!
!       psi(r-vec) = (1/r) * y(r) * Y_lm(r-hat)  .
!
!     For each total angular momentum L, we pseudize qfunc(r)
!     inside rc as:
!
!       qfunc(r) = r**(L+2) * [ a_1 + a_2*r**2 + a_3*r**4 ]
!
!     in such a way as to match qfunc and its 1'st derivative at
!     rc, and to preserve
!
!       integral dr r**L * qfunc(r)   ,
!
!     i.e., to preserve the L'th moment of the charge.  The array
!     qfunc has been set inside rc to correspond to this pseudized
!     version using the minimal L, namely L = | l1-l2 | (e.g., L=0
!     for diagonal elements).  The coefficients a_i (i=1,2,3)
!     are stored in the array qfcoef(i,L+1,j,k) for each L so that
!     the correctly pseudized versions of qfunc can be reconstructed
!     for each L.  (Note that for given l1 and l2, only the values
!     L = | l1-l2 | , | l1-l2 | +2 , ... , l1+l2 are ever used.)
!     ------------------------------------------------------
!
!
      use ion_parameters
      use ncprm
      use dft_mod
!
      implicit none
!
!    First the arguments passed to the subroutine
!
      integer                                                           &
     &      is,        &! The number of the pseudopotential
     &      ipp,       &! The type of pseudopotential
     &      iunps       ! The unit of the pseudo file
!
!   The local variables which are used to read the Vanderbilt file.
!   They are used only for the pseudopotential report. They are no 
!   longer used in the rest of the code.
!
      real(kind=8)                                                            &
     &       exfact,        &! index of the exchange and correlation used 
     &       etotpseu,      &! total pseudopotential energy
     &       wwnl(nvalx),   &! the occupation of the valence states
     &       ee(nvalx),     &! the energy of the valence states
     &       eloc,          &! energy of the local potential
     &       dummy,         &! dummy real variable
     &       rc(nvalx),     &! the cut-off radii of the pseudopotential
     &       eee(nbrx),     &! energies of the beta function
     &       ddd(nbrx,nbrx),&! the screened D_{\mu,\nu} parameters
     &       rcloc,         &! the cut-off radius of the local potential 
     &       vloc0(mmaxx),  &! the screened local potential
     &       rsatom(mmaxx), &! the charge density of pseudoatom
     &       chi(mmaxx,nvalx),&! atomic wavefunctions (not used)
     &       zp(nsx),       &! ionic charge of the pseudoatom
     &       z(nsx)          ! atomic charge
      integer                                                           &
     &       iver(3),       &! contains the version of the code
     &       idmy(3),       &! contains the date of creation of the pseudo
     &       ios,           &! integer variable for I/O control
     &       i,             &! dummy counter 
     &       nnlz(nvalx),   &! The nlm values of the valence states
     &       keyps,         &! the type of pseudopotential. Only US allowed
     &       irel,          &! it says if the pseudopotential is relativistic
     &       ifqopt(nsx),   &! level of Q optimization
     &       iptype(nbrx),  &! more recent parameters 
     &       npf,           &! as above
     &       nang,          &! number of angular momenta in pseudopotentials
     &       lloc,          &! angular momentum of the local part of PPs
     &       lmin, lmax,    &! min and max angular momentum in Q
     &       lp,            &! counter on Q angular momenta
     &       l,             &! counter on angular momenta
     &       jv,            &! beta function counter
     &       iv,            &! beta function counter
     &       nvales(nsx),   &! number of pseudo wavefunctions
     &       ir              ! mesh points counter
!
      character*20                                                      &
     &        title,        &! name of the chemical element
     &        xctype,       &! Name of the exchange-correlation
     &        fmt*60         ! format string
!
!     We first check the input variables
!
      if (is.lt.0 .or. is.gt.nsx)                                       &
     &   call error('readvan','Wrong is number', 1)
!
      if (ipp.lt.0 .or. ipp.gt.2)                                       &
     &   call error('readvan','Wrong pseudopotential type', 1)
!
      read(iunps, *, err=100 )                                          &
     &     (iver(i),i=1,3), (idmy(i),i=1,3)
      if ( iver(1).gt.7 .or. iver(1).lt.1 .or.                          &
     &     iver(2).gt.9 .or. iver(2).lt.0 .or.                          &
     &     iver(3).gt.9 .or. iver(3).lt.0      )                        &
     &   call error('readvan','wrong version numbers',1)
!
      read( iunps, '(a20,3f15.9)', err=100, iostat=ios )                &
     &     title, z(is), zp(is), exfact 
      if ( z(is).le.0.d0.or.z(is).gt.120.d0)                            &
     &     call error( 'readvan','wrong z', is )
      if ( zp(is).le.0.d0.or.zp(is).gt.25.d0)                           &
     &     call error( 'readvan','wrong zp', is )
      if ( exfact.lt.-6.or.exfact.gt.5)                                 &
     &     call error('readvan','Wrong xc in pseudopotential',1)
! convert from "our" conventions to Vanderbilt conventions
      if (exfact.eq.-5) exfact=bp88
      if (exfact.eq.-6) exfact=pw91
!
! check that the same dft level is used by all potentials
!
      if (is.eq.1) then
         dft=exfact
      else if (dft.ne.exfact) then
         call error('readvan','inconsistent DFT',is)
      end if
!
      read( iunps, '(2i5,1pe19.11)', err=100, iostat=ios )              &
     &     nvales(is), mesh(is), etotpseu
      if ( nvales(is).gt.nvalx )                                        &
     &     call error( 'readvan', 'increase nvalx', nvales(is) )
      if ( nvales(is).lt. 1 )                                           &
     &     call error( 'readvan', 'wrong nvales ', is )
      if ( mesh(is).gt.mmaxx .or. mesh(is).lt.0 )                       &
     &     call error( 'readvan','wrong mesh', is )
!
!     nnlz, wwnl, ee  give info (not used) on pseudo eigenstates
!
      read( iunps, '(i5,2f15.9)', err=100, iostat=ios ) ( nnlz(iv),     &
     &                      wwnl(iv), ee(iv), iv=1,nvales(is) )
      read( iunps, '(2i5,f15.9)', err=100, iostat=ios ) keyps,          &
     &                      ifpcor(is), rinner(1,is)
!
!     ifpcor   1 if "partial core correction" of louie, froyen,
!                 & cohen to be used; 0 otherwise
!
!     keyps= 0 --> standard hsc pseudopotential with exponent 4.0
!            1 --> standard hsc pseudopotential with exponent 3.5
!            2 --> vanderbilt modifications using defaults
!            3 --> new generalized eigenvalue pseudopotentials
!            4 --> frozen core all-electron case
      if ( keyps .lt. 0 .or. keyps .gt. 4 ) then
         call error('readvan','wrong keyps',keyps)
      else if (keyps.eq.4) then
         call error('readvan','keyps not implemented',keyps)
      else if (keyps.lt.3 .and. ipp.lt.2) then
         call error('readvan','HSC read, ultrasoft expected',keyps)
      else if (keyps.eq.3 .and. ipp.eq.2) then
         call error('readvan','ultrasoft read, HSC expected',keyps)
      end if
!
!     Read information on the angular momenta, and on Q pseudization
!     (version > 3.0)
!
      if (iver(1).ge.3) then
         read( iunps, '(2i5,f9.5,2i5,f9.5)', err=100, iostat=ios )      &
     &        nang, lloc, eloc, ifqopt(is), nqf(is), dummy
!
!    NB: In the Vanderbilt atomic code the angular momentum goes 
!        from 1 to nang
!
         if ( nang .gt. nvalx + 1 .or. nang.lt.0 )                      &
     &        call error(' readvan', 'Wrong nang', nang)
         if ( lloc .eq. -1 ) lloc = nang+1
         if ( lloc .gt. nang+1 .or. lloc .lt. 0 )                       &
     &        call error( 'readvan', 'wrong lloc', is )
         if ( nqf(is).gt.nqfx .or. nqf(is).lt.0 )                       &
     &        call error(' readvan', 'Wrong nqf', nqf(is))
         if ( ifqopt(is).lt.0 )                                         &
     &        call error( 'readvan', 'wrong ifqopt', is )
      end if
!
!     Reads and test the values of rinner (version > 5.1)
!     rinner = radius at which to cut off partial core or q_ij
!
      if (10*iver(1)+iver(2).ge.51) then
!
         read( iunps, *, err=100, iostat=ios )                          &
     &        (rinner(lp,is), lp=1,2*nang-1 )
!
         do lp = 1, 2*nang-1
            if (rinner(lp,is).lt.0.d0)                                  &
     &           call error('readvan','Wrong rinner', is )
         enddo
      else if (iver(1).gt.3) then
         do lp = 2, 2*nang-1
            rinner(lp,is)=rinner(1,is)
         end do
      end if
!
      if (iver(1).ge.4)                                                 &
     &     read( iunps, '(i5)',err=100, iostat=ios ) irel
!     
!       set the number of angular momentum terms in q_ij to read in
!
      if (iver(1).eq.1) then
         ipp  = 0
! old format: no distinction between nang and nvales
         nang = nvales(is)
! old format: no optimization of q_ij => 3-term taylor series
         nqf(is)=3
         nqlc(is)=5
      else if (iver(1).eq.2) then
         nang = nvales(is)
         nqf(is)=3
         nqlc(is) = 2*nang - 1
      else
         nqlc(is) = 2*nang - 1
      end if
!
      if ( nqlc(is).gt.lqx )                                            &
     &     call error(' readvan', 'Wrong  nqlc', nqlc(is) )
!
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                              ( rc(l), l=1,nang )
!
!     reads the number of beta functions 
!
      read( iunps, '(2i5)', err=100, iostat=ios ) nbeta(is), kkbeta(is)
!
      if( nbeta(is).gt.nbrx .or. nbeta(is).lt.0 )                       &
     &     call error( 'readvan','wrong nbeta', is )
      if( kkbeta(is).gt.mesh(is) .or. kkbeta(is).lt.0 )                 &
     &     call error( 'readvan','wrong kkbeta', is )
!
!    Now reads the main Vanderbilt parameters
!
      do iv=1,nbeta(is)
         read( iunps, '(i5)',err=100, iostat=ios ) lll(iv,is)
         read( iunps, '(1p4e19.11)',err=100, iostat=ios ) eee(iv),      &
     &                      ( betar(ir,iv,is), ir=1,kkbeta(is) )
         if ( lll(iv,is).gt.3 .or. lll(iv,is).lt.0 )                    &
     &      call error( 'readvan', ' wrong lll ? ', is )
         do jv=iv,nbeta(is)
            read( iunps, '(1p4e19.11)', err=100, iostat=ios )           &
     &           dion(iv,jv,is),                                        &
     &           ddd(iv,jv), qqq(iv,jv,is),                             &
     &           (qfunc(ir,iv,jv,is),ir=1,kkbeta(is)),                  &
     &           ((qfcoef(i,lp,iv,jv,is),i=1,nqf(is)),lp=1,nqlc(is))
!
!     Use the symmetry of the coefficients
!
            dion(jv,iv,is)=dion(iv,jv,is)
            qqq(jv,iv,is)=qqq(iv,jv,is)
!
            do ir = 1, kkbeta(is)
               qfunc(ir,jv,iv,is)=qfunc(ir,iv,jv,is)
            enddo
!
            do i = 1, nqf(is)
               do lp= 1, nqlc(is)
                  qfcoef(i,lp,jv,iv,is)=qfcoef(i,lp,iv,jv,is)
               enddo
            enddo
         enddo
      enddo
!
!    for versions later than 7.2
!
      if (10*iver(1)+iver(2).ge.72) then
         read( iunps, '(6i5)',err=100, iostat=ios )                     &
     &        (iptype(iv), iv=1,nbeta(is))
         read( iunps, '(i5,f15.9)',err=100, iostat=ios )                &
     &        npf, dummy
      end if
!
!   reads the local potential of the vanderbilt
!   Note:    For consistency with HSC pseudopotentials etc., the
!            program carries a bare local potential rucore for
!            each l value. For keyps=3 they are all the same.
!
      read( iunps, '(1p4e19.11)',err=100, iostat=ios ) rcloc,           &
     &                       ( rucore(ir,1,is), ir=1,mesh(is) )
!
!   If present reads the core charge
!
      if ( ifpcor(is).eq.1 ) then 
         if (iver(1).ge.7)                                              &
     &        read( iunps, '(1p4e19.11)', err=100, iostat=ios )         &
     &        dummy
         read( iunps, '(1p4e19.11)', err=100, iostat=ios )              &
     &                         ( rscore(ir,is), ir=1,mesh(is) )
      endif
!
!     Reads the screened local potential
!
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                             (vloc0(ir), ir=1,mesh(is))
!
!     Reads the valence atomic charge
!
      read( iunps, '(1p4e19.11)', err=100, iostat=ios )                 &
     &                             (rsatom(ir), ir=1,mesh(is))
!
!     Reads the logarithmic mesh (if version > 1)
!
      if (iver(1).gt.1) then
         read( iunps, '(1p4e19.11)',err=100, iostat=ios )               &
     &        (r(ir,is),ir=1,mesh(is))
         read( iunps, '(1p4e19.11)',err=100, iostat=ios )               &
     &        (rab(ir,is),ir=1,mesh(is))
      else
!
!     generate herman-skillman mesh (if version = 1)
!
         call herman_skillman_grid(mesh(is),z(is),cmesh(is),r(1,is))
      end if
!
!!!      if ( ifpcor(is).eq.1 ) then 
!
!!!    For compatibility with rho_core in the non-US case
!!!         do ir=1,mesh(is)
!!!            rscore(ir,is) = rscore(ir,is)/4.0/3.14159265/r(ir,is)**2
!!!         enddo
!!!      end if
!
!     Copy the local potential on all the channels
!
      do iv = 2, nbeta(is)
         do ir = 1, mesh(is)
            rucore(ir,iv,is)=rucore(ir,1,is)
         enddo
      enddo
!
!    Reads the wavefunctions of the atom
!      
      if (iver(1).ge.7) then
         read( iunps, *, err=100, iostat=ios ) i
         if (i.ne.nvales(is))                                           &
     &        call error('readvan','unexpected or unimplemented case',1)
      end if
!
      if (iver(1).ge.6)                                                 &
     &     read( iunps, *, err=100, iostat=ios )                        &
     &     ((chi(ir,iv),ir=1,mesh(is)),iv=1,nvales(is))
!      do iv = 1, nvales(is)
!            i = nnlz(iv) / 100
!            lchi(iv,is) = nnlz(iv)/10 - i * 10
!      enddo
!
      if (iver(1).eq.1) then
!
!   old version: read the q(r) here
!
         do iv=1,nbeta(is)
            do jv=iv,nbeta(is)
               lmin=lll(jv,is)-lll(iv,is)+1
               lmax=lmin+2*lll(iv,is)
               do l=lmin,lmax
                  read(iunps,*, err=100, iostat=ios)                    &
     &                 (qrl(ir,iv,jv,l,is),ir=1,kkbeta(is))
               end do
            end do
         end do
      else
!
!   new version: fill the q(r) here
!
         do iv=1,nbeta(is)
            do jv=iv,nbeta(is)
               lmin=lll(jv,is)-lll(iv,is)+1
               lmax=lmin+2*lll(iv,is)
               do l=lmin,lmax
                  do ir=1,kkbeta(is)
                     if (r(ir,is).ge.rinner(l,is)) then
                        qrl(ir,iv,jv,l,is)=qfunc(ir,iv,jv,is)
                     else
                        qrl(ir,iv,jv,l,is)=qfcoef(1,l,iv,jv,is)
                        do i = 2, nqf(is)
                           qrl(ir,iv,jv,l,is)=qrl(ir,iv,jv,l,is) +      &
     &                          qfcoef(i,l,iv,jv,is)*r(ir,is)**(2*i-2)
                        end do
                        qrl(ir,iv,jv,l,is) = qrl(ir,iv,jv,l,is) *       &
     &                                       r(ir,is)**(l+1)
                     end if
                  end do
               end do
            end do
         end do
      end if
!
!    Here we write on output informations on the read pseudopotential
!
      call dftname(xctype)
!
      write(6,200) is
200   format (/4x,60(1h=)/4x,'|  pseudopotential report',               &
     &        ' for atomic species:',i3,11x,'|')
      write(6,300) 'pseudo potential version', iver(1),                 &
     &     iver(2), iver(3)
300   format (4x,'|  ',1a30,3i4,13x,' |' /4x,60(1h-))
      write(6,400) title, xctype
400   format (4x,'|  ',2a20,' exchange-corr  |')
      write(6,500) z(is), is, zp(is), exfact
500   format (4x,'|  z =',f5.0,4x,'zv(',i2,') =',f5.0,4x,'exfact =',    &
     &     f10.5, 9x,'|')
      write(6,600) ifpcor(is), etotpseu
 600  format (4x,'|  ifpcor = ',i2,10x,' atomic energy =',f10.5,        &
     &     ' Ry',6x,'|')
      write(6,700)
700   format(4x,'|  index    orbital      occupation    energy',14x,'|')
      write(6,800) ( iv, nnlz(iv), wwnl(iv), ee(iv), iv=1,nvales(is) )
800   format(4x,'|',i5,i11,5x,f10.2,f12.2,15x,'|')
      if (iver(1).ge.3.and.nang.gt.0) then
         write(fmt,900) 2*nang-1, 40-8*(2*nang-2)
 900     format('(4x,''|  rinner ='',',i1,'f8.4,',i2,'x,''|'')')
         write(6,fmt)  (rinner(lp,is),lp=1,2*nang-1)
      end if
      write(6,1000)
1000  format(4x,'|    new generation scheme:',32x,'|')
      write(6,1100) nbeta(is),kkbeta(is),rcloc
1100  format(4x,'|    nbeta = ',i2,5x,'kkbeta =',i5,5x,                 &
     &     'rcloc =',f10.4,4x,'|'/                                      &
     &     4x,'|    ibeta    l     epsilon   rcut',25x,'|')
      do iv = 1, nbeta(is)
         lp=lll(iv,is)+1
         write(6,1200) iv,lll(iv,is),eee(iv),rc(lp)
1200      format(4x,'|',5x,i2,6x,i2,4x,2f7.2,25x,'|')
      enddo
      write(6,1300)
1300  format (4x,60(1h=))
!
      return
100   call error('readvan','error reading pseudo file', abs(ios) )
      end
!-----------------------------------------------------------------------
      subroutine recips(a,a1,a2,a3,b1,b2,b3)
!-----------------------------------------------------------------------
!   generates the reciprocal lattice vectors b1,b2,b3 given the real
!   space vectors a1,a2,a3.  the b's are units of 2pi/a.
!
      implicit none
      real(kind=8)  a1(3),a2(3),a3(3),b1(3),b2(3),b3(3)
      real(kind=8) den,s,a
      integer i,j,k,l,iperm, ir
!
      den=0.d0
      i=1
      j=2
      k=3
      s=1.d0
    1 continue
      do iperm=1,3
         den=den+s*a1(i)*a2(j)*a3(k)
         l=i
         i=j
         j=k
         k=l
      end do
!
      i=2
      j=1
      k=3
      s=-s
      if(s.lt.0.d0) go to 1
      i=1
      j=2
      k=3
      den=a/abs(den)
!
      do ir=1,3
         b1(ir)=den*(a2(j)*a3(k)-a2(k)*a3(j))
         b2(ir)=den*(a3(j)*a1(k)-a3(k)*a1(j))
         b3(ir)=den*(a1(j)*a2(k)-a1(k)*a2(j))
         l=i
         i=j
         j=k
         k=l
      end do
!
      return
      end
!-----------------------------------------------------------------------
      subroutine rhoset(cp,phi,bephi,qbecp,nss,ist,rho)
!-----------------------------------------------------------------------
!     input: cp (non-orthonormal), phi, bephi, qbecp
!     computes the matrix
!       rho = <s'c0|s cp> = <phi|s cp>
!     where  |phi> = s'|c0> = |c0> + sum q_ij |i><j|c0>
!     where s=s(r(t+dt)) and s'=s(r(t))  
!     routine makes use of  c(-q)=c*(q)
!
      use ion_parameters
      use cvan
      use elct
!
      implicit none
!
      integer nss, ist
      complex(kind=8)   cp(ngw,n), phi(ngw,n)
      real(kind=8)       bephi(nhsa,n), qbecp(nhsa,n), rho(nx,nx)
      integer i, j
      real(kind=8)    tmp1(nx,nx) ! automatic array
!
      call zero(nx*nss,rho)
!
!     <phi|cp>
!
      call MXMA(phi(1,ist),2*ngw,1,cp(1,ist),1,2*ngw,                   &
     &                              rho,1,nx,nss,2*ngw,nss)
!
!     q >= 0  components with weight 2.0
!
      do j=1,nss
         do i=1,nss
            rho(i,j)=2.*rho(i,j)
         end do
      end do
!
      if (ng0.eq.2) then
!
!     q = 0  components has weight 1.0
!
         do j=1,nss
            do i=1,nss
               rho(i,j) = rho(i,j) -                                    &
     &              real(phi(1,i+ist-1))*real(cp(1,j+ist-1))
            end do
         end do
      end if
#ifdef PARA
      call reduce(nx*nss,rho)
#endif
!
      if(nvb.gt.0)then
         call zero(nx*nss,tmp1)
!
         call MXMA(bephi(1,ist),nhsa,1,qbecp(1,ist),1,nhsa,               &
     &                                tmp1,1,nx,nss,nhsavb,nss)
!
         do j=1,n
            do i=1,n
               rho(i,j)=rho(i,j)+tmp1(i,j)
            end do
         end do
      endif
!
      return
      end
!-------------------------------------------------------------------------
      subroutine rhov(irb,eigrb,rhovan,rhog,rhor)
!-----------------------------------------------------------------------
!     this routine calculates arrays rhog 
!
!        n_v(g) = sum_i,ij rho_i,ij q_i,ji(g) e^-ig.r_i
!
!     routine makes use of c(-g)=c*(g)  and  beta(-g)=beta*(g)
!
      use ions
      use gvec
      use cvan
      use parm
      use elct
      use gvecb
      use parmb
      use control
      use qgb_mod
      use work1
      use work_box
#ifdef PARA
      use para_mod
#endif
!
      implicit none
!
      integer irb(3,nax,nsx)
      real(kind=8)    rhor(nnr,nspin), rhovan(nat,nhx,nhx,nspin)
      complex(kind=8) eigrb(ngb,nas,nsp), rhog(ng,nspin)
!
      integer isup, isdw, nc, ic, iv, jv, ig, ijv, is, iss,             &
     &     isa, ia, ism, ibig1, ibig2, ibig3, ir1, ir2, ir3, ir, ib
      real(kind=8) ra, ra1, ra2, qvr, qvi, sum, sum1, sum2, SSUM
      complex(kind=8) ci, fp, fm, CSUM
      complex(kind=8), allocatable::  qgbt(:,:)
      complex(kind=8), pointer:: v(:)
!
      call tictac(13,0)
      ci=(0.,1.)
!
!     localised vanderbilt charge for nspin
!
      v => wrk1
      call zero(2*nnr,v)
      allocate(qgbt(ngb,2))
!
      if(nspin.eq.1)then
!     ------------------------------------------------------------------
!     case of nspin=1
!     ------------------------------------------------------------------
         iss=1
         do is=1,nvb
            isa=0
            do ism=1,is-1
               isa=isa+na(ism)
            end do
            do ia=1,na(is),2
               call zero(2*ngb,qgbt(1,1))
               call zero(2*ngb,qgbt(1,2))
               nc=2
               if( ia.eq.na(is).and.na(is)/2.ne.real(na(is))/2. ) nc=1
               do ic=1,nc
                  isa=isa+1
                  ijv=0
                  do iv= 1,nh(is)
                     do jv=iv,nh(is)
                        ijv=ijv+1
                        sum=rhovan(isa,iv,jv,iss)
                        if(iv.ne.jv) sum=2.*sum
                        do ig=1,ngb
                           qgbt(ig,ic)=qgbt(ig,ic)+sum*qgb(ig,ijv,is)
                        end do
                     end do
                  end do
               end do
!     
               call zero(2*nnrb,qv)
!
               if(nc.eq.2)then
                  do ig=1,ngb
                     qv(npb(ig))=  eigrb(ig,ia  ,is)*qgbt(ig,1)         &
     &                  + ci*      eigrb(ig,ia+1,is)*qgbt(ig,2)
                     qv(nmb(ig))=                                       &
     &                       conjg(eigrb(ig,ia  ,is)*qgbt(ig,1))        &
     &                  + ci*conjg(eigrb(ig,ia+1,is)*qgbt(ig,2))
                  end do
               else
                  do ig=1,ngb
                     qv(npb(ig)) = eigrb(ig,ia,is)*qgbt(ig,1)
                     qv(nmb(ig)) = conjg(eigrb(ig,ia,is)*qgbt(ig,1))
                  end do
               endif
!
               call ivfftb(qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!     
               if(iprsta.gt.2) then
                  ra1= real(CSUM(nnrb,qv,1))
                  ra2=aimag(CSUM(nnrb,qv,1))
                  write(6,'(a,f12.8)') ' rhov: 1-atom g-sp = ',         &
     &                 omegab*real(qgbt(1,1))
                  write(6,'(a,f12.8)') ' rhov: 1-atom r-sp = ',         &
     &                 omegab*ra1/(nr1b*nr2b*nr3b)
                  write(6,'(a,f12.8)') ' rhov: 1-atom g-sp = ',         &
     &                 omegab*real(qgbt(1,2))
                  write(6,'(a,f12.8)') ' rhov: 1-atom r-sp = ',         &
     &                 omegab*ra2/(nr1b*nr2b*nr3b)
               endif
!
               do ir3=1,nr3b
                  ibig3=irb(3,ia,is)+ir3-1
                  ibig3=1+mod(ibig3-1,nr3)
                  if(ibig3.lt.1.or.ibig3.gt.nr3)                        &
     &                 call error('rhov','ibig3 wrong',ibig3)
#ifdef PARA
                  ibig3=ibig3-n3(me)
                  if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                  do ir2=1,nr2b
                     ibig2=irb(2,ia,is)+ir2-1
                     ibig2=1+mod(ibig2-1,nr2)
                     if(ibig2.lt.1.or.ibig2.gt.nr2)                     &
     &                 call error('rhov','ibig2 wrong',ibig2)
                     do ir1=1,nr1b
                        ibig1=irb(1,ia,is)+ir1-1
                        ibig1=1+mod(ibig1-1,nr1)
                        if(ibig1.lt.1.or.ibig1.gt.nr1)                  &
     &                       call error('rhov','ibig1 wrong',ibig1)
                        ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                        ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                        qvr=real(qv(ir))
                        v(ib) = v(ib)+cmplx(qvr,0.0)
                     end do
                  end do
#ifdef PARA
                  end if
#endif
               end do
!
               if(nc.eq.2)then
                  do ir3=1,nr3b
                     ibig3=irb(3,ia+1,is)+ir3-1
                     ibig3=1+mod(ibig3-1,nr3)
                     if(ibig3.lt.1) ibig3=ibig3+nr3
                     if(ibig3.lt.1.or.ibig3.gt.nr3)                     &
     &                    call error('rhov','ibig3 wrong',ibig3)
#ifdef PARA
                     ibig3=ibig3-n3(me)
                     if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                     do ir2=1,nr2b
                        ibig2=irb(2,ia+1,is)+ir2-1
                        ibig2=1+mod(ibig2-1,nr2)
                        if(ibig2.lt.1) ibig2=ibig2+nr2
                        if(ibig2.lt.1.or.ibig2.gt.nr2)                  &
     &                       call error('rhov','ibig2 wrong',ibig2)
                        do ir1=1,nr1b
                           ibig1=irb(1,ia+1,is)+ir1-1
                           ibig1=1+mod(ibig1-1,nr1)
                           if(ibig1.lt.1) ibig1=ibig1+nr1
                           if(ibig1.lt.1.or.ibig1.gt.nr1)               &
     &                          call error('rhov','ibig1 wrong',ibig1)
                           ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                           ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                           qvr=aimag(qv(ir))
                           v(ib) = v(ib)+cmplx(qvr,0.0)
                        end do
                     end do
#ifdef PARA
                     end if
#endif
                  end do
               endif
!
            end do
         end do
!
         do ir=1,nnr
            rhor(ir,iss)=rhor(ir,iss)+real(v(ir))        
         end do
!
         if(iprsta.gt.2) then
            ra =SSUM(2*nnr,v,1)
#ifdef PARA
            call reduce(1,ra)
#endif
            write(6,'(a,f12.8)')                                        &
     &           ' rhov: int  n_v(r)  dr = ',omega*ra/(nr1*nr2*nr3)
         endif
!
         call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!
         if(iprsta.gt.2) then
            write(6,*) ' rhov: smooth ',omega*rhog(1,iss)
            write(6,*) ' rhov: vander ',omega*v(1)
            write(6,*) ' rhov: all    ',omega*(rhog(1,iss)+v(1))
         endif
!
         do ig=1,ng
            rhog(ig,iss)=rhog(ig,iss)+v(np(ig))
         end do
!
         if(iprsta.gt.1) write(6,'(a,2f12.8)')                          &
     &        ' rhov: n_v(g=0) = ',omega*real(rhog(1,iss))
!
         if(iprsta.gt.4) then
            isa=0
            do is=1,nsp
               do ia=1,na(is)
                  isa=isa+1
                  write(6,*)
                  write(6,'(33x,a,i4)') '    rhovan(isa) ',isa
                  do iv=1,nh(is)
                     write(6,'(8f9.4)')                                 &
     &                    (rhovan(isa,iv,jv,iss),jv=1,nh(is))
                  end do
               end do
            end do
         endif
!
      else
!     ------------------------------------------------------------------
!     case of nspin=2
!     ------------------------------------------------------------------
         isup=1
         isdw=2
         do is=1,nvb
            isa=0
            do ism=1,is-1
               isa=isa+na(ism)
            end do
            do ia=1,na(is)
               call zero(2*ngb,qgbt(1,1))
               call zero(2*ngb,qgbt(1,2))
               isa=isa+1
               ijv=0
               do iv=1,nh(is)
                  do jv=iv,nh(is)
                     ijv=ijv+1
                     sum1=rhovan(isa,iv,jv,isup)
                     sum2=rhovan(isa,iv,jv,isdw)
                     if(iv.ne.jv) sum1=2.*sum1
                     if(iv.ne.jv) sum2=2.*sum2
                     do ig=1,ngb
                        qgbt(ig,1)=qgbt(ig,1)+sum1*qgb(ig,ijv,is)
                        qgbt(ig,2)=qgbt(ig,2)+sum2*qgb(ig,ijv,is)
                     end do
                  end do
               end do
!
               call zero(2*nnrb,qv)
!
               do ig=1,ngb
                  qv(npb(ig)) =    eigrb(ig,ia,is)*qgbt(ig,1)           &
     &                  + ci*      eigrb(ig,ia,is)*qgbt(ig,2)
                  qv(nmb(ig)) = conjg(eigrb(ig,ia,is)*qgbt(ig,1))       &
     &                  + ci*   conjg(eigrb(ig,ia,is)*qgbt(ig,2))
               end do
!
               call ivfftb(qv,nr1b,nr2b,nr3b,nr1bx,nr2bx,nr3bx)
!
               if(iprsta.gt.2) then
                  ra1= real(CSUM(nnrb,qv,1))
                  ra2=aimag(CSUM(nnrb,qv,1))
                  write(6,'(a,f12.8)') ' rhov: up   g-space = ',        &
     &                 omegab*real(qgbt(1,1))
                  write(6,'(a,f12.8)') ' rhov: up r-sp = ',             &
     &                 omegab*ra1/(nr1b*nr2b*nr3b)
                  write(6,'(a,f12.8)') ' rhov: dw g-space = ',          &
     &                 omegab*real(qgbt(1,2))
                  write(6,'(a,f12.8)') ' rhov: dw r-sp = ',             &
     &                 omegab*ra2/(nr1b*nr2b*nr3b)
               endif
!     
               do ir3=1,nr3b
                  ibig3=irb(3,ia,is)+ir3-1
                  ibig3=1+mod(ibig3-1,nr3)
                  if(ibig3.lt.1.or.ibig3.gt.nr3)                        &
     &                 call error('rhov','ibig3 wrong',ibig3)
#ifdef PARA
                  ibig3=ibig3-n3(me)
                  if (ibig3.gt.0.and.ibig3.le.npp(me)) then
#endif
                  do ir2=1,nr2b
                     ibig2=irb(2,ia,is)+ir2-1
                     ibig2=1+mod(ibig2-1,nr2)
                     if(ibig2.lt.1.or.ibig2.gt.nr2)                     &
     &                 call error('rhov','ibig2 wrong',ibig2)
                     do ir1=1,nr1b
                        ibig1=irb(1,ia,is)+ir1-1
                        ibig1=1+mod(ibig1-1,nr1)
                        if(ibig1.lt.1.or.ibig1.gt.nr1)                  &
     &                       call error('rhov','ibig1 wrong',ibig1)
                        ir=ir1+(ir2-1)*nr1bx+(ir3-1)*nr1bx*nr2bx
                        ib=ibig1+(ibig2-1)*nr1x+(ibig3-1)*nr1x*nr2x
                        qvr=real(qv(ir))
                        qvi=aimag(qv(ir))
                        v(ib) = v(ib)+cmplx(qvr,qvi)
                     end do
                  end do
#ifdef PARA
                  end if
#endif
               end do
!
            end do
         end do
!
         do ir=1,nnr
            rhor(ir,isup)=rhor(ir,isup)+real(v(ir)) 
            rhor(ir,isdw)=rhor(ir,isdw)+aimag(v(ir)) 
         end do
!
         if(iprsta.gt.2) then
            ra2 = SSUM(2*nnr,v,1)
#ifdef PARA
            call reduce(1,ra2)
#endif
            write(6,'(a,f12.8)') 'rhov:in n_v  ',omega*ra2/(nr1*nr2*nr3)
         endif
!
         call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!
         if(iprsta.gt.2) then
            write(6,*) 'rhov: smooth up',omega*rhog(1,isup)
            write(6,*) 'rhov: smooth dw',omega*rhog(1,isdw)
            write(6,*) 'rhov: vander up',omega*real(v(1))
            write(6,*) 'rhov: vander dw',omega*aimag(v(1))
            write(6,*) 'rhov: all up',                                  &
     &           omega*(rhog(1,isup)+real(v(1)))
            write(6,*) 'rhov: all dw',                                  &
     &           omega*(rhog(1,isdw)+aimag(v(1)))
         endif
!
         do ig=1,ng
            fp=  v(np(ig)) + v(nm(ig))
            fm=  v(np(ig)) - v(nm(ig))
            rhog(ig,isup)=rhog(ig,isup) + 0.5*cmplx(real(fp),aimag(fm))
            rhog(ig,isdw)=rhog(ig,isdw) + 0.5*cmplx(aimag(fp),-real(fm))
         end do
!
         if(iprsta.gt.2) write(6,'(a,2f12.8)')                          &
     &        ' rhov: n_v(g=0) up   = ',omega*real (rhog(1,isup))
         if(iprsta.gt.2) write(6,'(a,2f12.8)')                          &
     &        ' rhov: n_v(g=0) down = ',omega*real(rhog(1,isdw))
!
         if(iprsta.gt.4) then
            isa=0
            do is=1,nsp
               do ia=1,na(is)
                  isa=isa+1
                  write(6,*)
                  do iss = 1,2
                     write(6,'(33x,a,2i4)')                             &
     &                    '  rhovan(isa,iss)  isa iss ',isa,iss
                     do iv=1,nh(is)
                        write(6,'(8f9.4)')                              &
     &                       (rhovan(isa,iv,jv,iss),jv=1,nh(is))
                     end do
                  end do
               end do
            end do
         endif
!     
      endif
      deallocate(qgbt)
      call tictac(13,1)
!
      return
      end
!-----------------------------------------------------------------------
      subroutine set_fft_dim ( nr1, nr2, nr3, nr1s, nr2s, nr3s, nnr,    &
     &                         nr1x,nr2x,nr3x,nr1sx,nr2sx,nr3sx,nnrs ) 
!-----------------------------------------------------------------------
!
! Set the dimensions of fft arrays (for the scalar case).
! These may differ from fft transform lenghts for computational
! efficiency, thus avoiding or minimizing memory conflicts.
! See machine-dependent function "good_fft_dimension".
! fft arrays are effectively stored as one-dimensional arrays,
! but their memory layout is the same as that of a 3d array
! having dimensions (nr1x,nr2x,nr3x) and the like for smooth grid
!
      implicit none
! input
      integer nr1, nr2, nr3, nr1s, nr2s, nr3s
! output 
      integer nnr, nr1x, nr2x, nr3x, nr1sx, nr2sx, nr3sx, nnrs
!
!     dense grid
!
      integer good_fft_dimension
      external good_fft_dimension
!
      nr1x = good_fft_dimension(nr1 )
      nr2x = nr2
      nr3x = nr3
      nnr  = nr1x*nr2x*nr3x
!
!     smooth grid
!
      nr1sx= good_fft_dimension(nr1s)
      nr2sx= nr2s
      nr3sx= nr3s
      nnrs = nr1sx*nr2sx*nr3sx
!
      return
      end
!
!-------------------------------------------------------------------------
      subroutine set_fft_grid(a1,a2,a3,alat,gcut,nr1,nr2,nr3)
!-------------------------------------------------------------------------
!
      implicit none
! input
      real(kind=8) a1(3), a2(3), a3(3), alat, gcut
! input/output
      integer nr1, nr2, nr3
! local
      integer minr1, minr2, minr3, good_fft_order
      real(kind=8) a1n, a2n, a3n
      external good_fft_order
!
!
      a1n=sqrt(a1(1)**2+a1(2)**2+a1(3)**2)/alat
      a2n=sqrt(a2(1)**2+a2(2)**2+a2(3)**2)/alat
      a3n=sqrt(a3(1)**2+a3(2)**2+a3(3)**2)/alat
!
      minr1=int(2*sqrt(gcut)*a1n+1.)
      minr2=int(2*sqrt(gcut)*a2n+1.)
      minr3=int(2*sqrt(gcut)*a3n+1.)
!
      write(6,1010) gcut,minr1,minr2,minr3
1010  format(' set_fft_grid: gcut = ',f7.2,'  n1,n2,n3 min =',3i4)
      if (nr1.le.0) nr1=minr1
      if (nr2.le.0) nr2=minr2
      if (nr3.le.0) nr3=minr3
      nr1=good_fft_order(nr1)
      nr2=good_fft_order(nr2)
      nr3=good_fft_order(nr3)
      if (minr1-nr1.gt.2)                                               &
     &     call error('set_fft_grid','n1 too small ?',nr1)
      if (minr2-nr2.gt.2)                                               &
     &     call error('set_fft_grid','n2 too small ?',nr2)
      if (minr3-nr3.gt.2)                                               &
     &     call error('set_fft_grid','n3 too small ?',nr3)
!
      return
      end
!
!-------------------------------------------------------------------------
      subroutine sigset(cp,becp,qbecp,nss,ist,sig)
!-----------------------------------------------------------------------
!     input: cp (non-orthonormal), becp, qbecp
!     computes the matrix
!       sig = 1 - a ,  a = <cp|s|cp> = <cp|cp> + sum q_ij <cp|i><j|cp>
!     where s=s(r(t+dt)) 
!     routine makes use of c(-q)=c*(q)
!
      use ion_parameters
      use cvan
      use elct
!
      implicit none
!
      integer nss, ist
      complex(kind=8)  cp(ngw,n)
      real(kind=8) becp(nhsa,n), qbecp(nhsa,n), sig(nx,nx)
!
      integer i, j
      real(kind=8)    tmp1(nx,nx) ! automatic array
!
      call zero(nx*nss,sig )
      call MXMA(cp(1,ist),2*ngw,1,cp(1,ist),1,2*ngw,                    &
     &                                  sig,1,nx,nss,2*ngw,nss)
!
!     q >= 0  components with weight 2.0
!
      do j=1,nss
         do i=1,nss
            sig(i,j)=-2.*sig(i,j)
         end do
      end do
      if (ng0.eq.2) then
!
!     q = 0  components has weight 1.0
!
         do j=1,nss
            do i=1,nss
               sig(i,j) = sig(i,j) +                                    &
     &              real(cp(1,i+ist-1))*real(cp(1,j+ist-1))
            end do
         end do
      end if
#ifdef PARA
      call reduce(nx*nss,sig)
#endif
      do i=1,nss
         sig(i,i) = sig(i,i)+1.
      end do
!
      if(nvb.gt.0)then
         call zero(nx*nss,tmp1)
!
         call MXMA(becp(1,ist),nhsa,1,qbecp(1,ist),1,nhsa,                &
     &                              tmp1,1,nx,nss,nhsavb,nss)
!
         do j=1,nss
            do i=1,nss
               sig(i,j)=sig(i,j)-tmp1(i,j)
            end do
         end do
      endif
!
      return
      end
!-----------------------------------------------------------------------
      subroutine simpson( mesh, func, rab, intg )
!-----------------------------------------------------------------------
!
!    This routine computes the integral of a function defined on a
!    logaritmic mesh, by using the open simpson formula given on 
!    pag. 109 of Numerical Recipes. In principle it is used to
!    perform integrals from zero to infinity. The first point of
!    the function should be the closest to zero but not the value
!    in zero. The formula used here automatically includes the 
!    contribution from the zero point and no correction is required.
!
!    At least 8 integrating points are required.
!
!
!    last revised 12 May 1995 by Andrea Dal Corso
!
!
      implicit none
!
!    first the dummy variables
!
      integer                                                           &
     &       mesh                   ! input: the number of mesh points
      real(kind=8)                                                      &
     &       func(mesh),           &! input: the function to be integrated
     &       rab(mesh),            &! input: the derivative of the log mesh 
     &       intg                   ! output: the value of the integral
!
!    the parameters
!
      real(kind=8)                                                      &
     &      c(4)                    ! the coefficients of the expansion
!
!    here the local variables
!
      integer                                                           &
     &       i
!
      if ( mesh .lt. 8 )                                                &
     &   call error('simpson','few mesh points',8)
!
      c(1) = 109.0 / 48.d0
      c(2) = -5.d0 / 48.d0
      c(3) = 63.d0 / 48.d0
      c(4) = 49.d0 / 48.d0 
!
      intg = ( func(1)*rab(1) + func(mesh  )*rab(mesh  ) )*c(1)         &
     &     + ( func(2)*rab(2) + func(mesh-1)*rab(mesh-1) )*c(2)         &
     &     + ( func(3)*rab(3) + func(mesh-2)*rab(mesh-2) )*c(3)         &
     &     + ( func(4)*rab(4) + func(mesh-3)*rab(mesh-3) )*c(4)
      do i=5,mesh-4
         intg = intg + func(i)*rab(i)
      end do
!
      return
      end
!-------------------------------------------------------------------------
      subroutine strucf (ei1,ei2,ei3,sfac)
!-----------------------------------------------------------------------
! computes the structure factor sfac(ngs,nsp) and returns it in "pseu"
!
!
      use gvec
      use gvecs
      use parm
      use cnst
      use ions
!
      implicit none
      complex(kind=8) ei1(-nr1:nr1,nas,nsp), ei2(-nr2:nr2,nas,nsp),          &
     &           ei3(-nr3:nr3,nas,nsp), sfac(ngs,nsp)
      integer is, ig, ia
!
      call tictac(3,0)
!
      do is=1,nsp
         do ig=1,ngs
            sfac(ig,is)=(0.,0.)
         end do
         do ia=1,na(is)
            do ig=1,ngs
               sfac(ig,is)=sfac(ig,is) + ei1(in1p(ig),ia,is)            &
     &                                  *ei2(in2p(ig),ia,is)            &
     &                                  *ei3(in3p(ig),ia,is)
            end do
         end do
      end do
!
      call tictac(3,1)
!
      return
      end
!-------------------------------------------------------------------------
      subroutine tauset(phi,bephi,qbephi,nss,ist,tau)
!-----------------------------------------------------------------------
!     input: phi
!     computes the matrix
!        tau = <s'c0|s|s'c0> = <phi|s|phi>,  where  |phi> = s'|c0>
!     where s=s(r(t+dt)) and s'=s(r(t))  
!     routine makes use of c(-q)=c*(q)
!
      use ion_parameters
      use cvan
      use elct
!
      implicit none
      integer nss, ist
      complex(kind=8) phi(ngw,n)
      real(kind=8)  bephi(nhsa,n), qbephi(nhsa,n), tau(nx,nx)
      integer i, j
      real(kind=8)    tmp1(nx,nx) ! automatic array
!
      call zero(nx*nss,tau)
      call MXMA(phi(1,ist),2*ngw,1,phi(1,ist),1,2*ngw,                  &
     &                                   tau,1,nx,nss,2*ngw,nss)
!
!     q >= 0  components with weight 2.0
!
      do j=1,nss
         do i=1,nss
            tau(i,j)=2.*tau(i,j)
         end do
      end do
      if (ng0.eq.2) then
!
!     q = 0  components has weight 1.0
!
         do j=1,nss
            do i=1,nss
               tau(i,j) = tau(i,j) -                                    &
     &              real(phi(1,i+ist-1))*real(phi(1,j+ist-1))
            end do
         end do
      end if
#ifdef PARA
      call reduce(nx*nss,tau)
#endif
!
      if(nvb.gt.0)then
         call zero(nx*nss,tmp1)
!
         call MXMA(bephi(1,ist),nhsa,1,qbephi(1,ist),1,nhsa,              &
     &                                     tmp1,1,nx,nss,nhsavb,nss)
!
         do j=1,nss
            do i=1,nss
               tau(i,j)=tau(i,j)+tmp1(i,j)
            end do
         end do
      endif
!
      return
      end
!
!-------------------------------------------------------------------------
      subroutine updatc(ccc,x0,phi,bephi,becp,bec,cp)
!-----------------------------------------------------------------------
!     input ccc : dt**2/emass (unchanged in output)
!     input x0  : converged lambdas from ortho-loop (unchanged in output)
!     input cp  : non-orthonormal cp=c0+dh/dc*ccc
!     input bec : <cp|beta_i>
!     input phi 
!     output cp : orthonormal cp=cp+lambda*phi
!     output bec: bec=becp+lambda*bephi
!
      use ions
      use cvan
      use elct
      use work2
      use control
!
      implicit none
!
      complex(kind=8) cp(ngw,n), phi(ngw,n)
      real(kind=8)   bec(nhsa,n), x0(nx,nx), ccc
      real(kind=8)   bephi(nhsa,n), becp(nhsa,n)
! local variables
      integer i, j, ig, is, iv, ia, inl
      real(kind=8) wtemp(n,nhsa) ! automatic array
!
!     lagrange multipliers
!
      call zero(2*ngw*n,wrk2)
      do j=1,n
         call SCAL(n,ccc,x0(1,j),1)
      end do
!
!     wrk2 = sum_m lambda_nm s(r(t+dt))|m>
!
      call MXMA(phi,1,2*ngw,x0,nx,1,wrk2,1,2*ngw,2*ngw,n,n)
!
      do i=1,n
         do ig=1,ngw
            cp(ig,i)=cp(ig,i)+wrk2(ig,i)
         end do
      end do
!    
!     updating of the <beta|c(n,g)>
!
!     bec of vanderbilt species are updated 
!
      if(nvb.gt.0)then
         call MXMA(x0,1,nx,bephi,nhsa,1,wtemp,1,n,n,n,nhsavb)
!
         do i=1,n
            do inl=1,nhsavb
               bec(inl,i)=wtemp(i,inl)+becp(inl,i)
            end do
         end do
      endif
!
      if (iprsta.gt.2) then
         write(6,*)
         do is=1,nsp
            if(nsp.gt.1) then
               write(6,'(33x,a,i4)') ' updatc: bec (is)',is
               write(6,'(8f9.4)')                                       &
     &            ((bec(ish(is)+(iv-1)*na(is)+1,i),iv=1,nh(is)),i=1,n)
            else
               do ia=1,na(is)
                  write(6,'(33x,a,i4)') ' updatc: bec (ia)',ia
                  write(6,'(8f9.4)')                                    &
     &            ((bec(ish(is)+(iv-1)*na(is)+ia,i),iv=1,nh(is)),i=1,n)
               end do
            end if
            write(6,*)
         end do
      endif
!
      do j=1,n
         call SCAL(n,1.0/ccc,x0(1,j),1)
      end do
!
      return
      end
!______________________________________________________________________
      subroutine ggablyp4(nspin,rhog,gradr,rhor,exc)
!     _________________________________________________________________
!     becke-lee-yang-parr gga
!     
!     exchange: becke, pra 38, 3098 (1988) but derived from
!        pw91 exchange formula given in prb 48, 14944 (1993)
!        by setting "b3" and "b4" to 0.0
!     correlation: miehlich et al., cpl 157, 200 (1989) 
!     method by ja white & dm bird, prb 50, 4954 (1994) 
!
!     spin-polarized version by andras stirling 10/1998,
!     using original gga program of alfredo pasquarello 22/09/1994
!     and spin-unpolarized blyp routine of olivier parisel and 
!     alfredo pasquarello (02/1997)  
!
      use gvec
      use parm
      use cnst
!
      implicit none
! input
      integer nspin
      complex(kind=8) rhog(ng,nspin)
      real(kind=8) gradr(nnr,3,nspin), rhor(nnr,nspin)
! output
! on output: rhor contains the exchange-correlation potential
      real(kind=8)  exc
! local
      integer isdw, isup, isign, ir
!
      real(kind=8) abo, agdr, agdr2, agr, agr2, agur, agur2, arodw,     &
     &      arodw2, aroe, aroe2, aroup, aroup2, ax
      real(kind=8) byagdr, byagr, byagur, cden, cf, cl1, cl11, cl2,     &
     &      cl21, cl22, cl23, cl24, cl25, cl26, cl27, clyp, csum
      real(kind=8) dddn, dexcdg, dexcdgd, dexcdgu, df1d, df1u, df2d,    &
     &     df2u, dfd, dfnum1d, dfnum1u, dfnum2d, dfnum2u, dfs, dfu,     &
     &     dfxdd, dfxdg, dfxdgd, dfxdgu, dfxdu, dilta, dilta119, dl1dn, &
     &     dl1dnd, dl1dnu, dl2dd, dl2dg, dl2dgd, dl2dgu, dl2dn,         &
     &     dl2dnd, dl2dnd1, dl2dnu, dl2dnu1, dl2do, dlt, dodn,          &
     &     dsign, dwsign, dys, dysd, dysu
      real(kind=8) ex, excupdt, exd, exu, fac1, fac2, factor1, factor2, &
     &     fx, fxd, fxden, fxdend, fxdenu, fxnum, fxnumd, fxnumu, fxu
      real(kind=8) gkf, gkfd, gkfu, grdx, grdy, grdz, grux, gruy, gruz, &
     &     grx, gry, grz
      real(kind=8) omiga, pd, pi2, pider1, pider2, pider3, pider4, piexch, pu
      real(kind=8) rhodw, rhoup, roe, roedth, roeth, roeuth, rometh
      real(kind=8) s, s2, sd, sd2, sddw, sdup, su, su2, sysl, sysld, syslu
      real(kind=8) t113, upsign, usign
      real(kind=8) x1124, x113, x118, x13, x143, x16, x19, x23, x43,    &
     &     x4718, x524, x53, x672, x718, x76, x772, x83
      real(kind=8) ys, ysd, ysl, ysld, yslu, ysr, ysrd, ysru, ysu   
!===========================================================================
      real(kind=8) bb1, bb2, bb5, aa, bb, cc, dd, delt, eps
      parameter(bb1=0.19644797,bb2=0.2742931,bb5=7.79555418,            &
     &     aa=0.04918,                                                  &
     &     bb=0.132,cc=0.2533,dd=0.349,delt=1.0e-12,eps=1.0e-14)
      complex(kind=8) ci
!
!
      ci=(0.0,1.0)
      x13=1.0/3.0
      x19=1.0/9.0
      x23=2.0/3.0
      x43=4.0/3.0
      x53=5.0/3.0
      x16=1.0/6.0
      x76=7.0/6.0
      x83=8.0/3.0
      x113=11.0/3.0
      x4718=47.0/18.0
      x718=7.0/18.0
      x118=1.0/18.0
      x524=5.0/24.0
      x1124=11.0/24.0
      x143=14.0/3.0
      x772=7.0/72.0
      x672=6.0/72.0
!     
!     _________________________________________________________________
!     derived parameters from pi
!
      pi2=pi*pi
      ax=-0.75*(3.0/pi)**x13
      piexch=-0.75/pi 
      pider1=(0.75/pi)**x13
      pider2=(3.0*pi2)**x13
      pider3=(3.0*pi2/16.0)**x13
      pider4=-1.0/(6.0*pi2)
      cf=0.3*pider2*pider2
!     _________________________________________________________________
!     other parameters
!
      t113=2.0**x113
!
      rhodw=0.0
      grdx=0.0
      grdy=0.0
      grdz=0.0
!
      fac1=1.0
!     _________________________________________________________________
!     main loop
!
      isup=1
      isdw=2
      do ir=1,nnr
         rhoup=rhor(ir,isup)
         grux=gradr(ir,1,isup)
         gruy=gradr(ir,2,isup)
         gruz=gradr(ir,3,isup)
         if(nspin.eq.2) then
            rhodw=rhor(ir,isdw)
            grdx=gradr(ir,1,isdw)
            grdy=gradr(ir,2,isdw)
            grdz=gradr(ir,3,isdw)
         else
            rhodw=0.0
            grdx =0.0
            grdy =0.0
            grdz =0.0
         endif
         roe=rhoup+rhodw
         if(roe.eq.0.0) goto 100
         aroup=abs(rhoup)
         arodw=abs(rhodw)
         aroe=abs(roe)
         grx=grux + grdx
         gry=gruy + grdy
         grz=gruz + grdz
         agur2=grux*grux+gruy*gruy+gruz*gruz
         agur=sqrt(agur2)
         agdr2=grdx*grdx+grdy*grdy+grdz*grdz
         agdr=sqrt(agdr2)
         agr2=grx*grx+gry*gry+grz*grz
         agr=sqrt(agr2)
         roeth=aroe**x13
         rometh=1.0/roeth
         gkf=pider2*roeth
         sd=1.0/(2.0*gkf*aroe)
         s=agr*sd
         s2=s*s
!     _________________________________________________________________
!     exchange 
!
         if(nspin.eq.1) then
!
!
            ysr=sqrt(1.0+bb5*bb5*s2)
            ys=bb5*s+ysr
            ysl=log(ys)*bb1
            sysl=s*ysl
            fxnum=1.0+sysl+bb2*s2
            fxden=1.0/(1.0+sysl)
            fx=fxnum*fxden
!
            ex=ax*fx*roeth*aroe
!
!     ### potential contribution ###
!
            dys=bb5*(1.0+bb5*s/ysr)/ys
            dfs=-fxnum*(ysl+bb1*s*dys)*fxden*fxden                      &
     &           +(ysl+bb1*s*dys+2.0*s*bb2)*fxden
            dfxdu=(ax*roeth*x43)*(fx-dfs*s)
            dfxdg=ax*roeth*dfs*sd
!
!     ### end of potential contribution ###
!     
         else
!
            roeuth=(2.0*aroup)**x13
            roedth=(2.0*arodw)**x13
            gkfu=pider2*roeuth*aroup
            gkfd=pider2*roedth*arodw
            upsign=sign(1.d0,gkfu-eps)
            dwsign=sign(1.d0,gkfd-eps)
            factor1=0.5*(1+upsign)/(gkfu+(1-upsign)*eps)
            fac1=gkfu*factor1
            factor2=0.5*(1+dwsign)/(gkfd+(1-dwsign)*eps)
            fac2=gkfd*factor2
            sdup=1.0/2.0*factor1
            sddw=1.0/2.0*factor2
            su=agur*sdup
            su2=su*su
            sd=agdr*sddw
            sd2=sd*sd
!
            ysru=sqrt(1.0+bb5*bb5*su2)
            ysu=bb5*su+ysru
            yslu=log(ysu)*bb1
            syslu=su*yslu
            fxnumu=1.0+syslu+bb2*su2
            fxdenu=1.0/(1.0+syslu)
            fxu=fxnumu*fxdenu
            exu=piexch*2.0*gkfu*fxu*fac1
!
            ysrd=sqrt(1.0+bb5*bb5*sd2)
            ysd=bb5*sd+ysrd
            ysld=log(ysd)*bb1
            sysld=sd*ysld
            fxnumd=1.0+sysld+bb2*sd2
            fxdend=1.0/(1.0+sysld)
            fxd=fxnumd*fxdend
            exd=piexch*2.0*gkfd*fxd*fac2
!
            ex=0.5*(exu+exd)
!
!     ### potential contribution ###
!
            dysu=bb5*(1.0+bb5*su/ysru)/ysu
            pu=2.0*su*bb2
            dfnum1u=yslu+bb1*su*dysu+pu
            df1u=dfnum1u*fxdenu
            dfnum2u=fxnumu*(yslu+bb1*su*dysu)
            df2u=dfnum2u*fxdenu*fxdenu
            dfu=df1u-df2u
            dfxdu=ax*roeuth*x43*1.0*(fxu-dfu*su)*fac1
            dfxdgu=ax*aroup*roeuth*dfu*sdup*fac1
!
            dysd=bb5*(1.0+bb5*sd/ysrd)/ysd
            pd=2.0*sd*bb2
            dfnum1d=ysld+bb1*sd*dysd+pd
            df1d=dfnum1d*fxdend
            dfnum2d=fxnumd*(ysld+bb1*sd*dysd)
            df2d=dfnum2d*fxdend*fxdend
            dfd=df1d-df2d
            dfxdd=ax*roedth*x43*1.0*(fxd-dfd*sd)*fac2
            dfxdgd=ax*arodw*roedth*dfd*sddw*fac2
!     
!     ### end of potential contribution ###
!
         endif
!     _________________________________________________________________
!     correlation lyp(aroe,aroup,arodw,agr,agur,agdr)
!
         cden=1.0+dd*rometh
         cl1=-aa/cden
!
         omiga=exp(-cc*rometh)/cden/aroe**x113
         dilta=rometh*(cc+dd/cden)
         aroe2=aroe*aroe
         abo=aa*bb*omiga
!
         dodn=x13*omiga/aroe*(dilta-11.0)
         dddn=x13*(dd*dd*aroe**(-x53)/cden/cden-dilta/aroe)
!
         if(nspin.eq.1) then
!
            cl1=cl1*aroe
!
            cl21=4.0*cf*aroe**x83
            cl22=(x4718-x718*dilta)*agr2
            cl23=(2.5-x118*dilta)*agr2/2.0
            cl24=(dilta-11.0)/9.0*agr2/4.0
            cl25=x1124*agr2
!
            cl2=-abo*aroe2*(0.25*(cl21+cl22-cl23-cl24)-cl25)
!
!     ### potential contribution ###
!
            dl1dnu=-aa*(1/cden+x13*dd*rometh/cden/cden)
!
            dlt=x672+2.0*x772*dilta
            dl2dn=-abo*aroe*(cf*x143*aroe**x83-dlt*agr2)
            dl2do=cl2/omiga
            dl2dd=abo*aroe2*x772*agr2
            dl2dnu=dl2dn+dl2do*dodn+dl2dd*dddn
!     
            dl2dg=abo*aroe2*agr*dlt
!
!     ### end of potential contribution ###
!
         else
!
            cl11=cl1*4.0/aroe
            cl1=cl11*aroup*arodw
!
            aroup2=aroup*aroup
            arodw2=arodw*arodw
!
            cl21=t113*cf*(aroup**x83+arodw**x83)
            cl22=(x4718-x718*dilta)*agr2
            cl23=(2.5-x118*dilta)*(agur2+agdr2)
            dilta119=(dilta-11.0)/9.0
            cl24=dilta119/aroe*(aroup*agur2+arodw*agdr2)
            cl25=x23*aroe2*agr2
            cl26=(x23*aroe2-aroup2)*agdr2
            cl27=(x23*aroe2-arodw2)*agur2
!
            csum=cl21+cl22-cl23-cl24
            cl2=-abo*(aroup*arodw*csum-cl25+cl26+cl27)
!
!     ### potential contribution ###
!
!     *** cl1 has changed its form! ***
!
            dl1dn=cl1/aroe*(x13*dd/cden*rometh-1.0)
            dl1dnu=dl1dn+cl11*arodw
            dl1dnd=dl1dn+cl11*aroup
!     
            dl2dnu1=arodw*csum+                                         &
     &           arodw*aroup*(t113*cf*x83*aroup**x53-                   &
     &           dilta119*arodw/aroe2*(agur2-agdr2))-x43*aroe*agr2+     &
     &           x23*agdr2*(2.0*arodw-aroup)+x43*aroe*agur2
            dl2dnd1=aroup*csum+                                         &
     &           aroup*arodw*(t113*cf*x83*arodw**x53+                   &
     &           dilta119*aroup/aroe2*(agur2-agdr2))-x43*aroe*agr2+     &
     &           x23*agur2*(2.0*aroup-arodw)+x43*aroe*agdr2
!
            dl2do=cl2/omiga
            dl2dd=-abo*aroup*arodw*                                     &
     &           (-x718*agr2+x118*(agur2+agdr2)-                        &
     &           x19*(aroup*agur2+arodw*agdr2)/aroe)
!
            dl2dnu=-abo*dl2dnu1+dl2do*dodn+dl2dd*dddn
            dl2dnd=-abo*dl2dnd1+dl2do*dodn+dl2dd*dddn
!
            dl2dg=-abo*                                                 &
     &           (aroup*arodw*2.0*(x4718-x718*dilta)*agr-               &
     &           x43*aroe2*agr) 
            dl2dgu=-2.0*abo*agur*((x118*dilta-2.5-                      &
     &           dilta119*aroup/aroe)*aroup*arodw                       &
     &           +x23*aroe2-arodw2)
            dl2dgd=-2.0*abo*agdr*((x118*dilta-2.5-                      &
     &           dilta119*arodw/aroe)*aroup*arodw                       &
     &           +x23*aroe2-aroup2)
!
         endif
!
         clyp=cl1+cl2
!     _________________________________________________________________
!     updating of xc-energy
!
         excupdt=ex+clyp
!
         exc=exc+excupdt
!
!     _________________________________________________________________
!     first part xc-potential construction
!
!
         rhor(ir,isup)=dfxdu+(dl1dnu+dl2dnu)*fac1
         isign=sign(1.d0,agr-delt)
         byagr=0.5*(1+isign)/(agr+(1-isign)*delt)
!
         if(nspin.eq.1) then
!
            dexcdg=(dfxdg*aroe+dl2dg)*byagr
            gradr(ir,1,isup)=grx*dexcdg
            gradr(ir,2,isup)=gry*dexcdg
            gradr(ir,3,isup)=grz*dexcdg
!
         else
!
            rhor(ir,isdw)=dfxdd+(dl1dnd+dl2dnd)*fac2
!
            usign=sign(1.d0,agur-delt)
            dsign=sign(1.d0,agdr-delt)
            byagur=0.5*(1+usign)/(agur+(1-usign)*delt)
            byagdr=0.5*(1+dsign)/(agdr+(1-dsign)*delt)
!
            dexcdgu=(dfxdgu+dl2dgu)*byagur
            dexcdgd=(dfxdgd+dl2dgd)*byagdr
            dexcdg=dl2dg*byagr
!
            gradr(ir,1,isup)=(dexcdgu*grux+dexcdg*grx)*fac1
            gradr(ir,2,isup)=(dexcdgu*gruy+dexcdg*gry)*fac1
            gradr(ir,3,isup)=(dexcdgu*gruz+dexcdg*grz)*fac1
            gradr(ir,1,isdw)=(dexcdgd*grdx+dexcdg*grx)*fac2
            gradr(ir,2,isdw)=(dexcdgd*grdy+dexcdg*gry)*fac2
            gradr(ir,3,isdw)=(dexcdgd*grdz+dexcdg*grz)*fac2
!
         endif
!          
 100     continue
      end do
!
#ifdef PARA
      call reduce(1,exc)
#endif
      return
      end
!
!______________________________________________________________________
      subroutine ggapbe(nspin,rhog,gradr,rhor,excrho)
!     _________________________________________________________________
!     Perdew-Burke-Ernzerof gga
!     Perdew, et al. PRL 77, 3865, 1996
!
      use gvec
      use parm
      use cnst
!
      implicit none
! input
      integer nspin
      complex(kind=8) rhog(ng,nspin)
      real(kind=8)  gradr(nnr,3,nspin), rhor(nnr,nspin)
! output: excrho: exc * rho ;  E_xc = \int excrho(r) d_r
! output: rhor:   contains the exchange-correlation potential
      real(kind=8)  excrho
! local
      integer ir, icar, iss, isup, isdw, nspinx
      real(kind=8) lim1, lim2
      parameter ( lim1=1.d-8, lim2=1.d-8, nspinx=2 )
      real(kind=8) zet, arho(nspinx), grad(3,nspinx), agrad(nspinx),    &
     &     arhotot, gradtot(3), agradtot,                               &
     &     scl, scl1, wrkup, wrkdw,                                     &
     &     exrho(nspinx), dexdrho(nspinx), dexdg(nspinx),               &
     &     ecrho, decdrho(nspinx), decdg
!
!     main loop
!
      isup=1
      isdw=2
      do ir=1,nnr
!
         arho(isup) = abs(rhor(ir,isup))
         arhotot = arho(isup)
         zet = 0.d0
         do icar = 1, 3
            grad(icar,isup) = gradr(ir,icar,isup)
            gradtot(icar) = gradr(ir,icar,isup)
         enddo
!
         if (nspin.eq.2) then
            arho(isdw) = abs(rhor(ir,isdw))
            arhotot = abs(rhor(ir,isup)+rhor(ir,isdw))
            do icar = 1, 3
               grad(icar,isdw) = gradr(ir,icar,isdw)
               gradtot(icar) = gradr(ir,icar,isup)+gradr(ir,icar,isdw)
            enddo
            zet = (rhor(ir,isup) - rhor(ir,isdw)) / arhotot
            if (zet.ge. 1.d0) zet =  1.d0
            if (zet.le.-1.d0) zet = -1.d0
         endif
!
         do iss = 1, nspin
            agrad(iss) = sqrt( grad(1,iss)*grad(1,iss) +                &
     &                         grad(2,iss)*grad(2,iss) +                &
     &                         grad(3,iss)*grad(3,iss) )
            agradtot = sqrt( gradtot(1)*gradtot(1) +                    &
     &                       gradtot(2)*gradtot(2) +                    &
     &                       gradtot(3)*gradtot(3) )
         enddo
!
!     _________________________________________________________________
!     First it calculates the energy density excrho
!     exrho:  exchange term
!     ecrho:  correlation term
!
         if ( nspin.eq.2 ) then
            scl = 2.d0
            scl1 = 0.5d0
         else
            scl = 1.d0
            scl1 = 1.d0
         endif
         do iss = 1, nspin
            if ( arho(iss).gt.lim1) then
               call exchpbe( scl*arho(iss), scl*agrad(iss),             &
     &                       exrho(iss),dexdrho(iss),dexdg(iss))
               excrho = excrho + scl1*exrho(iss)
            else
               dexdrho(iss) = 0.d0
               dexdg(iss) = 0.d0
            endif
         enddo
         if ( arhotot.gt.lim1) then
            call ecorpbe( arhotot, agradtot, zet, ecrho,                &
     &                    decdrho(1), decdrho(2), decdg, nspin )
            excrho = excrho + ecrho
         else
            decdrho(isup) = 0.d0
            decdrho(isdw) = 0.d0
            decdg = 0.d0
         endif
!     _________________________________________________________________
!     Now it calculates the potential and writes it in rhor
!     it uses the following variables:
!     dexdrho = d ( ex*rho ) / d (rho)
!     decdrho = d ( ec*rho ) / d (rho)
!     dexdg   = (d ( ex*rho ) / d (grad(rho)_i)) * agrad / grad_i
!     decdg   = (d ( ec*rho ) / d (grad(rho)_i)) * agrad / grad_i
!     gradr  here is used as a working array
!
!     _________________________________________________________________
!    first part of the xc-potential : D(rho*exc)/D(rho)
!
         do iss = 1, nspin
            rhor(ir,iss) = dexdrho(iss) + decdrho(iss)
         enddo
!
!     gradr = D(rho*exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
         do iss = 1, nspin
            do icar = 1,3
               wrkup =0.d0
               wrkdw =0.d0
               if (agrad(iss).gt.lim2)                                  &
     &                    wrkup = dexdg(iss)*grad(icar,iss)/agrad(iss)
               if (agradtot.gt.lim2)                                    &
     &                    wrkdw = decdg*gradtot(icar)/agradtot
               gradr(ir,icar,iss) = wrkup + wrkdw
            enddo
         enddo
!
      end do
!
#ifdef PARA
      call reduce(1,excrho)
#endif
!
      return
      end
!
!______________________________________________________________________
      subroutine exchpbe(rho,agrad,ex,dexdrho,dexdg)
!     _________________________________________________________________
!
! Perdew-Burke-Ernzerof gga, Exchange term:
! Calculates the exchange energy density and the two functional derivative
! that will be used to calculate the potential
!
      implicit none
! input
! input rho:     charge density
! input agrad:   abs(grad rho)
      real(kind=8) rho, agrad
! ouput
! output ex: Ex[rho,grad_rho] = \int ex dr
! output dexdrho: d ex / d rho
! output dexdg:   d ex / d grad_rho(i) = dexdg*grad_rho(i)/abs(grad_rho)
      real(kind=8) ex, dexdrho, dexdg
! local
      real(kind=8) thrd, thrd4, pi, pi32td, ax, al, um, uk, ul
      parameter(thrd=.33333333333333333333d0,thrd4=4.d0/3.d0)
      parameter(pi=3.14159265358979323846264338327950d0)
      parameter(pi32td=3.09366772628014d0) ! pi32td=(3.d0*pi*pi)**0.333d0
      parameter(al=0.161620459673995d0)    ! al=1.0/(2.0*(pi32)**0.333d0)
      parameter(ax=-0.738558766382022405884230032680836d0)
      parameter(um=0.2195149727645171d0,uk=0.8040d0,ul=um/uk)
!
      real(kind=8) rhothrd, exunif, dexunif, kf, s, s2, p0, fxpbe, fs
!----------------------------------------------------------------------
! construct LDA exchange energy density
!
      rhothrd = rho**thrd
      dexunif = ax*rhothrd
      exunif  = rho*dexunif
!----------------------------------------------------------------------
! construct PBE enhancement factor
!
      kf = pi32td*rhothrd
      s = agrad/(2.d0*kf*rho)
      s2 = s*s
      p0 = 1.d0 + ul*s2
      fxpbe = 1.d0 + uk - uk/p0
      ex = exunif*fxpbe
!----------------------------------------------------------------------
! now calculates the potential terms
!
!  fs=(1/s)*d fxPBE/ ds
!
      fs=2.d0*uk*ul/(p0*p0)
      dexdrho = dexunif*thrd4*(fxpbe-s2*fs)
      dexdg = ax*al*s*fs
!
      return
      end
!----------------------------------------------------------------------
      subroutine ecorpbe(rho,agrad,zet,ectot,decup,decdn,decdg,nspin)
!     -----------------------------------------------------------------
!
!  Adapted from the Official PBE correlation code. K. Burke, May 14, 1996.
!
!   input: rho   = rho_up + rho_down; total  charge density
!   input: agrad = abs( grad(rho) )
!   input: zet   = (rho_up-rho_down)/rho
!   input: nspin
!  output: ectot = ec*rho       ---correlation energy density---
!  output: decup = d ( ec*rho ) / d (rho_up)
!  output: decdn = d ( ec*rho ) / d (rho_down)
!  output: decdg = (d ( ec*rho ) / d (grad(rho)_i)) * agrad / grad_i
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! References:
! [a] J.P.~Perdew, K.~Burke, and M.~Ernzerhof, 
!     {\sl Generalized gradient approximation made simple}, sub.
!     to Phys. Rev.Lett. May 1996.
! [b] J. P. Perdew, K. Burke, and Y. Wang, {\sl Real-space cutoff
!     construction of a generalized gradient approximation:  The PW91
!     density functional}, submitted to Phys. Rev. B, Feb. 1996.
! [c] J. P. Perdew and Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).
!----------------------------------------------------------------------
!----------------------------------------------------------------------
      implicit none
      real(kind=8) rho, agrad, zet, ectot, decup, decdn, decdg
      integer nspin
      real(kind=8) pi, pi32, alpha, thrd, thrdm, thrd2, sixthm, thrd4,  &
     &     gam, fzz, gamma, bet, delt, eta
! thrd*=various multiples of 1/3
! numbers for use in LSD energy spin-interpolation formula, [c](9).
!      gam= 2^(4/3)-2
!      fzz=f''(0)= 8/(9*gam)
! numbers for construction of PBE
!      gamma=(1-log(2))/pi^2
!      bet=coefficient in gradient expansion for correlation, [a](4).
!      eta=small number to stop d phi/ dzeta from blowing up at 
!          |zeta|=1.
      parameter(pi=3.1415926535897932384626433832795d0)
      parameter(pi32=29.608813203268075856503472999628d0)
      parameter(alpha=1.91915829267751300662482032624669d0)
      parameter(thrd=1.d0/3.d0,thrdm=-thrd,thrd2=2.d0*thrd)
      parameter(sixthm=thrdm/2.d0)
      parameter(thrd4=4.d0*thrd)
      parameter(gam=0.5198420997897463295344212145565d0)
      parameter(fzz=8.d0/(9.d0*gam))
      parameter(gamma=0.03109069086965489503494086371273d0)
      parameter(bet=0.06672455060314922d0,delt=bet/gamma)
      parameter(eta=1.d-12)
      real(kind=8) g, fk, rs, sk, twoksg, t
      real(kind=8) rtrs, eu, eurs, ep, eprs, alfm, alfrsm, z4, f, ec
      real(kind=8) ecrs, fz, eczet, comm, vcup, vcdn, g3, pon, b, b2, t2, t4
      real(kind=8) rs2, rs3, q4, q5, h, g4, t6, rsthrd, gz, fac
      real(kind=8) bg, bec, q8, q9, hb, hrs, hz, ht, pref
!----------------------------------------------------------------------
      if (nspin.eq.1) then
         g=1.d0
      else
         g=((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)*0.5d0
      endif
      fk=(pi32*rho)**thrd
      rs=alpha/fk
      sk=sqrt(4.d0*fk/pi)
      twoksg=2.d0*sk*g
      t=agrad/(twoksg*rho)
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! find LSD energy contributions, using [c](10) and Table I[c].
! eu=unpolarized LSD correlation energy
! eurs=deu/drs
! ep=fully polarized LSD correlation energy
! eprs=dep/drs
! alfm=-spin stiffness, [c](3).
! alfrsm=-dalpha/drs
! f=spin-scaling factor from [c](9).
! construct ec, using [c](8)
      rtrs=dsqrt(rs)
      call gcor2(0.0310907d0,0.21370d0,7.5957d0,3.5876d0,1.6382d0,      &
     &    0.49294d0,rtrs,eu,eurs)
      if (nspin.eq.2) then
         call gcor2(0.01554535d0,0.20548d0,14.1189d0,6.1977d0,3.3662d0, &
     &       0.62517d0,rtrs,ep,eprs)
         call gcor2(0.0168869d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,  &
     &       0.49671d0,rtrs,alfm,alfrsm)
         z4 = zet**4
         f=((1.d0+zet)**thrd4+(1.d0-zet)**thrd4-2.d0)/gam
         ec = eu*(1.d0-f*z4)+ep*f*z4-alfm*f*(1.d0-z4)/fzz
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! LSD potential from [c](A1)
! ecrs = dec/drs [c](A2)
! eczet=dec/dzeta [c](A3)
! fz = df/dzeta [c](A4)
         ecrs = eurs*(1.d0-f*z4)+eprs*f*z4-alfrsm*f*(1.d0-z4)/fzz
         fz = thrd4*((1.d0+zet)**thrd-(1.d0-zet)**thrd)/gam
         eczet = 4.d0*(zet**3)*f*(ep-eu+alfm/fzz)+fz*(z4*ep-z4*eu       &
     &           -(1.d0-z4)*alfm/fzz)
         comm = ec -rs*ecrs/3.d0-zet*eczet
         vcup = comm + eczet
         vcdn = comm - eczet
      else
         ecrs = eurs
         ec = eu
         vcup = ec -rs*ecrs/3.d0
      endif
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! PBE correlation energy
! g=phi(zeta), given after [a](3)
! delt=bet/gamma
! b=a of [a](8)
!      g=((1.d0+zet)**thrd2+(1.d0-zet)**thrd2)/2.d0
      g3 = g**3
      pon=-ec/(g3*gamma)
      b = delt/(dexp(pon)-1.d0)
      b2 = b*b
      t2 = t*t
      t4 = t2*t2
      rs2 = rs*rs
      rs3 = rs2*rs
      q4 = 1.d0+b*t2
      q5 = 1.d0+b*t2+b2*t4
      h = g3*(bet/delt)*dlog(1.d0+delt*Q4*t2/Q5)
      ectot = rho*(ec + h)
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! energy done. Now the potential, using appendix e of [b].
      t6 = t4*t2
      rsthrd = rs/3.d0
      fac = delt/b+1.d0
      bec = b2*fac/(bet*g3)
      q8 = q5*q5+delt*q4*q5*t2
      q9 = 1.d0+2.d0*b*t2
      hb = -bet*g3*b*t6*(2.d0+b*t2)/q8
      hrs = -rsthrd*hb*bec*ecrs
      ht = 2.d0*bet*g3*q9/q8
      comm = h+hrs-7.d0*t2*ht/6.d0
      if (nspin.eq.2) then
         g4 = g3*g
         bg = -3.d0*b2*ec*fac/(bet*g4)
         gz=(((1.d0+zet)**2+eta)**sixthm-                               &
     &   ((1.d0-zet)**2+eta)**sixthm)/3.d0
         hz = 3.d0*gz*h/g + hb*(bg*gz+bec*eczet)
         pref = hz-gz*t2*ht/g
         decup = vcup + comm + pref*(  1.d0 - zet)
         decdn = vcdn + comm + pref*( -1.d0 - zet)
      else
         decup = vcup + comm
      endif
      decdg = t*ht/twoksg
!
      return
      end
!______________________________________________________________________
      subroutine gcor2(a,a1,b1,b2,b3,b4,rtrs,gg,ggrs)
!     _________________________________________________________________
! slimmed down version of GCOR used in PW91 routines, to interpolate
! LSD correlation energy, as given by (10) of
! J. P. Perdew and Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).
! K. Burke, May 11, 1996.
!
      implicit none
      real(kind=8) a, a1, b1, b2, b3, b4, rtrs, gg, ggrs
      real(kind=8) q0, q1, q2, q3
!
      q0 = -2.d0*a*(1.d0+a1*rtrs*rtrs)
      q1 = 2.d0*a*rtrs*(b1+rtrs*(b2+rtrs*(b3+b4*rtrs)))
      q2 = dlog(1.d0+1.d0/q1)
      gg = q0*q2
      q3 = a*(b1/rtrs+2.d0*b2+rtrs*(3.d0*b3+4.d0*b4*rtrs))
      ggrs = -2.d0*a*a1*q2-q0*q3/(q1*(1.d0+q1))
!
      return
      end
!
      subroutine zero(n,a)
!
      integer n,i
      real(kind=8) a(n)
!
      do i=1,n
         a(i)=0.0
      end do
!
      return
      end
!______________________________________________________________________
      subroutine ggapw(nspin,rhog,gradr,rhor,exc)
!     _________________________________________________________________
!     perdew-wang gga (PW91)
!
      use gvec
      use parm
      use cnst
!
      implicit none
! input
      integer nspin
      complex(kind=8) rhog(ng,nspin)
      real(kind=8) gradr(nnr,3,nspin), rhor(nnr,nspin)
! output
      real(kind=8) exc
! local
      integer isup, isdw, ir
      real(kind=8) rhoup, rhodw, roe, aroup, arodw, aroe, rs, zeta
      real(kind=8) grxu, gryu, grzu, grhou, grxd, gryd, grzd, grhod, grho
      real(kind=8) ex, vx, ec,vc, sc, sx, v1x, v2x, v1c, v2c
      real(kind=8) ecrs, eczeta
      real(kind=8) exup, vcup, sxup, v1xup, v2xup, v1cup, v2cup
      real(kind=8) exdw, vcdw, sxdw, v1xdw, v2xdw, v1cdw, v2cdw
      real(kind=8) pi34, third, small
      parameter (pi34=0.75d0/3.141592653589793d+00,third=1.d0/3.d0)
      parameter(small=1.d-10)
!
!     _________________________________________________________________
!     main loop
!
      isup=1
      isdw=2
      exc=0.0
      do ir=1,nnr
         rhoup=rhor(ir,isup)
         if(nspin.eq.2) then
            rhodw=rhor(ir,isdw)
         else
            rhodw=0.0
         end if
         roe=rhoup+rhodw
         aroe=abs(roe)
         if (aroe.lt.small) then
            rhor(ir,isup)   =0.0
            gradr(ir,1,isup)=0.0
            gradr(ir,2,isup)=0.0
            gradr(ir,3,isup)=0.0
            if(nspin.eq.2) then
               rhor(ir,isdw)   =0.0
               gradr(ir,1,isdw)=0.0
               gradr(ir,2,isdw)=0.0
               gradr(ir,3,isdw)=0.0
            end if
            go to 100
         end if
         grxu =gradr(ir,1,isup)
         gryu =gradr(ir,2,isup)
         grzu =gradr(ir,3,isup)
         grhou=sqrt(grxu**2+gryu**2+grzu**2)
         if(nspin.eq.2) then
            grxd =gradr(ir,1,isdw)
            gryd =gradr(ir,2,isdw)
            grzd =gradr(ir,3,isdw)
            grhod=sqrt(grxd**2+gryd**2+grzd**2)
         else
            grxd =0.0
            gryd =0.0
            grzd =0.0
            grhod=0.0
         endif
         grho=sqrt((grxu+grxd)**2+(gryu+gryd)**2+(grzu+grzd)**2)
!
         rs=(pi34/aroe)**third
         if (nspin.eq.1) then
            call exchpw91(aroe,grho,ex,v1x,v2x)
            call pwlda(rs,ec,vc,ecrs)
            call corpw91ns(rs,grho,ec,ecrs,sc,v1c,v2c)
            exc = exc + roe*(ex+ec) + sc
            rhor(ir,isup) = vc + v1x + v1c
!
!     gradr = D(rho*exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
            gradr(ir,1,isup)=grxu*(v2x+v2c)
            gradr(ir,2,isup)=gryu*(v2x+v2c)
            gradr(ir,3,isup)=grzu*(v2x+v2c)
         else
            zeta=(rhoup-rhodw)/aroe        
            zeta=min(zeta, 1.d0)
            zeta=max(zeta,-1.d0)
            call exchpw91(2.d0*abs(rhoup),2.0*grhou,exup,v1xup,v2xup)
            call exchpw91(2.d0*abs(rhodw),2.0*grhod,exdw,v1xdw,v2xdw)
            call pwlsd(rs,zeta,ec,vcup,vcdw,ecrs,eczeta)
            call corpw91(rs,zeta,grho,ec,ecrs,eczeta,sc,v1cup,v1cdw,v2c)
            rhor(ir,isup) = vcup + v1xup + v1cup
            rhor(ir,isdw) = vcdw + v1xdw + v1cdw
            exc = exc+roe*(0.5*((1.d0+zeta)*exup+(1.d0-zeta)*exdw)+ec) &
    &             + sc
!
!     gradr = D(rho*exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
            gradr(ir,1,isup)=grxu*(2.0*v2xup+v2c)+grxd*v2c
            gradr(ir,2,isup)=gryu*(2.0*v2xup+v2c)+gryd*v2c
            gradr(ir,3,isup)=grzu*(2.0*v2xup+v2c)+grzd*v2c
            gradr(ir,1,isdw)=grxd*(2.0*v2xdw+v2c)+grxu*v2c
            gradr(ir,2,isdw)=gryd*(2.0*v2xdw+v2c)+gryu*v2c
            gradr(ir,3,isdw)=grzd*(2.0*v2xdw+v2c)+grzu*v2c
         end if
 100     continue
      end do
!
#ifdef PARA
      call reduce(1,exc)
#endif
!
      return
      end
!
!----------------------------------------------------------------------
      subroutine exchpw91(rho,grho,ex,v1x,v2x)
!----------------------------------------------------------------------
!
!  PW91 exchange for a spin-unpolarized electronic system
!  Modified from the "official" PBE code of Perdew, Burke et al.
!  input rho   : density
!  input grho:  abs(grad rho)
!  output:  exchange energy per electron (ex) and potentials
!          v1x = d(rho*exc)/drho
!          v2x = d(rho*exc)/d|grho| * (1/|grho|)
!
      implicit none
!  input
      real(kind=8) rho, grho
!  output
      real(kind=8) ex, v1x, v2x
! local
      real(kind=8) ex0, kf, s, s2, s3, s4, f, fs,                      &
     &       p0, p1, p2, p3, p4, p5, p6, p7
! parameters
      real(kind=8) a1, a2, a3, a4, a, b1, bx, pi34, thrd, thrd4
      parameter(a1=0.19645d0,a2=0.27430d0,a=7.7956d0,a4=100.d0)
! for becke exchange, set a3=b1=0
      parameter(a3=0.15084d0,b1=0.004d0)
! pi34=3/(4pi) ,  bx=(3pi^2)^(1/3)
      parameter(pi34=0.75d0/3.141592653589793d+00, bx=3.09366773d0,    &
     &          thrd=0.333333333333d0, thrd4=4.d0*thrd)
!
      if (rho.lt.1.d-10) then
         ex =0.0
         v1x=0.0
         v2x=0.0
      end if
!
!  kf=k_Fermi, ex0=Slater exchange energy
!
      kf = bx*(rho**thrd)
      ex0=-pi34*kf
      if (grho.lt.1.d-10) then
         ex =ex0
         v1x=ex0*thrd4
         v2x=0.0
      end if
      s  = grho/(2.d0*kf*rho)
      s2 = s*s
      s3 = s2*s
      s4 = s2*s2
      p0 = 1.d0/sqrt(1.d0+a*a*s2)
      p1 = log(a*s+1.d0/p0)
      p2 = exp(-a4*s2)
      p3 = 1.d0/(1.d0+a1*s*p1+b1*s4)
      p4 = 1.d0+a1*s*p1+(a2-a3*p2)*s2
!  f is the enhancement factor
      f = p3*p4
      ex = ex0*f
!  energy done. now the potential:
      p5 = b1*s2-(a2-a3*p2)
      p6 = a1*s*(p1+a*s*p0)
      p7 = 2.d0*(a2-a3*p2)+2.d0*a3*a4*s2*p2-4.d0*b1*s2*f
! fs = (1/s) dF(s)/ds
      fs = p3*(p3*p5*p6+p7)
      v1x = ex0*thrd4*(f-s2*fs)
      v2x = 0.5d0*ex0/kf*s*fs/grho
!
      return
      end
!
!----------------------------------------------------------------------
      subroutine corpw91ns(rs,grho,ec,ecrs,h,v1c,v2c)
!----------------------------------------------------------------------
!
!  PW91 correlation (gradient correction term) - no spin case
!  Modified from the "official" PBE code of Perdew, Burke et al.
!
!  input rs:   seitz radius
!  input zeta:  relative spin polarization
!  input grho: abs(grad rho)
!  input ec:   Perdew-Wang correlation energy
!  input ecrs:  d(rho*ec)/d r_s

!  output h  :  nonlocal part of correlation energy per electron
!  output v1c:  nonlocal parts of correlation potential
!         v1c = d(rho*exc)/drho
!         v2c = d(rho*exc)/d|grho|*(1/|grho|)
!
      implicit none
! input
      real(kind=8) rs, grho, ec, ecrs
! output
      real(kind=8) h, v1c, v2c
! local
      real(kind=8) rho, t, ks,  bet, delt, pon, b, b2, t2, t4, t6
      real(kind=8) q4, q5, q6, q7, q8, q9, r0, r1, r2, r3, r4, rs2, rs3
      real(kind=8) vcdn, ccrs, rsthrd, fac, bg, bec, coeff, cc
      real(kind=8) h0, h0b, h0rs, h0t, h1, h1t, h1rs, hrs, ht
! parameters
      real(kind=8) nu, cc0, cx, alf, c1, c2, c3, c4, c5, c6, c7, a4
      real(kind=8) ax, pi34
      parameter(nu=15.75592d0,cc0=0.004235d0,cx=-0.001667212d0)
      parameter(c1=0.002568d0,c2=0.023266d0,c3=7.389d-6,c4=8.723d0)
      parameter(c5=0.472d0,c6=7.389d-2,a4=100.d0, alf=0.09d0)
! ax=(4*1.9191583/pi)^(1/2), where k_F=1.9191583/r_s, k_s=boh*r_s^(1/2)
      parameter(ax=1.5631853d0, pi34=0.75d0/3.141592653589793d0)
!
!
      rs2 = rs*rs
      rs3 = rs2*rs
      rho=pi34/rs3
!  k_s=(4k_F/pi)^(1/2)
      ks=ax/sqrt(rs)
!  t=abs(grad rho)/(rho*2.*ks)
      t=grho/(2.d0*rho*ks)
      bet = nu*cc0
      delt = 2.d0*alf/bet
      pon = -delt*ec/bet
      b = delt/(exp(pon)-1.d0)
      b2 = b*b
      t2 = t*t
      t4 = t2*t2
      t6 = t4*t2
      q4 = 1.d0+b*t2
      q5 = 1.d0+b*t2+b2*t4
      q6 = c1+c2*rs+c3*rs2
      q7 = 1.d0+c4*rs+c5*rs2+c6*rs3
      cc = -cx + q6/q7
      r0 = 0.663436444d0*rs
      r1 = a4*r0
      coeff = cc-cc0-3.d0*cx/7.d0
      r2 = nu*coeff
      r3 = exp(-r1*t2)
      h0 = (bet/delt)*log(1.d0+delt*q4*t2/q5)
      h1 = r3*r2*t2
      h = (h0+h1)*rho
!  energy done. now the potential:
      ccrs = (c2+2.*c3*rs)/q7 - q6*(c4+2.*c5*rs+3.*c6*rs2)/q7**2
      rsthrd = rs/3.d0
      r4 = rsthrd*ccrs/coeff
      fac = delt/b+1.d0
      bg = -3.d0*b2*ec*fac/bet
      bec = b2*fac/bet
      q8 = q5*q5+delt*q4*q5*t2
      q9 = 1.d0+2.d0*b*t2
      h0b = -bet*b*t6*(2.d0+b*t2)/q8
      h0rs = -rsthrd*h0b*bec*ecrs
      h0t = 2.*bet*q9/q8
      h1rs = r3*r2*t2*(-r4+r1*t2/3.d0)
      h1t = 2.d0*r3*r2*(1.d0-r1*t2)
      hrs = h0rs+h1rs
      ht = h0t+h1t
      v1c = h0+h1+hrs-7.d0*t2*ht/6.d0
      v2c = t*ht/(2.d0*ks*grho)
!
      return
      end
!
!----------------------------------------------------------------------
      subroutine corpw91(rs,zeta,grho,ec,ecrs,eczeta,h,v1cup,v1cdn,v2c)
!----------------------------------------------------------------------
!
!  PW91 correlation (gradient correction term)
!  Modified from the "official" PBE code of Perdew, Burke et al.
!
!  input rs:   seitz radius
!  input zeta:  relative spin polarization
!  input grho: abs(grad rho)
!  input ec:   Perdew-Wang correlation energy
!  input ecrs:  d(rho*ec)/d r_s ?
!  input eczeta: d(rho*ec)/d zeta ?

!  output h: nonlocal part of correlation energy per electron
!  output v1cup,v1cdn:  nonlocal parts of correlation potentials
!         v1c** = d(rho*exc)/drho               (up and down components)
!         v2c   = d(rho*exc)/d|grho|*(1/|grho|) (same for up and down)
!
      implicit none
! input
      real(kind=8) rs, zeta, grho, ec, ecrs, eczeta
! output
      real(kind=8) h, v1cup, v1cdn, v2c
! local
      real(kind=8) rho, g, t, ks, gz, bet, delt, g3, g4, pon, b, b2, t2, t4, t6
      real(kind=8) q4, q5, q6, q7, q8, q9, r0, r1, r2, r3, r4, rs2, rs3
      real(kind=8) vcdn, ccrs, rsthrd, fac, bg, bec, coeff, cc
      real(kind=8) h0, h0b, h0rs, h0z, h0t, h1, h1t, h1rs, h1z
      real(kind=8) hz, hrs, ht, comm, pref
! parameters
      real(kind=8) nu, cc0, cx, alf, c1, c2, c3, c4, c5, c6, c7, a4
      real(kind=8) thrdm, thrd2, ax, pi34, eta
      parameter(nu=15.75592d0,cc0=0.004235d0,cx=-0.001667212d0)
      parameter(c1=0.002568d0,c2=0.023266d0,c3=7.389d-6,c4=8.723d0)
      parameter(c5=0.472d0,c6=7.389d-2,a4=100.d0, alf=0.09d0)
      parameter(thrdm=-0.333333333333d0,thrd2=0.666666666667d0)
! ax=(4*1.9191583/pi)^(1/2), where k_F=1.9191583/r_s, k_s=boh*r_s^(1/2)
      parameter(ax=1.5631853d0, pi34=0.75d0/3.141592653589793d0)
      parameter(eta=1.d-12)
!
!
      if (grho.lt.1.d-10) then
         h=0.0
         v1cup=0.0
         v1cdn=0.0
         v2c=0.0
      end if
      rs2 = rs*rs
      rs3 = rs2*rs
      rho=pi34/rs3
      g=((1.d0+zeta)**thrd2+(1.d0-zeta)**thrd2)/2.d0
!  k_s=(4k_F/pi)^(1/2)
      ks=ax/sqrt(rs)
!  t=abs(grad rho)/(rho*2.*ks*g)
      t=grho/(2.d0*rho*g*ks)
      bet = nu*cc0
      delt = 2.d0*alf/bet
      g3 = g**3
      g4 = g3*g
      pon = -delt*ec/(g3*bet)
      b = delt/(exp(pon)-1.d0)
      b2 = b*b
      t2 = t*t
      t4 = t2*t2
      t6 = t4*t2
      q4 = 1.d0+b*t2
      q5 = 1.d0+b*t2+b2*t4
      q6 = c1+c2*rs+c3*rs2
      q7 = 1.d0+c4*rs+c5*rs2+c6*rs3
      cc = -cx + q6/q7
      r0 = 0.663436444d0*rs
      r1 = a4*r0*g4
      coeff = cc-cc0-3.d0*cx/7.d0
      r2 = nu*coeff*g3
      r3 = dexp(-r1*t2)
      h0 = g3*(bet/delt)*log(1.d0+delt*q4*t2/q5)
      h1 = r3*r2*t2
      h = (h0+h1)*rho
!  energy done. now the potential:
      ccrs = (c2+2.*c3*rs)/q7 - q6*(c4+2.*c5*rs+3.*c6*rs2)/q7**2
      rsthrd = rs/3.d0
      r4 = rsthrd*ccrs/coeff
!  eta is a small quantity that avoids trouble if zeta=+1 or -1
      gz = ((1.d0+zeta+eta)**thrdm - (1.d0-zeta+eta)**thrdm)/3.d0
      fac = delt/b+1.d0
      bg = -3.d0*b2*ec*fac/(bet*g4)
      bec = b2*fac/(bet*g3)
      q8 = q5*q5+delt*q4*q5*t2
      q9 = 1.d0+2.d0*b*t2
      h0b = -bet*g3*b*t6*(2.d0+b*t2)/q8
      h0rs = -rsthrd*h0b*bec*ecrs
      h0z = 3.d0*gz*h0/g + h0b*(bg*gz+bec*eczeta)
      h0t = 2.*bet*g3*q9/q8
      h1rs = r3*r2*t2*(-r4+r1*t2/3.d0)
      h1z = gz*r3*r2*t2*(3.d0-4.d0*r1*t2)/g
      h1t = 2.d0*r3*r2*(1.d0-r1*t2)
      hrs = h0rs+h1rs
      ht = h0t+h1t
      hz = h0z+h1z
      comm = h0+h1+hrs-7.d0*t2*ht/6.d0
      pref = hz-gz*t2*ht/g
      comm = comm-pref*zeta
      v1cup = comm + pref 
      v1cdn = comm - pref
      v2c   = t*ht/(2.d0*ks*g*grho)
!
      return
      end
!----------------------------------------------------------------------
      subroutine pwlda(rs,ec,vc,ecrs)
!----------------------------------------------------------------------
!
!  uniform-gas, spin-unpolarised correlation of perdew and wang 1991
!  input:  rs   seitz radius
!  output: ec   correlation energy per electron
!          vc   potential
!          ecrs derivatives of ec wrt rs
!
      implicit none
! input
      real(kind=8) rs
! output
      real(kind=8) ec, vc, ecrs
!  local
      real(kind=8) q0, rs12, q1, q2, q3
! parameters
      real(kind=8) a, a1, b1, b2, b3, b4
      parameter(a =0.0310907d0, a1=0.21370d0, b1=7.5957d0,              &
     &          b2=3.5876d0,    b3=1.6382d0,  b4=0.49294d0)
!
      q0 = -2.d0*a*(1.d0+a1*rs)
      rs12 = sqrt(rs)
      q1 = 2.d0*a*rs12*(b1+rs12*(b2+rs12*(b3+b4*rs12)))
      q2 = log(1.d0+1.d0/q1)
      ec = q0*q2
      q3 = a*(b1/rs12+2.d0*b2+3.d0*b3*rs12+2.d0*b4*2.d0*rs)
      ecrs = -2.d0*a*a1*q2-q0*q3/(q1**2+q1)
      vc = ec - rs*ecrs/3.d0
!
      return
      end
!----------------------------------------------------------------------
      subroutine pwlsd(rs,zeta,ec,vcup,vcdn,ecrs,eczeta)
!----------------------------------------------------------------------
!
!  uniform-gas correlation of perdew and wang 1991
!  Modified from the "official" PBE code of Perdew, Burke et al.
!  input: seitz radius (rs), relative spin polarization (zeta)
!  output: correlation energy per electron (ec)
!          up- and down-spin potentials (vcup,vcdn)
!          derivatives of ec wrt rs (ecrs) & zeta (eczeta)
!
      implicit none
! input
      real(kind=8) rs, zeta
! output
      real(kind=8) ec, vcup, vcdn, ecrs, eczeta
! local
      real(kind=8) f, eu, ep, eurs, eprs, alfm, alfrsm, z4, fz, comm
      real(kind=8) rs12, q0, q1, q2, q3
! parameters
      real(kind=8) gam, fzz, thrd, thrd4
      parameter(gam=0.5198421d0,fzz=1.709921d0)
      parameter(thrd=0.333333333333d0,thrd4=1.333333333333d0)
!
      real(kind=8) au, au1, bu1, bu2, bu3, bu4
      parameter(au =0.0310907d0, au1=0.21370d0, bu1=7.5957d0,           &
     &          bu2=3.5876d0,    bu3=1.6382d0,  bu4=0.49294d0)
      real(kind=8) ap, ap1, bp1, bp2, bp3, bp4
      parameter(ap =0.01554535d0,ap1=0.20548d0, bp1=14.1189d0,          &
     &          bp2=6.1977d0,    bp3=3.3662d0,  bp4=0.62517d0 )
      real(kind=8) am, am1, bm1, bm2, bm3, bm4
      parameter(am =0.0168869d0, am1=0.11125d0, bm1=10.357d0,           &
     &          bm2=3.6231d0,    bm3=0.88026d0, bm4=0.49671d0 )
!
      rs12 = sqrt(rs)
!
      q0 = -2.d0*au*(1.d0+au1*rs)
      q1 = 2.d0*au*rs12*(bu1+rs12*(bu2+rs12*(bu3+bu4*rs12)))
      q2 = log(1.d0+1.d0/q1)
      eu = q0*q2
      q3 = au*(bu1/rs12+2.d0*bu2+3.d0*bu3*rs12+2.d0*bu4*2.d0*rs)
      eurs = -2.d0*au*au1*q2-q0*q3/(q1**2+q1)
!
      q0 = -2.d0*ap*(1.d0+ap1*rs)
      q1 = 2.d0*ap*rs12*(bp1+rs12*(bp2+rs12*(bp3+bp4*rs12)))
      q2 = log(1.d0+1.d0/q1)
      ep = q0*q2
      q3 = ap*(bp1/rs12+2.d0*bp2+3.d0*bp3*rs12+2.d0*bp4*2.d0*rs)
      eprs = -2.d0*ap*ap1*q2-q0*q3/(q1**2+q1)
!
      q0 = -2.d0*am*(1.d0+am1*rs)
      q1 = 2.d0*am*rs12*(bm1+rs12*(bm2+rs12*(bm3+bm4*rs12)))
      q2 = log(1.d0+1.d0/q1)
!  alfm is minus the spin stiffness alfc
      alfm=q0*q2
      q3 = am*(bm1/rs12+2.d0*bm2+3.d0*bm3*rs12+2.d0*bm4*2.d0*rs)
      alfrsm=-2.d0*am*am1*q2-q0*q3/(q1**2+q1)
!
      f = ((1.d0+zeta)**thrd4+(1.d0-zeta)**thrd4-2.d0)/gam
      z4 = zeta**4
      ec = eu*(1.d0-f*z4)+ep*f*z4-alfm*f*(1.d0-z4)/fzz
!  energy done. now the potential:
      ecrs = eurs*(1.d0-f*z4)+eprs*f*z4-alfrsm*f*(1.d0-z4)/fzz
      fz = thrd4*((1.d0+zeta)**thrd-(1.d0-zeta)**thrd)/gam
      eczeta = 4.d0*(zeta**3)*f*(ep-eu+alfm/fzz)+fz*(z4*ep-z4*eu        &
     &        -(1.d0-z4)*alfm/fzz)
      comm = ec -rs*ecrs/3.d0-zeta*eczeta
      vcup = comm + eczeta
      vcdn = comm - eczeta
!
      return
      end
!
!______________________________________________________________________
      subroutine fillgrad(nspin,rhog,gradr)
!     _________________________________________________________________
!
!     calculates gradient of charge density for gradient corrections
!     in: charge density on G-space    out: gradient in R-space
!
      use gvec
      use parm
      use work1
!
      implicit none
! input
      integer nspin
      complex(kind=8) rhog(ng,nspin)
! output
      real(kind=8)    gradr(nnr,3,nspin)
! local
      complex(kind=8), pointer:: v(:)
      complex(kind=8) ci
      integer iss, ig, ir
!
!
      v => wrk1
      ci=(0.0,1.0)
      do iss=1,nspin
         do ig=1,nnr
            v(ig)=(0.0,0.0)
         end do
         do ig=1,ng
            v(np(ig))=      ci*tpiba*gx(ig,1)*rhog(ig,iss)
            v(nm(ig))=conjg(ci*tpiba*gx(ig,1)*rhog(ig,iss))
         end do
         call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
         do ir=1,nnr
            gradr(ir,1,iss)=real(v(ir))
         end do
!
         do ig=1,nnr
            v(ig)=(0.0,0.0)
         end do
         do ig=1,ng
            v(np(ig))= tpiba*(      ci*gx(ig,2)*rhog(ig,iss)-           &
     &                                 gx(ig,3)*rhog(ig,iss) )
            v(nm(ig))= tpiba*(conjg(ci*gx(ig,2)*rhog(ig,iss)+           &
     &                                 gx(ig,3)*rhog(ig,iss)))
         end do
         call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
         do ir=1,nnr
            gradr(ir,2,iss)= real(v(ir))
            gradr(ir,3,iss)=aimag(v(ir))
         end do
      end do
!
      return
      end
!______________________________________________________________________
      subroutine grad2(nspin,gradr,rhor)
!     _________________________________________________________________
!
!     calculate the second part of gradient corrected xc potential
!     \sum_alpha (D / D r_alpha) ( D(rho*exc)/D(grad_alpha rho) )
!
      use gvec
      use parm
      use work1
!
      implicit none
! input
      integer nspin
      real(kind=8) gradr(nnr,3,nspin)
! input/output
      real(kind=8) rhor(nnr,nspin)
! local
      complex(kind=8), pointer:: v(:)
      complex(kind=8), allocatable:: x(:)
      complex(kind=8) ci, fp, fm
      integer iss, ig, ir
!
      v => wrk1
      allocate(x(ng))
      ci=(0.0,1.0)
      do iss=1, nspin
!
! x polarization
!
!     copy input gradr(r) into a complex array, v(r)...
!
         do ir=1,nnr
            v(ir)=cmplx(gradr(ir,1,iss),0.0)
         end do
!
!     bring v(r) to G-space, v(G)...
!
         call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!
!     multiply by (iG) to get x=(\grad_x gradr)(G)
!
         do ig=1,ng
            x(ig)=ci*tpiba*gx(ig,1)*v(np(ig))
         end do
!
! y and z polarizations: as above, two fft's together
!
         do ir=1,nnr
            v(ir)=cmplx(gradr(ir,2,iss),gradr(ir,3,iss))
         end do
         call fwfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
         do ig=1,ng
            fp=v(np(ig))+v(nm(ig))
            fm=v(np(ig))-v(nm(ig))
            x(ig) = x(ig) +                                             &
     &           ci*tpiba*gx(ig,2)*0.5*cmplx( real(fp),aimag(fm))
            x(ig) = x(ig) +                                             &
     &           ci*tpiba*gx(ig,3)*0.5*cmplx(aimag(fp),-real(fm))
         end do
!
!     x = \sum_alpha(\grad_alpha gradr)(G)
!     now bring back to R-space
!
         do ig=1,nnr
            v(ig)=(0.0,0.0)
         end do
         do ig=1,ng
            v(np(ig))=x(ig)
            v(nm(ig))=conjg(x(ig))
         end do
         call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)
!
!     v = \sum_alpha(\grad_x gradr)(r)
!
         do ir=1,nnr
            rhor(ir,iss)=rhor(ir,iss)-real(v(ir))
         end do
      end do
!
      deallocate(x)
!
      return
      end
!----------------------------------------------------------------------
      subroutine checkrho(nnr,nspin,rhor,rmin,rmax,rsum,rnegsum)
!----------------------------------------------------------------------
!
!     check \int rho(r)dr and the negative part of rho
!
      implicit none
      integer nnr, nspin
      real(kind=8) rhor(nnr,nspin), rmin, rmax, rsum, rnegsum
!
      real(kind=8) roe
      integer ir, iss
!
      rsum   =0.0
      rnegsum=0.0
      rmin   =100.
      rmax   =0.0
      do iss=1,nspin
         do ir=1,nnr
            roe=rhor(ir,iss)
            rsum=rsum+roe
            if (roe.lt.0.0) rnegsum=rnegsum+roe
            rmax=max(rmax,roe)
            rmin=min(rmin,roe)
         end do
      end do
#ifdef PARA
      call reduce(1,rsum)
      call reduce(1,rnegsum)
#endif
      return
      end
!______________________________________________________________________
      subroutine ggapwold(nspin,rhog,gradr,rhor,exc)
!     _________________________________________________________________
!     perdew-wang gga 
!     as given in y-m juan & e kaxiras, prb 48, 14944 (1993) 
!     method by ja white & dm bird, prb 50, 4954 (1994) 
!     non-spin polarized case only
!     _________________________________________________________________
!     by alfredo pasquarello 22/09/1994
!
      use gvec
      use parm
      use cnst
!
      implicit none
!
      integer nspin
      complex(kind=8) rhog(ng)
      real(kind=8) gradr(nnr,3), rhor(nnr), exc
!
      real(kind=8) bb1, bb2, bb3, bb4, bb5, alfa, beta, cc0, cc1, delt, &
     &     c1, c2, c3, c4, c5, c6, c7, a, alfa1, bt1, bt2, bt3, bt4
      parameter(bb1=0.19645,bb2=0.27430,bb3=-0.15084,bb4=0.004,         &
     & bb5=7.7956,alfa=0.09,beta=0.0667263212,cc0=15.75592,             &
     & cc1=0.003521,c1=0.001667,c2=0.002568,c3=0.023266,c4=7.389e-6,    &
     & c5=8.723,c6=0.472,c7=7.389e-2,a=0.0621814,alfa1=0.2137,          &
     & bt1=7.5957,bt2=3.5876,bt3=1.6382,bt4=0.49294,delt=1.0e-12) 
      real(kind=8) x13, x43, x76, pi2, ax, pider1, pider2, pider3,      &
     &     abder1, abder2, abder3
      complex(kind=8) ci
      integer isign, ir
      real(kind=8)                                                      &
     &     aexp, abig, abig2, agr, aroe, byagr, ccr, ccrnum, ccrden,    &
     &     dfxd, dfxdg, dys, dfs, dh1ds, dh1dg, dh1d, dh1dt, dexcdg,    &
     &     dexcd, dh1drs, dh0da, dadec, decdrs, decd, dh0dg, dcdrs,     &
     &     dh0d, dh0dt, eclog, ecr, ecden, fx, fxnum, fxden, fxexp,     &
     &     gkf, grx, gry, grz, h0, h1, h0den, h0arg, h0num,             &
     &     roeth, roe, rs, rs12, rs2, rs3, rs32, s, sd, s2, s3, s4,     &
     &     sysl, t, td, t2, t3, t4, xchge, ys, ysl, ysr
!
!
      if (nspin.ne.1) call error('ggapw','spin not implemented',nspin)
!
      ci=(0.0,1.0)
      x13=1.0/3.0
      x43=4.0/3.0
      x76=7.0/6.0
!     _________________________________________________________________
!     derived parameters from pi
!
      pi2=pi*pi
      ax=-0.75*(3.0/pi)**x13
      pider1=(0.75/pi)**x13
      pider2=(3.0*pi2)**x13
      pider3=(3.0*pi2/16.0)**x13
!     _________________________________________________________________
!     derived parameters from alfa and beta 
!
      abder1=beta*beta/(2.0*alfa)
      abder2=1.0/abder1
      abder3=2.0*alfa/beta
!     _________________________________________________________________
!     main loop
!
      do ir=1,nnr
         roe=rhor(ir)
         if(roe.eq.0.0) goto 100
         aroe=abs(roe)
         grx=gradr(ir,1)
         gry=gradr(ir,2)
         grz=gradr(ir,3)
         agr=sqrt(grx*grx+gry*gry+grz*grz)
         roeth=aroe**x13 
         rs= pider1/roeth
         gkf=pider2*roeth
         sd=1.0/(2.0*gkf*aroe)
         s=agr*sd
         s2=s*s
         s3=s*s2
         s4=s2*s2
!     _________________________________________________________________
!     exchange 
!
         ysr=sqrt(1.0+bb5*bb5*s2)
         ys=bb5*s+ysr
         ysl=log(ys)*bb1
         sysl=s*ysl
         fxexp=exp(-100.0*s2)
         fxnum=1.0+sysl+(bb2+bb3*fxexp)*s2
         fxden=1.0/(1.0+sysl+bb4*s4)
         fx=fxnum*fxden
         xchge=ax*fx*roeth
!     _________________________________________________________________
!     correlation ecr=ec(rho) 
!
         rs12=sqrt(rs)
         rs32=rs12*rs
         rs2=rs*rs
         rs3=rs*rs2
         ecden=a*(bt1*rs12+bt2*rs+bt3*rs32+bt4*rs2) 
         eclog=log(1.0+(1.0/ecden))
         ecr=-a*(1.0+alfa1*rs)*eclog
!     _________________________________________________________________
!     correlation h0(t,ecr)
!     
         td=pider3*sd/rs12
         t=agr*td
         t2=t*t
         t3=t*t2
         t4=t2*t2
         aexp=exp(-abder2*ecr)-1.0
         abig=abder3/aexp
         abig2=abig*abig
         h0num=t2+abig*t4
         h0den=1.0/(1.0+abig*t2+abig2*t4)
         h0arg=1.0+abder3*h0num*h0den
         h0=abder1*log(h0arg)
!     _________________________________________________________________
!     correlation h1(t,s,aroe)
!
         ccrnum=c2+c3*rs+c4*rs2
         ccrden=1.0/(1.0+c5*rs+c6*rs2+c7*rs3)
         ccr=c1+ccrnum*ccrden
         h1=cc0*(ccr-cc1)*t2*fxexp
!     _________________________________________________________________
!     updating of xc-energy
!
         exc=exc+(xchge+ecr+h0+h1)*aroe
!     _________________________________________________________________
!     first part xc-potential from exchange  
!     
         dys=bb5*(1.0+bb5*s/ysr)/ys
         dfs=-fxnum*(ysl+bb1*s*dys+4.0*bb4*s3)*fxden*fxden              &
     &        +(ysl+bb1*s*dys+2.0*s*(bb2+bb3*fxexp)                     &
     &        -200.0*s3*bb3*fxexp)*fxden
         dfxd=(ax*roeth*x43)*(fx-dfs*s)
         dfxdg=ax*roeth*dfs*sd
!     _________________________________________________________________
!     first part xc-potential from ecr 
!
         decdrs=-a*alfa1*eclog*rs + a*(1+alfa1*rs)                      &
     &        *a*(0.5*bt1*rs12+bt2*rs+1.5*bt3*rs32+2.0*bt4*rs2)         &
     &        /(ecden*ecden+ecden)
         decd=-x13*decdrs
!     _________________________________________________________________
!     first part xc-potential from h0 
!     
         dh0da=abder1/h0arg*abder3*h0den*                               &
     &        (t4-h0num*h0den*(t2+2.0*abig*t4))
         dadec=abder3*abder2*(aexp+1.0)/(aexp*aexp)
         dh0d=dh0da*dadec*decd 
         dh0dt=abder1/h0arg*abder3*h0den                                &
     &        *(2.0*t+4.0*abig*t3-h0num*h0den*(2.0*abig*t+4.0*abig2*t3))
         dh0d=dh0d-x76*t*dh0dt
         dh0dg=dh0dt*td
!     _________________________________________________________________
!     first part xc-potential from h1 
!
         dcdrs=(c3+2.0*c4*rs-ccrnum*ccrden*(c5+2.0*c6*rs+3.0*c7*rs2))   &
     &        *ccrden
         dh1drs=cc0*t2*fxexp*dcdrs
         dh1d=-x13*rs*dh1drs
         dh1dt=2.0*t*cc0*(ccr-cc1)*fxexp
         dh1d=dh1d-x76*t*dh1dt
         dh1ds=-200.0*s*cc0*(ccr-cc1)*t2*fxexp
         dh1d=dh1d-x43*s*dh1ds
         dh1dg=dh1dt*td+dh1ds*sd
!     _________________________________________________________________
!     first part of xc-potential: D(rho*exc)/D(rho)
!
         dexcd=dfxd+decd+dh0d+dh1d+ecr+h0+h1
         isign=sign(1.d0,agr-delt)
         byagr=0.5*(1+isign)/(agr+(1-isign)*delt)
         rhor(ir)=dexcd
!
!     gradr = D(rho*exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
         dexcdg=(dfxdg+dh0dg+dh1dg)*aroe*byagr
         gradr(ir,1)=gradr(ir,1)*dexcdg
         gradr(ir,2)=gradr(ir,2)*dexcdg
         gradr(ir,3)=gradr(ir,3)*dexcdg
 100     continue
      end do
!
#ifdef PARA
      call reduce(1,exc)
#endif
      return
      end
!
!-------------------------------------------------------------------------
      subroutine ylmr2b(l,ngy,ngb,gxnb,ylm)
!-----------------------------------------------------------------------
!     real spherical harmonics,  l is combined index for lm (l=1,2...25)
!     order:  s, p_x, p_z, p_y, d_xy, d_xz, d_z^2, d_yz, d_x^2-y^2  ....
!     the real spherical harmonics used here form bases for the
!     irriducible representations of the group O
!
!     see weissbluth  pages 128-130
!     errors in weissbluth have been corrected:
!        1.) elimination of the 7's from l=20 
!        2.) addition of the factor 1./sqrt(12.) to l=25
!
      use cnst
!
      implicit none
      integer l, ngy, ngb
      real(kind=8)  ylm(ngb), gxnb(ngb,3)
      real(kind=8)  c, gsq1, gsq2, gsq3
      integer ig
! 
      if (l.eq.1) then
         c=sqrt(1./fpi)
         do ig=1,ngy
            ylm(ig) = c
         end do
      else if (l.eq.2) then
!     x
         c=sqrt(3./fpi)
!
!       yml is undefined when  g=0 and l>0, so ylm(g=0)=0
!
         ylm(1)=0.         
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)
         end do
      else if (l.eq.3) then
!     z
         c=sqrt(3./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,3)
         end do
      else if (l.eq.4) then
!     y
         c=sqrt(3./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,2)
         end do
      else if (l.eq.5) then
!     x*y
         c=sqrt(15./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*gxnb(ig,2)
         end do
      else if (l.eq.6) then
!     x*z
         c=sqrt(15./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*gxnb(ig,3)
         end do
      else if (l.eq.7) then
!     (3.*z*z-1.0)
         c=sqrt(5.0/fpi/4.0)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*(3.*gxnb(ig,3)*gxnb(ig,3)-1.)
         end do
      else if (l.eq.8) then
!     y*z
         c=sqrt(15./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,2)*gxnb(ig,3)
         end do
      else if (l.eq.9) then
!     x*x-y*y
         c=sqrt(15./fpi/4.)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*(gxnb(ig,1)*gxnb(ig,1)-gxnb(ig,2)*gxnb(ig,2))
         end do
      else if (l.eq.10) then
!     x(x^2-3r^2/5)
         c=sqrt(7./fpi)*5./2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*(gxnb(ig,1)*gxnb(ig,1)-0.6)
         end do
      else if (l.eq.11) then
!     y(y^2-3r^2/5)
         c=sqrt(7./fpi)*5./2.
         ylm(1)=0.
        do ig=2,ngy
           ylm(ig) = c*gxnb(ig,2)*(gxnb(ig,2)*gxnb(ig,2)-0.6)
        end do
      else if (l.eq.12) then
!     xyz
         c=sqrt(7.*15./fpi)
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*gxnb(ig,2)*gxnb(ig,3)
         end do
      else if (l.eq.13) then
!     z(z^2-.6r^2)
         c=sqrt(7./fpi)*5./2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,3)*(gxnb(ig,3)*gxnb(ig,3)-0.6)
         end do
      else if (l.eq.14) then
!     z(x^2-y^2)
         c=sqrt(7.*15./fpi)/2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,3)*                                     &
     &           (gxnb(ig,1)*gxnb(ig,1)-gxnb(ig,2)*gxnb(ig,2))
         end do
      else if (l.eq.15) then
!     y(z^2-x^2)
         c=sqrt(7.*15./fpi)/2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,2)*                                     &
     &           (gxnb(ig,3)*gxnb(ig,3)-gxnb(ig,1)*gxnb(ig,1))
         end do
      else if (l.eq.16) then
!     x(y^2-z^2)
         c=sqrt(7.*15./fpi)/2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*                                     &
     &           (gxnb(ig,2)*gxnb(ig,2)-gxnb(ig,3)*gxnb(ig,3))
         end do
      else if (l.eq.17) then
!     a1
         c=sqrt(3.*7./fpi)*5./4.
         ylm(1)=0.
         do ig=2,ngy
            gsq1=gxnb(ig,1)*gxnb(ig,1)
            gsq2=gxnb(ig,2)*gxnb(ig,2)
            gsq3=gxnb(ig,3)*gxnb(ig,3)
            ylm(ig) = c*(gsq1*gsq1+gsq2*gsq2+gsq3*gsq3-0.6)
         end do
      else if (l.eq.18) then
!     yz(y^2-z^2)
         c=sqrt(9.*35./fpi)/2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,2)*gxnb(ig,3)*                          &
     &           (gxnb(ig,2)*gxnb(ig,2)-gxnb(ig,3)*gxnb(ig,3))
         end do
      else if (l.eq.19) then
!     zx(z^2-x^2)
        c=sqrt(9.*35./fpi)/2.
        ylm(1)=0.
        do ig=2,ngy
           ylm(ig) = c*gxnb(ig,1)*gxnb(ig,3)*                           &
     &          (gxnb(ig,3)*gxnb(ig,3)-gxnb(ig,1)*gxnb(ig,1))
        end do
      else if (l.eq.20) then
!     e\epsilon
         c=sqrt(9.*5./fpi)/4.
         ylm(1)=0.
         do ig=2,ngy
            gsq1=gxnb(ig,1)*gxnb(ig,1)
            gsq2=gxnb(ig,2)*gxnb(ig,2)
            gsq3=gxnb(ig,3)*gxnb(ig,3)
            ylm(ig) = c*( gsq1*gsq1-gsq2*gsq2 -6.*gsq3*(gsq1-gsq2) )
         end do
      else if (l.eq.21) then
!     xy(x^2-y^2)
         c=sqrt(9.*35./fpi)/2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) = c*gxnb(ig,1)*gxnb(ig,2)*                          &
     &           (gxnb(ig,1)*gxnb(ig,1)-gxnb(ig,2)*gxnb(ig,2))
         end do
      else if (l.eq.22) then
!     xy(z^2-1/7*r^2)
        c=sqrt(9.*5./fpi)*7./2.
        ylm(1)=0.
        do ig=2,ngy
           ylm(ig)=c*gxnb(ig,1)*gxnb(ig,2)*(gxnb(ig,3)*gxnb(ig,3)-1./7.)
        end do
      else if (l.eq.23) then
!     zx(y^2-1/7*r^2)
         c=sqrt(9.*5./fpi)*7./2.
         ylm(1)=0.
         do ig=2,ngy 
            ylm(ig) =                                                   &
     &           c*gxnb(ig,1)*gxnb(ig,3)*(gxnb(ig,2)*gxnb(ig,2)-1./7.)
         end do
      else if (l.eq.24) then
!     yz(x^2-1/7*r^2)
         c=sqrt(9.*5./fpi)*7./2.
         ylm(1)=0.
         do ig=2,ngy
            ylm(ig) =                                                   &
     &           c*gxnb(ig,2)*gxnb(ig,3)*(gxnb(ig,1)*gxnb(ig,1)-1./7.)
         end do
      else if (l.eq.25) then
!     e\theta
         c=sqrt(9.*5./fpi/3.)*7./2.
         ylm(1)=0.
         do ig=2,ngy
            gsq1=gxnb(ig,1)*gxnb(ig,1)
            gsq2=gxnb(ig,2)*gxnb(ig,2)
            gsq3=gxnb(ig,3)*gxnb(ig,3)
            ylm(ig) = c*( gsq3*gsq3-0.5*(gsq1*gsq1+gsq2*gsq2)-          &
     &           6./7.*(gsq3-0.5*(gsq1+gsq2)) )
         end do
      else if (l.ge.26) then
         call error(' ylmr2',' higher l not programmed  l=',l)
      endif
      return
      end