% Results from preliminary run.

clear;

% free energy of binding, in kcal/mol
yank.hydration.DG = [
    -2.563 0.003 % 1-methylpyrrole    
    -0.438 0.003 % 12-dichlorobenzene
%    -7.296 0.007 % 2-fluorobenzaldehyde
%    -3.944 0.004 % 23-benzofuran
    -1.269 0.000 % benzene
%    -1.917 0.017 % benzenedithiol
    -0.234 0.002 % ethylbenzene
%    -1.982 0.002% indene
%    -5.590 0.004 % indole
     0.319 0.003 % isobutylbenzene
    +0.414 0.002 % n-butylbenzene
    -4.474 0.009 % n-methylaniline
     0.089 0.003% n-propylbenzene
    -0.243 0.002 % o-xylene
    -0.056 0.002 % p-xylene
    -5.912 0.014 % phenol
%    -5.199 0.007%    % thieno_23c_pyridine
    -0.670 0.001 % toluene    
     ];

% enthalpy of hydration, in kcal/mol
yank.hydration.DH = [
    -2.620 0.168 % 1-methylpyrrole
    -0.467 0.194 % 12-dichlorobenzene
%    -7.593 0.196 % 2-fluorobenzaldehyde
%    -4.097 0.201 % 23-benzofuran
    -1.231 0.166 % benzene
%   -1.534 0.195 % benzenedithiol
    -0.080 0.209 % ethylbenzene    
%    -2.150 0.209 % indene
%    -5.506 0.235 % indole
    0.349 0.251% isobutylbenzene
    0.508 0.257 % n-butylbenzene
   -4.584 0.197 % n-methylaniline
    0.097 0.238 % n-propylbenzene
   -0.361 0.210 % o-xylene   
   -0.004 0.194 % p-xylene
   -5.834 0.201 % phenol
%  -5.109 0.190   % thieno_23c_pyridine
   -0.612 0.185 % toluene
    ];


% From David Mobley's supplementary material from the explicit solvent paper
experiment.hydration.DG = [
    -2.89 0.2 % 1-methylpyrrole
    -1.36 0.2 % 12-dichlorobenzene
	      % 2-fluorobenzaldehyde
              % 23-benzofuran
    -0.86 0.2 % benzene
              % benzenedithiol
    -0.79 0.2 % ethylbenzene
              % indene
              % indole
     0.16 0.2 % isobutylbenzene
    -0.40 0.2 % n-butylbenzene
    -4.69 0.2 % N-methylaniline
    -0.53 0.2 % n-propylbenzene
    -0.90 0.2 % o-xylene
    -0.80 0.2 % p-xylene    
    -6.61 0.2 % phenol
              % thieno_23c_pyridine
    -0.89 0.2 % toluene
		   ];

% David Mobley's numbers
explicit.hydration.DG = [
    -2.59 0.02 % 1-methylpyrrole
    -0.64 0.02 % 12-dichlorobenzene
	       % 2-fluorobenzaldehyde
               % 23-benzofuran
    -0.70 0.02 % benzene    
               % benzenedithiol    
    -0.59 0.02 % ethylbenzene
               % indene
               % indole
     0.00 0.02 % isobutylbenzene
    -0.30 0.02 % n-butylbenzene
    -5.74 0.03 % N-methylaniline
     0.01 0.02 % n-propylbenzene
    -0.52 0.02 % o-xylene    
    -0.67 0.02 % p-xylene    
    -5.67 0.03 % phenol
	       % thieno_23c_pyridine
    -0.71 0.02 % toluene
		   ];

