      program dl   ! in canonical ensemble   
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      dimension epxtot(ntotmx),epytot(ntotmx),epztot(ntotmx)
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz

      write(*,*)" Start MC Double layer"        

      open(2,file="in.dat")
      read(2,*)niter
      close(2)

      do iter=1,niter

         call readin    
         istep0=0
         call wrcf(istep0)
         call zero

         write(*,*)" Start equilibration"
         do istep=1,naver
            do ip=1,ntotmx
                call tryit
            enddo  
            if (mod(istep,isave).eq.0) then
              call wrcf(istep) 
              write(*,*)istep," I am equilibrating"
            endif    
         enddo  
     
         call zero
         write(*,*)" Start simulation"    
         do istep=1,nstep
             do ip=1,ntotmx
                call tryit
                call hist
             enddo 
             if (mod(istep,isave).eq.0) call saves(istep)
             if (mod(istep,isave).eq.0) call wrcf(istep)
         enddo  
        
         if (iter.eq.niter) then
               write(*,*)" End of simulation"
         else
               write(*,*)" End of iteration ",iter
               write(*,'(72("="))')
         endif

      enddo
      
      stop   
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine readin    
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      dimension epxtot(ntotmx),epytot(ntotmx),epztot(ntotmx)
      dimension     xn(nspcmx),xq(nspcmx)
      logical       runit
      logical       readit, ovrlap, targc, difdca
      character*2   fion
      character*4   fsalt
      character*9   fgz

c_______ CONSTANTS __________________________________________________
      mag=51328632
      pi=4.0*datan(1.d0) 
      twopi=2*pi
      sq2=dsqrt(2.d0)
      sq3=dsqrt(3.d0)

c_______ READIN PARAMETERS __________________________________________
c_______ CHECK FOR SIZES OF ARRAYS __________________________________
     
      open(2,file="in.dat")
     
c---read in first AND geometry---
      read(2,*)niter
      write(*,*)1
      read(2,*)targc,ntrg,ctrg
      read(2,*)xl,zl
      read(2,*)sigma(1),sigma(2)
      read(2,*)voltg
      read(2,*)
      read(2,*)

c---simulation parameters--
      read(2,*)drmax
      read(2,*)dchgmx
      read(2,*)dzgtrg  
      read(2,*)readit
      read(2,*)dsi
      read(2,*)tstar
      read(2,*)runit
      read(2,*)difdca
      read(2,*)ratchg
      read(2,*)nstep0
      read(2,*)naver 
      read(2,*)isave
     
      read(2,*)
      read(2,*)

c---species---
      read(2,*)nspec
      read(2,*)       
      write(*,'("  nspec = ",i6,"   max = ",i6)')nspec,nspcmx
      if (nspec.gt.nspcmx) stop "nspec too large"
      do ispec=1,nspec 
         read(2,*)eps(ispec),ratmov(ispec),
     :             z(ispec),d(ispec),dca(ispec),fion(ispec)
      if (.not.difdca) dca(ispec)=d(ispec)/2.0
      enddo      
      epsw=eps(1)

      read(2,*)
      read(2,*)
c---grand canonical---
      read(2,*)nsalt
      read(2,*)
      ratgrt=0.0
      ctot=0.0
      do i=1,nsalt
        read(2,*)isalt(i),ctarg(i),chex(i),ratgr(i),fsalt(i)
        ratgrt=ratgrt+ratgr(i)
        ctot=ctot+ctarg(i)
      enddo
      if (ratgrt.ge.1.0) stop "ratgrt>1"
      write(*,*)" Total cation conc = ",ctot

      close(2)

      do ispec=1,nspec
        alpha(ispec) =
     :      -((eps(ispec)-1)/(2*eps(ispec)+1))*(d(ispec)/2.0)**3
      enddo 

c specify constants
      
      tsi=300.0
      eps0=8.8542d-12
      boltz=1.3806d-23
      echg=1.6021917d-19
      beta=1./(boltz*tsi)
      avog=6.02214d23
      tomv=boltz*tsi/echg
      write(*,*)" tomv = ",tomv," Volts"
      
      if (runit) then
       tosi=1.0
       qstar=dsqrt(1.0/tstar)
       write(*,*)"-----------------------------------------"
       write(*,*)" I work in reduced units "
       write(*,*)" T* = ",tstar
       write(*,*)" q* = ",qstar
       write(*,*)" dsi = ",dsi
       write(*,*)"-----------------------------------------"
      else
       dsi=1.d-10
       tosi=1.0/(1.d-10)**3/avog/1000.
       epsi=eps(1)
       qstar=echg*dsqrt(beta/(4*pi*eps0*dsi))
       tstar=1.0/qstar**2
       write(*,*)"-----------------------------------------"
       write(*,*)" I work in real units "
       write(*,*)" T* = ",tstar
       write(*,*)" q* = ",qstar
       write(*,*)" dsi = ",dsi
       write(*,*)"-----------------------------------------"
      endif 

      voltg=voltg/qstar

c calc cell dimensions if targc=.true.

      if (targc) then
        rtrg=ctrg/tosi
        vbulkt=dble(ntrg)/rtrg
        vhalf=vbulkt/2.0
        xl=vhalf**(1./3.)
        zl=xl
      endif

      write(*,*)" zl = ",zl
      write(*,*)" starting xl = ",xl
      write(*,*)"-----------------------------------------"

c calc cation mole fractions

      do ispec=2,nspec
        xn(ispec)=ctarg(ispec-1)/ctot
        write(*,*)" cation mole fraction : ",ispec,xn(ispec)
      enddo  

c redefine ion numbers

      n(2)=ntrg
      write(*,*)" N(",2,") = ",n(2)
      if (nspec.ge.3) then
       ncation=n(2)
       do ispec=3,nspec
         n(ispec)=int(n(2)*xn(ispec)/xn(2))
         write(*,*)" N(",ispec,") = ",n(ispec)
         ncation=ncation+n(ispec)
       enddo
      endif 

c sigma in e unit

      if (.not.runit) then
        sigsi(1)=sigma(1)
        sigsi(2)=sigma(2)
        sigma(1)=(sigma(1)/echg)*dsi**2
        sigma(2)=(sigma(2)/echg)*dsi**2
      endif

c recalculate xl
      
      qwall=(sigma(1)+sigma(2))*xl**2
      write(*,*)" net charge on left wall /e  = ",sigma(1)*xl**2
      write(*,*)" net charge on right wall /e = ",sigma(2)*xl**2
      write(*,*)" net charge on walls /e (starting) = ",qwall
      nwall=int(qwall)
      if (dabs(qwall).gt.1.d-5) then
       xlsq=dble(nwall)/(sigma(1)+sigma(2))
       xl=dsqrt(xlsq)
       write(*,*)" xl = ",xl
       qwalle=(sigma(1)+sigma(2))*xl**2
       write(*,*)" net charge on walls /e = ",qwalle
      endif 

      sigma(1)=sigma(1)*qstar
      sigma(2)=sigma(2)*qstar
      
      write(*,*)" Initial sigma = ",sigma(1)

      qwall=sigma(1)+sigma(2)
      qwall=qwall*xl**2
      write(*,*)" net charge on walls = ",qwall

c calc N(Cl) from condition of electroneutrality

      nchg=nwall
      do ispec=2,nspec
        nchg=nchg+nint(z(ispec)*n(ispec))
      enddo
      n(1)=nchg
      write(*,*)" N(",1,") = ",n(1)
      write(*,*)"-----------------------------------------"
     
      ntry=0
      nimax=0
      ntot=0
      do ispec=1,nspec
        itry(ispec)=n(ispec)
        ntry=ntry+itry(ispec)
        ntot=ntot+n(ispec)
        if (n(ispec).gt.nimax) nimax=n(ispec)
      enddo  
      write(*,'("  nimax = ",i6,"   max = ",i6)')nimax,nionmx
      write(*,'("   ntot = ",i6,"   max = ",i6)')ntot,ntotmx
      write(*,'("   ntry = ",i6,"   max = ",i6)')ntry,ntrymx
      if (nimax.gt.nionmx) stop "nimax too large"
      if (ntot.gt.ntotmx) stop "ntot too large"
      if (ntry.gt.ntrymx) stop "ntry too large"
      write(*,*)"-----------------------------------------"

