function [freq] = phonon(k)
global dm;

IP = 1;
for I=1:6,
  for  J=1:6,
  A(I,J) = 0;
    maxij = dm(IP,1);
    IP = IP + 1;
      for M = 1:maxij,
       PHASE = k * dm(IP,1:3)' ;
       A(I,J) = A(I,J) + dm(IP,4)*exp(i*PHASE);
       IP = IP + 1;
      end;
 end;
end;

[V,D] = eig(A);
FREQ_K = diag(D);
for i=1:6
 if FREQ_K(i) < 0
   FREQ_K(i) = 0.;
 end;
FREQ_K(i)=sqrt(FREQ_K(i));
end;

freq = sort(FREQ_K);

return;


