Converts vibrational frequencies from atomic units to cm-1 for diatomics

convert_to_cm-1.py

#!/usr/bin/env python
"""Converts vibrational frequencies from atomic units to cm-1 for diatomics."""

import sys
from math import sqrt, pi

# Atomic masses obtained using
# import periodictable as pt
# for el in pt.elements: 
#   mass[el.symbol] = sorted([ (el[x].abundance,el[x].mass) for x in el.isotopes ])[-1][1] 

mass = {'H': 1.0078250321, 'He': 4.0026032497, 'Li': 7.016004, 'Be': 9.0121821, 'B': 11.0093055, 'C': 12.0, 'N': 14.0030740052, 'O': 15.9949146221, 'F': 18.9984032, 'Ne': 19.9924401759, 'Na': 22.98976967, 'Mg': 23.9850419, 'Al': 26.98153844, 'Si': 27.9769265327, 'P': 30.97376151, 'S': 31.97207069, 'Cl': 34.96885271, 'Ar': 39.962383123, 'K': 38.9637069, 'Ca': 39.9625912, 'Sc': 44.9559102, 'Ti': 47.9479471, 'V': 50.9439637, 'Cr': 51.9405119, 'Mn': 54.9380496, 'Fe': 55.9349421, 'Co': 58.9332002, 'Ni': 57.9353479, 'Cu': 62.9296011, 'Zn': 63.9291466, 'Ga': 68.925581, 'Ge': 73.9211782, 'As': 74.9215964, 'Se': 79.9165218, 'Br': 78.9183376, 'Kr': 83.911507, 'Rb': 84.9117893, 'Sr': 87.9056143, 'Y': 88.9058479, 'Zr': 89.9047037, 'Nb': 92.9063775, 'Mo': 97.9054078, 'Tc': 114.93828, 'Ru': 101.9043495, 'Rh': 102.905504, 'Pd': 105.903483, 'Ag': 106.905093, 'Cd': 113.9033581, 'In': 114.903878, 'Sn': 119.9021966, 'Sb': 120.903818, 'Te': 129.9062228, 'I': 126.904468, 'Xe': 131.9041545, 'Cs': 132.905447, 'Ba': 137.905241, 'La': 138.906348, 'Ce': 139.905434, 'Pr': 140.907648, 'Nd': 141.907719, 'Pm': 162.95352, 'Sm': 151.919728, 'Eu': 152.921226, 'Gd': 157.924101, 'Tb': 158.925343, 'Dy': 163.929171, 'Ho': 164.930319, 'Er': 165.93029, 'Tm': 168.934211, 'Yb': 173.9388581, 'Lu': 174.9407679, 'Hf': 179.9465488, 'Ta': 180.947996, 'W': 183.9509326, 'Re': 186.9557508, 'Os': 191.961479, 'Ir': 192.962924, 'Pt': 194.964774, 'Au': 196.966552, 'Hg': 201.970626, 'Tl': 204.974412, 'Pb': 207.976636, 'Bi': 208.980383, 'Po': 218.0089658, 'At': 223.02534, 'Rn': 228.03808, 'Fr': 232.04965, 'Ra': 234.05055, 'Ac': 236.05518, 'Th': 232.0380504, 'Pa': 231.0358789, 'U': 238.0507826, 'Np': 244.06785, 'Pu': 247.07407, 'Am': 249.07848, 'Cm': 252.08487, 'Bk': 254.0906, 'Cf': 256.09344, 'Es': 257.09598, 'Fm': 259.10059, 'Md': 260.10365, 'No': 262.10752, 'Lr': 263.11139, 'Rf': 264.11398, 'Db': 265.11866, 'Sg': 266.12193, 'Bh': 267.12774, 'Hs': 269.13411, 'Mt': 271.14123, 'Ds': 273.14925, 'Rg': 272.15348, 'Cn': 0, 'Nh': 0, 'Fl': 0, 'Mc': 0, 'Lv': 0, 'Ts': 0, 'Og': 0}

def convert(e1,e2,f):
  # Conversion factors
  hartree = 4.3597447222071e-18  # joules
  bohr    = 1./18897161646.321   # m
  amu     = 1.6605402e-27        # kg
  c       = 299792458.0          # m/s
  mole    = 6.02214076e23

  # Reduced mass in kg
  mu = mass[e1]*mass[e2] / (mass[e1]+mass[e2]) * amu

  # Frequency in reduced coordinates
  lam = (f * hartree / (bohr*bohr) ) / mu

  # Convert to wave numbers
  nu = sqrt(lam)/(2.*pi*c) * 0.01

  return nu

print("Element 1? "),
e1 = sys.stdin.readline().strip()
                                                                                              
print("Element 2? "),
e2 = sys.stdin.readline().strip()

print("Frequency (in hartree/bohr^2) ? "),
f = float(sys.stdin.readline().strip())

print(convert(e1,e2,f)," cm-1")