The Kitchin Research Group

import numpy as np
import matplotlib.pyplot as plt
from pycse import nlinfit

# data
V = np.array([row[0] for row in data]) 
E = np.array([row[1] for row in data])

plt.plot(V, E, '.')
plt.xlabel('Volume ($\AA^3$)')
plt.ylabel('Energy (Ha)')

def Murnaghan(vol, E0, B0, BP, V0):
    '''
    given a vector of parameters and volumes, return a vector of energies.
    equation From PRB 28,5480 (1983)
    '''
    
    E = E0 + B0*vol/BP*(((V0/vol)**BP)/(BP-1)+1) - V0*B0/(BP-1.)

    return E

guess = [-400, 0.5, 2, 210]
pars, pint, SE = nlinfit(Murnaghan, V, E, guess, alpha=0.05)
E0, B0, BP, V0 = pint

Vfit = np.linspace(V.min(), V.max())
plt.plot(Vfit, Murnaghan(Vfit, *pars))
plt.savefig('images/eos-uncertainty.png')

print '95% confidence intervals'
print 'V0 = {0} bohr**3'.format(V0)
print 'E0 = {0} Ha'.format(E0)
print 'B0 = {0} GPA'.format([x * 29421.010901602753 for x in B0])