% This file analyses the changes in the carbon back-bone 

load pdb1.xyz
load pdb2.xyz

X(:,1) = pdb1(:,2);
X(:,2) = pdb1(:,3);
X(:,3) = pdb1(:,4);
Y(:,1) = pdb2(:,2);
Y(:,2) = pdb2(:,3);
Y(:,3) = pdb2(:,4);

%  Bond lengths
length(size(X,1)) = 0;
for i=2:size(X)
     length(i) = norm(X(i,:)-X(i-1,:)) - norm(Y(i,:)-Y(i-1,:));
end

%  Bond angles and plane twist
angle(size(X,1)) = 0;
twist(size(X,1)) = 0;
for i=2:size(X,1)-2
alpha = X(i,:)-X(i-1,:);
beta = X(i+1,:) - X(i,:);
thetax = acos( dot(alpha,beta)/(norm(alpha)*norm(beta) ));
gamma = X(i+2,:) - X(i+1,:);
a = cross(alpha,beta);
b = cross(beta,gamma);
a = a /norm(a);
b = b / norm(b);
twistx = acos( dot(a,b));

alpha = Y(i,:)-Y(i-1,:);
beta = Y(i+1,:) - Y(i,:);
thetay = acos( dot(alpha,beta)/(norm(alpha)*norm(beta) ));
gamma = Y(i+2,:) - Y(i+1,:);
a = cross(alpha,beta);
b = cross(beta,gamma);
a = a /norm(a);
b = b / norm(b);
twisty = acos( dot(a,b));

angle(i) = thetay - thetax ;
twist(i) = twisty - twistx ;
end