h = figure(1);
set(h,'Position',[100 678 2000/3 2000/4]);

index = 0;
prefix     = 'Data/TwoPlatesTight_Charged_L25_index';
prefix_PMF = 'Data/TwoPlatesTight_Charged_PMF_L25_index';
fname      = sprintf([prefix, '%1.1d.mat'], index);
fname_PMF  = sprintf([prefix_PMF, '%1.1d.mat'], index);
load(fname);
geom_d = energy(1,:);
VdW_d  = energy(2,:);
elec_d = energy(3,:);
load(fname_PMF);
LJ_pair   = energy(1,:);
elec_pair = energy(2,:);
elec_d_outside = energy(3,:);
geom_infty = energy(4,1);
VdW_infty = energy(4,2);
elec_infty = energy(4,3);
geom  = geom_d - geom_infty;
VdW   = VdW_d - VdW_infty + LJ_pair;
elec  = elec_d + elec_d_outside - elec_infty + elec_pair;
total = geom + VdW + elec;
subplot(2,2,1);  plot(4:1:24,geom,'r-o');hold on;
subplot(2,2,2); plot(4:1:24,VdW,'r-o');hold on;
subplot(2,2,3); plot(4:1:24,elec,'r-o');hold on;
subplot(2,2,4); plot(4:1:24,total,'r-o');hold on;



index = 1;
prefix     = 'Data/TwoPlatesTight_Charged_L25_index';
prefix_PMF = 'Data/TwoPlatesTight_Charged_PMF_L25_index';
fname      = sprintf([prefix, '%1.1d.mat'], index);
fname_PMF  = sprintf([prefix_PMF, '%1.1d.mat'], index);
load(fname);
geom_d = energy(1,:);
VdW_d  = energy(2,:);
elec_d = energy(3,:);
load(fname_PMF);
LJ_pair   = energy(1,:);
elec_pair = energy(2,:);
elec_d_outside = energy(3,:);
geom_infty = energy(4,1);
VdW_infty = energy(4,2);
elec_infty = energy(4,3);
geom  = geom_d - geom_infty;
VdW   = VdW_d - VdW_infty + LJ_pair;
elec  = elec_d + elec_d_outside - elec_infty + elec_pair;
total = geom + VdW + elec;
subplot(2,2,1);  plot(4:1:24,geom,'k-v');hold on;
subplot(2,2,2); plot(4:1:24,VdW,'k-v');hold on;
subplot(2,2,3); plot(4:1:24,elec,'k-v');hold on;
subplot(2,2,4); plot(4:1:24,total,'k-v');hold on;


index = 2;
prefix     = 'Data/TwoPlatesTight_Charged_L25_index';
prefix_PMF = 'Data/TwoPlatesTight_Charged_PMF_L25_index';
fname      = sprintf([prefix, '%1.1d.mat'], index);
fname_PMF  = sprintf([prefix_PMF, '%1.1d.mat'], index);
load(fname);
geom_d = energy(1,:);
VdW_d  = energy(2,:);
elec_d = energy(3,:);
load(fname_PMF);
LJ_pair   = energy(1,:);
elec_pair = energy(2,:);
elec_d_outside = energy(3,:);
geom_infty = energy(4,1);
VdW_infty = energy(4,2);
elec_infty = energy(4,3);
geom  = geom_d - geom_infty;
VdW   = VdW_d - VdW_infty + LJ_pair;
elec  = elec_d + elec_d_outside - elec_infty + elec_pair;
total = geom + VdW + elec;
subplot(2,2,1);  plot(4:1:24,geom,'b-d');hold on;
subplot(2,2,2); plot(4:1:24,VdW,'b-d');hold on;
subplot(2,2,3); plot(4:1:24,elec,'b-d');hold on;
subplot(2,2,4); plot(4:1:24,total,'b-d');hold on;



