%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab code to extract atoms'configuration and Hessian matrix data from MD
% ++ and obtain dispersion relation curves for graphene sheet 
%
%   Yanming Wang, Saad Bahmla and Prof.Wei Cai
%   June 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize variables
lattice_const = 2.424;
%convert mass into eV/(A/ps)^2
m = 12.0107 * 1.6605e-27 * (1/1.6022e-19) / (1e10/1e12)^2; 
nsample = 101;
GAMMA = [0 0 0];
H = [sqrt(3)/2 1/2 0] * ( pi/lattice_const);
P = [sqrt(3)/3 1 0]  * 1.15 * ( pi/lattice_const);
% choose approperiate k-vector values as x-axis of the plot
ksample = [linspace(1,0,nsample)'* GAMMA    + linspace(0,1,nsample)'* H;
           linspace(1,0,nsample)'* H + linspace(0,1,nsample)'* P;
           linspace(1,0,nsample)'* P + linspace(0,1,nsample)'* GAMMA];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reading data from MD++ and generating K martix for the calculation
datafile='hessian.out';
data = load(datafile);
fid = fopen(['perf.cn']);
N = str2num(fgetl(fid));
s = zeros(N,3);
for i=1:N
    a = fgetl(fid);
    a2=str2num(a);
s(i,:)= a2;
end
lattice_matrix = zeros(3,3);
for i=1:3
    a = fgetl(fid);
    a2=str2num(a);
  lattice_matrix(i,:) = a2;    
end
fclose(fid);
K = reshape(data',N*3,N*3);
K = (K+K')/2;
eig(K);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Label the atoms into the right pairs and set reference atoms
flag = zeros(N,2);
Nx=12;
Ny=12;
for i=1:(N-N/Nx)/4
    id1 = 4*i - 1;
    id2 = 4*i + N/Nx;
    flag(id1,:) = [1 id1];
    flag(id2,:) = [2 id1];
end
for i=1:N/4
    id1 = 4*i - 3;
    id2 = 4*i - 2;
    flag(id1,:) = [1 id1];
    flag(id2,:) = [2 id1];
end
for i=1:(N-N/Nx*(Ny-1))/4
    id1 = 4*i - 1 + N/Nx*(Ny-1);
    id2 = 4*i;
    flag(id1,:) = [1 id1];
    flag(id2,:) = [2 id1];
end
disp(sprintf('number of atoms: N = %d',N));

% type flag to help you select j0 and j02
j0=1; j02=2;
if (flag(j0,1)~=1) || (flag(j0,2)~=j0) || (flag(j02,1)~=2) || (flag(j02,2)~=j0)
    disp('wrong choices of reference atoms!');
    flag(j0,:)
    flag(j02,:)
    return
end
ref1=s(j0,:);
ref2=s(j02,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop starts here

for p=1:length(ksample(:,1));
    k = ksample(p,:);
    Q=zeros(6,6);
    for i=1:N,
        L1 = [ K(3*j0-2, 3*i-2)  K(3*j0-2, 3*i-1)  K(3*j0-2, 3*i)
               K(3*j0-1, 3*i-2)  K(3*j0-1, 3*i-1)  K(3*j0-1, 3*i)
               K(3*j0,   3*i-2)  K(3*j0,   3*i-1)  K(3*j0,   3*i)];
        L2 = [ K(3*j02-2, 3*i-2) K(3*j02-2, 3*i-1) K(3*j02-2, 3*i)
               K(3*j02-1, 3*i-2) K(3*j02-1, 3*i-1) K(3*j02-1, 3*i)
               K(3*j02,   3*i-2) K(3*j02,   3*i-1) K(3*j02,   3*i)];
        if (flag(i,1) == 1)
            L = [L1, zeros(3); L2, zeros(3)];
        elseif (flag(i,1) == 2)
            L = [zeros(3), L1; zeros(3), L2];
        else
            disp('something is terribly wrong...');
        end
        id1 = flag(i,2);
        ds = s(id1,:)' - ref1';
        ds = ds - round(ds);
        dr = lattice_matrix * ds;
        Q = Q + exp(1i*k*dr) * L;       
    end
    
% solve eigenfunction and obtain six frequency values at each k.
    [v,d]=eig(Q);
    w1(p,1)=abs(sqrt(d(1,1)/m))/(2*pi);
    w2(p,1)=abs(sqrt(d(2,2)/m))/(2*pi);
    w3(p,1)=abs(sqrt(d(3,3)/m))/(2*pi);
    w4(p,1)=abs(sqrt(d(4,4)/m))/(2*pi);
    w5(p,1)=abs(sqrt(d(5,5)/m))/(2*pi);
    w6(p,1)=abs(sqrt(d(6,6)/m))/(2*pi);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
k_vec=[[0:100]*1,101+[0:100]*0.5,101*1.5+[0:100]*1.732/2];
plot(k_vec,w1,'.',k_vec,w2,'.',k_vec,w3,'.',k_vec,w4,'.',k_vec,w5,'.',k_vec,w6,'.');
set(gca,'XTick',[]);
set(gca,'XTick',[0 100 151 238.1]);
set(gca,'XtickLabel',{'','M','K','','N'});
xlim([0 k_vec(end)]);
peak = [max(w1),max(w2),max(w3),max(w4),max(w5),max(w6)];
height=max(peak)+5;
line([100;100],[0;height]);
line([151;151],[0;height]);
ylabel('Frenquency (THz)','fontsize',12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%