% 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);

symbols = 'sdv^<>ph';
markersize = 4;

h = plot_Fneq_panel('flexible', 1, timesteps, waters);
set(h, 'linestyle', '-');
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]);
for i = 1:length(h)
  set(h(i), 'marker', symbols(i));
  set(h(i), 'markersize', markersize);
end

h = plot_Fneq_panel('constrained', 1, timesteps, waters);
set(h, 'linestyle', '-.');
axis([0 6.1 0 0.1]);
xlabel('timestep / fs');
ylabel('\Delta F_{neq} / N_{wat}kT');
%box ons
set(gca, 'YTick', [0 0.02 0.04 0.06 0.08 0.1]);
set(gca, 'YTick', [0 0.025 0.05 0.075 0.1]);
for i = 1:length(h)
  set(h(i), 'marker', symbols(i));
  set(h(i), 'markersize', markersize);
end

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

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