module m_mixrho !  mixing routine of smrho, rh. T.kotani think this routine is too compicated to maintain. !It may be better to rewrite all
  ! with keeping the functionality
  use m_ll,only:ll
  use m_lgunit,only:stdo,stml !  integer,parameter,public:: kmxv=15
  public:: mixrho
  private
contains
  integer function ifile_handle() ! find unused file handle
    implicit none
    integer:: i
    logical:: nexist
    integer,save:: irem=2001
    character*256::nnn
    ifile_handle=-999999
    do i=irem,9999
       inquire(unit=i,opened=nexist,name=nnn)
       if( .NOT. nexist) then
          ifile_handle=i
          irem=i+1
          return
       endif
       !         print *,'mpipid i=',mpipid(1),i,trim(nnn)
    enddo
    do i=5001,irem
       inquire(unit=i,opened=nexist)
       if( .NOT. nexist) then
          ifile_handle=i
          irem=i
          return
       endif
    enddo
    call rx('ifile_handle: we did not find open file handle')
  end function ifile_handle
  subroutine mixrho(iter, qval,  sv_p_orhnew, sv_p_orhold, smrnew, smrho,rmsdel)
    use m_struc_def
    use m_supot,only: iv_a_okv,rv_a_ogv,k1,k2,k3
    use m_lmfinit,only:lat_alat,nbas,stdl,ispec,sspec=>v_sspec,nsp,&
         broyinit,nmixinit,betainit,killj,wtinit,wc,bexist,mix_nsave !elinls=>elinl,
    use m_lattic,only: lat_vol
    use m_supot,only:  lat_nabc, lat_ng
    use m_ext,only:    sname
    use m_ftox
    !!= Mix old and new charge densities =
    !!  Takao's version: real space mixing of smrho. It works OK. However,
    !!  we may need to fix it so that this is well-defined chi=|rho-f(rho)|**2 minimization mixing.
    !!     Warn. For Co case, I found broyden mixing works wrong.
    !!     ITER    NIT=30 MIX=B3,n=5,b=.7,w=1,0
    !!     However, as in copt, Broyden mixing is efficient.
    !! ----------------------------------------------------------------------
    ! i Inputs
    ! i iter
    ! i qval  :total valence charge, used to estimate Lindhard parameter
    ! i Read mixing parameter from m_lmfinit
    ! o   orhnew:On input, local parts of the density that generated the ham.
    ! o   orhold:Local parts of the density that the hamiltonian generated
    ! o         :On output, the mixed density
    ! o   smrho :On input, smooth density that generated the hamiltonian H
    ! o         :On output mixed smooth density
    ! o   smrnew:smooth density generated by the H(rhold,smrho) =>DESTROYED on output
    ! o Mixing parameters from m_lmfinit
    ! o         :(1)  broy  :mixing scheme: 0->Anderson 1->Broyden
    ! o         :(2)  beta  :linear mixing beta: nmix = (1-beta)nin + beta*nout
    ! o         :(3)  wc    :Broyden weights for prior iterations
    ! o         :(4,5) weights for n^+ + n^- and n^+ - n^- in spin pol cases
    ! o                Avg, diff weighted by wt(1) and wt(2) in the mixing.
    ! o                As special cases,
    ! o                wt(2)=0 freezes the moments and mixes the charge only
    ! o                wt(1)=0 freezes the charge and mixes the moments only
    ! o         :(6)  mxsav : maximum number of prior iter to save on disk
    ! o         :(7)  mmix  : maximum number of prior iter to mix
    ! o         :(8)  nkill : kill mixing file after multiples nkill
    ! o         :(9)  betv  : special potential mixing parameter
    ! o         :Outputs
    ! o         :(11) rmsdel : rms (output-input) density
    ! o         :(12) rmsl   : rms (output-input) local density
    ! o         :(13) nmix   : number of prior iter actually mixed
    ! o         :(14) broy   : actual mixing scheme used
    ! o         :(15) beta   : actual mixing beta used
    ! o         :(16-24) tj  : Anderson mixing coefficients
    ! o         :(25)        : (spin polarized case only)
    ! o                      : 1  if weight for n^+ + n^- is nonzero
    ! o                      : 10 if weight for n^+ - n^- is nonzero
    ! o                      : 11 if both nonzero

    !l Local variables
    !l   locmix :governs linear transformation of local densities for mixing
    !l          :Note see Bugs, below
    !l          :0  This was the only mode in versions 6.11 and earlier.
    !l          :   The smooth part of the density and the spherical part of
    !l          :   the local densities are mixed in an Anderson or Broyden scheme.
    !l          :   In the remaining options,
    !l          :   local densities are rotated into
    !l          :     rho1,rho2 -> w(r)*(rho1+rho2), rho1-rho2
    !l          :   where w(r) is defined in pvmix9, mode controlled by k9.
    !l          :   For now, k9=10.
    !l          :   This has no effect for linear mixing, but affects the
    !l          :   Anderson or Broyden weights in the spherical part.
    !l          :1  Same as locmix=0, but for the transformation
    !l              CAUTION: not tested in a long time
    !l          :2  Local densities are further mapped onto a G_kL expansion.
    !l          :   whose coefficients are included in the Anderson/Broyden
    !l          :   mixing schemes, including the nonspherical densities.
    !l          :   the residual (what is left of rho beyond the G_kL expansion)
    !l          :   is linearly mixed.
    !l              CAUTION: not tested in a long time
    !l          :3  Similar to locmix=2, but l=0 part is of the G_kL expansion
    !l          :   is projected out.
    !l          :
    !l   rms2   : <(rout-rin)**2>, where rin and rout are the input and
    !l          : screened output densities represented as a vector,
    !l          : including smoothed and local parts.  rmsdel changes with
    !l          : the choice of transformation of local densities;see locmix
    !l          : above.
    !l   rmsdel : (linked to dmxp(11)).
    !l          : On input, rmsdel=<(rout-rin)**2>, defined as rms2, above.
    !l          : On output, rmsdel is overwritten by rms2.
    !l   rmsuns : RMS difference in output-input sm. density, unscreened
    !l   rms    : RMS difference in output-input sm. density, screened
    !l   sumo   : integral of smrho
    !l   sums   : integral of screened smoothed output rho
    !l   kmxs   : k-cutoff for G_kL expansion of screening charge projected
    !l          : into local densities; see Remarks
    !l   kmxr   : k-cutoff for G_kL expansion used to smooth local densities
    !l          : for mixing
    !l   ng0    : condensed number of G vectors. That rho(G) is hermitian
    !l          : is exploited to reduce ng to ng0
    !l   qcell  : cell charge
    
    !b Bugs
    !b   Mixing scheme needs some revision.  For now, uses locmix=3 and k9=10
    !b   Problem: G-kL expansion doesn't represent rim density well; not
    !b          : good for rho1+rho2.
    !b   Better:  1. Use simpler trans: (rho1+rho2)/100 or so, and (rho1-rho2)
    !b            2. For nonspherical parts, mix coefficient to multipole in
    !b               the Anderson/Broyden scheme; linearly mix the rest.
    !r Remarks
    !r   To estimate the self-consistent density a model (Lindhard) dielectric
    !r   function is used.  A smooth "total density" is assembled on a mesh
    !r   from the smooth density + local contributions; the latter are
    !r   smoothed into a G_kL expansion of rho1-rho2 and added to the mesh
    !r   density.  This "mesh total density" is assembled for the input and
    !r   output densities smrho and smrnew, and the screening charge
    !r   (eps^-1 - 1) (smrnew(tot) - smrho(tot)) is obtained using the
    !r   Lindhard approximation for eps^-1.  The screening charge is then:
    !r
    !r     (1)  added to smrnew
    !r
    !r     (2)  added back onto local densities as a G_kL expansion
    !r          to both the true and smooth local mixed densities.
    !r
    !r   Finally, the (screened) output density is mixed with the input density
    !r   using an Anderson or Broyden mixing scheme.  Both mesh part and local
    !r   parts are mixed, the latter consists of the true, smooth,
    !r   and core parts.
    !r   How the mixing proceeds depends on locmix; see Local variables above.
    !r
    !r   The steps in this routine follow this path (see numbers in comments)
    !r   1.   smooth rho+,rho-  ->  smrho, smmag
    !r   2.   FT smrnew(r)-smrho(r) -> cn = smrho(G)
    !r   3.   rho1-rho2 -> smoothed, FT; added to smrho(G)
    !r   4.   cn = screening (eps^-1 - 1) smrho(G)
    !r   5.   Add screening to smrnew; smrnew(G) -> smrnew(r)
    !r   6.   smrnew = sm. screened nout, smrho,smrnew -> smooth rho+,rho-
    !r        At this stage, smrho = input smooth rho+,rho-
    !r                       smrnew = screened output rho+,rho-
    !r   7.   Add screening to rho1, rho2; scale if locmix > 0
    !r   8.   smrho,smrnew(r) -> smrho,smrnew(G)
    !r   9.   Pack in,out densities into large matrix a;
    !r        Read prior iterations into a
    !r   10.  Linear mixing of local densities
    !r   11.  Make modified a when wt(1) or wt(2) restricts q or spin
    !r   12.  Mix input, output densities
    !r   13.  Inverse of step 11.
    !r   14.  Poke contents of a into smrho,orhold
    !r   15.  Undo scaling of local rho (step 7)
    !r
    !m MPI
    !m   master process handles the mix files and broadcasts. All processes
    !m   then mix.
    implicit none
    real(8),save:: rmsdelsave,beta
    integer,save:: broy,nmix,mxsav
    real(8),save:: wt(3)
!    logical,save:: bexist
!    character,save:: fnaminit*8

    include "mpif.h"
    integer :: numprocs, ierr, status(MPI_STATUS_SIZE)
    integer :: MAX_PROCS
    parameter (MAX_PROCS = 100)
    integer :: resultlen
    character*(MPI_MAX_PROCESSOR_NAME) name
    character(10) :: shortname(0:MAX_PROCS-1)
    character(20) :: ext
    character(26) :: datim
    integer :: namelen(0:MAX_PROCS-1)
    double precision :: starttime, endtime
    integer :: iunit
    character(120) :: strn
    logical :: mlog,cmdopt
    integer::  iter , procid , master
    type(s_rv1) :: sv_p_orhold(3,1)
    type(s_rv1) :: sv_p_orhnew(3,1)
    double precision :: qval,elind=0d0
    double complex smrnew(k1,k2,k3,nsp),smrho(k1,k2,k3,nsp)
    integer :: i,i1,i2,i3,ib,ipl,ipr,is,k0,k9, &
         lmxl,n1,n2,n3,ng,nglob,nlml,nr,nmixr,nda,ifi,nlm0, &
         kkk,nnnew,nnmix,ngabc(3),igetss,nx, &
         nkill,isw,naa,kmxr,kmxs,locmix,offx,off2
    integer:: ng0 !,  oa , oaa !, ocn ,owk , oqkl
    integer ,allocatable :: ips0_iv(:)
    real(8),allocatable :: co_rv(:)
    complex(8) ,allocatable :: cg1_zv(:)
    complex(8) ,allocatable :: cg2_zv(:)
    complex(8) ,allocatable :: fkl_zv(:)
    real(8) ,allocatable :: rofi_rv(:)
    real(8) ,allocatable :: rwgt_rv(:)
    real(8) ,allocatable :: wk1_rv(:)
    real(8) ,allocatable :: wk2_rv(:)
    equivalence (n1,ngabc(1)),(n2,ngabc(2)),(n3,ngabc(3))
    parameter (nlm0=49)
    double precision :: a,beta0,dif,difx,difxu,fac,rms,rmt, &
         sumo,summ,sums,top,vol,alat,tpiba,pi,dquns,rmsuns,ddot,q1, &
         qin(2),qout(2),qscr(2),qmix(2),qcell,rms2,rms2f,rsmv,qmx, &
         dgets,rmsdel,srfpi,xx
    double complex xxc
    !      logical parmxp
    character sout*80,fnam*8
    integer ::iwdummy ,isp,nnnx,ng02,ng2, iprint!,ifile_handle
    real(8):: smmin,sss,wgtsmooth,wdummy(1,1,1,1)!,elinl
    complex(8),allocatable:: smrnewbk(:,:,:,:),w_owk(:)
    !      complex(8),allocatable:: w_oqkl(:)
    real(8),allocatable:: w_oqkl(:)
    real(8),allocatable:: w_oa(:),w_oaa(:)
    real(8),allocatable:: w_ocn(:)
    logical::       noelind
    logical:: mixrealsmooth
    logical:: init=.true., initd=.true.
    call tcn('mixrho')
!    if (mixmod == 'none') return
    !  Initial setup for mixrho (charge mixing). This is used in mixrho-parmax.
    call MPI_COMM_RANK( MPI_COMM_WORLD, procid, ierr )
    call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
    call MPI_GET_PROCESSOR_NAME(name, resultlen, ierr)
    call strcop(shortname(procid),name,10,'.',i)
    namelen(procid) = i-1
    master = 0
    mlog = .false. !cmdopt('--mlog',6,0,strn)
    call getpr(ipr)
    ipl = 1
    nx = 0
    kkk = k1*k2*k3
    pi = 4d0*datan(1d0)
    srfpi = dsqrt(4d0*pi)
    locmix = 3 !original lmv7 by mark.
    qmix(2) = 0
    ngabc=lat_nabc
    ng=lat_ng
    vol=lat_vol
    alat = lat_alat
!    elinl = elinls
    kmxs = 3
    kmxr = 8
    if(initd) then
       mxsav = mix_nsave !out 
       rmsdelsave=0d0
       beta=1d0
       initd=.false.
    endif
    rms2 = 0
    fnam = 'mixm' !fnaminit
    wt(1:2)   = wtinit
    !     In case parmxp doesn't touch wt, unset flag 
    wt(3) = 0 !2023feb
    nmix = nmixinit
    broy = broyinit
    if(bexist) beta=betainit
    
!    write(stdo,ftox)'wwwwwwmix2',broy,nmix,ftof(wt),'beta=',ftof(beta),ftof(elinl),ftof(wc),killj!,trim(fnam)
    !r Periodic file deletion (nkill):  parmxp returns nkill as -nkill when
    !r   mod(iter,nkill) is zero, as a signal that the current mixing
    !r   file is to be deleted after it is used.  Here nitj is the sum of
    !r   the number of iterations prior to the current group.
    nkill = max(killj,0)
    if (nkill > 1) then
       if (mod(iter,nkill) == 0) nkill=-nkill
    endif
    if (nkill == 1) nkill=-nkill
    
    rmsdel = rmsdelsave
    if (nmix == 0) broy = 0
    ! ... Initial printout
    !      call awrit7(' Mix density with beta=%;4d'//
    !     .  '%?;n;, Lindhard E=%;3d;;'//
    !     .  '%?;(n==1);, Broyden scheme using n=%,1i;%j;'//
    !     .  '%?;(n==0);, Anderson scheme using n=%,1i;%j;',
    !     .  ' ',80,stdo,beta,isw(elinl.ne.0),elinl,broy,nmix,
    !     .  broy,nmix)
    !      call rx('done')
    if (nsp == 1) wt(2)=0d0
    if (sum(wt**2)==0d0) call rx('MIXRHO: bad mixing weights wt=0')
    call dscal(3,1/dsqrt(wt(1)**2+wt(2)**2+wt(3)**2),wt,1)
    if (nmix < 0) nmix = mxsav
    ! ... Initial charges
    call rhoqm(smrho,k1,k2,k3,n1,n2,n3,nsp,vol,qin)
    call rhoqm(smrnew,k1,k2,k3,n1,n2,n3,nsp,vol,qout)
    fac = vol/(n1*n2*n3)
    !!== RMS differences screened and input sm density; count neg points ==
    sumo = sum(dble(smrho))
    sums = sum(dble(smrnew))
    summ = 0d0
    top = 0d0
    rms = 0d0
    nnnew = 0
    nnmix = 0
    do  i3 = 1, n3
       do  i2 = 1, n2
          do  i1 = 1, n1
             dif = dble(sum(smrnew(i1,i2,i3,1:nsp))-sum(smrho(i1,i2,i3,1:nsp)))
             rms = rms + dif*dif
             top = dmax1(top,dabs(dif))
             xxc = beta*sum(smrnew(i1,i2,i3,1:nsp)) + (1d0-beta)*sum(smrho(i1,i2,i3,1:nsp))
             if( sum(dble(smrnew(i1,i2,i3,1:nsp))) < 0) nnnew = nnnew+1
             if( dble(xxc) < 0) nnmix = nnmix+1
             summ = summ + dble(xxc)
          enddo
       enddo
    enddo
    rms = dsqrt(fac*rms/vol)
    !!== elind mode ==> removed
    do  ib = 1, nbas
       is = ispec(ib) 
       a   =sspec(is)%a
       nr  =sspec(is)%nr
       rmt =sspec(is)%rmt
       lmxl=sspec(is)%lmxl
       rsmv=sspec(is)%rsmv
       if (lmxl < 0) cycle
       nlml = (lmxl+1)**2
       if (nlml > nlm0) call rxi('mixrho: nlml > nlm0, need',nlml)
       allocate(rofi_rv(nr))
       allocate(rwgt_rv(nr))
       call radmsh ( rmt , a , nr , rofi_rv)
       call radwgt ( rmt , a , nr , rwgt_rv)
       call pvmix9 ( 1 , 0 , nr , nlml * nsp , 0 , 0d0 , rofi_rv , sv_p_orhold( 1 , ib )%v &
            , sv_p_orhold( 2 , ib )%v )
       call pvmix9 ( 1 , 0 , nr , nlml * nsp , 0 , 0d0 , rofi_rv , sv_p_orhnew( 1 , ib )%v &
            , sv_p_orhnew( 2 , ib )%v )
       if (allocated(rwgt_rv)) deallocate(rwgt_rv)
       if (allocated(rofi_rv)) deallocate(rofi_rv)
    enddo                     ! Loop over sites
    if (allocated(fkl_zv)) deallocate(fkl_zv)
    if(allocated(w_ocn)) deallocate(w_ocn)
    ! --- Fancy mixing of smoothed + some representation of local rho ---
    ! ... Count number of elts from local densities for fancy mixing scheme
    nda = 0
    do  ib = 1, nbas
       is = ispec(ib)!int(ssite(ib)%spec)
       nr = sspec(is)%nr
       lmxl = sspec(is)%lmxl
       nlml = (lmxl+1)**2
       if (lmxl == -1) goto 20
       !       include spherical part of local densities only
       if (locmix <= 1) then
          nda = nda + 2*nr
          !       also include nonspherical part, G_kL expansion
       elseif (locmix >= 2) then
          nda = nda + 2*(kmxr+1)*nlml
          !         and spherical part of rho1-rho2
          if (locmix == 3) nda = nda + 2*nr
       else
          call rxi('bad locmix,',locmix)
       endif
