mopac606_nbo/src/thermo.f

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