del[f_] := {D[f,bx],D[f,by],D[f,bz]}
equilibrium = { bx -> 0 , by -> 0, bz -> 0}
Vext = a(E^(2 bx) + E^(2 by) + E^(2 bz) )
Tad = b E^(-2 (bx+by+bz)/3 )
U = (xi -1 ) Tad
ForceA = D[ Tad + U +Vext, bx] /. equilibrium
bvalue = Simplify[ Solve[ForceA == 0 , b]][[1]]
Ead = Simplify[ Tad +U +Vext /. bvalue  ]
KA = del[del[ Ead ]] /. equilibrium
ihalf = { {1,0,0},{0,1,0},{0,0,lam}}
mat = ihalf.KA.ihalf / ( a 2)
omega2A = Eigenvalues[mat]
Tdia = b (E^(-2 bx) + E^(-2 by) +E^( -2 bz) ) /3
Edia = Tdia + U + Vext /. bvalue
Kdia = del[ del [ Edia ]] /. equilibrium
mat = Simplify[ ihalf.Kdia.ihalf / ( 2 a) ]
omega2dia = Eigenvalues[mat]