index = 3;
prefix     = 'Data/TwoPlatesTight_Charged_L25_index';
prefix_PMF = 'Data/TwoPlatesTight_Charged_PMF_L25_index';
fname      = sprintf([prefix, '%1.1d.mat'], index);
fname_PMF  = sprintf([prefix_PMF, '%1.1d.mat'], index);
load(fname);
geom_d = energy(1,:);
VdW_d  = energy(2,:);
elec_d = energy(3,:);
load(fname_PMF);
LJ_pair   = energy(1,:);
elec_pair = energy(2,:);
elec_d_outside = energy(3,:);
geom_infty = energy(4,1);
VdW_infty = energy(4,2);
elec_infty = energy(4,3);
geom  = geom_d - geom_infty;
VdW   = VdW_d - VdW_infty + LJ_pair;
elec  = elec_d + elec_d_outside - elec_infty + elec_pair;
total = geom + VdW + elec;
subplot(2,2,1);  plot(4:1:24,geom,'m-s');hold on;
subplot(2,2,2); plot(4:1:24,VdW,'m-s');hold on;
subplot(2,2,3); plot(4:1:24,elec,'m-s');hold on;
subplot(2,2,4); plot(4:1:24,total,'m-s');hold on;


index = 4;
prefix     = 'Data/TwoPlatesTight_Charged_L25_index';
prefix_PMF = 'Data/TwoPlatesTight_Charged_PMF_L25_index';
fname      = sprintf([prefix, '%1.1d.mat'], index);
fname_PMF  = sprintf([prefix_PMF, '%1.1d.mat'], index);
load(fname);
geom_d = energy(1,:);
VdW_d  = energy(2,:);
elec_d = energy(3,:);
load(fname_PMF);
LJ_pair   = energy(1,:);
elec_pair = energy(2,:);
elec_d_outside = energy(3,:);
geom_infty = energy(4,1);
VdW_infty = energy(4,2);
elec_infty = energy(4,3);
geom  = geom_d - geom_infty;
VdW   = VdW_d - VdW_infty + LJ_pair;
elec  = elec_d + elec_d_outside - elec_infty + elec_pair;
total = geom + VdW + elec;
%set(0,'defaulttextinterpreter','latex');
subplot(2,2,1);  plot(4:1:24,geom,'g-*');hold on;
h_legend = legend('  0.0e/  0.0e','  0.1e/ 0.1e','- 0.1e/ 0.1e','  0.2e/ 0.2e','- 0.2e/ 0.2e');
set(h_legend,'Location','SouthEast','FontSize',12);
xlim([4,24]); xlabel('d (\AA)','fontsize',15,'interpreter', 'latex');
ylim([-300,100]); ylabel('G_{geom}^{PMF}( k_BT )','fontsize',15,'rot',90);
set(gca,'xtick',[4,6,8,10,12,14,16,18,20,22,24],'fontsize',15);



subplot(2,2,2); plot(4:1:24,VdW,'g-*');hold on;
xlim([4,24]); xlabel('d (\AA)','fontsize',15,'interpreter', 'latex');
ylim([-20,80]); ylabel('G_{VdW}^{PMF}( k_BT )','fontsize',15,'rot',90);
h_legend = legend('  0.0e/  0.0e','  0.1e/ 0.1e','- 0.1e/ 0.1e','  0.2e/ 0.2e','- 0.2e/ 0.2e');
set(h_legend,'Location','NorthEast','FontSize',12);
set(gca,'xtick',[4,6,8,10,12,14,16,18,20,22,24],'fontsize',15);



subplot(2,2,3); plot(4:1:24,elec,'g-*');hold on;
h_legend = legend('  0.0e/  0.0e','  0.1e/ 0.1e','- 0.1e/ 0.1e','  0.2e/ 0.2e','- 0.2e/ 0.2e');
set(h_legend,'Location','NorthEast','FontSize',12);
xlim([4,24]); xlabel('d (\AA)','fontsize',15,'interpreter', 'latex');
ylim([-300,600]); ylabel('G_{elec}^{PMF}( k_BT )','fontsize',15,'rot',90);
set(gca,'xtick',[4,6,8,10,12,14,16,18,20,22,24],'fontsize',15);



subplot(2,2,4); plot(4:1:24,total,'g-*');hold on;
h_legend = legend('  0.0e/  0.0e','  0.1e/ 0.1e','- 0.1e/ 0.1e','  0.2e/ 0.2e','- 0.2e/ 0.2e');
set(h_legend,'Location','NorthEast','FontSize',12);
xlim([4,24]); xlabel('d (\AA)','fontsize',15,'interpreter', 'latex');
ylim([-300,600]); ylabel('G_{tot}^{PMF}( k_BT )','fontsize',15,'rot',90);
set(gca,'xtick',[4,6,8,10,12,14,16,18,20,22,24],'fontsize',15);