% Plot nonequilibrium free energy.

clear;

timesteps = {0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00}; % timesteps to plot
timesteps = {0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.50, 3.00, 3.50, 4.0, 5.0, 6.0, 7.0}; % timesteps to plot
timesteps = {0.50, 1.00, 1.50, 2.00, 3.00, 4.0, 5.0, 6.0, 7.0}; % timesteps to plot
waters = {220, 480, 734, 1094, 1288, 1575}; % number of waters to plot

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 = {3520, 1760, 880, 440, 220}; % number of waters to plot (reverse order)

ntimesteps = length(timesteps);
colors = jet(ntimesteps);
nwaters = length(waters);
nsteps = [1 2 4 8 16 32 64 128 256 512 1024 2048];

constraints = 'constrained';
%constraints = 'flexible';

for water_index = 1:nwaters
	  nwat = waters{water_index};
          prefix = sprintf('TIP3P-%d-%s', nwat, constraints);          

clf;
hold on;
for i = 1:ntimesteps
  try
    timestep = timesteps{i};
    filename = sprintf('analysis3/%s-timestep-%04.2f.out', prefix, timestep)
    x = textread(filename);
    color = colors(i);
    h = errorbar(x(:,1), x(:,9), 2*x(:,10), '.-')
    set(h, 'color', colors(i,:));	 
  catch
    % Stability limit reached.
    ntimesteps = i-1;
    break;
  end
end

timesteps_labels = cell(1,ntimesteps);
for i = 1:ntimesteps
	  timestep = timesteps{i};
          timesteps_labels{i} = sprintf('%4.2f fs', timestep);
end

h = legend(timesteps_labels);
set(h, 'location', 'northwest');

xlabel('number of timesteps');
ylabel('nonequilibrium free energy / kT');
axis([-10 2100 0 10]);
plot([-10 2100], [0 0], 'k-');
set(gca, 'xtick', nsteps);
title(prefix);

filename = sprintf('%s-convergence.eps', prefix);
exportfig(gcf, filename, 'width', 8, 'height', 6, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));

end
