clear;

datadir = 'pyopenmm/processed-data/'; % source directory for analyzed data (.txt files)

xmin = 2.5; 
xmax = 9;
maxits = 10000;
nbins = 50;
%bins = linspace(xmin, xmax, nbins);
bins = xmin + (xmax - xmin)/nbins * ((1:nbins) - 0.5);
barcolor = 'b';

resolution = 450; % dpi for figures
scale = 1.4; % figure scaling

npoints = 1000;
rpoints = linspace(xmin, xmax, npoints);
sigma    = 3.4;
r_WCA    = 2.^(1./6.) * sigma;
h = 5.0;
r0 = r_WCA;
w = 0.5 * r_WCA;
upoints = h*(1-((rpoints-r0-w)/w).^2).^2;
Ppoints = exp(-upoints) .* rpoints.^2;
Ppoints = Ppoints / max(Ppoints);

%% UMBRELLA SAMPLING

filename = sprintf('%s/umbrella-solvent.txt', datadir);
x = textread(filename); 
r = x(:,1);
tau = integrated_autocorrelation_time(r);
disp(sprintf('integrated autocorrelation time of umbrella trajectory: %.1f', tau));
du = x(:,2);
dr = (xmax - xmin) / nbins * 1.0001;
bin_index = floor((r - xmin) / dr) + 1;
log_Pr = zeros(1,nbins);
for i = 1:nbins
  indices = find(bin_index == i);
  if length(indices) == 0
    log_Pr(i) = min(-du);
  else
    log_Pr(i) = logsum(-du(indices));
  end
end
log_Pr(:) = log_Pr(:) - max(log_Pr(:));
Pr = exp(log_Pr);

figure;

subplot('position', [0.07 0.20 0.73 0.67]); 
tau = integrated_autocorrelation_time(r);
plot(r/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0]);
set(gca,'YTick',[1 2]);
xlabel('MC iteration');
ylabel('r / r_0');
%title(sprintf('umbrella sampling (\\tau = %.1f)', tau));
subplot('position', [0.8 0.20 0.17 0.67]); 
[N,X] = hist(r, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Pr * max(N)/max(Pr), bins, 'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

% Generate PDF.
filename = 'wca-dimer-umbrella.eps';
%print('-depsc', filename);
%set(gcf, 'renderer', renderer); 
%exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*3.5, 'resolution', resolution);
%system(sprintf('epstopdf %s', filename));
% Generate separate text and figure
addpath matlab/exportfig
addpath matlab/epscombine
exportfig(gcf, 'zbuffer', 'color', 'cmyk', 'width', scale*3.5, 'height', scale*1.0, 'separatetext', 1, 'renderer', 'zbuffer', 'resolution', resolution);
exportfig(gcf, 'painter', 'color', 'cmyk', 'width', scale*3.5, 'height', scale*1.0, 'separatetext', 1, 'renderer', 'painters');
epscombine('painter_t.eps','zbuffer.eps',filename);
system(sprintf('epstopdf %s', filename));

%% MAKE TRAJECTORY PLOT
figure;

subplot('position', [0.07 0.83 0.73 0.13]); 
filename = sprintf('%s/md-vacuum.txt', datadir);
x = textread(filename); 
tau = integrated_autocorrelation_time(x);
plot(x/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0]); 
set(gca,'YTick',[1 2]);
set(gca,'XTick',[]);
ylabel('r / r_0');
title(sprintf('vacuum MD (\\tau = %.1f)', tau));

subplot('position', [0.8 0.83 0.17 0.13]); 
[N,X] = hist(x, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Ppoints * max(N),rpoints,'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

subplot('position', [0.07 0.64 0.73 0.13]); 
filename = sprintf('%s/mc-vacuum.txt', datadir);
x = textread(filename); 
tau = integrated_autocorrelation_time(x);
plot(x/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0])
set(gca,'YTick',[1 2]);
set(gca,'XTick',[]);
ylabel('r / r_0');
title(sprintf('vacuum MD + MC (\\tau = %.1f)', tau));

subplot('position', [0.8 0.64 0.17 0.13]); 
[N,X] = hist(x, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Ppoints * max(N),rpoints,'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

subplot('position', [0.07 0.45 0.73 0.13]); 
filename = sprintf('%s/md-solvent.txt', datadir);
x = textread(filename); 
tau = integrated_autocorrelation_time(x);
plot(x/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0])
set(gca,'YTick',[1 2]);
set(gca,'XTick',[]);
ylabel('r / r_0');
title(sprintf('solvent MD (\\tau = %.1f)', tau));

subplot('position', [0.8 0.45 0.17 0.13]); 
[N,X] = hist(x, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Pr * max(N), bins,'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

subplot('position', [0.07 0.26 0.73 0.13]); 
filename = sprintf('%s/mc-solvent.txt', datadir);
x = textread(filename); 
tau = integrated_autocorrelation_time(x);
plot(x/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0])
set(gca,'YTick',[1 2]);
set(gca,'XTick',[]);
ylabel('r / r_0');
title(sprintf('solvent MD + MC (\\tau = %.1f)', tau));

subplot('position', [0.8 0.26 0.17 0.13]); 
[N,X] = hist(x, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Pr * max(N), bins,'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

subplot('position', [0.07 0.07 0.73 0.13]); 
filename = sprintf('%s/ncmc-solvent.txt', datadir);
x = textread(filename); 
tau = integrated_autocorrelation_time(x);
plot(x/r0,'k.'); 
hold on; 
plot([0 maxits], 5.5/r0 * [1 1], 'k:'); 
axis([0 maxits xmin/r0 xmax/r0])
set(gca,'YTick',[1 2]);
xlabel('MC iteration');
ylabel('r / r_0');
title(sprintf('solvent MD + 2048-step NCMC (\\tau = %.1f)', tau));

subplot('position', [0.8 0.07 0.17 0.13]); 
[N,X] = hist(x, bins);
h = plot(N,X,'k-');
hold on;
h = plot(Pr * max(N), bins,'r-');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
axis([0 max(N) xmin xmax]);
hold on; 
plot([1 max(N)], 5.5/r0 * [1 1], 'k:'); 

% Generate PDF.
filename = 'wca-dimer-trajectory.eps';
%print('-depsc', filename);
%exportfig(gcf, filename, 'color', 'cmyk', 'width', scale*3.5, 'height', scale*3.5, 'resolution', resolution);
%system(sprintf('epstopdf %s', filename));
% Generate separate text and figure
addpath matlab/exportfig
addpath matlab/epscombine
exportfig(gcf, 'zbuffer', 'color', 'cmyk', 'width', scale*3.5, 'height', scale*3.5, 'separatetext', 1, 'renderer', 'zbuffer', 'resolution', resolution);
exportfig(gcf, 'painter', 'color', 'cmyk', 'width', scale*3.5, 'height', scale*3.5, 'separatetext', 1, 'renderer', 'painters');
epscombine('painter_t.eps','zbuffer.eps',filename);
system(sprintf('epstopdf %s', filename));
