import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scm.plams import Units

file_p = "C_diamond_3x3x3_p_coul.xy"
file_0 = "C_diamond_3x3x3_0_coul.xy"
xy_p = np.loadtxt(file_p, usecols=(0, 1), skiprows=1)
xy_0 = np.loadtxt(file_0, usecols=(0, 1), skiprows=1)

to_eV = Units.convert(1, "hartree", "eV")

# Coul
x = xy_p[:, 0]
y_p = xy_p[:, 1]
y_0 = xy_0[:, 1]
y_p0 = (y_0 - y_p) * to_eV  # eV

# Get dV at half away from the max
idV = np.argmax(y_p0) - int(len(y_p0) / 2)
dV_0p = y_p0[idV]

print(f"dV_0/p {dV_0p:6.3f} eV")

fig, ax = plt.subplots(1)
ax.plot(x, y_p0, color="b", lw=2, label="DFT coulomb V_0/p = V_0 - V_p")
ax.plot(x, [0.0] * len(x), color="r", lw=1, label="V = 0.0 eV", ls="--")
ax.plot(x[idV], y_p0[idV], marker="+", ms=10, color="k", label=f"dV_0/p = {dV_0p:6.3f} eV")
ax.legend()

ax.set_xlabel("x (Å)")
ax.set_ylabel("Potential (eV)")
plt.tight_layout()
plt.savefig("dV_0p.png")
plt.show()