I tried to solve Hamiltonian system ($Q$ is a vector of all generalized coordinates, $P$ - of generalized momentum)
$$
\frac{\mathrm{d} Q}{\mathrm{d} t}=\frac{\partial H}{\partial P} \\
\frac{\mathrm{d} P}{\mathrm{d} t}=-\frac{\partial H}{\partial Q}
$$
where
$$
H=\frac{1}{2} \sum_{i=1}^{n}\overrightarrow{p}_{i}^{2} + \frac{k}{2} \sum_{i=2}^{n}(\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|-rest)^2 + \frac{\kappa}{2} \sum_{i=2}^{n-1}\arccos^2\frac{(\overrightarrow{q}_{i}-\overrightarrow{q}_{i-1})\cdot (\overrightarrow{q}_{i+1}-\overrightarrow{q}_{i})}{\left \| \overrightarrow{q}_{i}-\overrightarrow{q}_{i-1} \right \|\left \| \overrightarrow{q}_{i+1}-\overrightarrow{q}_{i} \right \|}
$$
with "SymplecticPartitionedRungeKutta" method numerically using Mathematica (description could be found here https://en.wikipedia.org/wiki/Symplectic_integrator http://www.sciencedirect.com/science/article/pii/S0377042701004927)

The code is below

 ϰ = 20;
 k = 20;
 rest = Sqrt[5];
 n = 3;
 dim = 3;
 Q = Table[Table[Subscript[q, j][i][t], {j, 1, dim}], {i, 1, n}];
 P = Table[Table[Subscript[p, j][i][t], {j, 1, dim}], {i, 1, n}];
 icsQ = {{-1, -1, 0}, {0, 1, 0}, {1, -1, 0}};
 icsP = Table[Table[0, {j, 1, dim}], {i, 1, n}];
 
 H = 1/2 (Sum[P[[i]].P[[i]], {i, 1, n}] + 
 Sum[(Norm[Q[[i]] - Q[[i - 1]]] - rest)^2, {i, 2, n}]*k + 
 Sum[(VectorAngle[Q[[i]] - Q[[i - 1]], 
 Q[[i + 1]] - Q[[i]]])^2, {i, 2, n - 1}]*ϰ);
 Q' = Flatten[D[#, t] & /@ Q];
 P' = Flatten[D[#, t] & /@ P];

 Conjugate'[y_[x_][t_]] := 1;
 Abs'[x_] := Sign[x];

 var = Flatten[Q~Join~P];
 ics = Flatten[icsQ~Join~icsP];
 eqn = Table[Q'[[i]] == D[H, var[[n*dim + i]]], {i, 1, n*dim}] ~Join~
 Table[P'[[i]] == -D[H, var[[i]]], {i, 1, n*dim}] ~Join~
 Table[(var[[i]] /. t -> 0) == ics[[i]], {i, Length[var]}];

 method = {"SymplecticPartitionedRungeKutta", "PositionVariables" -> var[[1 ;; dim*n]]};
 sol = NDSolve[eqn, var, {t, 0, 10}, 
 Method -> method, WorkingPrecision -> 10, 
 MaxStepSize -> 0.001, MaxSteps -> ∞]

However, error pops out: 

> NDSolve::sprksep: "The Hamiltonian of the differential system in the
> method does not appear to be in separable form. Try using the method
> ImplicitRungeKutta with coefficients
> ImplicitRungeKuttaGaussCoefficients."

But it's easy to see that the Hamiltonian of the system definately has separate form H(p,q)=T(p)+V(q).

Can anyone help?

P.S. In order to prevent questions. You could notice that I defined two "strange" functions: 

 Conjugate'[y_[x_][t_]] := 1;
 Abs'[x_] := Sign[x];

Idea: Mathematica doesn't take the derivative from non-analytical functions. And when you take derivative in form `D[H,q[t]]`, terms of type `D[Conjugate[q[t]],q[t]]` appear. If I don't fix it, I will get an error:

> NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0