% Generate main figure of nonequilibrium free energy as a function of timestep and number of waters.

clear;
clf;

timesteps = {0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}; % timesteps to plot
%waters = {7040, 3520, 1760, 880, 440, 220}; % number of waters to plot (reverse order)
waters = {3520, 1760, 880, 440, 220}; % number of waters to plot (reverse order)
nwaters = length(waters);

%subplot(2,2,1)
subplot('position', [0.1 0.50 0.38 0.4])
plot_Fneq_panel('flexible', 0, timesteps, waters)
axis([0 6.1 0 30]);
ylabel('\Delta F_{neq} / kT');
title('unconstrained TIP3P');
set(gca, 'XTick', []);
%box on
set(gca, 'YTick', [0 10 20 30]);

%subplot(2,2,2)
subplot('position', [0.55 0.50 0.38 0.4])
plot_Fneq_panel('constrained', 0, timesteps, waters)
axis([0 6.1 0 30]);
ylabel('\Delta F_{neq} / kT');
title('constrained TIP3P');
set(gca, 'XTick', []);
%box on
set(gca, 'YTick', [0 10 20 30]);

%subplot(2,2,3)
subplot('position', [0.1 0.05 0.38 0.4])
plot_Fneq_panel('flexible', 1, timesteps, waters)
axis([0 6.1 0 0.1]);
xlabel('timestep / fs');
ylabel('\Delta F_{neq} / N_{wat}kT');
%box on
set(gca, 'YTick', [0 0.025 0.05 0.075 0.1]);

%subplot(2,2,4)
subplot('position', [0.55 0.05 0.38 0.4])
plot_Fneq_panel('constrained', 1, timesteps, waters)
axis([0 6.1 0 0.1]);
xlabel('timestep / fs');
ylabel('\Delta F_{neq} / N_{wat}kT');
%box ons
set(gca, 'YTick', [0 0.025 0.05 0.075 0.1]);

legend_labels = cell(nwaters,1);
for j = 1:nwaters
   legend_labels{j} = sprintf('%d', waters{j});
end
h = legend(legend_labels);
set(h, 'location', 'northwest');

filename = sprintf('Fneq-main-figure.eps');
exportfig(gcf, filename, 'width', 7.5, 'height', 4, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));
