program ionch ! the original double layer code implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz write(*,*)" Start MC Double layer" open(2,file="in.dat") read(2,*)niter close(2) do iter=1,niter call readin istep0=0 call wrcf(istep0) call zero write(*,*)" Start equilibration" do istep=1,naver do ip=1,ntotmx call tryit enddo if (mod(istep,isave).eq.0) then call wrcf(istep) write(*,*)istep," I am equilibrating" endif enddo call zero write(*,*)" Start simulation" do istep=1,nstep do ip=1,ntotmx call tryit call hist enddo if (mod(istep,isave).eq.0) call saves(istep) if (mod(istep,isave).eq.0) call wrcf(istep) enddo if (iter.eq.niter) then write(*,*)" End of simulation" else write(*,*)" End of iteration ",iter write(*,'(72("="))') endif enddo stop end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine readin implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical readit, ovrlap, targc, difdca dimension xn(nspcmx),xq(nspcmx) c_______ CONSTANTS __________________________________________________ mag=51328632 pi=4.0*datan(1.d0) twopi=2*pi sq2=dsqrt(2.d0) sq3=dsqrt(3.d0) c_______ READIN PARAMETERS __________________________________________ c_______ CHECK FOR SIZES OF ARRAYS __________________________________ open(2,file="in.dat") c---read in first AND geometry--- read(2,*)niter write(*,*)1 read(2,*)targc,ntrg,ctrg read(2,*)xl,zl read(2,*)sigma(1),sigma(2) read(2,*)voltg read(2,*) read(2,*) c---simulation parameters-- read(2,*)drmax read(2,*)dchgmx read(2,*)dzgtrg read(2,*)readit read(2,*)dsi read(2,*)tstar read(2,*)runit read(2,*)difdca read(2,*)ratchg read(2,*)nstep0 read(2,*)naver read(2,*)isave read(2,*) read(2,*) c---species--- read(2,*)nspec read(2,*) write(*,'(" nspec = ",i6," max = ",i6)')nspec,nspcmx if (nspec.gt.nspcmx) stop "nspec too large" do ispec=1,nspec read(2,*)eps(ispec),ratmov(ispec), : z(ispec),d(ispec),dca(ispec),fion(ispec) if (.not.difdca) dca(ispec)=d(ispec)/2.0 enddo 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) c specify constants tsi=300.0 eps0=8.8542d-12 boltz=1.3806d-23 echg=1.6021917d-19 beta=1./(boltz*tsi) avog=6.02214d23 tomv=boltz*tsi/echg write(*,*)" tomv = ",tomv," Volts" if (runit) then tosi=1.0 qstar=dsqrt(1.0/tstar) write(*,*)"-----------------------------------------" write(*,*)" I work in reduced units " write(*,*)" T* = ",tstar write(*,*)" q* = ",qstar write(*,*)" dsi = ",dsi write(*,*)"-----------------------------------------" else dsi=1.d-10 tosi=1.0/(1.d-10)**3/avog/1000. epsi=eps(1) qstar=echg*dsqrt(beta/(4*pi*eps0*epsi*dsi)) tstar=1.0/qstar**2 write(*,*)"-----------------------------------------" write(*,*)" I work in real units " write(*,*)" T* = ",tstar write(*,*)" q* = ",qstar write(*,*)" dsi = ",dsi write(*,*)"-----------------------------------------" endif voltg=voltg/qstar c calc cell dimensions if targc=.true. if (targc) then rtrg=ctrg/tosi vbulkt=dble(ntrg)/rtrg vhalf=vbulkt/2.0 xl=vhalf**(1./3.) zl=xl endif write(*,*)" zl = ",zl write(*,*)" starting xl = ",xl write(*,*)"-----------------------------------------" c calc cation mole fractions do ispec=2,nspec xn(ispec)=ctarg(ispec-1)/ctot write(*,*)" cation mole fraction : ",ispec,xn(ispec) enddo c redefine ion numbers n(2)=ntrg write(*,*)" N(",2,") = ",n(2) if (nspec.ge.3) then ncation=n(2) do ispec=3,nspec n(ispec)=int(n(2)*xn(ispec)/xn(2)) write(*,*)" N(",ispec,") = ",n(ispec) ncation=ncation+n(ispec) enddo endif c sigma in e unit if (.not.runit) then sigsi(1)=sigma(1) sigsi(2)=sigma(2) sigma(1)=(sigma(1)/echg)*dsi**2 sigma(2)=(sigma(2)/echg)*dsi**2 endif c recalculate xl qwall=(sigma(1)+sigma(2))*xl**2 write(*,*)" net charge on left wall /e = ",sigma(1)*xl**2 write(*,*)" net charge on right wall /e = ",sigma(2)*xl**2 write(*,*)" net charge on walls /e (starting) = ",qwall nwall=int(qwall) if (dabs(qwall).gt.1.d-5) then xlsq=dble(nwall)/(sigma(1)+sigma(2)) xl=dsqrt(xlsq) write(*,*)" xl = ",xl qwalle=(sigma(1)+sigma(2))*xl**2 write(*,*)" net charge on walls /e = ",qwalle endif sigma(1)=sigma(1)*qstar sigma(2)=sigma(2)*qstar write(*,*)" Initial sigma = ",sigma(1) qwall=sigma(1)+sigma(2) qwall=qwall*xl**2 write(*,*)" net charge on walls = ",qwall c calc N(Cl) from condition of electroneutrality nchg=nwall do ispec=2,nspec nchg=nchg+nint(z(ispec)*n(ispec)) enddo n(1)=nchg write(*,*)" N(",1,") = ",n(1) write(*,*)"-----------------------------------------" ntry=0 nimax=0 ntot=0 do ispec=1,nspec itry(ispec)=n(ispec) ntry=ntry+itry(ispec) ntot=ntot+n(ispec) if (n(ispec).gt.nimax) nimax=n(ispec) enddo write(*,'(" nimax = ",i6," max = ",i6)')nimax,nionmx write(*,'(" ntot = ",i6," max = ",i6)')ntot,ntotmx write(*,'(" ntry = ",i6," max = ",i6)')ntry,ntrymx if (nimax.gt.nionmx) stop "nimax too large" if (ntot.gt.ntotmx) stop "ntot too large" if (ntry.gt.ntrymx) stop "ntry too large" write(*,*)"-----------------------------------------" c constract output filenames open(2,file="fnames.dat") do ispec=1,nspec write(2,'("gz-",a2,".dat")')fion(ispec) enddo close(2) open(2,file="fnames.dat") do ispec=1,nspec read(2,*)fgz(ispec) enddo close(2) if (iter.eq.niter) then nstep=nstep0 else nstep=nstep0/10 endif do igc=1,nsalt ispec=isalt(igc) igcval(igc)=nint(z(ispec)) iv=igcval(igc) ispcin(igc,1)=isalt(igc) do icount=2,iv+1 ispcin(igc,icount)=1 enddo vgc(igc) = 1.0 do i=1,iv+1 ispec = ispcin(igc,i) ri = d(ispec)/2.0 vgc(igc)=vgc(igc)*xl**2*2*(zl-dca(ispec)) enddo write(*,*)" vgc = ",igc,vgc(igc) enddo ctrgcl=0.0 do igc=1,nsalt iv=igcval(igc) ctrgcl=ctrgcl+iv*ctarg(igc) enddo write(*,*)" Cl conc = ",ctrgcl write(*,*)"-----------------------------------------" if (iter.eq.1) then write(*,'(" iter = ",i2)')iter do igc=1,nsalt iv=igcval(igc) rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv chemp(igc)=dlog(rtargi)+chex(igc) write(*,*)igc,chex(igc),chemp(igc),iv enddo else write(*,*) write(*,'(" iter = ",i2)')iter do igc=1,nsalt iv=igcval(igc) rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv ispec=isalt(igc) rbulki=cbulk(ispec)/tosi rblkcl=cbulk(1)/tosi chexo=chex(igc) chex(igc)=chemp(igc)-dlog(rbulki*rblkcl**iv) chempo=chemp(igc) chemp(igc)=dlog(rtargi)+chex(igc) write(*,*) igc,ctarg(igc),cbulk(ispec) write(*,*) igc,chempo,chemp(igc) write(*,*) igc,chexo,chex(igc) enddo endif c_______ CALCULATE REDUCED CHARGES __________________________________ c_______ CHECK FOR ELECTRONEUTRALITY ________________________________ c_______ FILL UP DD, D2, AND NWH ARRAYS _____________________________ charge=qwalle ii=0 do i=1,nspec charge=charge+z(i)*n(i) q(i)=z(i)*qstar ri=d(i)/2.0 do j=1,nspec dd(i,j)=(d(i)+d(j))/2.0 d2(i,j)=dd(i,j)**2 enddo do j=1,itry(i) ii=ii+1 nwh(ii)=i enddo enddo if (charge.gt.0.001) stop "Cell is no neutral" if (readit) then c_______ READIN CONFIGURATION _______________________________________ c_______ FILL UP INDEX ARRAYS AND LOOKUP TABLES _____________________ open (3,file="c0") read(3,*)nspec read(3,*) do ispec=1,nspec read(3,*)n(ispec) read(3,*) read(3,*) do i=1,n(ispec) ii=(ispec-1)*nionmx+i read(3,*)rxi,ryi,rzi rx(ii)=rxi ry(ii)=ryi rz(ii)=rzi indspc(ispec,i)=ii ispcbk(ii)=i enddo read(3,*) enddo close(3) else c_______ MAKE CONFIGURATION _________________________________________ call initc endif call lookup c_______ Geometry for gz distributions ______________________________ do i=1,nspec nzg(i) = int(2*(zl-dca(i))/dzgtrg)+1 dzg(i) = 2*(zl-dca(i))/dble(nzg(i)) if (nzg(i).gt.nzgmx) then write(*,*) "nzg is too large" write(*,'(" nzg = ",i6," nzgmx = ",i6)')nzg(i),nzgmx stop endif write(*,'(40("_"))') enddo return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine tryit implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz ireg=1 vszam=ranff(mag) ratgr2=0.0 do igc=1,nsalt ratgr2=ratgr2+ratgr(igc) ratgr1=ratgr2-ratgr(igc) if (vszam.gt.ratgr1.and.vszam.lt.ratgr2) then if (ranff(mag).lt.0.5) then call create(igc) else call destr(igc) endif endif enddo if (vszam.gt.ratgr2.and.vszam.lt.ratgr2+ratchg) then call rechg else iwh = int(ranff(mag)*ntry)+1 ispec = nwh(iwh) if (n(ispec).eq.0) goto 1 i = int(dble(n(ispec))*ranff(mag)) + 1 call move(ispec,i) 1 continue endif return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine lookup implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz do ispec=1,nspec do i=1,n(ispec) ii=indspc(ispec,i) rxi=rx(ii) ryi=ry(ii) rzi=rz(ii) do jspec=1,nspec do j=1,n(jspec) jj=indspc(jspec,j) if (ii.ne.jj) then rxj=rx(jj) ryj=ry(jj) rzj=rz(jj) rxij=rxi-rxj ryij=ryi-ryj rzij=rzi-rzj rxij = rxij-dnint(rxij/xl)*xl ryij = ryij-dnint(ryij/xl)*xl rijsq= rxij*rxij+ryij*ryij+rzij*rzij if (rijsq.lt.d2(ispec,jspec)) stop : "lookup: Overlap in initial configuration" rii(ii,jj)=dsqrt(rijsq) else rii(ii,jj)=0.0 endif enddo enddo enddo enddo return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ function ranff ( seed ) implicit doubleprecision (a-h,o-z) integer a14,p31,seed,b15,b16,xhi,xalo,leftlo,fhi,kf15 data a14/16807/,b15/32768/,b16/65536/,p31/2147483647/ xhi=seed/b16 xalo=(seed-xhi*b16)*a14 leftlo =xalo/b16 fhi=xhi*a14+leftlo kf15=fhi/b15 seed=(((xalo-leftlo*b16 )-p31) +(fhi-kf15*b15)*b16) +kf15 if (seed .lt. 0) seed=seed+p31 xhi=seed/b16 xdr=(float(xhi)*65536.0d0)+float(seed-xhi*b16) ranff=xdr*4.6566128752457969241d-10 seed=xdr return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine hist implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz athist=athist+1.0 do ispec=1,nspec do i=1,n(ispec) ii=indspc(ispec,i) rzi=rz(ii) k=int((rzi+zl-dca(ispec))/dzg(ispec))+1 gz(ispec,k)=gz(ispec,k)+1.0 if (dabs(rzi).lt.zl/2.0) an(ispec)=an(ispec)+1.0 enddo enddo asigma(1)=asigma(1)+sigma(1) asigma(2)=asigma(2)+sigma(2) return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine initc implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. write(*,*)" Making initial configuration" do ispec=1,nspec do i=1,nionmx indspc(ispec,i)=0 enddo enddo 11 continue ii=0 do 10 ispec=1,nspec ntrg=n(ispec) ni=0 30 ni=ni+1 20 rxinew = (ranff(mag)-0.5)*xl ryinew = (ranff(mag)-0.5)*xl rzinew = (2*ranff(mag)-1)*(zl-dca(ispec)) ii=ii+1 if (ii.gt.1000000) goto 11 call overlp(ispec,ni,rxinew,ryinew,rzinew,ovrlap) if (ovrlap) then goto 20 else ii=(ispec-1)*nionmx+ni rx(ii)=rxinew ry(ii)=ryinew rz(ii)=rzinew indspc(ispec,ni)=ii ispcbk(ii)=ni c write(*,*)ii,ispec,ni,rxinew,ryinew,rzinew endif if (ni.lt.ntrg) goto 30 10 continue write(*,*)" Initial configuration made" return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine overlp(ispec,i,rxi,ryi,rzi,ovrlap) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. do jspec=1,ispec if (jspec.eq.ispec) then nj=i-1 else nj=n(jspec) endif if (nj.gt.0) then dda=d2(ispec,jspec) do j=1,nj jj=indspc(jspec,j) rxa=rx(jj) rya=ry(jj) rza=rz(jj) rxia=rxi-rxa ryia=ryi-rya rzia=rzi-rza rxia=rxia-dnint(rxia/xl)*xl ryia=ryia-dnint(ryia/xl)*xl riasq=rxia*rxia+ryia*ryia+rzia*rzia if (riasq.lt.dda) then ovrlap=.true. return endif enddo endif enddo return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine wrcf(istep) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz c write(*,*)' begin wrcf' open(2,file="c.dat") write(2,*) nspec,istep write(2,*) do ispec=1,nspec write(2,*)n(ispec),d(ispec),dca(ispec) write(2,*)nzg(ispec),dzg(ispec) write(2,*) do i=1,n(ispec) ii=indspc(ispec,i) write(2,*)rx(ii),ry(ii),rz(ii) enddo write(2,*) enddo close(2) return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine saves(istep) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz character*8 fngz dimension cz(nspcmx,nzgmx),dinteg(nspcmx),gamma1(nspcmx), : gamma2(nspcmx) open(2,file="out.dat") write(*,'(72("#"))') write(*,*) write(2,*) write(*,*)"--- Save ----, Iteration: ",iter," ---- Step: ", : istep,"------" write(*,'(72("_"))') write(2,'(72("#"))') write(2,*)"--- Save ----, Iteration: ",iter," ---- Step: ", : istep,"------" write(2,'(72("_"))') write(2,'(" zl [A] = ",f10.5)')zl write(2,'(" xl [A] = ",f10.5)')xl write(2,'(72("-"))') write(*,'(72("-"))') if (runit) then write(*,*)" Left electrode charge = ",sigma(1) write(*,*)" Right electrode charge = ",sigma(2) write(2,*)" Left electrode charge = ",sigma(1) write(2,*)" Right electrode charge = ",sigma(2) else write(*,*)" Left electrode charge = ",sigma(1) write(*,*)" Right electrode charge = ",sigma(2) write(2,*)" Left electrode charge = ",sigma(1) write(2,*)" Right electrode charge = ",sigma(2) write(2,'(72("-"))') write(*,'(72("-"))') c write(*,*)" Left electrode charge [C/m2] = ",sigsi(1) c write(*,*)" Right electrode charge [C/m2] = ",sigsi(2) c write(2,*)" Left electrode charge [C/m2] = ",sigsi(1) c write(2,*)" Right electrode charge [C/m2] = ",sigsi(2) endif write(2,'(72("-"))') write(*,'(72("-"))') write(2,*) do igr=1,nsalt write(2,'(a4,": ctarg = ",f6.3,", chemp = ",f10.5, : ", excess = ",f10.6)') : fsalt(igr),ctarg(igr),chemp(igr),chex(igr) enddo write(*,'(12x," N "," q [e] "," d "," dca" )') write(2,'(12x," N "," q [e] "," d "," dca" )') do i=1,nspec write(2,'(a2," ion: ",1x,i7,2x,f5.2,2x,2(2x,f8.4))') : fion(i),n(i),z(i),d(i),dca(i) write(*,'(a2," ion: ",1x,i7,2x,f5.2,2x,2(2x,f8.4))') : fion(i),n(i),z(i),d(i),dca(i) enddo c ##### Calculate and write out concentration profiles ############### do ispec=1,nspec fngz = fgz(ispec) ri=d(ispec)/2.0 open(22,file=fngz) do i=1,nzg(ispec) vslab = xl**2*dzg(ispec) cz(ispec,i) = tosi*gz(ispec,i)/athist/vslab rgz=-zl+ri+(i-0.5)*dzg(ispec) write(22,*)rgz,cz(ispec,i),dzg(ispec) enddo close(22) enddo c ##### Calculate bulk concentration and ############################# c ##### number of ions in the filter ################################# write(*,'(72("_"))') write(2,'(72("_"))') vbulk=xl**2*zl do ispec=1,nspec cbulk(ispec)=tosi*an(ispec)/athist/vbulk write(*,'(" Bulk concentration of ",a2," = ",f10.5)') ' fion(ispec),cbulk(ispec) write(2,'(" Bulk concentration of ",a2," = ",f10.5)') ' fion(ispec),cbulk(ispec) enddo write(*,'(72("_"))') write(2,'(72("_"))') write(*,*)" Average electrode charge 1 = ",asigma(1)/athist write(*,*)" Average electrode charge 2 = ",asigma(2)/athist c #### Salt data #################################################### write(*,'(72("_"))') write(2,'(72("_"))') write(2,'(10x,1x," mu OLD",2x,3x,"mu NEW", : 5x,"Excess OLD",3x,"excess NEW ")') write(*,'(10x,1x," mu OLD",2x,3x,"mu NEW", : 5x,"Excess OLD",3x,"excess NEW ")') do igc=1,nsalt iv=igcval(igc) rtargi=(ctarg(igc)/tosi)*(ctrgcl/tosi)**iv ispec=isalt(igc) rbulki=cbulk(ispec)/tosi rblkcl=cbulk(1)/tosi chexnw=chemp(igc)-dlog(rbulki*rblkcl**iv) chempn=dlog(rtargi)+chexnw write(2,'(4x,a4,": ",4(1x,f10.6))')fsalt(igc),chemp(igc), : chempn,chex(igc),chexnw write(*,'(4x,a4,": ",4(1x,f10.6))')fsalt(igc),chemp(igc), : chempn,chex(igc),chexnw enddo c ##### Acceptance ratios ############################################ write(*,'(72("_"))') write(2,'(72("_"))') write(*,*)" Electrode charge exchange: ",achg(2)/achg(1) write(2,*)" Electrode charge exchange: ",achg(2)/achg(1) write(*,'(3x," Move "," Jump ")') write(2,'(3x," Move "," Jump ")') do i=1,nspec write(*,'(a2,":",2(2x,f7.5))') : fion(i), amove(2,i)/amove(1,i),ajump(2,i)/ajump(1,i) write(2,'(a2,":",2(2x,f7.5))') : fion(i), amove(2,i)/amove(1,i),ajump(2,i)/ajump(1,i) enddo do icount=1,nsalt write(*,'(72("_"))') write(*,'(" Create ",a4," : " : ,f10.0," / ",f10.0," = ",f7.5)') : fsalt(icount), : acreat(icount,2),acreat(icount,1), : acreat(icount,2)/acreat(icount,1) write(*,'(" Destroy ",a4," : " : ,f10.0," / ",f10.0," = ",f7.5)') : fsalt(icount), : adest(icount,2),adest(icount,1), : adest(icount,2)/adest(icount,1) write(2,'(72("_"))') write(2,'(" Create ",a4," : " : ,f10.0," / ",f10.0," = ",f7.5)') : fsalt(icount), : acreat(icount,2),acreat(icount,1), : acreat(icount,2)/acreat(icount,1) write(2,'(" Destroy ",a4," : " : ,f10.0," / ",f10.0," = ",f7.5)') : fsalt(icount), : adest(icount,2),adest(icount,1), : adest(icount,2)/adest(icount,1) enddo write(2,'(72("_"))') write(*,'(72("_"))') pot1=0.0 do ispec=1,nspec nzg2=nzg(ispec)/2 dinteg(ispec)=0.0 gamma1(ispec)=0.0 do i=1,nzg2 zg1=dca(ispec)+(i-1)*dzg(ispec) zg2=zg1+dzg(ispec) dinteg(ispec)=dinteg(ispec) : +(cz(ispec,i)/tosi)*(zg2**2-zg1**2)/2.0 gamma1(ispec)=gamma1(ispec) : +(cz(ispec,i)-cbulk(ispec))*dzg(ispec) enddo pot1=pot1+z(ispec)*dinteg(ispec) gamma1(ispec)=(gamma1(ispec)/tosi)*xl**2 enddo pot1=-pot1*2*twopi*qstar**2 write(*,*)" Potential at left wall /kT = ",pot1 write(2,*)" Potential at left wall /kT = ",pot1 pot2=0.0 do ispec=1,nspec nzg2=nzg(ispec)/2 dinteg(ispec)=0.0 gamma2(ispec)=0.0 do i=1,nzg2 zg1=dca(ispec)+(i-1)*dzg(ispec) zg2=zg1+dzg(ispec) ii=nzg(ispec)-i+1 dinteg(ispec)=dinteg(ispec) : +(cz(ispec,ii)/tosi)*(zg2**2-zg1**2)/2.0 gamma2(ispec)=gamma2(ispec) : +(cz(ispec,ii)-cbulk(ispec))*dzg(ispec) enddo pot2=pot2+z(ispec)*dinteg(ispec) gamma2(ispec)=(gamma2(ispec)/tosi)*xl**2 enddo pot2=-pot2*2*twopi*qstar**2 write(*,*)" Potential at right wall /kT = ",pot2 write(2,*)" Potential at right wall /kT = ",pot2 if (.not.runit) then write(2,'(72("-"))') write(*,'(72("-"))') write(*,*)" Potential at left wall /mV = ",pot1*tomv*1000 write(2,*)" Potential at left wall /mV = ",pot1*tomv*1000 write(*,*)" Potential at right wall /mV = ",pot2*tomv*1000 write(2,*)" Potential at right wall /mV = ",pot2*tomv*1000 endif write(2,'(72("_"))') write(*,'(72("_"))') c1=sigma(1)/(tstar*pot1) c2=sigma(2)/(tstar*pot2) write(*,*)" Capacitance of left DL = ",c1 write(2,*)" Capacitance of left DL = ",c1 write(*,*)" Capacitance of right DL = ",c2 write(2,*)" Capacitance of right DL = ",c2 if (.not.runit) then write(2,'(72("-"))') write(*,'(72("-"))') write(*,*)" Capacitance of left DL [F/m2] = ", : sigsi(1)/(tomv*pot1) write(2,*)" Capacitance of left DL [F/m2] = ", : sigsi(1)/(tomv*pot1) write(*,*)" Capacitance of right DL [F/m2] = ", : sigsi(2)/(tomv*pot2) write(2,*)" Capacitance of right DL [F/m2]= ", : sigsi(2)/(tomv*pot2) endif write(2,'(72("_"))') write(*,'(72("_"))') chgn1=0.0 chgn2=0.0 do ispec=1,nspec write(*,*)" Adsorption of ion ",ispec, : " on left wall = ", : gamma1(ispec) write(*,*)" Adsorption of ion ",ispec, : " on right wall = ", : gamma2(ispec) chgn1=chgn1+z(ispec)*gamma1(ispec) chgn2=chgn2+z(ispec)*gamma2(ispec) write(2,*)" Adsorption of ion ",ispec, : " on left wall = ", : gamma1(ispec) write(2,*)" Adsorption of ion ",ispec, : " on right wall = ", : gamma2(ispec) enddo write(2,'(72("-"))') write(*,'(72("-"))') chgn1t=(sigma(1)/qstar)*xl**2 chgn2t=(sigma(2)/qstar)*xl**2 write(*,*)" net charge in left DL = ",chgn1,chgn1t write(*,*)" net charge in right DL = ",chgn2,chgn2t write(2,*)" net charge in left DL = ",chgn1,chgn1t write(2,*)" net charge in right DL = ",chgn2,chgn2t write(2,'(72("-"))') write(*,'(72("-"))') do ispec=1,nspec write(*,*)" Relative adsorption of ion ",ispec, : " on left wall = ", : gamma1(ispec)/dabs(chgn1t) write(*,*)" Relative adsorption of ion ",ispec, : " on right wall = ", : gamma2(ispec)/dabs(chgn2t) write(2,*)" Relative adsorption of ion ",ispec, : " on left wall = ", : gamma1(ispec)/dabs(chgn1t) write(2,*)" Relative adsorption of ion ",ispec, : " on right wall = ", : gamma2(ispec)/dabs(chgn2t) enddo write(2,'(72("_"))') write(*,'(72("_"))') return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine zero implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz c_______ ZERO OUT STAFF _____________________________________________ asigma(1)=0.0 asigma(2)=0.0 do ispec=1,nspec amove(1,ispec)=0.0 amove(2,ispec)=0.0 ajump(1,ispec)=0.0 ajump(2,ispec)=0.0 an(ispec)=0.0 do i=1,nzg(ispec) gz(ispec,i) = 0.0 enddo do icount=1,3 acreat(icount,1)=0.0 acreat(icount,2)=0.0 adest(icount,1)=0.0 adest(icount,2)=0.0 enddo enddo athist = 0.0 return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine move(ispec,i) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. c write(*,*)"move" ii=indspc(ispec,i) rxiold = rx(ii) ryiold = ry(ii) rziold = rz(ii) rzold = dabs(rziold) di=d(ispec) ri=di/2.0 if (ranff(mag).lt.ratmov(ispec)) then amove(1,ispec) = amove(1,ispec)+1 rxinew = rxiold + (2.0*ranff(mag)-1.0) * drmax ryinew = ryiold + (2.0*ranff(mag)-1.0) * drmax rzinew = rziold + (2.0*ranff(mag)-1.0) * drmax rxinew = rxinew-dnint(rxinew/xl)*xl ryinew = ryinew-dnint(ryinew/xl)*xl drznew = dabs(rzinew) if (drznew.gt.zl-dca(ispec)) return imove=1 else ajump(1,ispec) = ajump(1,ispec)+1 c %%%%%%%%% Generate new position in box %%%%%%%%%%%%%%%%%%%%% rxinew = (ranff(mag)-0.5)*xl ryinew = (ranff(mag)-0.5)*xl rzinew = (2*ranff(mag)-1)*(zl-dca(ispec)) drznew = dabs(rzinew) imove=0 endif call ergi(ispec,i,rxinew,ryinew,rzinew,deltu,ovrlap) if (ovrlap) return 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 do jspec=1,nspec do j=1,n(jspec) jj=indspc(jspec,j) rii(ii,jj)=riinew(ispec,jj) rii(jj,ii)=riinew(ispec,jj) enddo enddo endif return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine ergi(ispec,i,rxi,ryi,rzi,uii,ovrlap) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. uii=0.0 qi=q(ispec) di=d(ispec) ri=di/2.0 ii=indspc(ispec,i) c ################################################################### c_______ SUM OVER OTHER FREE CHARGES ________________________________ c ################################################################### do jspec=1,nspec ddi=d2(ispec,jspec) qj=q(jspec) nj=n(jspec) do j=1,nj if (ispec.eq.jspec.and.i.eq.j) goto 1 jj = indspc(jspec,j) rxj = rx(jj) ryj = ry(jj) rzj = rz(jj) c_______ DISTANCE IS CALCULATED FOR IONS IN NEW CONFIGURATION ______ c_______ STORED IN RIINEW __________________________________________ rxij = rxi-rxj ryij = ryi-ryj rzij = rzi-rzj rxij = rxij-dnint(rxij/xl)*xl ryij = ryij-dnint(ryij/xl)*xl rijsq= rxij*rxij+ryij*ryij+rzij*rzij if (rijsq.lt.ddi) then ovrlap = .true. return endif rij= dsqrt(rijsq) riinew(ispec,jj)=rij uiinew = qi*qj/rij uii=uii + uiinew zi = dabs(rzij) uii = uii + qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl)) c_______ DISTANCE IS FROM LOOKUP TABLE ______________________________ c_______ FOR OLD CONFIGURATION ______________________________________ rij=rii(ii,jj) uiiold = qi*qj/rij uii=uii - uiiold zi = dabs(rz(ii)-rzj) uii = uii - qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl)) 1 continue enddo enddo uii = uii - qi*sigma(1)*(-twopi*(rz(ii)+zl)) : - qi*sigma(2)*(-twopi*(zl-rz(ii))) uii = uii + qi*sigma(1)*(-twopi*(rzi+zl)) : + qi*sigma(2)*(-twopi*(zl-rzi)) return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine rechg implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. c write(*,*)"move" sigold(1)=sigma(1) sigold(2)=sigma(2) achg(1) = achg(1)+1.0 dsigma = (2.0*ranff(mag)-1.0) * dchgmx sigma(1) = sigold(1) + dsigma sigma(2) = sigold(2) - dsigma call erchg(deltu) if ((deltu.le.0.0).or.(dexp(-deltu).gt.ranff(mag))) then achg(2) = achg(2) + 1.0 else sigma(1)=sigold(1) sigma(2)=sigold(2) endif return end c_______________________________________________________________________ c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ subroutine erchg(uii) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. uii=0.0 c ################################################################### c_______ SUM OVER OTHER FREE CHARGES ________________________________ c ################################################################### do jspec=1,nspec qj=q(jspec) nj=n(jspec) do j=1,nj jj = indspc(jspec,j) rzj = rz(jj) uii = uii - qj*sigold(1)*(-twopi*(rzj+zl)) : - qj*sigold(2)*(-twopi*(zl-rzj)) uii = uii + qj*sigma(1)*(-twopi*(rzj+zl)) : + qj*sigma(2)*(-twopi*(zl-rzj)) enddo enddo uii = uii - (-twopi*2*zl*sigold(1)*sigold(2)*xl**2) uii = uii + (-twopi*2*zl*sigma(1)*sigma(2)*xl**2) uii = uii + voltg*sigold(1)*xl**2 uii = uii - voltg*sigma(1)*xl**2 return end c ################################################################### c ################################################################### c ################################################################### c ################################################################### subroutine create(igc) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. acreat(igc,1) = acreat(igc,1)+1.0 iv=igcval(igc) ntrial = ntot + iv+1 if ( ntrial .ge. ntotmx ) stop 'maximum number of atoms in box' ispec1=ispcin(igc,1) nterm=n(ispec1)+1 do icount=2,iv+1 ispec=ispcin(igc,icount) nterm=nterm*(n(ispec)+icount-1) enddo do icount=1,iv+1 ispec = ispcin(igc,icount) rxinew = (ranff(mag)-0.5)*xl ryinew = (ranff(mag)-0.5)*xl rzinew = (2*ranff(mag)-1)*(zl-dca(ispec)) drznew = dabs(rzinew) rznew(icount) = rzinew rxnew(icount) = rxinew rynew(icount) = ryinew enddo call encrea(igc,deltu,ovrlap) if (ovrlap) then return endif bfcr = dexp(-deltu+chemp(igc))*vgc(igc)/dble(nterm) if (bfcr.gt.ranff(mag)) then call add(igc) acreat(igc,2) = acreat(igc,2)+1.0 endif return end c ################################################################### c ################################################################### c ################################################################### c ################################################################### subroutine destr(igc) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz adest(igc,1) = adest(igc,1) + 1.0 iv=igcval(igc) ispec1=ispcin(igc,1) nterm=n(ispec1) if (nterm.eq.0) return do icount=2,iv+1 ispec=ispcin(igc,icount) nterm=nterm*(n(ispec)-icount+2) enddo do icount=1,iv+1 ispec = ispcin(igc,icount) 1 continue i = int(n(ispec)*ranff(mag)) + 1 ii = indspc(ispec,i) if (icount.gt.2) then do icold=2,icount-1 ispeco = ispcin(igc,icold) iold = iout(icold) iiold = indspc(ispeco,iold) if (ii.eq.iiold) goto 1 enddo endif iout(icount) = i enddo if (iv.gt.1) then if (iout(2).lt.iout(3)) then ipuf=iout(2) iout(2)=iout(3) iout(3)=ipuf endif if (iv.eq.3) then if (iout(3).lt.iout(4)) then ipuf=iout(4) iout(4)=iout(3) iout(3)=ipuf endif if (iout(2).lt.iout(3)) then ipuf=iout(2) iout(2)=iout(3) iout(3)=ipuf endif endif endif do icount=1,iv+1 ispec = ispcin(igc,icount) i = iout(icount) ii = indspc(ispec,i) rxnew(icount) = rx(ii) rynew(icount) = ry(ii) rznew(icount) = rz(ii) enddo call endest(igc,deltu,ovrlap) bfdr = dexp(-deltu-chemp(igc))*dble(nterm)/vgc(igc) if (bfdr.gt.ranff(mag)) then call remove(igc) adest(igc,2) = adest(igc,2)+1.0 endif return end c################################################################### c################################################################### c################################################################### subroutine remove(igc) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz iv=igcval(igc) do icount=1,iv+1 ispec=ispcin(igc,icount) i=iout(icount) ii = indspc(ispec,i) if (i.lt.n(ispec)) then do k = i+1,n(ispec) indspc(ispec,k-1) = indspc(ispec,k) ispcbk(indspc(ispec,k-1))=k-1 enddo indspc(ispec,n(ispec)) = ii ispcbk(indspc(ispec,n(ispec)))=n(ispec) endif n(ispec) = n(ispec)-1 enddo ntot=ntot-iv-1 return end c################################################################### c################################################################### c################################################################### subroutine add(igc) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz dimension inew(4) iv=igcval(igc) do icount=1,iv+1 ispec=ispcin(igc,icount) inew(icount) = n(ispec) + 1 ii = indspc(ispec,inew(icount)) if (ii.eq.0) then ii = (ispec-1)*nionmx+inew(icount) indspc(ispec,inew(icount)) = ii ispcbk(ii)=inew(icount) endif n(ispec) = inew(icount) rx(ii)=rxnew(icount) ry(ii)=rynew(icount) rz(ii)=rznew(icount) enddo ntot=ntot+iv+1 do icount=1,iv+1 ispec=ispcin(igc,icount) i=inew(icount) ii = indspc(ispec,i) do jspec=1,nspec do j=1,n(jspec) jj=indspc(jspec,j) if (ii.ne.jj) then rii(ii,jj)=riinew(icount,jj) rii(jj,ii)=riinew(icount,jj) endif enddo enddo enddo do icount=1,iv+1 ispec=ispcin(igc,icount) ii = indspc(ispec,inew(icount)) do jcount=1,iv+1 if (icount.ne.jcount) then jspec=ispcin(igc,jcount) jj = indspc(jspec,inew(jcount)) rii(ii,jj)=riicrn(icount,jcount) endif enddo enddo return end c##################################################################### c##################################################################### c##################################################################### subroutine encrea(igc,uii,ovrlap) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz logical ovrlap ovrlap=.false. uii=0.0 iv=igcval(igc) do icount1=1,iv ispec1=ispcin(igc,icount1) q1=q(ispec1) do icount2=icount1+1,iv+1 ispec2=ispcin(igc,icount2) q2=q(ispec2) ddi=d2(ispec1,ispec2) rxij = rxnew(icount1)-rxnew(icount2) ryij = rynew(icount1)-rynew(icount2) rzij = rznew(icount1)-rznew(icount2) rxij = rxij-dnint(rxij/xl)*xl ryij = ryij-dnint(ryij/xl)*xl rijsq= rxij*rxij+ryij*ryij+rzij*rzij if (rijsq.lt.ddi) then ovrlap = .true. return endif rij=dsqrt(rijsq) riicrn(icount1,icount2)=rij riicrn(icount2,icount1)=rij uii=uii + q1*q2/rij zi = dabs(rzij) uii = uii + q1*(q2/xl**2) : *(-twopi*zi-wpot(zi,xl)) enddo c uii = uii + q1*sigma(1)*(-twopi*(rznew(icount1)+zl)) c : + q1*sigma(2)*(-twopi*(zl-rznew(icount1))) enddo do icount=1,iv+1 ispec=ispcin(igc,icount) qi=q(ispec) do jspec=1,nspec qj=q(jspec) ddi=d2(ispec,jspec) do j=1,n(jspec) jj=indspc(jspec,j) rxij = rxnew(icount)-rx(jj) ryij = rynew(icount)-ry(jj) rzij = rznew(icount)-rz(jj) rxij = rxij-dnint(rxij/xl)*xl ryij = ryij-dnint(ryij/xl)*xl rijsq= rxij*rxij+ryij*ryij+rzij*rzij if (rijsq.lt.ddi) then ovrlap = .true. return endif rij=dsqrt(rijsq) riinew(icount,jj)=rij uii=uii + qi*qj/rij zi = dabs(rzij) uii = uii + qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl)) enddo enddo uii = uii + qi*sigma(1)*(-twopi*(rznew(icount)+zl)) : + qi*sigma(2)*(-twopi*(zl-rznew(icount))) enddo return end c################################################################### c##################################################################### c##################################################################### subroutine endest(igc,uii,ovrlap) implicit doubleprecision (a-h,o-z) parameter(nionmx=2000,ntotmx=12000,ntrymx=12000, : nspcmx=6,nzgmx=4000) 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),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), : nzg(nspcmx),fgz(nspcmx) common/accu/ an(nspcmx),cbulk(nspcmx),asigma(2) common/adj/ amove(2,nspcmx),ajump(2,nspcmx), : acreat(3,2),adest(3,2),achg(2) common/bufer/ riinew(nspcmx,ntotmx) common/simp/ ratchg,dchgmx,nstep,naver,isave,runit logical runit character*2 fion character*4 fsalt character*8 fgz uii=0.0 iv=igcval(igc) do icount=1,iv+1 ispec=ispcin(igc,icount) i=iout(icount) ii=indspc(ispec,i) qi=q(ispec) do jspec=1,nspec qj=q(jspec) ddi=d2(ispec,jspec) do j=1,n(jspec) jj=indspc(jspec,j) if (ii.ne.jj) then rij=rii(ii,jj) uii=uii - qi*qj/rij zi = dabs(rz(ii)-rz(jj)) uii = uii - qi*(qj/xl**2)*(-twopi*zi-wpot(zi,xl)) endif enddo enddo uii = uii - qi*sigma(1)*(-twopi*(rz(ii)+zl)) : - qi*sigma(2)*(-twopi*(zl-rz(ii))) enddo do icount1=1,iv ispec1=ispcin(igc,icount1) i1=iout(icount1) ii1=indspc(ispec1,i1) q1=q(ispec1) do icount2=icount1+1,iv+1 ispec2=ispcin(igc,icount2) i2=iout(icount2) ii2=indspc(ispec2,i2) q2=q(ispec2) rij=rii(ii1,ii2) uii=uii + q1*q2/rij zi = dabs(rz(ii1)-rz(ii2)) uii = uii + q1*(q2/xl**2) : *(-twopi*zi-wpot(zi,xl)) enddo c uii = uii - q1*sigma(1)*(-twopi*(rz(ii1)+zl)) c : - q1*sigma(2)*(-twopi*(zl-rz(ii1))) enddo return end c ################################################################### c ################################################################### function wpot(z,wl) implicit doubleprecision (a-h,o-z) common/const/ pi,twopi,sq2,sq3,tosi,tomv,mag x=z/wl x2=x*x r1=dsqrt(0.5+x2) r2=dsqrt(0.25+x2) r3=(0.5+r1)/r2 wpot = 4*dlog(r3) - 4*x*( pi/2.0 - datan(4*x*r1) ) wpot = wpot*wl return end