Ph = Pi/5;
Th = ArcCos[Cos[Ph]/(1+ Cos[Ph])];

(* vertices of icosaherdon, distorted from equilibrium positions, except for north pole *)
V0 = {0,0,1};
V1 = {Sin[Th+y1 t]Cos[x1 /Sin[Th]t], Sin[Th+y1 t]Sin[x1 /Sin[Th]t], Cos[Th+y1 t]};
V2 = {Sin[Th+y2 t]Cos[2Ph+x2 /Sin[Th] t], Sin[Th+y2 t]Sin[2Ph+x2 /Sin[Th]t], Cos[Th+y2 t]};
V3 = {Sin[Th+y3 t]Cos[4Ph+x3 /Sin[Th] t], Sin[Th+y3 t]Sin[4Ph+x3 /Sin[Th]t], Cos[Th+y3 t]};
V4 = {Sin[Th+y4 t]Cos[6Ph+x4 /Sin[Th]t], Sin[Th+y4 t]Sin[6Ph+x4 /Sin[Th]t], Cos[Th+y4 t]};
V5 = {Sin[Th+y5 t]Cos[8Ph+x5 /Sin[Th]t], Sin[Th+y5 t]Sin[8Ph+x5 /Sin[Th]t], Cos[Th+y5 t]};

(* angular deviations for displaced vertices*)
d01 = ArcCos[V0.V1]-Th;
d02 = ArcCos[V0.V2]-Th;
d03 = ArcCos[V0.V3]-Th;
d04 = ArcCos[V0.V4]-Th;
d05 = ArcCos[V0.V5]-Th;

d12 = ArcCos[V1.V2]-Th;
d13 = ArcCos[-V1.V3]-Th;
d14 = ArcCos[-V1.V4]-Th;
d15 = ArcCos[V1.V5]-Th;

d23 = ArcCos[V2.V3]-Th;
d24 = ArcCos[-V2.V4]-Th;
d25 = ArcCos[-V2.V5]-Th;

d34 = ArcCos[V3.V4]-Th;
d35 = ArcCos[-V3.V5]-Th;

d45 = ArcCos[V4.V5]-Th;

(*sum of angular deviations that enters first order expansion term in energy *)
Sumd = d01 + d02 + d03 + d04 + d05 + d12 + d13 + d14 + d15 + d23 + d24 + d25 + d34 + d35 + d45;



(* second order expansion of angle, needed to perform harmonic analysis*)

Simplify[N[Series[Sumd, {t, 0,2}]]]

(* more nicely rewritten as a function *)

SumdE[x1_, x2_,x3_, x4_, x5_, y1_, y2_, y3_,y4_,y5_ ]:=0.5` ((x1^2)+(x2^2)+(x3^2)+(x4^2)+(x5^2))+0.08541019662496852`(x1 x2+ x2 x3+ x3 x4 + x4 x5 + x5 x1)-0.5854101966249684` (x3 x5 + x1 x3 + x2 x4 + x2 x5 + x1 x4)+0.4755282581475768`(- x1 y4+ x1 y3- x2 y5+ x2 y4- x3 y1+ x3 y5- x4 y2+ x4 y1- x5 y3+ x5 y2)+0.6249999999999999`( (y1^2)+ (y2^2)+ (y3^2)+ (y4^2)+ (y5^2))-1.0112712429686843`( y1 y2+y2 y3 + y3 y4 +y4 y5 + y5 y1)+0.3862712429686843`( y1 y3 + y2 y4 + y3 y5 + y4 y1 + y5 y2)+0.2938926261462365`(x1 y5-x1 y2+x2  y1-x2 y3 + x3 y2-x3 y4+ x4 y3- x4 y5+x5 y4  - x5 y1)

(* uniform displacements in both NS and EW direction does not generate quadratic terms *)

+0.6249999999999999`+0.3862712429686843`-1.0112712429686843`
0.5` +0.08541019662496852`-0.5854101966249684`

(* sum of squares of angles, enters 2nd order expansion*)

SumdSq = d01^2 + d02^2  + d03^2  + d04^2  + d05^2  + d12^2  + d13^2  + d14^2  + d15^2  + d23^2  + d24^2  + d25^2  + d34^2  + d35^2  + d45^2 ;

Simplify[N[Series[SumdSq, {t, 0,2}]]]

