#!/usr/bin/env python
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

# Both are from the top site and the reaction coordinate is taken from the elbow plot script that uses lines with varying angles to determine the MEP
NN = np.array(
[17.59249722,
18.82288479,
20.16782288,
21.63550304,
23.24490764,
25.02567856,
27.01809095,
29.22516413,
31.6842327,
34.45387154,
37.54747069,
41.01383819,
44.80109416,
48.87550096,
53.27145261,
57.94636353,
62.84194763,
67.70922886,
71.84129161,
75.8361983,
79.97619006,
83.96644962,
87.3314466,
89.94352924,
92.01426026,
93.65220404,
94.97388021,
96.04891701,
96.92438963,
97.63918172,
98.22558691,
98.69822923,
99.07601227,
99.37077704,
99.59265796,
99.74929522,
99.84775143,
99.89348763,
99.89044205,
99.8404788,
99.75102805,
99.62321876,
99.45970269,
99.26237019]
)

DFT = np.array(
[20.8987176,
22.21091779,
23.63890034,
25.18266526,
26.85186108,
28.69473045,
30.74021897,
32.98832663,
35.4872961,
38.28537002,
41.39219692,
44.59550913,
47.82776693,
51.4266689,
55.39221504,
59.7437024,
64.43288835,
69.28609932,
73.71477493,
78.29782704,
82.75544824,
86.67275173,
89.7699301,
92.08557748,
93.88985273,
95.30818675,
96.46601044,
97.41156645,
98.17380038,
98.82025194,
99.38951525,
99.82369914,
100.1710462,
100.4605022,
100.6824184,
100.8560919,
100.9718743,
101.0587111,
101.0876567,
101.039414,
100.9911714,
100.9043346,
100.7885522,
100.6438243]
)

mpl.rc("text", usetex=True)
fig = plt.figure(figsize=(3.69,3))
#fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(3.69,3.))

#plt.subplot(2,1,1)

x = np.arange(0, 44, 1)
plt.plot( x, NN, label='HD-NNP' )
plt.plot( x, DFT, label='DFT' )

plt.fill_between(x, DFT - 4.2, DFT + 4.2, facecolor='grey', alpha=0.5, zorder=-1)

plt.plot( [38,38], [-10., 120.], c='k', linestyle='--' )

plt.legend(loc='best', numpoints=1, frameon=True)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(0,43)
plt.xticks([0, 38], ['0.0', '1.0'])
plt.ylim(0, 110.)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel(r'E (kJ/mol)')
plt.xlabel('Reaction coordinate')

plt.tight_layout()
#plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('reactioncoordinate.pdf')
#plt.show()