c constract output filenames      

      open(2,file="fnames.dat")
      do ispec=1,nspec
        write(2,'("gz-",a2,".dat   pz-",a2,".dat")')fion(ispec),
     :               fion(ispec)
      enddo
      close(2)

      open(2,file="fnames.dat")
      do ispec=1,nspec
        read(2,*)fgz(ispec,1),fgz(ispec,2)
      enddo
      close(2)

      if (iter.eq.niter) then
        nstep=nstep0
      else
        nstep=nstep0/10
      endif

      do igc=1,nsalt
          ispec=isalt(igc)
          igcval(igc)=nint(z(ispec))
          iv=igcval(igc)
          ispcin(igc,1)=isalt(igc)
          do icount=2,iv+1
              ispcin(igc,icount)=1
          enddo
          vgc(igc) = 1.0
          do i=1,iv+1
               ispec = ispcin(igc,i)
               ri = d(ispec)/2.0
               vgc(igc)=vgc(igc)*xl**2*2*(zl-dca(ispec))
          enddo
          write(*,*)" vgc = ",igc,vgc(igc)
      enddo 

      ctrgcl=0.0
      do igc=1,nsalt
        iv=igcval(igc)
        ctrgcl=ctrgcl+iv*ctarg(igc)
      enddo  
      write(*,*)" Cl conc = ",ctrgcl
        
      write(*,*)"-----------------------------------------"
      if (iter.eq.1) then
        write(*,'("     iter = ",i2)')iter
        do igc=1,nsalt
          iv=igcval(igc)
          rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv
          chemp(igc)=dlog(rtargi)+chex(igc)
          write(*,*)igc,chex(igc),chemp(igc),iv
        enddo  
      else
        write(*,*)                           
        write(*,'("     iter = ",i2)')iter
        do igc=1,nsalt
          iv=igcval(igc)
          rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv
          ispec=isalt(igc)
          rbulki=cbulk(ispec)/tosi
          rblkcl=cbulk(1)/tosi
          chexo=chex(igc)
          chex(igc)=chemp(igc)-dlog(rbulki*rblkcl**iv)
          chempo=chemp(igc)
          chemp(igc)=dlog(rtargi)+chex(igc)
          write(*,*) igc,ctarg(igc),cbulk(ispec)
          write(*,*) igc,chempo,chemp(igc)
          write(*,*) igc,chexo,chex(igc)
        enddo  
      endif
              
c_______ CALCULATE REDUCED CHARGES __________________________________
c_______ CHECK FOR ELECTRONEUTRALITY ________________________________
c_______ FILL UP DD, D2, AND NWH ARRAYS _____________________________

      charge=qwalle
      ii=0
      do i=1,nspec
            charge=charge+z(i)*n(i)
            q(i)=z(i)*qstar
            ri=d(i)/2.0
            do j=1,nspec
                  dd(i,j)=(d(i)+d(j))/2.0
                  d2(i,j)=dd(i,j)**2
            enddo
            do j=1,itry(i)
                  ii=ii+1
                  nwh(ii)=i
            enddo      
      enddo  
      if (charge.gt.0.001) stop "Cell is no neutral"

      if (readit) then
 
c_______ READIN CONFIGURATION _______________________________________
c_______ FILL UP INDEX ARRAYS AND LOOKUP TABLES _____________________
              
            open (3,file="c0")
            read(3,*)nspec
            read(3,*)
           
            do ispec=1,nspec
                  read(3,*)n(ispec)                      
                  read(3,*)                      
                  read(3,*)
                  do i=1,n(ispec)
                        ii=(ispec-1)*nionmx+i
                        read(3,*)rxi,ryi,rzi
                        rx(ii)=rxi
                        ry(ii)=ryi
                        rz(ii)=rzi
                        indspc(ispec,i)=ii
                        ispcbk(ii)=i
                  enddo      
            read(3,*)
            enddo
            close(3)

      else

c_______ MAKE CONFIGURATION _________________________________________
            call initc

      endif

      call lookup

c_______ Geometry for gz distributions ______________________________ 
      do i=1,nspec

             nzg(i)  = int(2*(zl-dca(i))/dzgtrg)+1
             dzg(i)  = 2*(zl-dca(i))/dble(nzg(i))
             if (nzg(i).gt.nzgmx) then 
               write(*,*) "nzg is too large"
               write(*,'(" nzg = ",i6,"   nzgmx = ",i6)')nzg(i),nzgmx
               stop
             endif  
         
         write(*,'(40("_"))')

      enddo

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine tryit
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      
      ireg=1
      vszam=ranff(mag)
      ratgr2=0.0
      do igc=1,nsalt
        ratgr2=ratgr2+ratgr(igc)
        ratgr1=ratgr2-ratgr(igc)
        if (vszam.gt.ratgr1.and.vszam.lt.ratgr2) then   
           if (ranff(mag).lt.0.5) then
             call create(igc)
           else
             call destr(igc)
           endif
        endif
      enddo
      
      if (vszam.gt.ratgr2.and.vszam.lt.ratgr2+ratchg) then
          call rechg
      else 
          iwh   = int(ranff(mag)*ntry)+1
          ispec = nwh(iwh)
          if (n(ispec).eq.0) goto 1
          i = int(dble(n(ispec))*ranff(mag)) + 1
          call move(ispec,i)
1         continue
      endif

      auc=auc+uc
      aup=aup+up

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine lookup
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      dimension epxtot(ntotmx),epytot(ntotmx),epztot(ntotmx)
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
     
      do ispec=1,nspec
         do i=1,n(ispec)
             ii=indspc(ispec,i)
             rxi=rx(ii)
             ryi=ry(ii)
             rzi=rz(ii)
      
             do jspec=1,nspec
                  do j=1,n(jspec)
                      jj=indspc(jspec,j)
                      if (ii.ne.jj) then
                         rxj=rx(jj)
                         ryj=ry(jj)
                         rzj=rz(jj)
                         rxij=rxi-rxj
                         ryij=ryi-ryj
                         rzij=rzi-rzj
                         rxij = rxij-dnint(rxij/xl)*xl
                         ryij = ryij-dnint(ryij/xl)*xl
                         rijsq= rxij*rxij+ryij*ryij+rzij*rzij
                         if (rijsq.lt.d2(ispec,jspec)) stop 
     :                      "lookup: Overlap in initial configuration"
                         rii(ii,jj)=dsqrt(rijsq)
                      else
                         rii(ii,jj)=0.0          
                      endif
                  enddo
             enddo

         enddo
      enddo
      
      uc=0.0

      do ispec=1,nspec
        qi=q(ispec)
        di=d(ispec)
        ri=di/2.0
        do i=1,n(ispec)
          ii=indspc(ispec,i)
          do jspec=1,nspec
            ddi=d2(ispec,jspec)
            qj=q(jspec)
            do j=1,n(jspec)
               if (ispec.eq.jspec.and.i.eq.j) goto 1
               jj  = indspc(jspec,j)
               rij=rii(ii,jj)         
               uiiold = qi*qj/(epsw*rij)
               uc=uc + uiiold   
               zi = dabs(rz(ii)-rz(jj))
               ushi = qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl))/epsw
               uc = uc + ushi
c               write(*,*)uc,uiiold,ushi
1              continue               
            enddo
      enddo
      enddo
      enddo
      uc=uc/2.0
      
      uc = uc + qi*sigma(1)*(-twopi*(rz(ii)+zl))/epsw 
     :        + qi*sigma(2)*(-twopi*(zl-rz(ii)))/epsw

      call fullit

      up=0.0

      do jspec = 1,nspec
         qj=q(jspec)
         do j = 1,n(jspec)
           jj = indspc(jspec,j)
           rxj=rx(jj)
           ryj=ry(jj)
           rzj=rz(jj)
           phitot=0.0
           do kspec = 1,nspec
             do k = 1,n(kspec)
               kk = indspc(kspec,k)
               rxk=rx(kk)
               ryk=ry(kk)
               rzk=rz(kk)
               if (jj.ne.kk) then
                 rij=rii(jj,kk)
                 rxij = rxj-rxk     
                 ryij = ryj-ryk   
                 rzij = rzj-rzk    
                 rxij = rxij-dnint(rxij/xl)*xl
                 ryij = ryij-dnint(ryij/xl)*xl
c                 rzij = rzij-dnint(rzij/xl)*xl
                 phitot=phitot+(dipx(kk)*rxij
     :                         +dipy(kk)*ryij
     :                         +dipz(kk)*rzij)/rij**3
c                 write(*,*)jj,kk,rij,phitot
               endif
             enddo    
           enddo 
           up = up+qj*phitot
           up = up+twopi*(-sigma(1)+sigma(2))*dipz(jj)
