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