20     continue
    enddo

    if(mixrealsmooth()) then
       ng02 = k1*k2*k3
       ng2 = ng02
       nda = nda + ng02
       allocate(w_ocn(ng02*nsp),co_rv(ng02*nsp))
       call dcopy(ng02*nsp,dreal(smrnew), 1, w_ocn,1)
       call dcopy(ng02*nsp,dreal(smrho),  1, co_rv,1)
       ! cccccccccccccccccccccc
       wgtsmooth=1d0/ng02**0.5d0 !try this relative weight. Not correct weight charge mixing.
       if(iprint()>10) print *,'wgtsmooth=',wgtsmooth
       ! cccccccccccccccccccccc
       w_ocn= w_ocn*wgtsmooth !
       co_rv= co_rv*wgtsmooth !
    else
       ! ... 8. FFT smooth densities -> rho(G) for mixing
       ! ... Find and include number of smooth mesh points.
       allocate(ips0_iv(ng))
       call lgstar ( 0 , ng , 1 , rv_a_ogv , ng0 , ips0_iv , iwdummy )
       ng02= ng0*2
       ng2=ng*2
       nda = nda + ng02
       allocate(w_ocn(2*ng*nsp),co_rv(ng*2*nsp))
       allocate(w_owk(kkk*nsp))
       call dpcopy(smrnew,w_owk,1,kkk*2*nsp,1d0)
       call fftz3(w_owk,n1,n2,n3,k1,k2,k3,nsp,0,-1)
       call gvgetf ( ng , nsp , iv_a_okv , k1 , k2 , k3 , w_owk, w_ocn)
       call lgstar ( 1 , ng , nsp , rv_a_ogv , i , ips0_iv , w_ocn)
       call dpcopy(smrho,w_owk,1,kkk*2*nsp,1d0)
       call fftz3(w_owk,n1,n2,n3,k1,k2,k3,nsp,0,-1)
       call gvgetf ( ng , nsp , iv_a_okv , k1 , k2 , k3 , w_owk, co_rv)
       call lgstar ( 1 , ng , nsp , rv_a_ogv , i , ips0_iv , co_rv)
       deallocate(w_owk)
    endif
    allocate(w_oa(nda*nsp*(mxsav+2)*2))
    w_oa=0d0
    if(locmix >= 2) then
       allocate(w_oqkl(2*(kmxr+1)*nlm0*nsp*4*nbas))
       w_oqkl=0d0
    else
       allocate(w_oqkl(1))
    endif
    ! --- 9. Read prior iterations from disk; update with current iter ---
    if (procid == master) then
       ifi=ifile_handle()
       open(ifi,file=trim(fnam)//'.'//trim(sname),form='unformatted')
    endif
    call MPI_BCAST(ifi,1,MPI_INTEGER,master,MPI_COMM_WORLD,ierr)
!    if (mlog) then
!       call gettime(datim)
!       call awrit3(' mixrho '//datim//' Process %i of %i on ' &
!            //shortname(procid)(1:namelen(procid))// &
!            ' bcast ifi = %i',' ',256,stml,procid,numprocs,ifi)
!    endif
    k9 = 10
    call pvmix5 ( nmix , mxsav , fnam , ifi , rmsdel , locmix &
         , k9 , nbas , kmxr , nlm0 ,  nsp , sv_p_orhold &
         , sv_p_orhnew , co_rv , w_ocn , ng2 , ng02 , (/0/) , nda , w_oa , w_oqkl &
         , rms2 , nmixr )
    rmsdel = rms2
    nmix = min(nmix,nmixr)
    ! ... Write this and prior iterations onto disk
    call pvmix5 ( nmix , mxsav , fnam , - ifi ,  rmsdel , &
         locmix , k9 , nbas , kmxr , nlm0 ,   nsp , sv_p_orhold &
         , sv_p_orhnew , co_rv , w_ocn , ng , ng02 , (/0/) , nda , w_oa , w_oqkl &
         , rms2 , nmixr )
    if (nkill < 0) then
       if(ipr>=20) then
          sout = ' ...  deleting file '//fnam
          write(stdo,"(a,i5)") trim(sout)//" nkill=",nkill
       endif
       close(ifi,status='delete')
    else
       close(ifi)
    endif

    ! --- Printout smooth density mixing data ---
!    if (nnnew+nnmix > 0 .AND. ipr >= 20) call awrit2( &
!         ' mixrho: (warning) scr. and lin-mixed densities had'// &
!         ' %i and %i negative points',' ',80,stdo,nnnew,nnmix)
    call rhoqm(smrnew,k1,k2,k3,n1,n2,n3,nsp,vol,qscr)
    !     if (ipr .gt. 30) write(stdo,100)
    !    .  sumo*fac,sumu*fac,sums*fac,rms,summ*fac
    if (ipr >45) then
       write(stdo,100) &
            qin(1),qout(1),qscr(1),rms,beta*qscr(1)+(1-beta)*qin(1)
       if (nsp == 2) write(stdo,101) &
            qin(2),qout(2),beta*qscr(2)+(1-beta)*qin(2)
    endif
100 format(' charges:',7x,'old',11x,'new',9x,'screened', &
         6x,'rms diff',7x,'lin mix'/' smooth ',5f14.6)
101 format(' mmom   ',2f14.6,28x,f14.6)
    ! --- 10. Linear mixing of local densities  ---
    call pvmix3 ( nbas , nsp , beta , locmix , wt &
         , kmxr , nlm0 , k9 , w_oqkl , sv_p_orhold , sv_p_orhnew , difx )
    difxu = difx

    !!== Main Mixing part ==
    !!=== 11. Spin polarized case: separate weighting for spin channels ===
    if (nsp == 2 .OR. nx > 0) then
       if (wt(2) /= 0 .AND. nsp == 1) &
            call rx1('mixrho wt(2)=%d invalid when nsp=1',wt(2))
       !       Map to ASA call:
       !             dim. of a     extra    #-mix   #-elts/spin  dim of aa
       !       ASA    nda           nx      naa      npq            naa
       !       FP     nda*nsp       0       naa      nda            naa
       if (nx > 0) call rx('mixrho not ready for nx=0')
       naa = 0
       if (wt(1) /= 0) naa = naa+nda
       if (wt(2) /= 0) naa = naa+nda
       if (wt(3) /= 0) naa = naa+nx
       offx = 0                !offset to extra elements (none now)
       off2 = (nsp-1)*nda      !offset to spin down part of a
       allocate(w_oaa(naa*(mxsav+2)*2))
       w_oaa=0d0
       call pqsclf(0,nda*nsp,nda,offx,off2,nx,naa,mxsav,wt, &
            w_oa,w_oaa,rms2f)
       if (wt(1)*wt(2) /= 0) then
          !          rms2f = rms2f * 2
       else
          write(stdo,ftox)' Constrained spin mixing wt =',ftof(wt), &
               ' Constrained rms DQ=',ftod(rms2f)
          rmsdel = rms2f
       endif
       !!=== 12. Mix the soup of densities ===
       beta0 = beta
       call pvmix6(broy,nmix,nmixr,mxsav,beta,wc,naa,w_oaa)
       !!=== 13. Restore matrix a to rho+, rho===
       call pqsclb(nda*nsp,nda,offx,off2,nx,naa,mxsav,wt,w_oa,w_oaa)
       call pqsclc(nda*nsp,nda*nsp,nx,mxsav,w_oa)
       deallocate(w_oaa)
    else
       naa = nda
       !!=== 12. Mix the soup of densities ===
       beta0 = beta
       call pvmix6(broy,nmix,nmixr,mxsav,beta,wc,naa,w_oa)
    endif
    ! ... 14. Poke mixed smooth and local densities into smrho,rhoold
    allocate(w_owk(kkk))
    call pvmix7 (   nbas , nsp , nda , w_oa , n1 , &
         n2 , n3 , k1 , k2 , k3 , locmix , wt , k9 , kmxr , nlm0 , w_oqkl &
         , ng , ng2, ng02 , iv_a_okv , ips0_iv , rv_a_ogv , co_rv , w_owk , &
         sv_p_orhold , smrho, wgtsmooth )
    deallocate(w_oqkl)
    deallocate(w_owk,w_oa)
    if (allocated(co_rv)) deallocate(co_rv)
    deallocate(w_ocn)
    if (allocated(ips0_iv)) deallocate(ips0_iv)
    ! ... 15. Restore local densities: rho+ +/- rho-  -> rho+, rho-'
    call dpzero(qmix,2)
    do  ib = 1, nbas
       is = ispec(ib)!int(ssite(ib)%spec)
       a=sspec(is)%a
       nr=sspec(is)%nr
       rmt=sspec(is)%rmt
       lmxl=sspec(is)%lmxl
       rsmv=sspec(is)%rsmv
       if (lmxl < 0) cycle
       nlml = (lmxl+1)**2
       if (nlml > nlm0) call rxi('mixrho: nlml > nlm0, need',nlml)
       allocate(rofi_rv(nr))
       allocate(rwgt_rv(nr))
       call radmsh ( rmt , a , nr , rofi_rv)
       call radwgt ( rmt , a , nr , rwgt_rv)
       call pvmix9 ( 1 , - 1 , nr , nlml * nsp , 0 , 0d0 , rofi_rv , &
            sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v )
       call pvmix9 ( 1 , - 1 , nr , nlml * nsp , 0 , 0d0 , rofi_rv , &
            sv_p_orhnew( 1 , ib )%v , sv_p_orhnew( 2 , ib )%v )
       !   ... Add net local charge to qmix
       do  i = 1, nsp
          off2 = 1 + nr*nlml*(i-1)
          allocate(wk1_rv(nr))
          allocate(wk2_rv(nr))
          call dpscop ( sv_p_orhold( 1 , ib )%v , wk1_rv , nr , off2 , 1 , 1d0 )
          call dpscop ( sv_p_orhold( 2 , ib )%v , wk2_rv , nr , off2 , 1 , 1d0 )
          qmx = srfpi*(ddot(nr,wk1_rv,1,rwgt_rv,1)- ddot ( nr , wk2_rv , 1 , rwgt_rv , 1))
          !         print *, 'spin, qmx',i,qmx
          if (i == 2) qmix(2) = qmix(2) + qmx-q1
          q1 = qmx
          qmix(1) = qmix(1) + qmx
          if (allocated(wk2_rv)) deallocate(wk2_rv)
          if (allocated(wk1_rv)) deallocate(wk1_rv)
       enddo
       if (allocated(rwgt_rv)) deallocate(rwgt_rv)
       if (allocated(rofi_rv)) deallocate(rofi_rv)
    enddo                     !Loop over sites

    ! ... Force density positive
    !     call rhopos(smrho,k1,k2,k3,n1,n2,n3)

    ! ... Net interstitial charge and cell charge
    summ = 0
    do    i = 1, nsp
       do    i3 = 1, n3
          do    i2 = 1, n2
             do    i1 = 1, n1
                !       smrho(i1,i2,i3,i) = smrho(i1,i2,i3,i) + .001
                summ = summ + dble(smrho(i1,i2,i3,i))
             enddo
          enddo
       enddo
    enddo
    !     Net system charge
    qcell  = qval - summ*fac - qmix(1)

    !     Add constant to smrho to force charge neutrality
    xx = qcell/(vol*nsp)
    summ = 0

    ! cccccccccccccccccccccccc
    if(iprint()>10) write(6,"(a,d13.5,d13.5)")' mixrho: add corrections to qcell smrho =',qcell,xx
    !       print *,'vvv no xx test'
    !       xx=0d0
    ! ccccccccccccccccccccccc
    do    i  = 1, nsp
       do    i3 = 1, n3
          do    i2 = 1, n2
             do    i1 = 1, n1
                smrho(i1,i2,i3,i) = smrho(i1,i2,i3,i) + xx
             enddo
          enddo
       enddo
    enddo
    !       summ = summ + dble(smrho(i1,i2,i3,i))
    !     qcell  = qval - summ*fac - qmix(1)
    !     print *, 'qcell after shift',qcell; stop

    ! --- Clean up ---
!    if (nmix < 0) dmxp(7) = -nmix
    !    if (beta0 /= beta) betakeep = beta
    rmsdelsave = rmsdel
    ! dmxp(12) = difx
    ! dmxp(13) = iabs(nmix)
    ! dmxp(14) = broy
    ! dmxp(15) = beta
    ! if (broy == 1) dmxp(15) = 1
    ! dmxp(25) = 0
    ! if (wt(1) /= 0 .OR. wt(2) /= 0) dmxp(25) = 1
    ! if (wt(3) /= 0) dmxp(25) = dmxp(25) + 10

    ! ... Printout
    if (ipr >= 10 .AND. abs(qcell) > 1d-6) write(stdo, &
         '('' add q='',f10.6,'' to preserve neutrality'')') qcell

!     if (ipr >= 30) then
!        if (elinl /= 0) then
!           write(stdo,450) ' unscreened ',rmsuns,difxu
!           write(stdo,450) '   screened ',rms,difx,rmsdel
! 450       format(a,'rms difference:  smooth',f10.6,'   local',f10.6: &
!                '   tot',f10.6)
!        endif
!    
!    elseif (ipr >= 20) then
!       call awrit3(' rms smooth dq=%;3g  max local dq=%;3g  dq=%;3g', &
!            ' ',80,stdo,rms,difx,rmsdel)
!    endif
!    if (ipl > 0 .AND. ipr > 0) then
!       call awrit6('%xfp %?;n;elind %;4g  ;%j;'// &
!            '%?#(n==1)#Broyden n %1i wc %;3g#%2j#%-3j'// &
!            '%?#(n==0)#Anderson n %i beta %j%;3g#%3j#', &
!            sout,80,0,isw(elinl.ne.0),elinl,broy,nmix,wc,beta)
!       call awrit3('%a  sm-dq %;3g  mx loc %;3g  dq %;3g',sout,80, &
!            -stdl,rms,difx,rmsdel)
!    endif
    call tcx('mixrho')
    nnnew = 0
    nnmix = 0
    do  i3 = 1, n3
       do  i2 = 1, n2
          do  i1 = 1, n1
             sumo = sumo + dble(smrho(i1,i2,i3,1))
             sums = sums + dble(smrnew(i1,i2,i3,1))
             dif = dble(smrnew(i1,i2,i3,1)-smrho(i1,i2,i3,1))
             rms = rms + dif*dif
             top = dmax1(top,dabs(dif))
             xxc = beta*smrnew(i1,i2,i3,1) + (1d0-beta)*smrho(i1,i2,i3,1)
             if (dble(smrnew(i1,i2,i3,1)) < 0) nnnew = nnnew+1
             if (dble(xxc) < 0) nnmix = nnmix+1
             summ = summ + dble(xxc)
          enddo
       enddo
    enddo
    ! cccccccccccccccccccccccccccccccccc
    do isp=1,nsp
       nnnx = 0
       smmin=1d0
       do i1=1,k1
          do i2=1,k2
             do i3=1,k3
                sss=dreal(smrho(i1,i2,i3,isp))
                if(sss<0d0) then
                   nnnx=nnnx+1
                   if(sss<smmin) then
                      smmin=sss
                   endif
                   !            smrho(i1,i2,i3,isp)=minimumrho
                endif
             enddo
          enddo
       enddo
       if(nnnx>0) then
          write(6,"(a,2i8,d13.5)") ' mixrho: warning. negative smrho; isp number min=',isp,nnnx,smmin
          !          print *,'enforce smrho >',minimumrho
       else
          if(iprint()>45) write(6,"(a,i3)") ' mixrho: all smrho are positive for isp=',isp
       endif
    enddo
  end subroutine mixrho


  subroutine pvmix3 ( nbas , nsp , beta , locmix &
       , wt , kmxr , nlm0 , k9 , qkl , sv_p_orhold , sv_p_orhnew , difx  )
    use m_struc_def  
    use m_lmfinit,only: ispec,sspec=>v_sspec
    !- Linearly mix local densities, possibly subtracting G_kL expansion
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   ssite :struct for site-specific information; see routine usite
    !i     Elts read: spec
    !i     Stored:    *
    !i     Passed to: *
    !i   sspec :struct for species-specific information; see routine uspec
    !i     Elts read: a nr rmt rg lmxl
    !i     Stored:    *
    !i     Passed to: *
    !i   nbas  :size of basis
    !i   nsp   :2 for spin-polarized case, otherwise 1
    !i   beta  :linear mixing parameter
    !l Local variables
    !l   locmix :governs linear transformation of local densities for mixing
    !l          :see description in Local variables, routine mixrho.
    !i   kmxr   :k-cutoff for G_kL expansion when locmix=2; also dim. qkl
    !i   nlm0   :second dimension of qkl
    !i   qkl    :(locmix>=2) G_kL expansion of local densities
    !i          : qkl(:,:,isp,1,ib) = rho1(rhold), G_kL expansion
    !i          : qkl(:,:,isp,2,ib) = rho2(rhold), G_kL expansion
    !i          : qkl(:,:,isp,3,ib) = rho1(rhnew), G_kL expansion
    !i          : qkl(:,:,isp,4,ib) = rho2(rhnew), G_kL expansion
    ! o Inputs/Outputs
    ! o  orhold :On input, local densities generating hamiltonian
    ! o         :transformed into unscaled rho1+rho2 and rho1-rho2 (pvmix9)
    ! o         :If a G_kL expansion has been generated (if locmix>=2)
    ! o         :this expansion is first subtracted from w(orhold).
    ! o         :On output, w(orhold) is overwritten by the linear
    ! o         :combination (1-beta)*w(orhold) + beta*w(orhnew)
    ! o  orhnew :On input, local densities gen. by ham. (maybe screened)
    ! o         :transformed into unscaled rho1+rho2 and rho1-rho2 (pvmix9)
    ! o         :If a G_kL expansion has been generated (if locmix>=2)
    ! o         :this expansion is first subtracted from w(orhnew).
    !o Outputs
    !o   difx   :maximum rms difference in rhonew-rhoold
    !r Remarks
    !r
    !u Updates
    !u   21 Nov 01 First created
    ! ----------------------------------------------------------------------
    !     implicit none
    ! ... Passed parameters
    integer :: nbas,nsp,kmxr,nlm0,locmix,k9
    type(s_rv1) :: sv_p_orhold(3,nbas)

    type(s_rv1) :: sv_p_orhnew(3,nbas)

    ! ino wt(3)      real(8):: difx , beta , wt(2) , qkl(0:kmxr,nlm0,nsp,4,nbas)
    real(8):: difx , beta , wt(3) , qkl(0:kmxr,nlm0,nsp,4,nbas)
!    type(s_site)::ssite(*)
!    type(s_spec)::sspec(*)

    ! ... Local parameters
    integer :: ib,is,igetss,nr,nlml,m,lmxl,k9l
    integer::  i !orsm(4) ,
    real(8) ,allocatable :: ri_rv(:)
    real(8) ,allocatable :: rwgt_rv(:)
    double precision :: a,rmt,rg,difa,rf
    ! ... Heap
    real(8),allocatable:: w_orsm(:,:)
    real(8):: wdummy
    complex(8):: cdummy(1,1,1,1)
    real(8):: rdummy(1,1,1)

    difx  = 0
    do  ib = 1, nbas
       is = ispec(ib)!int(ssite(ib)%spec)

       a=sspec(is)%a
       nr=sspec(is)%nr
       rmt=sspec(is)%rmt
       rg=sspec(is)%rg

       lmxl = int(sspec(is)%lmxl)

       if (lmxl == -1) goto 10
       nlml = (lmxl+1)**2
       allocate(ri_rv(nr))
       allocate(rwgt_rv(nr))
       call radmsh ( rmt , a , nr , ri_rv)
       call radwgt ( rmt , a , nr , rwgt_rv)
       rf = rmt/3

       if (locmix >= 2) then
          allocate(w_orsm(nr*nlml*nsp,4))
       else
          allocate(w_orsm(1,4)) !dummy
       endif

       do  m = 1, 4
          if (locmix >= 2) then
             !           Generate on a mesh the smooth density to be subtracted
             !            call defrr(orsm(m),nr*nlml*nsp)
             i = 1001
             if (wt(1) == 0) i = 11001
             if (wt(2) == 0) i = 21001
             call pkl2ro ( i , 1 , rg , kmxr , nr , nlml , nsp , ri_rv &
                  , rwgt_rv , kmxr , nlm0 , cdummy , qkl ( 0 , 1 , 1 , m , ib) , &
                  w_orsm(1,m) , rdummy , difa)

             !           Undo scaling of rho1+rho2 for linear mix
             if (locmix == 3 .AND. (m == 2 .OR. m == 4)) then
                !             qlk(1,2) -> scaled rho1+rho2, rho1-rho2
                !             qkl(3,4) -> scaled rhn1+rhn2, rhn1-rhn2
                !             if (m .eq. 2) print *, m, 'old'
                !             if (m .eq. 4) print *, m, 'new'
                !             call prrmsh('scaled rhosm1+2 from qkl',w(ori),w(orsm(m-1)),
                !     .         nr,nr,nlml*nsp)
                !             call prrmsh('scaled rhosm1-2 from qkl',w(ori),w(orsm(m)),
                !     .         nr,nr,nlml*nsp)
                k9l = k9 - mod(k9,10)
                call pvmix9 ( k9l , - 1 , nr , nlml * nsp , 0 , rf , ri_rv &
                     , w_orsm(1,m-1) , w_orsm(1,m))
                !             call prrmsh('unscaled rhosm1+2',w(ori),w(orsm(m-1)),
                !     .         nr,nr,nlml*nsp)
                !             call prrmsh('unscaled rhosm1-2',w(ori),w(orsm(m)),
                !     .         nr,nr,nlml*nsp)
             endif
          else
             w_orsm(1,m) = 1
          endif
       enddo

       !   ... Constraints on rnew, smoothed rho
       !        i = 10
       !        if (wt(1) .eq. 0) i = 110
       !        if (wt(2) .eq. 0) i = 210
       !        call pvmix8(i,nr,nlml,nsp,w(orhnew(1,ib)),w(orhnew(2,ib)))
       !        call pvmix8(i,nr,nlml,nsp,w(orsm(1)),w(orsm(2)))
       !        call pvmix8(i,nr,nlml,nsp,w(orsm(3)),w(orsm(4)))

       !   ... Linear mix, this site
       call pvmix4 ( nr , nlml , nsp , ri_rv , locmix , wt , beta , &
            ib , w_orsm ( 1 , 1 ) , w_orsm ( 1 , 2 ) , w_orsm ( 1 , 3 ) , &
            w_orsm ( 1 , 4 ) , sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v &
            , sv_p_orhold( 3 , ib )%v , sv_p_orhnew ( 1 , ib ) %v , sv_p_orhnew &
            ( 2 , ib ) %v , sv_p_orhnew ( 3 , ib ) %v , rwgt_rv , difa )
       difx  = dmax1(difx,difa)
       if (allocated(rwgt_rv)) deallocate(rwgt_rv)
       if (allocated(ri_rv)) deallocate(ri_rv)
       ! i#error, rlse?  orsm = ['orsm(4)', 'integer']
       !        call rlse(orsm(1))
       deallocate(w_orsm)
10     continue
    enddo
  end subroutine pvmix3
  subroutine pvmix4(nr,nlml,nsp,ri,locmix,wt,beta,ib, &
       rhos1,rhos2,rhns1,rhns2, &
       rho1,rho2,rho3,rhn1,rhn2,rhn3,rwgt,dif)
    !- Linear mixing of local densities for one site
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   nr    :number of radial mesh points
    !i   nlml  :L-cutoff for rho1,rho2
    !i   nsp   :2 for spin-polarized case, otherwise 1
    !i   ri    :radial mesh points
    !i   locmix:>2 => work with difference rho-rhos
    !i   beta  :Linear mixing beta
    !i   ib    :site index (printout only)
    !i   rhos1 :input local true density, G_kL expansion
    !i   rhos2 :input local smoothed density, G_kL expansion
    !i   rhns1 :output local true density, G_kL expansion
    !i   rhns2 :output local smoothed density, G_kL expansion
    !i   rho1  :input local true density (generates output rho)
    !i   rho2  :input local smoothed density (generates output rho)
    !i   rho3  :input core charge (generates output rho)
    !i   rhn1  :output local true density (generated by input rho)
    !i   rhn2  :output local smoothed density (generated by input rho)
    !i   rhn3  :output core charge (generated by input rho)
    !i   rwgt  :radial mesh weights
    !o Outputs
    !o   dif   :integral (rho_out - rho_in)**2
    !o   rho1,rho2,rho3 are overwritten by the linearly mixed densities
    !r Remarks
    !u Updates
    !u   09 Dec 08 Suppress q,moment mixing when wt(1) or wt(2) is zero
    ! ----------------------------------------------------------------------
    !     implicit none
    ! ... Passed parameters
    integer :: nr,nlml,nsp,ib,locmix
    double precision :: dif,beta,wt(3)
    double precision :: ri(nr),rwgt(nr), &
         rho1(nr,nlml,nsp),rho2(nr,nlml,nsp),rho3(nr,nsp), &
         rhn1(nr,nlml,nsp),rhn2(nr,nlml,nsp),rhn3(nr,nsp), &
         rhos1(nr,nlml,nsp),rhos2(nr,nlml,nsp), &
         rhns1(nr,nlml,nsp),rhns2(nr,nlml,nsp)
    ! ... Local parameters
    integer :: i,ir,iprint,np,npc
    double precision :: pi,srfpi,vsph,ddot,qnew,qold,qmix,rmt, &
         amnew,amold,amix

    !      print *, 'start of pvmix4'
    !      call prrmsh('rho1',ri,rho1,nr,nr,nlml*nsp) ! rho1+rho2
    !      call prrmsh('rho2',ri,rho2,nr,nr,nlml*nsp) ! rho1-rho2
    !      call prrmsh('rhn1',ri,rhn1,nr,nr,nlml*nsp) ! rhscr1+rhscr2
    !      call prrmsh('rhn2',ri,rhn2,nr,nr,nlml*nsp) ! rhscr1-rhscr2

    !      call prrmsh('rhos1',ri,rhos1,nr,nr,nlml*nsp) ! rhos1+rhos2
    !      call prrmsh('rhos2',ri,rhos2,nr,nr,nlml*nsp) ! rhos1-rhos2
    !      call prrmsh('rhns1',ri,rhns1,nr,nr,nlml*nsp) ! rhscrs1+rhscrs2
    !      call prrmsh('rhns2',ri,rhns2,nr,nr,nlml*nsp) ! rhscrs1-rhscrs2

    ! ... Setup
    !      stdo = lgunit(1)
    pi = 4d0*datan(1d0)
    srfpi = dsqrt(4d0*pi)
    rmt = ri(nr)
    ! ... rho+,rho- -> rho,amag  and some integrated quantities
    call splrho(0,nsp,nr,nlml,rho1,rho2,rho3)
    call splrho(0,nsp,nr,nlml,rhn1,rhn2,rhn3)
    !      qnew = srfpi*(ddot(nr,rhn1,1,rwgt,1) - ddot(nr,rhn2,1,rwgt,1))
    !      qold = srfpi*(ddot(nr,rho1,1,rwgt,1) - ddot(nr,rho2,1,rwgt,1))
    !      amnew = srfpi*(ddot(nr,rhn1(1,1,nsp),1,rwgt,1) -
    !     .               ddot(nr,rhn2(1,1,nsp),1,rwgt,1))
    !      amold = srfpi*(ddot(nr,rho1(1,1,nsp),1,rwgt,1) -
    !     .               ddot(nr,rho2(1,1,nsp),1,rwgt,1))
    qnew = srfpi*ddot(nr,rhn2,1,rwgt,1)
    qold = srfpi*ddot(nr,rho2,1,rwgt,1)
    amnew = srfpi*ddot(nr,rhn2(1,1,nsp),1,rwgt,1)
    amold = srfpi*ddot(nr,rho2(1,1,nsp),1,rwgt,1)

    dif = 0
    do  i  = 1, nsp
       if (wt(1) == 0 .AND. i == 1) cycle
       if (wt(2) == 0 .AND. i == 2) cycle
       do  ir = 1, nr
          dif  = dif + rwgt(ir)*(rhn2(ir,1,i)-rho2(ir,1,i))**2
       enddo
    enddo
    vsph = (4d0*pi/3d0)*rmt**3
    dif = dsqrt(dif/vsph)
    ! --- Overwrite rho with (1-beta)*rho + beta*rhn ---
    !     wt(1)=0 => only mix spin part (set i=2)
    !     wt(2)=0 => only mix charge part (set i=1)
    np = nr*nlml*nsp
    i = 1
    if (wt(1) == 0) then
       np = nr*nlml
       i = 2
    elseif (wt(2) == 0) then
       np = nr*nlml
       i = 1
    endif
    call dscal(np,1-beta,rho1(1,1,i),1)
    call daxpy(np,beta,rhn1(1,1,i),1,rho1(1,1,i),1)
    call dscal(np,1-beta,rho2(1,1,i),1)
    call daxpy(np,beta,rhn2(1,1,i),1,rho2(1,1,i),1)
    npc = nr*nsp
    if (wt(1) == 0 .OR. wt(2) == 0) npc = nr
    call dscal(npc,1-beta,rho3(1,i),1)
    call daxpy(npc,beta,rhn3(1,i),1,rho3(1,i),1)
    ! ... Charges and mag. moments of mixed density
    qmix = srfpi*(ddot(nr,rho2,1,rwgt,1))
    if (nsp == 2) amix = srfpi*ddot(nr,rho2(1,1,nsp),1,rwgt,1)
    ! --- Subtract smoothed (1-beta)*rhos + beta*rhns ---
    if (locmix >= 2) then
       call splrho(20,nsp,nr,nlml,rhos1,rhos2,rho3)
       call splrho(20,nsp,nr,nlml,rhns1,rhns2,rhn3)
       np = nr*nlml*nsp
       call daxpy(np,-(1-beta),rhos1,1,rho1,1)
       call daxpy(np,-beta,rhns1,1,rho1,1)
       call daxpy(np,-(1-beta),rhos2,1,rho2,1)
       call daxpy(np,-beta,rhns2,1,rho2,1)
       call splrho(21,nsp,nr,nlml,rhos1,rhos2,rho3)
       call splrho(21,nsp,nr,nlml,rhns1,rhns2,rhn3)
    endif
    ! ... Restore rho,amag -> rho+,rho-
    call splrho(1,nsp,nr,nlml,rho1,rho2,rho3)
    call splrho(1,nsp,nr,nlml,rhn1,rhn2,rhn3)
    if (iprint() >45 ) then
       write(stdo,100) ib,qold,qnew,qnew,dif,qmix
100    format(' site ',i4,f12.6,5f14.6)
       !  100   format(' site ',i4,f12.6,f14.6,14x,2f14.6)
       if (nsp == 2) write(stdo,101) amold,amnew,amix
101    format(' mmom   ',2f14.6,28x,f14.6)
    endif
    !       print *, 'end of pvmix4'
    !       call prrmsh('rho1 after lin mix',ri,rho1,nr,nr,nlml*nsp)
    !       call prrmsh('rho2 after lin mix',ri,rho2,nr,nr,nlml*nsp)
  end subroutine pvmix4
  subroutine pvmix5 ( nmix , mxsav , fnam , ifi , rmsdel &
       , locmix , k9 , nbas , kmxr , nlm0 , nsp , sv_p_orhold &
       , sv_p_orhnew , co , cn , ng2 , ng02 , cnst , nda , a , qkl , rms2 &
       , nmixr )
    use m_lmfinit,only:ispec,sspec=>v_sspec
    use m_struc_def
    use m_ftox
    !- Copy rho into holding array, read prior iterations from disk
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   nmix  :number of prior iterations sought (for printout)
    !i   mxsav :max no. prior iteration to read, also dimensions a
    !i   fnam  :file name (for printout)
    !i   ifi   :file logical unit; sign delineates whether copying into
    !i         :or out of holding array
    !i   lbin=T
    !i  rmsdel :Same as rms2 (see Outputs), from prior iteration.
    !i         :If no prior iteration, rsmdel=0.  For printout only.
    !i  locmix :switch governing linear transformation of local densities
    !i          for mixing; see Local variables in subroutine mixrho.
    !i  nbas   :size of basis
    !i  ssite  :struct for site-specific information; see routine usite
    !i     Elts read: spec
    !i     Stored:    *
    !i     Passed to: rhogkl
    !i   sspec :struct for species-specific information; see routine uspec
    !i     Elts read: a nr rmt lmxl
    !i     Stored:    *
    !i     Passed to: rhogkl
    !i   nsp   :2 for spin-polarized case, otherwise 1
    !i   orhold:input local density this iteration,
    !i         :transformed into unscaled rho1+rho2 and rho1-rho2 (pvmix9)
    !i   orhnew:output local density this iteration,
    !i         :transformed into unscaled rho1+rho2 and rho1-rho2 (pvmix9)
    !i   co    :input smooth density this iteration, in FT form
    !i   cn    :output smooth density this iteration, in FT form
    !i   ng    :leading dimension of co,cn
    !i   ng0   :number of Fourier coefficients
    !    ng02=ng0*2 in the mixrealsmooth()=F (original case)
    !i   cnst  :site-dependent vector marking whether to constrain
    !i         :selected mixing (not used here)
    !o Inputs/Outputs
    ! o  nmixr :(output for ifi>0) number of prior iter read from file
    ! o        :(input  for ifi<0) number of prior iter written to file
    ! o        :on output number of prior iterations read from file
    ! o  a,nda :work array a and leading dimension (nda must be >= na)
    ! o         a(*,i,1) holds f(xi) for iteration i-1
    ! o         a(*,i,2) holds   xi  for iteration i-1
    !o Outputs
    !o  qkl    :(locmix>=2) G_kL expansion of local densities
    !o         : qkl(:,:,isp,1,ib) = rho1(rhold), G_kL expansion
    !o         : qkl(:,:,isp,2,ib) = rho2(rhold), G_kL expansion
    !o         : qkl(:,:,isp,3,ib) = rho1(rhnew), G_kL expansion
    !o         : qkl(:,:,isp,4,ib) = rho2(rhnew), G_kL expansion
    !o  rms2   :rms difference in f-x, where f=output density and x the
    !o          input density.  f and x include smoothed and local
    !o          contributions, with the latter depending on choice of
    !o          linear transformation; see locmix above.  Thus, there
    !o          is some arbitrariness in rms2.
    !r Remarks
    !r   Content of a(*) = a(*,isp,:,:):
    !r   Case locmix = 0:
    !r     a(1:2*ng0) = Fourier coffs of sm density (hermitian => ng0<ng)
    !r     a(1+2*ng0:2*ng0+2*nr) = l=0 part of rho(1), rho(2)
    !r     ... nda = 2*ng0 + 2*nr
    !r   Case locmix = 1:
    !r     a(1:2*ng0) = Fourier coffs of sm density (hermitian => ng0<ng)
    !r     a(1+2*ng0:2*ng0+2*nr) = Weighted l=0 part of rho(1), rho(2)
    !r     ... nda = 2*ng0 + 2*nr
    !r   Case locmix = 2
    !r     a(1:2*ng0) = Fourier coffs of sm density (hermitian => ng0<ng)
    !r     a(1+2*ng0:2*ng0+(1+kmxr)*nlml) = P_kL expansion of rho1-rho2
    !r     ... nda = 2*ng0 + (1+kmxr)*nlml
    !r   Case locmix = 3
    !r     a(1:2*ng0) = Fourier coffs of sm density (hermitian => ng0<ng)
    !r     a(1+2*ng0:2*ng0+2*nr) = Weighted l=0 part of rho(1), rho(2)
    !r     a(1+k:2*ng0:2*ng0+2*nr+(1+kmxr)*nlml) = P_kL exp. of rho1-rho2
    !r                                             w/ l=0 part zeroed out
    !r     ... nda = 2*ng0 + 2*nr + (1+kmxr)*nlml
    ! -----------------------------------------------------------------
    implicit none
    ! ... Passed parameters
    include "mpif.h"
    integer :: procid, master, numprocs, ierr, status(MPI_STATUS_SIZE)
    integer :: MAX_PROCS
    parameter (MAX_PROCS = 100)
    integer :: resultlen
    character*(MPI_MAX_PROCESSOR_NAME) name
    character(10) :: shortname(0:MAX_PROCS-1)
    character(20) :: ext
    character(26) :: datim
    integer :: namelen(0:MAX_PROCS-1)
    double precision :: starttime, endtime
    character(120) :: strn
    logical :: mlog,cmdopt
    integer :: iunit
    logical :: readerror
    logical :: lbin
    integer :: ng2,ng02,nda,nmix,mxsav,ifi,nbas,nr,nsp,cnst(0:*),locmix, &
         k9,kmxr,nlm0
    type(s_rv1) :: sv_p_orhold(3,1)

    type(s_rv1) :: sv_p_orhnew(3,1)

    double precision :: a(nda,nsp,mxsav+2,2),rms2,rmsdel
    real(8):: co(ng2,nsp) , cn(ng2,nsp) , qkl(0:kmxr,nlm0,nsp,4,1)
!    type(s_site)::ssite(*)
!    type(s_spec)::sspec(*)

    character fnam*8
    ! ... Local parameters
    integer :: ib,na,i,j,k,m,np,iprint,nmixr,is,igetss, off,nlml,lmxl,k9l !awrite,
    real(8) ,allocatable :: rofi_rv(:)

    logical :: lcnst,lfdmp,lddump
    double precision :: ddot,rmt,aat,rf

    character outs*80
    ! ... Heap

    call MPI_COMM_RANK( MPI_COMM_WORLD, procid, ierr)
    call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr)
    call MPI_GET_PROCESSOR_NAME(name, resultlen, ierr)
    call strcop(shortname(procid),name,10,'.',i)
    namelen(procid) = i-1
    master = 0
