      program doublelayer
      implicit doubleprecision (a-h,o-z)

      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      common/logi/   readit,contin,yesmem
      common/integ/  r(3,20000),c(3,20000),cbulk(3),cbr(10),daver 
      common/wrout/  qp1si,qp2si,qnsi,chg11,chg22,tstar,eps,zlr(10),
     :               sigsi,sleft,ctarg
      common/sub/    difpot(10000),ksub
      
      dimension x(3)
      logical readit, ovrlap, yesmem, contin
      integer step

      open(2,file='in')
      read(2,*)niter
      close(2)

      do iter=1,niter

      call readin(iter)    

      ntot = n(1)+n(2)+n(3)
      x(1) = dble(n(1))/dble(ntot)
      x(2) = dble(n(2))/dble(ntot)
      x(3) = dble(n(3))/dble(ntot)
      xp = x(1)+x(2)

      call zero
 
      write(*,*)' Start equilibration'
      isim=0
      do 1 step=1,naver
       do 2 ip=1,ntot
           which = ranff(mag)
           if (which.lt.x(1)) then
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(1)
              else
                 call jump(1)
              endif 
           elseif (which.gt.xp) then 
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(3) 
              else
                 call jump(3)
              endif 
           else
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(2)
              else
                 call jump(2)
              endif 
           endif
2      continue       
        if (mod(step,isave).eq.0) call wrcf 
1     continue
 
      call zero
      write(*,*)' Start simulation'    
      isim=1
      do 11 step=1,nstep
       do 12 ip=1,ntot
           which = ranff(mag)
           if (which.lt.x(1)) then
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(1)
              else
                 call jump(1)
              endif 
           elseif (which.gt.xp) then 
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(3) 
              else
                 call jump(3)
              endif 
           else
              flmv=ranff(mag)
              if (flmv.gt.0.5) then
                 call move(2)
              else
                 call jump(2)
              endif 
           endif
12     continue       
        if (mod(step,ihist).eq.0) call hist
        if (mod(step,isave).eq.0) call saves(iter,step)
        if (mod(step,isave).eq.0) call wrcf 
11     continue

      enddo

      write(*,*)' End of simulation'
      
      stop   
      end

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

      subroutine readin(iter)    
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      common/logi/   readit,contin,yesmem
      common/integ/  r(3,20000),c(3,20000),cbulk(3),cbr(10),daver 
      common/wrout/  qp1si,qp2si,qnsi,chg11,chg22,tstar,eps,zlr(10),
     :               sigsi,sleft,ctarg
      logical readit, ovrlap, yesmem, contin
      dimension chg1(2001),chg2(2001),chg3(2001)

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

c_______ READIN PARAMETERS __________________________________________
      open (2,file='in')
      read(2,*)niter
      read(2,*)readit
      read(2,*)contin
      read(2,*)
      read(2,*)tsi
      read(2,*)eps
      read(2,*)ctarg
      read(2,*)xl   
      read(2,*)sleft
      read(2,*)zl0
      read(2,*)       
      read(2,*)n(1)
      read(2,*)n(2)
      read(2,*)n(3)
      read(2,*)
      read(2,*)qp1si
      read(2,*)qp2si
      read(2,*)qnsi
      read(2,*)
      read(2,*)dp1
      read(2,*)dp2
      read(2,*)dn
      read(2,*)
      read(2,*)dz  
      read(2,*)drmxp1
      read(2,*)drmxp2
      read(2,*)drmxn
      read(2,*)
      read(2,*)nstep0 
      read(2,*)naver 
      read(2,*)isave
      read(2,*)ihist
      close(2)

      avog=6.02214d23
      tosi=1.0/(1.d-10)**3/avog/1000.
      dsi=1.d-10
      eps0=8.8542d-12
      boltz=1.3806d-23
      echg=1.6021917d-19
      beta=1./(boltz*tsi)

      sigsi=0.0
      sigma=(sigsi/echg)*dsi**2

      z(1)=qp1si
      z(2)=qp2si
      z(3)=qnsi

      rtarg=ctarg/tosi
c     v=dble(n(3))/rtarg
c     xl=v**(1./3.)
      ni=n(1)+n(2)+n(3)
      dq=-(n(1)*qp1si+n(2)*qp2si+n(3)*qnsi)
c      xl=dsqrt(dq/sigma)
      if (niter.ne.1) then
            if (iter.eq.1) then
                  zl=dble(n(3))/(rtarg*xl**2)
            else
                  zl=zl*cbulk(3)/ctarg
            endif
            if (iter.ne.niter) then
                  nstep=20000    
            else
                  nstep=nstep0
            endif
      else
            nstep = nstep0
            zl    = zl0
      endif      
      v=xl*xl*zl
      zlr(iter)=zl

