1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
SUBROUTINE THERMO(A,B,C,LINEAR,SYM,WT,VIBS,NVIBS,ESCF)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION VIBS(*)
LOGICAL LINEAR
CHARACTER KEYWRD*241, KOMENT*81, TITLE*81, TMPKEY*241
COMMON /KEYWRD/ KEYWRD
COMMON /TITLES/ KOMENT,TITLE
C
C
C THERMO CALCULATES THE VARIOUS THERMODYNAMIC QUANTITIES FOR A
C SPECIFIED TEMPERATURE GIVEN THE VIBRATIONAL FREQUENCIES, MOMENTS OF
C INERTIA, MOLECULAR WEIGHT AND SYMMETRY NUMBER.
C
C REFERENCE: G.HERZBERG MOLECULAR SPECTRA AND MOLECULAR STRUCTURE
C VOL 2, CHAP. 5
C
C ---- TABLE OF SYMMETRY NUMBERS ----
C
C C1 CI CS 1 D2 D2D D2H 4 C(INF)V 1
C C2 C2V C2H 2 D3 D3D D3H 6 D(INF)H 2
C C3 C3V C3H 3 D4 D4D D4H 8 T TD 12
C C4 C4V C4H 4 D6 D6D D6H 12 OH 24
C C6 C6V C6H 6 S6 3
C
C
C PROGRAM LIMITATIONS: THE EQUATIONS USED ARE APPROPRIATE TO THE
C HIGH TEMPERATURE LIMIT AND WILL BEGIN TO BE INADEQUATE AT TEMPERA-
C TURES BELOW ABOUT 100 K. SECONDLY THIS PROGRAM IS ONLY APPROPRIATE
C IN THE CASE OF MOLECULES IN WHICH THERE IS NO FREE ROTATION
C
C
C
C
*******************************************************************
*
* THE FOLLOWING CONSTANTS ARE NOW DEFINED:
* PI = CIRCUMFERENCE TO DIAMETER OF A CIRCLE
* R = GAS CONSTANT IN CALORIES/MOLE
* H = PLANCK'S CONSTANT IN ERG-SECONDS
* AK = BOLTZMANN CONSTANT IN ERG/DEGREE
* AC = SPEED OF LIGHT IN CM/SEC
*******************************************************************
SAVE PI, R, H, AK, AC
DIMENSION TRANGE(300)
DATA PI /3.14159D0 /
DATA R/1.98726D0/
DATA H/6.626D-27/
DATA AK/1.3807D-16/
DATA AC/2.99776D+10/
*******************************************************************
IT1=200
IT2=400
ISTEP=10
TMPKEY=KEYWRD
I=INDEX(TMPKEY,'THERMO(')
IF(I.NE.0) THEN
C
C ERASE ALL TEXT FROM TMPKEY EXCEPT THERMO DATA
C
TMPKEY(:I)=' '
TMPKEY(INDEX(TMPKEY,')'):)=' '
IT1=READA(TMPKEY,I)
IF(IT1.LT.100) THEN
WRITE(6,'(//10X,''TEMPERATURE RANGE STARTS TOO LOW,'',
1'' LOWER BOUND IS RESET TO 30K'')')
IT1=100
ENDIF
I=INDEX(TMPKEY,',')
IF(I.NE.0) THEN
TMPKEY(I:I)=' '
IT2=READA(TMPKEY,I)
IF(IT2.LT.IT1) THEN
IT2=IT1+200
ISTEP=10
GOTO 10
ENDIF
I=INDEX(TMPKEY,',')
IF(I.NE.0) THEN
TMPKEY(I:I)=' '
ISTEP=READA(TMPKEY,I)
IF(ISTEP.LT.1)ISTEP=1
ELSE
ISTEP=(IT2-IT1)/20
IF(ISTEP.EQ.0)ISTEP=1
IF(ISTEP.GE.2.AND. ISTEP.LT.5)ISTEP=2
IF(ISTEP.GE.5.AND. ISTEP.LT.10)ISTEP=5
IF(ISTEP.GE.10.AND. ISTEP.LT.20)ISTEP=10
IF(ISTEP.GT.20.AND. ISTEP.LT.50)ISTEP=20
IF(ISTEP.GT.50.AND. ISTEP.LT.100)ISTEP=50
IF(ISTEP.GT.100)ISTEP=100
ENDIF
ELSE
IT2=IT1+200
ENDIF
ENDIF
10 CONTINUE
WRITE(6,'(//,A)')TITLE
WRITE(6,'(A)')KOMENT
IF(LINEAR) THEN
WRITE(6,'(//10X,''MOLECULE IS LINEAR'')')
ELSE
WRITE(6,'(//10X,''MOLECULE IS NOT LINEAR'')')
ENDIF
WRITE(6,'(/10X,''THERE ARE'',I3,'' GENUINE VIBRATIONS IN THIS '',
1''SYSTEM'')')NVIBS
WRITE(6,20)
20 FORMAT(10X,'THIS THERMODYNAMICS CALCULATION IS LIMITED TO',/
110X,'MOLECULES WHICH HAVE NO INTERNAL ROTATIONS'//)
WRITE(6,'(//20X,''CALCULATED THERMODYNAMIC PROPERTIES'')')
WRITE(6,'(42X,''*'')')
WRITE(6,'('' TEMP. (K) PARTITION FUNCTION H.O.F.'',
1'' ENTHALPY HEAT CAPACITY ENTROPY'')')
WRITE(6,'( '' KCAL/MOL'',
1'' CAL/MOLE CAL/K/MOL CAL/K/MOL'',/)')
DO 30 I=1,NVIBS
30 VIBS(I)=ABS(VIBS(I))
ILIM=1
DO 40 ITEMP=IT1,IT2,ISTEP
ILIM=ILIM+1
40 TRANGE(ILIM)=ITEMP
TRANGE(1)=298.D0
DO 80 IR=1,ILIM
ITEMP=TRANGE(IR)
T=ITEMP
C *** INITIALISE SOME VARIABLES ***
C1=H*AC/AK/T
QV=1.0D0
HV=0.0D0
E0=0.0D0
CPV=0.0D0
SV1=0.0D0
SV2=0.0D0
C *** CONSTRUCT THE FREQUENCY DEPENDENT PARTS OF PARTITION FUNCTION
DO 50 I=1,NVIBS
WI=VIBS(I)
EWJ=EXP(-WI*C1)
QV=QV/(1-EWJ)
HV=HV+WI*EWJ/(1-EWJ)
E0=E0+WI
CPV=CPV+WI*WI*EWJ/(1-EWJ)/(1-EWJ)
SV1=SV1+LOG(1.0D0-EWJ)
50 SV2=SV2+WI*EWJ/(1-EWJ)
C *** FINISH CALCULATION OF VIBRATIONAL PARTS ***
HV=HV*R*H*AC/AK
E0=E0*1.4295D0
CPV=CPV*R*C1*C1
SV=SV2*R*C1-R*SV1
C *** NOW CALCULATE THE ROTATIONAL PARTS (FIRST LINEAR MOLECULES
IF(.NOT.LINEAR) GOTO 60
QR=1/(C1*A*SYM)
HR=R*T
CPR=R
SR=R*(LOG(T*AK/(H*AC*A*SYM)))+R
GOTO 70
60 QR=SQRT(PI/(A*B*C*C1*C1*C1))/SYM
HR=3.0D0*R*T/2.0D0
CPR=3.0D0*R/2.0D0
SR=0.5D0*R*(3.D0*LOG(T*AK/(H*AC))
1-2.D0*LOG(SYM)+LOG(PI/(A*B*C))+3.D0)
70 CONTINUE
C *** CALCULATE INTERNAL CONTRIBUTIONS ***
QINT=QV*QR
HINT=HV+HR
CPINT=CPV+CPR
SINT=SV+SR
C *** CONSTRUCT TRANSLATION CONTRIBUTIONS ***
QTR=(SQRT(2.D0*PI*WT*T*AK*1.6606D-24)/H)**3
HTR=5.0D0*R*T/2.0D0
CPTR=5.0D0*R/2.0D0
STR=2.2868D0*(5.0D0*LOG10(T)+3.0D0*LOG10(WT))-2.3135D0
C *** CONSTRUCT TOTALS ***
CPTOT=CPTR+CPINT
STOT=STR+SINT
HTOT=HTR+HINT
C *** OUTPUT SECTION ***
IF(IR.EQ.1)THEN
H298=HTOT
ELSE
WRITE(6,'(/,I7,'' VIB.'',G18.4
1 ,13X,3F11.5 )')ITEMP,QV, HV, CPV, SV
WRITE(6,'(7X,'' ROT.'',G13.3
1 ,16X,3F11.3 )') QR, HR, CPR, SR
WRITE(6,'(7X,'' INT.'',G13.3
1 ,16X,3F11.3 )') QINT,HINT,CPINT,SINT
WRITE(6,'(7X,'' TRA.'',G13.3
1 ,16X,3F11.3)')
2 QTR, HTR, CPTR, STR
WRITE(6,'(7X,'' TOT.'',13X,F17.3,F11.4,2F11.4)')
1 ESCF+(HTOT-H298)/1000.D0,HTOT,CPTOT,STOT
ENDIF
80 CONTINUE
WRITE(6,'(/3X,'' * NOTE: HEATS OF FORMATION ARE RELATIVE TO THE'',
1/12X,'' ELEMENTS IN THEIR STANDARD STATE AT 298K'')')
END
|
