clear
simname = ['NM_1fs_1ns'];
load (['allPotE.' simname '.out'])
load (['allKinE.' simname '.out'])
load ('temperature.set')
temp = num2str(temperature);
numtemps = size(temperature);
numtemps = numtemps(:,1);
allTotE=allPotE+allKinE;
eqsteps = 10

PEmeans=zeros(1,numtemps);
PEstdvs=zeros(1,numtemps);
PEvars=zeros(1,numtemps);
KEmeans=zeros(1,numtemps);
KEvars=zeros(1,numtemps);
TEmeans=zeros(1,numtemps);
TEvars=zeros(1,numtemps);

figure(1)
i=1;
hold all
while i<numtemps + 1
[histVals histCent] = hist(allPotE(eqsteps:end,i));
plot(histCent, histVals)
PEmeans(i)=mean(allPotE(eqsteps:end,i));
PEstdvs(i)=std(allPotE(eqsteps:end,i));
PEvars(i)=var(allPotE(eqsteps:end,i));
i=i+1;
end
title('Potential Energy Dist for REM Replicas','FontSize',14,'FontWeight','Bold')
xlabel('Potential Energy')
ylabel('Counts')
legend(temp,'Location','NorthEastOutside')

print -depsc PotEhist;

figure(2)
i=1;
hold all
while i<numtemps + 1
[histVals histCent] = hist(allKinE(eqsteps:end,i));
plot(histCent, histVals)
KEmeans(i)=mean(allKinE(eqsteps:end,i));
KEvars(i)=var(allKinE(eqsteps:end,i));
i=i+1;
end
title('Kinetic Energy Dist for REM Replicas','FontSize',14,'FontWeight','Bold')
xlabel('Kinetic Energy')
ylabel('Counts')
legend(temp,'Location','NorthEastOutside')

print -depsc KinEhist;

figure(3)
i=1;
hold all
while i<numtemps + 1
[histVals histCent] = hist(allTotE(eqsteps:end,i));
plot(histCent, histVals)
TEmeans(i)=mean(allTotE(eqsteps:end,i));
TEvars(i)=var(allTotE(eqsteps:end,i));
i=i+1;
end
title('Total Energy Dist for REM Replicas','FontSize',14,'FontWeight','Bold')
xlabel('Total Energy')
ylabel('Counts')
legend(temp,'Location','NorthEastOutside')

print -depsc TotEhist;

clear histCent
clear histVals

%save REMdist_10Reps_NM.mat