#!/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;