Gaussian: recalculating thermochem
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