http://server.ccl.net/cca/software/SOURCES/FORTRAN/unimol/brw.for.shtml
CCL brw.for
```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')
789  FORMAT(20A4)
WRITE(6,789) TITLE
DELR = 0.0005
DELR = DELR*1E-10
300  CONTINUE
IF(EPS.LT.0.) GO TO 987
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
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)
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
```
 Modified: Fri Dec 15 17:00:00 1995 GMT Page accessed 5230 times since Sat Apr 17 21:35:23 1999 GMT