import os
import sys
import numpy as np
import matplotlib.pyplot as plt

E_vacuum, E_reqm = [], []
E_vacuum_06, E_vacuum_07, E_vacuum_08, E_vacuum_09= [], [], [], []
E_reqm_06, E_reqm_07, E_reqm_08, E_reqm_09= [], [], [], []
d1_vacuum, d1_reqm = [], []
d1_vacuum_06, d1_vacuum_07, d1_vacuum_08, d1_vacuum_09= [], [], [], []
d1_reqm_06, d1_reqm_07, d1_reqm_08, d1_reqm_09= [], [], [], []
d2_vacuum, d2_reqm = [], []
d2_vacuum_06, d2_vacuum_07, d2_vacuum_08, d2_vacuum_09= [], [], [], []
d2_reqm_06, d2_reqm_07, d2_reqm_08, d2_reqm_09= [], [], [], []
angle_vacuum, angle_reqm = [], []
angle_vacuum_06, angle_vacuum_07, angle_vacuum_08, angle_vacuum_09= [], [], [], []
angle_reqm_06, angle_reqm_07, angle_reqm_08, angle_reqm_09= [], [], [], []

f_read = open('Final_AuCO_Summary.csv', 'r')
for i in range(11):
    f_read.readline()
for j in range(1385):
    line = f_read.readline()
    index = line.split()[0]
    coordination = int(line.split()[1])
    Ead_vacuum = float(line.split()[2])
    Ead_reqm = float(line.split()[3])
    AuC_vacuum = float(line.split()[4])
    AuC_reqm = float(line.split()[5])
    com_vacuum = float(line.split()[6])
    com_reqm = float(line.split()[7])
    AuCO_angle_vacuum = float(line.split()[8])
    AuCO_angle_reqm = float(line.split()[9])
    E_vacuum.append(Ead_vacuum)
    E_reqm.append(Ead_reqm)
    d1_vacuum.append(AuC_vacuum)
    d1_reqm.append(AuC_reqm)
    d2_vacuum.append(com_vacuum)
    d2_reqm.append(com_reqm)
    angle_vacuum.append(AuCO_angle_vacuum)
    angle_reqm.append(AuCO_angle_reqm)
    if coordination <= 6:
        E_vacuum_06.append(Ead_vacuum)
        E_reqm_06.append(Ead_reqm)
        d1_vacuum_06.append(AuC_vacuum)
        d1_reqm_06.append(AuC_reqm)
        d2_vacuum_06.append(com_vacuum)
        d2_reqm_06.append(com_reqm)
        angle_vacuum_06.append(AuCO_angle_vacuum)
        angle_reqm_06.append(AuCO_angle_reqm)
    if coordination == 7:
        E_vacuum_07.append(Ead_vacuum)
        E_reqm_07.append(Ead_reqm)
        d1_vacuum_07.append(AuC_vacuum)
        d1_reqm_07.append(AuC_reqm)
        d2_vacuum_07.append(com_vacuum)
        d2_reqm_07.append(com_reqm)
        angle_vacuum_07.append(AuCO_angle_vacuum)
        angle_reqm_07.append(AuCO_angle_reqm)
    if coordination == 8:
        E_vacuum_08.append(Ead_vacuum)
        E_reqm_08.append(Ead_reqm)
        d1_vacuum_08.append(AuC_vacuum)
        d1_reqm_08.append(AuC_reqm)
        d2_vacuum_08.append(com_vacuum)
        d2_reqm_08.append(com_reqm)
        angle_vacuum_08.append(AuCO_angle_vacuum)
        angle_reqm_08.append(AuCO_angle_reqm)
    if coordination >= 9:
        E_vacuum_09.append(Ead_vacuum)
        E_reqm_09.append(Ead_reqm)
        d1_vacuum_09.append(AuC_vacuum)
        d1_reqm_09.append(AuC_reqm)
        d2_vacuum_09.append(com_vacuum)
        d2_reqm_09.append(com_reqm)
        angle_vacuum_09.append(AuCO_angle_vacuum)
        angle_reqm_09.append(AuCO_angle_reqm)
f_read.close()    

rmse_E = np.mean(np.sqrt((np.array(E_vacuum)-np.array(E_reqm))**2))
rmse_d1 = np.mean(np.sqrt((np.array(d1_vacuum)-np.array(d1_reqm))**2))
rmse_d2 = np.mean(np.sqrt((np.array(d2_vacuum)-np.array(d2_reqm))**2))
rmse_angle = np.mean(np.sqrt((np.array(angle_vacuum)-np.array(angle_reqm))**2))

rmse_E_06 = np.mean(np.sqrt((np.array(E_vacuum_06)-np.array(E_reqm_06))**2))
rmse_d1_06 = np.mean(np.sqrt((np.array(d1_vacuum_06)-np.array(d1_reqm_06))**2))
rmse_d2_06 = np.mean(np.sqrt((np.array(d2_vacuum_06)-np.array(d2_reqm_06))**2))
rmse_angle_06 = np.mean(np.sqrt((np.array(angle_vacuum_06)-np.array(angle_reqm_06))**2))

