**Bug introduced in 9.0 or earlier and persisting through 10.4.1 or later**

----

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