#!/usr/bin/python

import numpy
import matplotlib
import math

fig_size = (3.4,3.4*(numpy.sqrt(5)-1.0)/2.0)
matplotlib.rcParams.update(
        {'legend.borderpad':0.3,
         'legend.labelspacing':0.3,
         'legend.fontsize':'medium',
         'text.usetex':True,
         'axes.labelsize': 8,
         'font.size': 8,
         'legend.fontsize': 8,
         'xtick.labelsize': 8,
         'ytick.labelsize': 8,
         'figure.figsize': fig_size}
)

import pylab

theta, u_xc, u_xc_err = numpy.genfromtxt("dmqmc/rs_1_N_14.dat", usecols=(0,5,6), unpack=True)
x_i = numpy.genfromtxt("madelung.dat", usecols=(2), unpack=True)
E_hf = 13.60355734
r_s = 1.0
E_M = x_i[1] * 0.5 * 14.0

T_0 = 6.0*math.pow((((2.0*math.pi*math.pow(3.0,(1.0/3.0)))/(r_s*math.pow((4.0*math.pi*14.0),(1.0/3.0))))),2.0)
print "T_0", T_0
print "E_M", E_M
# Data for E_c and E_c_err was published in
# V. A. Neufeld, A. J. W. Thom (2017), "A study of the dense uniform electron gas with high orders of coupled cluster", JCP, 147, 194105. 
# (Data repository: Neufeld, V., & Thom, A. J. (2018). Research data and further information supporting "A study of the dense uniform electron gas with high orders of coupled cluster" [Dataset]. https://doi.org/10.17863/CAM.14336).
E_c = numpy.asarray([-0.514495863928, -0.530688601627, -0.531082596101, -0.531126436765, -0.5313])
E_c_err = numpy.asarray([8.80695138074e-05, 0.000167256507004, 0.000183888618594, 0.000187439694048, 0.0002])
E_xc = (E_c + E_hf + E_M - T_0)/14.0
E_xc_err = (E_c_err)/14.0
print "E_xc", E_xc
print "E_xc_err", E_xc_err

f = pylab.gcf()
ax = pylab.gca()

pylab.axhline(E_xc[0], color='g', linewidth=0.4, label=r'$\textrm{CCSD}$')
pylab.axhline(E_xc[1], color='r', linewidth=0.4, label=r'$\textrm{CCSDT}$')
#pylab.axhline(E_xc[2], color='c', linewidth=0.4, label=r'$\textrm{CCSDTQ}$')
#pylab.axhline(E_xc[3], color='m', linewidth=0.4, label=r'$\textrm{CCSDTQ5}$')
pylab.axhline(E_xc[4], color='k', linewidth=0.4, label=r'$\textrm{FCIQMC}$')

(_, caps, _) = pylab.errorbar(theta, u_xc, yerr=u_xc_err, color='blue', marker='x', markerfacecolor='blue', markeredgecolor='blue', markersize=3.0, mew=0.5, linewidth=0.5, capsize=1.5, label=r'$\textrm{DMQMC}$')
for cap in caps:
    cap.set_color('blue')
    cap.set_markeredgewidth(0.6)

ax.set_xlim(0.1, 10)
ax.set_ylim(-0.58, -0.42)
ax.set_xscale('log')
ax.set_xlabel(r'$\Theta$')
ax.set_ylabel(r'$E_{\textrm{xc}}/E_{\textrm{h}}$')
ax.legend(loc='lower right', numpoints=1, fontsize=7)

# inset zoomed in
axi = pylab.axes((0.32,0.62,0.3,0.3))
axi.axhline(E_xc[0], color='g', linewidth=0.4, label=r'$\textrm{CCSD}$')
axi.axhline(E_xc[1], color='r', linewidth=0.4, label=r'$\textrm{CCSDT}$')
#axi.axhline(E_xc[2], color='c', linewidth=0.4, label=r'$\textrm{CCSDTQ}$')
#axi.axhline(E_xc[3], color='m', linewidth=0.4, label=r'$\textrm{CCSDTQ5}$')
axi.axhline(E_xc[4], color='k', linewidth=0.4, label=r'$\textrm{FCIQMC}$')

(_, caps, _) = axi.errorbar(theta, u_xc, yerr=u_xc_err, color='blue', marker='x', markerfacecolor='blue', markeredgecolor='blue', markersize=3.0, mew=0.5, linewidth=0.5, capsize=1.5, label=r'$\textrm{DMQMC}$')
for cap in caps:
    cap.set_color('blue')
    cap.set_markeredgewidth(0.6)

axi.set_xlim(0.12,0.13)
axi.set_ylim(-0.558,-0.55)
axi.locator_params(nbins=3)

f.subplots_adjust(left=0.18, bottom=0.18, right=0.96, top=0.97, hspace=0)

pylab.savefig('rs1.eps', dpi=1200, bbox_inches='tight')