rmse_E_07 = np.mean(np.sqrt((np.array(E_vacuum_07)-np.array(E_reqm_07))**2))
rmse_d1_07 = np.mean(np.sqrt((np.array(d1_vacuum_07)-np.array(d1_reqm_07))**2))
rmse_d2_07 = np.mean(np.sqrt((np.array(d2_vacuum_07)-np.array(d2_reqm_07))**2))
rmse_angle_07 = np.mean(np.sqrt((np.array(angle_vacuum_07)-np.array(angle_reqm_07))**2))

rmse_E_08 = np.mean(np.sqrt((np.array(E_vacuum_08)-np.array(E_reqm_08))**2))
rmse_d1_08 = np.mean(np.sqrt((np.array(d1_vacuum_08)-np.array(d1_reqm_08))**2))
rmse_d2_08 = np.mean(np.sqrt((np.array(d2_vacuum_08)-np.array(d2_reqm_08))**2))
rmse_angle_08 = np.mean(np.sqrt((np.array(angle_vacuum_08)-np.array(angle_reqm_08))**2))

rmse_E_09 = np.mean(np.sqrt((np.array(E_vacuum_09)-np.array(E_reqm_09))**2))
rmse_d1_09 = np.mean(np.sqrt((np.array(d1_vacuum_09)-np.array(d1_reqm_09))**2))
rmse_d2_09 = np.mean(np.sqrt((np.array(d2_vacuum_09)-np.array(d2_reqm_09))**2))
rmse_angle_09 = np.mean(np.sqrt((np.array(angle_vacuum_09)-np.array(angle_reqm_09))**2))

fig1 = plt.figure()
fig1.set_size_inches(6, 6)
x_range=np.arange(-1.3, -0.1, 0.01)
plt.xlim(-1.4,-0.05)
plt.ylim(-1.4,-0.05)
plt.plot(E_vacuum, E_reqm, 'r.', label='RMSE = {:.3f} eV'.format(rmse_E))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('E_Vacuum/eV')
plt.ylabel('E_ReQM/eV')
plt.title('CO Adsorption Enegyer Comparison ')
plt.legend()
fig1.savefig('CO_AdE_1.pdf')

fig2 = plt.figure()
fig2.set_size_inches(6, 6)
x_range=np.arange(-1.3, -0.1, 0.01)
plt.xlim(-1.4,-0.05)
plt.ylim(-1.4,-0.05)
plt.plot(E_vacuum_06, E_reqm_06, 'r.', label='Coordination = 6, RMSE = {:.3f} eV'.format(rmse_E_06))
plt.plot(E_vacuum_07, E_reqm_07, 'g.', label='Coordination = 7, RMSE = {:.3f} eV'.format(rmse_E_07))
plt.plot(E_vacuum_08, E_reqm_08, 'b.', label='Coordination = 8, RMSE = {:.3f} eV'.format(rmse_E_08))
plt.plot(E_vacuum_09, E_reqm_09, 'm.', label='Coordination = 9, RMSE = {:.3f} eV'.format(rmse_E_09))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('E_Vacuum/eV')
plt.ylabel('E_ReQM/eV')
plt.title('CO Adsorption Enegyer Comparison ')
plt.legend(fontsize=8)
fig2.savefig('CO_AdE_2.pdf')

fig3 = plt.figure()
fig3.set_size_inches(6, 6)
x_range=np.arange(1.9, 2.6, 0.01)
plt.xlim(1.9, 2.6)
plt.ylim(1.9, 2.6)
plt.plot(d1_vacuum, d1_reqm, 'r.', label='RMSE = {:.3f} A'.format(rmse_d1))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('d_Vacuum/Angstram')
plt.ylabel('d_ReQM/Angstram')
plt.title('CO Adsorption Distance (Au-C) Comparison ')
plt.legend()
fig3.savefig('CO_AuC_1.pdf')

fig4 = plt.figure()
fig4.set_size_inches(6, 6)
x_range=np.arange(1.9, 2.6, 0.01)
plt.xlim(1.9, 2.6)
plt.ylim(1.9, 2.6)
plt.plot(d1_vacuum_06, d1_reqm_06, 'r.', label='Coordination = 6, RMSE = {:.3f} A'.format(rmse_d1_06))
plt.plot(d1_vacuum_07, d1_reqm_07, 'g.', label='Coordination = 7, RMSE = {:.3f} A'.format(rmse_d1_07))
plt.plot(d1_vacuum_08, d1_reqm_08, 'b.', label='Coordination = 8, RMSE = {:.3f} A'.format(rmse_d1_08))
plt.plot(d1_vacuum_09, d1_reqm_09, 'm.', label='Coordination = 9, RMSE = {:.3f} A'.format(rmse_d1_09))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('d_Vacuum/Angstram')
plt.ylabel('d_ReQM/Angstram')
plt.title('CO Adsorption Distance (Au-C) Comparison ')
plt.legend(fontsize=8)
fig4.savefig('CO_AuC_2.pdf')