(2. (x1^2)+2.` (x2^2)+2.` (x3^2)-1.4472135954999579` x3 x4+2.` (x4^2)-0.552786404500042` x3 x5-1.4472135954999583` x4 x5+2.` (x5^2)+0.8506508083520399` x3 y1-0.8506508083520399` x4 y1-0.5257311121191335` x5 y1+2.5` (y1^2)+0.525731 x3 y2+0.8506508083520399` x4 y2-0.8506508083520398` x5 y2+0.19098300562505255` y1 y2+2.5` (y2^2)+0.` x3 y3+0.5257311121191335` x4 y3+0.8506508083520399` x5 y3+1.3090169943749475` y1 y3+0.1909830056250526` y2 y3+2.5` (y3^2)-0.5257311121191335` x3 y4-3.3306690738754696`*^-16 x4 y4+0.5257311121191337` x5 y4+1.3090169943749475` y1 y4+1.3090169943749472` y2 y4+0.19098300562505255` y3 y4+2.5` (y4^2)+x1 (-1.4472135954999579` x2-0.5527864045000421` x3-0.5527864045000421` x4-1.4472135954999579` x5+0.` y1-0.5257311121191335` y2-0.8506508083520401` y3+0.8506508083520401` y4+0.5257311121191335` y5)+x2 (-1.4472135954999583` x3-0.552786404500042` x4-0.552786404500042` x5+0.5257311121191335` y1-3.3306690738754696`*^-16 y2-0.5257311121191336` y3-0.8506508083520398` y4+0.8506508083520398` y5)-0.8506508083520398` x3 y5-0.5257311121191336` x4 y5+1.1102230246251565`*^-16 x5 y5+0.19098300562505255` y1 y5+1.3090169943749472` y2 y5+1.3090169943749472` y3 y5+0.1909830056250526` y4 y5+2.4999999999999996` (y5^2)) (((t+0.`))^2)+((O[t+0.`])^3)

(* rewritten more nicely as a function *)

SumdSqE[x1_, x2_,x3_, x4_, x5_, y1_, y2_, y3_,y4_,y5_ ]:=2((x1^2)+(x2^2)+(x3^2)+(x4^2)+(x5^2))-1.4472135954999579`(x1 x2+ x2 x3+ x3 x4 + x4 x5 + x5 x1)-0.552786404500042` (x3 x5 + x1 x3 + x2 x4 + x2 x5 + x1 x4)+0.5257311121191335`(x1 y5-x1 y2+x2  y1-x2 y3 + x3 y2-x3 y4+ x4 y3- x4 y5+x5 y4  - x5 y1)-0.8506508083520401`(- x1 y4+ x1 y3- x2 y5+ x2 y4- x3 y1+ x3 y5- x4 y2+ x4 y1- x5 y3+ x5 y2)+2.4999999999999996`( (y1^2)+ (y2^2)+ (y3^2)+ (y4^2)+ (y5^2))+0.1909830056250526`( y1 y2+y2 y3 + y3 y4 +y4 y5 + y5 y1)+1.3090169943749472`( y1 y3 + y2 y4 + y3 y5 + y4 y1 + y5 y2)

(*define a 10x10 matrix for quadratic form and diagonalize to determine the stable and unstable modes. V1 is the derivative of the two-body potential, it multiplies Sumd; V2 is the second derivative, it multiplies SumdSq/2. Actually here I dropped 2, so need to redefine the results appropriately, i.e. V2 has to be 2 times bigger, or V1 two times Smaller to see the critial behavior. Turns out that when V1 is -4/3 of V2 (i.e. first derivative is 2/3 or the second one) then there is an instbility, there are 5 zero modes, which correspond to costless distortions with all K vectors 2Pi k/5 *)

M = Table[0, {m,10}, {n,10}];
V1 = 1; V2 =2;
For[i=1, i<11, i++, 
For[j = 1, j < 11, j++,
ind1 = {0,0,0,0,0,0,0,0,0,0};
ind2= {0,0,0,0,0,0,0,0,0,0};
ind3 = {0,0,0,0,0,0,0,0,0,0};
If[j!=i,  ind1[[i]] = 1; ind1[[j]] = 1; ind2[[i]]=1; ind3[[j]]=1; M[[i,j]] = V1*( SumdE[ind1[[1]],ind1[[2]],ind1[[3]],ind1[[4]],ind1[[5]],ind1[[6]],ind1[[7]],ind1[[8]],ind1[[9]],ind1[[10]] ] - SumdE[ind2[[1]],ind2[[2]],ind2[[3]],ind2[[4]],ind2[[5]],ind2[[6]],ind2[[7]],ind2[[8]],ind2[[9]],ind2[[10]] ] - SumdE[ind3[[1]],ind3[[2]],ind3[[3]],ind3[[4]],ind3[[5]],ind3[[6]],ind3[[7]],ind3[[8]],ind3[[9]],ind3[[10]] ] )/2    
+V2*( SumdSqE[ind1[[1]],ind1[[2]],ind1[[3]],ind1[[4]],ind1[[5]],ind1[[6]],ind1[[7]],ind1[[8]],ind1[[9]],ind1[[10]] ] - SumdSqE[ind2[[1]],ind2[[2]],ind2[[3]],ind2[[4]],ind2[[5]],ind2[[6]],ind2[[7]],ind2[[8]],ind2[[9]],ind2[[10]] ] - SumdSqE[ind3[[1]],ind3[[2]],ind3[[3]],ind3[[4]],ind3[[5]],ind3[[6]],ind3[[7]],ind3[[8]],ind3[[9]],ind3[[10]] ] )/2    ];
If[j==i,  ind1[[i]] = 1; 
M[[i,i]] = V1*( SumdE[ind1[[1]],ind1[[2]],ind1[[3]],ind1[[4]],ind1[[5]],ind1[[6]],ind1[[7]],ind1[[8]],ind1[[9]],ind1[[10]] ])
+V2*( SumdSqE[ind1[[1]],ind1[[2]],ind1[[3]],ind1[[4]],ind1[[5]],ind1[[6]],ind1[[7]],ind1[[8]],ind1[[9]],ind1[[10]] ])];

]
]
Eigenvalues[M]

(*collecting terms in order to go into momentum space*)

(2V1+ .5V2)((x1^2)+(x2^2)+(x3^2)+(x4^2)+(x5^2))+
(-1.4472135954999579`V1 +0.08541019662496852`V2) (x1 x2+ x2 x3+ x3 x4 + x4 x5 + x5 x1)+
(-0.552786404500042` V1  -0.5854101966249684`V2)(x3 x5 + x1 x3 + x2 x4 + x2 x5 + x1 x4)+
(0.5257311121191335`V1 + 0.2938926261462365`V2)(x1 y5-x1 y2+x2  y1-x2 y3 + x3 y2-x3 y4+ x4 y3- x4 y5+x5 y4  - x5 y1)+(-0.8506508083520401`V1 + 0.4755282581475768`V2)(- x1 y4+ x1 y3- x2 y5+ x2 y4- x3 y1+ x3 y5- x4 y2+ x4 y1- x5 y3+ x5 y2)+
(2.4999999999999996`V1 + +0.6249999999999999`V2)( (y1^2)+ (y2^2)+ (y3^2)+ (y4^2)+ (y5^2))+
(0.1909830056250526`V1 -1.0112712429686843`V2) ( y1 y2+y2 y3 + y3 y4 +y4 y5 + y5 y1)+
(1.3090169943749472`V1 + 0.3862712429686843`V2)( y1 y3 + y2 y4 + y3 y5 + y4 y1 + y5 y2)

(* construction of 2x2 matrix in momentum space (x_q, y_q). K takes values 2Pik/5. This approach is in perfect agreement with 10x10. There is mode coupling between x and y but only at nonzero K. At critical value of V1/V2 = -4/3 there is a flat band at 0 energy *)

Nm = Table[0, {m,2}, {n,2}]
V1 = -2.; V2 = 1; K = -0 Pi/5;
Nm[[1,1]] = (2V2+ .5V1)+ (-1.4472135954999579`V2 +0.08541019662496852`V1)Cos[K] + (-0.552786404500042` V2  -0.5854101966249684`V1)Cos[2K];
Nm[[1,2]] = (0.5257311121191335`V2 + 0.2938926261462365`V1)(-I)  Sin[K] + (-0.8506508083520401`V2 + 0.4755282581475768`V1)(I) Sin[2K] 
Nm[[2,2]]= (2.4999999999999996`V2 + +0.6249999999999999`V1) + (0.1909830056250526`V2 -1.0112712429686843`V1)Cos[K] + (1.3090169943749472`V2 + 0.3862712429686843`V1)Cos[2K];
Nm[[2,1]] = -Nm[[1,2]];
Eigenvalues[Nm]