C calculate BRW value of average energy transfer C AUTHOR: GILBERT. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION TITLE(20) DOUBLE PRECISION M,NU,K,MLIGHT,MAV,MBG,MM,KTWVN C THE FOLLOWING DATA ARE LOCAL EPSILON AND SIGMA VALUES C FOR INDIVIDUAL ATOMS DATA SIGH/3./,SIGC/3.2/,SIGO/3.2/,SIGF/3.2/,SIGCL/3.4/, 1 SIGS/3.4/,SIGBR/3.6/,SIGI/4.0/,SIGN/3.2/, 2 EPSH/6.5/,EPSC/20.3/,EPSO/20.3/,EPSF/20.3/,EPSCL/120./, 3 EPSS/120./,EPSBR/190./,EPSI/230./,EPSN/20.3/ C OPEN(UNIT=5,STATUS='OLD',FILE='BRW.DAT') C OPEN(UNIT=6,STATUS='NEW',FILE='BRW.OUT') WRITE(6,1000) 1000 FORMAT(10X,20(2H *),/, 1 ' BIASED RANDOM WALK MODEL B. VERSION: June 29 1991') READ(5,789) TITLE 789 FORMAT(20A4) WRITE(6,789) TITLE DELR = 0.0005 DELR = DELR*1E-10 300 CONTINUE READ(5,*) EPS,SIGMA,M,T,NATOMS,MLIGHT IF(EPS.LT.0.) GO TO 987 READ(5,*) EPSBG,SIGBG,MBG WRITE(6,1) EPSBG,SIGBG,MBG 1 FORMAT(///,' BATH GAS: EPSILON (K) =',F10.2,/, 1 ' SIGMA (A) =',F10.2,/, 2 ' MASS (AMU) =',F10.2) C FIND KT IN CM-1 KTWVN = T/1.439 WRITE(6,2) NATOMS 2 FORMAT(' NUMBER OF DISTINCT ATOMS:',I5) MAV=0. EPSRL = 0. SIGRL=0. NTOT = 0 WRITE(6,700) 700 FORMAT(/,' INPUT ATOMS:',/, 1 4X,' MASS STOICH.CFT. LOCAL EPS. LOCAL SIGMA') DO 4 I=1,NATOMS READ(5,*) MM,NST,EPSL,SIGL IF(EPSL.GT.0.) GO TO 704 C ASSIGN DEFAULT VALUES TO LOCAL EPSILON, SIGMA IF(MM.GT.3.5) GO TO 500 C HYDROGEN OR AN ISOTOPE EPSL = EPSH SIGL = SIGH GO TO 704 500 IF(MM.GT.13.5) GO TO 501 C CARBON EPSL = EPSC SIGL = SIGC GO TO 704 501 IF(MM.GT.15.) GO TO 502 C NITROGEN EPSL = EPSN SIGL = SIGN GO TO 704 502 IF(MM.GT.18.5) GO TO 503 C OXYGEN EPSL = EPSO SIGL = SIGO GO TO 704 503 IF(MM.GT.20.) GO TO 504 C FLUORINE EPSL = EPSF SIGL = SIGF GO TO 704 504 IF(MM.GT.34.) GO TO 505 C SULFUR EPSL = EPSS SIGL = SIGS GO TO 704 505 IF(MM.GT.38.) GO TO 506 C CHLORINE EPSL = EPSCL SIGL = SIGCL GO TO 704 506 IF(MM.GT.81.) GO TO 507 C BROMINE EPSL = EPSBR SIGL = SIGBR GO TO 704 C ONLY IODINE IS LEFT 507 EPSL = EPSI SIGL = SIGI 704 CONTINUE WRITE(6,703) MM,NST,EPSL,SIGL 703 FORMAT(F10.2,I10,5X,F10.2,5X,F10.2) FNST = NST MAV = MAV + FNST*MM NTOT = NTOT+NST EPSRL = EPSRL + FNST*EPSL SIGRL = SIGRL + FNST*SIGL 4 CONTINUE TOT = NTOT MAV = MAV/TOT C MODEL A MAV = MAV*MBG/(MAV+MBG) C MODEL B: MAV = 1./((1./(MAV*TOT-MAV)) + (1./MAV)) ELOC = EPSRL/TOT SLOC = SIGRL/TOT ELOC = DSQRT(ELOC*EPSBG) SLOC = 0.5*(SLOC+SIGBG) READ(5,*) ECOLL,NVIB,NU TSTAR = T/EPS OM22 = 1./(0.636 + 0.567*DLOG10(TSTAR)) WRITE(6,3) EPS,SIGMA,M,T,OM22,MAV,MLIGHT ERAT = ECOLL*1.439/(T*NVIB) 3 FORMAT(/,' EPSILON (K) =',T23,F10.0,/,' SIGMA (A) =',T23, 1 F10.2,/,' RED. MASS (A.M.U.) =',T23,F10.2,/, 2 ' TEMPERATURE (K) =',T23,F10.0,/, 3 ' OMEGA 22 =',T23,F10.3, 4 /,' AVERAGE ATOM/ATOM MASS =',F10.2,' AMU', 5 /,' MASS LIGHTEST ATOM =',F10.2,' AMU') WRITE(6,101) ECOLL,NVIB,NU 101 FORMAT(' COLLISION ENERGY (CM-1):',F10.0,/, 1 ' (NOT USED IN MODEL B)',//, 2 ' NO. VIBRATIONS:',I5, ' (NOT USED IN MODEL B)', 3 /,' HIGHEST FREQUENCY (CM-1):', F10.0) NU = NU*3E10 K = 39.5*(MLIGHT*1.66E-27)*NU*NU WRITE(6,31) SLOC,ELOC 31 FORMAT(' LOCAL SIGMA, EPSILON',2F8.2) SLOC=SLOC*1E-10 ELOC=ELOC*1.38E-23 C MODEL A EBAR = ECOLL*1.986E-23/NVIB C THIS PROGRAM IS FOR MODEL B C MODEL B: EBAR = 2.*KTWVN*1.986E-23 EPS = EPS*1.38E-23 SIGMA = SIGMA*1E-10 DELTAX = DSQRT(2.*EBAR/K) C = 2.*3.14159*NU C WRITE(6,78) C C 78 FORMAT(' C =',1P,E10.2,' S-1') M = M*1.66E-27 E = 2.*1.38E-23*T D= SIGMA*OM22 D = DMAX1(D,SIGMA) BOND=0.666667 B = BOND*D C B IS AVERAGE IMPACT PARAMETER C WRITE(6,51) B*1E10 C 51 FORMAT(//,' IMPACT PARAMETER (A):',F5.3) N=0 FAVSLW=0. TAUC = 0. R = D 5 SONR = SIGMA/R SONR6 = SONR**6 SONR12 = SONR6*SONR6 VEFF = 4.*EPS*(SONR12-SONR6) + E*((B/R)**2) IF(VEFF.GE.E) GO TO 10 TAUC = TAUC +(1./SQRT(E-VEFF)) R = R-DELR GO TO 5 10 TAUC = TAUC*DELR*SQRT(2.*M) C WRITE(6,52) R*1E10 C 52 FORMAT(' TURNING PT (A) =',F6.3) WRITE(6,9) TAUC 9 FORMAT(/,' COLLISION TIME (S) =',1P,E9.2) D= SLOC*OM22 D = DMAX1(D,SLOC) BOND=0.666667 B = BOND*D N=0 FAVFST=0. E=EBAR R = D 6 SONR = SLOC/R SONR6 = SONR**6 SONR12 = SONR6*SONR6 VEFF = 4.*ELOC*(SONR12-SONR6) + E*((B/R)**2) IF(VEFF.GE.E) GO TO 11 FORCE = 4.*ELOC*(-12.*SONR12+6.*SONR6) -2.*E*((B/R)**2) FORCE = FORCE/R FAVFST=FAVFST+DABS(FORCE) N=N+1 R = R-DELR GO TO 6 11 CONTINUE C WRITE(6,52) R*1E10 FORCE = 4.*ELOC*(-12.*SONR12+6.*SONR6) -2.*E*((B/R)**2) FORCE = FORCE/R C WRITE(6,21) FORCE C 21 FORMAT(/,' FAST FORCE AT TURNING POINT (NEWTON) =',1P,E9.2) FFAST=DABS(FORCE) FAVFST = FAVFST/N F = 1.20 C WRITE(6,1002) F C1002 FORMAT(' SIGMA/R FOR :',F10.3) FTERM = DABS((4.)*(-12.*(F**12) + 6.*(F**6))) X = SLOC/F DELTAV = FTERM*ELOC*DELTAX/X IF(DELTAV.GT.0.5*EBAR) DELTAV = 0.5*EBAR EIDOT = (EBAR-DELTAV)*NU ECM = EIDOT/1.986E-23 ECM2=ECM*ECM DVCM = DELTAV/1.986E-23 C WRITE(6,86) DVCM C86 FORMAT(' DELTA V (CM-1):',F10.1) C WRITE(6,120) ECM2 C120 FORMAT(' WITH MODIFIED DELTA V',1P,E10.2, C 1 ' CM-2 S-2') EIDOT2=EIDOT*EIDOT C WRITE(6,71) FAVFST C 71 FORMAT(' AVERAGE FAST FORCE =',1P,E9.2) FAV = FFAST FAV = DABS(FAV) C WRITE(6,127) FAV C127 FORMAT(' AVERAGE FORCE USED:',1P,E9.2,' NEWTON') A = FAV/SQRT(MAV*1.66E-27*EBAR*0.5) C WRITE(6,2233) A C2233 FORMAT(' CALC. A =',1P,E9.2,' S-1') S = 2.*A*EIDOT2*TAUC/(A*A+C*C) S = DSQRT(S)/1.986E-23 C THIS VALUE OF S IS THAT WITHOUT THE SEMI-EMPIRICAL CORRECTION WRITE(6,246) S 246 FORMAT(' S = ', 1 F10.1,' CM-1',/,10x,'**********') IF(S.GT.1000.) WRITE(6,387) 387 FORMAT(' WARNING: A VERY HIGH VALUE, METHOD MAY BE INVALID' 1 ,/,' FOR THESE PARAMETERS') GO TO 300 987 CONTINUE C CLOSE(5,STATUS='KEEP') C CLOSE(6,STATUS='KEEP') STOP END