c_______ CALCULATE PARTICLE DIAMETERS _______________________________
      daver=(n(1)*dp1+n(2)*dp2+n(3)*dn)/dble(ni)

      dp1p2=(dp1+dp2)/2.0
      dp1n=(dp1+dn)/2.0
      dp2n=(dp2+dn)/2.0
      
      drmxp1=drmxp1*dp1
      drmxp2=drmxp2*dp2
      drmxn =drmxn*dn

      d(1,1)=dp1
      d(1,2)=dp1p2
      d(1,3)=dp1n
      d(2,1)=dp1p2
      d(2,2)=dp2 
      d(2,3)=dp2n
      d(3,1)=dp1n
      d(3,2)=dp2n
      d(3,3)=dn

      d2(1,1)=dp1**2
      d2(1,2)=dp1p2**2
      d2(1,3)=dp1n**2
      d2(2,1)=dp1p2**2
      d2(2,2)=dp2**2 
      d2(2,3)=dp2n**2
      d2(3,1)=dp1n**2
      d2(3,2)=dp2n**2
      d2(3,3)=dn**2

c_______ REDUCE CHARGE AND DIPOLE MOMENT ____________________________ 
c      dp1si=dp1si*1.d-10
c      dp2si=dp2si*1.d-10
c      dnsi=dnsi*1.d-10
      
      
      qp1=qp1si*echg
      qp2=qp2si*echg
      qn=qnsi*echg

      qp1sq=beta*qp1*qp1/(4*pi*eps0*eps*dsi)
      qp1=dsqrt(qp1sq)
      qp2sq=beta*qp2*qp2/(4*pi*eps0*eps*dsi)
      qp2=dsqrt(qp2sq)
      qnsq=beta*qn*qn/(4*pi*eps0*eps*dsi)
      qn=-dsqrt(qnsq)
      q(1)=qp1
      q(2)=qp2
      q(3)=qn 
      tstar=4*pi*eps0*eps*dsi/(beta*echg**2)

      qleft=sleft*dq
      qright=dq-qleft
     
      qlsq=beta/(4*pi*eps0*eps*dsi)
      qleft=qleft*dsqrt(qlsq)*echg
      qrsq=beta/(4*pi*eps0*eps*dsi)
      qright=qright*dsqrt(qrsq)*echg
      
      chg(1)=qleft/xl**2
      chg(2)=qright/xl**2

