% Program computes frequencies and participation ratios for vibrational normal modes of 
% ordered 4-atom chain through diagonalization of its dynamical matrix

%PARAMETERS:
n_atoms=4; %number of atoms
K0=1; %spring constant normalization factor (when plotting reintroduce the units for 
%the frequency as sqrt(Ksi/Msi) where Ksi=10.6 eV/A^2 and Msi is sillicon atom mass) 
random_vector=zeros(1,n_atoms); %no disorder
K=K0*(ones(1,n_atoms)+random_vector); %spring constants are K0*(1+random_vector(i))

%VARIABLES:
dmatrix=zeros(n_atoms,n_atoms); %set all elements of the dynamical matrix to zero

%Construct dynamical matrix using spring constant vector K
dmatrix(1,n_atoms)=-K(n_atoms); %appears due to periodic boundary conditions
dmatrix(n_atoms,1)=-K(n_atoms); %also due to periodic bc (dmatrix is symmetric!)
dmatrix(1,1)=K(n_atoms)+K(1);
dmatrix(1,2)=-K(1);
dmatrix(n_atoms,n_atoms) = K(n_atoms-1) + K(n_atoms);
dmatrix(n_atoms,n_atoms-1) = -K(n_atoms-1);

dmatrix(2,1)=-K(1);
dmatrix(2,2)=K(1)+K(2);
dmatrix(2,3)=-K(2);
dmatrix(3,2)=dmatrix(2,3);
dmatrix(3,3)=K(2)+K(3);
dmatrix(3,4)=-K(3);


%Find Eigenvalues and Eigenstates of the dynamical matrix
[EigVec,EigVal]=eig(dmatrix);
EigVal=diag(EigVal);

freq=sqrt(EigVal); %frequencies of normal modes are square root of eigenvalues

%MATALAB shortcut to find participation ratio
EigVec=EigVec.^4; %take forth power of each element of matrix EigVec
p=sum(EigVec).^(-1); %for matrix X, sum(X) is a row vector whose components are sums over each column of X


pr_result=[freq.'; p]; %data stored into a two row matrix before saving on the disk

fid=fopen('p_vs_freq.dat','wt');
fprintf(fid,'%g  %g\n',pr_result); %frequency vs. participation ratio saved on the disk in two column data file
fclose(fid);
