boltzmann.py

#!/usr/bin/env python3

"""
ABOUT

 Calculation Boltzmann population of conformers at 298K.
 Relative energies (kcal/mol) are given as arguments to the script.

 Output:
 The percent abundance associated with each relative free energy
 is displayed on stdout.

INVOCATION~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Example:
 python3 ./boltzmann.py 0.00 0.05 1.00
 or
 ./boltzmann.py 0.00 0.05 1.00

 Output will be:

           0.00: 47.53%
           0.05: 43.69%
           1.00: 8.78%

AUTHOR ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 Billy Wayne McCann
 email : [email protected]
 license : ItsYours (BSD-like)
"""

from sys import exit, argv
from math import exp, log

if len(argv[1:]) == 0:
    print("Input the deltaG values as arguments to this script.")
    exit(0)

delta_Gs = map(float, argv[1:])

# gas constant times 298K
RT = 0.5921


def exponential(delta_G):
    return exp(-delta_G / RT)


def sum_exponentials(energies):
    return sum(map(exponential, energies))


# partition function
distribution = sum_exponentials(delta_Gs)

print("partition function     Q = {0:.2f}".format(distribution))
print("Relative Abundances")
for delta_G in delta_Gs:
    percent_abundance = exponential(delta_G) / distribution * 100
    print("\t%.2f: %.2f%%" % (delta_G, percent_abundance))