import numpy.polynomial.polynomial as poly
import matplotlib.pyplot as plt

filename = 'e-cut-conv.hdf5'
e_cut_list = range(20,110,10)
energies = []

for i,e_cut in enumerate(e_cut_list): 
    energy = nlread(filename, TotalEnergy)[i].evaluate()
    energies.append(energy.inUnitsOf(eV))
    print(e_cut)
    print(energy)

for i,e_cut in enumerate(e_cut_list):
    print(e_cut)
    print(energies[i]-energies[-1])

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('E$_{cut}$ (Hartree)')
ax.set_ylabel('E$_{tot}$ (eV)')
ax.plot(e_cut_list, energies-energies[-1], linewidth=2,label='Relative total energy of bulk copper')
ax.legend()
plt.savefig('e-cut-conv-7k.png')

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('E$_{cut}$ (Hartree)')
ax.set_ylabel('E$_{tot}$ (meV)')
ax.axhline(y=10, xmin=0, xmax=1, color='k', ls='--', label='10 meV')
ax.axhline(y=1, xmin=0, xmax=1, color='k', ls=':', label='1 meV')
ax.scatter(e_cut_list[3:], (energies[3:]-energies[-1])*1000)
ax.plot(e_cut_list[3:], (energies[3:]-energies[-1])*1000, linewidth=2,label='Relative total energy of bulk copper')
ax.legend()
ax.set_ylim(0,12)
ax.set_xlim(40,100)
plt.savefig('e-cut-conv-7k-zoom.png')