subroutine translate(IL) INCLUDE 'global.inc' ! move the molecule (translational move) ! XCT,YCT,ZCT are the trial positions ! UDIF = H(OLD) - H(NEW) ! GDIF = EGG(old) - EGG(new) XCT = XC(IL) + TD - TD2 * RNDM() YCT = YC(IL) + TD - TD2 * RNDM() ZCT = ZC(IL) + TD - TD2 * RNDM() ! *** Calculate UDIF --------------------------------------------------- call period_check(XCT,YCT,xcb,ycb,iper) ! check if in box (x,y) CALL POTgr(XCb,YCb,ZCT,EG) ! calculate interaction with the wall call ener1(XCT,YCT,ZCT,entot,IL) ! calculates interaction with others GDIF = EGG(IL)-EG EDIF = (ENMOL(IL)-EGG(IL)) - entot UDIF = GDIF+ EDIF ! accumulate the difference between new and old ENMOLT(IL) = entot + EG max = NABL1(IL) * if(max.ne.0) then * DO LL= 1,max * L = LIST1(IL,LL,4) * UDIF = UDIF + UU(L,IL) - UTT(L,IL) * enddo * * endif ! *** UDIF calculated ; check acceptence ------------------------- IF (UDIF.le.0.0) then XR = RNDM() RTEST = TEMP * LOG(XR) IF (UDIF.LT.RTEST) RETURN IACT = IACT + 1 endif IAT = IAT + 1 ! *** store the new accepted position *************************** XC(IL) = XCb YC(IL) = YCb ZC(IL) = ZCT ENTH = ENTH - EDIF - GDIF UGT = UGT - GDIF EGG(IL) = EG ENMOL(IL) = ENMOLT(IL) ! DO I = 1,Npart ! UU(I,IL) = UTT(I,IL) ! UU(IL,I) = UU(I,IL) ! enddo if(iper.eq.1) call LSTMAKE ! accepted move shifted outside the zbox ! if(ENMOL(IL).gt.0.0) then ! write(11,*) ' positive, accepted tranlations move' ! write(11,*)entot,EG,il ! call prints(il) ! stop 'Energy ENMOL detected positive !!!' ! endif END ! *** ************************************************ subroutine rotate(IL) INCLUDE 'global.inc' ! move the molecule (translational move) ! XCT,YCT,ZCT are the trial positions ! UDIF = H(OLD) - H(NEW) ! GDIF = EGG(old) - EGG(new) 700 XRU = 1. - 2.0*RNDM() XRV = 1. - 2.0*RNDM() XRW = 1. - 2.0*RNDM() DS2 = XRU * XRU + XRV * XRV + XRW * XRW IF (DS2.GT.1.0) GOTO 700 UCT = U(IL) + OD * XRU VCT = V(IL) + OD * XRV WCT = W(IL) + OD * XRW DS = SQRT( UCT * UCT + VCT * VCT + WCT * WCT ) UCT = UCT / DS VCT = VCT / DS WCT = WCT / DS ! CALL ENER1 (XC(IL),YC(IL),ZC(IL),UCT,VCT,WCT) ! CALL POTgr(XCb,YCb,ZCT,EG) ! calculate interaction with the wall ! call ener1(XCb,YCb,ZCT,entot,IL) ! calculates interaction with others IF (EDIF.lt.0.0) then XR = RNDM() RTEST = RB * LOG(XR) IF (EDIF.LT.RTEST) RETURN IACO = IACO + 1 endif IAO=IAO + 1 U(IL) = UCT V(IL) = VCT W(IL) = WCT ENTH = ENTH - UDIF UGT = UGT - GDIF EGG(IL) = EG ENMOL(IL) = entot+EG ! DO I = 1,Npart ! UU(I,IL) = UTT(I,IL) ! UU(IL,I) = UU(I,IL) ! enddo end ! *** ************************************************ subroutine period_check(XCT,YCT,xcb,ycb,iper) INCLUDE 'global.inc' real XCT,YCT,xcb,ycb iper = 0 xcb = xct if(xct.gt.xpore) xcb=xct - xpore if(xct.lt.0.) xcb=xct + xpore if(xcb.ne.xct)iper = 1 ycb = yct if(yct.gt.ypore) ycb=yct - ypore if(yct.lt.0.) ycb=yct + ypore if(ycb.ne.yct)iper = 1 end ! *** ************************************************ subroutine exchange INCLUDE 'global.inc' volch = vol - 2.*1.5*xpore*ypore if ((rndm().lt.0.5).and.(npart.gt.2)) then ! remove !!! ip=INT(npart*rndm())+1 if(ip.gt.npart) ip=npart eno = ENmol(ip) if(eno.gt.10.) then ! this is automatic remove when molecular ! call prints(ip) ! energy positive (although theoretically, IAR = IAR + 1 ! not possible) call update_remove(ip) return endif coef = npart/(zz*volch) xr = rndm() ! print *,xr, log(xr) rrnd = log(xr) test = log(coef) + eno/TEMP if (rrnd.lt.test) then IAR = IAR + 1 call update_remove(ip) endif else ! insert !!! call pos_gen(xct,yct,zct,xpore,ypore,zpore) ! geomtry of ! inserting if(npart.ge.M) print *, 'Npart out of limit !!!' If(npart.gt.0) then call ener0(xct,yct,zct,enn) else enn = 0.0 endif ! print *, enn ,npart call potgr(xct,yct,zct,EG) ent = enn + EG coefi = zz*volch/(npart+1) xr = rndm() rrnd = log(xr) test = log(coefi) - ent/TEMP if (rrnd.lt.test) then IAA = IAA +1 xc(npart+1) = xct yc(npart+1) = yct zc(npart+1) = zct call update_add endif endif N = Npart 101 format(/,a,i7,5x,a,i7) end *************************************************************************** * subroutine update_remove(il) INCLUDE 'global.inc' xc(il) = xc(npart) yc(il) = yc(npart) zc(il) = zc(npart) Npart = Npart -1 if(Npart.lt.0) Npart=0 call update end !************************************************************************** subroutine update_add INCLUDE 'global.inc' Npart = Npart + 1 call update end !************************************************************************** function rndm() ! call random_number(rndm) rndm = rand_rec(1) ! look at initial to start generator end ! *** ****************************** function rand_rec(idum) dimension r(97) parameter (m1=259200, ia1=7141, ic1=54773, rm1=1./m1) parameter (m2=134456, ia2=8121, ic2=28411, rm2=1./m2) parameter (m3=243000, ia3=4561, ic3=51349) data iff /0/ save if(idum.lt.0.or.iff.eq.0)then iff = 1 ix1 = mod(ic1-idum,m1) ix1 = mod(ia1*ix1+ic1,m1) ix2 = mod(ix1,m2) ix1 = mod(ia1*ix1+ic1,m1) ix3 = mod(ix1,m3) do j = 1,97 ix1 = mod(ia1*ix1+ic1,m1) ix2 = mod(ia2*ix2+ic2,m2) r(j) = (float(ix1)+float(ix2)*rm2)*rm1 enddo endif ix1 = mod(ia1*ix1+ic1,m1) ix2 = mod(ia2*ix2+ic2,m2) ix3 = mod(ia3*ix3+ic3,m3) j = 1 + (97*ix3)/m3 if(j.lt.1.or.j.gt.97)j = mod((iabs(j)+1),97) rand_rec = r(j) r(j) = (float(ix1) + float(ix2)*rm2)*rm1 end ! *** ! ************************************************************ ! *** subroutine pos_gen(x,y,z,a,b,c) ! This subroutine is only for call of different structures of pores call gen_in_slit(x,y,z,a,b,c) end ! *** ******************************************************** subroutine gen_in_slit(x,y,z,a,b,d) real x,y,z,a,b,d, dd data dd/1.5/ ! molecules generated 1.5 A z = dd + rndm()*(d-2.*dd) ! from the wall y = rndm()*b x = rndm()*a end ! *** ! ************************************************************ ! *** subroutine prints(ip) INCLUDE 'global.inc' write(11,*) '*******' write(11,*) ib,ip ! do j = 1,npart ! write(11,'(i4,3(2x,f15.4))')j,egg(j),enmol(j),uu(j,ip) ! enddo end