      program integ
      implicit doubleprecision (a-h,o-z)
      parameter(nionmx=500,ntotmx=1500,ntrymx=1500, 
     :          nspcmx=3,nzgmx=400)
      common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag
      common/conf/  rx(ntotmx),ry(ntotmx),rz(ntotmx),
     :              indspc(nspcmx,nionmx),rii(ntotmx,ntotmx), 
     :              n(nspcmx),ispcbk(ntotmx)
      common/spec/  d(nspcmx),dd(nspcmx,nspcmx),d2(nspcmx,nspcmx),
     :              dca(nspcmx),z(nspcmx),q(nspcmx),sigma(2),sigold(2),
     :              ratmov(nspcmx),sigsi(2),tstar,voltg, 
     :              eps(nspcmx),epsw,qstar,itry(nspcmx),
     :              nspec,nwh(ntrymx),ntry
      common/geom/  xl,zl,drmax,fion(nspcmx),fsalt(nspcmx)
      common/gc/    chemp(3),ctarg(3),vgc(3),riicrn(4,4),
     :              rxnew(4),rynew(4),rznew(4),r2new(4),
     :              chex(3),ctrgcl,ratgr(3),ratgrt,
     :              iout(4),ispcin(3,4),igcval(3),
     :              isalt(3),ication(3),nsalt,ntot,niter,iter
      common/profl/ athist,dzg(nspcmx),gz(nspcmx,nzgmx),
     :              pol(nspcmx,nzgmx),nzg(nspcmx),fgz(nspcmx,2)
      common/accu/  an(nspcmx),cbulk(nspcmx),asigma(2)
      common/adj/   amove(2,nspcmx),ajump(2,nspcmx),
     :              acreat(3,2),adest(3,2),achg(2)
      common/bufer/ riinew(nspcmx,ntotmx)
      common/simp/  ratchg,dchgmx,nstep,naver,isave,runit
      common/erg/   uc,auc,ucav,up,aup,upav,adip(nspcmx),adipsq(nspcmx)
      common/itern/ dipxn(ntotmx),dipyn(ntotmx),dipzn(ntotmx),
     :              dipx(ntotmx),dipy(ntotmx),dipz(ntotmx),
     :              dipx0n(ntotmx),dipy0n(ntotmx),dipz0n(ntotmx),
     :              dipx0(ntotmx),dipy0(ntotmx),dipz0(ntotmx),
     :              eqxn(ntotmx,ntotmx),eqyn(ntotmx,ntotmx),
     :              eqzn(ntotmx,ntotmx),
     :              eqx(ntotmx,ntotmx),eqy(ntotmx,ntotmx),
     :              eqz(ntotmx,ntotmx),
     :              epxn(ntotmx,ntotmx),epyn(ntotmx,ntotmx),
     :              epzn(ntotmx,ntotmx),
     :              epx(ntotmx,ntotmx),epy(ntotmx,ntotmx),
     :              epz(ntotmx,ntotmx),
     :              alpha(nspcmx),ialpha
      logical       runit
      logical       readit, ovrlap, targc, difdca
      character*2   fion
      character*4   fsalt
      character*9   fgz
      character*9   fngz,fnpz
      dimension     cz(nspcmx,nzgmx),dinteg(nspcmx),gamma1(nspcmx),
     :              gamma2(nspcmx),pz(nspcmx,nzgmx),pinteg(nspcmx),
     :              einteg(nspcmx)
      
      mag=51328632
      pi=4.0*datan(1.d0) 
      twopi=2*pi
      sq2=dsqrt(2.d0)
      sq3=dsqrt(3.d0)
      tsi=300.0
      eps0=8.8542d-12
      boltz=1.3806d-23
      echg=1.6021917d-19
      beta=1./(boltz*tsi)
      avog=6.02214d23
      tomv=boltz*tsi/echg
      dsi=1.d-10
      tosi=1.0/(1.d-10)**3/avog/1000.

      open(2,file="in.dat")
     
