% This Matlab Script Calculates a Ramachandran Distribution for all
% backbone dihedrals
function [Bins]=RamachandranAll(fn,maxfe)

%clear

dihedraldata = load(fn);%'allDihedrals.out');
numDihedrals = 55 % Number of Dihedral Pairs 57 residues - 2 incomplete terminals


for i=1:numDihedrals
    % Offset 3 columns (col1 = step, col2 = PSI Res1, col3 = Omega Res1)
    % Only want Phi and Psi values, skip omegas
    Phis(:,i) = dihedraldata(:,i*3+1); % dihedral index 10 in ProtoMol
    Psis(:,i) = dihedraldata(:,i*3+2); % dihedral index 17 in ProtoMol
end

dataSize = size(Phis(:,1))
binRes =  5 % Bin resolution in degrees

PhiBins = -180:binRes:180;
PsiBins = -180:binRes:180;
tbinCount = size(PhiBins)
binCount = tbinCount(2)
Bins = zeros(binCount,binCount);

for j=1:dataSize
    for k=1:numDihedrals
        
        % added to shift axis to match published results
        while Phis(j,k) >= 180
            Phis(j,k) = Phis(j,k) - 360;
        end
        while Psis(j,k) >= 180
            Psis(j,k) = Psis(j,k) - 360;
        end
        while Phis(j,k) < -180
            Phis(j,k) = Phis(j,k) + 360;
        end
        while Psis(j,k) < -180
            Psis(j,k) = Psis(j,k) + 360;
        end
        % end shift
        
        curPhiBin = floor((Phis(j,k)+180)/binRes)+1;
        curPsiBin = floor((Psis(j,k)+180)/binRes)+1;
        Bins(curPsiBin,curPhiBin) = Bins(curPsiBin,curPhiBin) + 1;
    end
end

Bins;
Bins = Bins ./ sum(sum(Bins));
for i=1:binCount
    for j=1:binCount
        if Bins(i,j) > 0
            Bins(i,j) = -log(Bins(i,j));%0.6*
        end
    end
end

maxb = max(max(Bins))
if exist('maxfe') == 1
    maxb = maxfe;
end

for i=1:binCount
    for j=1:binCount
        if Bins(i,j) == 0
            Bins(i,j) = maxb;
        end
        if Bins(i,j) > maxb;
            Bins(i,j) = maxb;
        end
    end
end

%flip y data if imagesc
for i=1:(binCount/2)
    for j=1:binCount
       temp = Bins(i,j);
       Bins(i,j) = Bins(binCount-i+1,j);
       Bins(binCount-i+1,j) = temp;;
    end
end

figure(3);
load('MyColormaps','mycmap');
set(gcf,'Colormap',mycmap);
imagesc(PhiBins,PsiBins,Bins)
%mesh(PhiBins,PsiBins,Bins);
hold all
set(gca,'FontSize',17,...
  'Layer','top',...
  'XTick',[-150 -100 -50 0 50 100 150 ],...
  'XTickLabel',{'-150','-100','-50','0','50','100','150'},...
  'YTick',[-150 -100 -50 0 50 100 150 ],...
  'YTickLabel',{'150','100','50','0','-50','-100','-150'})
colorbar1 = colorbar('peer',...
  gca,[0.9225 0.131 0.04771 0.794],...
  'Box','on',...
  'Location','manual');
xlabel('\Phi')
xlim([-182.5 177.5])
ylim([-177.5 182.5])
ylabel('\Psi')
hold off

% figure(5);
% mesh(PhiBins,PsiBins,Bins);
% hold all
% xlabel('\Phi')
% xlim([-180 180])
% ylim([-180 180])
% ylabel('\Psi')
% zlabel('Free energy')
% title('Imp Solv BPTI - 1ns - All Backbone Dihedrals - 300K')
% hold off;
%print -depsc '3DimpSolWW_ramaNM.eps'
%print -dpng '3DimpSolWW_rama.png'

% figure(6);
% mb = max(max(Bins))
% minb = min(min(Bins))
% %mb = 0.0036;
% contVals = minb:(mb-minb)/100:mb;
% %contVals = 1:1:50
% %Bins = Bins / 6.15385
% contour(PhiBins,PsiBins,Bins,contVals);
% hold all
% xlabel('\Phi')
% xlim([-180 180])
% ylim([-180 180])
% ylabel('\Psi')
% title('Imp Solv BPTI - 0.5ns - All Backbone Dihedrals - 300K')
% hold off;
%print -depsc '2DimpSolWW_ramaNM.eps'
%print -depsc '2DimpSolWW_rama.png'