Subject: Gaussian: recalculating thermochem

From: Wai-To Chan <[email-masked]>
Date: Sun, 5 Nov 2000 22:01:57 -0500 (EST)

To Jeremey Greenwood,

  Appended to this message is a piece of FORTRAN code (thermo.f) that 
may serve your purpose. It can evaluate thermochemical quantities
including vibrational and rotational entropies as well as total 
enthalpy (H) and free energy (G) at a temperature specified within
the program using quantities from a gaussian 98 harmonic vibrations
output file. I guess the program you already received from other 
would do this but this program may still help you with developing
an automatic procedure.
You only have to specify the name of the g98 file and 
the program would read off the harmonic vibrations from 
the file. 

  To get the program to run follow these steps 

Compile and link  the program 

	  f77 thermo.f 

  I compiled the program on a SGI machine but it should run on
AIX and other major unix stations.

Submit a g98 harmonic vibration file done with default
settings xxx.log to a test run. You should get 
a listing of thermal quantities which duplicates 
those in your g98 output file.
So if you want to repeat your old g98 jobs
redone at different temperature simply repeat the run with 
the setting of temp from 298.15 K changed to the desired value. 

   Obviously writing such a tiny program doesn't require much 
imagination and everyone who finds it useful are welcomed to
download it. But since I lifted the two major subroutines
for parsing the g98 file 
>from external source please retain the comments where
I acknowledged their original source in any future modifications.  

Wai-To Chan

C****************thermal.f*********************
c234567890123456789012345678901234567890123456789012345678901234567890
      program read
C
      integer ncard, nome, num, zero
      integer index, NUMBER, ILEN
      integer i,j, lpst, nflag, nrot
      real*8 dnum, omega(100), vthe(100), therm
      real*8 mass, mmass, II(4), I1, I2, I3, Rconv
      real*8 Ezpe, Esum, temp, Kb, R, sum, autoj
      real*8 hP, pie, Nav, qtran, qrot, qvib
      real*8 Stran, Srot, Svib, Evib, Etrn, Erot
      real*8 utokg, jtoc, molV, sigma, totS, totE, totG
      real*8 one, oneh, half, two, hc, term, totH
      real*8 sumvE, sumvS, theT, ext
c   Size of Gaussian 98 output file is set to 10000 lines
      character*80 line(10000)
      character*14 freq /'Frequencies --'/ 
      character*15 molM /'Molecular mass:'/
      character*26 rotcs /'ROTATIONAL CONSTANTS (GHZ)'/
      character*25 rotc  /'ROTATIONAL CONSTANT (GHZ)'/
      character*26 rots  /'ROTATIONAL SYMMETRY NUMBER'/ 
      character*22 zcorr /'Zero-point correction='/
      character*32 sumof /'Sum of electronic and zero-point'/
      character*4 LEXT /'.log'/
      character*40 RNT1 
      parameter(zero=0,Kb=1.381D-23,jtoc=1.D0/4.184D0)
      parameter(hP=6.626D-34,pie=22.D0/7.D0,Nav=6.022D23)
      parameter(R=8.314D0,molV=24.79D-3,utokg=1.661D-27)
c     parameter(R=8.314D0,molV=1.D-4,utokg=1.661D-27)
c     Change the value of temp set below if you 
c   want to calculate H and G at different temperature  
      parameter(temp=298.15D0,hc=6.626D-34 * 2.998D10)
      parameter(Rconv=16.86D0*1.661D-27*1.D-20,autoj=627.5D0*4.184D3)
      parameter(one=1.D0,oneh=1.5D0,half=0.5D0,two=2.D0)
c  Use .log as extension of your g98 output file as in xxx.log
      call MAKNAME(1,RNT1,ILEN,LEXT)
      if (ILEN .EQ. 0) stop 'usage: a.out file1'
         open (unit=11, file=RNT1, form='formatted')
      write(6,'(a,a/)') 'Gaussian output file: ', RNT1
30    format(a80)
	  ncard=1
      nome=1
80    continue
      read(11,30,end=90) line(ncard)
       ncard=ncard+1
       goto 80
c
90    continue
c
       ncard = ncard - 1
       nrot=1
       nome=1
      do 100 i = 1, ncard 