c---read in first AND geometry---
      read(2,*)niter
      write(*,*)1
      read(2,*)targc,ntrg,ctrg
      read(2,*)xl,zl
      read(2,*)sigma(1),sigma(2)
      read(2,*)voltg
      read(2,*)
      read(2,*)

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

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

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

      close(2)
     
      open(2,file="c.dat")  
      read(2,*) nspec,istep                             
      read(2,*)
      do ispec=1,nspec
            read(2,*)n(ispec),d(ispec),dca(ispec)
            read(2,*)nzg(ispec),dzg(ispec)
            read(2,*)
            do i=1,n(ispec)
                  ii=indspc(ispec,i)
                  read(2,*)rx(ii),ry(ii),rz(ii)
            enddo      
      read(2,*)
      enddo
      close(2)

      qstar=echg*dsqrt(beta/(4*pi*eps0*eps(1)*dsi))
      write(*,*)" qstar = ",qstar,echg,beta,pi,eps0,eps(1),dsi
      nzg2=nzg(1)/2

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

        do ispec=1,nspec
          fngz  = fgz(ispec,1)
          fnpz  = fgz(ispec,2)
          ri=d(ispec)/2.0
          open(12,file=fngz)
          open(22,file=fnpz)
           do i=1,nzg(ispec)
             read(12,*)rgz,cz(ispec,i),dzg(ispec)
             read(22,*)rgz,pz(ispec,i)
           enddo   
          close(12)
          close(22)
        enddo
      
      
      open(22,file="pot-left.dat")

      do j=1,nzg2-2
         zgj=dca(1)+(j-1)*dzg(1)
         pot1=0.0
         potp1=0.0
         do ispec=1,nspec
           dinteg(ispec)=0.0
           einteg(ispec)=0.0
           pinteg(ispec)=0.0
           gamma1(ispec)=0.0
           do i=j,nzg2
             zg1=dca(ispec)+(i-1)*dzg(ispec)
             zg2=zg1+dzg(ispec)
             dinteg(ispec)=dinteg(ispec)
     :           +(cz(ispec,i)/tosi)*(zg2**2-zg1**2)/2.0
             einteg(ispec)=einteg(ispec)
     :           +(cz(ispec,i)/tosi)* dzg(ispec)
             pinteg(ispec)=pinteg(ispec)
     :           +pz(ispec,i)* dzg(ispec)
             gamma1(ispec)=gamma1(ispec)
     :                     +(cz(ispec,i)-cbulk(ispec))*dzg(ispec)
           enddo 
           pot1=pot1+z(ispec)*(zgj*einteg(ispec)-dinteg(ispec))
           potp1=potp1+pinteg(ispec)
           gamma1(ispec)=(gamma1(ispec)/tosi)*xl**2
         enddo
         pot1=pot1*2*twopi*qstar**2
         potp1=-potp1*2*twopi
         write(22,*)zgj,pot1+potp1,pot1,potp1
      enddo
      close(22)
      
c      stop

c         write(*,*)"  Ion potential at left wall /kT  = ",pot1
c         write(2,*)"  Ion potential at left wall /kT  = ",pot1
c         write(*,*)" Pol. potential at left wall /kT  = ",potp1
c         write(2,*)" Pol. potential at left wall /kT  = ",potp1
c         write(*,*)" Tot. potential at left wall /kT  = ",potp1+pot1
c         write(2,*)" Tot. potential at left wall /kT  = ",potp1+pot1
        
      open(23,file="pot-right.dat")
      do j = 1, nzg2-2
       zgj=dca(1)+(j-1)*dzg(1)
       pot2=0.0
       potp2=0.0
       do ispec=1,nspec
c        nzg2=nzg(ispec)/2
        dinteg(ispec)=0.0
        einteg(ispec)=0.0
        pinteg(ispec)=0.0
        gamma2(ispec)=0.0
        do i=j,nzg2
          zg1=dca(ispec)+(i-1)*dzg(ispec)
          zg2=zg1+dzg(ispec)
          ii=nzg(ispec)-i+1
          dinteg(ispec)=dinteg(ispec)
     :        +(cz(ispec,ii)/tosi)*(zg2**2-zg1**2)/2.0
          einteg(ispec)=einteg(ispec)
     :     +(cz(ispec,ii)/tosi)*dzg(ispec)
          pinteg(ispec)=pinteg(ispec)
     :        -pz(ispec,ii)*dzg(ispec)
          gamma2(ispec)=gamma2(ispec)
     :                  +(cz(ispec,ii)-cbulk(ispec))*dzg(ispec)
       enddo 
        pot2=pot2+z(ispec)*(zgj*einteg(ispec)-dinteg(ispec))
        potp2=potp2+pinteg(ispec)
        gamma2(ispec)=(gamma2(ispec)/tosi)*xl**2
      enddo
      pot2=-pot2*2*twopi*qstar**2
      potp2=-potp2*2*twopi
      write(23,*)zgj,pot2+potp2,pot2,potp2
c      write(*,*)"  Ion potential at right wall /kT  = ",pot2
c      write(2,*)"  Ion potential at right wall /kT  = ",pot2
c      write(*,*)" Pol. potential at right wall /kT  = ",potp2
c      write(2,*)" Pol. potential at right wall /kT  = ",potp2
c      write(*,*)" Tot. potential at right wall /kT  = ",potp2+pot2
c      write(2,*)" Tot. potential at right wall /kT  = ",potp2+pot2
      end do  
      close(23)
c     stop
      end
