clear all; clc;

%pdbfile = input('Please enter pbd file (for example, 1YXT)','s');
pdbfile = '1YXT';
pmin = 90; pmax = 120; %define P-loop residues

fileIn = [pdbfile,'.pdb'];
%fileOut = [pdbfile,'_network.dat'];
fileOutRes = [pdbfile,'_residues.dat'];

xray = pdbread(fileIn);
n = length(xray.Model.Atom); %number of atoms in pdb file
j = 1;

for i=1:n
   if ((xray.Model.Atom(1,i).resSeq>=pmin) && (xray.Model.Atom(1,i).resSeq<=pmax))
        patoms(j)=i; %indeces of atoms of P-loop residues
        res(j)=xray.Model.Atom(1,i).resSeq;
        j = j+1;
    end
end

%% extract residue coordinates
crd = zeros(3,length(res));
for i = 1:length(patoms)
    crd(1,i) = xray.Model.Atom(1,patoms(i)).X;
    crd(2,i) = xray.Model.Atom(1,patoms(i)).Y;
    crd(3,i) = xray.Model.Atom(1,patoms(i)).Z;
end

%% extract and include ligand coordinates
%ligand is treated as an additional residue

ligandid = xray.Heterogen.hetID;
ligandid = 'ANP';
%p = length(xray.Model.HeterogenAtom);
for i = 1:length(xray.Model.HeterogenAtom)
    a = xray.Model.HeterogenAtom(1,i).resName;
    if size(xray.Model.HeterogenAtom(1,i).resName) == size(ligandid)
       if (xray.Model.HeterogenAtom(1,i).resName == ligandid)
         crd(1,end+1) = xray.Model.HeterogenAtom(1,i).X;
         crd(2,end) = xray.Model.HeterogenAtom(1,i).Y;
         crd(3,end) = xray.Model.HeterogenAtom(1,i).Z;
         res(end+1) = pmax+1;
       end
   end
end

pmax = pmax+1; %added ligand

%% calculate distances between each atom
%p1 = patoms(1); p2 = patoms(end); %first and last atoms of P-loop
d = zeros(length(res));
cutoff = 4.5;

for i=1:length(res) %first atoms
    for j=1:length(res) %second atom
        d(i,j) = sqrt((crd(1,j)-crd(1,i))^2+(crd(2,j)-crd(2,i))^2+(crd(3,j)-crd(3,i))^2);
        
        %assign 1 if residues 'interact' and 0 if not
        if d(i,j)<=cutoff
            d(i,j)=1;
              else d(i,j)=0;
        end
        
         %eliminate 'self' interactions
        if  (res(i)==res(j))
            d(i,j)=0;
        end
    end
end

%% From atoms to residue ineractions
%index[first, last] - a set of atoms in the residue
for i = pmin:pmax
    k = find(res==i);
    index(i-pmin+1,:)=[k(1),k(end)];
end

resmatrix = zeros(pmax-pmin+1);
for i = 1:pmax-pmin+1
    for j = 1:pmax-pmin+1
        if find(d(index(i,1):index(i,2),index(j,1):index(j,2))==true)
            resmatrix(i,j)=1;
        end
    end
end

%% results --> files

%atom interactions
%[k,l,val] = find(d);
%fid = fopen(fileOut,'w');
%for i = 1:length(k)
%    fprintf(fid,'%d %d %f\n', k(i),l(i),val(i));
%end
%fclose(fid);

%residue interactions
[k,l,val] = find(resmatrix);
fid = fopen(fileOutRes,'w');
for i = 1:length(k)
    fprintf(fid,'%d %d %f\n', k(i),l(i),val(i));
end
fclose(fid);