fig5 = plt.figure()
fig5.set_size_inches(6, 6)
x_range=np.arange(2.5, 3.0, 0.01)
plt.xlim(2.5, 3.0)
plt.ylim(2.5, 3.0)
plt.plot(d2_vacuum, d2_reqm, 'r.', label='RMSE = {:.3f} A'.format(rmse_d2))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('d_Vacuum/Angstram')
plt.ylabel('d_ReQM/Angstram')
plt.title('CO Center of Mass Distance to Au Comparison ')
plt.legend()
fig5.savefig('CO_com_dis_1.pdf')

fig6 = plt.figure()
fig6.set_size_inches(6, 6)
x_range=np.arange(2.5, 3.0, 0.01)
plt.xlim(2.5, 3.0)
plt.ylim(2.5, 3.0)
plt.plot(d2_vacuum_06, d2_reqm_06, 'r.', label='Coordination = 6, RMSE = {:.3f} A'.format(rmse_d1_06))
plt.plot(d2_vacuum_07, d2_reqm_07, 'g.', label='Coordination = 7, RMSE = {:.3f} A'.format(rmse_d1_07))
plt.plot(d2_vacuum_08, d2_reqm_08, 'b.', label='Coordination = 8, RMSE = {:.3f} A'.format(rmse_d1_08))
plt.plot(d2_vacuum_09, d2_reqm_09, 'm.', label='Coordination = 9, RMSE = {:.3f} A'.format(rmse_d1_09))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference Line y=x')
plt.xlabel('d_Vacuum/Angstram')
plt.ylabel('d_ReQM/Angstram')
plt.title('CO Center of Mass Distance to Au Comparison  ')
plt.legend(fontsize=8)
fig6.savefig('CO_com_dis_2.pdf')

fig7 = plt.figure()
fig7.set_size_inches(6, 6)
x_range=np.arange(120, 190, 0.01)
plt.xlim(120, 190)
plt.ylim(120, 190)
plt.plot(angle_vacuum, angle_reqm, 'r.', label='RMSE = {:.3f} A'.format(rmse_d2))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference y=x')
plt.xlabel('angle_Vacuum/degree')
plt.ylabel('angle_ReQM/degree')
plt.title('Adsorption Angle Comparison ')
plt.legend()
fig7.savefig('CO_Angle_1.pdf')

fig8 = plt.figure()
fig8.set_size_inches(6, 6)
x_range=np.arange(120, 190, 0.01)
plt.xlim(120, 190)
plt.ylim(120, 190)
plt.plot(angle_vacuum_06, angle_reqm_06, 'r.', label='Coordination = 6, RMSE = {:.3f} A'.format(rmse_d1_06))
plt.plot(angle_vacuum_07, angle_reqm_07, 'g.', label='Coordination = 7, RMSE = {:.3f} A'.format(rmse_d1_07))
plt.plot(angle_vacuum_08, angle_reqm_08, 'b.', label='Coordination = 8, RMSE = {:.3f} A'.format(rmse_d1_08))
plt.plot(angle_vacuum_09, angle_reqm_09, 'm.', label='Coordination = 9, RMSE = {:.3f} A'.format(rmse_d1_09))
plt.plot(x_range, x_range, 'k-', linewidth=2, label='Reference y=x')
plt.xlabel('d_Vacuum/degree')
plt.ylabel('d_ReQM/degree')
plt.title('Adsorption Angle Comparison ')
plt.legend(fontsize=8)
fig8.savefig('CO_Angle_2.pdf')

fig9 = plt.figure()
plt.hist(E_vacuum_06, bins=50, label = 'Coordination = 6')
plt.hist(E_vacuum_07, bins=50, label = 'Coordination = 7')
plt.hist(E_vacuum_08, bins=50, label = 'Coordination = 8')
plt.hist(E_vacuum_09, bins=50, label = 'Coordination = 9')
plt.xlim(-1.4, -0.1)
plt.ylim(0, 60)
plt.xlabel('Adsorption Energy/eV')
plt.ylabel('Counts')
plt.title('CO Adsorption Energy Distribution (Vacuum)')
plt.legend()
fig9.savefig('CO_hist_Vacuum.pdf')

fig10 = plt.figure()
plt.hist(E_reqm_06, bins=50, label = 'Coordination = 6')
plt.hist(E_reqm_07, bins=50, label = 'Coordination = 7')
plt.hist(E_reqm_08, bins=50, label = 'Coordination = 8')
plt.hist(E_reqm_09, bins=50, label = 'Coordination = 9')
plt.xlim(-1.4, -0.1)
plt.ylim(0, 60)
plt.xlabel('Adsorption Energy/eV')
plt.ylabel('Counts')
plt.title('CO Adsorption Adsorption Energy (ReQM)')
plt.legend()
fig10.savefig('CO_hist_reqm.pdf')


