# NBodySimulation: custom pairwise potential

I'm trying to define a Lennard-Jones potential in 2D for NBodySimulation

\[Epsilon] = 1; rm = 1;

LJ[b1_, b2_] :=
With[{rd = rm/EuclideanDistance[b1["Position"], b2["Position"]]},
\[Epsilon]*(rd^12 - 2*rd^6)
]

data = NBodySimulation[
<|
"PairwisePotential" -> LJ
|>,
{
<|"Mass" -> 1, "Position" -> {-0.6, 1}, "Velocity" -> {0., 0}|>,
<|"Mass" -> 2, "Position" -> {0.6, -1}, "Velocity" -> {0., 0}|>
},
10
];


This gives the following error:

NDSolveProcessEquations::ndnum: Encountered non-numerical value for a derivative at t == 0..


I've seen similar questions on this site relating to NDSolve. Accord to these to issue could be either with the initial conditions, or the function given to NDSolve.

Here, the issue is clearly not with the initial conditions, as replacing LJ with "Harmonic" works perfectly:

data = NBodySimulation[
<|
"PairwisePotential" -> "Harmonic"
|>,
{
<|"Mass" -> 1, "Position" -> {-0.6, 1}, "Velocity" -> {0., 0}|>,
<|"Mass" -> 2, "Position" -> {0.6, -1},
"Velocity" -> {0., 0}|>
},
10
];


So maybe NDSolve doesn't like EuclideanDistance. However when I calculate the distance manually, I still run into error messages:

LJ[b1_, b2_] :=
With[{rd =
1/((b1["Position"][[1]] -
b2["Position"][[1]])^2 + (b1["Position"][[2]] -
b2["Position"][[2]])^2)},
\[Epsilon]*(rd^6 - 2*rd^3)
]


In this case NBodySimulation[...] gives the following error:

Part::partw: Part 2 of q1[t] does not exist.
Part::partw: Part 2 of q2[t] does not exist.


Based on this I think that b1["Position"] passed to my potential function is not really a vector, or is not evaluated when I'm trying to reference its components. But then how possibly can I implement a 2D or 3D potential for NBodySimulation?

• FWIW: Using Norm to compute the distance in the denominator of rd works fine...i.e.: rd = rm/Norm[b1["Position"] - b2["Position"]] Apr 14 '20 at 19:44
• @JoshuaSchrier You're correct, thank you so much! Bit off-topic, but could you explain why Norm works and my other two versions don't, more specifically, why I can't reference components of the position vectors? Apr 14 '20 at 19:59
• Try replacing your function with LJ[b1_, b2_] := With[{rd = rm/Norm[b1["Position"] - b2["Position"]]}, Print[b1]; \[Epsilon]*(rd^12 - 2*rd^6)] You'll notice that b1 is a Symbol, not a vector. Although I'm not sure why EuclideanDistance should fail to work in that case. Apr 14 '20 at 20:06