% Plot nonequilibrium free energy as a function of timestep and number of waters.

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
timesteps = {0.50, 1.00, 1.50, 2.00, 3.00, 4.0, 5.0, 6.0, 7.0}; % timesteps 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 = {220, 480, 734, 1094, 1288, 1575}; % number of waters to plot
waters = {1575, 1288, 1094, 734, 480, 220}; % number of waters to plot (reverse order)

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

max_timestep = 7.0;
ntimesteps = length(timesteps);
nwaters = length(waters);
colors = jet(ntimesteps);

dt_ni = zeros(nwaters, ntimesteps);
Fneq_ni = zeros(nwaters, ntimesteps);
dFneq_ni = zeros(nwaters, ntimesteps);

clf;
hold on;
for i = 1:ntimesteps
  try
    for j = 1:nwaters
      timestep = timesteps{i};
      nwat = waters{j};
      filename = sprintf('analysis2/TIP3P-%d-%s-timestep-%04.2f.out', nwat, constraints, timestep)
      x = textread(filename);
      dt_ni(j,i) = timestep;
      
      %index = 9; % which ntimesteps to use (256)
      %index = 10; % which ntimesteps to use (512)
      index = 11; % which ntimesteps to use (1024)
      %index = size(x,1); % which ntimesteps to use
      Fneq_ni(j,i) = x(index,9);
      dFneq_ni(j,i) = x(index,10);
    end
  catch
    % Could not read file---the timestep stabiity limit was reached.
    disp(sprintf('stability limit reached at i = %d', i));
    disp(filename);
    ntimesteps = i-1
    dt_ni = dt_ni(:, 1:ntimesteps);
    Fneq_ni = Fneq_ni(:, 1:ntimesteps);
    dFneq_ni = dFneq_ni(:, 1:ntimesteps);
    max_timestep = timesteps{ntimesteps};
    break
  end
end

%h = errorbar(dt_ni', Fneq_ni', 2*dFneq_ni', '.-');
h = semilogy(dt_ni', Fneq_ni', '.-');
set(h, 'markersize', 10);

xlabel('timestep / fs');
ylabel('nonequilibrium free energy / kT');
xmin = 0;
xmax = max_timestep + 0.1;
%axis([xmin xmax -1 30]);
%plot([xmin xmax], [0 0], 'r-');
nsteps = 2^(index-1);
title(sprintf('nonequilibrium free energy vs. timestep for %s TIP3P water box (after %d steps to reach steady-state)', constraints, nsteps));
set(gca, 'XTick', dt_ni(1,:));

% Fit a polynomial.
%npoints = 1000;
%N = 5; % order of polynomial
%M = ntimesteps; % number of timesteps to include
%objective = @(a) sum(((Fneq_i(1:M) - a*dt_i(1:M).^N) ./ dFneq_i(1:M)).^2);
%a = fminsearch(objective, 0);
%X = linspace(xmin, xmax, npoints);
%Y = a * X.^N;
%plot(X, Y, 'b-');

%h = legend('simulation', sprintf('%.3f * dt^%d', a, N));
%set(h, 'location', 'northwest');

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

filename = sprintf('nonequilibrium-free-energy-vs-timestep-and-waters-%s.eps', constraints);
exportfig(gcf, filename, 'width', 8, 'height', 6, 'color', 'cmyk');
system(sprintf('epstopdf %s', filename));