c_______ WRITE PARAMETERS ___________________________________________
      write(*,*)
      write(*,*) ' SIMULATION PARAMETERS IN SI UNITS'
      write(*,*)
      write(*,'('' Reduced temperature              = '',f10.5)')
     :      tstar
      write(*,*)
      write(*,'('' L [A]                            = '',f10.5)')
     :      xl
      write(*,'('' H [A]                            = '',f10.5)')
     :      zl
      write(*,*)
      
      write(*,'('' Average diameter [A]             = '',f10.5)')
     :      daver
      write(*,'('' L [d]                            = '',f10.5)')
     :      xl/daver
      write(*,'('' H [d]                            = '',f10.5)')
     :      zl/daver
      write(*,*)
      
      write(*,'(12x,''   N   '',''  q [e]  '',''  d [A]  '')')
      write(*,'(''  cation 1  '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(1),qp1si,dp1
      if (n(2).gt.0) then
      write(*,'(''  cation 2  '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(2),qp2si,dp2
      endif
      write(*,'(''  anion     '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(3),qnsi,dn
      write(*,*)
      
      write(*,'('' Surface charge left [C/m2]       = '',f10.5)')
     :      sleft*sigsi
      write(*,'('' Surface charge right [C/m2]      = '',f10.5)')
     :      sigsi-sleft*sigsi
      write(*,*)
      
      write(*,'('' Surface charge left [e]          = '',f10.5)')
     :      sleft*dq
      write(*,'('' Surface charge right [e]         = '',f10.5)')
     :      dq-sleft*dq
      write(*,*)
      
      write(*,'('' Surface charge left [d]          = '',f10.5)')
     :      (sleft*dq)/(xl/daver)**2
      write(*,'('' Surface charge right [d]         = '',f10.5)')
     :      (dq-sleft*dq)/(xl/daver)**2
      write(*,*)
      
      write(*,'('' Surface charge left [reduced]    = '',f10.5)')
     :      chg(1)  
      write(*,'('' Surface charge right [reduced]   = '',f10.5)')
     :      chg(2)        
      write(*,*)
      
      chg11 = (sleft*dq)/xl**2
      chg22 = (dq-sleft*dq)/xl**2
      
      if (readit) then
 
c_______ READIN CONFIGURATION _______________________________________
      if (contin) then
       open (3,file='c')
      else
       open (3,file='c0')
      endif
      read(3,*)n(1),n(2),n(3)   
      read(3,*)d(1,1),d(2,2),d(3,3)    
      read(3,*)
      do ispec=1,3
      ni=n(ispec)
      if (ni.gt.0) then
      do 29 i=1,ni
29    read(3,*)rx(ispec,i),ry(ispec,i),rz(ispec,i)
      read(3,*)
      endif
      enddo
      close(3)

      else

c_______ MAKE CONFIGURATION _________________________________________
      call grand

      endif

c_______ DEFINE GRID FOR DISTRIBUTION FUNCTIONS _____________________ 
      
c      if (dz.gt.(dabs(dn-dp1)/2.0)) then
c            write(*,*)"dz is too large"
c            write(*,*)"it should be dz<(dn-dp1)/2"
c            stop
c      endif      
c      nzg1=int(dabs(dn-dp1)/dz/2.0)
      dzg=dz
      nzg1=nint((dp2/2.0-dp1/2.0)/dzg)
      write(*,'('' nzg1          = '',i8)')nzg1  
      nzg2=nint((dn/2.0-dp2/2.0)/dzg)
      write(*,'('' nzg2          = '',i8)')nzg2  
      nzg(1)=int((zl-dp1/2.0)/dzg)+1
      nzg(2)=int((zl-dp2/2.0)/dzg)+1
      nzg(3)=int((zl-dn/2.0)/dzg)+1

      write(*,'(''            dz = '',f8.5)')dzg
      write(*,'('' nzg (cation1) = '',i8)')nzg(1)
      if (n(2).gt.0) then
      write(*,'('' nzg (cation2) = '',i8)')nzg(2)
      endif
      write(*,'('' nzg (anion)   = '',i8)')nzg(3)
      if (nzg(1).gt.20000) stop 'nzg(1) is too big, should be < 20000'
      if (nzg(2).gt.20000) stop 'nzg(2) is too big, should be < 20000'
      if (nzg(3).gt.20000) stop 'nzg(3) is too big, should be < 20000'

      do ispec=1,3
      do 21 i=1,nzg(ispec)
21     gz(ispec,i) = 0.0
      enddo
       athist = 0.0
       
      return
      end

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

      subroutine move(ispec)
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      logical ovrlap

c      write(*,*)'move p1'
      i = int(dble(n(ispec))*ranff(mag)) + 1

      rxiold = rx(ispec,i)
      ryiold = ry(ispec,i)
      rziold = rz(ispec,i)
      atm(ispec)    = atm(ispec)+1

      rxinew = rxiold + ( 2.0 * ranff( mag ) - 1.0 ) * drmax(ispec)
      ryinew = ryiold + ( 2.0 * ranff( mag ) - 1.0 ) * drmax(ispec)
      rzinew = rziold + ( 2.0 * ranff( mag ) - 1.0 ) * drmax(ispec)
      dhalf=d(ispec,ispec)/2.0
      if (rzinew.lt.dhalf.or.rzinew.gt.(zl)) then
            au=au+u
            return
      endif      
      rxinew = rxinew - xl*dnint(rxinew/xl)
      ryinew = ryinew - xl*dnint(ryinew/xl)

      call ergi(ispec,i,rxiold,ryiold,rziold,uold,ucsho)
      if (ovrlap) stop 'Overlap in old move'
    
      call ergi(ispec,i,rxinew,ryinew,rzinew,unew,ucshn)
      if (ovrlap) then
       au = au + u
       return
      endif
      deltu  = unew - uold

      if ( deltu .lt. 75.0 ) then
         if ( deltu .le. 0.0 ) then
            u      = u + deltu 
            acm(ispec) = acm(ispec) + 1.0
            rx(ispec,i) = rxinew
            ry(ispec,i) = ryinew
            rz(ispec,i) = rzinew
         elseif ( dexp ( - deltu ) .gt. ranff( mag ) ) then
            u      = u + deltu 
            acm(ispec) = acm(ispec) + 1.0
            rx(ispec,i) = rxinew
            ry(ispec,i) = ryinew
            rz(ispec,i) = rzinew
         endif
        endif

       au =au + u
        
        return
        end

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

      subroutine jump(ispec)
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      logical ovrlap, move

c      write(*,*)'jump p1'
      i = int(dble(n(ispec))*ranff(mag)) + 1

      rxiold = rx(ispec,i)
      ryiold = ry(ispec,i)
      rziold = rz(ispec,i)
      atf(ispec) = atf(ispec)+1

      rxinew = ranff(mag)*xl-xl/2.0      
      ryinew = ranff(mag)*xl-xl/2.0 
      dhalf=d(ispec,ispec)/2.0
      rzinew = dhalf+ranff(mag)*(zl-dhalf)
      
      call ergi(ispec,i,rxiold,ryiold,rziold,uold,ucsho)
      if (ovrlap) stop 'Overlap in old jump'
    
      call ergi(ispec,i,rxinew,ryinew,rzinew,unew,ucshn)
      if (ovrlap) then
       au = au + u
       return
      endif

      deltu  = unew - uold

      if ( deltu .lt. 75.0 ) then
         if ( deltu .le. 0.0 ) then
            u      = u + deltu 
            acf(ispec) = acf(ispec) + 1.0
            rx(ispec,i) = rxinew
            ry(ispec,i) = ryinew
            rz(ispec,i) = rzinew
         elseif ( dexp ( - deltu ) .gt. ranff( mag ) ) then
            u      = u + deltu 
            acf(ispec) = acf(ispec) + 1.0
            rx(ispec,i) = rxinew
            ry(ispec,i) = ryinew
            rz(ispec,i) = rzinew
         endif
        endif

       au =au + u
        
        return
        end

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

      subroutine ergi(ipart,i,rxi,ryi,rzi,ui,ucshi)
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      logical ovrlap
      ovrlap=.false.
     
      ui=0.0
      ucshi=0.0
      qi=q(ipart)
      
      do 2 jpart=1,3
            ddi=d2(ipart,jpart)
            di=d(ipart,jpart)
            qj=q(jpart)
            nj=n(jpart)
            if (nj.eq.0) goto 2
            do 1 j=1,nj
                  if (ipart.eq.jpart.and.i.eq.j) goto 1
                  rxij=rxi-rx(jpart,j)
                  ryij=ryi-ry(jpart,j)
                  rzij=rzi-rz(jpart,j)
                  rxij=rxij-xl*dnint(rxij/xl)
                  ryij=ryij-xl*dnint(ryij/xl)
                  rijsq=rxij*rxij+ryij*ryij+rzij*rzij
                  if (rijsq.lt.ddi) then
                        ovrlap = .true.
                        return
                  endif 
                  rij=dsqrt(rijsq)
                  ui=ui+qi*qj/rij
                  zi=dabs(rzij)
                  wl=xl
                  ucshi = ucshi + qi*(qj/wl**2)*(-twopi*zi-wpot(zi,wl))
1     continue
2     continue      
      
      ui = ui + qi*chg(1)*(-twopi*rzi) + qi*chg(2)*(-twopi*(zl-rzi))
c      ui = ui + qi*chg(1)*wpot(rzi,xl) + qi*chg(2)*wpot(zl-rzi,xl)
      
      ui = ui + ucshi
      
      return
      end

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

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

c    ###################################################################
c    ###################################################################
      
      subroutine zero
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      common/logi/   readit,contin,yesmem
      common/sub/    difpot(10000),ksub
      logical readit, ovrlap, yesmem, contin
      do ispec=1,3
            acm(ispec)=0.0
            atm(ispec)=0.0
            acf(ispec)=0.0
            atf(ispec)=0.0
            if (.not.contin) then
            do 1 j=1,nzg(ispec)
1           gz(ispec,j)=0.0
            endif
      enddo
      au=0.0
      u=0.0
      athist=0.0
      ksub=0
      
      return
      end

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

      subroutine hist
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      
c      write(*,*)' begin hist'
      athist=athist+1.0
      
      do ispec=1,3
      
      ni=n(ispec)
      dhalf=d(ispec,ispec)/2.0
      
      do 1 i=1,ni
       rzi=rz(ispec,i)
       k0i=int((rzi-dhalf)/dzg)+1
       gz(ispec,k0i)=gz(ispec,k0i)+1.0
1     continue 
      
      enddo
     
      return
      end

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    #################################################################

      subroutine grand
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      logical  ovrlap

      write(*,*)' Making initial configuration'

c ************ cations **************************************** 

      do 1 ispec=1,3

      np1trg=n(ispec)
      np1=0
      ip1=0

30    ip1=ip1+1
20    rxinew=xl*ranff(mag)-xl/2.0
      ryinew=xl*ranff(mag)-xl/2.0
      dhalf=d(ispec,ispec)/2.0
      rzinew=dhalf+(zl-dhalf)*ranff(mag)
      call overlp(ispec,ip1,rxinew,ryinew,rzinew)
 
      if (ovrlap) then
       goto 20
      else
       rx(ispec,ip1)=rxinew
       ry(ispec,ip1)=ryinew
       rz(ispec,ip1)=rzinew
       np1=ip1
      endif 

      if (np1.lt.np1trg) goto 30

1     continue

      write(*,*)' Initial configuration made'

      open(2,file='cstart')  
      write(2,*)n(1),n(2),n(3)
      write(2,*)d(1,1),d(2,2),d(3,3)
      write(2,*)
      do 2 ispec=1,3
      ni=n(ispec)
      if (ni.gt.0) then
      do 18 i=1,n(ispec)
18    write(2,'(3(2x,f10.6))')rx(ispec,i),ry(ispec,i),rz(ispec,i)
      write(2,*)
      endif
2     continue      
      close(2)
     
      return
      end
       
c    ###################################################################
c    ###################################################################

      subroutine overlp(ispec,i,rxi,ryi,rzi)
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/ovr/    ovrlap
      logical ovrlap
      ovrlap=.false.

      do 2 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 1 j=1,nj
                  if (ispec.eq.jspec.and.i.eq.j) goto3
                  rxa=rx(jspec,j)
                  rya=ry(jspec,j)
                  rza=rz(jspec,j)
                  rxia=rxi-rxa
                  ryia=ryi-rya
                  rzia=rzi-rza
                  rxia=rxia-xl*dnint(rxia/xl)
                  ryia=ryia-xl*dnint(ryia/xl)
                  riasq=rxia*rxia+ryia*ryia+rzia*rzia
                  if (riasq.lt.dda) then
                        ovrlap=.true.
                        return
                  endif
3           continue
1           continue
            endif
2     continue            

      return
      end

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

      subroutine wrcf
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/wrout/  qp1si,qp2si,qnsi,chg11,chg22,tstar,eps,zlr(10),
     :               sigsi,sleft,ctarg
c      write(*,*)' begin wrcf'
c      write(*,*)' u =',u
      open(2,file='c')  
      write(2,*)n(1),n(2),n(3)
      write(2,*)z(1),z(2),z(3)
      write(2,*)d(1,1),d(2,2),d(3,3)
      write(2,*)nzg(1),nzg(2),nzg(3) 
      write(2,*)nzg1,nzg2 
      write(2,*)chg(1)              
      write(2,*)dzg,tsi,eps                
      write(2,*)zl,xl                 
      write(2,*)
      do 2 ispec=1,3
      ni=n(ispec)
      if (ni.gt.0) then
      do 18 i=1,n(ispec)
18    write(2,'(3(2x,f10.6))')rx(ispec,i),ry(ispec,i),rz(ispec,i)
      write(2,*)
      endif
2     continue      
      close(2)
      
      return
      end

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

      subroutine saves(iter,step)
      implicit doubleprecision (a-h,o-z)
      common/pos/    rx(3,1000),ry(3,1000),rz(3,1000)
      common/mol/    xl,zl,d2(3,3),d(4,4),q(3),chg(2),z(3),n(3),ntot
      common/const/  pi,twopi,sq2,sq3,tosi,tsi,mag
      common/erg/    u,au
      common/gz/     gz(3,20000),dzg,athist,nzg(3),nzg1,nzg2
      common/run/    drmax(3)
      common/sim/    nstep,naver,isave,iprint,ihist,iratio
      common/adj/    acm(3),atm(3),acf(3),atf(3)
      common/integ/  r(3,20000),c(3,20000),cbulk(3),cbr(10),daver 
      common/wrout/  qp1si,qp2si,qnsi,chg11,chg22,tstar,eps,zlr(10),
     :               sigsi,sleft,ctarg
      common/sub/    difpot(10000),ksub
      integer        step
      dimension densp(3,20000),densn(3,20000),
     :          charge(20000),rhoc(20000),rhocx(20000),poti(20000),
     :          total(20000)
      dimension rnew(2001),dpnew1(2001),dpnew2(2001),dnnew(2001)
      dimension chg1(2001),chg2(2001),chg3(2001)
      dimension x(20000),f(20000)
            
      rleft  = zl/3.0
      rright = 5.0*zl/6.0
      nright = int(rright/dzg)+1

      open(22,file='l')
       write(22,*)athist
      close(22)
      

      dnorm=(xl**2*dzg)*(athist)

      open(22,file='gp1')
      dhalf=d(1,1)/2.0
      do i=1,nzg(1)
            rzi=dhalf+(i-1)*dzg+dzg/2.0
            r(1,i)=rzi
            gzpp0=tosi*gz(1,i)/dnorm
            c(1,i)=gzpp0
            write(22,*)rzi,gzpp0
      enddo   
      close(22)

      if (n(2).gt.0) then
      open(22,file='gp2')
      dhalf=d(2,2)/2.0
      do i=1,nzg(2)
            rzi=dhalf+(i-1)*dzg+dzg/2.0
            r(2,i)=rzi
            gzpp0=tosi*gz(2,i)/dnorm
            c(2,i)=gzpp0
            write(22,*)rzi,gzpp0
      enddo   
      close(22)
      endif

      open(22,file='gn')
      dhalf=d(3,3)/2.0
      do i=1,nzg(3)
            rzi=dhalf+(i-1)*dzg+dzg/2.0
            r(3,i)=rzi
            gzpp0=tosi*gz(3,i)/dnorm
            c(3,i)=gzpp0
            write(22,*)rzi,gzpp0
      enddo   
      close(22)

      write(*,'(72(''_''))')
      write(*,*)step,'   Save '
      write(*,*)acm(1)/atm(1),acm(3)/atm(3),acf(1)/atf(1),acf(3)/atf(3)
c      write(*,*)' end saves'
      
      c1bulk=0.0
      a1=0.0
      do i=1,nzg(1)
            rzi=r(1,i)
            if (rzi.gt.rleft.and.rzi.lt.rright) then
                  c1bulk=c1bulk+c(1,i)
                  a1=a1+1.0
            endif      
      enddo   
      c1bulk=c1bulk/a1
      
      if (n(2).gt.0) then
      c2bulk=0.0
      a2=0.0
      do i=1,nzg(2)
            rzi=r(2,i)
            if (rzi.gt.rleft.and.rzi.lt.rright) then
                  c2bulk=c2bulk+c(2,i)
                  a2=a2+1.0
            endif      
      enddo   
      c2bulk=c2bulk/a2
      endif
     
      c3bulk=0.0
      a3=0.0
      do i=1,nzg(3)
            rzi=r(3,i)
            if (rzi.gt.rleft.and.rzi.lt.rright) then
                  c3bulk=c3bulk+c(3,i)
                  a3=a3+1.0
            endif      
      enddo   
      c3bulk=c3bulk/a3
      
      cbulk(1)=c1bulk
      cbulk(2)=c2bulk
      cbulk(3)=c3bulk
      
      write(*,*)' Bulk concentration of cation(1) = ',c1bulk
      if (n(2).gt.0) then
      write(*,*)' Bulk concentration of cation(2) = ',c2bulk
      endif
      write(*,*)' Bulk concentration of cation(3) = ',c3bulk

      open(2,file='r')
      
      write(2,*)" Different ion diameters " 
c      write(2,*)" No second cation (n2=0)" 
      write(2,*)" Charge only on left wall (sleft=1)"
      write(2,*)
      
      dbp1=cbulk(1)
      dbp2=cbulk(2)
      dbn =cbulk(3)
     
      write(2,'(12x,''   N   '',''  q [e]  '',''  d [A]  '')')
      write(2,'(''  cation 1  '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(1),qp1si,d(1,1)
c      if (n(2).gt.0) then
      write(2,'(''  cation 2  '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(2),qp2si,d(2,2)
c      endif
      write(2,'(''  anion     '',i7,2x,f5.2,2x,2x,f8.4)')
     :      n(3),qnsi,d(3,3)
      write(2,*)
      
      write(2,'('' Reduced quantities are calculated by the average'',
     :          '' diameter '',f5.3,'' A'')')daver
      write(2,*)

      write(2,'('' Surface charge [C/m2]  = '',f10.5)')
     :      sleft*sigsi
      write(2,'('' Surface charge reduced = '',f10.5)')
     :      chg11*daver**2
      write(2,'('' Reduced temperature    = '',f10.5)')
     :      tstar

      write(2,*)
      write(2,'('' Geometry'')')
      write(2,'('' L [A] = '',f8.3)')xl
      write(2,'('' Targeted c(anion) = '',f8.5,'' M'')')ctarg
      write(2,*)
      
      cbr(iter) = cbulk(3)
      write(2,'('' History '')')
            write(2,'('' Iter '','' H [A]  '', ''  c(anion) '')')
      do it=1,iter
            write(2,'(2x,i2,2x,f8.3,1x,f8.5)')it,zlr(it),cbr(it)
      enddo
      write(2,*)
      write(2,*)' c(cation1) = ',cbulk(1)
      write(2,*)' c(cation2) = ',cbulk(2)
      write(2,*)' c(anion) =  ',cbulk(3)
      write(2,*)
c      write(*,'(40(''_''))')
      
c ******************************************************************
c The calculations from here are for 
c     equal ion diameters 
c     no second cation (n2=0) 
c     charge only on left wall (sleft=1)
c ******************************************************************
      
      write(2,'(''Number of MC cycles = '',i6)')step
      write(2,*)
      
c ***** Contact values on the left *********************************

      dsi2=d(1,1)/2.0
      call fitl(1,nzg(1),3,r,c,al,bl,cl)
      gpl1=dsi2**2*al+dsi2*bl+cl
      write(2,'('' Contact left - cations1  = '',f10.5)') gpl1  
      dsi2=d(2,2)/2.0
      call fitl(2,nzg(2),3,r,c,al,bl,cl)
      gpl2=dsi2**2*al+dsi2*bl+cl
      write(2,'('' Contact left - cations2  = '',f10.5)') gpl2  
      dsi2=d(3,3)/2.0
      call fitl(3,nzg(3),3,r,c,al,bl,cl)
      gnl=dsi2**2*al+dsi2*bl+cl
      write(2,'('' Contact left - anions   = '',f10.5)') gnl  

c      call fitr(1,nzg,3,r,c,al,bl,cl)
c      gpr=dsi2**2*al+dsi2*bl+cl
c      write(2,'('' Contact right - cations = '',f10.5)') gpr  
c      call fitr(3,nzg,3,r,c,al,bl,cl)
c      gnr=dsi2**2*al+dsi2*bl+cl
c      write(2,'('' Contact right - anions  = '',f10.5)') gnr  
      write(2,*)                                              

c ***** New curve **************************************************
      
c     if (d(1,1).gt.d(3,3)) then
c           write(*,*)" d(anion) should be larger than d(cation)"
c           stop
c     endif

      open(22,file='rhop1.n')
      open(24,file='rhop2.n')
      open(23,file='rhon.n')
      open(12,file='charge')

      dz2=dzg/2.0
      rnew(1)  = r(1,1) - dz2
      dpnew1(1) = gpl1
      chg1(1)   = gpl1
      ii=1
       write(22,*)rnew(1),0.0       
       write(22,*)rnew(1),dpnew1(1)
       write(24,*)rnew(1),0.0
       write(23,*)rnew(1),0.0
       write(12,*)rnew(1),chg1(1)
      
      do i=1,nzg1
      
       ii=ii+1
       rnew(ii) = r(1,i)+dz2
       x(1) = r(1,i)
       x(2) = r(1,i+1)
       x(3) = r(1,i+2)
       
       f(1) = c(1,i)
       f(2) = c(1,i+1)
       f(3) = c(1,i+2)
       call fit(3,x,f,a,b,c0)
       dpnew1(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       chg1(ii)=z(1)*dpnew1(ii)
       write(12,*)rnew(ii),chg1(ii)
       write(22,*)rnew(ii),dpnew1(ii)

       if (i.eq.nzg1) then
            dpnew2(ii) = gpl2
            write(24,*)rnew(ii),0.0
            write(24,*)rnew(ii),dpnew2(ii)
            chg2(ii)=z(1)*dpnew1(ii)+z(2)*dpnew2(ii)
            write(12,*)rnew(ii),chg2(ii)
       endif     

      enddo   

      
      do i=1,nzg2
      
       ii=ii+1
       rnew(ii) = r(1,nzg1+i)+dz2
c     cation1       
       x(1) = r(1,nzg1+i)
       x(2) = r(1,nzg1+i+1)
       x(3) = r(1,nzg1+i+2)
       f(1) = c(1,nzg1+i)
       f(2) = c(1,nzg1+i+1)
       f(3) = c(1,nzg1+i+2)
       call fit(3,x,f,a,b,c0)
       dpnew1(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       write(22,*)rnew(ii),dpnew1(ii)
c     cation2       
       x(1) = r(2,i)
       x(2) = r(2,i+1)
       x(3) = r(2,i+2)
       f(1) = c(2,i)
       f(2) = c(2,i+1)
       f(3) = c(2,i+2)
       call fit(3,x,f,a,b,c0)
       dpnew2(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       write(24,*)rnew(ii),dpnew2(ii)
      
       chg2(ii)=z(1)*dpnew1(ii)+z(2)*dpnew2(ii)
       write(12,*)rnew(ii),chg2(ii)
      
       if (i.eq.nzg2) then
            dnnew(ii) = gnl
            write(23,*)rnew(ii),0.0
            write(23,*)rnew(ii),dnnew(ii)
            chg3(ii)=z(1)*dpnew1(ii)+z(2)*dpnew2(ii)+z(3)*dnnew(ii)
            write(12,*)rnew(ii),chg3(ii)
       endif     
      
      enddo  


      do i=1,nright
      
       ii=ii+1
       rnew(ii) = r(1,nzg1+nzg2+i)+dz2
c     cation1       
       x(1) = r(1,nzg1+nzg2+i)
       x(2) = r(1,nzg1+nzg2+i+1)
       x(3) = r(1,nzg1+nzg2+i+2)
       f(1) = c(1,nzg1+nzg2+i)
       f(2) = c(1,nzg1+nzg2+i+1)
       f(3) = c(1,nzg1+nzg2+i+2)
       call fit(3,x,f,a,b,c0)
       dpnew1(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       write(22,*)rnew(ii),dpnew1(ii)
c     cation2       
       x(1) = r(2,nzg2+i)
       x(2) = r(2,nzg2+i+1)
       x(3) = r(2,nzg2+i+2)
       f(1) = c(2,nzg2+i)
       f(2) = c(2,nzg2+i+1)
       f(3) = c(2,nzg2+i+2)
       call fit(3,x,f,a,b,c0)
       dpnew2(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       write(24,*)rnew(ii),dpnew2(ii)
c     anion         
       x(1) = r(3,i)
       x(2) = r(3,i+1)
       x(3) = r(3,i+2)
       f(1) = c(3,i)
       f(2) = c(3,i+1)
       f(3) = c(3,i+2)
       call fit(3,x,f,a,b,c0)
       dnnew(ii) = a*rnew(ii)**2 + b*rnew(ii) + c0
       write(23,*)rnew(ii),dnnew(ii)

       chg3(ii)=z(1)*dpnew1(ii)+z(2)*dpnew2(ii)+z(3)*dnnew(ii)
       write(12,*)rnew(ii),chg3(ii)

      enddo   
      
      close(22)
      close(24)
      close(23)
      close(12)
     
      return
      end
      
c    ###################################################################
c    ###################################################################

       subroutine simp2(n1,n2,r,f,ff)
       implicit doubleprecision (a-h,o-z)
       dimension r(20000),f(20000)
       m=n2-n1
       ba=r(n2)-r(n1)
       dx=ba/dble(m)
       if (m.eq.1) then
        ff=dx*(f(n1)+f(n2))/2.0
        return
       endif
       if (m.eq.2) then 
        ff=dx*(f(n1)/3.0+4*f(n1+1)/3.0+f(n2)/3.0)
        return
       endif
       if (m.eq.3) then
        ff=dx*(3*(f(n1)+f(n2))+9*(f(n1+1)+f(n2-1)))/8.0
        return
       endif
       if (m.eq.4) then
        ff=dx*(14*(f(n1)+f(n2))+64*(f(n1+1)+f(n2-1))+24*f(n1+2))/45.0
        return
       endif
       ff=3*(f(n1)+f(n2))/8.0+7*(f(n1+1)+f(n2-1))/6.0+
     :    23*(f(n1+2)+f(n2-2))/24.0 
       do 1 i=n1+3,n2-3
1      ff=ff+f(i)
       ff=dx*ff
       return
       end
                            
c    ###################################################################
c    ###################################################################

      subroutine fit(n,x,f,a,b,c)
      implicit doubleprecision (a-h,o-z)
      dimension x(20000),f(20000)
      x0=0.0
      x1=0.0
      x2=0.0
      x3=0.0
      x4=0.0
      fx0=0.0
      fx1=0.0
      fx2=0.0
      do 1 i=1,n
      x0=x0+1.0
      x1=x1+x(i)
      x2=x2+x(i)**2
      x3=x3+x(i)**3
      x4=x4+x(i)**4
      fx0=fx0+f(i)
      fx1=fx1+f(i)*x(i)
      fx2=fx2+f(i)*x(i)**2
1     continue
      y1=x1-x2*x2/x3
      y2=x0-x2*x1/x3
      y3=fx0-x2*fx1/x3
      y4=x3-x4*x2/x3
      y5=x2-x4*x1/x3
      y6=fx2-x4*fx1/x3
      z1=y2-y1*y5/y4
      z2=y3-y1*y6/y4
      c=z2/z1
      b=(y6-c*y5)/y4
      a=(fx0-c*x0-b*x1)/x2
      return
      end

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

      subroutine fitl(k,nzg,n,x,f,a,b,c)
      implicit doubleprecision (a-h,o-z)
      dimension x(3,20000),f(3,20000)
      x0=0.0
      x1=0.0
      x2=0.0
      x3=0.0
      x4=0.0
      fx0=0.0
      fx1=0.0
      fx2=0.0
      do 1 i=1,n
      x0=x0+1.0
      x1=x1+x(k,i)
      x2=x2+x(k,i)**2
      x3=x3+x(k,i)**3
      x4=x4+x(k,i)**4
      fx0=fx0+f(k,i)
      fx1=fx1+f(k,i)*x(k,i)
      fx2=fx2+f(k,i)*x(k,i)**2
1     continue
      y1=x1-x2*x2/x3
      y2=x0-x2*x1/x3
      y3=fx0-x2*fx1/x3
      y4=x3-x4*x2/x3
      y5=x2-x4*x1/x3
      y6=fx2-x4*fx1/x3
      z1=y2-y1*y5/y4
      z2=y3-y1*y6/y4
      c=z2/z1
      b=(y6-c*y5)/y4
      a=(fx0-c*x0-b*x1)/x2
      return
      end

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

      subroutine fitr(k,nzg,n,x,f,a,b,c)
      implicit doubleprecision (a-h,o-z)
      dimension x(3,20000),f(3,20000)
      x0=0.0
      x1=0.0
      x2=0.0
      x3=0.0
      x4=0.0
      fx0=0.0
      fx1=0.0
      fx2=0.0
      do 1 i=1,n
      x0=x0+1.0
      x1=x1+x(k,i)
      x2=x2+x(k,i)**2
      x3=x3+x(k,i)**3
      x4=x4+x(k,i)**4
      fx0=fx0+f(k,nzg-i+1)
      fx1=fx1+f(k,nzg-i+1)*x(k,i)
      fx2=fx2+f(k,nzg-i+1)*x(k,i)**2
1     continue
      y1=x1-x2*x2/x3
      y2=x0-x2*x1/x3
      y3=fx0-x2*fx1/x3
      y4=x3-x4*x2/x3
      y5=x2-x4*x1/x3
      y6=fx2-x4*fx1/x3
      z1=y2-y1*y5/y4
      z2=y3-y1*y6/y4
      c=z2/z1
      b=(y6-c*y5)/y4
      a=(fx0-c*x0-b*x1)/x2
      return
      end