!    mlog = cmdopt('--mlog',6,0,strn)
    mlog=.false.
    !      stdo = lgunit(1)
    lcnst = cnst(0) .gt. 0

    ! --- Copy rho to a; calc. RMS error; read prior iter from ifi ---
    if (ifi > 0) then

       !   ... FT Mesh density
       na = 1
       do  i = 1, nsp
          call dpcopy(co(1,i),a(na+00,i,1,2),1,ng02,1d0)
          call dpcopy(cn(1,i),a(na+00,i,1,1),1,ng02,1d0)
       enddo
       na = 1 + ng02

       !   ... Site densities, depending on locmix
       do  ib = 1, nbas
          is = ispec(ib)!int(ssite(ib)%spec)


          aat=sspec(is)%a
          nr=sspec(is)%nr
          rmt=sspec(is)%rmt
          lmxl=sspec(is)%lmxl

          if (lmxl == -1) cycle
          rf = rmt/3
          nlml = (lmxl+1)**2
          allocate(rofi_rv(nr))

          call radmsh ( rmt , aat , nr , rofi_rv)


          !   ... If site is constrained, freeze density at rhold (not checked)
          if (lcnst) then
             if (cnst(ib) /= 0) then
                do  i = 1, nsp
                   off = 1+nr*nlml*(i-1)
                   call dpscop ( sv_p_orhold( 1 , ib )%v , sv_p_orhnew ( 1 , ib &
                        ) %v , nr , off , off , 1d0 )


                   call dpscop ( sv_p_orhold( 2 , ib )%v , sv_p_orhnew ( 2 , ib &
                        ) %v , nr , off , off , 1d0 )


                enddo
             endif
          endif

          !   --- Copy local densities to mix matrix a, depending on locmix ---

          !   ... Overwrite rho1,rho2 with linear transf. according to k9
          ! X     mch tmp/rho1 tmp/rho2 -+ out.fe -- -coll 2,18  -px
          !        call prrmsh('rhold1',w(orofi),w(orhold(1,1)),nr,nr,nlml*nsp)
          !        call prrmsh('rhold2',w(orofi),w(orhold(2,1)),nr,nr,nlml*nsp)

          !   ... locmix=0,1,3 : copy spherical part of local rho to a
          !       NB: Increment na after all actions takend w/ local densities
          if (locmix <= 1 .OR. locmix == 3) then

             !         Scaling rho1+rho2, old and new densities
             k9l = k9 - mod(k9,10)
             call pvmix9 ( k9l , 0 , nr , nlml * nsp , 0 , rf , rofi_rv , &
                  sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v )


             call pvmix9 ( k9l , 0 , nr , nlml * nsp , 0 , rf , rofi_rv , &
                  sv_p_orhnew( 1 , ib )%v , sv_p_orhnew( 2 , ib )%v )


             !         call prrmsh('1 post',w(orofi),w(orhold(1,1)),nr,nr,nlml*nsp)
             !         call prrmsh('2 post',w(orofi),w(orhold(2,1)),nr,nr,nlml*nsp)
             !         Copy (scaled) spherical part of rho1+rho2, rho1-rho2 to a
             do  i = 1, nsp
                off = 1+nr*nlml*(i-1)
                call dpscop ( sv_p_orhold( 1 , ib )%v , a ( na + 00 , i , 1 , &
                     2 ) , nr , off , 1 , 1d0 )

                call dpscop ( sv_p_orhold( 2 , ib )%v , a ( na + nr , i , 1 , &
                     2 ) , nr , off , 1 , 1d0 )

                call dpscop ( sv_p_orhnew( 1 , ib )%v , a ( na + 00 , i , 1 , &
                     1 ) , nr , off , 1 , 1d0 )

                call dpscop ( sv_p_orhnew( 2 , ib )%v , a ( na + nr , i , 1 , &
                     1 ) , nr , off , 1 , 1d0 )

             enddo

             !         Undo scaling transformation of rho1+rho2
             !          call pvmix9(k9l,-1,nr,nlml*nsp,0,rf,w(orofi),w(orhold(1,ib)),
             !     .      w(orhold(2,ib)))
             !          call pvmix9(k9l,-1,nr,nlml*nsp,0,rf,w(orofi),w(orhnew(1,ib)),
             !     .      w(orhnew(2,ib)))
             !         call prrmsh('1 rest',w(orofi),w(orhold(1,1)),nr,nr,nlml*nsp)
             !         call prrmsh('2 rest',w(orofi),w(orhold(2,1)),nr,nr,nlml*nsp)

          endif
          !       print *, 'a(nr)'
          !       print *, a(na+nr-1,:,1,:)

          !   ... locmix>=2 : Make G_kL expansion of rho; copy coffs to a
          if (locmix >= 2) then
             do  i = 1, nsp
                call pshpr(iprint()-30)
                call rhogkl ( ib , ib , i , 2 , sv_p_orhold , &
                     kmxr , qkl ( 0 , 1 , i , 1 , ib ) )

                call rhogkl ( ib , ib , i , 3 , sv_p_orhold , &
                     kmxr , qkl ( 0 , 1 , i , 2 , ib ) )

                call rhogkl ( ib , ib , i , 2 , sv_p_orhnew , &
                     kmxr , qkl ( 0 , 1 , i , 3 , ib ) )

                call rhogkl ( ib , ib , i , 3 , sv_p_orhnew , &
                     kmxr , qkl ( 0 , 1 , i , 4 , ib ) )

                call poppr
                !           Convert qkl for 2nd spin channel into 2nd channel only,
                !           since rhogkl makes expansion combining channels 1+2
                if (i == 2) then
                   j = (1+kmxr)*nlm0
                   !             qkl(nsp=2) -> qkl(nsp=1) - qkl(nsp=2)
                   do  m = 1, 4
                      call dpsadd(qkl(0,1,1,m,ib),qkl(0,1,1,m,ib),j,1+j,1,-1d0)
                   enddo
                endif

                !           locmix=3 : zero out spherical part of G_kL
                if (locmix == 3) then
                   !             print *, '!! skip zero qkl'
                   do  m = 1, 4
                      call dpzero(qkl(0,1,i,m,ib),kmxr+1)
                   enddo
                endif

                !           locmix=3 : increment index by 2*nr
                if (locmix == 3) na = na + 2*nr

                !           Copy qkl into a
                np = (1+kmxr)*nlml
                call dpscop(qkl(0,1,i,1,ib),a(na+00,i,1,2),np,1,1,1d0)
                call dpscop(qkl(0,1,i,2,ib),a(na+np,i,1,2),np,1,1,1d0)
                call dpscop(qkl(0,1,i,3,ib),a(na+00,i,1,1),np,1,1,1d0)
                call dpscop(qkl(0,1,i,4,ib),a(na+np,i,1,1),np,1,1,1d0)

                !           locmix=3, : undo increment unless final spin
                if (locmix == 3 .AND. i < nsp) na = na - 2*nr

             enddo                 ! end of spin loop
             !         call prmx('qkl',qkl,1+kmxr,1+kmxr,nlm0*nsp)
          endif

          !       Undo scaling transformation of rho1+rho2
          if (locmix <= 1 .OR. locmix == 3) then
             call pvmix9 ( k9l , - 1 , nr , nlml * nsp , 0 , rf , rofi_rv &
                  , sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v )


             call pvmix9 ( k9l , - 1 , nr , nlml * nsp , 0 , rf , rofi_rv &
                  , sv_p_orhnew( 1 , ib )%v , sv_p_orhnew( 2 , ib )%v )


             !         call prrmsh('1 rest',w(orofi),w(orhold(1,1)),nr,nr,nlml*nsp)
             !         call prrmsh('2 rest',w(orofi),w(orhold(2,1)),nr,nr,nlml*nsp)
          endif

          !   ... Increment index na to a
          if (locmix <= 1) then
             na = na + 2*nr
          elseif (locmix >= 2) then
             na = na + 2*np
          endif

          if (allocated(rofi_rv)) deallocate(rofi_rv)

       enddo                     !Loop over sites
       na = na-1
       if (nda /= na) call rx('mixrho: bug in pvmix5')
       na = na*nsp
       !       call prmx('a',a,nda,na,1)
       rms2 = dsqrt(dabs(ddot(na,a,1,a,1)     - &
            2*ddot(na,a,1,a(1,1,1,2),1) + &
            ddot(na,a(1,1,1,2),1,a(1,1,1,2),1))/(na-0))*nsp
       if (mxsav == 0) return
       !   ... File I/O of previous iterations; nmixr is number actually read
       nmixr = 0
       readerror = .true.
       if (procid == master) then
          readerror = .false.
          !if (lbin)
          read (ifi,err=311,end=311) nmixr, na
          !if ( .NOT. lbin) read (ifi,*,err=311,end=311) nmixr, na
          goto 312