c           write(*,*)jj,upn,qj*phitot
         enddo 
      enddo     
      up=up/2.0

      write(*,*)" Total Coulomb = ",uc,uc/dble(ntot)
      write(*,*)" Total Polarization = ",up,up/dble(ntot)

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

       function ranff ( seed )
       implicit doubleprecision (a-h,o-z)
       integer a14,p31,seed,b15,b16,xhi,xalo,leftlo,fhi,kf15
       data a14/16807/,b15/32768/,b16/65536/,p31/2147483647/
       xhi=seed/b16
       xalo=(seed-xhi*b16)*a14
       leftlo =xalo/b16
       fhi=xhi*a14+leftlo
       kf15=fhi/b15
       seed=(((xalo-leftlo*b16 )-p31) +(fhi-kf15*b15)*b16) +kf15
       if (seed .lt. 0) seed=seed+p31
       xhi=seed/b16
       xdr=(float(xhi)*65536.0d0)+float(seed-xhi*b16)
       ranff=xdr*4.6566128752457969241d-10
       seed=xdr
       return
       end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine hist 
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      
      athist=athist+1.0
      
      do ispec=1,nspec
        adip0=0.0
        adip20=0.0
        do i=1,n(ispec)
          ii=indspc(ispec,i)
          dip0=dipx(ii)**2+dipy(ii)**2+dipz(ii)**2
          adip0=adip0+dsqrt(dip0)
          adip20=adip20+dip0
        enddo
        adip0=adip0/dble(n(ispec))
        adip20=adip20/dble(n(ispec))
        adip(ispec)=adip(ispec)+adip0
        adipsq(ispec)=adipsq(ispec)+adip20
      enddo

      do ispec=1,nspec
        do i=1,n(ispec)
          ii=indspc(ispec,i)
          rzi=rz(ii)
          k=int((rzi+zl-dca(ispec))/dzg(ispec))+1
          gz(ispec,k)=gz(ispec,k)+1.0
          pol(ispec,k)=pol(ispec,k)+dipz(ii) 
          if (dabs(rzi).lt.zl/2.0) an(ispec)=an(ispec)+1.0
        enddo    
      enddo
  
      asigma(1)=asigma(1)+sigma(1)
      asigma(2)=asigma(2)+sigma(2)

      return
      end
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine initc
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
      
      write(*,*)" Making initial configuration"

      do ispec=1,nspec
       do i=1,nionmx
        indspc(ispec,i)=0
       enddo
      enddo
      
11    continue

      ii=0
      do 10 ispec=1,nspec

      ntrg=n(ispec)
      ni=0

30    ni=ni+1
20    rxinew = (ranff(mag)-0.5)*xl
      ryinew = (ranff(mag)-0.5)*xl
      rzinew = (2*ranff(mag)-1)*(zl-dca(ispec))

      ii=ii+1
      if (ii.gt.1000000) goto 11      
      
      call overlp(ispec,ni,rxinew,ryinew,rzinew,ovrlap)
      if (ovrlap) then
       goto 20
      else
       ii=(ispec-1)*nionmx+ni
       rx(ii)=rxinew
       ry(ii)=ryinew
       rz(ii)=rzinew
       indspc(ispec,ni)=ii
       ispcbk(ii)=ni
c       write(*,*)ii,ispec,ni,rxinew,ryinew,rzinew
      endif 

      if (ni.lt.ntrg) goto 30

