#!/usr/bin/perl
# Evaluate a vibrational partition function
# C. David Sherrill, Georgia Institute of Technology
# March 14, 2006
# Notes:
#
# This is NOT the usual way to do this, which involves closed formulas
# instead of an explicit sum over states. Our goal was to explicitly
# account for anharmonicity, which is usually neglected in vibrational
# partition functions using the harmonic oscillator model. If you
# want something that strictly follows harmonic oscillator, this is
# not the best (or even a correct) way to do this. Summary: I wrote
# this script for something I needed, not to be useful to you. Sorry!
# I hope, however, it is useful as a PERL example for beginners.
#
# This script figures out the vibrational partition function for each
# mode, using a simple expression for the vibrational energy of mode,
#
# E(n) = (n+1/2) \omega_e - (n+1/2)^2 \omega_e x_e
#
# The script is dumb and will just add up contributions from each
# quantum number n until it reaches some maximum quantum number defined
# for each mode in the array @v. Typically contributions die off quickly
# with increasing n, so this is ok with modest (~5-20) values of n.
#
# *** However, note that the above energy formula will break down for
# large n, so don't go too crazy. Look at the contributions printed
# out by each quantum number, and break it off at some appropriate value.
# Do NOT let n get so large that the vibrational energies start decreasing
# again! ***
#
# The total vibrational partition function is simply the product of the
# partition functions for each mode.
#
# Ok, enough disclaimers. Here's the simple script.
#
# Temperature in K
$T = 298;
# Harmonic vibrational frequencies, in wavenumbers
@harms = (3364, 1708, 522, 522);
# How many levels to go up in each mode
@v = (20, 30, 20, 20);
# Boltzmann constant in cm^-1 K^-1
$Kb = 0.69509;
# Anharmonicities, estimated as 2pct of harmonics
$i = 0;
foreach $harm(@harms)
{
$anharm = $harm * 0.02;
@anharms[$i++] = $anharm;
}
print "Temperature: $T\n";
print "Harmonics: @harms\n";
print "Anharmonicities: @anharms\n";
$i = 0;
foreach $harm(@harms)
{
$anharm = @anharms[$i];
print "\n\n Contribution for harmonic frequency ", $i+1, " = $harm, ";
print "anharmonicity = $anharm\n";
$q[$i] = 0.0;
for ($v = 0; $v < $v[$i]; $v++) {
$E = ($v + 0.5) * $harm - ($v + 0.5) * ($v + 0.5) * $anharm;
print "E (v=$v): $E\n";
$x = exp(-$E / ($Kb * $T));
$q[$i] += $x;
print "contrib to q: $x\n";
print "q[$i] = $q[$i]\n";
}
$i++;
}
$qvib = 1.0;
foreach $x(@q) {
$qvib *= $x;
}
print "\n\nSummary\n";
print "Vibrational partition functions (per mode):\n";
foreach $x(@q) {
printf "%8.5e\n", $x;
}
print "Total vibrational partition function: ";
printf "%8.5e\n", $qvib;