% 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).
%
% For simplicity, consider a polymer where successive 
% bond-angles are fixed, equal to >> theta <<, while torsion
% angles phi_j are entirely free. 
%                                               7
%        ____________i+2                       /
%        |            ^                       /
%        |            |                      /\theta=pi-theta'
%        |    \ b(i)  |____n(i+1) --------->----------> x-axis
%        i-----\----->i+1                   \    b(i;0)
%       / ^   phi(i)   \                     \
%      /   \            \                     \
%     n(i)  \            \                     V z-axis
%           i-1------------         
%                                 q(i-1) * p(i)* b(i;0) = b(i)
%      q(i)   = q(i-1)   * p(i;0)
%      p(i;0) = d(i-1;0) * a(i;0)
%      d(i;0) = [cos(phi(i)/2),sin(phi(i)/2),0,0]   Rx(phi(i))
%      a(i;0) = [cos(theta/2),0,0,-sin(theta/2)]    Rz(-theta)
%      n(i)   = q(i-1)*d(i-1;0)*e3
%      b(i)   = q(i)*e1 = q(i-1)*d(i-1;0)*a(i;0)*e1
%      b(0)   = e1 (unused arm on Atom 1 afixed to the x-axis;
%                   other arm at a(1)*e1, d(0) = I )
%      phi(1), phi(N) may be defined but have no effect                        
%                                  
%                                 |sin(phi(i))| = ||n(i)Xn(i+1)||
%    R(i+1) = R(i) + b(i)
%    n(i+1) = b(i) X b(i+1)
%
% CLOSURE CONDITIONS
%
%  I.               R(1)=R(N+1)
% II.               b(N)=e1
%      
%    
% To place the chain given the dihedral angles phi(1:n)
% proceed from atoms 1 to N:
%   
%   
% 
clear all
[a,b,c]=sphere(10);
N = 5;
omega = .035123*pi;
theta = pi/3;
r1=0;r2=0;r3=0;
R=eye(4);
figure
% Build Initial Configuration as a Spiral Staircase
for i = 1:N
   phi(i) = omega;
   Q = rotmat(phi(i),theta,1);
   R=R*Q;
   X(i) = r1 + R(1,4);Y(i) = r2 + R(2,4); Z(i) = r3 + R(3,4);
   r1   = X(i); r2 = Y(i); r3 = Z(i);
end
hold on
for i = 1:N
   surf(.1*a+X(i),.1*b+Y(i),.1*c+Z(i));
end
plot3(X,Y,Z,'y-')
hold off
for kk = 1:30
   r1=0;r2=0;r3=0;
   R=eye(4);
   for i = 1:N
      phi(i) = omega*kk/2;
      if i==22, phi(i) = 3*omega*kk;end
      if i==30, phi(i) = i*i*omega*kk;end
      Q = rotmat(phi(i),theta,1);
      R=R*Q;
      X(i) = r1 + R(1,4);Y(i) = r2 + R(2,4); Z(i) = r3 + R(3,4);
      r1   = X(i); r2 = Y(i); r3 = Z(i);
   end
   figure
   view(6*kk,50)
   hold on
   for i = 1:N
      surf(.1*a+X(i),.1*b+Y(i),.1*c+Z(i));
   end
   plot3(X,Y,Z,'y-')
   pause(1)
   hold off
end

hold off

% Idea: state can be expressed in terms of local quaternions. Each
% atom and bond can be reconstructed, given the torsion angles.
% To ensure that the localization conditions are enforced, assume
% the backbone quaternions are numbered consecutively,
% q1,q2,...,qm; let the total quaternions be Qn=q1*q2*...*qn.
% I want a move that fixes QN1, QN2, with N2=N1+L; L is the 
% length of a (small) segment of the backbone that will be
% perturbed, subject to the constraint that QN1, QN2 remain
% fixed (to be more precise, RN1 and RN2 must also remain
% fixed; this just adds a few more constraints, but the 
% analysis is similar, so let's leave it out of the 
% discussion!). Define q(N1+k):=p(k) (then p(L)=q(N2))
% Since QN2 = QN1*p(1)*...*p(L) let's assume that
% all the intermediate quaternions, p(k) are perturbed
% slightly:  p(k) -->p(k) + d(k), d(k) small. Then,
% QN2 = QN1*(p(1)+d(1))*...*(p(L)+d(L)) = QN1*p(1)*...*p(L)
% and for small d(k) expand keeping only the linear terms.
% This gives a rank=4 constraint on the d(k); it's null
% vectors are all the allowable deformations that leave
% the two ends of the chain fixed (if we also include the
% position vectors, then we get 6 conditions which are
% again linear for small deformations). 


%function Q = rotmat(phi,theta)
% st=sin(theta);ct=cos(theta);sp=sin(phi);cp=cos(phi);
% Q(1,1)=-ct;   Q(1,2)=-st;    Q(1,3)=0;   Q(1,4)=-b*ct;
% Q(2,1)=st*cp; Q(2,2)=-ct*cp; Q(2,3)=-sp; Q(2,4)=b*st*cp;
% Q(3,1)=st*sp; Q(3,2)=-ct*sp; Q(3,3)=cp;  Q(3,4)=b*st*sp;
% Q(4,1)=0;     Q(4,2)=2;      Q(4,3)=0;   Q(4,4)=1;
