      program ionch   
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/names/ fion(10),fsalt(10),fgz(10),fgin(10),fc,fout,fh,
     :              fnaca(10,4),fdata,focc,frdf(10,10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      common/simp/  nstep,naver,isave
      common/bulk/  zbulk1,zbulk2,rbulk,vbulk,anaca(nspcmx,4,20000),
     :              asuba(nspcmx,4,20000),ksub
      common/occ/   aocc(3,0:4,0:4,0:4,0:4),nfree
      common/radf/  ardf(3,nspcmx,nspcmx),
     :              rdf(3,nspcmx,nspcmx,nrgmx),
     :              drdf,nrdf,calrdf
      character*11  fgz,fgin,fc,fout,fnaca,fh,fdata,focc
      character*15  frdf
      character*2   fion
      character*4   fsalt
      logical       chonly
      logical       calrdf,calgin

      write(*,*)" Start MC for calcium channel"

      open(2,file="in.dat")
      read(2,*)     
      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 accum
             enddo 
             call hist
             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=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/names/ fion(10),fsalt(10),fgz(10),fgin(10),fc,fout,fh,
     :              fnaca(10,4),fdata,focc,frdf(10,10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      common/simp/  nstep,naver,isave
      common/bulk/  zbulk1,zbulk2,rbulk,vbulk,anaca(nspcmx,4,20000),
     :              asuba(nspcmx,4,20000),ksub
      common/occ/   aocc(3,0:4,0:4,0:4,0:4),nfree
      common/radf/  ardf(3,nspcmx,nspcmx),
     :              rdf(3,nspcmx,nspcmx,nrgmx),
     :              drdf,nrdf,calrdf
      character*11  fgz,fgin,fc,fout,fnaca,fh,fdata,focc
      character*15  frdf
      character*2   fion
      character*4   fsalt
      logical       chonly
      logical       readit, ovrlap, targc
      logical       calrdf,calgin

c_______ CONSTANTS __________________________________________________
c      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 __________________________________
 
      call readi(mag,fgz,fgin,fc,fout,fnaca,fh,fdata,focc,frdf)
     
      open(2,file="in.dat")
     
c---read in first AND geometry---
      do i=1,16
        read(2,*) 
      enddo

c---simulation parameters--
      read(2,*)drmax(1),drmax(2)
      read(2,*)drg,dzgtrg  
      read(2,*)readit
      read(2,*)nstep0
      read(2,*)naver 
      read(2,*)isave
      read(2,*)calgin
      read(2,*)calrdf
      
      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"
      ntot=0
      ntry=0
      nimax=0
      istr=0
      do ispec=1,nspec 
         read(2,*)eps(ispec),n(ispec),itry(ispec),ratmov(ispec),
     :             ratexc(ispec),z(ispec),d(ispec),chonly(ispec),
     :             fion(ispec)
         ntot=ntot+n(ispec)
         ntry=ntry+itry(ispec)
         if (n(ispec).gt.nimax) nimax=n(ispec)
         if (chonly(ispec)) istr=istr+1
      enddo      
      nfree=nspec-istr
      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"

      read(2,*)
      read(2,*)
c---grand canonical---
      read(2,*)nsalt
      read(2,*)
      ratgrt=0.0
      do i=1,nsalt
        read(2,*)isalt(i),ctarg(i),chex(i),ratgr(i),fsalt(i)
        ratgrt=ratgrt+ratgr(i)
      enddo
      read(2,*)
      read(2,*)ratreg(1),ratreg(2),ratreg(3),ratreg(4)

      close(2)

      open(2,file="geometry.dat")
       read(2,*)zl(1),zl(2),zl(3)                        
       read(2,*)zl(4)                                
       read(2,*)rl(1),rl(2),rl(3)                        
       read(2,*)rl(4),rl(5)                               
       read(2,*)rlvest,rlcurv                              
       read(2,*)epsw,epspr,deps                               
      close(2)

c_______ FIX ZL AND NSTEP ACCORDING TO ITERATION STATUS _____________
c_______ VERY SPECIFIC TO THE ACTUAL GEOMETRY _______________________ 
      
      tsi=300.0
      avog=6.02214d23
      tosi=1.0/(1.d-10)**3/avog/1000.
      write(*,*)" tosi = ",tosi
      dsi=1.d-10
      eps0=8.8542d-12
      boltz=1.3806d-23
      echg=1.6021917d-19
      beta=1./(boltz*tsi)
     
      zocc=2.5

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

      write(*,'(72("_"))')
      write(*,'(72("_"))')
      write(*,*)" Regions for particle insertion/deletion"
      do ispec=istr+1,nspec
       ri=d(ispec)/2.0
       if (ispec.eq.istr+1) then
          zcr1(1,ispec)=-zl(4)
          zcr2(1,ispec)= zl(4)
          rcr(1,ispec) = rl(5)    
          zcr1(2,ispec)=-zl(4)
          zcr2(2,ispec)= zl(4)
          rcr(2,ispec) = rl(5)    
          zcr1(3,ispec)=-zl(4)
          zcr2(3,ispec)= zl(4)
          rcr(3,ispec) = rl(5)    
          zcr1(4,ispec)=-zl(4)
          zcr2(4,ispec)= zl(4)
          rcr(4,ispec) = rl(5)    
       else
          zcr1(1,ispec)=-zocc
          zcr2(1,ispec)= zocc
          rcr(1,ispec) = rl(1)-ri 
          zcr1(2,ispec)=-zl(1)
          zcr2(2,ispec)= zl(1)
          rcr(2,ispec) = rl(1)-ri 
          zcr1(3,ispec)=-zl(2)   
          zcr2(3,ispec)= zl(2)   
          rcr(3,ispec) = rl(2)      
          zcr1(4,ispec)=-zl(4)
          zcr2(4,ispec)= zl(4)
          rcr(4,ispec) = rl(5)    
       endif
       do ireg=1,4 
        write(*,*)fion(ispec)," in region ",ireg,": z1-z2 = ", 
     :      zcr1(ireg,ispec),zcr2(ireg,ispec)
        write(*,*)fion(ispec)," in region ",ireg,": r = ", 
     :      rcr(ireg,ispec)
       enddo 
       write(*,'(72("_"))')
      enddo 

      write(*,'(72("_"))')
      write(*,'(72("_"))')
      write(*,*)" Regions for calculating average ion numbers"
      do ispec=istr+1,nspec
          ri=d(ispec)/2.0
          zreg1(1,ispec)=-zocc       
          zreg2(1,ispec)= zocc       
          rreg(1,ispec) = rl(1)-ri    
          zreg1(2,ispec)=-zl(1)      
          zreg2(2,ispec)= zl(1)      
          rreg(2,ispec) = rl(1)-ri   
          zreg1(3,ispec)=-zl(2)   
          zreg2(3,ispec)= zl(2)   
          rreg(3,ispec) = rl(2)        
          zreg1(4,ispec)=-zl(4)    
          zreg2(4,ispec)= zl(4)      
          rreg(4,ispec) = rl(5)      
          do ireg=1,4
            write(*,*)fion(ispec)," in region ",ireg,": z1-z2 = ", 
     :         zreg1(ireg,ispec),zreg2(ireg,ispec)
            write(*,*)fion(ispec)," in region ",ireg,": r = ", 
     :         rreg(ireg,ispec)
          enddo
          write(*,'(72("_"))')
      enddo 
      write(*,'(72("_"))')
      write(*,'(72("_"))')

      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)=istr+1
          enddo
          do ireg=1,4  
            vgc(ireg,igc) = 1.0
            do i=1,iv+1
               ispec = ispcin(igc,i)
               ri = d(ispec)/2.0
               rmaxsq = rcr(ireg,ispec)**2
               vgc(ireg,igc)=vgc(ireg,igc)*rmaxsq*pi
     :                       *(zcr2(ireg,ispec)-zcr1(ireg,ispec))
            enddo
            write(*,*)" vgc = ",ireg,igc,vgc(ireg,igc)
          enddo  
      enddo 

      ctrgcl=0.0
      do igc=1,nsalt
        iv=igcval(igc)
        ctrgcl=ctrgcl+iv*ctarg(igc)
      enddo  
        
      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(istr+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
        
      zbulk1 = zl(2) + 0.3*(zl(4)-zl(2))
      zbulk2 = zl(4) - 0.3*(zl(4)-zl(2))
      rbulk  = 2.0*rl(5)/3.0
      vbulk  = rbulk**2*pi*(zbulk2-zbulk1)*2
              
c_______ CALCULATE REDUCED CHARGES __________________________________
c_______ CHECK FOR ELECTRONEUTRALITY ________________________________
c_______ FILL UP DD, D2, AND NWH ARRAYS _____________________________

      charge=0.0
      ii=0
      qstar=echg*dsqrt(beta/(4*pi*eps0*dsi))
      tstar=4*pi*eps0*dsi/(beta*echg**2)
      write(*,'(" q*          = ",f10.5)')qstar 
      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"

c_______ WRITE PARAMETERS ___________________________________________
      write(*,*)
      write(*,*) " SIMULATION PARAMETERS IN SI UNITS"
      write(*,*)
      write(*,'(" T [K]       = ",f10.5)')tsi  
      write(*,*)             
      write(*,'(72("_"))')
      write(*,*)"           filter radius rl(1) = ", rl(1)
      write(*,*)"  vestibule outer radius rl(2) = ", rl(2)
      write(*,*)"  flat wall outer radius rl(3) = ", rl(3)
      write(*,*)"   outer cylinder radius rl(4) = ", rl(4)
      write(*,*)"    cell cylinder radius rl(5) = ", rl(5)
      write(*,'(72("_"))')
      write(*,*)"           filter length zl(1) = ", zl(1)
      write(*,*)"          channel length zl(2) = ", zl(2)
      write(*,*)"   outer cylinder length zl(3) = ", zl(3)
      write(*,*)"             cell length zl(4) = ", zl(4)
      write(*,'(72("_"))')
      write(*,*)"              vestibule radius = ", rlvest
      write(*,*)"              curveture radius = ", rlcurv
      write(*,'(72("_"))')
      write(*,*)"    eps solvent = ", epsw
      write(*,*)"    eps protein = ", epspr
      write(*,'(72("_"))')
      write(*,*)      
      write(*,*)" Number of structural ion species = ", istr
      
      write(*,'(12x,"   N   ","  q [e]  ","  d [A]  ")')
      do i=1,nspec
      write(*,'("  ion   ",i1,1x,i7,2x,f5.2,2x,2x,f8.4)')
     :      i,n(i),z(i),d(i)
      enddo
     
      write(*,*)
      
      if (readit) then
 
c_______ READIN CONFIGURATION _______________________________________
c_______ FILL UP INDEX ARRAYS AND LOOKUP TABLES _____________________
              
            open (3,file="c0")
            read(3,*)nspec
            read(3,*)
            read(3,*)
            read(3,*)
            read(3,*)
            read(3,*)
           
            do ispec=1,nspec
                  read(3,*)n(ispec),d(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
                        r2(ii)=dsqrt(rx(ii)**2+ry(ii)**2)
                        indspc(ispec,i)=ii
                        ispcbk(ii)=i
                  enddo      
            read(3,*)
            enddo
            close(3)

      else

c_______ MAKE CONFIGURATION _________________________________________
            call initc

      endif

      call lookup

      do ireg=1,3
        do ispec=1,nspec
          nin(ireg,ispec)=0
          do i=1,n(ispec)
            ii=indspc(ispec,i)
            zi =rz(ii)
            if (zi.gt.zreg1(ireg,ispec).and.
     :          zi.lt.zreg2(ireg,ispec)) then
                nin(ireg,ispec)=nin(ireg,ispec)+1
                indreg(ireg,ispec,nin(ireg,ispec))=ii 
            endif  
          enddo   
        enddo   
      enddo   

c_______ VERY SPECIFIC TO THE ACTUAL GEOMETRY _______________________ 
     
c_______ Geometry for jmpin and jmpout ______________________________ 
      do ispec=1,nspec
         rljpin(ispec) = rl(2)
         zljpin(ispec) = zl(2)+d(ispec)/2.0
         vin(ispec) =rljpin(ispec)**2*zljpin(ispec)
         vout(ispec)=rl(5)**2*(zl(4)-zljpin(ispec))
      enddo

c_______ Geometry for gin 2D distributions __________________________ 
      rinup  = rl(4)+10.0
      zinlft = zl(2)+10.0
      nrgr   = int(rinup)/drg+1
      rinup  = nrgr*drg
      nrgz   = int(zinlft)/drg+1
      zinlft = -nrgz*drg
      nrgz   = 2*nrgz
       write(*,'(i1," nrg in r dimension= ",i8," drg = ",f8.4)')
     :      i,nrgr,drg
       write(*,'(i1," nrg in z dimension= ",i8," drg = ",f8.4)')
     :      i,nrgz,drg
      if (nrgr.gt.nrgz) then
       nrg=nrgr
      else
       nrg=nrgz
      endif 
       if (nrg.gt.nrgmx) then 
         write(*,*) "nrg is too large"
         write(*,'(" nrg = ",i6,"   nrgmx = ",i6)')nrg,nrgmx
         stop
       endif  
     
c_______ Geometry for gz distributions ______________________________ 
      do i=1,nspec

             nzreg(i)=5

             ri=d(i)/2.0
             zzreg(0,i)=-zl(4)
             zzreg(1,i)=-zl(2)-ri
             if (chonly(ispec)) then
              zzreg(2,i)=-zl(1)+ri    
              zzreg(3,i)=zl(1)-ri
             else
              zzreg(2,i)=-zl(1)    
              zzreg(3,i)= zl(1)  
             endif
             zzreg(4,i)=zl(2)+ri
             zzreg(5,i)=zl(4)
       
             nzgtot(i)=0
             do j=1,5
               nzg(j,i)  = int((zzreg(j,i)-zzreg(j-1,i))/dzgtrg)+1
               nzgtot(i)=nzgtot(i)+nzg(j,i)
               dzg(j,i)  = (zzreg(j,i)-zzreg(j-1,i))/dble(nzg(j,i))
               write(*,'(a2,", region ",i2,": From ",f10.3," to ",f10.3,
     :            ": nzg = ",i8," dzg = ",f8.4)')
     :            fion(i),j,zzreg(j-1,i),zzreg(j,i),nzg(j,i),dzg(j,i)
             enddo
             write(*,*)" nzgtot = ",nzgtot(i)
             if (nzgtot(i).gt.nzgmx) then 
               write(*,*) "nzg is too large"
               write(*,'(" nzg = ",i6,"   nzgmx = ",i6)')nzgtot(i),nzgmx
               stop
             endif  
         
         write(*,'(40("_"))')

      enddo

c_______ Geometry for rdf ___________________________________________ 
      drdf=drg
      nrdf=int(13.0/drdf)+1

      call calv

      return
      end

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

      subroutine tryit
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical       calrdf,calgin
         
      vszam=ranff(mag)
      ratrg2=0.0
      do irg=1,4    
        ratrg2=ratrg2+ratreg(irg)
        ratrg1=ratrg2-ratreg(irg)
        if (vszam.gt.ratrg1.and.vszam.lt.ratrg2) then   
           ireg = irg
        endif   
      enddo
c        write(*,*)vszam,ratreg(1),ratreg(2),ratreg(3),ratreg(4),ireg
      
      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(ireg,igc)
           else
             call destr(ireg,igc)
           endif
        endif
      enddo
      
      if (vszam.gt.ratgr2) then
          iwh   = int(ranff(mag)*ntry)+1
          ispec = nwh(iwh)
          if (n(ispec).eq.0) goto 1
          i = int(dble(n(ispec))*ranff(mag)) + 1
          ii=indspc(ispec,i)
          if (chonly(ispec)) then
             call move(ispec,i)
          else
              if (ranff(mag).gt.ratexc(ispec)) then
                 call move(ispec,i)
              else
                 rzold = dabs(rz(ii))
                 if (rzold.gt.zljpin(ispec)) then
                    call jmpin(ispec,i)
                 else                          
                    call jmpout(ispec,i)
                 endif
              endif
          endif
1         continue
      endif

      return
      end

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

      subroutine wall(ispec,drznew,r2inew,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      logical       chonly
      logical ovrlap
      ovrlap=.false.

      di=d(ispec)
      ri=di/2.0
      rl10=rl(1)-ri
      zl10=zl(1)-ri
      zl20=zl(2)+ri
      rl50=rl(5)
      zl40=zl(4)
      rl4=(rlvest+ri)*(rlvest+ri)
      
      if (chonly(ispec)) then
              if (r2inew.gt.rl10.or.drznew.gt.zl10) then
                      ovrlap=.true.
                      return
              endif        
      else
              if (drznew.lt.zl(1)) then
                      if (r2inew.gt.rl10) then
                        ovrlap=.true.
                        return
                      endif  
              elseif (drznew.gt.zl20) then
                      if (r2inew.gt.rl50.or.drznew.gt.zl40) then  
                        ovrlap=.true.
                        return
                      endif  
              else
                      rsq=(rl(2)-r2inew)*(rl(2)-r2inew)+
     :                    (drznew-zl(1))*(drznew-zl(1))
                      if (rsq.lt.rl4.or.r2inew.gt.rl(2)) then  
                        ovrlap=.true.
                        return
                      endif  
              endif        
      endif
      return
      end

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

      subroutine lookup
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      logical       chonly

      if (dabs(deps).gt.1.d-5) then
      
         write(*,*)" I am reading in A (aold), LU decomposed A (aa),",
     :             " and indx "
         open(12,file="amx.dat")  
         open(22,file="indx.dat")
          read(12,*)npatch
          read(12,*)  
          do i=1,npatch
            read(22,*)ii,indx(i)
            do j=1,npatch
               read(12,*)ii,jj,amx(i,j)
            enddo
            read(12,*)  
          enddo
         close(12)
         close(22)
         write(*,*)" stuff read in"
        
c         open(2,file="npatch.dat")
c         read(2,*)npatch
c         close(2)
       
         open(2,file="centers.dat")
         open(12,file="normal.dat")
         open(22,file="area.dat")
         do i=1,npatch
           read(2,*)px(i),py(i),pz(i)
           read(12,*)ux(i),uy(i),uz(i)
           read(22,*)a(i)            
         enddo
         close(2)
         close(12)
         close(22)
      
      endif
     
      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
                         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
              
             if (dabs(deps).gt.1.d-5) then
               do jj=1,npatch        
                  rxj=px(jj)
                  ryj=py(jj)
                  rzj=pz(jj)
                  rxij=rxi-rxj
                  ryij=ryi-ryj
                  rzij=rzi-rzj
                  rijsq= rxij*rxij+ryij*ryij+rzij*rzij
                  rip(ii,jj)=dsqrt(rijsq)
               enddo     
             endif  

         enddo
      enddo

      if (dabs(deps).gt.1.d-5) then
        
        do k=1,npatch
          rxk=px(k)
          ryk=py(k)
          rzk=pz(k)
          c(k)=0.0
          do jspec=1,nspec
              qj   = q(jspec)
              epsj=eps(jspec)
              do j=1,n(jspec)
                 jj=indspc(jspec,j)
                 rxj=rx(jj)
                 ryj=ry(jj)
                 rzj=rz(jj)
                 rij=rip(jj,k)
                 rzkj = (rxk-rxj)*ux(k)+(ryk-ryj)*uy(k)+(rzk-rzj)*uz(k)
                 c(k)=c(k)-deps*qj*rzkj/(4*pi*epsj*rij*rij*rij)
              enddo
          enddo  
          h(k)=c(k)
        enddo

c_______ CALCULATE INITIAL H ________________________________________

        call lubksb(amx,npatch,indx,h)
       
        areat=0.0
        chgt=0.0
        do i=1,npatch
         chgt=chgt+a(i)*h(i)
         areat=areat+a(i)
        enddo
        chgt=chgt/qstar
        write(*,*)" Total initial induced charge / e = ",chgt
        write(*,*)" Total area = ",areat
      
      endif
      
      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 accum
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/occ/   aocc(3,0:4,0:4,0:4,0:4),nfree
      logical       chonly
      logical       calrdf,calgin
  
      do ireg=1,3
        i1=nin(ireg,istr+1)
        i2=nin(ireg,istr+2)
        i3=nin(ireg,istr+3)
        i4=nin(ireg,istr+4)
        aocc(ireg,i1,i2,i3,i4)=aocc(ireg,i1,i2,i3,i4)+1.0
      enddo

      do jspec=1,nspec
        anin(1,jspec)=anin(1,jspec)+nin(1,jspec)
        anin(2,jspec)=anin(2,jspec)+nin(2,jspec)
        anin(3,jspec)=anin(3,jspec)+nin(3,jspec)
        an(jspec)=an(jspec)+dble(n(jspec))
      enddo
    
      do i=1,npatch
        ah(i)=ah(i)+h(i)
      enddo
    
      ataccu=ataccu+1.0
      
      return
      end
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine hist 
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/bulk/  zbulk1,zbulk2,rbulk,vbulk,anaca(nspcmx,4,20000),
     :              asuba(nspcmx,4,20000),ksub
      common/radf/  ardf(3,nspcmx,nspcmx),
     :              rdf(3,nspcmx,nspcmx,nrgmx),
     :              drdf,nrdf,calrdf
      logical       chonly
      logical       calrdf,calgin
    
      athist=athist+1.0
      
      do ispec=1,nspec

         ni=n(ispec)
         ri=d(ispec)/2.0
         do i=1,ni
            
            ii=indspc(ispec,i)
            rzi=rz(ii)
            r2i=r2(ii)
            
            nsofar=0
            do j=1,nzreg(ispec)
              if (rzi.gt.zzreg(j-1,ispec).and.rzi.le.zzreg(j,ispec))
     :          then
                    k=int((rzi-zzreg(j-1,ispec))/dzg(j,ispec))+1
                    gz(ispec,k+nsofar)=gz(ispec,k+nsofar)+1.0
              endif
              nsofar=nsofar+nzg(j,ispec)
            enddo
            
            if (calgin) then 
               if (dabs(rzi).lt.-zinlft.and.r2i.lt.rinup) then
                    k=int((rzi-zinlft)/drg)+1
                    l=int(r2i/drg)+1
                    gin(nspec+1,k,l)=gin(nspec+1,k,l)+z(ispec)
                    gin(ispec,k,l)=gin(ispec,k,l)+1.0     
               endif        
            endif   
           
            if (dabs(rzi).gt.zbulk1.and.dabs(rzi).lt.zbulk2
     :            .and.r2i.lt.rbulk) 
     :        anin(4,ispec)=anin(4,ispec)+1.0

         enddo 
      
      enddo
   
      if (calrdf) then

         do ispec=1,nspec
           
           do i=1,nin(1,ispec)
             ii=indreg(1,ispec,i)
             rxi=rx(ii)
             ryi=ry(ii)
             rzi=rz(ii)
             do jspec=1,nspec
               do j=1,n(jspec)
                 jj=indspc(jspec,j)
                 rzj=rz(jj)
                 if (dabs(rzi-rzj).lt.20.0) then
                   rxj=rx(jj)
                   ryj=ry(jj)
                   rijsq=(rxi-rxj)**2+(ryi-ryj)**2+(rzi-rzj)**2
                   if (rijsq.lt.225.0) then
                     rij=dsqrt(rijsq)
                     k=int((rij-dd(ispec,jspec))/drdf)+1
                     if (j.le.nrdf) then
                       rdf(1,ispec,jspec,k)=rdf(1,ispec,jspec,k)+1.0
                       ardf(1,ispec,jspec)=ardf(1,ispec,jspec)+1.0
                     endif  
                   endif  
                 endif
               enddo
             enddo  
           enddo
         
           do i=1,nin(2,ispec)
             ii=indreg(2,ispec,i)
             rxi=rx(ii)
             ryi=ry(ii)
             rzi=rz(ii)
             do jspec=1,nspec
               do j=1,n(jspec)
                 jj=indspc(jspec,j)
                 rzj=rz(jj)
                 if (dabs(rzi-rzj).lt.20.0) then
                   rxj=rx(jj)
                   ryj=ry(jj)
                   rijsq=(rxi-rxj)**2+(ryi-ryj)**2+(rzi-rzj)**2
                   if (rijsq.lt.225.0) then
                     rij=dsqrt(rijsq)
                     k=int((rij-dd(ispec,jspec))/drdf)+1
                     if (j.le.nrdf) then
                       rdf(2,ispec,jspec,k)=rdf(2,ispec,jspec,k)+1.0
                       ardf(2,ispec,jspec)=ardf(2,ispec,jspec)+1.0
                     endif  
                   endif  
                 endif
               enddo
             enddo  
           enddo
         
         enddo  

      endif

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

      subroutine initc
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      logical       chonly
      logical ovrlap
      ovrlap=.false.
      
      write(*,*)" Making initial configuration"

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

      open(22,file="in-structural.dat")

      ii=0
      do 10 ispec=1,nspec

      if (chonly(ispec)) then
        
        do i=1,n(ispec)
          read(22,*)rxi,ryi,rzi
          ii=(ispec-1)*nionmx+i
          rx(ii)=rxi
          ry(ii)=ryi
          rz(ii)=rzi
          r2i=dsqrt(rxi**2+ryi**2)
          r2(ii)=r2i
          indspc(ispec,i)=ii
          ispcbk(ii)=i
        enddo

      else

        ntrg=n(ispec)
        ni=0
        ri=d(ispec)/2.0
        
30      ni=ni+1
        if (chonly(ispec)) then
              rlisq=(rl(1)-ri)**2
        else      
              rlisq=rl(5)**2
        endif     
20      rsqnew = ranff(mag) * rlisq
        r2inew = dsqrt(rsqnew)
        
        ii=ii+1
        if (ii.gt.1000000) goto 11      
        
        if (chonly(ispec)) then
              rzinew = (2*ranff(mag)-1)*(zl(1)-ri)
        else
              rzinew = (2*ranff(mag)-1)*zl(4)
              drznew = dabs(rzinew)
              call wall(ispec,drznew,r2inew,ovrlap)
              if (ovrlap) goto 20
        endif
        
        finew  = ranff(mag) * twopi
        rxinew = r2inew * dcos(finew)      
        ryinew = r2inew * dsin(finew)
        
        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
         r2(ii)=r2inew
         indspc(ispec,ni)=ii
         ispcbk(ii)=ni
c         write(*,*)ii,ispec,ni,rxinew,ryinew,rzinew
        endif 
        
        if (ni.lt.ntrg) goto 30

      endif

10    continue

      close(22)
      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=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      logical       chonly
      logical       calrdf,calgin
      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
                     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=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/names/ fion(10),fsalt(10),fgz(10),fgin(10),fc,fout,fh,
     :              fnaca(10,4),fdata,focc,frdf(10,10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      character*11  fgz,fgin,fc,fout,fnaca,fh,fdata,focc
      character*15  frdf
      character*2   fion
      character*4   fsalt
      logical       chonly
      logical       calrdf,calgin
c      write(*,*)' begin wrcf'
      open(2,file=fc)  
      write(2,*) nspec,istep                             
      write(2,*) zl(1),zl(2),zl(3)           
      write(2,*) rl(1),rl(2),rl(3),rl(4)           
      write(2,*) rl(5),zl(4)
      write(2,*) rlvest,rlcurv 
      write(2,*)
      do ispec=1,nspec
            write(2,*) n(ispec),d(ispec)
            write(2,*) nzgtot(ispec),nrgz,nrgr
            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=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/names/ fion(10),fsalt(10),fgz(10),fgin(10),fc,fout,fh,
     :              fnaca(10,4),fdata,focc,frdf(10,10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      common/simp/  nstep,naver,isave
      common/bulk/  zbulk1,zbulk2,rbulk,vbulk,anaca(nspcmx,4,20000),
     :              asuba(nspcmx,4,20000),ksub
      common/occ/   aocc(3,0:4,0:4,0:4,0:4),nfree
      common/radf/  ardf(3,nspcmx,nspcmx),
     :              rdf(3,nspcmx,nspcmx,nrgmx),
     :              drdf,nrdf,calrdf
      logical       chonly
      logical       calrdf,calgin
      character*11  fgz,fgin,fngz,fngin,fc,fout,fnaca,fh,focc,fdata
      character*15  frdf
      character*2   fion
      character*4   fsalt
      dimension     cz(nspcmx,nzgmx), cin(nspcmx,nrgmx,nrgmx),
     :              cav(nzgmx), rzi(nzgmx), r2i(nrgmx) 
      dimension     cfilt(nspcmx),cchan(nspcmx),ctot(nspcmx),
     :              czocc(nspcmx)   
     
      write(*,*)" Begin saves"	 
      open(2,file=fout)
        
      write(2,'(" zl4 [A]     = ",f10.5)')zl(4)
      write(2,'(" rl5 [A]     = ",f10.5)')rl(5)
      write(2,*)
      write(2,'(" rl1 [A]     = ",f10.5)')rl(1)
      write(2,'(" epspr       = ",f10.5)')epspr   
      write(2,*)
      do igr=1,nsalt
        write(2,'(a4,": ctarg = ",f11.8,", chemp = ",f10.5,
     :       ", excess = ",f10.6)')
     :       fsalt(igr),ctarg(igr),chemp(igr),chex(igr) 
      enddo 
      write(2,*)
      write(2,'(12x,"   N   ","  q [e]  ","  d [A]  ")')
      do i=1,nspec
      write(2,'(a2," ion: ",1x,i7,2x,f5.2,2x,2x,f8.4)')
     :      fion(i),n(i),z(i),d(i)
      enddo
      ksub=ksub+1
      
c ##### Calculate and write out concentration profiles ###############

      if (calgin) then
        do ispec=1,nspec+1
            fngin = fgin(ispec)
            open(12,file=fngin)
            do i=1,nrgz 
               do j=1,nrgr         
                  cin(ispec,i,j) = 
     :                 tosi*gin(ispec,i,j)/athist/vjin(ispec,i,j)
                  zri=zinlft+(i-1)*drg+drg/2.0
                  rrj=(j-1)*drg+drg/2.0
                  write(12,*)zri,rrj,cin(ispec,i,j),vjin(ispec,i,j)
               enddo
               write(12,*)
            enddo   
            close(12)
        enddo
      endif
        
        do ispec=1,nspec
          fngz  = fgz(ispec)
          open(22,file=fngz)
           do i=1,nzgtot(ispec)
             anz = tosi*gz(ispec,i)/athist
             cz(ispec,i) = tosi*gz(ispec,i)/athist/vjz(ispec,i)
             write(22,*)rgz(ispec,i),cz(ispec,i),anz,vjz(ispec,i),
     :                  dzgi(ispec,i)
           enddo   
          close(22)
        enddo

      if (calrdf) then
        do ispec=1,nspec
          do jspec=1,nspec
            ddi=dd(ispec,jspec)
            open(22,file=frdf(ispec,jspec))
            do j=1,nrdf
              rlow=ddi+(j-1)*drdf
              rhigh=rlow+drdf
              ri=rlow+drdf/2.0
              vshell=4*pi*(rhigh**3-rlow**3)/3.0
              if (ardf(1,ispec,jspec).gt.0.01) then
                rdfi1=rdf(1,ispec,jspec,j)/ardf(1,ispec,jspec)
              else
                rdfi1=0.0
              endif  
              if (ardf(2,ispec,jspec).gt.0.01) then
                rdfi2=rdf(2,ispec,jspec,j)/ardf(2,ispec,jspec)
              else
                rdfi2=0.0
              endif  
              write(22,*)ri,rdfi1,rdfi2,rdfi1/vshell,rdfi2/vshell
            enddo
            close(22)
          enddo
        enddo
      endif

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

      chfilt=0.0
      chzocc=0.0
      chchan=0.0
      chbulk=0.0
      do ispec=1,nspec
        czocc(ispec)=anin(1,ispec)/ataccu
        cfilt(ispec)=anin(2,ispec)/ataccu
        cchan(ispec)=anin(3,ispec)/ataccu
        cbulk(ispec)=tosi*anin(4,ispec)/athist/vbulk
        ctot(ispec)=an(ispec)/ataccu
        chfilt=chfilt+z(ispec)*cfilt(ispec)
        chzocc=chzocc+z(ispec)*czocc(ispec)
        chchan=chchan+z(ispec)*cchan(ispec)
        if (.not.chonly(ispec)) chbulk=chbulk+z(ispec)*cbulk(ispec)
      enddo
      chfilt=chfilt
      chzocc=chzocc
      chchan=chchan
      
      ina=istr+2
      ica=istr+3
      do ispec=1,nspec
         anaca(ispec,1,ksub)=czocc(ispec)
         anaca(ispec,2,ksub)=cfilt(ispec)
         anaca(ispec,3,ksub)=cchan(ispec)
         anaca(ispec,4,ksub)=cbulk(ispec)
         if (ksub.gt.1) then
           do j=1,4
            asuba(ispec,j,ksub)=ksub*anaca(ispec,j,ksub)
     :                         -(ksub-1)*anaca(ispec,j,ksub-1)
           enddo
         else
           do j=1,4
            asuba(ispec,j,ksub)=anaca(ispec,j,ksub)
           enddo
         endif  
      enddo
     
      do ispec=istr+2,nspec
        do j=1,4
          open(12,file=fnaca(ispec,j))
          do k=1,ksub
            write(12,*)k,anaca(ispec,j,k),asuba(ispec,j,k)
          enddo  
          close(12)
        enddo
      enddo  
      
        write(*,'(72("#"))')
        write(*,'(i9,". cycle -------- ",i2,". iteration")')istep,iter
        write(*,'(72("_"))')
        write(*,'(10x,"|",4x,"BULK",4x,"|",
     :        3x,"CHANNEL",2x,"|",3x,"FILTER",2x,"|",
     :        3x,"CENTER ",                        
     :        3x,"|",3x,"TOTAL  ")')
        
        write(2,'(72("#"))')
        write(2,'(i9,". cycle -------- ",i2,". iteration")')istep,iter
        write(2,'(72("_"))')
        write(2,'(10x,"|",4x,"BULK",4x,"|",
     :        3x,"CHANNEL",2x,"|",3x,"FILTER",2x,"|",
     :        3x,"CENTER ",                       
     :        3x,"|",3x," TOTAL  ")')
        
      do ispec=1,nspec
        write(*,'(1x,a2,"  ion: ",1("|",1x,e10.5,1x),
     :  3("|",1x,f10.7,1x),"|",f8.3,1x)')
     ;    fion(ispec),cbulk(ispec),cchan(ispec),cfilt(ispec),
     :    czocc(ispec),ctot(ispec)
        write(2,'(1x,a2,"  ion: ",1("|",1x,e10.5,1x),
     :  3("|",1x,f10.7,1x),"|",f8.3,1x)')
     ;    fion(ispec),cbulk(ispec),cchan(ispec),cfilt(ispec),
     :    czocc(ispec),ctot(ispec)
      enddo      
      write(*,'(72("_"))')
      write(*,'("  Charge: ",4("|",1x,f10.5,1x))')
     ;    chbulk,chchan,chfilt,chzocc
      write(2,'(72("_"))')
      write(2,'("  Charge: ",4("|",f10.5,1x))')
     ;    chbulk,chchan,chfilt,chzocc
      
      write(*,'(72("_"))')
      write(2,'(72("_"))')
      
c #### Salt data ####################################################

        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(istr+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(*,'(3x,"    Move ","    Jump ","      In ",
     :             "   Out ","  Accepted In  ",
     :             " Accepted Out ")')
      write(2,'(72("_"))')
      write(2,'(3x,"    Move ","    Jump ","      In ",
     :             "   Out ","  Accepted In  ",
     :             " Accepted Out ")')
      do i=1,nspec
         if (chonly(i)) then
            write(*,'(a2,":",1(2x,f7.5))')fion(i), amove(2,i)/amove(1,i)
            write(2,'(a2,":",1(2x,f7.5))')fion(i), amove(2,i)/amove(1,i)
         else
            write(*,'(a2,":",4(2x,f7.5),2(2x,i9,3x))')
     :      fion(i), amove(2,i)/amove(1,i),
     :      ajump(2,i)/ajump(1,i),ajin(2,i)/ajin(1,i),
     :      ajout(2,i)/ajout(1,i),int(ajin(2,i)),int(ajout(2,i))
            write(2,'(a2,":",4(2x,f7.5),2(2x,i9,3x))')
     :      fion(i), amove(2,i)/amove(1,i),
     :      ajump(2,i)/ajump(1,i),ajin(2,i)/ajin(1,i),
     :      ajout(2,i)/ajout(1,i),int(ajin(2,i)),int(ajout(2,i))
         endif
      enddo

      do icount=1,nsalt         
        do ireg=1,4
          if (ratreg(ireg).gt.1.d-5) then
          write(*,'(72("_"))')
          write(*,'("  Create ",a4," in region ",i1," :  "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),ireg,
     :      acreat(ireg,icount,2),acreat(ireg,icount,1),
     :      acreat(ireg,icount,2)/acreat(ireg,icount,1)
          write(*,'(" Destroy ",a4," in region ",i1," :  "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),ireg,
     :      adest(ireg,icount,2),adest(ireg,icount,1),
     :      adest(ireg,icount,2)/adest(ireg,icount,1)
          write(2,'(72("_"))')
          write(2,'("  Create ",a4," in region ",i1," :  "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),ireg,
     :      acreat(ireg,icount,2),acreat(ireg,icount,1),
     :      acreat(ireg,icount,2)/acreat(ireg,icount,1)
          write(2,'(" Destroy ",a4," in region ",i1," :  "
     :            ,f10.0," / ",f10.0," = ",f7.5)')
     :      fsalt(icount),ireg,
     :      adest(ireg,icount,2),adest(ireg,icount,1),
     :      adest(ireg,icount,2)/adest(ireg,icount,1)
          endif
        enddo
      enddo
      
      write(2,'(72("_"))')
      write(*,'(72("_"))')
        
c ##### Occupancies ##################################################

      open(12,file=focc)
         if (nfree.eq.2) then

            do ireg=1,3
              do i=0,4
              do j=0,4
                 k=0  
                 m=0  
                 aocci=aocc(ireg,i,j,k,m)
                 aocci=aocci/ataccu
                 write(12,'(2(i1,1x),f15.10)')i,j,aocci
              enddo
              enddo
              write(12,*)
            enddo     
            write(12,'(72("_"))')
                 
            do ireg=1,3      
              do j=0,4
                azcl=0.0
                do i=0,4
                   k=0  
                   m=0  
                   azcl=azcl+aocc(ireg,i,j,k,m)/ataccu
                enddo
                write(12,'(1(i1,1x),f15.10)')j,azcl  
              enddo
              write(12,*)
            enddo      

         elseif (nfree.eq.3) then

            do ireg=1,3
              do i=0,4
              do j=0,4
              do k=0,4
                 m=0  
                 aocci=aocc(ireg,i,j,k,m)
                 aocci=aocci/ataccu
                 write(12,'(3(i1,1x),f15.10)')i,j,k,aocci
              enddo
              enddo
              enddo
              write(12,*)
            enddo     
            write(12,'(72("_"))')
      
            do ireg=1,3      
              do j=0,4
              do k=0,4
                azcl=0.0
                do i=0,4
                   m=0  
                   azcl=azcl+aocc(ireg,i,j,k,m)/ataccu
                enddo
                write(12,'(2(i1,1x),f15.10)')j,k,azcl  
              enddo
              enddo
              write(12,*)
            enddo      
                 
         elseif (nfree.eq.4) then

            do ireg=1,3
              do i=0,4
              do j=0,4
              do k=0,4
              do m=0,4
                 aocci=aocc(ireg,i,j,k,m)
                 aocci=aocci/ataccu
                 write(12,'(4(i1,1x),f15.10)')i,j,k,m,aocci
              enddo
              enddo
              enddo
              enddo
              write(12,*)
            enddo     
            write(12,'(72("_"))')
      
            do ireg=1,3      
              do j=0,4
              do k=0,4
              do m=0,4
                azcl=0.0
                do i=0,4
                   azcl=azcl+aocc(ireg,i,j,k,m)/ataccu
                enddo
                write(12,'(3(i1,1x),f15.10)')j,k,m,azcl  
              enddo
              enddo
              enddo
              write(12,*)
            enddo      

         endif
      close(12)
      
      close(2)

      open(2,file=fh)
      do i=1,npatch
        write(2,*)i,ah(i)/ataccu
      enddo
      close(2)
      close(32)

      open(2,file=fdata,form='unformatted')
      write(2)d,dd,d2,z,q,ratmov,ratexc,eps,qstar,itry,chonly,
     :        nspec,istr,fion
      write(2)zl,rl,rlvest,rlcurv,epsw,epspr,deps,npatch
      write(2)drmax,vin,vout 
      write(2)chemp,ctarg,vgc,chex,ctrgcl,ratgr,ispcin,igcval,isalt,
     :        ication,nsalt,niter,fsalt
      write(2)athist,ataccu,dzg,drg,nzg,nrg,nrgr,nrgz
      write(2)an,ah,anin,cbulk,zocc,amove,ajump,ajin,ajout,acreat,adest
      write(2)nstep,naver,isave
      write(2)zbulk1,zbulk2,rbulk,vbulk,anaca,ksub
      write(2)nfree,drdf,nrdf,calrdf
      write(2)cfilt,cchan,ctot,czocc   
      write(2)chbulk,chchan,chfilt,chzocc
      close(2)

      return
      end
      
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      
      subroutine zero
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bulk/  zbulk1,zbulk2,rbulk,vbulk,anaca(nspcmx,4,20000),
     :              asuba(nspcmx,4,20000),ksub
      common/occ/   aocc(3,0:4,0:4,0:4,0:4),nfree
      common/radf/  ardf(3,nspcmx,nspcmx),
     :              rdf(3,nspcmx,nspcmx,nrgmx),
     :              drdf,nrdf,calrdf
      logical       chonly
      logical       calrdf,calgin

c_______ ZERO OUT STAFF _____________________________________________ 
     
      do i=1,3
       do ispec=1,nspec
        do jspec=1,nspec
          ardf(i,ispec,jspec)=0.0
          do j=1,nrdf
            rdf(i,ispec,jspec,j)=0.0
          enddo
        enddo
       enddo
      enddo  

      do i=0,4
       do j=0,4
        do k=0,4
         do m=0,4
           do ireg=1,3
             aocc(ireg,i,j,k,m)=0.0
           enddo  
         enddo
        enddo
       enddo
      enddo

      do j=1,4
       do ispec=nspec+1,istr+4
        nin(j,ispec)=0
       enddo
      enddo

      do ispec=1,nspec
       
       amove(1,ispec)=0.0
       amove(2,ispec)=0.0
       ajump(1,ispec)=0.0
       ajump(2,ispec)=0.0
       ajin(1,ispec)=0.0
       ajin(2,ispec)=0.0
       ajout(1,ispec)=0.0
       ajout(2,ispec)=0.0
       anin(1,ispec)=0.0
       anin(2,ispec)=0.0
       anin(3,ispec)=0.0
       anin(4,ispec)=0.0
     
       an(ispec)=0.0
     
         do i=1,nrgz
            do j=1,nrgr   
             gin(ispec,i,j) = 0.0
            enddo
         enddo   
       
         ii=0      
         do j=1,nzreg(ispec)
           do i=1,nzg(j,ispec)
             ii=ii+1
             gz(ispec,ii) = 0.0
           enddo
         enddo  

      enddo
     
      do ireg=1,4
        do icount=1,3
         acreat(ireg,icount,1)=0.0
         acreat(ireg,icount,2)=0.0
         adest(ireg,icount,1)=0.0
         adest(ireg,icount,2)=0.0
        enddo 
      enddo 
      
      athist = 0.0
      ataccu = 0.0
      ksub   = 0
      
      do i=1,npatch
        ah(i)=0.0
      enddo
     
      return
      end
      
c ####################################################################  
c ####################################################################

      
      subroutine ludcmp(a,n,indx,d)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      parameter(tiny0=1.0d-20)
      dimension a(npchmx,npchmx), vv(npchmx), indx(npchmx)
      d=1.0
      do i=1,n
            aamax=0.0
            do j=1,n
              if (dabs(a(i,j)).gt.aamax) aamax=dabs(a(i,j))
            enddo
            if (aamax.eq.0.0) pause 'singular matrix in ludcmp'
            vv(i)=1.0/aamax
      enddo
      do j=1,n
        do i=1,j-1
          sum0=a(i,j)
          do k=1,i-1
            sum0=sum0-a(i,k)*a(k,j)
          enddo
          a(i,j)=sum0
        enddo
        aamax=0.0
        do i=j,n
          sum0=a(i,j)
          do k=1,j-1
            sum0=sum0-a(i,k)*a(k,j)
          enddo
          a(i,j)=sum0
          dum=vv(i)*dabs(sum0)
          if (dum.ge.aamax) then
             imax=i
             aamax=dum
          endif
        enddo
        if (j.ne.imax) then
           do k=1,n
             dum=a(imax,k)
             a(j,k)=dum
           enddo
           d=-d
           vv(imax)=vv(j)
        endif
        indx(j)=imax
        if (a(j,j).eq.0) a(j,j)=tiny0
        if (j.ne.n) then
           dum=1.0/a(j,j)
           do i=j+1,n
             a(i,j)=a(i,j)*dum
           enddo
        endif
      enddo
      return
      end

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

      subroutine lubksb(a,n,indx,b)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      parameter(tiny0=1.0d-20)
      dimension a(npchmx,npchmx), b(npchmx), indx(npchmx)
      ii=0
      do i=1,n
        ll=indx(i)
        sum0=b(ll)
        b(ll)=b(i)
        if (ii.ne.0) then
           do j=ii,i-1
             sum0=sum0-a(i,j)*b(j)
           enddo
        elseif (sum0.ne.0) then
           ii=i
        endif
        b(i)=sum0
      enddo
      do i=n,1,-1
        sum0=b(i)
        do j=i+1,n
          sum0=sum0-a(i,j)*b(j)
        enddo
        b(i)=sum0/a(i,i)
      enddo
      return
      end


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

      subroutine calv
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/profl/ athist,ataccu,dzg(10,nspcmx),drg, 
     :              vjz(nspcmx,nzgmx),vjin(nspcmx,nrgmx,nrgmx),
     :              rgz(nspcmx,nzgmx),rgin(nspcmx,nrgmx),
     :              gz(nspcmx,nzgmx),gin(nspcmx,nrgmx,nrgmx),
     :              zinlft,rinup,zzreg(0:10,nspcmx),dzgi(nspcmx,nzgmx),
     :              nzg(10,nspcmx),nrg,nrgr,nrgz,nzgtot(nspcmx),
     :              nzreg(nspcmx),calgin
      logical       chonly,choold
      logical       calrdf,calgin
      dimension     rzi(nzgmx), r2i(nrgmx) 
      logical       ovrlap

      if (calgin) then 
        do ispec=1,nspec
          
          do i=1,nrgz
            do j=1,nrgr
               rlow   = (j-1)*drg                            
               rup    =     j*drg                             
               vjin(ispec,i,j) = drg*pi*(rup**2-rlow**2)
            enddo
          enddo  
        
        enddo
      endif  
        
        open(2,file="vjz.dat")
        
        do ispec=1,nspec
         
         choold=chonly(ispec)
         chonly(ispec)=.false.

         ri=d(ispec)/2.0
         ii=0

c________ region 1 ___________________________________________________         

         do i=1,nzg(1,ispec)
            ii=ii+1
            vjz(ispec,ii) = dzg(1,ispec)*pi*rl(5)**2
            rgz(ispec,ii) = zzreg(0,ispec)+(i-1)*dzg(1,ispec)
     :                     +dzg(1,ispec)/2.0
            write(2,*)rgz(ispec,ii),vjz(ispec,ii)
            dzgi(ispec,ii)=dzg(1,ispec)
         enddo   
         
c________ region 2 ___________________________________________________         

          do i=1,nzg(2,ispec)
             ii=ii+1
             rgz(ispec,ii) = 
     :               zzreg(1,ispec)+(i-1)*dzg(2,ispec)+dzg(2,ispec)/2.0
             zlow = rgz(ispec,ii)-dzg(2,ispec)/2.0
             zup  = zlow+dzg(2,ispec)
c--------- MC -------------------------------------------------------- 
             vjtot=dzg(2,ispec)*pi*rl(2)**2
             accmc=0.0
             do imc=1,10000 
                 rsqnew = ranff(mag) * rl(2)*rl(2)
                 r2inew = dsqrt(rsqnew)
                 finew  = ranff(mag) * twopi
                 rxinew = r2inew * dcos(finew)      
                 ryinew = r2inew * dsin(finew)
                 rannum = ranff(mag)
                 rzinew = zlow+ranff(mag)*(zup-zlow)
                 drznew = dabs(rzinew)
                 call wall(ispec,drznew,r2inew,ovrlap)
                 if (.not.ovrlap) accmc=accmc+1.0
                 ratio = accmc/dble(imc)
             enddo
             vjz(ispec,ii) = ratio*vjtot
             write(2,*)rgz(ispec,ii),vjz(ispec,ii)
             dzgi(ispec,ii)=dzg(2,ispec)
          enddo   
         
c________ region 3 ___________________________________________________         

              do i=1,nzg(3,ispec)
                 ii=ii+1
                 rgz(ispec,ii) = 
     :               zzreg(2,ispec)+(i-1)*dzg(3,ispec)+dzg(3,ispec)/2.0
                 vjz(ispec,ii) = dzg(3,ispec)*pi*(rl(1)-ri)**2
                 write(2,*)rgz(ispec,ii),vjz(ispec,ii)
                 dzgi(ispec,ii)=dzg(3,ispec)
              enddo   
         
c________ region 4 ___________________________________________________         

          do i=1,nzg(4,ispec)
             ii=ii+1
             rgz(ispec,ii) = 
     :               zzreg(3,ispec)+(i-1)*dzg(4,ispec)+dzg(4,ispec)/2.0
             zlow = rgz(ispec,ii)-dzg(4,ispec)/2.0
             zup  = zlow+dzg(4,ispec)
c--------- MC -------------------------------------------------------- 
             vjtot=dzg(4,ispec)*pi*rl(2)**2
             accmc=0.0
             do imc=1,10000 
                 rsqnew = ranff(mag) * rl(2)*rl(2)
                 r2inew = dsqrt(rsqnew)
                 finew  = ranff(mag) * twopi
                 rxinew = r2inew * dcos(finew)      
                 ryinew = r2inew * dsin(finew)
                 rannum = ranff(mag)
                 rzinew = zlow+ranff(mag)*(zup-zlow)
                 call wall(ispec,rzinew,r2inew,ovrlap)
                 if (.not.ovrlap) accmc=accmc+1.0
                 ratio = accmc/dble(imc)
             enddo
             vjz(ispec,ii) = ratio*vjtot
             write(2,*)rgz(ispec,ii),vjz(ispec,ii)
             dzgi(ispec,ii)=dzg(4,ispec)
          enddo   
         
c________ region 5 ___________________________________________________         

         do i=1,nzg(5,ispec)
            ii=ii+1
            vjz(ispec,ii) = dzg(5,ispec)*pi*rl(5)**2
            rgz(ispec,ii) = zzreg(4,ispec)+(i-1)*dzg(5,ispec)
     :                     +dzg(5,ispec)/2.0
            write(2,*)rgz(ispec,ii),vjz(ispec,ii)
            dzgi(ispec,ii)=dzg(5,ispec)
         enddo   
        
        write(2,*)
       
       chonly(ispec)=choold

       enddo
        
      close(2)
         
      return
      end
c_______________________________________________________________________
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      subroutine move(ispec,i)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical ovrlap
      ovrlap=.false.
c      write(*,*)"move"	

      ii=indspc(ispec,i)
      rxiold = rx(ii)
      ryiold = ry(ii)
      rziold = rz(ii)
      r2iold = r2(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
           if (rzold.lt.zl(2)+ri) then
            drmaxi=drmax(1)*di
           else 
            drmaxi=drmax(2)*di
           endif
           rxinew = rxiold + (2.0*ranff(mag)-1.0) * drmaxi
           ryinew = ryiold + (2.0*ranff(mag)-1.0) * drmaxi
           rzinew = rziold + (2.0*ranff(mag)-1.0) * drmaxi
           drznew = dabs(rzinew)
           r2inew = dsqrt(rxinew*rxinew+ryinew*ryinew)
           call wall(ispec,drznew,r2inew,ovrlap)
           if (ovrlap) return  
           imove=1

      else
      
           ajump(1,ispec) = ajump(1,ispec)+1
c  %%%%%%%%% Generate new position in a cylinder %%%%%%%%%%%%%%
           if (chonly(ispec)) then
               rl1=rl(1)-ri
               zl1=zl(1)-ri
           else
               rl1=rl(5)
               zl1=zl(4)
           endif
           rsqnew = ranff(mag) * rl1*rl1
           r2inew = dsqrt(rsqnew)
           finew  = ranff(mag) * twopi
           rxinew = r2inew * dcos(finew)      
           ryinew = r2inew * dsin(finew)
           rzinew = (2*ranff(mag)-1)*zl1
           drznew = dabs(rzinew)
           if (.not.chonly(ispec)) then
              call wall(ispec,drznew,r2inew,ovrlap)
              if (ovrlap) return   
           endif
           imove=0 

      endif
     
      call ergi(ispec,i,rxinew,ryinew,rzinew,r2inew,dui,dup,ovrlap)
      if (ovrlap) return 
      
      deltu=dui+dup

      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
            r2(ii) = r2inew
            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

            do ireg=1,3
               if ((rziold.lt.zreg1(ireg,ispec).or.
     :              rziold.gt.zreg2(ireg,ispec)).and.
     :             (rzinew.gt.zreg1(ireg,ispec).and.
     :              rzinew.lt.zreg2(ireg,ispec)))
     :               call addreg(ireg,ispec,ii)
               if ((rzinew.lt.zreg1(ireg,ispec).or.
     :              rzinew.gt.zreg2(ireg,ispec)).and.
     :             (rziold.gt.zreg1(ireg,ispec).and.
     :              rziold.lt.zreg2(ireg,ispec)))
     :               call delreg(ireg,ispec,ii)
            enddo

            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
               rip(ii,j)=ripnew(ispec,j)
              enddo
            endif

      else

            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
                h(j)=hold(j)
                c(j)=cold(j)
              enddo  
            endif  

      endif

      return
      end

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

      subroutine jmpin(ispec,i)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical ovrlap
      ovrlap=.false.
c      write(*,*)" jmpin"

      ii=indspc(ispec,i)
      ajin(1,ispec) = ajin(1,ispec)+1

      di=d(ispec)
      ri=di/2.0

      rsqnew = ranff(mag) * rljpin(ispec)*rljpin(ispec)
      r2inew = dsqrt(rsqnew)
      finew  = ranff(mag) * twopi
      rxinew = r2inew * dcos(finew)      
      ryinew = r2inew * dsin(finew)
      rzinew = (2*ranff(mag)-1)*zljpin(ispec)
      drznew = dabs(rzinew)
      call wall(ispec,drznew,r2inew,ovrlap)
      if (ovrlap) return 
      
      call ergi(ispec,i,rxinew,ryinew,rzinew,r2inew,dui,dup,ovrlap)
      if (ovrlap) return 
      
      deltu=dui+dup
     
c      write(*,*)"deltu in jmpin = ",dui,dup
     
      expbu = vin(ispec)*dexp(-deltu)/vout(ispec)
      if (expbu.gt.ranff(mag)) then
            ajin(2,ispec) = ajin(2,ispec) + 1.0
            rziold = rz(ii)
            rx(ii) = rxinew
            ry(ii) = ryinew
            rz(ii) = rzinew
            r2(ii) = r2inew
            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
            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
                rip(ii,j)=ripnew(ispec,j)
              enddo
            endif  

            do ireg=1,3
               if ((rziold.lt.zreg1(ireg,ispec).or.
     :              rziold.gt.zreg2(ireg,ispec)).and.
     :             (rzinew.gt.zreg1(ireg,ispec).and.
     :              rzinew.lt.zreg2(ireg,ispec)))
     :               call addreg(ireg,ispec,ii)
               if ((rzinew.lt.zreg1(ireg,ispec).or.
     :              rzinew.gt.zreg2(ireg,ispec)).and.
     :             (rziold.gt.zreg1(ireg,ispec).and.
     :              rziold.lt.zreg2(ireg,ispec)))
     :               call delreg(ireg,ispec,ii)
            enddo

      else
            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
                h(j)=hold(j)
                c(j)=cold(j)
              enddo  
            endif  
      endif

      return
      end

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

      subroutine jmpout(ispec,i)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/moove/ drmax(2),vin(nspcmx),vout(nspcmx),rljpin(nspcmx),
     :              zljpin(nspcmx) 
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical ovrlap
      ovrlap=.false.
c      write(*,*)" jmpout"

      ii=indspc(ispec,i)
      ajout(1,ispec) = ajout(1,ispec)+1

      di=d(ispec)
      ri=di/2.0
      rziold=rz(ii)

      rsqnew = ranff(mag) * rl(5)*rl(5)
      r2inew = dsqrt(rsqnew)
      finew  = ranff(mag) * twopi
      rxinew = r2inew * dcos(finew)      
      ryinew = r2inew * dsin(finew)
      rzinew = (2*ranff(mag)-1)*(zl(4)-zljpin(ispec))
      if (rzinew.lt.0.0) then
           rzinew = -zljpin(ispec)+rzinew
      else
           rzinew =  zljpin(ispec)+rzinew
      endif     
      
      call ergi(ispec,i,rxinew,ryinew,rzinew,r2inew,dui,dup,ovrlap)
      if (ovrlap) return 
      
      deltu=dui+dup
     
      expbu = vout(ispec)*dexp(-deltu)/vin(ispec)
      if (expbu.gt.ranff(mag)) then
            ajout(2,ispec) = ajout(2,ispec) + 1.0
            rzold  = dabs(rz(ii))
            rx(ii) = rxinew
            ry(ii) = ryinew
            rz(ii) = rzinew
            drznew = dabs(rzinew)
            r2(ii) = r2inew
            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
            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
                rip(ii,j)=ripnew(ispec,j)
              enddo
            endif  

            do ireg=1,3
               if ((rziold.lt.zreg1(ireg,ispec).or.
     :              rziold.gt.zreg2(ireg,ispec)).and.
     :             (rzinew.gt.zreg1(ireg,ispec).and.
     :              rzinew.lt.zreg2(ireg,ispec)))
     :               call addreg(ireg,ispec,ii)
               if ((rzinew.lt.zreg1(ireg,ispec).or.
     :              rzinew.gt.zreg2(ireg,ispec)).and.
     :             (rziold.gt.zreg1(ireg,ispec).and.
     :              rziold.lt.zreg2(ireg,ispec)))
     :               call delreg(ireg,ispec,ii)
            enddo

      else
            if (dabs(deps).gt.1.d-5) then
              do j=1,npatch
                h(j)=hold(j)
                c(j)=cold(j)
              enddo  
            endif  
      endif

      return
      end

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

      subroutine ergi(ispec,i,rxi,ryi,rzi,r2i,uii,uip,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical ovrlap
      ovrlap=.false.
     
      uii=0.0
      uip=0.0
      
      qi=q(ispec)
      di=d(ispec)
      ri=di/2.0
      ii=indspc(ispec,i)
      epsi=eps(ispec)
      
c    ###################################################################
c_______ SUM OVER OTHER FREE CHARGES ________________________________
c    ###################################################################
      
      do jspec=1,nspec
            ddi=d2(ispec,jspec)
            qj=q(jspec)
            nj=n(jspec)
            epsj = eps(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
               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*(1.0/epsj+1.0/epsi)/(2*rij)   
               uii=uii + uiinew
c_______ DISTANCE IS FROM LOOKUP TABLE ______________________________
c_______ FOR OLD CONFIGURATION ______________________________________
               rij=rii(ii,jj)         
               uiiold = qi*qj*(1.0/epsj+1.0/epsi)/(2*rij)   
               uii=uii - uiiold   
1              continue               
            enddo
      enddo         
   
      if (dabs(deps).lt.1.d-5) return
     
c    ###################################################################
c_______ SUM OVER PATCHES ___________________________________________
c    ###################################################################
      
       call calch(ispec,i,rxi,ryi,rzi)
       do k=1,npatch
           do jspec=1,nspec
              nj=n(jspec)
              qj=q(jspec)
              do j=1,nj
                jj=indspc(jspec,j)
                if (ii.eq.jj) then
c_______ NEW CONFIGURATION: FROM RNEW CALCULATED IN CALCH ___________
                   rij = ripnew(ispec,k)      
                   uip = uip + qj*h(k)*a(k)/(2*rij)   
c_______ OLD CONFIGURATION: LOOKUP TABLE ____________________________
                   rij = rip(jj,k)         
                   uip = uip - qj*hold(k)*a(k)/(2*rij)
                else  
c_______ IONS NOT MOVED: LOOKUP TABLE _______________________________
                   rij = rip(jj,k)
                   uip = uip + qj*dh(k)*a(k)/(2*rij)
                endif 
              enddo
           enddo
       enddo
      
      return
      end

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

      subroutine calch(ispec,i,rxi,ryi,rzi)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
       
c_______ CALCULATE AND STORE DISTANCES BETWEEN _____________________
c_______ PATCHES AND PARTICLE II (RETURN IF OVERLAP) _______________
        
      qi=q(ispec)
      di=d(ispec)
      ri=di/2.0
      epsi=eps(ispec)
      
c_______ STORE OLD H AND C __________________________________________
       
      do k=1,npatch
       hold(k)=h(k)
       cold(k)=c(k)
      enddo
      
c_______ CALCULATE VECTOR C _________________________________________
      
      ii=indspc(ispec,i)
      rxiold=rx(ii)
      ryiold=ry(ii)
      rziold=rz(ii)
      
      do k=1,npatch
        
        rxk = px(k)
        ryk = py(k)
        rzk = pz(k)
c ##### NEW ########################################################## 
        rxki = rxk-rxi
        ryki = ryk-ryi
        rzki = rzk-rzi
        rijsq = rxki*rxki+ryki*ryki+rzki*rzki
        rij = dsqrt(rijsq)
        ripnew(ispec,k)=rij
        rzki  = rxki*ux(k)+ryki*uy(k)+rzki*uz(k)
        dc(k) = -deps*qi*rzki/(4*pi*epsi*rij*rij*rij)
c ##### OLD ########################################################## 
        rij   = rip(ii,k)
        rzki  = (rxk-rxiold)*ux(k)+(ryk-ryiold)*uy(k)+(rzk-rziold)*uz(k)
        dc(k) = dc(k)+deps*qi*rzki/(4*pi*epsi*rij*rij*rij)
        c(k)  = cold(k)+dc(k)
        h(k)  = c(k)
      enddo
      
c_______ CALCULATE VECTOR H _________________________________________

      call lubksb(amx,npatch,indx,h)
     
      do k=1,npatch
       dh(k)=h(k)-hold(k)
      enddo
      
      return
      end

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

      subroutine create(ireg,igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      character*4   fsalt
      logical       chonly
      logical ovrlap
      ovrlap=.false.
c      write(*,*)"create ",igc, " in region ",ireg
      
      acreat(ireg,igc,1) = acreat(ireg,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)
      if (ireg.eq.4) then
        dnterm=dble(n(ispec1)+1)
      else
        dnterm=dble(nin(ireg,ispec1)+1)
      endif
      do icount=2,iv+1
        ispec=ispcin(igc,icount)
        dnterm=dnterm*dble(n(ispec)+icount-1)
      enddo  

      do icount=1,iv+1
         ispec  = ispcin(igc,icount)
         ri     = d(ispec)/2.0
         rmaxsq = rcr(ireg,ispec)**2
         rsqnew = ranff(mag) * rmaxsq 
         r2inew = dsqrt(rsqnew)
         finew  = ranff(mag) * twopi
         rxinew = r2inew * dcos(finew)      
         ryinew = r2inew * dsin(finew)
         rzinew = zcr1(ireg,ispec)+ranff(mag)
     :            *(zcr2(ireg,ispec)-zcr1(ireg,ispec))
         drznew = dabs(rzinew)
         call wall(ispec,drznew,r2inew,ovrlap)
         if (ovrlap) return   
         rznew(icount) = rzinew                
         rxnew(icount) = rxinew            
         rynew(icount) = ryinew           
         r2new(icount) = r2inew
      enddo
    
      call encrea(igc,uii,uip,ovrlap)
      if (ovrlap) then
        return
      endif
      
      deltu=uii+uip
      bfcr = dexp(-deltu+chemp(igc))*vgc(ireg,igc)/dnterm
      
      if (bfcr.gt.ranff(mag)) then
         call add(igc)
         acreat(ireg,igc,2) = acreat(ireg,igc,2)+1.0
      else
          if (dabs(deps).gt.1.d-5) then
            do j=1,npatch
              h(j)=hold(j)
              c(j)=cold(j)
            enddo
          endif  
      endif

      return
      end

c    ###################################################################
c    ###################################################################
c    ###################################################################
c    ###################################################################
      
      subroutine destr(ireg,igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),ajin(2,nspcmx), 
     :              ajout(2,nspcmx),acreat(4,3,2),adest(4,3,2)
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      character*4   fsalt
      logical       chonly
c      write(*,*)"delete ",igc, " in region ",ireg
      
      adest(ireg,igc,1) = adest(ireg,igc,1) + 1.0
      
      iv=igcval(igc)
      
      ispec1=ispcin(igc,1)
      if (ireg.eq.4) then
        dnterm=dble(n(ispec1))
      else
        dnterm=dble(nin(ireg,ispec1))
      endif
      if (nterm.eq.0) return
      do icount=2,iv+1
        ispec=ispcin(igc,icount)
        dnterm=dnterm*dble(n(ispec)-icount+2)
      enddo  
      
      do icount=1,iv+1
        ispec  = ispcin(igc,icount)
1       continue
        if (ireg.eq.4) then
          i   = int(n(ispec)*ranff(mag)) + 1
          ii  = indspc(ispec,i)                 
        else
          if (ispec.ne.istr+1) then
              j   = int(nin(ireg,ispec)*ranff(mag)) + 1
c              write(*,*)" j in destroy",j,nin(ireg,ispec)
              ii  = indreg(ireg,ispec,j)
              i   = ispcbk(ii)
              ii2 = indspc(ispec,i)
              if (ii.ne.ii2) stop "Gebasz in create"
          else
              i   = int(n(ispec)*ranff(mag)) + 1
              ii  = indspc(ispec,i)                 
          endif
        endif
        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)
        r2new(icount) = r2(ii)
      enddo
      
      call endest(igc,uii,uip,ovrlap)
      
      deltu=uii+uip
      bfdr = dexp(-deltu-chemp(igc))*dnterm/vgc(ireg,igc)
         
      if (bfdr.gt.ranff(mag)) then
         call remove(igc)
         adest(ireg,igc,2) = adest(ireg,igc,2)+1.0
      else
          if (dabs(deps).gt.1.d-5) then
            do j=1,npatch
              h(j)=hold(j)
              c(j)=cold(j)
            enddo
          endif  
      endif
      
      return
      end

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

      subroutine remove(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly

      iv=igcval(igc)
  
      do icount=1,iv+1
        ispec=ispcin(igc,icount)  
        i=iout(icount)
        ii = indspc(ispec,i)
        rzinew=rznew(icount)
        drznew=dabs(rznew(icount))
        ri=d(ispec)/2.0
        do ireg=1,3
           if (rzinew.gt.zreg1(ireg,ispec).and.
     :         rzinew.lt.zreg2(ireg,ispec))
     :         call delreg(ireg,ispec,ii)
        enddo
        if (i.lt.n(ispec)) then
          do k = i+1,n(ispec)
c             iprev  = indspc(ispec,k-1)
             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=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/geom/  zl(6),rl(6),rlvest,rlcurv,epsw,epspr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      dimension     inew(3)
   
      iv=igcval(igc)
   
      do icount=1,iv+1
          ispec=ispcin(igc,icount)  
          inew(icount) = n(ispec) + 1
          ii = indspc(ispec,inew(icount))
          rzinew=rznew(icount)
          drznew=dabs(rznew(icount))
          ri=d(ispec)/2.0
          if (ii.eq.0) then
            ii = (ispec-1)*nionmx+inew(icount)
            indspc(ispec,inew(icount)) = ii 
            ispcbk(ii)=inew(icount)
          endif
          do ireg=1,3
               if (rzinew.gt.zreg1(ireg,ispec).and.
     :             rzinew.lt.zreg2(ireg,ispec))
     :             call addreg(ireg,ispec,ii)
          enddo
          n(ispec) = inew(icount)
          rx(ii)=rxnew(icount)
          ry(ii)=rynew(icount)
          rz(ii)=rznew(icount)
          r2(ii)=r2new(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
        if (dabs(deps).gt.1.d-5) then
          do k=1,npatch
            rip(ii,k)=ripnew(icount,k)
          enddo
        endif  
      enddo
      
      return
      end

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

      subroutine encrea(igc,uii,uip,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
      logical ovrlap
      ovrlap=.false.
      
      uii=0.0
      uip=0.0
   
      iv=igcval(igc)
   
      do icount1=1,iv
        ispec1=ispcin(igc,icount1)
        do icount2=icount1+1,iv+1
          ispec2=ispcin(igc,icount2)
          ddi=d2(ispec1,ispec2)
          rxij = rxnew(icount1)-rxnew(icount2)
          ryij = rynew(icount1)-rynew(icount2)
          rzij = rznew(icount1)-rznew(icount2)
          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 + q(ispec1)*q(ispec2)*(1.0/eps(ispec1) 
     :                                  +1.0/eps(ispec2))/(2.0*rij) 
        enddo     
      enddo     
      
      do icount=1,iv+1
        ispec=ispcin(igc,icount)
        qi=q(ispec)
        epsi=eps(ispec)
        do jspec=1,nspec
          qj=q(jspec)
          epsj=eps(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)
             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*(1.0/epsi+1.0/epsj)/(2.0*rij)
          enddo 
        enddo
      enddo 
      
      if (dabs(deps).lt.1.d-5) return
     
      call chcrea(igc)
      do k=1,npatch
        do jspec=1,nspec
          nj=n(jspec)
          qj=q(jspec)
          do j=1,nj
             jj=indspc(jspec,j)
             rij = rip(jj,k)
             uip = uip + qj*dh(k)*a(k)/(2.0*rij)
          enddo
        enddo
        do icount=1,iv+1
          ispec=ispcin(igc,icount)
          qj  = q(ispec)
          rij = ripnew(icount,k)
          uip = uip + qj*h(k)*a(k)/(2.0*rij)
        enddo
      enddo

      return
      end

c###################################################################
c#####################################################################
c#####################################################################
      
      subroutine endest(igc,uii,uip,ovrlap)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly

      uii=0.0
      uip=0.0
      
      iv=igcval(igc)
      
      do icount=1,iv+1
        ispec=ispcin(igc,icount)
        i=iout(icount)
        ii=indspc(ispec,i)
        qi=q(ispec)
        epsi=eps(ispec)
        do jspec=1,nspec
          qj=q(jspec)
          epsj=eps(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*(1.0/epsi+1.0/epsj)/(2.0*rij) 
             endif
          enddo 
        enddo 
      enddo 
      
      do icount1=1,iv
        ispec1=ispcin(igc,icount1)
        i1=iout(icount1)
        ii1=indspc(ispec1,i1)
        do icount2=icount1+1,iv+1
          ispec2=ispcin(igc,icount2)
          i2=iout(icount2)
          ii2=indspc(ispec2,i2)
          rij=rii(ii1,ii2)
          uii=uii + q(ispec1)*q(ispec2)*(1.0/eps(ispec1) 
     :                                  +1.0/eps(ispec2))/(2.0*rij) 
        enddo     
      enddo     

      if (dabs(deps).lt.1.d-5) return
     
      call chdest(igc)
      
      do k=1,npatch
        do jspec=1,nspec
          nj=n(jspec)
          qj=q(jspec)
          do j=1,nj
             jj=indspc(jspec,j)
             rij = rip(jj,k)
             ieqj=0
             do icount=1,iv+1
               ispec=ispcin(igc,icount)
               i=iout(icount)
               ii=indspc(ispec,i)
               if (ii.eq.jj) then
                  uip = uip - qj*hold(k)*a(k)/(2.0*rij)
                  ieqj=1
               endif
             enddo  
             if (ieqj.eq.0) uip = uip + qj*dh(k)*a(k)/(2.0*rij)
          enddo
        enddo
      enddo
     
      return
      end

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

      subroutine chcrea(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
       
      iv=igcval(igc)

      do k=1,npatch
        hold(k)=h(k)
        cold(k)=c(k)
        dc(k)=0.0
      enddo

      do icount=1,iv+1
        ispec=ispcin(igc,icount)
        qi=q(ispec)
        di=d(ispec)
        epsi=eps(ispec)
        do j=1,npatch
          rxj   = px(j)
          ryj   = py(j)
          rzj   = pz(j)
          rxij  = rxj-rxnew(icount)
          ryij  = ryj-rynew(icount)
          rzij  = rzj-rznew(icount)
          rijsq = rxij*rxij+ryij*ryij+rzij*rzij
          rij = dsqrt(rijsq)
          ripnew(icount,j)=rij
          rk  = rxij*ux(j) + ryij*uy(j) + rzij*uz(j)
          dc(j)= dc(j)-deps*qi*rk/(4.0*pi*epsi*rij*rij*rij)
          c(j)=cold(j)+dc(j)
          h(j)=c(j)
        enddo
      enddo

c_______ CALCULATE VECTOR H _________________________________________

      call lubksb(amx,npatch,indx,h)
     
      do k=1,npatch
       dh(k)=h(k)-hold(k)
      enddo
      return
      end
      
c###################################################################
c#####################################################################
c#####################################################################

      subroutine chdest(igc)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/const/ pi,twopi,sq2,sq3,tosi,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),r2(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              z(nspcmx),q(nspcmx),
     :              ratmov(nspcmx),ratexc(nspcmx), 
     :              eps(nspcmx),qstar,itry(nspcmx),
     :              chonly(nspcmx),
     :              nspec,nwh(ntrymx),ntry,istr
      common/patch/ px(npchmx),py(npchmx),pz(npchmx),h(npchmx),
     :              ux(npchmx),uy(npchmx),uz(npchmx),a(npchmx),
     :              amx(npchmx,npchmx),c(npchmx),deps,
     :              rip(ntotmx,npchmx),
     :              npatch,indx(npchmx),npch(10)
      common/gc/    chemp(3),ctarg(3),vgc(4,3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,ratreg(4),
     :              zcr1(4,nspcmx),zcr2(4,nspcmx),rcr(4,nspcmx),
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/bufer/ hold(npchmx), cold(npchmx),ripnew(nspcmx,npchmx),
     :              riinew(nspcmx,ntotmx),dc(npchmx),dh(npchmx)
      logical       chonly
       
      iv=igcval(igc)
     
      do k=1,npatch
        hold(k)=h(k)
        cold(k)=c(k)
        dc(k)=0.0
      enddo

      do icount=1,iv+1
        ispec = ispcin(igc,icount)
        i     = iout(icount)
        ii    = indspc(ispec,i)
        qi=q(ispec)
        di=d(ispec)
        epsi=eps(ispec)
        do j=1,npatch
          rxj   = px(j)
          ryj   = py(j)
          rzj   = pz(j)
          rxij  = rxj-rxnew(icount)
          ryij  = ryj-rynew(icount)
          rzij  = rzj-rznew(icount)
          rij   = rip(ii,j)
          rk    = rxij*ux(j) + ryij*uy(j) + rzij*uz(j)
          dc(j) = dc(j)+deps*qi*rk/(4.0*pi*epsi*rij*rij*rij)
          c(j)  = cold(j)+dc(j)
          h(j)  = c(j)
        enddo
      enddo

c_______ CALCULATE VECTOR H _________________________________________

      call lubksb(amx,npatch,indx,h)
     
      do k=1,npatch
       dh(k)=h(k)-hold(k)
      enddo

      return
      end

c###################################################################
c###################################################################
c###################################################################
      
      subroutine addreg(ireg,ispec,i)
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
      
      inew = nin(ireg,ispec)+1
      nin(ireg,ispec) = inew
      indreg(ireg,ispec,inew) = i
c      write(*,*)"addreg ",ireg,ispec,inew,i

      return
      end

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

      subroutine delreg(ireg,ispec,ii)  
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=2000,ntotmx=12000,ntrymx=12000,npchmx=2000, 
     :          nspcmx=6,nzgmx=4000,nrgmx=500)
      common/accu/  an(nspcmx),ah(npchmx),anin(4,nspcmx),
     :              zreg1(4,nspcmx),zreg2(4,nspcmx),rreg(4,nspcmx), 
     :              cbulk(nspcmx),zocc,indreg(3,nspcmx,nionmx),
     :              nin(4,nspcmx) 
    
c      write(*,*)"delreg ",ireg,ispec,nin(ireg,ispec),ii
      do j=1,nin(ireg,ispec)
        i1=indreg(ireg,ispec,j)
c        write(*,*)ireg,ispec,nin(ireg,ispec),i,j,i1
        if (ii.eq.i1) then 
          j1=j
          goto 1
        endif
      enddo
      stop "Gebasz in delreg"
1     continue      

      if (j1.lt.nin(ireg,ispec)) then
          do k = j1+1,nin(ireg,ispec)
             indreg(ireg,ispec,k-1) = indreg(ireg,ispec,k) 
          enddo
      endif
     
      nin(ireg,ispec) = nin(ireg,ispec)-1

      return
      end
