program MCGC !********************************************************************* ! MONTE-CARLO STUDY OF pores ! ! Version 2008 -H2 in carbon pores ! ! An MC step is the movement of all N (not constant) particles !******************************************************************** include 'global.inc' DIMENSION F(NBmax),FBAR(NBmax) DATA Z049/0.45/,Z051/0.55/,SMALL/0.01/ ! *** READ-IN INPUT procedures ----------------------------------- open(6,file='mc.lst',status='unknown') call input print *, 'input finished' !****************************************************************** ! INITIALIZATION OF CALCULATIONS !****************************************************************** ! ! UU(I,J) IS THE STORED INTERACTION ENERGY BETWEEN ! MOLECULES OF TYPE I,J. THIS IS A SYMMETRIC ARRAY. ! WHEN MOLECULE I IS MOVED , ONLY THE (N-1) INTERACTIONS ! (I.NE.J) NEED TO BE UPDATED. ! !****************************************************************** call initial print *, 'initial finished' !****************************************************************** ! Verification of the reading: write data into "mc.lst" file ! and write positions into "pos.ini" call outControl print *, 'control_out finished' ! *** *************************************************************** ! create the file 'mc.pos' with box and molecules coordinates ! after each bin open(7,file='mc.pos',status='unknown') ! *** ************************************************************** ! create a file with (after each bin) IB,UGt/N,ESM/N,EPAIR/N,UTOT/N open(9,file='mc_ene.dat',status='unknown') ! *** ************************************************************** ! create a file with IB,N,HINT/n,HBAR/n,f(ib)/n,vbar/n open(10,file='mc_ent.dat',status='unknown') ! *** create a file of problem reports ------------------- ! open(11,file='mc_irregul.out',status='unknown') ! *** re-build list of neighbors ------------------------ CALL LSTMAKE ! already done in 'initial', but any way .... ! *** *********** MAIN LOOP ************************************ ! over the bins HT = 0.0 ! total energy WT = 0.0 ! wall energy VT = 0.0 ! volume print *, 'Begin NBIN loop' do IB = 1,NBIN CALL LSTMAKE ! *** ------- intermediate saving ----------------------------- ind_res = mod(IB,1000) if(ind_res.eq.0) then call restart print *, IB ,Npart,UGT/N,ENTH/Npart,float(ITT)/(IAT+1) endif ! *** ------------ bin accumulates ---------------------------- WB = 0.0 HB = 0.0 VB = 0.0 ibn = 0 ! index of steps: ! becomes =0 after UPDATENB steps ! *** here starts the loop over the steps --------------------- DO IS = 1,NSTEP !*** STEP LOOP STARTS *** ! reactualize the list of neighbors every UPDATENB steps ibn=ibn + 1 if (ibn.eq.UPDATENB)then call update ibn = 0 endif ! move the particles do I = 1,ncycle !*** CYCLE LOOP STARTS *** random=rndm() ifnum=INT(random*(npart+nexc))+1 ! Select a target molecule if((ifnum.le.npart).and.(npart.gt.1))then ITT = ITT + 1 CALL translate(ifnum) ! translational move if((n_Molec.ge.2).and.(N_if_orn.ne.0)) then CALL rotate(ifnum) ! orientational move endif else if (MX.ne.0) then ITE = ITE + 1 nold = npart ! save number of particles call exchange ! exchange particle endif endif ! print *, 'End of cycle:',I,IS,IB enddo ! *** ncycle loop ends **************************** ! *** ****************************************************************** ! compute the running sum of H,ULJ and V in a bin HB = HB + ENTH/N ! /N - last change to avoid sum WB = WB + UGT/N ! of energy of different number of particles VB = VB + VOL N = Npart ! this is important: the number of particles ! might has changed !!!!!!! enddo ! *** END of the loop over the STEPS ***************** ! *** ***************************************************************** ! compute the average value of H and V in the bin ! IB is the number of bins already done WBB = WB / NSTEP HBB = HB / NSTEP VBB = VB / NSTEP F(IB) = HBB ! cumulate H, and V and their averaged values ! over the bins already done HT = HT + HB WT = WT + WB VT = VT + VB IT2 = IB * NSTEP FBAR(IB) = HT / IT2 VBAR = VT / IT2 ! write the state of the system after IB bins WRITE (6,610) IB ,N 610 FORMAT (' AFTER',I5,' BINS, THE STATE OF ',i5,' ATOMS IS') WRITE (6,611) FBAR(IB) 611 FORMAT (5X,'HCUM=',E16.10,5X) WRITE (6,641) F(IB),VBB 641 FORMAT (5X,'HBIN=',E16.10,5X,'VBIN=',E16.10,5X) ! WRITE (6,612) ENTH,VOL ! 612 FORMAT (5X,'ENTH=',E16.10,5X,'VOLUM= ',E16.10,5X) write(7,'(2(2x,i6),a)') IB, N,' Nbin Npart' do i = 1,N write(7,'(6f12.6,i3)')xc(i),yc(i),zc(i) enddo HINT = ENTH/N VINT = VOL * 0.602252/N WRITE (6,6912) HINT,VINT 6912 FORMAT (5X,'H/N= ',E16.10,5X,'V/N= ',E16.10) HBAR = FBAR(IB) VBAR = VBAR * 0.602252 / N WRITE (6,615) ITT,IAT,IACT,ITO,IAO,IACO,ITE,IAA,IAR 615 FORMAT (3(1X,I9)) IF (STEPTST.GT.SMALL) then TESTT = FLOAT(IAT) / FLOAT(ITT) TDSTEP = TD / 10. IF (TESTT.GE.Z051) TD = TD + TDSTEP IF (TESTT.LE.Z049) TD = TD - TDSTEP TD2 = TD + TD if(n_Molec.ge.2) then TESTO = FLOAT(IAO) / FLOAT(ITO) ODSTEP = OD / 10. IF (TESTT.GE.Z051) OD = OD + ODSTEP IF (TESTT.LE.Z049) OD = OD - ODSTEP OD2 = OD + OD endif endif WRITE (6,'(3(2x,f12.6))') TD,OD WRITE (6,616) 616 FORMAT (1X,43('-*-')) WRITE (9,8000) IB,N,UGT/N,ENTH/N WRITE (10,8000) IB,N,ENTH/N,HBB,HBAR 8000 FORMAT (1X,2(2x,I5),8(1X,E16.10)) enddo ! end of the loop over the bins **************** * *** ******* MAIN LOOP FINISHED ************************** WRITE (6,617) 617 FORMAT (1X,'BIN NO.',3X,'BIN ENTHALPY',4X,'AV. ENTHALPY') DO I = 1,NBIN WRITE (6,618) I,F(I),FBAR(I) enddo 618 FORMAT(1X,I6,2E17.6) close(6) close(7) close(9) close(10) close(11) ! *** output data for further calculations to the file "restart.res" call restart END * *** ***************************************************************