% Results from 5 ns/replica run for FKBP.

clear;

% From Holt et al., JACS 115:9925, 1993.
experiment.binding.Ki = [
%    1600 130 % L1
    2000 300 % L2
    660 60 % L3
%    186 24 % L4
    110 20 % L5
    12 5 % L6
%    250 30 % L7
    10 1 % L8
    7 2 % L9
%   450 150 % L0 (listed as 300-600)
%    300 40 % L11
    30 10 % L12
%    100 20 % L13
    1 0.4 % L14    
    0.4 0.1 % L20 / FK506
] * 1e-9; % M

experiment.binding.DG = 0.562 * log(experiment.binding.Ki); % kcal/mol at 283 K
experiment.binding.DG(:,2) = 0.562 ./ experiment.binding.Ki(:,1) .* experiment.binding.Ki(:,2); % error bars

yank.binding.DG = [
    -7.880 0.050 % LG2
    -9.275 0.072 % LG3
    -12.586 0.101 % LG5
    -11.025 0.265 % LG6
    -13.438 0.108 % LG8
    -13.597 0.165 % LG9
    -11.748 0.265 % L12
    -13.278 0.099 % L14
    -12.512 0.063 % L20
    ]; % kcal/mol
%yank.binding.DG(:,1) = yank.binding.DG(:,1) + 0.562 * log(8 * pi^2); % neglected correction

 molecules.binding = {
     'L2',
     'L3',
     'L5',
     'L6',
     'L8',
     'L9',
     'L12',
     'L14',
     'L20/FK506'
     };
 
% YANK vs experiment
figure(1);
clf;
plot(experiment.binding.DG(:,1), yank.binding.DG(:,1), 'k.', 'markersize', 8);
xlabel('experiment DG (kcal/mol)');
ylabel('YANK DG (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([experiment.binding.DG(:,1); yank.binding.DG(:,1)])-1;
Fmax = max([experiment.binding.DG(:,1); yank.binding.DG(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(experiment.binding.DG,1)
  DG1 = experiment.binding.DG(i,1);
  dDG1 = experiment.binding.DG(i,2);  
  DG2 = yank.binding.DG(i,1);
  dDG2 = yank.binding.DG(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.binding{i}); % label
end
title('FKBP');
print -depsc yank-vs-experiment-fkbp.eps
system('epstopdf yank-vs-experiment-fkbp.eps');

% STOP
stop

% YANK versus explicit solvent
figure(2);
clf;
plot(explicit.binding.DG(:,1), yank.binding.DG(:,1), 'k.', 'markersize', 8);
xlabel('explicit DG (kcal/mol)');
ylabel('YANK DG (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([explicit.binding.DG(:,1); yank.binding.DG(:,1)])-1;
Fmax = max([explicit.binding.DG(:,1); yank.binding.DG(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(explicit.binding.DG,1)
  DG1 = explicit.binding.DG(i,1);
  dDG1 = explicit.binding.DG(i,2);  
  DG2 = yank.binding.DG(i,1);
  dDG2 = yank.binding.DG(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.binding{i}); % label
end
print -depsc yank-vs-explicit-fkbp.eps
system('epstopdf yank-vs-explicit-fkbp.eps')

% YANK DG vs DH
figure(3);
clf;
plot(yank.binding.DH(:,1), yank.binding.DG(:,1), 'k.', 'markersize', 8);
xlabel('YANK DH (kcal/mol)');
ylabel('YANK DG (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([yank.binding.DH(:,1); yank.binding.DG(:,1)])-1;
Fmax = max([yank.binding.DH(:,1); yank.binding.DG(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(yank.binding.DH,1)
  DG1 = yank.binding.DH(i,1);
  dDG1 = yank.binding.DH(i,2);  
  DG2 = yank.binding.DG(i,1);
  dDG2 = yank.binding.DG(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.binding{i}); % label
end
print -depsc dg-vs-dh-fkbp.eps
system('epstopdf dg-vs-dh-fkbp.eps')

% explicit vs experiment
figure(4);
clf;
plot(experiment.binding.DG(:,1), explicit.binding.DG(:,1), 'k.', 'markersize', 8);
xlabel('experiment DG (kcal/mol)');
ylabel('explicit DG (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([experiment.binding.DG(:,1); explicit.binding.DG(:,1)])-1;
Fmax = max([experiment.binding.DG(:,1); explicit.binding.DG(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(experiment.binding.DG,1)
  DG1 = experiment.binding.DG(i,1);
  dDG1 = experiment.binding.DG(i,2);  
  DG2 = explicit.binding.DG(i,1);
  dDG2 = explicit.binding.DG(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.binding{i}); % label
end
print -depsc explicit-vs-experiment-fkbp.eps
system('epstopdf explicit-vs-experiment-fkbp.eps');

% hydration: yank DG vs experiment
figure(5);
clf;
plot(experiment.hydration.DG(:,1), yank.hydration.DG(:,1), 'k.', 'markersize', 8);
xlabel('experiment hydration DG (kcal/mol)');
ylabel('yank hydration DG (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([experiment.hydration.DG(:,1); yank.hydration.DG(:,1)])-1;
Fmax = max([experiment.hydration.DG(:,1); yank.hydration.DG(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(experiment.hydration.DG,1)
  DG1 = experiment.hydration.DG(i,1);
  dDG1 = experiment.hydration.DG(i,2);  
  DG2 = yank.hydration.DG(i,1);
  dDG2 = yank.hydration.DG(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.hydration{i}); % label
end
print -depsc yank-hydration-vs-experiment-fkbp.eps
system('epstopdf yank-hydration-vs-experiment-fkbp.eps');

% hydration: yank DH vs experiment
figure(6);
clf;
plot(experiment.hydration.DG(:,1), yank.hydration.DH(:,1), 'k.', 'markersize', 8);
xlabel('experiment hydration DG (kcal/mol)');
ylabel('yank hydration DH (kcal/mol)');
axis square;
oldaxis = axis;
Fmin = min([experiment.hydration.DG(:,1); yank.hydration.DH(:,1)])-1;
Fmax = max([experiment.hydration.DG(:,1); yank.hydration.DH(:,1)])+1;
axis([Fmin Fmax Fmin Fmax]);
hold on;
plot([Fmin Fmax], [Fmin Fmax], 'r-'); % x = y line
% error bars
for i = 1:size(experiment.hydration.DG,1)
  DG1 = experiment.hydration.DG(i,1);
  dDG1 = experiment.hydration.DG(i,2);  
  DG2 = yank.hydration.DH(i,1);
  dDG2 = yank.hydration.DH(i,2);    
  plot(DG1 + [-dDG1 +dDG1], DG2 * [1 1], 'k-');
  plot(DG1 * [1 1], DG2 + [-dDG2, +dDG2], 'k-');  
  text(DG1+0.1,DG2+0.1,molecules.hydration{i}); % label
end
print -depsc yank-hydration-DH-vs-experiment-fkbp.eps
system('epstopdf yank-hydration-DH-vs-experiment-fkbp.eps');


