% 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
max_timestep = 7.0;
ntimesteps = length(timesteps);
colors = jet(ntimesteps);

dt_i = zeros(1, ntimesteps);
Fneq_i = zeros(1, ntimesteps);
dFneq_i = zeros(1, ntimesteps);

clf;
hold on;
for i = 1:ntimesteps
	  timestep = timesteps{i};
	  filename = sprintf('timestep-%04.2f.out', timestep)
	  x = textread(filename);
          dt_i(i) = timestep;
	  Fneq_i(i) = x(end,9);
          dFneq_i(i) = x(end,10);
end

h = errorbar(dt_i, Fneq_i, 2*dFneq_i, 'k.');
set(h, 'markersize', 10);

xlabel('timestep / fs');
ylabel('nonequilibrium free energy / kT');
xmin = 0;
xmax = max_timestep + 0.1;
axis([xmin xmax -1 12]);
plot([xmin xmax], [0 0], 'r-');
title('nonequilibrium free energy vs. timestep');
set(gca, 'XTick', dt_i);

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


exportfig(gcf, 'nonequilibrium-free-energy-vs-timestep.eps', 'width', 8, 'height', 6, 'color', 'cmyk');
system('epstopdf nonequilibrium-free-energy-vs-timestep.eps');
