clear all
close all
fclose('all');
clc

font=16;
num=40;
x=0:9:360;
init=zeros(1,num+1);
ES1=init;
ET1=init;
SOCS1T1=init;
SOCT2T1=init;

%%
tickx=0:45:360;
%Gas
yyaxis left
for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Gas/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
    
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Gas/T2/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'T2:'));
    SOCT2T1(j)=str2double(C{Index(2)+1});
    
    fclose(fid);
    
end

SOCS1T1(end)=SOCS1T1(1);
SOCT2T1(end)=SOCT2T1(1);

SOCT2T1g=SOCT2T1;
SOCS1T1g=SOCS1T1;

SOCplot(x,SOCS1T1,font,'','-o')
hold on
SOCplot(x,SOCT2T1,font,'Gas','-s')
xlim([40 310])
legend('$<S_{1}|H_{SO}|T_{1}>$','$<T_{2}|H_{SO}|T_{1}>$','interpreter','latex')

figure
%%

%Hexane

for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Hexane/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
    
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/Hexane/T2/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'T2:'));
    SOCT2T1(j)=str2double(C{Index(2)+1});
    
    fclose(fid);
    
end

SOCS1T1(end)=SOCS1T1(1);
SOCT2T1(end)=SOCT2T1(1);

SOCT2T1h=SOCT2T1;
SOCS1T1h=SOCS1T1;

SOCplot(x,SOCS1T1,font,'','-o')
hold on
SOCplot(x,SOCT2T1,font,'Hexane','-s')
xlim([40 310])
legend('$<S_{1}|H_{SO}|T_{1}>$','$<T_{2}|H_{SO}|T_{1}>$','interpreter','latex')


%%

figure
%ACN

for j=1:num
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/ACN/S1/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'S1:'));
    SOCS1T1(j)=str2double(C{Index(2)+1});
    
    filename = sprintf("C:/Users/Gebruiker/Desktop/Internship/plots/SOC-Pyr-DMA-dihedral-data/SOC-Pyr-DMA-dihedral-data/ACN/T2/angle(%d)/slurm.out",j);
    fid=fopen(filename);
    C = textscan(fid, '%s');
    C=C{1};
    Index = find(contains(C,'T2:'));
    SOCT2T1(j)=str2double(C{Index(2)+1});
    
    fclose(fid);
    
end

SOCS1T1(end)=SOCS1T1(1);
SOCT2T1(end)=SOCT2T1(1);

SOCT2T1a=SOCT2T1;
SOCS1T1a=SOCS1T1;

SOCplot(x,SOCS1T1,font,'','-o')
hold on
SOCplot(x,SOCT2T1,font,'ACN','-s')
xlim([40 310])
legend('$<S_{1}|H_{SO}|T_{1}>$','$<T_{2}|H_{SO}|T_{1}>$','interpreter','latex')

%%
figure
SOCplot(x,SOCS1T1g,font,'','-o')
hold on
SOCplot(x,SOCS1T1h,font,'','-s')
SOCplot(x,SOCS1T1a,font,'','-d')
thS1=69.9; % S1
xline(thS1)% S1
thS1=65.7; % S1
xline(thS1)
thS1=71.6;  % S1
xline(thS1)

xlim([40 310])
legend('$Gas$','$NHX$','$ACN$','interpreter','latex')
title('$<S_{1}|H_{SO}|T_{1}>$','interpreter','latex')

xticks(tickx)
% figure
% SOCplot(x,SOCT2T1g,font,'','-o')
% hold on
% SOCplot(x,SOCT2T1h,font,'','-s')
% SOCplot(x,SOCT2T1a,font,'','-d')
% xlim([40 310])
% legend('$Gas$','$Hexane$','$ACN$','interpreter','latex')
% title('$<T_{2}|H_{SO}|T_{1}>$','interpreter','latex')


%% 

E0G =[9.0414    4.0142    1.6597    0.6723    0.2751    0.1053    0.0283    0.0013    0.0067    0.0387    0.0968    0.1878    0.3383    0.6164    1.1432    2.0870    3.4888    4.1963    3.2059 .....
    2.0130    1.1613    0.6240    0.3331    0.1989    0.1448    0.1258    0.1166    0.1077    0.0991    0.0988    0.1317    0.2583    0.6408    1.6857    4.3766   10.7489   22.2016   32.7855 ....
    31.6133   19.3090    9.0414];


E0H =[8.5801    3.7999    1.5747    0.6440    0.2674    0.1037    0.0283    0.0015    0.0060    0.0361    0.0917    0.1804    0.3290    0.6064    1.1322    2.0703    3.4348    4.0591    3.0741....
    1.9246    1.1041    0.5892    0.3128    0.1868    0.1366    0.1191    0.1110    0.1019    0.0925    0.0916    0.1246    0.2540    0.6470    1.7239    4.4962   11.0050   22.4568   32.5843....
    30.5036   18.3237    8.5801];

E0A =[4.8437    2.1615    0.9589    0.4500    0.2156    0.0912    0.0261    0.0016    0.0053    0.0339    0.0889    0.1838    0.3567    0.6734    1.2330    2.0711    2.7650    2.5283    1.7579....
    1.0688    0.5877    0.3048    0.1646    0.1072    0.0903    0.0890    0.0902    0.0874    0.0821    0.0881    0.1438    0.3442    0.9404    2.5604    6.5495   14.7502   26.9509   25.3640....
    19.7474   10.2653    4.8437];

yyaxis right
% CT_plot(x,E0G,font','','-d')
% hold on
% CT_plot(x,E0A,font','','-d')
% hold on
CT_plot(x,E0H,font','','-d')
ylim([0 1])
ylabel('$E_0 \: (eV)$','interpreter','latex')
yline([0.026])
xlim([50 300])
legend('$Gas$','$ACN$','$NHX$','$E0$','$k_B T$' ,'interpreter','latex')