c
      if (index(line(i), freq) .ne. zero) then
        lpst = 16
      do 50 j = 1, 3
        nflag = NUMBER(line(i), lpst, num, omega(nome))
      if (nflag .gt. 0) goto 100
        nome = nome + 1
50    continue
      endif
c
       if (index(line(i), molM) .ne. zero) then
        lpst=17
        nflag = NUMBER(line(i), lpst, num, mmass)
       endif
       if (index(line(i), zcorr) .ne. zero) then
         lpst=24
         nflag = NUMBER(line(i), lpst, num, Ezpe)
       endif
       if (index(line(i), sumof) .ne. zero) then
        lpst=44
        nflag = NUMBER(line(i), lpst, num, Esum)
	   endif
c
       if (index(line(i), rots) .ne. zero) then
        lpst=28
        nflag = NUMBER(line(i), lpst, num, sigma)
        endif
       if (index(line(i), rotc) .ne. zero) then
        lpst=28
        nrot=1
        nflag = NUMBER(line(i), lpst, num, II(nrot)) 
       endif
       if (index(line(i), rotcs) .ne. zero) then
	   nrot=3
       do 51 j = 1, 3 
       nflag = NUMBER(line(i), lpst, num, II(j))
51     continue
       endif 
c
100   continue
       close(11)
      write(6,41) 'Molecular mass/amu', mmass
      write(6,'(a,3f12.6)') 
     * 'Rotational Constants: ', (II(j), j = 1, nrot)  
      write(6,41) 'Zero-point correction: ', Ezpe
      write(6,41) 'Sum of electronic and zero-point energies: ', Esum
      write(6,*) 'number of vibrations: ', nome-1
41    format(a,g12.6/)
c     Evaluate translational partition function qtran
c234567
c     qtran=dsqrt((2.*pie*mmass*Nav*Kb*temp)/(hP**2))
      mass = mmass * utokg
      therm = hP * dsqrt(1.D0/(2.D0 * pie * mass * Kb*temp))
      qtran = molV/(therm**3)
      Stran = 2.5D0*R + R*dlog(qtran/Nav)
c     Stran = 2.5D0*R + R*dlog(qtran)
      Etrn = oneh*R*temp
      write(6,*) 'Translational     Q          E          S   '
      write(6,'(11x,3g12.4/)') qtran, Etrn*jtoc/1.D3, Stran*jtoc
c     Evaluate rotational partion function qrot
       qrot = 8.D0 * pie**2 * Kb * temp / hP**2
      if (nrot .eq. 1) then
       I1 =  Rconv/(.0334D0*II(1)) 
       qrot = qrot * I1 / sigma
       Erot = R*temp
       Srot = one*R + R*dlog(qrot)
      else 
       I1 = Rconv/(.0334D0*II(1))
       I2 = Rconv/(.0334D0*II(2))
       I3 = Rconv/(.0334D0*II(3))
       qrot = dsqrt(pie) * qrot**oneh * dsqrt(I1*I2*I3) / sigma 
c23456789012345678901234567890123456789012345678901234567890123456789012
       Erot = oneh*R*temp
       Srot = oneh*R + R*dlog(qrot)
      endif
      write(6,*) 'Rotational   Q         E         S       '
      write(6,'(8x,3g12.4/)') qrot, Erot*jtoc/1.D3, Srot*jtoc
c     Evaluate vibrational partition function
          qvib=one
          sumvE = 0.D0
          sumvS = 0.D0
c  omega(j) contains valueso harmonic frequencies read off your g98 output file
       do 42 j = 1, nome-1
        vthe(j) = hc*omega(j)/Kb
        theT = vthe(j)/temp
        ext = dexp(theT)
        sumvE = vthe(j)/two + vthe(j) * (one / (ext - 1) ) + sumvE
        sumvS = (theT/(ext - one) - dlog(one - (one/ext)) ) + sumvS 
c
c       if (j .eq. 1) write(6,*) 'vib1 S: ', R*sumvS/jtoc
c       qvib = qvib * (dexp(-vthe(j)/(two*temp)) /
c    1                (1 - dexp(-vthe(j)/temp)))
c       term = (one/dsqrt(ext))/(one - (one/ext))
        term = one/(one - one/ext)
        qvib = qvib*term 
