function [] = constructSeqParms(seq)

%Example usage: constructSeqParms('GCTATATATATATATAGC');
% 
%This function constructs a sequence dependent shape vector and stiffness
%matrix for a rigid base DNA of a given sequence 'seq'. It uses a dimer and
%monomer dependent parameter set, loaded from the file 'paramset1.mat'.
%The resulting stiffness matrix (of the dimersion NxN, where N=12*nbp-6 and 
%nbp is the length of the sequence) is saved in the variable 'stiff'
%and the shape vector (of the dimention Nx1) - in the variable 'shapes'.
%The parameters in the variables 'shapes' and 'stiff' are ordered by basepairs 
%along the sequence, and for every basepair the order is:
%Intra rotations (Buckle, Propeller, Opening), Intra translations 
%(Shear, Stretch, Stagger), Inter rotations (Tilt, Roll,Twist) 
%and Inter translations (Shift, Slide, Rise). E.g. the entry shapes((i-1)*12+1) 
%corresponds to Buckle of the i-th basepair for i=1...nbp, and the entry 
%shapes((i-1)*12+7) contains the Tilt value for the i-th basepair
%step, i=1..nbp-1. The shapes are also saved in arrays reordered by a 
%parameter: the variable eta contains all Intra rotations along the
%sequence, w - Intra translations, u - Inter rotations, v - Inter
%translations.
%
%References: O. Gonzalez, D. Petkeviciute and J.H. Maddocks, A sequence-dependent
%rigid-base model of DNA, J. Chem. Phys 138 (5), 2013.

load paramset1

nbp = size(seq,2);
stiff = zeros(nbp*12-6);
c = zeros(nbp*12-6,1);

for i = 1:(nbp-1) 
    
  k = (i-1)*12 + 1;  
  stiff(k:k+17,k:k+17) = stiff(k:k+17,k:k+17) + dimer(fsi(dimer,seq(i:i+1))).b18;
  c(k:k+17) = c(k:k+17) + dimer(fsi(dimer,seq(i:i+1))).c18;
  
end; 

for i = 1:nbp
                      
  k = (i-1)*12 + 1;  
  stiff(k:k+5,k:k+5) = stiff(k:k+5,k:k+5) + base(fsi(base,seq(i))).b;   
  c(k:k+5) = c(k:k+5) + base(fsi(base,seq(i))).c;   
 
end; 

shapes = stiff\c;
[eta,w,u,v] = vector2shapes(shapes);

filename = strcat('constructed_seq_parms.mat');
save(filename,'stiff','shapes','eta','w','u','v');

end

%-------------------------------------------------------------

function i = fsi(struc, s)
 
n = size(struc);

for j = 1:n(2)
    if(strcmpi(struc(j).S,s)) 
        i=j; 
        return 
    end;    
end    

end

%-------------------------------------------------------------

function [eta,w,u,v] = vector2shapes(y)

    n = (size(y,1)+6)/12;
    eta = zeros(n,3);
    w = zeros(n,3);
    u = zeros(n-1,3);
    v = zeros(n-1,3);
 
    for i = 1:n
        j = (i-1)*12;
        eta(i,:) = y(j+1:j+3)'; 
        w(i,:) = y(j+4:j+6)';
        if (i<n)
            u(i,:) = y(j+7:j+9)';
            v(i,:) = y(j+10:j+12)';
        end;    
    end;    

end