311       continue
          readerror = .true.
312       continue
       endif
       call MPI_BCAST(readerror,1,MPI_LOGICAL, &
            master,MPI_COMM_WORLD,ierr)
       ! if (mlog) then
       !    call gettime(datim)
       !    call awrit3(' mixrho '//datim//' Process %i of %i on ' &
       !         //shortname(procid)(1:namelen(procid))// &
       !         ' bcast readerror = %l',' ', &
       !         256,stml,procid,numprocs,readerror)
       ! endif
       if (readerror) goto 31
       call MPI_BCAST(nmixr,1,MPI_INTEGER,master,MPI_COMM_WORLD,ierr)
       ! if (mlog) then
       !    call gettime(datim)
       !    call awrit3(' mixrho '//datim//' Process %i of %i on ' &
       !         //shortname(procid)(1:namelen(procid))// &
       !         ' bcast nmixr = %i',' ', &
       !         256,stml,procid,numprocs,nmixr)
       ! endif
       call MPI_BCAST(na,1,MPI_INTEGER,master,MPI_COMM_WORLD,ierr)
       ! if (mlog) then
       !    call gettime(datim)
       !    call awrit3(' mixrho '//datim//' Process %i of %i on ' &
       !         //shortname(procid)(1:namelen(procid))// &
       !         ' bcast na = %i',' ', &
       !         256,stml,procid,numprocs,na)
       ! endif
       if (nda*nsp /= na) then
          if (procid == master) then
             write(stdo,ftox)' mixrho:  expecting',nda*nsp,'elements but found',na,'discarding file'