molecules.hydration = {
    '1-methylpyrrole',
    '1,2-dichlorobenzene',
%    '2-fluorobenzaldehyde',
%    '2,3-benzofuran',
    'benzene',
%    'benzenedithiol',
    'ethylbenzene',
%    'indene',
%    'indole'
    'isobutylbenzene',
    'n-butylbenzene',
    'N-methylaniline',
    'n-propylbenzene',
    'o-xylene',
    'p-xylene',
    'phenol',
%    'thieno-2,3c-pyridine',
    'toluene'
    };
		    
 yank.binding.DG = [
     -3.231  0.092 % 1-methylpyrrole
     -9.136  0.328 % 12-dichlorobenzene
     -5.920  0.314 % 2-fluorobenzaldehyde
     -5.285  0.203 % 23-benzofuran
%     -6.335  0.080 % benzene
%     -8.106  0.078 % benzenedithiol
     -9.102  0.083 % ethylbenzene
     -9.428  0.081 % indene
     -7.107  0.055 % indole
     -9.932  0.076 % isobutylbenzene
    -11.979  0.204 % n-butylbenzene
     -6.559  0.022 % n-methylaniline
     -6.758  0.172 % n-propylbenzene
     -8.039  0.085 % o-xylene
     -7.068  0.147 % p-xylene
     -5.532  0.012 % phenol
%     -8.015  0.017 % thieno_23c_pyridine
     -7.307  0.052 % toluene
     ];
 yank.binding.DG(:,1) = yank.binding.DG(:,1) + 0.6 * log(8 * pi^2); % neglected correction
 
 yank.binding.DH = [
     -1.557  2.971 % 1-methylpyrrole
    -10.750  7.409 % 12-dichlorobenzene
     -8.428  8.290 % 2-fluorobenzaldehyde
     -1.695  2.526 % 23-benzofuran
%     -8.500 2.237 % benzene
%     21.663 1.819 % benzenedithiol
     -10.726 1.911 % ethylbenzene
      -9.119 1.690 % indene
      -6.845 1.300 % indole
     -19.078 1.469 % isobutylbenzene
     -25.702 4.082 % n-butylbenzene
     -11.050 2.573 % n-methylaniline
     -26.368 1.630 % n-propylbenzene
     -19.166 1.939 % o-xylene
      -8.155 2.742 % p-xyene
     -17.353 2.154 % phenol
%     -11.442 1.491 % thieno_23c_pyridine
     -11.170 1.293 % toluene
		   ];

 experiment.binding.DG = [
             -4.44 0.2  % 1-methylpyrrole
             -6.37 0.2  % 12-dichlorobenzene
              0.00 2.74  % 2-fluorobenzaldehyde (NB)
             -5.46 0.03  % 23-benzofuran
%             -5.19 0.16  % benzene
%                        % benzenedithiol (ND)
             -5.76 0.07 % ethylbenzene
             -5.13 0.01 % indene
             -4.89 0.06 % indole
             -6.51 0.06 % isobutylbenzene
             -6.70 0.02 % n-butylbenzene
             -4.7  0.2 % n-methylaniline
             -6.55 0.02 % n-propylbenzene
             -4.60 0.06 % o-xylene
             -4.67 0.06 % p-xylene
              0.00 2.74 % phenol (NB)
         % thieno_23c_pyridine (NB)
             -5.52 0.06 % toluene
]; 

 % David Mobley's explicit solvent results
 explicit.binding.DG = [
             -4.32 0.08 % 1-methylpyrrole
             -5.66 0.15 % 12-dichlorobenzene
             -3.14 0.13 % 2-fluorobenzaldehyde 
             -3.66 0.06 % 23-benzofuran
%             -3.95 0.20 % benzene
%             -2.79 0.13 % benzenedithiol       
             -5.82 0.14 % ethylbenzene
             -1.63 0.07 % indene
             -1.37 0.10 % indole
             -8.09 0.18 % isobutylbenzene
             -5.70 0.20 % n-butylbenzene
             -5.37 0.11 % n-methylaniline
             -5.44 0.10 % n-propylbenzene
             -3.23 0.25 % o-xylene
             -3.59 0.12 % p-xylene
             -1.07 0.20 % phenol
             % -2.56 0.07 % thieno_23c_pyridine
             -4.07 0.10 % toluene
     ];

 molecules.binding = {
     '1-methylpyrrole',
     '1,2-dichlorobenzene',
     '2-fluorobenzaldehyde',
     '2,3-benzofuran',
%     'benzene',
     'ethylbenzene',
     'indene',
     'indole',
     'isobutylbenzene',
     'n-butylbenzene',
     'N-methylaniline',
     'n-propylbenzene',
     'o-xylene',
     'p-xylene',
     'phenol',
     'toluene'
     }; 
 
% 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
print -depsc yank-vs-experiment-lysozyme.eps
system('epstopdf yank-vs-experiment-lysozyme.eps');

% 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-lysozyme.eps
system('epstopdf yank-vs-explicit-lysozyme.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-lysozyme.eps
system('epstopdf dg-vs-dh-lysozyme.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-lysozyme.eps
system('epstopdf explicit-vs-experiment-lysozyme.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-lysozyme.eps
system('epstopdf yank-hydration-vs-experiment-lysozyme.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-lysozyme.eps
system('epstopdf yank-hydration-DH-vs-experiment-lysozyme.eps');


