%-----------------------------------------------------------------------
function proline(chi5)
%chi5 = 0;
%-----------------------------------------------------------------------
% written by Vageli Coutsias, Dept. of Math&Stats, UNM  (C) 2008
%            vageli@math.unm.edu                                
% generates conformations for proline sidechain, following procedure in
% "The flexibility in the proline ring couples to the protein backbone",
% Bosco Ho, E.A. Coutsias, C. Seok, K. Dill, 
% Protein Science, 2005, 14:1011-1018
% (several typos in the paper are corrected here-see forthcoming Errata)
% internal coordinates are consistent with any imposed value of the 
% chi5 torsion, where it is assumed that 
%               chi5 = phi + 60 degrees
% the bond lengths are fixed to canonical values (averages for TRANS-PRO):
% Ca-Cb = 1.530, Cb-Cg = 1.501, Cg-Cd = 1.513, Cd-N = 1.476, N-Ca = 1.464A
% The bond angles N -Ca-Cb = 103.7, Ca-Cb-Cg = 103.9, Cg-N -Ca = 111.5 deg
% are also fixed, while the remaining bond angles Cb-Cg-Cd, Cg-Cd-N
% are free to vary, as are the torsions chi1,..., chi4. Since resulting
% structures will be slightly strained, energy minimization is recommended
% to distribute strain among all the bond angles. To give a sense of the
% strain, the deviation of Cb-Cg-Cd, Cg-Cd-N from "canonical" values
% (resp. 104.5, 102.7) is also output. In general, there are two pucker 
% states, the up- and down-pucker forms, defined by the convention:
% Up:        180 > chi2 >  10 degrees
% down:     -180 < chi2 < -10 degrees
% We also identify a planar conformation by
% planar:   -10 < chi2 < 10 degrees
% for testing: +/- the value below gives planar limiting conformations
%chi5 =  31.65090679922;%92383308949865750037133693
% This program may be used and modified freely, please acknowledge use
% by citing this website and the above publication. Thanks!
%-----------------------------------------------------------------------
tol = 10^(-6);           % tolerance for differentiating structures
deg2rad = pi/180;
alpha = 103.7;           % N -Ca-Cb
beta  = 103.9;           % Ca-Cb-Cg
gamma = 104.5;           % Cb-Cg-Cd  not imposed
delta = 102.7;           % Cg-Cd-N   not imposed
nu    = 111.5;           % Cg-N -Ca
AB    = 1.530;           % Ca-Cb
BG    = 1.501;           % Cb-Cg
GD    = 1.513;           % Cg-Cd
DN    = 1.476;           % Cd-N
NA    = 1.464;           % N -Ca
% 1. Solution of geometrical problem: given chi5, find value of chi1
%    so that all bond lengths and bond angles N-Ca-Cb, Ca-Cb-Cg, Cg-N-Ca
%    are held fixed to (arbitrary) prescribed values
%
%           SETUP FOR TETRAHEDRAL EQUATION
%
%                            |.   
%                           |        .        |\
%                          |                 D .\
%                 <-------O<---------------------O 
%                       N| nu             .       \
%                       |             .            \     ^
%                      |           .                \   /
%               sigma |        .                     \ /
%                    | xi   .                         O G
%                   |     .   theta         .        / .
%                  |    .            .              /  .
%                 |  .       .                     /   .
%                v.  .          eta         beta  / B  .
%               O------------------------------->O---------------->
%              |A     tau      
%             |pi-alpha
%            v
  ar = alpha*deg2rad; br = beta*deg2rad; 
  nr = nu*deg2rad; sigma = chi5*deg2rad; 
  [ca, sa] = COS_SIN(ar); 
  [cb, sb] = COS_SIN(br);
  [cn, sn] = COS_SIN(nr);
%  trigonometry: determine AC, AD, xi, eta and theta
%  AG^2 =      AB^2 + BG^2 - 2 AB BG cos(beta)
   AG   = sqrt(AB^2 + BG^2 - 2*AB*BG*cb);
%  AD^2 =      NA^2 + DN^2 - 2 NA DN cos(nu)
   AD   = sqrt(NA^2 + DN^2 - 2*NA*DN*cn);