!             call awrit2 &
!                  ( %i elements but found %i ...' &
!                  //' discarding file',' ',80,stdo,nda*nsp,na)
          endif
          nmixr = 0
          goto 31
       endif
       k = min(mxsav,nmixr)
       if (procid == master) then
          readerror = .false.
          do  30  j = 1, k
!             if (lbin) then
                if ( .NOT. lddump(a(1,1,j+1,1),na,ifi)) readerror = .TRUE. 
                if ( .NOT. lddump(a(1,1,j+1,2),na,ifi)) readerror = .TRUE. 
!             else
!                if ( .NOT. lfdmp(a(1,1,j+1,1),na,ifi)) readerror = .TRUE. 
!                if ( .NOT. lfdmp(a(1,1,j+1,2),na,ifi)) readerror = .TRUE. 
!             endif
             if (readerror) goto 313
             nmixr = j
30        enddo
       endif
313    continue
       if (nmixr > 0) then
          call MPI_BCAST(a,nda*nsp*(mxsav+2)*2,MPI_DOUBLE_PRECISION, &
               master,MPI_COMM_WORLD,ierr)
          ! if (mlog) then
          !    call gettime(datim)
          !    call awrit3(' mixrho '//datim//' Process %i of %i on ' &
          !         //shortname(procid)(1:namelen(procid))// &
          !         ' bcast (mix) a %i d.p. numbers',' ', &
          !         256,stml,procid,numprocs,nda*nsp*(mxsav+2)*2)
          ! endif
       endif
31     continue
       if (iprint() >= 20) then !this is needed for test
          write(stdo,fmt=ftox, advance='no')' mixrho: sought',nmix,'iter from file '//trim(fnam)
          write(stdo,fmt='("; read ",g0," RMS DQ=",es8.2e1)',advance='no') nmixr, rms2
          if (rmsdel/= 0) write(stdo,"('  last it=',es8.2e1)",advance='no')rmsdel
          write(stdo,*)
       endif
       ! --- Save iterations on disk ---
    else
       if (procid == master) then
          rewind (-ifi)
          write (-ifi) min(nmix+1,mxsav), nda*nsp
          do  j = 1, min(nmix+1,mxsav)
             write(-ifi)a(:,:,j,1)
             write(-ifi)a(:,:,j,2)
          enddo
       endif
    endif

  end subroutine pvmix5


  subroutine pvmix6(broy,nmix,mmix,mxsav,beta,wc,nda,a)
    use m_amix,only: amix
    !- Mixing of the total density
    ! ------------------------------------------------------------------
    !i  broy   : 0 for Anderson mixing
    !i         : 1 for Broyden mixing
    !i  nmix   : nmix: number of iter to try and mix
    !i  mmix   : number of iterates available to mix
    !i  mxsav  : second dimension of a
    ! o Inputs/Outputs
    ! o  a     :On input:
    ! o         (*,i,1)  output vector for prev. iteration i
    ! o         (*,i,2)  input  vector for prev. iteration i
    ! o        :On output:
    ! o         (*,0,1) and (*,0,2) estimated new vector
    ! o         (*,i,1), i>1 destroyed
    !r  Remarks
    !r    Broyden mixing adapted from Duane Johnson.  Notation:
    !r    x^(m): input vector for iteration m
    !r    F^(m): difference between output and input vector in iteration m
    ! ------------------------------------------------------------------
    !     implicit none
    ! ... Passed parameters
    integer :: nda,nmix,mmix,mxsav,broy
    double precision :: beta,a(nda,0:mxsav+1,2)
    ! ... Local parameters
    integer :: j,iprint!,amix
    double precision :: tjmax,tj(10),norm(100),kpvt(10),wc,wctrue,rms2
    !     integer idamax
    !     double precision dval,rmskm(100),rmsmx
    parameter (tjmax = 5d0)
    ! ... This is the result of linear mixing, q=0 smooth rho
    !     sqmix = a(1,0,1)*beta + a(1,0,2)*(1-beta)
    real(8),allocatable:: aaa(:)
    !$$$ccccccccccccccccccccccccccccccccccccccccc
    !$$$c takao simple mixing test
    !$$$      print *,' vvv: takao simple mixing test'
    !$$$      allocate(aaa(nda))
    !$$$c      aaa = .2d0*a(:,0,1)+.8d0*a(:,0,2)
    !$$$c      aaa = .40d0*a(:,0,1)+.6d0*a(:,0,2)
    !$$$c      aaa = .9d0*a(:,0,1)+.1d0*a(:,0,2)
    !$$$      aaa = 1.d0*a(:,0,1)+.0d0*a(:,0,2)
    !$$$c       aaa = .5d0*a(:,0,1)+.5d0*a(:,0,2)
    !$$$c      aaa = .8d0*a(:,0,1)+.2d0*a(:,0,2)
    !$$$      call dcopy(nda,aaa,1,a(1,0,1),1)
    !$$$      call dcopy(nda,aaa,1,a(1,0,2),1)
    !$$$      deallocate(aaa)
    !$$$      return
    !$$$cccccccccccccccccccccccccccccccccccccccccc
    ! --- Anderson mixing ---
    if (nmix == 0 .OR. nmix == 1 .OR. broy == 0) then
       ! ... amix needs f-x for prior iterations
       do  30  j = 1, nmix
          call daxpy(nda,-1d0,a(1,j,2),1,a(1,j,1),1)
30     enddo
       !        call query('beta',4,beta)
       nmix = amix(nda,nmix,mxsav,0,beta,iprint(),tjmax,a,tj,  rms2)
       !      if (iprint() .ge. 20) write(stdo,101) strn,rms2,(tj(j), j=1,nmix)
       !  101 format(1x,a4,':     rms diff',f10.6:'  tj=',5f8.3)
       !     call prm('a out',a,nda,nda,(mxsav+2)*2)
       ! ... Copy x*=a(*,0,2) to a(*,0,1)
       call dcopy(nda,a(1,0,2),1,a(1,0,1),1)
       ! --- Broyden mixing, Duane Johnson's approach ---
    elseif (broy == 1) then
       call pqmixb(nda,nmix,mmix,mxsav,beta,wc,rms2,a,wctrue)
       !   ... Copy x*=a(*,0,2) to a(*,0,1)
       call dcopy(nda,a(1,0,2),1,a(1,0,1),1)
    else
       call rx('pvmix6: bad value for broy')
    endif
    ! ... Add the change relative to linear mixing
    !     sqmix = (a(1,0,1) - sqmix)*vol
  end subroutine pvmix6

  subroutine pvmix7 ( nbas , nsp , nda , a , n1 &
       , n2 , n3 , k1 , k2 , k3 , locmix , wt , k9 , kmxr , nlm0 , qkl &
       , ng , ng2, ng02 , kv , ips0 , gv , crho , wk , sv_p_orhold , smrho &
       ,wgtsmooth)
    use m_lmfinit,only: ispec,sspec=>v_sspec
    use m_struc_def 
    !- Restore mixed density into specific arrays
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   ssite :struct for site-specific information; see routine usite
    !i     Elts read: spec
    !i     Stored:    *
    !i     Passed to: *
    !i   sspec :struct for species-specific information; see routine uspec
    !i     Elts read: a nr rmt lmxl rg
    !i     Stored:    *
    !i     Passed to: *
    !i   nbas  :size of basis
    !i   nsp   :2 for spin-polarized case, otherwise 1
    !i   nda   :leading dimension of a
    !i   a     :mixed density from pvmix6, smoothed + local densities
    !i         :For local parts, a contains SCALED rho1+rho2, rho1-rho2
    !i         :If locmix =  0,1, or 3, a contains spherical part on mesh
    !i         :If locmix =  2,3 a full rho, GkL expansion
    !i   n1..3 :FT mesh
    !i   k1..3 :dimensions smrho
    !i   locmix:mixing mode for local densities
    !i   ng    :number of FT G-vectors
    !i   ng0   :condensed number of G vector (excluding hermitian equiv)
    !i   kv    :indices for gather/scatter operations (gvlist.f)
    !i   ips0  :permutation array mapping ng0 vectors to ng vectors (lgstar)
    !i   gv    :list of reciprocal lattice vectors G (gvlist.f)
    !i   crho  :FT coefficients of smrho(G)
    !i   wk    :complex work array of dimension (k1,k2,k3)
    !i   smrho :smooth density that generated the hamiltonian
    !i   orhold:local  density that generated the hamiltonian,
    !i         :a portion of which which has been linearly mixed (pvmix4)
    !i         :orhold(1) contains rho1+rho2;
    !i         :orhold(2) contains rho1-rho2; see pvmix9
    !o Outputs
    !o   smrho :overwritten by mixed smooth density
    !o   orhold:mixed local density is returned, in rho1+rho2,rho10rho2 form
    !r Remarks
    !u Updates
    ! ----------------------------------------------------------------------
    implicit none
    ! ... Passed parameters
    integer :: nsp,n1,n2,n3,k1,k2,k3,ng,ng02,nda,na,nr,nbas,locmix,k9, &
         kv(ng,3),ips0(ng),kmxr,nlm0,ng2
    type(s_rv1) :: sv_p_orhold(3,1)

    real(8):: gv(ng,3) , a(nda,nsp) , qkl(0:kmxr,nlm0,nsp,4,nbas) &
         , rf , wt(2), wgtsmooth
!    type(s_site)::ssite(*)
!    type(s_spec)::sspec(*)

    double complex smrho(k1,k2,k3,nsp),wk(k1,k2,k3)
    real(8):: crho(ng2,nsp)
    ! ... Local parameters
    real(8) ,allocatable :: rofi_rv(:)
    real(8) ,allocatable :: rwgt_rv(:)

    integer :: ib,is,igetss,i,m,lmxl,nlml,off,np,k9l,orsm(2)
    double precision :: aat,rmt,rg,xx
    ! ... Heap
    integer ::iwdummy
    complex(8):: cdummy(1,1,1,1)
    real(8):: rdummy(1,1,1)
    real(8),allocatable:: w_orsm(:,:)
    logical:: mixrealsmooth

    ! ... Restore mixed smoothed density
    if(mixrealsmooth()) then
       do  i = 1, nsp
          call dscal(ng02,-1d0/wgtsmooth,crho(1,i),1)
          call daxpy(ng02,1d0/wgtsmooth,a(1,i),1,crho(1,i),1)
          call daxpy(k1*k2*k3,1d0,crho(1,i),1,smrho(1,1,1,i),2)
       enddo
    else
       do  i = 1, nsp
          !       Add difference of mixed, old to minimize perturbation
          ! ccccccccccccccc
          !        print *,' vvv pvmix7 TEST!!! a=0 test'
          !        a=0d0
          ! ccccccccccccc
          call dscal(ng02,-1d0,crho(1,i),1)
          call daxpy(ng02,1d0,a(1,i),1,crho(1,i),1)
          call lgstar(2,ng,1,gv,m,ips0,crho(1,i))
          call gvputf(ng,1,kv,k1,k2,k3,crho(1,i),wk)
          call fftz3(wk,n1,n2,n3,k1,k2,k3,1,0,1)
          call daxpy(2*k1*k2*k3,1d0,wk,1,smrho(1,1,1,i),1)
       enddo
    endif

    ! ... Update local densities
    na = 1 + ng02
    k9l = k9 - mod(k9,10)
    do  ib = 1, nbas
       is = ispec(ib)!int(ssite(ib)%spec)

       aat=sspec(is)%a
       nr=sspec(is)%nr
       rmt=sspec(is)%rmt
       lmxl=sspec(is)%lmxl

       if (lmxl == -1) cycle
       rf = rmt/3
       rg = sspec(is)%rg
       nlml = (lmxl+1)**2
       allocate(rofi_rv(nr))
       allocate(rwgt_rv(nr))
       call radmsh ( rmt , aat , nr , rofi_rv)
       call radwgt ( rmt , aat , nr , rwgt_rv)

       !       Scale rho1+rho2 to match scaled spherical parts and sm parts
       !       call prrmsh('1 mix7a',w(orofi),w(orhold(1,ib)),nr,nr,nlml*nsp)
       call pvmix9 ( k9l , 0 , nr , nlml * nsp , 0 , rf , rofi_rv , &
            sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v )


       !       call prrmsh('1 mix7b',w(orofi),w(orhold(1,ib)),nr,nr,nlml*nsp)

       !       Overwrite sph. rhold with scaled mixed rho; do not unscale
       if (locmix <= 1 .OR. locmix == 3) then
          do  i = 1, nsp
             off = 1+nr*nlml*(i-1)

             call dpscop ( a ( na + 00 , i ) , sv_p_orhold( 1 , ib )%v , nr &
                  , 1 , off , 1d0 )

             call dpscop ( a ( na + nr , i ) , sv_p_orhold( 2 , ib )%v , nr &
                  , 1 , off , 1d0 )

             !            call pvmix9(k9l,-1,nr,1,off-1,rf,w(orofi),w(orhold(1,ib)),
             !     .        w(orhold(2,ib)))
             !          call prrmsh('1 mix7c',w(orofi),w(orhold(1,1)),nr,nr,nlml*nsp)
             !          call prrmsh('2 mix7c',w(orofi),w(orhold(2,1)),nr,nr,nlml*nsp)
          enddo
          na = na + 2*nr
       endif

       !   ... Mixing mode 2,3: add mixed sm density to existing rhold
       if (locmix >= 2) then

          !         call prrmsh('1 scale',w(orofi),w(orhold(1,ib)),nr,nr,nlml*nsp)
          np = (1+kmxr)*nlml
          allocate(w_orsm(nr*nlml*nsp,2))
          do  m = 1, 2
             !            call defrr(orsm(m),nr*nlml*nsp)
             do  i = 1, nsp
                call dpscop(a(na+np*(m-1),i),qkl(0,1,i,m,ib),np,1,1,1d0)
             enddo
             i = 1001
             if (wt(1) == 0) i = 11001
             if (wt(2) == 0) i = 21001
             call pkl2ro ( i , 1 , rg , kmxr , nr , nlml , nsp , rofi_rv , &
                  rwgt_rv , kmxr , nlm0 , cdummy , qkl ( 0 , 1 , 1 , m , ib ) &
                  , w_orsm ( 1 , m ) , rdummy , xx )

          enddo
          !         Don't undo scaling of smoothed density since rho1+rho2 scaled
          !          call pvmix9(k9l,-1,nr,nlml*nsp,0,rf,w(orofi),
          !     .      w(orsm(1)),w(orsm(2)))
          !         Add to rhold(1,2)
          call daxpy ( nr * nlml * nsp , 1d0 , w_orsm ( 1 , 1 ) , 1 , sv_p_orhold( 1 , ib )%v &
               , 1 )

          call daxpy ( nr * nlml * nsp , 1d0 , w_orsm ( 1 , 2 ) , 1 , sv_p_orhold( 2 , ib )%v &
               , 1 )

          deallocate(w_orsm)
          !          call rlse(orsm(2))
          !          call rlse(orsm(1))
          na = na + 2*np
       endif
       !       call prrmsh('1 mix7d',w(orofi),w(orhold(1,ib)),nr,nr,nlml*nsp)

       call pvmix9 ( k9l , - 1 , nr , nlml * nsp , 0 , rf , rofi_rv &
            , sv_p_orhold( 1 , ib )%v , sv_p_orhold( 2 , ib )%v )



       !      call prrmsh('rho1,mixed',w(orofi),w(orhold(1,ib)),nr,nr,nlml*nsp)
       !      call prrmsh('rho2,mixed',w(orofi),w(orhold(2,ib)),nr,nr,nlml*nsp)

       if (allocated(rwgt_rv)) deallocate(rwgt_rv)
       if (allocated(rofi_rv)) deallocate(rofi_rv)
       ! i  bug fix,  unecessary
       ! i         call rlse(orsm(1))

    enddo                     !Loop over sites
    na = na-1
    if (nda /= na) call rx('mixrho: bug in pvmix7')
  end subroutine pvmix7
  subroutine pvmix9(mode,linv,nr,nlml,off,rf,rofi,rho1,rho2)
    !- Transformation of local densities rho1,rho2 for mixing
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   mode  :transformation mode; see Remarks
    !i         : mode  transform
    !i         : 1s digit
    !i         : <=0    do nothing
    !i         :  1     rho1 <= rho1+rho2  rho2 <= rho1-rho2
    !i         : 10s digit
    !i         : 10     scale rho1 by gaussian exp(-(r/rf)**2)
    !i         : 20     scale rho1 by gaussian 1/(1+(r/rf)**2)
    !i   linv  :0  forward transform
    !i         :-1 inverse transform
    !i   nr    :number of radial mesh points
    !i   nlml  :number of L channels
    !i   off   :offset to start of first point (0 for first point)
    !i   rf    :factor used in scaling; see mode
    ! o Inputs/Outputs
    !i   rho1  :On input, local true density, or the transformation
    !i         :of it if the inverse is to be made
    !i         :On output, some linear combination of rho1,rho2.
    !i   rho2  :On input, local smoothed density, or the transformation
    !i         :of it if the inverse is to be made
    !i         :On output, some linear combination of rho1,rho2
    !r Remarks
    !r   This routine takes as input rho1,rho2 and returns in
    !r   rho1,rho2 some linear combination of them, depending on mode
    !r   If linv=1, the inverse transform is done.
    !u Updates
    !u   18 Nov 01 First created
    ! ----------------------------------------------------------------------
    !     implicit none
    ! ... Passed parameters
    integer :: mode,nr,nlml,linv,off
    double precision :: rho1(nr,nlml),rho2(nr,nlml),rofi(nr),rf
    ! ... Local parameters
    integer :: i,ilm,mode0,mode1
    double precision :: ri
    logical:: l_dummy_isanrg,isanrg

    if (mode <= 0) return
    mode0 = mod(mode,10)
    mode1 = mod(mode/10,10)
    ! ino isanrg is logical function,       call isanrg(linv,-1,0,'pvmix9:','linv', .true.)
    l_dummy_isanrg=isanrg(linv,-1,0,'pvmix9:','linv', .true.)
    ! ino isanrg is logical function,       call isanrg(mode0,0,1,'pvmix9:','1s digit mode', .true.)
    l_dummy_isanrg=isanrg(mode0,0,1,'pvmix9:','1s digit mode', .true.)
    ! ino isanrg is logical function,       call isanrg(mode1,0,2,'pvmix9:','10s digit mode', .true.)
    l_dummy_isanrg=isanrg(mode1,0,2,'pvmix9:','10s digit mode', .true.)

    ! ... Forward
    if (linv == 0) then
       if (mode0 == 1) then
          call dsumdf(nr*nlml,1d0,rho1,off,1,rho2,off,1)
       endif

       if (mode1 == 1) then
          do  ilm = 1, nlml
             do  i = 1, nr
                ri = rofi(i)/rf
                rho1(i+off,ilm) = rho1(i+off,ilm) * exp(-ri**2)
             enddo
          enddo

       elseif (mode1 == 2) then
          do  ilm = 1, nlml
             do  i = 1, nr
                ri = rofi(i)/rf
                rho1(i+off,ilm) = rho1(i+off,ilm) / (1 + ri**2)
             enddo
          enddo
       endif

       ! ... Inverse
    else
       if (mode1 == 1) then
          do  ilm = 1, nlml
             do  i = 1, nr
                ri = rofi(i)/rf
                rho1(i+off,ilm) = rho1(i+off,ilm) * exp(ri**2)
             enddo
          enddo

       elseif (mode1 == 2) then
          do  ilm = 1, nlml
             do  i = 1, nr
                ri = rofi(i)/rf
                rho1(i+off,ilm) = rho1(i+off,ilm) * (1 + ri**2)
             enddo
          enddo
       endif

       if (mode0 == 1) then
          call dsumdf(nr*nlml,0.5d0,rho1,off,1,rho2,off,1)
       endif

    endif

    !      call prrmsh('pvmix9 rho1',rofi,rho1,nr,nr,nlml)
    !      call prrmsh('pvmix9,rho2',rofi,rho2,nr,nr,nlml)

  end subroutine pvmix9



  subroutine pqsclf(mode,nda,npq,offx,off2,nx,na,mxsav,wt,a,a2,rms2)
    !- Split into (a+ + a-) and (a+ - a-); include extra data
    ! ------------------------------------------------------------------
    !i Inputs
    !i   mode  :1 a is already split into (a+ + a-) and (a+ - a-)
    !i   nda   :leading dimension of a
    !i   npq   :number of elements to spin-split with wt(1),wt(2)
    !i   offx  :(nx>0)offset to location in a of extra elements
    !i   off2  :offset to spin-down (should be 0 for nsp=1)
    !i   nx    :number of extra data to weight with w3
    !i   na    :dimension of a2 and number of data with nonzero weight
    !i   wt    :1,2 for P,Q (up+dn) and (up-dn), and wt(3) for extra
    !i   a     :(na,1:2,i,1) output vector for prev. iteration i
    !i         :(na,1:2,i,2) input  vector for prev. iteration i
    !o Outputs
    !o   a2    :(*,i,1) output charge and moments, scaled by weights
    !o         :(*,i,2) input  charge and moments, scaled by weights
    !o   rms2  :root mean square of a2(:,0,1)-a2(:,0,2)
    !r Remarks
    !r   q = a2(1) = wt1*(a(1)+a(2)) and mom = a2(2) = wt2*(a(1)-a(2))
    !r   If wt(1) or wt(2) is zero, a2 holds q or mom only.
    !r   na: no. elts to mix: 2*na if wt1,wt2 ne 0, otherwise na
    ! ------------------------------------------------------------------
    !     implicit none
    integer :: mode,nda,npq,nx,mxsav,na,offx,off2
    double precision :: wt(3),a(nda,0:mxsav+1,2),a2(na,0:mxsav+1,2), &
         rms2,ddot
    integer :: is,ia,ja

    !      print *, '!!'; wt(1) = 1 ; wt(2) = 1
    ja = 0
    if (wt(1) == 0 .AND. wt(2) == 0) goto 11
    do  10  is = 0, mxsav+1
       ja = 0
       do  12  ia = 1, npq

          !         if ((ia .eq. -3221 .or. ia .eq. 1) .and. is .eq. 0) then
          !           print *, 'hi',ia,wt(1:2)
          !           print 987, 'rhold+,-',a(ia,is,2),a(ia+off2,is,2)
          !           print 987, 'rhnew+,-',a(ia,is,1),a(ia+off2,is,1)
          !  987      format(a,2x,4f12.6)
          !C          pause
          !         endif

          if (wt(1) /= 0 .AND. mode == 0) then
             ja = ja+1
             !           Given (rhnew+ + rhnew-)*wt(1) ,  (rhold+ + rhold-)*wt(1)
             a2(ja,is,1) = (a(ia,is,1) + a(ia+off2,is,1))*wt(1)
             a2(ja,is,2) = (a(ia,is,2) + a(ia+off2,is,2))*wt(1)
          elseif (wt(1) /= 0) then
             ja = ja+1
             a2(ja,is,1) = a(ia,is,1)*wt(1)
             a2(ja,is,2) = a(ia,is,2)*wt(1)
          endif
          if (wt(2) /= 0 .AND. mode == 0) then
             ja = ja+1
             !           Given (rhnew+ - rhnew-)*wt(2) ,  (rhold+ - rhold-)*wt(2)
             a2(ja,is,1) = (a(ia,is,1) - a(ia+off2,is,1))*wt(2)
             a2(ja,is,2) = (a(ia,is,2) - a(ia+off2,is,2))*wt(2)
          elseif (wt(2) /= 0) then
             ja = ja+1
             a2(ja,is,1) = a(ia+off2,is,1)*wt(2)
             a2(ja,is,2) = a(ia+off2,is,2)*wt(2)
          endif
12     enddo
10  enddo
11  continue

    if (wt(3) /= 0) then
       do  20  ia = 1, nx
          ja = ja+1
          do  22  is = 0, mxsav+1
             a2(ja,is,1) = a(ia+offx,is,1)*wt(3)
             a2(ja,is,2) = a(ia+offx,is,2)*wt(3)
22        enddo
20     enddo
    endif

    rms2 = dsqrt(dabs(ddot(na,a2,1,a2,1)     - &
         2*ddot(na,a2,1,a2(1,0,2),1) + &
         ddot(na,a2(1,0,2),1,a2(1,0,2),1))/(na-0))

    if (ja /= na) call rx('pqsclf: element mismatch')

  end subroutine pqsclf
  subroutine pqsclb(nda,npq,offx,off2,nx,na,mxsav,wt,a,a2)
    !- Undo split into wt1*q and wt2*mom done by pqsclf
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   nda   :leading dimension of a
    !i   npq   :number of P,Q
    !i   offx  :offset to location in a of extra elements
    !i   off2  :offset to spin-down (should be 0 for nsp=1)
    !i   nx    :number of extra data
    !i   na    :dimensions a2; number of data with nonzero weight
    !i   wt    :1,2 for P,Q (up+dn) and (up-dn), and wt(3) for extra
    !i   a2    :mixed q,mom scaled by weights, or portion of a with nonzero weights
    !o Outputs
    !o   a     :a2 is unscaled and restored into a
    !r Remarks
    !u Updates
    !u   10 Dec 08
    ! ----------------------------------------------------------------------
    !     implicit none
    integer :: nda,na,npq,nx,mxsav,offx,off2
    double precision :: wt(3),a(nda,0:mxsav+1,2), &
         a2(na,0:mxsav+1,2),sum,diff
    integer :: is,ia,ja

    if (wt(1) == 0 .AND. wt(2) == 0) then
       ja = 0
    elseif (wt(1) /= 0 .AND. wt(2) /= 0) then
       do  is = 0, mxsav+1
          ja = 1
          do  ia = 1, npq
             !          if ((ia .eq. 3221 .or. ia .eq. -3) .and. is .eq. 0) then
             !            print *, 'hi',ia,wt(1:2)
             !            print 987, 'mxold+,-',
             !     .        (a2(ja,is,2)/wt(1) + a2(ja+1,is,2)/wt(2))/2,
             !     .        (a2(ja,is,2)/wt(1) - a2(ja+1,is,2)/wt(2))/2
             !            print 987, 'mxnew+,-',
             !     .        (a2(ja,is,1)/wt(1) + a2(ja+1,is,1)/wt(2))/2,
             !     .        (a2(ja,is,1)/wt(1) - a2(ja+1,is,1)/wt(2))/2
             !C 987       format(a,2x,4f12.6)
             !C           pause
             !          endif
             !         mixed  rhonew+,rhold+
             a(ia,is,1)     = (a2(ja,is,1)/wt(1) + a2(ja+1,is,1)/wt(2))/2
             a(ia,is,2)     = (a2(ja,is,2)/wt(1) + a2(ja+1,is,2)/wt(2))/2
             !         mixed  rhonew-,rhold-
             a(ia+npq,is,1) = (a2(ja,is,1)/wt(1) - a2(ja+1,is,1)/wt(2))/2
             a(ia+npq,is,2) = (a2(ja,is,2)/wt(1) - a2(ja+1,is,2)/wt(2))/2
             ja = ja+2
          enddo
       enddo
       ja = 2*npq
    elseif (wt(1) /= 0) then
       diff = 0
       do  ia = 1, npq
          do  is = 0, mxsav+1
             ja = ia+off2
             sum  = a2(ia,is,1)/wt(1)
             diff = (a(ia,is,1) - a(ja,is,1))
             a(ia,is,1) = (sum + diff)/2
             a(ja,is,1) = (sum - diff)/2
             sum  = a2(ia,is,2)/wt(1)
             diff = (a(ia,is,2) - a(ja,is,2))
             a(ia,is,2) = (sum + diff)/2
             a(ja,is,2) = (sum - diff)/2
          enddo
       enddo
       ja = npq
    elseif (wt(2) /= 0) then
       do  ia = 1, npq
          do  is = 0, mxsav+1
             ja = ia+npq

             !         if ((ia .eq. -3221 .or. ia .eq. 1) .and. is .eq. 0) then
             !            print *, 'hi',ia,wt(1:2)
             !            print 987, 'given qold, mixed (+)-(-)',
             !     .        (a(ia,is,2) + a(ja,is,2)),a2(ia,is,2)
             !            sum  = (a(ia,is,2) + a(ja,is,2))
             !            diff = a2(ia,is,2)/wt(2)
             !            print 987, 'mixed +,-                ',
             !     .        (sum + diff)/2, (sum - diff)/2
             !  987       format(a,2x,4f12.6)
             !          endif

             !         given rhnew+ + rhnew- ; mixed  rhnew+ - rhnew-
             sum  = (a(ia,is,1) + a(ja,is,1))
             diff = a2(ia,is,1)/wt(2)
             a(ia,is,1) = (sum + diff)/2
             a(ja,is,1) = (sum - diff)/2
             !         sum = given rhold+ + rhold- ; diff = mixed  rhold+ - rhold-
             sum  = (a(ia,is,2) + a(ja,is,2))
             diff = a2(ia,is,2)/wt(2)
             a(ia,is,2) = (sum + diff)/2
             a(ja,is,2) = (sum - diff)/2
          enddo
       enddo
       ja = npq
    endif

    ! --- Extra mixing data ---
    if (wt(3) /= 0) then
       do  is = 0, mxsav+1
          do  ia = 1, nx
             a(ia+offx,is,1) = a2(ia+ja,is,1)/wt(3)
             a(ia+offx,is,2) = a2(ia+ja,is,2)/wt(3)
          enddo
       enddo
    endif

  end subroutine pqsclb

  subroutine pqsclc(nda,npq,nx,mxsav,a)
    !- Copy a(:,0,2) into a(:,0,1)
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   nda   :leading dimension of a
    !i   npq   :number of P,Q
    !i   nx    :number of extra data
    ! o Inputs/Outputs
    ! o  a     :On input
    ! o        :(*,0,1) mixed output vector
    ! o        :(*,0,2) mixed input vector
    ! o        :On output
    ! o        :(*,0,2) -> (*,0,1)
    !r Remarks
    !u Updates
    !u   10 Dec 08 First created
    ! ----------------------------------------------------------------------
    implicit none
    integer :: nda,npq,nx,mxsav
    double precision :: a(nda,0:mxsav+1,2)
    integer :: is,ia

    is = 0
    do  ia = 1, npq+nx
       !        if (a(ia,is,1) .ne. a(ia,is,2)) then
       !          call rx('pqsclc: oops')
       !        endif
       a(ia,is,1) = a(ia,is,2)
    enddo
    !     stop 'here'
  end subroutine pqsclc

  subroutine pqmixb(nda,nmix,mmix,mxsav,beta,wc,rms2,a,wctrue)
    !      use m_lgunit,only:stdo
    use m_ftox
    !- Broyden mixing of a vector, Duane Johnson's approach
    ! ------------------------------------------------------------------
    !i  mmix: number of iterates available to mix
    !i  a:    (*,i,1)  output values for prev. iteration i
    !i        (*,i,2)  input  values for prev. iteration i
    ! o nmix: nmix > 0: number of iter to try and mix
    !i        nmix < 0: use mmix instead of nmix.
    !o  nmix: (abs)  number of iter actually mixed.
    !o        (sign) <0, intended that caller update nmix for next call.
    !r  Notations:
    !r  x^(m): input vector for iteration m
    !r  F^(m): difference between output and input vector in iteration m
    ! ------------------------------------------------------------------
    implicit none
    integer :: nda,nmix,mmix,mxsav
    double precision :: beta,rms2,wctrue,a(nda,0:mxsav+1,2)
    double precision :: ddot,dval,wc
    integer:: im , km , i , iprint , i1mach , imix , jmix
    real(8) ,allocatable :: xmp1_rv(:)
!    integer:: broyj
    real(8) ,allocatable :: f_rv(:)
    real(8) ,allocatable :: ui_rv(:)
    real(8) ,allocatable :: vti_rv(:)
    real(8) ,allocatable :: xold_rv(:)
    real(8) ,allocatable :: df_rv(:)
    real(8) ,allocatable :: dx_rv(:)
    real(8) ,allocatable :: wk_rv(:)
    ! --- Allocate some arrays ---
    allocate(xmp1_rv(nda))
    allocate(f_rv(nda))
    allocate(ui_rv(nda))
    allocate(vti_rv(nda))
    allocate(xold_rv(nda))
    allocate(df_rv(nda))
    allocate(dx_rv(nda))
    ! ... imix is a local copy of nmix
    imix = nmix
    if (imix < 0) imix = mmix
    ! --- Starting from iteration mmix, build the Jacobian matrix ---
1   jmix = min(mmix,iabs(imix))
    allocate(wk_rv(nda*2*(jmix+2)))

    do  10  km = 1, jmix
       !   ... this loops from most-distant to most-recent
       im = jmix-km+1
       !   ... this loops from most-recent to most-distant
       !       im = km

       call dcopy ( nda , a ( 1 , im - 1 , 1 ) , 1 , dx_rv , 1 )

       call daxpy ( nda , - 1d0 , a ( 1 , im - 1 , 2 ) , 1 , dx_rv &
            , 1 )

       rms2 = dsqrt ( ddot ( nda , dx_rv , 1 , dx_rv , 1 ) / ( &
            nda - 0 ) )


       ! ---   Determine wc_true if wc < 0 ---
       if (wc < 0) then
          wctrue = -wc/100/dsqrt(nda*rms2**2)
          wctrue = min(max(wctrue,1d0),1d4)
       else
          wctrue = wc
       endif
       if (km == 1) wctrue = .01d0

       i = iprint()
       if (km /= jmix) i = i-20
       i = broyj ( nda , a ( 1 , im - 1 , 2 ) , dx_rv , km , 0 , &
            i , beta , 0d0 , 0d0 , 0d0 , wctrue , wk_rv , nda , xmp1_rv &
            )

10  enddo

    ! --- Check for interactive change of nmix ---
    ! NB negative sign signals request for permanent change in nmix
    im = imix
    !      if (iprint() .gt. 30) call query('redo, nmix=',2,imix)
    if (iabs(imix) > mmix .AND. imix /= im) &
         write(stdo,ftox)' (warning) only ',mmix,' iter available'
    if (im /= imix) goto 1
    nmix = imix
    ! ... If no prior iter allowed, give up on nmix
    ! i#error, have return with len(w_varlist)>0 at line 769
    if ( nmix == 0 ) then
       if (allocated(wk_rv)) deallocate(wk_rv)
       if (allocated(dx_rv)) deallocate(dx_rv)
       if (allocated(df_rv)) deallocate(df_rv)
       if (allocated(xold_rv)) deallocate(xold_rv)
       if (allocated(vti_rv)) deallocate(vti_rv)
       if (allocated(ui_rv)) deallocate(ui_rv)
       if (allocated(f_rv)) deallocate(f_rv)
       if (allocated(xmp1_rv)) deallocate(xmp1_rv)
       return
    endif


    ! --- Printout ---
    if (iprint() > 60 .OR. (iprint() >= 40 .AND. nda <= 100)) then
       print 310
       do  12  i = 1, nda
          if ( dabs ( a ( i , 0 , 1 ) - a ( i , 0 , 2 ) ) >= 5d-9 ) &
               print 311 , i , a ( i , 0 , 2 ) , a ( i , 0 , 1 ) , a ( i , 0 &
               , 1 ) - a ( i , 0 , 2 ) , xmp1_rv(i) 
12     enddo
311    format(i5,4f14.6)
310    format(14x,'Old',11X,' New',9X,'Diff',10X,'Mixed')
    endif
    ! --- Save x^(m+2) into a(*,0,2) and exit ---
    call dcopy ( nda , xmp1_rv , 1 , a ( 1 , 0 , 2 ) , 1 )
    if (allocated(xmp1_rv)) deallocate(xmp1_rv)
    if (allocated(f_rv)) deallocate(f_rv)
    if (allocated(ui_rv)) deallocate(ui_rv)
    if (allocated(vti_rv)) deallocate(vti_rv)
    if (allocated(xold_rv)) deallocate(xold_rv)
    if (allocated(df_rv)) deallocate(df_rv)
    if (allocated(dx_rv)) deallocate(dx_rv)
    if (allocated(wk_rv)) deallocate(wk_rv)
  end subroutine pqmixb
  subroutine pqmxup(na,mxsav,nclass,nl,nsp,nx,lmx,pnu,qnu,xnew,pold,qold,xold,cnst,nda,a,rms2)
    !- Copy from holding array into P,Q
    implicit none
    integer :: nda,mxsav,nclass,nl,nsp,nx,lmx(nclass),cnst(0:*)
    double precision :: a(nda,0:mxsav+1,2),rms2, &
         pnu(nl,nsp,nclass), qnu(3,nl,nsp,nclass), xnew(nx), &
         pold(nl,nsp,nclass),qold(3,nl,nsp,nclass),xold(nx),ddot
    integer :: ic,isp,l,na,i
    logical :: lcnst
    lcnst = cnst(0) .gt. 0
    na = 1
    rms2 = 0
    do  110  isp = 1, nsp
       do  11  ic = 1, nclass
          if (lcnst) then
             if (cnst(ic) /= 0) goto 11
          endif
          do  12  l = 0, lmx(ic)
             !       print *, pold(l+1,isp,ic) - a(na,0,2)
             pnu(l+1,isp,ic) = a(na,0,2)
             rms2 = rms2 + (qold(1,l+1,isp,ic) - a(na+1,0,2))**2
             do  21  i = 1, 3
                !       print *, i, qold(i,l+1,isp,ic) -  a(na+i,0,2)
                qnu(i,l+1,isp,ic) = a(na+i,0,2)
21           enddo
             na = na+4
12        enddo
11     enddo
110 enddo
    call dcopy(nx,a(na,0,2),1,xnew,1)
    rms2 =  rms2 + ddot(nx,xnew,1,xnew,1) - &
         2*ddot(nx,xnew,1,xold,1) + &
         ddot(nx,xold,1,xold,1)
    na = na + nx
    rms2 = dsqrt(4*rms2/(na-0))
  end subroutine pqmxup
  subroutine dpsadd(adest,asrc,nel,n1,n2,fac)
    integer :: n1,n2,i,iadd,ntop,nel
    double precision :: asrc(1),adest(1),fac
    call daxpy(nel,fac,asrc(n2),1,adest(n1),1)
  end subroutine dpsadd
  subroutine pkl2ro(mode,ib,rsm,kmax,nr,nlml,nsp,rofi,rwgt,k0,nlm0,fklc,fklr,rho1,rho2,qmx)
    !- Put PkL or GkL expansion of a function on a radial mesh for one site
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   mode  :a compound of digits :
    !i         :1s digit
    !i         :  0 add P_kL expansion of rho to rho1 (and possibly rho2)
    !i         :    Here, fkl are coffs to P_kL expansion
    !i         :  1 add G_kL expansion of rho to rho1 (and possibly rho2)
    !i         :    Here, fkl are coffs to G_kL expansion
    !i         :10s digit
    !i         :  0 add expansion to rho1 only; rho2 is not touched
    !i         :  1 add expansion to both rho1 and rho2
    !i         :100s digit
    !i         :  0 coefficients fkl are real (uses fklr)
    !i         :  1 coefficients fkl are complex (uses fklc)
    !i         :1000s digit
    !i         :  0 do not initialize rho1,rho2 before adding expansion
    !i         :  1 initialize rho1,rho2 before adding expansion
    !i         :10000s digit for spin polarized case
    !i         :  1 zero out charge, rho1+ + rho1- (and rho2+ + rho2-)
    !i         :  2 zero out spin, rho1+ - rho1- (and rho2+ + rho2-)
    !i   ib    :site index (used in addressing fkl expansion)
    !i   rsm   :smoothing radius for P_kL (or G_kL) expansion
    !i   kmax  :k-cutoff for P_kL (or G_kL) expansion
    !i   nr    :number of radial mesh points
    !i   nlml  :L-cutoff for rho1,rho2
    !i   nsp   :2 if fkl is spin-pol; 1 if not
    !i   rofi  :radial mesh points
    !i   rwgt  :radial mesh weights
    !i   k0    :leading dimension of fkl
    !i   nlm0  :second dimension of fkl
    !i   fklc  :complex coefficients to P_kL (or G_kL) expansion for this site
    !i         :(only one of fklc or fklr is used; see mode)
    !i   fklr  :real coefficients to P_kL (or G_kL) expansion for this site
    !i         :(only one of fklc or fklr is used; see mode)
    !o Outputs
    !o  rho1  :fkL P_kL (or G_kL) added to local true density rho1
    !o  rho2  :fkL P_kL (or G_kL) added to local smoothed density rho2
    !i   qmx   :charge in rho1 after fkl PkL is added
    !r Remarks
    !u Updates
    !u   21 Nov 01 First created
    ! ----------------------------------------------------------------------
    !     implicit none
    ! ... Passed parameters
    integer :: mode,ib,k0,kmax,nlm0,nr,nlml,nsp
    double precision :: qmx,rsm
    double precision :: rofi(nr),rwgt(nr)
    double precision :: rho1(nr,nlml,nsp),rho2(nr,nlml,nsp)
    double complex   fklc(0:k0,nlm0,nsp,ib)
    double precision :: fklr(0:k0,nlm0,nsp,ib)
    ! ... Local parameters
    integer :: lmx,i,ilm,isp,k,l,lmxl,mode0,mode1,mode2,mode3,mode4,np
    double precision :: add,pi,r,rl,srfpi,sum1
    parameter(lmx=10)
    double precision :: fkl(0:kmax,nlm0,2),pkl(0:kmax,0:lmx)

    pi = 4d0*datan(1d0)
    srfpi = dsqrt(4*pi)
    lmxl = ll(nlml)
    mode0 = mod(mode,10)
    mode1 = mod(mode/10,10)
    mode2 = mod(mode/100,10)
    mode3 = mod(mode/1000,10)
    mode4 = mod(mode/10000,10)

    if (lmxl > lmx) call rxi('pklr2o: increase lmx, need',lmxl)
    if (nlml > nlm0) call rxi('pklr2o: increase nlm0, need',nlml)

    do  isp = 1, nsp
       do  ilm = 1, nlml
          do  k = 0, kmax
             if (mode2 == 0) then
                fkl(k,ilm,isp) = fklr(k,ilm,isp,ib)
             else
                fkl(k,ilm,isp) = dble(fklc(k,ilm,isp,ib))
             endif
          enddo
       enddo
    enddo

    if (mode3 == 1) call dpzero(rho1,nr*nlml*nsp)
    if (mode3 /= 0 .AND. mode1 /= 0) call dpzero(rho2,nr*nlml*nsp)

    do  i = 2, nr
       r = rofi(i)
       if (mode0 == 0) then
          call radpkl(r,rsm,kmax,lmxl,kmax,pkl)
       else
          call radgkl(r,rsm,kmax,lmxl,kmax,pkl)
       endif
       do  ilm = 1, nlml
          l = ll(ilm)
          rl = r**l
          do  isp = 1, nsp
             if (mode1 /= 0) then
                do  k = 0, kmax
                   add = fkl(k,ilm,isp)*pkl(k,l)*r*r*rl
                   rho1(i,ilm,isp) = rho1(i,ilm,isp) + add
                   rho2(i,ilm,isp) = rho2(i,ilm,isp) + add
                enddo
             else
                do  k = 0, kmax
                   add = fkl(k,ilm,isp)*pkl(k,l)*r*r*rl
                   rho1(i,ilm,isp) = rho1(i,ilm,isp) + add
                enddo
             endif
          enddo
       enddo
    enddo
    sum1 = 0d0
    do  isp = 1, nsp
       do  i = 1, nr
          sum1 = sum1 + rwgt(i)*rho1(i,1,isp)
       enddo
    enddo
    qmx = sum1*srfpi

    ! zero out spin or charge
    if (mode4 == 0 .OR. nsp /= 2) return
    i = 20                   ! No core
    if (mode1 == 0) i = 30 ! No rho2
    call splrho(i,nsp,nr,nlml,rho1,rho2,sum1)
    np = nr*nlml
    if (mode4 == 1) then   ! Zero density
       call dpzero(rho1(1,1,1),np)
       if (mode1 /= 0) then ! Including rho2
          call dpzero(rho2(1,1,1),np)
       endif
    endif
    if (mode4 == 2) then   ! Zero spin
       call dpzero(rho1(1,1,nsp),np)
       if (mode1 /= 0) then ! Including rho2
          call dpzero(rho2(1,1,nsp),np)
       endif
    endif
    call splrho(i+1,nsp,nr,nlml,rho1,rho2,sum1)

  end subroutine pkl2ro

  subroutine rhoqm(smrho,k1,k2,k3,n1,n2,n3,nsp,vol,qsum)
    !- Return charge, magnetic moment of smooth density
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   smrho :smooth density on uniform mesh
    !i   k1..k3:
    !i   n1..n3:
    !i   nsp   :2 for spin-polarized case, otherwise 1
    !i   vol   :cell volume
    !o Outputs
    !o   qsum  :qsum(1) = smrho(+) + smrho(-)
    !o         :qsum(2) = smrho(+) - smrho(-) (for nsp=2 only)
    !l Local variables
    !l         :
    !r Remarks
    !r   Input smrho is assumed to be (rho1, rho2)
    !r   If instead smrho=(rho1+rho2,rho1-rho2) => qsum(1,2) = q+amom, q-amom
    !u Updates
    !u   13 Dec 08 First created
    ! ----------------------------------------------------------------------
    implicit none
    integer :: k1,k2,k3,n1,n2,n3,nsp
    double complex smrho(k1,k2,k3,nsp)
    double precision :: vol,qsum(2)
    integer :: i,i1,i2,i3
    double precision :: sumi,q1,fac
    qsum(1) = 0
    qsum(2) = 0
    fac = vol/(n1*n2*n3)
    q1 = 0
    do  i = 1, nsp
       sumi = 0
       do  i3 = 1, n3
          do  i2 = 1, n2
             do  i1 = 1, n1
                sumi = sumi + dble(smrho(i1,i2,i3,i))
             enddo
          enddo
       enddo
       if (i == 2) qsum(2) = qsum(2) + q1-sumi
       q1 = sumi
       qsum(1) = qsum(1) + sumi
    enddo
    qsum(1) = fac*qsum(1)
    qsum(2) = fac*qsum(2)
    !     write(*,333) qsum ! 333 format(' rhoqm : istl charge, moment = ',2f13.7)
  end subroutine rhoqm

  integer function broyj(n,xin,gin,ir,isw,ipr,beta,dxmx,xtol,gtol, wc,wk,ndw,xnew)
    !- One Broyden step in finding gin = f[xin]-xin = 0
    ! ----------------------------------------------------------------------
    !i Inputs
    !i   n:     number of variables
    !i   ir:    Number of iterations of x and g.
    !i          1 initiates a new sequence of mixing;
    !i          broyj uses linear mixing for this iteration.
    !i   isw    1s digit (not implemented)
    !i          0  find minimum
    !i          1  find maximum
    !i         10s digit not used
    !i        100s digit not used
    !i       1000s digit governs convergence criterion:
    !i          1 return when |grad| < gtol
    !i          2 return when max dx < xtol
    !i          3 return when either (1) or (2) satisfied
    !i          4 return when both (1) and (2) are satisfied
    !i   beta:  linear mixing parameter (ir=1 only)
    !i   xin:   input vector, this iteration
    !i   gin:   output-input vector, f[xin]-xin, this iteration
    !i   wc:    weighting for this iteration
    ! o  wk     workspace of 2*ndw*(ir+2), ndw>=n
    ! o         wk must be preserved between calls to broyj.
    ! o         (*,1,0) x of the prior iteration.
    ! o         (*,2,0) g of the prior iteration.
    ! o         (*,1..2,1..ir-1) u and vt of this and prior iterations
    ! o         (*,1,ir) g(this iter) - g (prior iter).
    !o Outputs
    !o   xnew   estimate of x
    !o   broyj
    !r Remarks
    !r   Adapted from Duane Johnson
    ! ----------------------------------------------------------------------
    implicit none
    integer :: isw,ir,n,ipr,ndw
    double precision :: beta,dxmx,wc,xin(n),gin(n),xnew(n),xtol,gtol, wk(ndw,2,0:ir)
    ! Local variables
    integer :: i,ip,j,k,irm1,irm2,lm,ln,nn,i1mach,isw1,isw2,isw3 !dinv
    integer :: ierr
    double precision :: aij,cmj,dfnorm,fac1,fac2,gmi,one,zero,ddot,w0
    parameter (zero=0d0,one=1d0,nn=20)
    double precision :: a(nn,nn),cm(nn),w(nn),d(nn,nn)
    double precision :: betx,diff,gmax,xmax
    !     double precision wl(nn,3),u(nn,nn),v(nn,nn)
    save w,cm,a,w0

    isw1 = mod(isw/10,10)
    isw2 = mod(isw/100,10)
    isw3 = mod(isw/1000,10)
    if (ir > nn) call rxi('broyj: increase nn, need',ir)

    ! --- First iteration: simple mixing ---
    if (ir == 1) then
       w0 = wc
       betx = beta
       gmax = 0
       do  k = 1, n
          gmax = max(gmax,abs(gin(k)))
       enddo
       If (dxmx > 0d0 .AND. gmax > dxmx) then
          betx = beta*dxmx/gmax
       endif
       do    k = 1, n
          xnew(k) = xin(k) + betx*gin(k)
       enddo
       ! --- Subsequent iterations: Broyden mixing ---
    else
       !   ... Make xold, gold
       do  20  k = 1, n
          wk(k,1,0) = xin(k) - wk(k,1,0)
          wk(k,1,ir) = gin(k) - wk(k,2,0)
20     enddo

       !   --- Coefficient matrices and the sum for corrections ---
       !   ... dfnorm = |g(i)-g(i-1)|, used for normalization
       dfnorm = dsqrt(ddot(n,wk(1,1,ir),1,wk(1,1,ir),1))
       fac2 = one/(dfnorm+1d-12)
       fac1 = beta*fac2
       !   ... Shuffle each prior u,vt to prior+1 iteration
       irm1 = ir-1
       irm2 = ir-2
       do  30  j = irm2, 1, -1
          call dcopy(n,wk(1,1,j),1,wk(1,1,j+1),1)
          call dcopy(n,wk(1,2,j),1,wk(1,2,j+1),1)
30     enddo
       !   ... Make u,vt for this iteration
       do    k = 1, n
          wk(k,1,1) = fac1*wk(k,1,ir) + fac2*wk(k,1,0)
          wk(k,2,1) = fac2*wk(k,1,ir)
       enddo

       !   --- Make  a and b = ( w0**2 I + a )^-1 (symmetric) ---
       do  42  j = 1, irm2
          aij = zero
          cmj = zero
          do    k = 1, n
             cmj = cmj + wk(k,2,ir-j)*gin(k)
             aij = aij + wk(k,2,ir-j)*wk(k,2,1)
          enddo
          a(irm1,j) = aij
          a(j,irm1) = aij
          cm(j) = cmj
42     enddo
       aij = zero
       cmj = zero
       do  k = 1, n
          cmj = cmj + wk(k,2,1)*gin(k)
          aij = aij + wk(k,2,1)*wk(k,2,1)
       enddo
       a(irm1,irm1) = aij
       cm(irm1) = cmj
       w(irm1) = wc

       !   ... Set up and calculate beta matrix
       do   lm = 1, irm1
          do  ln = 1, irm1
             d(ln,lm) = a(ln,lm)*w(ln)*w(lm)
          enddo
          d(lm,lm) = w0**2 + a(lm,lm)*w(lm)*w(lm)
       enddo

       !   --- Invert to make d ---
       !        if (dinv(' ',irm1,nn,d) .ne. 0) then
       !          call rx('broyj: matrix singular')
       !       endif
       call matinv2(irm1,d(1:irm1,1:irm1),ierr)

       !   ... Invert with singular value decomposition
       !        call svd(nn,irm1,irm1,d,wl(1,2),.true.,u,.true.,v,ierr,wl)
       !        call dpzero(d,nn**2)
       !        do  60  ln = 1, irm1
       !   60   d(ln,ln) = 1
       !        call svbksb(nn,irm1,irm1,irm1,wl(1,2),u,v,d,d,wl)
       !    ... This one sometimes hangs up
       !        call rs(nn,irm1,d,wl(1,3),1,v,wl,wl(1,2),ierr)
       !        call dpzero(d,nn**2)
       !        do  60  ln = 1, irm1
       !          print *, 'evl',ln, wl(ln,3)
       !   60   d(ln,ln) = 1
       !        call svbksb(nn,irm1,irm1,irm1,wl(1,3),v,v,d,d,wl)

       !   --- xnew <- vector for the new iteration ---
       do   k = 1, n
          xnew(k) = xin(k) + beta*gin(k)
       enddo
       do  70  i = 1, irm1
          gmi = zero
          do  ip = 1, irm1
             gmi = gmi + cm(ip)*d(ip,i)*w(ip)
          enddo
          do  k = 1, n
             xnew(k) = xnew(k) - gmi*wk(k,1,ir-i)*w(i)
          enddo
70     enddo
       !   ... Cap to maximum allowed shift xnew-xin
       if (dxmx > 0d0) then
          diff = 0
          do  k = 1, n
             diff = max(diff,abs(xnew(k)-xin(k)))
          enddo
          if (diff > dxmx) then
             betx = dxmx/diff
             do   k = 1, n
                xnew(k) = xin(k) + betx*(xnew(k)-xin(k))
             enddo
          endif
       endif
    endif
    ! --- Cleanup, setup for next call ---
    xmax = 0
    gmax = 0
    diff = 0
    do  110  k = 1, n
       xmax = max(xmax,abs(xnew(k)-xin(k)))
       gmax = max(gmax,dabs(gin(k)))
       diff = diff + (xnew(k)-xin(k))**2
       wk(k,2,0) = gin(k)
       wk(k,1,0) = xin(k)
110 enddo
    diff = dsqrt(diff/n)
    j = ir+1
    if (isw3 /= 0 .AND. (gmax == 0 .OR. &
         gmax < gtol .AND. xmax < xtol .AND. isw3 == 4 .OR. &
         gmax < gtol .AND. (isw3 == 1 .OR. isw3 == 3)  .OR. &
         xmax < xtol .AND. (isw3 == 2 .OR. isw3 == 3)  .OR. &
         gmax < gtol .AND. (isw3 == 1 .OR. isw3 == 3))) j = 0
    broyj = j
  end function broyj
subroutine rhogkl ( ib1 , ib2 , nsp , mode , sv_p_orhoat , kmax , qkl )!- G_kL expansion of valence sphere densities
  use m_lgunit,only:stdo
  use m_struc_def  
  use m_lmfinit,only: ispec,sspec=>v_sspec
  use m_hansr,only:corprm
  !i Inputs
  !i  ib1,ib2: compute expansion coefficents for sites ib1..ib2
  !i   nsp   :1 make qkl for first spin (possibly the only one)
  !i         :2 make qkl combining spins 1 and 2
  !i   mode  : a compound of digits specifying what is to be included
  !i         : in the expansion coefficients
  !i         : 1s   digit = 1 include local density rho1-rho2
  !i         :              2 include local density rho1
  !i         :              3 include local density rho2
  !i         : 10s  digit = 1 include core density rhoc
  !i                        2 include -1 * core density from sm-hankel
  !i                        3 combination 1+2
  !i         : 100s digit = 1 add -1 * nuclear density Y0 Z delta(r)
  !i   kmax  :make expansion coffs to polynomial cutoff kmax
  !i   ssite :struct for site-specific information; see routine usite
  !i     Elts read: spec
  !i     Stored:    *
  !i     Passed to: *
  !i   sspec :struct for species-specific information; see routine uspec
  !i     Elts read: lmxl z qc a nr rmt rg
  !i     Stored:    *
  !i     Passed to: corprm
  !i   orhoat:vector of offsets containing site density
  !o Outputs
  !o   qkl  :Expansion coefficients, stored as a single long vector.
  !o        := integral pkl Y_L integrand
  !o        :where integrand is according to mode
  !r Remarks
  !r    In the spin-polarized case, up- and down- spins are combined.
  !u Updates
  !u   19 Oct 01 Adapted from rhomom.f
  ! ----------------------------------------------------------------------
  implicit none
  integer:: ib1 , ib2 , nsp , mode , kmax
  type(s_rv1) :: sv_p_orhoat(3,1)
  real(8):: qkl(0:kmax,1)
  integer:: ipr , iprint, nrmx , j1 , ib , is &
       , igetss , lmxl , nr , nlml , ilm , j , lfoc , k , l , m
  real(8) ,allocatable :: rofi_rv(:)
  real(8) ,allocatable :: rwgt_rv(:)
  real(8) ,allocatable :: h_rv(:)
  parameter( nrmx=1501)
  double precision :: z,qc,a,rmt,qcorg,qcorh,qsc,cofg,cofh,rg, ceh,rfoc,df(0:20)
  ipr  = iprint()
  allocate(rofi_rv(nrmx))
  allocate(rwgt_rv(nrmx))
  call stdfac(20,df)
  if (ipr >= 40) write(stdo,221)
  j1 = 1
  do  ib = ib1, ib2
     is = ispec(ib)
     lmxl=sspec(is)%lmxl
     z=sspec(is)%z
     qc=sspec(is)%qc
     a=sspec(is)%a
     nr=sspec(is)%nr
     rmt=sspec(is)%rmt
     rg=sspec(is)%rg

     if (lmxl == -1) goto 10
     call corprm(is,qcorg,qcorh,qsc,cofg,cofh,ceh,lfoc,rfoc,z)
     qc = qcorg+qcorh
     nlml = (lmxl+1)**2
     call radmsh ( rmt , a , nr , rofi_rv )
     call radwgt ( rmt , a , nr , rwgt_rv )
     allocate(h_rv(nr*(kmax+1)*(lmxl+1)))
     call pvrgkl ( mode , kmax , nlml , nr , nsp , rofi_rv , rwgt_rv &
          , sv_p_orhoat( 1 , ib )%v , sv_p_orhoat( 2 , ib )%v , sv_p_orhoat( 3 , ib )%v &
          , h_rv , cofh , rg , ceh , rfoc , z , qkl ( 0 , j1 ) )
     deallocate(h_rv)

     if (ipr >= 40) then
        write(stdo,222) ib,0,1,(qkl(k,j1), k=0,kmax)
        ilm = 1
        do  l = 1, lmxl
           do  m = -l, l
              ilm = ilm+1
              j = j1+ilm-1
              if (dabs(qkl(0,j))*df(2*l+1) > 1d-6) write(stdo,220) 0,ilm, &
                   (qkl(k,j)*df(2*l+1),k=0,kmax)

           enddo
        enddo
     endif
222  format(2x,'ib=',i3,i5,i6,10f12.6)
220  format(9x,i4,i6,f12.6,10f12.6)
221  format(/' rhogkl:    k   ilm      qkl (2l+1)!! ...')
     !$$$#if DEBUG
     !$$$        if (ib .eq. 1) then
     !$$$          print *, 'ib=',ib
     !$$$          call prtrkl ( mode , kmax , rg , nr , nlml , nsp , rofi_rv ,
     !$$$     .     sv_p_orhoat( 1 , ib )%v , sv_p_orhoat( 2 , ib )%v , sv_p_orhoat( 3 , ib )%v
     !$$$     .     , qkl ( 0 , j1 ) )
     !$$$
     !$$$
     !$$$        endif
     !$$$#endif
     j1 = j1+nlml
10   continue
  enddo
  if (allocated(rwgt_rv)) deallocate(rwgt_rv)
  if (allocated(rofi_rv)) deallocate(rofi_rv)
end subroutine rhogkl
subroutine pvrgkl(mode,kmax,nlml,nr,nsp,rofi,rwgt,rho1,rho2,rhoc, pkl,cofh,rg,ceh,rfoc,z,qkl)
  use m_hansr,only:hansmr
  use m_ll,only:ll
  !- Multipole moments for one site
  ! ----------------------------------------------------------------------
  !i Inputs
  !i   mode  : a compound of digits specifying what is to be included
  !i         : in the expansion coefficients
  !i         : 1s   digit = 1 include local density rho1-rho2
  !i         :              2 include local density rho1
  !i         :              3 include local density rho2
  !i         : 10s  digit = 1 include core density rhoc
  !i                        2 include -1 * core density from sm-hankel
  !i                        3 combination 1+2
  !i         : 100s digit = 1 add -1 * nuclear density Y0 Z delta(r)
  !i   kmax  :k-cutoff for polynomial expansion of radial part
  !i   nlml  :L-cutoff for charge
  !i   nr    :number of radial mesh points
  !i   nsp   :number of spins
  !i   rofi  :radial mesh points
  !i   rwgt  :radial integration weights
  !i   rho1  :local true density*r**2, tabulated on a radial mesh
  !i   rho2  :local smoothed density*r**2, tabulated on a radial mesh
  !i   rhoc  :core density
  !i   cofh  :coefficient to Hankel part of pseudocore density (corprm)
  !i   rg    :smoothing radius for compensating gaussians
  !i   ceh   :energy of hankel function to fit core tail
  !i   rfoc  :smoothing radius for hankel head fitted to core tail
  !i   z     :nuclear charge
  !o Outputs
  !o   qkl  :expansion coefficients for rho
  !w Workarea:
  !w   pkl:
  !r Remarks
  !r   Q_kL = integral p_kl (rho1-rho2) + l=0 contr. from core spillout
  !r   The core spillout term is:
  !r      qcore(rhoc)-z  - sm_qcore-sm_qnuc
  !r   pvrgkl makes this Q_kL when mode=131; partial contr for other modes
  !r   NB: p0l = a**l and scaling factor for k=0 is 4*pi/(a**l * (2l+1)!!)
  !r       => q0l = 4*pi/(2l+1)!! q_l, where q_l is the multipole moment
  ! ----------------------------------------------------------------------
  !     implicit none
  ! ... Passed parameters
  integer :: mode,kmax,nlml,nr,nsp
  double precision :: ceh,cofh,rfoc,rg,z
  double precision :: rofi(1),rwgt(1),qkl(0:kmax,nlml), &
       rhoc(nr,nsp),rho1(nr,nlml,nsp),rho2(nr,nlml,nsp), &
       pkl(nr,0:kmax,0:*)
  ! ... Local parameters
  integer :: n0,i,ilm,l,m,lmxl,isp,k
  parameter (n0=10)
  double precision :: ag,fac,y0,xi(0:n0),fpi,factk,dfact, &
       df(0:20),wk(nr),smrch,f1,f2
  !     double precision gkl(0:kmax,0:nlml),wk2(0:kmax,0:nlml),rl
  !     double precision sumh,samh,y0

  !     call prrmsh('rho1',rofi,rho1,nr,nr,nlml*nsp)
  !     call prrmsh('rho2',rofi,rho2,nr,nr,nlml*nsp)

  fpi  = 16d0*datan(1d0)
  y0   = 1d0/dsqrt(fpi)
  lmxl = ll(nlml)
  call stdfac(20,df)
  call vecpkl(rofi,rg,nr,kmax,lmxl,nr,kmax,wk,1,pkl,pkl)
  !     Non-vectorized form ... should be able to integrate with
  !     either pkl -> gkl exp, or with gkl -> pkl exp, but
  !     something is wrong ...  doesn't work
  !      do  i = 1, nr
  !        call radgkl(rofi(i),rg,kmax,lmxl,kmax,wk2)
  !C       call radpkl(rofi(i),rg,kmax,lmxl,kmax,wk2)
  !        rl = 1
  !        do  l = 0, lmxl
  !          pkl(i,0:kmax,l) = wk2(0:kmax,l)*rl
  !          rl = rl * rofi(i)
  !        enddo
  !      enddo
  call dpzero(qkl,nlml*(kmax+1))

  ! ... rho1-rho2 contribution (or rho1, or rho2, depending on mode)
  if (mod(mode,10) > 0) then
     if (mod(mode,10) == 1) then
        f1 = 1
        f2 = -1
     elseif (mod(mode,10) == 2) then
        f1 = 1
        f2 = 0
     elseif (mod(mode,10) == 3) then
        f1 = 0
        f2 = 1
     else
        call rx('rhogkl: bad mode')
     endif
     ilm = 0
     do  l = 0, lmxl
        do  m = -l, l
           ilm = ilm+1
           do  k = 0, kmax
              do  i = 1, nr
                 !             If rg is small enough, these should all integrate to 1
                 !             call radgkl(rofi(i),rg,kmax,lmxl,kmax,gkl)
                 !              if (m .eq. -l) qkl(k,ilm) = qkl(k,ilm) +
                 !     .          rwgt(i)*rofi(i)**(2+l) * gkl(k,l) * pkl(i,k,l)
                 do  isp = 1, nsp
                    qkl(k,ilm) = qkl(k,ilm) + rwgt(i) * pkl(i,k,l) * &
                         (f1*rho1(i,ilm,isp) + f2*rho2(i,ilm,isp))
                 enddo
              enddo
           enddo
        enddo
     enddo
  endif

  ! ... Core part (spec'd by 10s digit mode)
  if (mod(mode/10,10) > 0) then
     do  k = 0, kmax
        do  isp = 1, nsp
           !         Case 1 or 3: add core density
           if (mod(mod(mode/10,10),2) /= 0) then
              do  i = 1, nr
                 qkl(k,1) = qkl(k,1) + y0*rwgt(i)*rhoc(i,isp)*pkl(i,k,0)
              enddo
           endif
           !         Case 2 or 3: subtract core density from sm. Hankel
           if (mod(mode/10,10) >= 2) then
              do  i = 1, nr
                 call hansmr(rofi(i),ceh,1/rfoc,xi,1)
                 smrch = cofh*xi(0)*rofi(i)**2
                 qkl(k,1) = qkl(k,1) - rwgt(i)*smrch*pkl(i,k,0)
              enddo
           endif
        enddo
     enddo
  endif
  ! ... Nuclear part (spec'd by 100s digit mode)
  if (mod(mode/100,10) == 1) then
     do  k = 0, kmax
        qkl(k,1) = qkl(k,1) - y0*z*pkl(1,k,0)
     enddo
  endif
  ! ... Scale to get coefficients of the G_kL; see radpkl
  ag = 1/rg
  ilm = 0
  dfact = 1
  do  l = 0, lmxl
     do  m = -l, l
        ilm = ilm+1
        factk = 1d0
        do  k = 0, kmax
           fac = fpi / ((4*ag*ag)**k * ag**l * factk * dfact)
           qkl(k,ilm) = qkl(k,ilm) * fac
           factk = factk*(k+1)
        enddo
     enddo
     dfact = dfact*(2*l+3)
  enddo
end subroutine pvrgkl
subroutine vecpkl(r,rsm,nr,kmax,lmax,nrx,k0,wk,lrl,p,gp)! Vector of p_kl polynomials, or r^l p_kl
  ! ----------------------------------------------------------------------
  !i Inputs
  !i   r     :vector of points
  !i   rsm   :smoothing radius
  !i   nr    :number of points
  !i   kmax  :make p for k=0..kmax
  !i   lmax  :make p for l=0..lmax
  !i   nrx   :leading dimension of p
  !i   k0    :second dimension of p
  !i   wk    :work array of length nr
  !i   lrl   :if 1s digit = 0, returns p_kl; otherwise returns r^l p_kl
  !i         :if 10s digit nonzero, returns gp; otherwise gp is not touched.
  !o Outputs
  !o   p     :radial part of spherical polynomials P_kL; see Remarks
  !o   gp    :radial derivative of p from l=0..lmax-1 (depending on lrl).
  !r Remarks
  !r   P_kL are polyonomials orthogonal in the following sense:
  !r                                          (4a^2)^k a^l k! (2l+1)!!
  !r    int P_kL G_k'L' = delta_kk'*delta_ll'  ----------------------
  !r                                                    4pi
  !r   and are defined in J. Math. Phys. 39, 3393 (1988).
  !r   Combining eqns 12.7 and 5.19 in that paper, we obtain
  !r    p_kl = a**l / (2a**2)^(k+l) (2l+1)!! / (2k+2l+1)!! phi_kl
  !r    p_0l = a**l
  !r    p_1l = a**l (2*(ar)**2/(2l+3) - 1)
  !r    p_kl = [(2*(ar)**2 - (4k+2l-1))p_k-1,l - 2(k-1)p_k-2,l]
  !r           / (2k+2l+1)
  !u Updates
  !u   22 Aug 01 bug fix for gp when kmax=0
  !u   25 Jan 00 veckl generates gp as well as p.
  ! ----------------------------------------------------------------------
  !     implicit none
  integer :: nr,kmax,lmax,nrx,k0,lrl
  double precision :: r(nrx),wk(nr),rsm,p(nrx,0:k0,0:*), &
       gp(nrx,0:k0,0:*)
  integer :: i,l,k
  double precision :: a,xx,xx2,xx3

  if (kmax < 0 .OR. lmax < 0) return
  if (kmax > k0) call rx('vecpkl: kmax gt k0')
  if (rsm <= 0) call rx('vecpkl: rsm <= 0')
  a = 1d0/rsm

  ! --- Set wk = 2*a**2*r**2 ---
  xx = 2*a*a
  do  6  i = 1, nr
     wk(i) = xx*r(i)**2
6 enddo

  ! --- Do explicitly for k=0,1 ---
  do    l = 0, lmax
     xx = a**l
     do    i = 1, nr
        p(i,0,l) = xx
     enddo
  enddo

  if (kmax > 0) then
     do   l = 0, lmax
        xx = a**l
        xx2 = 1/dble(2*l+3)
        do   i = 1, nr
           p(i,1,l)=xx*(wk(i)*xx2-1d0)
        enddo
     enddo
  endif

  ! --- Recursion for higher k ---
  do    k = 2, kmax
     xx3 = 2*(k-1)
     do    l = 0, lmax
        xx2 = (4*k+2*l-1)
        xx = 1/dble(2*k+2*l+1)
        do    i = 1, nr
           p(i,k,l) = xx*((wk(i)-xx2)*p(i,k-1,l) - xx3*p(i,k-2,l))
        enddo
     enddo
  enddo

  ! --- Radial derivative of p ---
  if (mod(lrl/10,10) /= 0) then

     !  ... Set wk = 2*a**2*r**2
     xx = 2*a*a
     do  16  i = 1, nr
        wk(i) = xx*r(i)
16   enddo

     do    k = 0, kmax
        do    l = 0, lmax-1
           xx2 = dble(2*k+2*l+3)/(a*(2*l+3))
           do    i = 1, nr
              gp(i,k,l) = wk(i)*(p(i,k,l) - xx2*p(i,k,l+1))
           enddo
        enddo
     enddo
  endif

  ! --- Scale by r^l if lrl nonzero ---
  if (mod(lrl,10) == 0) return
  do  50  i = 1, nr
     wk(i) = 1
50 enddo

  do  52  l = 1, lmax

     !   ... gP scales as  r*l gP +  l*r^l-1 P
     if (mod(lrl/10,10) /= 0 .AND. l < lmax) then
        do   k = 0, kmax
           do   i = 1, nr
              gp(i,k,l) = wk(i)*r(i)*gp(i,k,l) + l*wk(i)*p(i,k,l)
           enddo
        enddo
     endif

     do  54  i = 1, nr
        wk(i) = wk(i)*r(i)
54   enddo
     do   k = 0, kmax
        do   i = 1, nr
           p(i,k,l) = p(i,k,l)*wk(i)
        enddo
     enddo
52 enddo
end subroutine vecpkl
end module m_mixrho

subroutine splrho(mode,nsp,nr,nlml,rho1,rho2,rhoc)  !- Overwrite spin pol local rho+,rho- with rho,rho+ - rho-, or reverse
  ! ----------------------------------------------------------------------
  !i Inputs
  !i   mode  :1s digit
  !i         :0 input (rho+,rho-) -> (rho+ + rho-, rho+ - rho-)
  !i         :1 input (rho+ + rho-, rho+ - rho-) -> (rho+,rho-)
  !i         :10s digit
  !i         :1 suppress splitting of rho2
  !i         :2 suppress splitting of rhoc
  !i         :3 suppress both
  !i   nsp   :2 for spin-polarized case, otherwise 1
  !i   nr    :number of radial mesh points
  !i   nlml  :L-cutoff
  !i   rho1  :local true density, tabulated on a radial mesh
  !i   rho2  :local smoothed density, tabulated on a radial mesh
  !i   rhoc  :core density
  !r Remarks
  !u Updates
  ! ----------------------------------------------------------------------
  implicit none
  integer :: mode,nsp,nr,nlml
  double precision :: rho1(nr,nlml,nsp),rho2(nr,nlml,nsp),rhoc(nr,nsp)
  double precision :: fac
  if (nsp == 1) return
  fac = 1
  if (mod(mode,10) /= 0) fac = .5d0
  call dsumdf(nr*nlml,fac,rho1,0,1,rho1(1,1,2),0,1)
  if (mod(mod(mode/10,10),2)  == 0)  call dsumdf(nr*nlml,fac,rho2,0,1,rho2(1,1,2),0,1)
  if (mod(mod(mode/10,10)/2,2) == 0) call dsumdf(nr,fac,rhoc,0,1,rhoc(1,2),0,1)
end subroutine splrho
subroutine lgstar(mode,ng,n,gv,ng0,ips0,cg)  !- Compresses F.T. of a real function, using fact it is hermitian
  ! ----------------------------------------------------------------------
  !i Inputs
  !i   mode  :0, count number of inequivalent points ng0,
  !i             and make ips0.  cg is not used.
  !i         :1  same as mode 0, but also compress cg
  !i         :2  use ips0 to undo compression of cg
  !i   ng    :number of G-vectors
  !i   n     :cg array holds n functions;
  !i   gv    :list of reciprocal lattice vectors G (gvlist.f)
  !o   ng0   :(mode=2) number of inequivalent points
  !i   ips0  :(mode=2) permutation array.
  !i   cg    :list of g-vectors for each of n functions
  !o Outputs
  !o   ng0   :(mode=0,1) number of inequivalent points
  !o   ips0  :(mode=0,1) array of permutation indices.  A negative value
  !o         :signifies the point's hermitian counterpart falls earlier in
  !o         :the list, and points to that element.
  !r Remarks
  !r   Hacked from svgsym, using only one symmetry operation
  !r   to reduce (G,-G) pairs to a single element.
  !u Updates
  ! ----------------------------------------------------------------------
  !     implicit none
  integer :: mode,ng,n,ips0(ng),ng0,iprint
  double precision :: gv(ng,3)
  double complex cg(ng,n)
  integer :: i,i0,i00,irep,k,j,j0,lwarn,m
  double precision :: v(3),df
  if (mode < 0 .OR. mode > 2) call rxi('lgstar, bad mode',mode)
  lwarn = 0
  if (mode == 2) goto 200
  ! --- mode = 0,1 ---
  ng0 = 0
  do  5  i = 1, ng
     ips0(i) = 0
5 enddo
  ! --- Main loop: look for next unclassified vector ---
  i00 = 1
  do  10  irep = 1, ng+1
     i0 = 0
     do  12  i = i00, ng
        i0 = i
        if (ips0(i) == 0) goto 80
12   enddo
     goto 81
80   continue
     !   ... Apply all point ops, find in list, add to phase sum
     ng0 = ng0 + 1
     ips0(i0) = ng0
     if (mode == 1)  then
        do  21  m = 1, n
           cg(ng0,m) = cg(i0,m)
21      enddo
     endif
     do  20  k = 1, 1
        v(1) = gv(i0,1)
        v(2) = gv(i0,2)
        v(3) = gv(i0,3)
        do  22  j = i0+1,ng
           df = (v(1)+gv(j,1))**2+(v(2)+gv(j,2))**2+(v(3)+gv(j,3))**2
           j0 = j
           if (df < 1d-8) goto 70
22      enddo
        !     ... No matching vector here ... should only happen for G=0
        i00 = i0
        goto 10
70      continue
        ips0(j0) = -i0
        if (mode == 1) then
           if (abs(cg(i0,1)-dconjg(cg(j0,1))) > 1d-9) lwarn = lwarn+1
        endif
20   enddo
     i00 = i0
10 enddo
  call rxi('bug in lgstar, irep=',irep)
81 continue
  if (lwarn > 1 .AND. iprint() >= 10) print 345, lwarn
345 format(' lgstar (warning):',i6, ' points not hermitian')
  return
  ! --- mode = 2 ---
200 continue
  ! ... Unpack original points first
  do    m = 1, n
     do    i = ng, 1, -1
        k = ips0(i)
        if (k > 0) cg(i,m) = cg(k,m)
     enddo
  enddo
  ! ... Unpack hermitian points
  do    m = 1, n
     do    i = 1, ng
        k = -ips0(i)
        if (k > 0) cg(i,m) = dconjg(cg(k,m))
     enddo
  enddo
end subroutine lgstar
subroutine dsumdf(n,scal,a1,ofa1,l1,a2,ofa2,l2)! Returns scaled sum and difference of two vectors
  ! ----------------------------------------------------------------
  !i Inputs
  !i   n    :number elements to scale and combine
  !i   scal :scale sum and difference by scal; see Outputs
  !i   a1   :first vector
  !i   ofa1 :offset to first entry in a1
  !i   l1   :skip length in a1
  !i   a2   :second vector
  !i   ofa2 :offset to first entry in a2
  !i   l2   :skip length in a2
  !o Outputs
  !o   a1   :a1 <- scal*(a1+a2)
  !o   a2   :a1 <- scal*(a1-a2)
  ! ----------------------------------------------------------------
  implicit none
  integer :: n,l1,l2,ofa1,ofa2
  double precision :: scal, a1(1), a2(1)
  real(8) ,allocatable :: a_rv(:)   ! --- a1-a2-> temp;  a1+a2 -> a1;  temp -> a2 ---
  allocate(a_rv(n))
  call dcopy ( n , a1 ( 1 + ofa1 ) , l1 , a_rv , 1 )
  call daxpy ( n , - 1d0 , a2 ( 1 + ofa2 ) , l2 , a_rv , 1 )
  call daxpy (n,1d0,a2(1+ofa2),l2,a1(1+ofa1),l1)
  call dcopy ( n , a_rv , 1 , a2 ( 1 + ofa2 ) , l2 )
  deallocate(a_rv)
  if (scal == 1) return
  call dscal(n,scal,a1(1+ofa1),l1)
  call dscal(n,scal,a2(1+ofa2),l1)
end subroutine dsumdf
subroutine dpdump(array,length,ifile)! Binary I/O of an array
  integer:: length,ifile
  double precision :: array(length)
  if (ifile > 0) read(ifile) array
  if (ifile < 0) write(-ifile) array
end subroutine dpdump
logical function lddump(array,length,ifile)! Binary I/O of an array, returning T if I/O without error or EOF
  integer :: length,ifile
  double precision :: array(length),xx,yy
  lddump = .true.
  if (ifile > 0) then
     yy = array(length)
     xx = -1.9283746d0*datan(1d0) !       (some random number)
     array(length) = xx
     read(ifile,end=90,err=91) array
     if (xx /= array(length)) return
     array(length) = yy
     goto 90
90   continue
91   continue
     lddump = .false.
  else
     write(-ifile) array
  endif
end function lddump