10    continue

      write(*,*)" Initial configuration made"

      return
      end
       
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine overlp(ispec,i,rxi,ryi,rzi,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
      
      do jspec=1,ispec
            if (jspec.eq.ispec) then
                  nj=i-1
            else      
                  nj=n(jspec)
            endif
            
            if (nj.gt.0) then
               dda=d2(ispec,jspec)
               do j=1,nj
                     jj=indspc(jspec,j)
                     rxa=rx(jj)
                     rya=ry(jj)
                     rza=rz(jj)
                     rxia=rxi-rxa
                     ryia=ryi-rya
                     rzia=rzi-rza
                     rxia=rxia-dnint(rxia/xl)*xl
                     ryia=ryia-dnint(ryia/xl)*xl
                     riasq=rxia*rxia+ryia*ryia+rzia*rzia
                     if (riasq.lt.dda) then
                           ovrlap=.true.
                           return
                     endif
               enddo     
            endif
      enddo               

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine wrcf(istep)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
c      write(*,*)' begin wrcf'
      open(2,file="c.dat")  
      write(2,*) nspec,istep                             
      write(2,*)
      do ispec=1,nspec
            write(2,*)n(ispec),d(ispec),dca(ispec)
            write(2,*)nzg(ispec),dzg(ispec)
            write(2,*)
            do i=1,n(ispec)
                  ii=indspc(ispec,i)
                  write(2,*)rx(ii),ry(ii),rz(ii)
            enddo      
      write(2,*)
      enddo
      close(2)

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine saves(istep)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      character*9   fngz,fnpz
      dimension     cz(nspcmx,nzgmx),dinteg(nspcmx),gamma1(nspcmx),
     :              gamma2(nspcmx),pz(nspcmx,nzgmx),pinteg(nspcmx) 
     
      open(2,file="out.dat")
      write(*,'(72("#"))')
      write(*,*)
      write(2,*)
      write(*,*)"--- Save ----, Iteration: ",iter," ---- Step: ",
     :          istep,"------"
      write(*,'(72("_"))')
      write(2,'(72("#"))')
      write(2,*)"--- Save ----, Iteration: ",iter," ---- Step: ",
     :          istep,"------"
      write(2,'(72("_"))')
        
      write(2,'(" zl [A]     = ",f10.5)')zl
      write(2,'(" xl [A]     = ",f10.5)')xl
      write(2,'(72("-"))')
      write(*,'(72("-"))')

      if (runit) then
        write(*,*)"  Left electrode charge = ",sigma(1) 
        write(*,*)" Right electrode charge = ",sigma(2) 
        write(2,*)"  Left electrode charge = ",sigma(1) 
        write(2,*)" Right electrode charge = ",sigma(2) 
      else
        write(*,*)"  Left electrode charge = ",sigma(1) 
        write(*,*)" Right electrode charge = ",sigma(2) 
        write(2,*)"  Left electrode charge = ",sigma(1) 
        write(2,*)" Right electrode charge = ",sigma(2) 
        write(2,'(72("-"))')
        write(*,'(72("-"))')
c        write(*,*)"  Left electrode charge [C/m2] = ",sigsi(1)
c        write(*,*)" Right electrode charge [C/m2] = ",sigsi(2) 
c        write(2,*)"  Left electrode charge [C/m2] = ",sigsi(1)
c        write(2,*)" Right electrode charge [C/m2] = ",sigsi(2) 
      endif
      write(2,'(72("-"))')
      write(*,'(72("-"))')

      write(2,*)
      do igr=1,nsalt
        write(2,'(a4,": ctarg = ",f6.3,", chemp = ",f10.5,
     :       ", excess = ",f10.6)')
     :       fsalt(igr),ctarg(igr),chemp(igr),chex(igr) 
      enddo 
      write(*,'(12x,"   N   ","  q [e]  ","  d      "," dca" )')
      write(2,'(12x,"   N   ","  q [e]  ","  d      "," dca" )')
      do i=1,nspec
        write(2,'(a2," ion: ",1x,i7,2x,f5.2,2x,2(2x,f8.4))')
     :      fion(i),n(i),z(i),d(i),dca(i)
        write(*,'(a2," ion: ",1x,i7,2x,f5.2,2x,2(2x,f8.4))')
     :      fion(i),n(i),z(i),d(i),dca(i)
      enddo
      
c ##### Calculate and write out concentration profiles ###############
        
        do ispec=1,nspec
          fngz  = fgz(ispec,1)
          fnpz  = fgz(ispec,2)
          ri=d(ispec)/2.0
          open(12,file=fngz)
          open(22,file=fnpz)
           do i=1,nzg(ispec)
             vslab = xl**2*dzg(ispec)
             cz(ispec,i) = tosi*gz(ispec,i)/athist/vslab
             pz(ispec,i) = pol(ispec,i)/athist/vslab
             rgz=-zl+ri+(i-0.5)*dzg(ispec)
             write(12,*)rgz,cz(ispec,i),dzg(ispec)
             write(22,*)rgz,pz(ispec,i)
           enddo   
          close(12)
          close(22)
        enddo

c ##### Calculate bulk concentration and #############################
c ##### number of ions in the filter #################################

      write(*,'(72("_"))')
      write(2,'(72("_"))')
      vbulk=xl**2*zl
      do ispec=1,nspec
        cbulk(ispec)=tosi*an(ispec)/athist/vbulk
        write(*,'(" Bulk concentration of ",a2," = ",f10.5)')
     '      fion(ispec),cbulk(ispec)
        write(2,'(" Bulk concentration of ",a2," = ",f10.5)')
     '      fion(ispec),cbulk(ispec)
      enddo
      write(*,'(72("_"))')
      write(2,'(72("_"))')
      write(*,*)" Average electrode charge 1 = ",asigma(1)/athist
      write(*,*)" Average electrode charge 2 = ",asigma(2)/athist
      
c #### Salt data ####################################################

        dii=dble(istep*ntotmx)
        ucav=auc/dii
        upav=aup/dii
      
        write(*,'(72("_"))')
        write(2,'(72("_"))')
        write(*,*)" Average energies = ",ucav/dble(ntot),
     :            upav/dble(ntot),
     :           (ucav+upav)/dble(ntot)
        write(*,*)" Energies = ",uc/dble(ntot),up/dble(ntot)
        write(2,*)" Energies = ",ucav/dble(ntot),upav/dble(ntot),
     :           (ucav+upav)/dble(ntot)
        write(*,'(72("_"))')
        write(2,'(72("_"))')

c       call lookup

        write(*,'(72("_"))')
        write(2,'(72("_"))')
        do ispec=1,nspec
        write(*,*)fion(ispec)," Dipole moment = ",adip(ispec)/athist,
     :           dsqrt(adipsq(ispec)/athist)
        enddo
        write(*,'(72("_"))')
        write(2,'(72("_"))')
        write(*,'(72("_"))')
        write(2,'(72("_"))')

        write(2,'(10x,1x," mu OLD",2x,3x,"mu NEW",
     :        5x,"Excess OLD",3x,"excess NEW ")')
        write(*,'(10x,1x," mu OLD",2x,3x,"mu NEW",
     :        5x,"Excess OLD",3x,"excess NEW ")')
        do igc=1,nsalt
          iv=igcval(igc)
          rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv
          ispec=isalt(igc)
          rbulki=cbulk(ispec)/tosi
          rblkcl=cbulk(1)/tosi
          chexnw=chemp(igc)-dlog(rbulki*rblkcl**iv)
          chempn=dlog(rtargi)+chexnw
          write(2,'(4x,a4,": ",4(1x,f10.6))')fsalt(igc),chemp(igc),
     :        chempn,chex(igc),chexnw     
          write(*,'(4x,a4,": ",4(1x,f10.6))')fsalt(igc),chemp(igc),
     :        chempn,chex(igc),chexnw     
        enddo  

c ##### Acceptance ratios ############################################

      write(*,'(72("_"))')
      write(2,'(72("_"))')
      write(*,*)" Electrode charge exchange: ",achg(2)/achg(1)
      write(2,*)" Electrode charge exchange: ",achg(2)/achg(1)
      
      write(*,'(3x,"    Move ","    Jump ")')
      write(2,'(3x,"    Move ","    Jump ")')
      
      do i=1,nspec
            write(*,'(a2,":",2(2x,f7.5))')
     :      fion(i), amove(2,i)/amove(1,i),ajump(2,i)/ajump(1,i)
            write(2,'(a2,":",2(2x,f7.5))')
     :      fion(i), amove(2,i)/amove(1,i),ajump(2,i)/ajump(1,i)
      enddo

      do icount=1,nsalt         
          write(*,'(72("_"))')
          write(*,'("  Create ",a4," : "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),
     :      acreat(icount,2),acreat(icount,1),
     :      acreat(icount,2)/acreat(icount,1)
          write(*,'(" Destroy ",a4," : "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),
     :      adest(icount,2),adest(icount,1),
     :      adest(icount,2)/adest(icount,1)
          write(2,'(72("_"))')
          write(2,'("  Create ",a4," : "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),
     :      acreat(icount,2),acreat(icount,1),
     :      acreat(icount,2)/acreat(icount,1)
          write(2,'(" Destroy ",a4," : "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),
     :      adest(icount,2),adest(icount,1),
     :      adest(icount,2)/adest(icount,1)
      enddo
      
      write(2,'(72("_"))')
      write(*,'(72("_"))')

      pot1=0.0
      potp1=0.0
      do ispec=1,nspec
        nzg2=nzg(ispec)/2
        dinteg(ispec)=0.0
        pinteg(ispec)=0.0
        gamma1(ispec)=0.0
        do i=1,nzg2
          zg1=dca(ispec)+(i-1)*dzg(ispec)
          zg2=zg1+dzg(ispec)
          dinteg(ispec)=dinteg(ispec)
     :        +(cz(ispec,i)/tosi)*(zg2**2-zg1**2)/2.0
          pinteg(ispec)=pinteg(ispec)
     :        +pz(ispec,i)* dzg(ispec)
          gamma1(ispec)=gamma1(ispec)
     :                  +(cz(ispec,i)-cbulk(ispec))*dzg(ispec)
        enddo 
        pot1=pot1+z(ispec)*dinteg(ispec)
        potp1=potp1+pinteg(ispec)
        gamma1(ispec)=(gamma1(ispec)/tosi)*xl**2
      enddo
      pot1=-pot1*2*twopi*qstar**2
      potp1=-potp1*2*twopi
      write(*,*)"  Ion potential at left wall /kT  = ",pot1
      write(2,*)"  Ion potential at left wall /kT  = ",pot1
      write(*,*)" Pol. potential at left wall /kT  = ",potp1
      write(2,*)" Pol. potential at left wall /kT  = ",potp1
      write(*,*)" Tot. potential at left wall /kT  = ",potp1+pot1
      write(2,*)" Tot. potential at left wall /kT  = ",potp1+pot1
        
      pot2=0.0
      potp2=0.0
      do ispec=1,nspec
        nzg2=nzg(ispec)/2
        dinteg(ispec)=0.0
        pinteg(ispec)=0.0
        gamma2(ispec)=0.0
        do i=1,nzg2
          zg1=dca(ispec)+(i-1)*dzg(ispec)
          zg2=zg1+dzg(ispec)
          ii=nzg(ispec)-i+1
          dinteg(ispec)=dinteg(ispec)
     :        +(cz(ispec,ii)/tosi)*(zg2**2-zg1**2)/2.0
          pinteg(ispec)=pinteg(ispec)
     :        -pz(ispec,ii)*dzg(ispec)
          gamma2(ispec)=gamma2(ispec)
     :                  +(cz(ispec,ii)-cbulk(ispec))*dzg(ispec)
        enddo 
        pot2=pot2+z(ispec)*dinteg(ispec)
        potp2=potp2+pinteg(ispec)
        gamma2(ispec)=(gamma2(ispec)/tosi)*xl**2
      enddo
      pot2=-pot2*2*twopi*qstar**2
      potp2=-potp2*2*twopi
      write(*,*)"  Ion potential at right wall /kT  = ",pot2
      write(2,*)"  Ion potential at right wall /kT  = ",pot2
      write(*,*)" Pol. potential at right wall /kT  = ",potp2
      write(2,*)" Pol. potential at right wall /kT  = ",potp2
      write(*,*)" Tot. potential at right wall /kT  = ",potp2+pot2
      write(2,*)" Tot. potential at right wall /kT  = ",potp2+pot2
      
      if (.not.runit) then
        write(2,'(72("-"))')
        write(*,'(72("-"))')
        write(*,*)" Potential at left wall /mV  = ",pot1*tomv*1000
        write(2,*)" Potential at left wall /mV  = ",pot1*tomv*1000
        write(*,*)" Potential at right wall /mV = ",pot2*tomv*1000
        write(2,*)" Potential at right wall /mV = ",pot2*tomv*1000
      endif

      write(2,'(72("_"))')
      write(*,'(72("_"))')
      
      c1=sigma(1)/(tstar*pot1)
      c2=sigma(2)/(tstar*pot2)
      write(*,*)"  Capacitance of left DL = ",c1
      write(2,*)"  Capacitance of left DL = ",c1
      write(*,*)" Capacitance of right DL = ",c2
      write(2,*)" Capacitance of right DL = ",c2      

      if (.not.runit) then
        write(2,'(72("-"))')
        write(*,'(72("-"))')
        write(*,*)"  Capacitance of left DL [F/m2] = ",
     :            sigsi(1)/(tomv*pot1)
        write(2,*)"  Capacitance of left DL [F/m2] = ",
     :            sigsi(1)/(tomv*pot1)
        write(*,*)" Capacitance of right DL [F/m2] = ",
     :            sigsi(2)/(tomv*pot2)
        write(2,*)" Capacitance of right DL [F/m2]= ",
     :            sigsi(2)/(tomv*pot2)
      endif
      
      write(2,'(72("_"))')
      write(*,'(72("_"))')
   
      chgn1=0.0
      chgn2=0.0
      do ispec=1,nspec
        write(*,*)" Adsorption of ion ",ispec,
     :            " on left wall  = ",
     :            gamma1(ispec)
        write(*,*)" Adsorption of ion ",ispec,
     :            " on right wall = ",
     :            gamma2(ispec)
        chgn1=chgn1+z(ispec)*gamma1(ispec)
        chgn2=chgn2+z(ispec)*gamma2(ispec)
        write(2,*)" Adsorption of ion ",ispec,
     :            " on left wall  = ",
     :            gamma1(ispec)
        write(2,*)" Adsorption of ion ",ispec,
     :            " on right wall = ",
     :            gamma2(ispec)
      enddo

      write(2,'(72("-"))')
      write(*,'(72("-"))')

      chgn1t=(sigma(1)/qstar)*xl**2
      chgn2t=(sigma(2)/qstar)*xl**2
      write(*,*)"  net charge in left DL = ",chgn1,chgn1t
      write(*,*)" net charge in right DL = ",chgn2,chgn2t
      write(2,*)"  net charge in left DL = ",chgn1,chgn1t
      write(2,*)" net charge in right DL = ",chgn2,chgn2t
     
      write(2,'(72("-"))')
      write(*,'(72("-"))')

      do ispec=1,nspec
        write(*,*)" Relative adsorption of ion ",ispec,
     :            " on left wall  = ",
     :            gamma1(ispec)/dabs(chgn1t)
        write(*,*)" Relative adsorption of ion ",ispec,
     :            " on right wall = ",
     :            gamma2(ispec)/dabs(chgn2t)
        write(2,*)" Relative adsorption of ion ",ispec,
     :            " on left wall  = ",
     :            gamma1(ispec)/dabs(chgn1t)
        write(2,*)" Relative adsorption of ion ",ispec,
     :            " on right wall = ",
     :            gamma2(ispec)/dabs(chgn2t)
      enddo

      write(2,'(72("_"))')
      write(*,'(72("_"))')

      return
      end
      
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      subroutine zero
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz

c_______ ZERO OUT STAFF _____________________________________________ 

      auc=0.0
      aup=0.0

      asigma(1)=0.0
      asigma(2)=0.0

      do ispec=1,nspec
       
       adip(ispec)=0.0
       adipsq(ispec)=0.0

       amove(1,ispec)=0.0
       amove(2,ispec)=0.0
       ajump(1,ispec)=0.0
       ajump(2,ispec)=0.0
     
       an(ispec)=0.0
     
       do i=1,nzg(ispec)
          gz(ispec,i) = 0.0
          pol(ispec,i) = 0.0
       enddo  
     
       do icount=1,3
         acreat(icount,1)=0.0
         acreat(icount,2)=0.0
         adest(icount,1)=0.0
         adest(icount,2)=0.0
       enddo 
      
      enddo 
      
      athist = 0.0
     
      return
      end
      
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine move(ispec,i)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
c      write(*,*)"move"	

      ii=indspc(ispec,i)
      rxiold = rx(ii)
      ryiold = ry(ii)
      rziold = rz(ii)
      rzold  = dabs(rziold)
      di=d(ispec)
      ri=di/2.0
      
      if (ranff(mag).lt.ratmov(ispec)) then
      
           amove(1,ispec) = amove(1,ispec)+1
           rxinew = rxiold + (2.0*ranff(mag)-1.0) * drmax
           ryinew = ryiold + (2.0*ranff(mag)-1.0) * drmax
           rzinew = rziold + (2.0*ranff(mag)-1.0) * drmax
           rxinew = rxinew-dnint(rxinew/xl)*xl
           ryinew = ryinew-dnint(ryinew/xl)*xl
           drznew = dabs(rzinew)
           if (drznew.gt.zl-dca(ispec)) return
           imove=1

      else
      
           ajump(1,ispec) = ajump(1,ispec)+1
c  %%%%%%%%% Generate new position in box %%%%%%%%%%%%%%%%%%%%%
           rxinew = (ranff(mag)-0.5)*xl
           ryinew = (ranff(mag)-0.5)*xl
           rzinew = (2*ranff(mag)-1)*(zl-dca(ispec))
           drznew = dabs(rzinew)
           imove=0 

      endif
     
      call ergi(ispec,i,rxinew,ryinew,rzinew,ucold,ucnew,upn,ovrlap)
      if (ovrlap) then
c        write(*,*)"end move"
        return 
      endif  
      deltu=ucnew-ucold+upn-up

      if ((deltu.le.0.0).or.(dexp(-deltu).gt.ranff(mag))) then
            if (imove.eq.1) then
               amove(2,ispec) = amove(2,ispec) + 1.0
            else
               ajump(2,ispec) = ajump(2,ispec) + 1.0
            endif   
            rx(ii) = rxinew
            ry(ii) = ryinew
            rz(ii) = rzinew
            uc=uc+ucnew-ucold
            up=upn
            do jspec=1,nspec
             do j=1,n(jspec)
              jj=indspc(jspec,j)
              dipx(jj)=dipxn(jj)
              dipy(jj)=dipyn(jj)
              dipz(jj)=dipzn(jj)
              dipx0(jj)=dipx0n(jj)
              dipy0(jj)=dipy0n(jj)
              dipz0(jj)=dipz0n(jj)
              eqx(ii,jj)=eqxn(ii,jj)
              eqy(ii,jj)=eqyn(ii,jj)
              eqz(ii,jj)=eqzn(ii,jj)
              eqx(jj,ii)=eqxn(jj,ii)
              eqy(jj,ii)=eqyn(jj,ii)
              eqz(jj,ii)=eqzn(jj,ii)
             enddo
            enddo 
            do jspec=1,nspec
             do j=1,n(jspec)
               jj=indspc(jspec,j)
               rii(ii,jj)=riinew(ispec,jj)
               rii(jj,ii)=riinew(ispec,jj)
             enddo
            enddo
      endif

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine ergi(ispec,i,rxi,ryi,rzi,ucold,ucnew,upn,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
     
      ucold=0.0
      ucnew=0.0
      
      qi=q(ispec)
      di=d(ispec)
      ri=di/2.0
      ii=indspc(ispec,i)
      
c    ###################################################################
c_______ SUM OVER OTHER FREE CHARGES ________________________________
c    ###################################################################
      
      do jspec=1,nspec
            ddi=d2(ispec,jspec)
            qj=q(jspec)
            nj=n(jspec)
            do j=1,nj
               if (ispec.eq.jspec.and.i.eq.j) goto 1
               jj  = indspc(jspec,j)
               rxj = rx(jj)
               ryj = ry(jj)
               rzj = rz(jj)
c_______ DISTANCE IS CALCULATED FOR IONS IN NEW CONFIGURATION ______ 
c_______ STORED IN RIINEW __________________________________________ 
               rxij = rxi-rxj
               ryij = ryi-ryj
               rzij = rzi-rzj
               rxij = rxij-dnint(rxij/xl)*xl
               ryij = ryij-dnint(ryij/xl)*xl
               rijsq= rxij*rxij+ryij*ryij+rzij*rzij
               if (rijsq.lt.ddi) then
                     ovrlap = .true.
                     return
               endif 
               rij= dsqrt(rijsq)
               riinew(ispec,jj)=rij
               uiinew = qi*qj/(epsw*rij)
               ucnew=ucnew + uiinew
               zi = dabs(rzij)
               ucnew = ucnew 
     :            + qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl))/epsw
c_______ DISTANCE IS FROM LOOKUP TABLE ______________________________
c_______ FOR OLD CONFIGURATION ______________________________________
               rij=rii(ii,jj)         
               uiiold = qi*qj/(epsw*rij)   
               ucold=ucold + uiiold   
               zi = dabs(rz(ii)-rzj)
               ucold = ucold 
     :             + qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl))/epsw
1              continue               
            enddo
      enddo         
      
      ucold = ucold + qi*sigma(1)*(-twopi*(rz(ii)+zl))/epsw
     :          + qi*sigma(2)*(-twopi*(zl-rz(ii)))/epsw
      ucnew = ucnew + qi*sigma(1)*(-twopi*(rzi+zl))/epsw 
     :          + qi*sigma(2)*(-twopi*(zl-rzi))/epsw
    
      upn=0.0

c      up=0.0
      ialpha=1
      if (ialpha.eq.0) return

      call iterate(ispec,i,rxi,ryi,rzi)

      do jspec = 1,nspec
         qj=q(jspec)
         do j = 1,n(jspec)
           jj = indspc(jspec,j)
           if (jj.eq.ii) then
             rxj=rxi
             ryj=ryi
             rzj=rzi
           else  
             rxj=rx(jj)
             ryj=ry(jj)
             rzj=rz(jj)
           endif
           phitot=0.0
           do kspec = 1,nspec
             do k = 1,n(kspec)
               kk = indspc(kspec,k)
               if (kk.eq.ii) then
                 rxk=rxi
                 ryk=ryi
                 rzk=rzi
               else  
                 rxk=rx(kk)
                 ryk=ry(kk)
                 rzk=rz(kk)
               endif
               if (jj.ne.kk) then
                 if (jj.eq.ii) then
                    rij=riinew(ispec,kk)
                 elseif (kk.eq.ii) then
                    rij=riinew(ispec,jj)
                 else
                    rij=rii(jj,kk)
                 endif
                 rxij = rxj-rxk     
                 ryij = ryj-ryk   
                 rzij = rzj-rzk    
                 rxij = rxij-dnint(rxij/xl)*xl
                 ryij = ryij-dnint(ryij/xl)*xl
c                 rzij = rzij-dnint(rzij/xl)*xl
                 phitot=phitot+(dipxn(kk)*rxij
     :                         +dipyn(kk)*ryij
     :                         +dipzn(kk)*rzij)/rij**3
c                 write(*,*)jj,kk,rij,phitot
               endif
             enddo    
           enddo 
           upn = upn+qj*phitot
           upn = upn+twopi*(-sigma(1)+sigma(2))*dipz(jj)
c           write(*,*)jj,upn,qj*phitot
         enddo 
      enddo     
      upn=upn/2.0

c      write(*,*)" end erg"
      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine rechg
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
c      write(*,*)"move"	

      sigold(1)=sigma(1)
      sigold(2)=sigma(2)
      achg(1) = achg(1)+1.0
     
      dsigma = (2.0*ranff(mag)-1.0) * dchgmx
      sigma(1) = sigold(1) + dsigma 
      sigma(2) = sigold(2) - dsigma 
     
      call erchg(deltu)

      if ((deltu.le.0.0).or.(dexp(-deltu).gt.ranff(mag))) then
            achg(2) = achg(2) + 1.0
      else
            sigma(1)=sigold(1)
            sigma(2)=sigold(2)
      endif

      return
      end

c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine erchg(uii)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
     
      uii=0.0
      
c    ###################################################################
c_______ SUM OVER OTHER FREE CHARGES ________________________________
c    ###################################################################
      
      do jspec=1,nspec
            qj=q(jspec)
            nj=n(jspec)
            do j=1,nj
               jj  = indspc(jspec,j)
               rzj = rz(jj)
               uii = uii - qj*sigold(1)*(-twopi*(rzj+zl)) 
     :               - qj*sigold(2)*(-twopi*(zl-rzj))
               uii = uii + qj*sigma(1)*(-twopi*(rzj+zl)) 
     :               + qj*sigma(2)*(-twopi*(zl-rzj))
            enddo
      enddo         
    
      uii = uii - (-twopi*2*zl*sigold(1)*sigold(2)*xl**2)
      uii = uii + (-twopi*2*zl*sigma(1)*sigma(2)*xl**2)
      
      uii = uii + voltg*sigold(1)*xl**2
      uii = uii - voltg*sigma(1)*xl**2
      
      return
      end

c    ###################################################################
c    ###################################################################
c    ###################################################################
c    ###################################################################

      subroutine create(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
      
      acreat(igc,1) = acreat(igc,1)+1.0

      iv=igcval(igc)
      
      ntrial = ntot + iv+1
      if ( ntrial .ge. ntotmx ) stop 'maximum number of atoms in box'
   
      ispec1=ispcin(igc,1)
      nterm=n(ispec1)+1
      do icount=2,iv+1
        ispec=ispcin(igc,icount)
        nterm=nterm*(n(ispec)+icount-1)
      enddo  

      do icount=1,iv+1
         ispec  = ispcin(igc,icount)
         rxinew = (ranff(mag)-0.5)*xl
         ryinew = (ranff(mag)-0.5)*xl
         rzinew = (2*ranff(mag)-1)*(zl-dca(ispec))
         drznew = dabs(rzinew)
         rznew(icount) = rzinew                
         rxnew(icount) = rxinew            
         rynew(icount) = ryinew           
      enddo
    
      call encrea(igc,deltu,ovrlap)
      if (ovrlap) then
        return
      endif
      
      bfcr = dexp(-deltu+chemp(igc))*vgc(igc)/dble(nterm)
      
      if (bfcr.gt.ranff(mag)) then
         call add(igc)
         acreat(igc,2) = acreat(igc,2)+1.0
      endif

      return
      end

c    ###################################################################
c    ###################################################################
c    ###################################################################
c    ###################################################################
      
      subroutine destr(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      
      adest(igc,1) = adest(igc,1) + 1.0
      
      iv=igcval(igc)
      
      ispec1=ispcin(igc,1)
      nterm=n(ispec1)
      if (nterm.eq.0) return
      do icount=2,iv+1
        ispec=ispcin(igc,icount)
        nterm=nterm*(n(ispec)-icount+2)
      enddo  
      
      do icount=1,iv+1
        ispec  = ispcin(igc,icount)
1       continue
        i   = int(n(ispec)*ranff(mag)) + 1
        ii  = indspc(ispec,i)                 
        if (icount.gt.2) then
          do icold=2,icount-1
             ispeco = ispcin(igc,icold)
             iold   = iout(icold)
             iiold  = indspc(ispeco,iold)
             if (ii.eq.iiold) goto 1
          enddo
        endif  
        iout(icount)  = i
      enddo
     
      if (iv.gt.1) then
           if (iout(2).lt.iout(3)) then
             ipuf=iout(2)
             iout(2)=iout(3)
             iout(3)=ipuf
           endif  
           if (iv.eq.3) then
               if (iout(3).lt.iout(4)) then
                 ipuf=iout(4)
                 iout(4)=iout(3)
                 iout(3)=ipuf
               endif  
               if (iout(2).lt.iout(3)) then
                 ipuf=iout(2)
                 iout(2)=iout(3)
                 iout(3)=ipuf
               endif  
           endif        
      endif

      do icount=1,iv+1
        ispec  = ispcin(igc,icount)
        i      = iout(icount)
        ii     = indspc(ispec,i)                 
        rxnew(icount) = rx(ii)
        rynew(icount) = ry(ii)
        rznew(icount) = rz(ii)
      enddo
      
      call endest(igc,deltu,ovrlap)
      
      bfdr = dexp(-deltu-chemp(igc))*dble(nterm)/vgc(igc)
         
      if (bfdr.gt.ranff(mag)) then
         call remove(igc)
         adest(igc,2) = adest(igc,2)+1.0
      endif
      
      return
      end

c###################################################################
c###################################################################
c###################################################################

      subroutine remove(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz

      iv=igcval(igc)
  
      do icount=1,iv+1
        ispec=ispcin(igc,icount)  
        i=iout(icount)
        ii = indspc(ispec,i)
        if (i.lt.n(ispec)) then
          do k = i+1,n(ispec)
             indspc(ispec,k-1) = indspc(ispec,k)
             ispcbk(indspc(ispec,k-1))=k-1
          enddo
          indspc(ispec,n(ispec)) = ii
          ispcbk(indspc(ispec,n(ispec)))=n(ispec)
        endif
        n(ispec) = n(ispec)-1
      enddo  
        
      ntot=ntot-iv-1
      
      return
      end
      
c###################################################################
c###################################################################
c###################################################################
      
      subroutine add(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      dimension     inew(4)
   
      iv=igcval(igc)
   
      do icount=1,iv+1
          ispec=ispcin(igc,icount)  
          inew(icount) = n(ispec) + 1
          ii = indspc(ispec,inew(icount))
          if (ii.eq.0) then
            ii = (ispec-1)*nionmx+inew(icount)
            indspc(ispec,inew(icount)) = ii 
            ispcbk(ii)=inew(icount)
          endif
          n(ispec) = inew(icount)
          rx(ii)=rxnew(icount)
          ry(ii)=rynew(icount)
          rz(ii)=rznew(icount)
      enddo
      ntot=ntot+iv+1
      
      do icount=1,iv+1
        ispec=ispcin(igc,icount)  
        i=inew(icount)
        ii = indspc(ispec,i)
        do jspec=1,nspec
          do j=1,n(jspec)
              jj=indspc(jspec,j)
              if (ii.ne.jj) then
                rii(ii,jj)=riinew(icount,jj)
                rii(jj,ii)=riinew(icount,jj)
              endif
          enddo
        enddo
      enddo

      do icount=1,iv+1
        ispec=ispcin(igc,icount)  
        ii = indspc(ispec,inew(icount))
        do jcount=1,iv+1
          if (icount.ne.jcount) then
              jspec=ispcin(igc,jcount)  
              jj = indspc(jspec,inew(jcount))
              rii(ii,jj)=riicrn(icount,jcount)
          endif  
        enddo
      enddo
      
      return
      end

c#####################################################################
c#####################################################################
c#####################################################################

      subroutine encrea(igc,uii,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
      logical ovrlap
      ovrlap=.false.
      
      uii=0.0
   
      iv=igcval(igc)
   
      do icount1=1,iv
        ispec1=ispcin(igc,icount1)
        q1=q(ispec1)
        do icount2=icount1+1,iv+1
          ispec2=ispcin(igc,icount2)
          q2=q(ispec2)
          ddi=d2(ispec1,ispec2)
          rxij = rxnew(icount1)-rxnew(icount2)
          ryij = rynew(icount1)-rynew(icount2)
          rzij = rznew(icount1)-rznew(icount2)
          rxij = rxij-dnint(rxij/xl)*xl
          ryij = ryij-dnint(ryij/xl)*xl
          rijsq= rxij*rxij+ryij*ryij+rzij*rzij
          if (rijsq.lt.ddi) then
              ovrlap = .true.
              return
          endif 
          rij=dsqrt(rijsq)
          riicrn(icount1,icount2)=rij
          riicrn(icount2,icount1)=rij
          uii=uii + q1*q2/rij 
          zi = dabs(rzij)
          uii = uii + q1*(q2/xl**2)
     :              *(-twopi*zi-wpot(zi,xl))
        enddo     
c        uii = uii + q1*sigma(1)*(-twopi*(rznew(icount1)+zl)) 
c     :            + q1*sigma(2)*(-twopi*(zl-rznew(icount1)))
      enddo     
      
      do icount=1,iv+1
        ispec=ispcin(igc,icount)
        qi=q(ispec)
        do jspec=1,nspec
          qj=q(jspec)
          ddi=d2(ispec,jspec)
          do j=1,n(jspec)
             jj=indspc(jspec,j)
             rxij = rxnew(icount)-rx(jj)
             ryij = rynew(icount)-ry(jj)
             rzij = rznew(icount)-rz(jj)
             rxij = rxij-dnint(rxij/xl)*xl
             ryij = ryij-dnint(ryij/xl)*xl
             rijsq= rxij*rxij+ryij*ryij+rzij*rzij
             if (rijsq.lt.ddi) then
                ovrlap = .true.
                return
             endif 
             rij=dsqrt(rijsq)
             riinew(icount,jj)=rij
             uii=uii + qi*qj/rij
             zi = dabs(rzij)
             uii = uii + qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl))
          enddo 
        enddo
        uii = uii + qi*sigma(1)*(-twopi*(rznew(icount)+zl)) 
     :            + qi*sigma(2)*(-twopi*(zl-rznew(icount)))
      enddo 

      return
      end

c###################################################################
c#####################################################################
c#####################################################################
      
      subroutine endest(igc,uii,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      logical       runit
      character*2   fion
      character*4   fsalt
      character*9   fgz

      uii=0.0
      
      iv=igcval(igc)
      
      do icount=1,iv+1
        ispec=ispcin(igc,icount)
        i=iout(icount)
        ii=indspc(ispec,i)
        qi=q(ispec)
        do jspec=1,nspec
          qj=q(jspec)
          ddi=d2(ispec,jspec)
          do j=1,n(jspec)
             jj=indspc(jspec,j)
             if (ii.ne.jj) then
                 rij=rii(ii,jj)
                 uii=uii - qi*qj/rij 
                 zi = dabs(rz(ii)-rz(jj))
                 uii = uii - qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl))
             endif
          enddo 
        enddo 
        uii = uii - qi*sigma(1)*(-twopi*(rz(ii)+zl)) 
     :            - qi*sigma(2)*(-twopi*(zl-rz(ii)))
      enddo 
      
      do icount1=1,iv
        ispec1=ispcin(igc,icount1)
        i1=iout(icount1)
        ii1=indspc(ispec1,i1)
        q1=q(ispec1)
        do icount2=icount1+1,iv+1
          ispec2=ispcin(igc,icount2)
          i2=iout(icount2)
          ii2=indspc(ispec2,i2)
          q2=q(ispec2)
          rij=rii(ii1,ii2)
          uii=uii + q1*q2/rij 
          zi = dabs(rz(ii1)-rz(ii2))
          uii = uii + q1*(q2/xl**2)
     :              *(-twopi*zi-wpot(zi,xl))
        enddo     
c        uii = uii - q1*sigma(1)*(-twopi*(rz(ii1)+zl)) 
c     :            - q1*sigma(2)*(-twopi*(zl-rz(ii1)))
      enddo     

      return
      end

c    ###################################################################
c    ###################################################################

      function wpot(z,wl)
      implicit doubleprecision (a-h,o-z)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      x=z/wl
      x2=x*x
      r1=dsqrt(0.5+x2)
      r2=dsqrt(0.25+x2)
      r3=(0.5+r1)/r2
      wpot = 4*dlog(r3) - 4*x*( pi/2.0 - datan(4*x*r1) )
      wpot = wpot*wl
      return
      end

c    ###################################################################
c    ###################################################################

      subroutine iterate(ispec,i,rxi,ryi,rzi)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ riinew(nspcmx,ntotmx)
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      dimension epxtot(ntotmx),epytot(ntotmx),epztot(ntotmx)
      logical runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
c      write(*,*)"start iterate"
c-------------------        
         qi = q(ispec)
         ii = indspc(ispec,i)

c --- calculate E acting between ion ii and other ions 

c      open(56,file="iterate.dat")

      do jspec=1,nspec
       do j=1,n(jspec)
         jj=indspc(jspec,j)
         if (jj.eq.ii) then
           rxj=rxi
           ryj=ryi
           rzj=rzi
         else  
           rxj=rx(jj)
           ryj=ry(jj)
           rzj=rz(jj)
         endif
         eqxtot=0.0
         eqytot=0.0
         eqztot=0.0
         do kspec = 1,nspec
             qk = q(kspec)
             do k = 1,n(kspec)
               kk = indspc(kspec,k)
               if (kk.eq.ii) then
                 rxk=rxi
                 ryk=ryi
                 rzk=rzi
               else  
                 rxk=rx(kk)
                 ryk=ry(kk)
                 rzk=rz(kk)
               endif
               if (jj.ne.kk) then
                 if (jj.eq.ii) then
                    rij=riinew(ispec,kk)
                 elseif (kk.eq.ii) then
                    rij=riinew(ispec,jj)
                 else
                    rij=rii(jj,kk)
                 endif   
                 rxij = rxj-rxk
                 ryij = ryj-ryk
                 rzij = rzj-rzk
                 rxij = rxij-dnint(rxij/xl)*xl
                 ryij = ryij-dnint(ryij/xl)*xl
c                 rzij = rzij-dnint(rzij/xl)*xl
                 eqxn(jj,kk) = rxij*qk/(epsw*rij**3)
                 eqyn(jj,kk) = ryij*qk/(epsw*rij**3)
                 eqzn(jj,kk) = rzij*qk/(epsw*rij**3)
                 eqxtot = eqxtot+eqxn(jj,kk)
                 eqytot = eqytot+eqyn(jj,kk)
                 eqztot = eqztot+eqzn(jj,kk)
              endif  
             enddo    
         enddo 
         dipx0n(jj) = alpha(jspec)*eqxtot
         dipy0n(jj) = alpha(jspec)*eqytot
         dipz0n(jj) = alpha(jspec)*eqztot
       enddo
      enddo 
c if you have to print values,can print species by species individually!

      do iter0=1,4 
c--------------------------------      
        do jspec = 1,nspec  
           do j = 1,n(jspec)
             jj = indspc(jspec,j)
             if (jj.eq.ii) then
               rxj=rxi
               ryj=ryi
               rzj=rzi
             else  
               rxj=rx(jj)
               ryj=ry(jj)
               rzj=rz(jj)
             endif
             epxtot(jj)=0.0
             epytot(jj)=0.0
             epztot(jj)=0.0
             do kspec = 1,nspec
               do k = 1,n(kspec)
                 kk = indspc(kspec,k)
                 if (kk.eq.ii) then
                   rxk=rxi
                   ryk=ryi
                   rzk=rzi
                 else  
                   rxk=rx(kk)
                   ryk=ry(kk)
                   rzk=rz(kk)
                 endif
                 if (jj.ne.kk) then
                   if (jj.eq.ii) then
                      rij=riinew(ispec,kk)
                   elseif (kk.eq.ii) then
                      rij=riinew(ispec,jj)
                   else
                      rij=rii(jj,kk)
                   endif   
                   rxij = rxj-rxk     
                   ryij = ryj-ryk   
                   rzij = rzj-rzk    
                   rxij = rxij-dnint(rxij/xl)*xl
                   ryij = ryij-dnint(ryij/xl)*xl
c                   rzij = rzij-dnint(rzij/xl)*xl
c                   rij = dsqrt(rxij**2+ryij**2+rzij**2)
                   scal = dipxn(kk)*rxij
     :                   +dipyn(kk)*ryij
     :                   +dipzn(kk)*rzij
                   epxn(jj,kk) = 3*rxij*scal/rij**5-
     :                           dipxn(kk)/rij**3    
                   epyn(jj,kk) = 3*ryij*scal/rij**5-
     :                           dipyn(kk)/rij**3    
                   epzn(jj,kk) = 3*rzij*scal/rij**5-
     :                           dipzn(kk)/rij**3    
                   epxtot(jj)=epxtot(jj)+epxn(jj,kk)
                   epytot(jj)=epytot(jj)+epyn(jj,kk)
                   epztot(jj)=epztot(jj)+epzn(jj,kk)
                 endif
               enddo    
             enddo 
           enddo 
        enddo     
c------------------------------------------
        err0=0.0
        do jspec = 1,nspec
           do j = 1,n(jspec)
             jj = indspc(jspec,j)
             dipxo=dipxn(jj)
             dipyo=dipyn(jj)
             dipzo=dipzn(jj)
             dipxn(jj) = dipx0n(jj) + alpha(jspec)*epxtot(jj)
             dipyn(jj) = dipy0n(jj) + alpha(jspec)*epytot(jj)
             dipzn(jj) = dipz0n(jj) + alpha(jspec)*epztot(jj)
c             write(*,*)iter,jj,dipxo,dipxn(jj),epxtot(jj)
c             write(*,*)iter,jj,dipyo,dipyn(jj),alpha(jspec)
c             write(*,*)iter,jj,dipzo,dipzn(jj)
             err0=err0+(dipxn(jj)-dipxo)**2
     :                +(dipyn(jj)-dipyo)**2
     :                +(dipzn(jj)-dipzo)**2
           enddo 
        enddo     
        err0=dsqrt(err0)

      enddo

c      do jspec=1,nspec
c       do j=1,n(jspec) 
c         jj=indspc(jspec,j) 
c           write(56,*)epxtot(jj)
c       enddo
c      enddo
c      close(56)

C      stop
      
c      write(*,*)" subr. iterate err = ",err0

      return
      end

c    ###################################################################
c    ###################################################################

      subroutine fullit
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ riinew(nspcmx,ntotmx)
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/  dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :               dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :               dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :               dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :               eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :               eqzn(ntotmx,ntotmx),
     :               eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :               eqz(ntotmx,ntotmx),
     :               epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :               epzn(ntotmx,ntotmx),
     :               epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :               epz(ntotmx,ntotmx),
     :               alpha(nspcmx),ialpha
      dimension epxtot(ntotmx),epytot(ntotmx),epztot(ntotmx)
      logical runit
      character*2   fion
      character*4   fsalt
      character*9   fgz
c-------------------        

c --- calculate E acting between ion ii and other ions 

c     open(56,file="fullit.dat")
       do ispec=1,nspec
         epsi = eps(ispec)
         qi = q(ispec)
         do i=1,n(ispec)
          ii = indspc(ispec,i)
          eqxtot=0.0
          eqytot=0.0
          eqztot=0.0
          do jspec = 1,nspec
              qj = q(jspec)
              do j = 1,n(jspec)
                jj = indspc(jspec,j)
                if (ii.ne.jj) then
                  rij=rii(ii,jj)
                  rxij = rx(ii)-rx(jj)
                  ryij = ry(ii)-ry(jj)
                  rzij = rz(ii)-rz(jj)
                  rxij = rxij-dnint(rxij/xl)*xl  ! min. image convention
                  ryij = ryij-dnint(ryij/xl)*xl
c                  rzij = rzij-dnint(rzij/xl)*xl
                  eqx(ii,jj) = rxij*qj/(epsw*rij**3)!e.f.point-charges
                  eqy(ii,jj) = ryij*qj/(epsw*rij**3)
                  eqz(ii,jj) = rzij*qj/(epsw*rij**3)
                  eqxtot=eqxtot+eqx(ii,jj)
                  eqytot=eqytot+eqy(ii,jj)
                  eqztot=eqztot+eqz(ii,jj)
c                  write(56,*)ii,jj,eqx(ii,jj)
               endif  
              enddo    
          enddo 
c----- iteration 0 ------------------------------------        
          dipx0(ii) = alpha(ispec)*eqxtot
          dipy0(ii) = alpha(ispec)*eqytot
          dipz0(ii) = alpha(ispec)*eqztot
c          write(56,*)ii,dipx0(ii)
c          write(*,*)eqxtot
c          write(*,*)dipx0n(ii),dipy0n(ii),dipz0n(ii)
        enddo
       enddo
c       stop
      
c      write(*,*)"A"

        do jspec = 1,nspec
           do j = 1,n(jspec)
             jj = indspc(jspec,j)
             dipx(jj) = dipx0(jj)
             dipy(jj) = dipy0(jj)
             dipz(jj) = dipz0(jj)
c             write(*,*)dipxn(jj),dipyn(jj),dipzn(jj)
           enddo
        enddo

      do iter0=1,10  
c--------------------------------      
        do jspec = 1,nspec  
           do j = 1,n(jspec)
             jj = indspc(jspec,j)
             rxj=rx(jj)
             ryj=ry(jj)
             rzj=rz(jj)
             epxtot(jj)=0.0
             epytot(jj)=0.0
             epztot(jj)=0.0
             do kspec = 1,nspec
               do k = 1,n(kspec)
                 kk = indspc(kspec,k)
                 rxk=rx(kk)
                 ryk=ry(kk)
                 rzk=rz(kk)
                 if (jj.ne.kk) then
                   rij=rii(jj,kk)
                   rxij = rxj-rxk     
                   ryij = ryj-ryk   
                   rzij = rzj-rzk    
                   rxij = rxij-dnint(rxij/xl)*xl
                   ryij = ryij-dnint(ryij/xl)*xl
c                   rzij = rzij-dnint(rzij/xl)*xl
                   scal = dipx(kk)*rxij
     :                   +dipy(kk)*ryij
     :                   +dipz(kk)*rzij
                   epx(jj,kk) = 3*rxij*scal/rij**5-   ! e.f.-dipole
     :                           dipx(kk)/rij**3
                   epy(jj,kk) = 3*ryij*scal/rij**5-
     :                           dipy(kk)/rij**3    
                   epz(jj,kk) = 3*rzij*scal/rij**5-
     :                           dipz(kk)/rij**3    
                   epxtot(jj)=epxtot(jj)+epx(jj,kk) !total-point+dipole
                   epytot(jj)=epytot(jj)+epy(jj,kk)
                   epztot(jj)=epztot(jj)+epz(jj,kk)
                 endif
               enddo    
             enddo
c             if (jj.eq.1) write(*,*)epxtot(jj),epytot(jj),epztot(jj)
           enddo 
        enddo     
c------------------------------------------
        err0=0.0
        do jspec = 1,nspec
           do j = 1,n(jspec)
             jj = indspc(jspec,j)
             dipxo=dipx(jj)  ! for first iteration same as for point
          !   charges
             dipyo=dipy(jj)
             dipzo=dipz(jj)
             dipx(jj) = dipx0(jj) + alpha(jspec)*epxtot(jj)
             ! updating dip.moments
             dipy(jj) = dipy0(jj) + alpha(jspec)*epytot(jj)
             dipz(jj) = dipz0(jj) + alpha(jspec)*epztot(jj)
c             write(*,*)iter,jj,dipxo,dipxn(jj),epxtot(jj)
c             write(*,*)iter,jj,dipyo,dipyn(jj),alpha(jspec)
c             write(*,*)iter,jj,dipzo,dipzn(jj)
             err0=err0+(dipx(jj)-dipxo)**2
     :                +(dipy(jj)-dipyo)**2
     :                +(dipz(jj)-dipzo)**2
           enddo 
        enddo     
        err0=dsqrt(err0)

      enddo

c      open(56,file="fieldp-fu.dat")
c      do jspec=1,nspec
c       do j=1,n(jspec) 
c         jj=indspc(jspec,j) 
c           write(56,*)epxtot(jj)
c       enddo
c      enddo
c      close(56)

C      stop
      
      write(*,*)" subr. iterate err = ",err0

      return
      end