%  sin(xi)  = DN*sin(nu)/AD
   sxi      = DN*sn/AD;
   cxi      = sqrt(1-sxi^2);  % acute angle, use positive square root
%  sin(eta) = BC*sin(beta)/AC
   seta     = BG*sb/AG;
   ceta     = sqrt(1-seta^2); % acute angle,  ""  ""
   [cs, ss] = COS_SIN(sigma);
%  GD^2 = AG^2 + AD^2 - 2*AG*AD*cos(theta)
   ct = (AG^2 + AD^2 - GD^2)/(2*AG*AD);
%Tetrahedral equation; following Bricard, JMathPurAppl, 1897, 3, 113-150
%english translation: http://www.math.unm.edu/~vageli/papers/bricard.pdf
%     sin(xi)*sin(eta)*cos(alpha)*cos(sigma)*cos(tau)
%    -cos(xi)*sin(eta)*sin(alpha)           *cos(tau)
%    +sin(xi)*sin(eta)           *sin(sigma)*sin(tau)
%    -sin(xi)*cos(eta)*sin(alpha)*cos(sigma)
%    -cos(xi)*cos(eta)*cos(alpha)                     + cos(theta) = 0
%(sign change from Bricard's form due to different definition of torsion
% tau, here taken as the dihedral( N-A-B-G) =  - Bricard's definition)
% solve for tau:
%  ( sin(xi)*cos(alpha)*cos(sigma)-cos(xi)*sin(alpha))*cos(tau)
%+(sin(xi)                 *sin(sigma)                   )*sin(tau)
%=(cos(eta)*(sin(xi)*sin(alpha)*cos(sigma)+cos(xi)*cos(alpha)) 
%- cos(theta))/sin(eta)
%    A * cos(tau) + B * sin(tau) = C
     A =  sxi*ca*cs-cxi*sa;
     B =  sxi*ss;
     C = (ceta*(sxi*sa*cs + cxi*ca)-ct)/seta;
     SAB = sqrt(A^2+B^2);
     tau0 = atan2(B,A);
     if (abs(C/SAB)>1)
         fprintf('\n impossible configuration \n')
         return
     else
         tau1 = acos(C/SAB);
     end
     if (abs(tau1) <= tol | abs(abs(tau1)-pi) <= tol )
     %if ( tau1 == 0 | abs(tau1) == pi)
         tau = tau0+abs(tau1); 
         flag = 0;    % unique solution
     else
         taup = tau0 + abs(tau1);
         taum = tau0 - abs(tau1);
         flag = 1;    % two solutions
     end
% 2. Cartesian and Internal coordinates for the compatible structures
     b_ang(1) = 0; % ignored
     b_ang(2) = nr;
     b_ang(3) = ar;
     b_ang(4) = br;
     b_ang(5) = 0; % ignored
     b_len(1) = DN;
     b_len(2) = NA;
     b_len(3) = AB;
     b_len(4) = BG;
     b_len(5) = GD; % ignored
     t_ang(1) = 0;  % ignored
     t_ang(2) = sigma;
     t_ang(4) = 0;  % ignored
     t_ang(5) = 0;  % ignored
     R0 = zeros(3,1); Q = eye(3); %R0(1) = -NA;
     if (flag == 0) % only one structure returned
fprintf('\n---------------Single conformation------------------------\n');
         t_ang(3) = tau;
         atoms = CHAIN_XYZ ( 5, b_len, b_ang, t_ang, 0, R0, Q );
         atoms0 = [atoms(1:3,2:5),atoms(1:3,1)];
         [t_ang0, b_ang0, b_len0, R0, Q] = CHAIN_TORS( 5,  atoms0 );
         print_internals(b_len0, b_ang0, t_ang0, gamma, delta)
         R0 = zeros(3,1); Q = eye(3); R0(1) = -NA;
         atoms0 = CHAIN_XYZ(5,b_len0,b_ang0,t_ang0,1,R0,Q);
         print_cartesians(atoms0)
     else
fprintf('\n---------------first conformation-------------------------\n');
         t_ang(3) = taup;
         atoms = CHAIN_XYZ ( 5, b_len, b_ang, t_ang, 0, R0, Q );
         atomsp = [atoms(1:3,2:5),atoms(1:3,1)];
         [t_angp, b_angp, b_lenp, R0, Q] = CHAIN_TORS( 5,  atomsp );
         print_internals(b_lenp, b_angp, t_angp, gamma, delta)
                  R0 = zeros(3,1); Q = eye(3); R0(1) = -NA;
         atomsp = CHAIN_XYZ(5,b_lenp,b_angp,t_angp,1,R0,Q);
         print_cartesians(atomsp)
fprintf('\n-------------second conformation--------------------------\n');     
         t_ang(3) = taum;
         atoms = CHAIN_XYZ ( 5, b_len, b_ang, t_ang, 0, R0, Q );
         atomsm = [atoms(1:3,2:5),atoms(1:3,1)];
         [t_angm, b_angm, b_lenm, R0, Q] = CHAIN_TORS( 5,  atomsm );
         print_internals(b_lenm, b_angm, t_angm, gamma, delta)
         R0 = zeros(3,1); Q = eye(3); R0(1) = -NA;
         atomsm = CHAIN_XYZ(5,b_lenp,b_angm,t_angm,1,R0,Q);
         print_cartesians(atomsm)
     end
fprintf('\n----------------------------------------------------------\n');
return
%-----------------------------------------------------------------------
function print_internals(b_len, b_ang, t_ang, gamma, delta)
% convert all angles to degrees
         rad2deg = 180/pi;
         b_ang = b_ang*rad2deg;
         t_ang = t_ang*rad2deg;
         % identify Pucker
         if (t_ang(3)<10 & t_ang(3)> -10)
             pucker = 'PLANAR';
         else if (t_ang(3)>10 & t_ang(3)<180)
                 pucker = '  UP  ';
             else
                 pucker = ' DOWN ';
             end
         end
         fprintf('\n  %s PUCKER  \n', pucker);
         fprintf('\n  INTERNAL COORDINATES \n');
fprintf('\n   bond length         bond angle                torsion\n');
  fprintf(' Ca-Cb  %6.3f     N-Ca-Cb  %6.3f   N-Ca-Cb-Cg   %6.3f \n',...
         b_len(2),b_ang(2),t_ang(2));
  fprintf(' Cb-Cg  %6.3f    Ca-Cb-Cg  %6.3f   Ca-Cb-Cg-Cd  %6.3f \n',...
         b_len(3),b_ang(3),t_ang(3));
  fprintf(' Cg-Cd  %6.3f    Cb-Cg-Cd  %6.3f   Cb-Cg-Cd-N   %6.3f \n',...
         b_len(4),b_ang(4),t_ang(4));
  fprintf(' Cd-N   %6.3f    Cg-Cd-N   %6.3f   Cg-Cd-N-Ca   %6.3f \n',...
         b_len(5),b_ang(5),t_ang(5));
  fprintf('  N-Ca  %6.3f    Cd-N-Ca   %6.3f   Cd-N-Ca-Cb   %6.3f \n',...
         b_len(1),b_ang(1),t_ang(1));
         fprintf('\n bond angle variation \n')
         fprintf(' D<Cb-Cg-Cd> =  %6.3f   D<Cg-Cd-N> =   %6.3f   \n',...
                   b_ang(4)- gamma,b_ang(5)-delta);
%-----------------------------------------------------------------------
function print_cartesians(atoms)
    fprintf('\n  CARTESIAN  COORDINATES \n');
    fprintf('\n                X            Y           Z   \n');
    fprintf('  Ca        %6.3f       %6.3f      %6.3f     \n',...
                             atoms(1,2),atoms(2,2),atoms(3,2));
    fprintf('  Cb        %6.3f       %6.3f      %6.3f     \n',...
                             atoms(1,3),atoms(2,3),atoms(3,3));   
    fprintf('  Cg        %6.3f       %6.3f      %6.3f     \n',...
                             atoms(1,4),atoms(2,4),atoms(3,4));
    fprintf('  Cd        %6.3f       %6.3f      %6.3f     \n',...
                             atoms(1,5),atoms(2,5),atoms(3,5));
    fprintf('  N         %6.3f       %6.3f      %6.3f     \n',...
                             atoms(1,1),atoms(2,1),atoms(3,1));
%-----------------------------------------------------------------------
% MODULE GEOM
%-----------------------------------------------------------------------
% CONTAINS
%-----------------------------------------------------------------------
%
    function atoms = CHAIN_XYZ ( n, b_len, b_ang, t_ang, space, R0, Q )
%
%   calculates cartesian from internal coordinates  for a linear 
%   chain optimally using backwards rotation. For a discussion see 
%   "Efficiency of Rotational Operators for Geometric Manipulation 
%   of Chain Molecules", Chaok Seok and E.A. Coutsias, Bulletin of 
%   the Korean Chemical Society, 28(10), p.1705-1708, (2007).
%-----------------------------------------------------------------
%   total cost:
%   (2 mul + 1 add)                           first 3 atoms
%   (27 mul + 15 add )*(n-4)   mid-chain
%     4 mul +  3 add + 2 integer add          terminal atom
%   2*(n-3) + 1    cos_sin calls (each:  5 mul + 2 add + 1 dtan ) 
%-----------------------------------------------------------------
%
%   INPUT:
%     n - number of atoms
%     t_ang - torsion angles, b_ang - bond angles, b_len - bond lengths
%     angles are input in radians
%
%   OUTPUT: atoms in BODY frame of first three atoms
%           R1 = 0
%           R2 positive x-axis
%           R3 positive xy-halfplane
%
% atom #1
       atoms(1:3,1:3) = 0;
% atom #2
       atoms(1,2) = b_len(1);
% atom #3
       [ca, sa] =  COS_SIN(b_ang(2));
       atoms(1,3) = b_len(1) - b_len(2) * ca;
       atoms(2,3) =            b_len(2) * sa;
% initialize rotation matrix
       U(1,1)   = -ca;
       U(1,2)   = -sa;
       U(2,1)   =  sa;
       U(2,2)   = -ca;
       U(1:2,3) = 0;
       U(3,1:2) = 0; 
       U(3,3)   = 1;
% all other atoms:
       n1 = n-1;
       n2 = n-2;
       j1 = 3;
       j2 = 2;
       for j = 4:n1
          [ca, sa] =  COS_SIN(b_ang(j1));
          [ct, st] =  COS_SIN(t_ang(j2));
          s(1:3,1)  = U(1:3,2);
          U(1:3,2) = s(1:3)*ct + U(1:3,3)*st;
          U(1:3,3) =-s(1:3)*st + U(1:3,3)*ct;
          s(1:3,1)  = U(1:3,2);
          U(1:3,2) =-s(1:3)*ca - U(1:3,1)*sa;
          U(1:3,1) = s(1:3)*sa - U(1:3,1)*ca;
          atoms(1:3,j) = atoms(1:3,j1) + b_len(j1) * U(1:3,1);
          j2 = j1;
          j1 = j;
       end
          [ca, sa] =  COS_SIN(b_ang(n1));
          [ct, st] =  COS_SIN(t_ang(n2));
          v(1) =-b_len(n1) * ca;
          ba   = b_len(n1) * sa;
          v(2) = ba*ct;
          v(3) = ba*st;
          s = U*v';
          atoms(1:3,n) = atoms(1:3,n1) + s(1:3);
          % now set to frame (default is R0 = 0, Q = Identity)
       if (space == 1)
           for i = 1:n
               atoms(:,i) = R0 + Q*atoms(:,i);
           end
       end
%
%      end function CHAIN_XYZ
%-----------------------------------------------------------------------
   function [t_ang, b_ang, b_len, R0, Q] = CHAIN_TORS(n, atoms )
%
%       (j-1) --------  (j)
%                b_ang(j) \
%                 t_ang(j) \b_len(j)
%                           \
%                           (j+1) ------  (j+2)
%
%       calculates torsions, bond angles, bond lengths from atoms
%       for a linear chain
%       returns angles in radians
%
%       INPUT:
%       n - number of atoms
%       atoms - atoms in chain (always treated as a ring for
%                               purposes of loop closure)
%       OUTPUT: t_ang, b_ang, b_len (if not a ring, then the values
%                                    computed outside are irrelevant)
%      The chain can be reconstructed in Cartesian form without these
%      six additional DOFs.
       R0 = atoms(:,1);
       Q = FRAME(atoms(1:3,1:3)); 
%      An origin and orientation are also calculated, so the exact 
%      cartesian configuration may be reconstructed if so desired
 for j = 1:n-1
     v = atoms(:,j)- atoms(:,j+1);
     b_len(j) = sqrt ( dot(v,v) );
 end
     v = atoms(:,n)- atoms(:,1);
     b_len(n) = sqrt ( dot(v,v) );
 for j = 2:n-1
     b_ang(j) = BOND_ANGLE( atoms(:,j-1), atoms(:,j), atoms(:,j+1) );
 end
     b_ang(n) = BOND_ANGLE( atoms(:,n-1), atoms(:,n), atoms(:,1) );
     b_ang(1) = BOND_ANGLE( atoms(:,n  ), atoms(:,1), atoms(:,2) );
 for j = 2:n-2
    t_ang(j)=TORSION(atoms(:,j-1),atoms(:,j),atoms(:,j+1),atoms(:,j+2));
 end
    t_ang(n-1)=TORSION(atoms(:,n-2),atoms(:,n-1),atoms(:,n),atoms(:,1));
    t_ang(n  )=TORSION(atoms(:,n-1),atoms(:,n  ),atoms(:,1),atoms(:,2));
    t_ang(1  )=TORSION(atoms(:,n  ),atoms(:,  1),atoms(:,2),atoms(:,3));
%  end function CHAIN_TORS
%-----------------------------------------------------------------------
       function U = FRAME(R)
% construct body frame for r1  origin
%                          r2  r1-->r2: positive x axis
%                          r3  on pos. y side of xy plane
%-----------------------------------------------------------------------
     U = zeros(3, 3);
      U(:, 1) = R(:, 2) - R(:, 1);
      U(:, 1) = U(:, 1) / norm(U(:, 1));

      U(:, 3) = CROSS(U(:, 1), R(:, 3) - R(:, 2));
      U(:, 3) = U(:, 3) / norm(U(:, 3));

      U(:, 2) = CROSS(U(:, 3), U(:, 1));

%     end function FRAME
% -----------------------------------------------------------------------
function [c, s] = COS_SIN(theta)
% Expects angle theta in radians
% Returns the cosine and sine of theta via tan(theta / 2).
% -----------------------------------------------------------------------                                                                             
   t = tan(theta / 2);
   t2   = t * t;
   t2p1 = t2 + 1;
   c = (1 - t2)/t2p1;
   s = 2 * t/t2p1;
   % end SUBROUTINE COS_SIN
% -----------------------------------------------------------------------
      function angle = BOND_ANGLE ( a, b, c )
% calculate bond angle in radians; interior angle, 3 atoms
      p = a-b;   q = c-b;
      r = dot(p,q)/sqrt(dot(p,p)*dot(q,q));
      angle = atan2 ( sqrt(1.-r*r), r);
      %      end function bond_angle
%-----------------------------------------------------------------------
      function chi = TORSION( a, b, c, d )
% calculate torsion angle in radians (-pi<= angle < pi)
% fails if a=b; gives 0 if c=d
       if (b == c)
           chi = pi;
       else
           r   = d - b;
           s(:,1) = c - b;
           s(:,2) = a - b;
           s(:,3) =  CROSS ( s(:,1), s(:,2));
           s(:,2) =  CROSS ( s(:,3), s(:,1));
           for i = 2:3
               f = sqrt( s(1,i)^2 + s(2,i)^2 + s(3,i)^2 );
               s(1:3,i) = s(1:3,i) / f;
           end
           y = r(1)*s(1,2) + r(2)*s(2,2) + r(3)*s(3,2);
           z = r(1)*s(1,3) + r(2)*s(2,3) + r(3)*s(3,3);
           chi = atan2 ( z, y );
       end
%      end function TORSION
%-----------------------------------------------------------------------
      function  w = CROSS ( u, v )
% vector product:       c  =  a  X  b
  w = [ u(2)*v(3)-u(3)*v(2); u(3)*v(1)-u(1)*v(3); u(1)*v(2)-u(2)*v(1) ];
%      end function cross
%-----------------------------------------------------------------------
%end MODULE GEOM
%-----------------------------------------------------------------------