c234567
42     continue
       Evib = R*sumvE
       Svib = R*sumvS
       write(6,*) 'Vibrational     Q           E          S  '
       write(6,'(13x,3g12.4/)')  qvib, Evib*jtoc/1.D3, Svib*jtoc 
       totS = Stran + Srot + Svib
       totE = Esum + (Etrn + Erot + Evib)/autoj - Ezpe
       totH = totE + R*temp/autoj
       totG = totH - temp*totS/autoj
       write(6,111) 'Total partition function: ', qtran*qvib*qrot
       write(6,111) 'Total entropy:          ', totS*jtoc
       write(6,111) 'Total zero-temp energy: ', Esum
       write(6,111) 'Thermal correction: ', (Etrn + Erot + Evib)/autoj 
     1 - Ezpe 
       write(6,111) 'Total thermal energy:  ', totE
       write(6,111) 'Total enthalpy:        ', totH
       write(6,111) 'Total free energy:     ', totG
111    format (a,g14.8)
       stop
       end
c
C W.-T. Chan notes: This subroutine is taken from the McMaster AIMPAC package 
      SUBROUTINE MAKNAME(I,STRING,L,EXT)
      CHARACTER*(*) STRING,EXT
      INTEGER I,J,N,L
      CALL GETARG(I,STRING)
      J = LEN(STRING)
      DO 10 N = 1,J
        IF(STRING(N:N) .EQ. ' ') THEN
          L = N - 1
          STRING = STRING(1:L)//EXT
          RETURN
        ENDIF
10    CONTINUE
      STOP ' FAILED TO MAKE A FILE NAME '
      END
c
cWai-To Chan note: This routine is taken from the McMaster AIMPAC package.  
C SKK ================================================================== SKK
C
        FUNCTION        NUMBER	(LINE, LPST, NUM, DEC)
C
C CONVERTS A CHARACTER STRING OF NUMBERS INTO ACTUAL NUMBERS EITHER
C INTEGERS OR DECIMAL MAY BE READ.
C NUMBER = 1 IF ALL THE REMAINING CHARACTERS ARE BLANK
C        = 2 IF CHARACTERS ARE NOT RECOGNISED AS A NUMBER, LPST IS RESET
C SKK ================================================================== SKK

        DOUBLE PRECISION DEC, TEN
        CHARACTER*(*)   LINE
        CHARACTER       BLANK, COMMA, DOT, MINUS, L
        CHARACTER       CTEN(0:9)
        DATA    CTEN    /'0','1','2','3','4','5','6','7','8','9'/
        PARAMETER       (BLANK = ' ', COMMA = ',')
        PARAMETER       (DOT   = '.', MINUS = '-')
        INTEGER         ITEN
        PARAMETER       (ITEN = 10, TEN = 10.0D0)
        NUM     = 0
        DEC     = 0.0D0
        NP      = 0
        ND      = 0
        MS      = 0
        NUMBER	= 0
        LPEND   = LEN (LINE)
5       IF (LINE(LPST:LPST) .EQ. BLANK) THEN
                LPST    = LPST + 1
                IF (LPST .GT. LPEND) THEN
                        NUMBER	= 1
                        RETURN
                END IF
                GOTO 5
        END IF
        LBEFOR  = LPST

        DO 1 I  = LBEFOR, LPEND
        LPST    = I
        L       = LINE(I:I)
        IF (L .EQ. BLANK .OR. L .EQ. COMMA) THEN
                GOTO 2
        ELSE IF (L .EQ. MINUS) THEN
                MS      = 1
                GOTO 1
        ELSE IF (L .EQ. DOT) THEN
                NP      = 1
                GOTO 1
        END IF
        DO 3 J  = 0, 9
        IF (L .EQ. CTEN(J)) THEN
                N       = J
                GOTO 4
        END IF
3       CONTINUE
        NUMBER	= 2
        LPST    = LBEFOR
        RETURN

4       IF (NP .EQ. 1) THEN
                ND      = ND + 1
                DEC     = DEC + DFLOAT(N)/TEN**ND
        ELSE
                NUM     = NUM*ITEN + N
        END IF
1       CONTINUE

2       DEC     = DFLOAT(NUM) + DEC
        IF (MS .EQ. 0) RETURN
        DEC     = -DEC
        NUM     = -NUM
        RETURN
        END