function [w2,V]=dk(k)
global npa mass cn fc;
fa=exp(i*2*pi*(fc(:,1:3)*k'));
co=[fa.*fc(:,4) fa.*fc(:,5) fa.*fc(:,6) fa.*fc(:,7) fa.*fc(:,8) fa.*fc(:,9) fa.*fc(:,10) fa.*fc(:,11) fa.*fc(:,12)];
D=zeros(3*npa,3*npa);
p=0;l=1;
for n=1:npa,
for m=1:n-1,
q=p+cn(l);
A=reshape(sum(co(p+1:q,:)),3,3)/sqrt(mass(n)*mass(m));
D(3*m-2:3*m,3*n-2:3*n)=A;
D(3*n-2:3*n,3*m-2:3*m)=A';
p=q;l=l+1;
end;
q=p+cn(l);
D(3*n-2:3*n,3*n-2:3*n)=reshape(sum(co(p+1:q,:)),3,3)/mass(n);
p=q;l=l+1;
end;
[V,W2]=eig(D);
w2=real(diag(W2));
