% script Local_Moves; 
%                    perform localized moves in a polymer chain.
% >> Perform moves of the chain that involve only 
%     a subset of atoms. Augment to avoid steric clashes.  <<
% PROBLEM 1: given chain, deform torsion angles so first and last atoms
%            approach to distance 1
% PROBLEM 2: further deform torsion angles to allow joining of first and
%            last atoms with correct bond angle
% PROBLEM 3: introduce arbitrary moves of closed chain, centered about any
%            given torsion angle. Find maximally distorting infinitesimal 
%            deformation (leading singular vector).
%
% In a polymer successive bond-angles are fixed, equal to >> theta <<, 
% while torsion angles phi_j are entirely free. In a polypeptide chain
% (where the sequence N-Ca-C' repeats, with the dihedral Ca-C'-N-Ca (W)
% fixed, the angles Phi, Psi, of rotation about the two axes off the Ca
% nodes are rather awkward degrees of freedom, due to their disconforming
% with the overall chain geometry. Here we perform moves of the chain
% which respect the local steric constraints on the backbone level (i.e.
% sidechain clashes are not addressed, nor are long-range clashes).
%  In pursuing general answers to the 3 problems above, consider first all
% possible states of a tripeptide with two end-bonds fixed in space. The system
% is composed of 9 nodes, with bonds 1-2 and 8-9 fixed. As a peptide, the
% dihedrals around 3-4, 6-7 are fixed in the trans value Pi.
clear all
close all
view(10,180)
[a,b,c]=sphere(10);                 %           .
bb = 1;                             %        . theta. 
alpha = pi-109*pi/180;              % bb/2 .            .  bb
ct=cos(alpha); st=sin(alpha);       %    .     .5*dd        .
dd = bb*(5+4*ct)^(1/2);             %  *-----------------------*
omega = pi/120;                     % dd = bb*sqrt(5+4*cos(theta))
N = 9; %
phi0 = asin(st*bb/(2*dd));
d0   = bb*sin(phi0); 
RR   = 3;%input('enter distance Ca1-Ca3');%Ca2
hh   = (dd^2-RR^2/4)^(1/2);%            /|\
%                                      / | \
%                                     /  |  \    
R=zeros(3,N);%                    dd /   |   \ dd
B=zeros(3,N);%                      /    |    \
% define Cai in body frame         /     |hh   \
R(:,8) = [RR;0;0];%               /      |      \         
%Calculate position of Ca2,3:  Ca1-------|-------Ca3
R(1,8) = RR;%                          RR
R(1,5) = RR/2;
R(2,5) = hh;
% Now define parameters for the remaining atoms
% Cone angle and radius of motion for atom i
  phi = zeros(N,1)+phi0;
  phi(1) = .5;%input('enter angle of left end');
  phi(9) = .8;%input('enter angle of right end');
  phi(2) = 0; phi(5) = 0; phi(8) = 0;
  %-----------------------------------------
  dm = zeros(N,1)+d0;
  dm(1) = bb*sin(phi(1));
  dm(9) = bb*sin(phi(9));
  dm(2) = 0; dm(5) = 0; dm(8) = 0;
  %-----------------------------------------
  % Define the axes of rotation for each atom
  %    Length
  L0 = (1-dm.^2).^(1/2);
  %    Vectors
  L(:,1) = [-L0(1); 0; 0];
  L(:,2) = [0;0;0];
  L(:,3) = L0(3)*[RR/(2*dd); hh/dd;0];
  L(:,4) = -L(:,3);
  L(:,5) = [0;0;0];
  L(:,6) = L0(6)*[RR/(2*dd);-hh/dd;0];
  L(:,7) = -L(:,6);
  L(:,8) = [0;0;0];
  L(:,9) = [L0(9);0;0];
  %-----------------------------------------
  % define biorthogonal vector used for circle
  %r0 = zeros(3,N);
  r2 = zeros(3,N);
  for i = 1:N
     r2(:,i) = cross(L(:,i),[0;0;1])/L0(i);
  end
  %-----------------------------------------
  % Define rotation angles
  pi2 = pi/2;
  tspan = [0 12];
  [theta(4,:),theta(6,:),ttt] = wheels(dm(4),dm(6),alpha,tspan);
  tt = length(ttt);
  theta(1,1:tt) = pi2;
  theta(2,1:tt) = 0;
  theta(3,:) = -pi-theta(4,:);
  theta(5,1:tt) = 0;
  theta(7,:) = pi-theta(6,:);
  theta(8,1:tt) = 0;
  theta(9,1:tt) = -pi2;
 
  %-----------------------------------------
  % Define atom positions
  for ss = 1:tt
     for i = [2,5,8]
        for j = [-1,1]
           ij = i+j;
        R(:,ij) = R(:,i)+L(:,ij)+[0;0;cos(theta(ij,ss))]+r2(:,ij)*sin(theta(ij,ss));
        end
     end
     %plot(1:N,a0(1:N))
     %figure
     %clf
     %hold on
     %axis([-1,4,-1,4,-2.5,2.5])
     %for i = 1:N
        %surf(.05*a+R(1,i),.05*b+R(2,i),.05*c+R(3,i));
     %end
     %view((ss/tt)*360,30)
    % plot3(R(1,:),R(2,:),R(3,:),'y-',R(1,4:2:6),R(2,4:2:6),R(3,4:2:6),'r-')
    %plot3(R(1,4:2:6),R(2,4:2:6),R(3,4:2:6),'r-')
    distance(ss) = norm(R(:,4)-R(:,6));
     %hold off
     %pause(1)
  end
  plot(ttt,distance)
%for kk = 1:40
%   pause(1)
%  view(6*kk,50)
%end
