Clear[f, diag]
f[k_, l_, r1_, r2_] := If[r1 == r2, {0, 0}, -k ((r1 - r2)/Norm[r1 - r2]*(Norm[r1 - r2] - l)) ](*r1 is the point feeling the force due to r2*)
diag[n_Integer] := Join[Table[i <-> i + n + 1, {i, Select[Range[1, n^2 - n], Mod[#, n] != 0 &]}], Table[i <-> i + n - 1, {i, Select[Range[2, n^2 - n + 1], Mod[#, n] != 1 &]}]] (*create a graph of only diagonals*)
l = 1; ld = Sqrt[2] l; k1 = 5; kd1 = k1/2; k2 = 1; kd2 = k2/2; n = 10;
pos = Flatten[Table[{x, y}, {x, 1, n, 1}, {y, 1, n, 1}], 1];
pos1 = pos; pos2 = pos;
m = AdjacencyMatrix[GridGraph[{n, n}]];
md = AdjacencyMatrix[Graph[Range[n^2], diag[n] ]];
dim = Dimensions[pos][[1]];
v = Partition[
   Flatten[Table[
     If[i > 2, {0, 0}, RandomReal[{-1, 1}, 2]], {i, 1, n}, {j, 1, n}]], 2];
v = Chop[(# - {Mean[v[[All, 1]] ], Mean[v[[All, 2]] ]}) & /@ v];
T = Total[Norm[#]^2 & /@ v];
v = 1 v/Sqrt[T];

acc1 = N@Table[
    Sum[ f[k1, l, pos[[i]], pos[[j]] ]*m[[i, j]] + f[kd1, ld, pos[[i]], pos[[j]] ]*md[[i, j]] , {j, 1, dim}] , {i, 1, dim}];
acc2 = N@Table[
    Sum[ f[k2, l, pos[[i]], pos[[j]] ]*m[[i, j]] + f[kd2, ld, pos[[i]], pos[[j]] ]*md[[i, j]] , {j, 1, dim}] , {i, 1, dim}];
dt = 0.03;
posold1 = pos1; pos1 = posold1 + v dt + acc1/2 dt^2;
posold2 = pos2; pos2 = posold2 + v dt + acc2/2 dt^2;
pos1[[1]] = pos[[1]];  pos1[[n]] = pos[[n]]; pos2[[1]] = pos[[1]]; pos2[[n]] = pos[[n]]; 

nstep = 500;
evo1 = Reap[Do[
      acc1 = Table[Sum[ f[k1, l, pos1[[i]], pos1[[j]] ]*m[[i, j]] + f[kd1, Sqrt[2] l, pos1[[i]], pos1[[j]] ]*md[[i, j]] , {j, 1, dim}] , {i, 1, dim}];
      posoldold1 = posold1;
      posold1 = pos1;
      pos1 = 2 posold1 - posoldold1 + acc1 dt^2;
      pos1[[1]] = pos[[1]]; pos1[[n]] = pos[[n]];
      Sow[pos1];
      , {nstep}];][[2, 1]];
Print["evo1 done!"];
evo2 = Reap[Do[
      acc2 = Table[Sum[ f[k2, l, pos2[[i]], pos2[[j]] ]*m[[i, j]] + f[kd2, Sqrt[2] l, pos2[[i]], pos2[[j]] ]*md[[i, j]] , {j, 1, dim}] , {i, 1, dim}];
      posoldold2 = posold2;
      posold2 = pos2;
      pos2 = 2 posold2 - posoldold2 + acc2 dt^2;
      pos2[[1]] = pos[[1]]; pos2[[n]] = pos[[n]];
      Sow[pos2];
      , {nstep}];][[2, 1]];

p0 = Table[
   Grid[{{
      Show[
       Graphics[{Gray, Disk[#, Offset[2]] & /@ pos }]
       ,
       AdjacencyGraph[m + md, VertexCoordinates -> evo1[[j]], VertexStyle -> Black, VertexShapeFunction -> Function[{pt, v, size}, Disk[pt, Offset[2]]], EdgeStyle -> Directive[Gray]]
        , PlotRange -> {{1 - 0.5, n + 0.5}, {1 - 0.5, n + 0.5}}, PlotLabel -> "Stiff crystal", LabelStyle -> {Black, Bold}, ImageSize -> 300]
      ,
      Show[
       Graphics[{Gray, Disk[#, Offset[2]] & /@ pos }]
       ,
       AdjacencyGraph[m + md, VertexCoordinates -> evo2[[j]], VertexStyle -> Black, VertexShapeFunction -> Function[{pt, v, size}, Disk[pt, Offset[2]]], EdgeStyle -> Directive[Gray]]
        , PlotRange -> {{1 - 0.5, n + 0.5}, {1 - 0.5, n + 0.5}}, PlotLabel -> "Pliable crystal", LabelStyle -> {Black, Bold}, ImageSize -> 300]
      }}]
   , {j, 2, nstep, 3}];
ListAnimate[p0]