%Import the files
%assume you are in ~/protomol/examples/alanDipeptideSol
clear;
DELIMITER = ' ';
HEADERLINES = 0;
ad = importdata('./alanine.EW.forceD3.xyz0', DELIMITER, HEADERLINES);  % explicit solvent, Hassan05 charges
c  = importdata('../alanDipeptideVac/alanine.EW.forceD.xyz0', DELIMITER, HEADERLINES); % dielectric D = 53.1,S=0.25
c1  = importdata('../alanDipeptideVac/alanine.EW.forceD1.xyz0', DELIMITER, HEADERLINES); % dielectric D=53,S=0.20
c2  = importdata('../alanDipeptideVac/alanine.EW.forceD2.xyz0', DELIMITER, HEADERLINES); % dielectric D=78,S=0.25
  c3  = importdata('../alanDipeptideVac/alanine.EW.forceD3.xyz0', DELIMITER, HEADERLINES); % dielectric D=78,S=0.25, Hassan05 charges
  h1  = importdata('../alanDipeptideVac/alanine.EW.forceH1.xyz0', DELIMITER, HEADERLINES); % scpism D=78, Hassan05 charges
  h2  = importdata('../alanDipeptideVac/alanine.EW.forceH2.xyz0', DELIMITER, HEADERLINES); % scpism D=78, Hassan05 charges
  h3  = importdata('../alanDipeptideVac/alanine.EW.forceH3.xyz0', DELIMITER, HEADERLINES); % unmodified Coulomb, Hassan05 charges
    h4  = importdata('../alanDipeptideVac/alanine.EW.forceH4.xyz0', DELIMITER, HEADERLINES); % D = 53.1, Hassan05 charges (WORST)
m  = ad(:,4); 
ad = ad(:,1:3);
c = c(:,1:3);
c1 = c1(:,1:3);
c2 = c2(:,1:3);
c3 = c3(:,1:3);
h1 = h1(:,1:3);
h2 = h2(:,1:3);
h3 = h3(:,1:3);
h4 = h4(:,1:3);
l = length(m);

% denominator
d=zeros(l,1);
for i=1:l
d(i)=norm(ad(i,:)) * m(i);
end
n=sum(d)/l;

% numerator
diff = c - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfc = max(diffn) / n;

%avg error for coulomb dielec
avgfc = mean(diffn) / n;

% numerator
diff = c1 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfc1 = max(diffn) / n;

%avg error for coulomb dielec
avgfc1 = mean(diffn) / n;

% numerator
diff = c2 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfc2 = max(diffn) / n;

%avg error for coulomb dielec
avgfc2 = mean(diffn) / n;

% numerator
diff = c3 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfc3 = max(diffn) / n;

%avg error for coulomb dielec
avgfc3 = mean(diffn) / n;

% numerator
diff = h1 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfh1 = max(diffn) / n;

%avg error for coulomb dielec
avgfh1 = mean(diffn) / n;

% numerator
diff = h2 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfh2 = max(diffn) / n;

%avg error for coulomb dielec
avgfh2 = mean(diffn) / n;

% numerator
diff = h3 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfh3 = max(diffn) / n;

%avg error for coulomb dielec
avgfh3 = mean(diffn) / n;

% numerator
diff = h4 - ad;
diffn = zeros(l,1);
for i=1:l
diffn(i)=norm(diff(i,:)) * m(i);
end

% max error for coulomb dielec
maxfh4 = max(diffn) / n;

%avg error for coulomb dielec
avgfh4 = mean(diffn) / n;


