""" Molecular Orbitals for CO """
from numpy import *
from numpy.linalg import *

d=1.13
c_s=-24.4
c_p=-11.3
o_s=-35.1
o_p=-13.6

# to calculate a different molecule, the easiest way is just to use existing formalism, but change the numbers. Parameters for various elements can be faoun in the paper of Harrison. For example, the following set of parameters will perform a calculation for N2 moelcule (although in the output the first N atom will be labeled as "C" and the second as "O". 
#d=1.10
#c_s=-24.4
#c_p=-11.3
#o_s=-35.1
#o_p=-13.6


# values from W. A. Harrison, Phys. Rev. B vol 24, page 5835 (1981)
vss=-1.32*27.2*(0.529/d)**2
vsp=1.42*27.2*(0.529/d)**2
vpp=2.22*27.2*(0.529/d)**2
vppp=-0.63*27.2*(0.529/d)**2


basis=( '|C, s>', '|C,pz>', '|O, s>', '|O,pz>', '|C, px>', '|C,py>', '|O,px>', '|O,py>') 
h=array([[   c_s,       0.,      vss,      vsp,        0.,       0.,      0.,       0.],  # |C, s>
         [    0.,      c_p,     -vsp,      vpp,        0.,       0.,      0.,       0.],  # |C,pz>
         [   vss,     -vsp,      o_s,       0.,        0.,       0.,      0.,       0.],  # |O, s>
         [   vsp,      vpp,       0.,      o_p,        0.,       0.,      0.,       0.],  # |O,pz>
         [    0.,       0.,       0.,       0.,       c_p,       0.,    vppp,       0.],  # |C,px>
         [    0.,       0.,       0.,       0.,        0.,      c_p,      0.,     vppp],  # |C,py>
         [    0.,       0.,       0.,       0.,      vppp,       0.,     o_p,       0.],  # |O,px>
         [    0.,       0.,       0.,       0.,        0.,     vppp,      0.,      o_p]]) # |O,py>


# diagonalize h to find energies and wavevectors
e,psi=eigh(h)

# print the results
for level in range(8):
   print 'energy=%6.2f      |psi>=' % (e[level]),
   for orbital in range(8):
      c=psi[orbital][level]
      if (c >= 5.0e-3):
         print '+%4.2f%6s' % (c, basis[orbital]),
      if (c <= -5.0e-3):
         print '-%4.2f%6s' % (abs(c), basis[orbital]), 
   print
   print 'probabilities per cent:  ',
   for orbital in range(8):
      c=psi[orbital][level]
      if (abs(c) >= 5.0e-3):
         print ' %4.0f%6s' % (100*c*c,basis[orbital]), 
   print
   print