Re: Vectorized molecular dynamics with Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg117382] Re: Vectorized molecular dynamics with Mathematica?
- From: Oliver Ruebenkoenig <ruebenko at wolfram.com>
- Date: Wed, 16 Mar 2011 06:30:57 -0500 (EST)
On Tue, 15 Mar 2011, DmitryG wrote:
[...]
Hi Dmitry,
>
> Thank you, Bobby!
>
> There is a mistake in the vectorized version, I write vx instead of v
> (an artefact of rewriting 2d version of the problem to a 1d version).
> Upon correction the code works and the results are correct. However,
> the execution time in the vectorized version increases dramatically
> with the number of particles NPart and largely exceeds the execution
> time of the non-vectorized version. This is not what we expect from
> vectorization;-)). Try NPart=5 and you'll see the effect.
>
> Can anybody suggest how to make the vectorized version fast? Maybe my
> initial idea to use Outer[..] was not good and there are better
> programming means? Below is the corrected and working vectorized
> version.
>
I think, the idea was quite right. There, is however, a caveat, for Outer.
e.g.:
Compile[{{xIn, _Real, 1}, {yIn, _Real, 1}}, Outer[Sqrt, xIn, yIn]]
The message say, that Outer can only be compiled for a function argument
that is either Time, Plus or List and
Function[{x, y}, (-2*(-x + y))/E^(-x + y)^2]
can thus not be compiled. Outer is tricky function to compile, one reason
is that the function argument can contract the tensor and the compiler can
then not figure out what rank the result to be returned needs to have.
Now, you have several options.
> Thank you for your help,
>
> Dmitry
>
> *************************************************
> NPart = 5; (* Number of particles *)
> tMax = 1000; (* Maximal time of the calculation *)
> m = 1; (* Mass of a particle *)
>
> (* Pair interaction *)
> UInt[x_] = -Exp[-x^2]; (* Pair interaction energy *)
> FInt[x_, y_] = -UInt'[x] /. x -> y - x (* Force between two particles
> (screened attraction) *)
>
>
> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces between all
> particles *)
>
> (* Demonstrating pair interaction *)
> xx = {x1, x2, x3, x4}; (* Position vector for a system of 4 particles
> *)
> MatrixForm[Total[Outer[FInt, xx, xx]]] (* Vector of forces in the
> RHS of Newton's equations *)
>
> (* Confinement *)
> PowConf = 20; xConf = 30.;
> UConf[x_] := (x/xConf)^PowConf; (* Confinement potential *)
> FConf[x_] = -UConf'[x] (* Confinement force *)
>
> (* Setting up and solving equations *)
> Equations = {
> m v'[t] == F[x[t], x[t]] + FConf[x[t]],
> x'[t] == v[t]
> }
>
> aIni = 0.5 xConf; vv0 = 0.1;
> x0 = Table[RandomReal[{-aIni, aIni}], {i, 1, NPart}];
> v0 = Table[RandomReal[{-vv0, vv0}], {i, 1, NPart}];
>
> IniConds = {x[0] == x0, v[0] == v0};
>
> Vars = {x, v};
>
> Timing[
> Sol = NDSolve[Join[Equations, IniConds], Vars, {t, 0, tMax}, MaxSteps
> -> Infinity];]
(*
Because Outer[f,...] can not be compiled, NDSolve uses the uncompiled
version. Nevertheless, you can for NDSolve to compile what ever is
compilable with the option
Sol = NDSolve[Join[Equations, IniConds], Vars, {t, 0, tMax}, MaxSteps
-> Infinity,
Compiled->True
];]
this will give a small speedup.
*)
> xt[t_] := x[t] /. Sol[[1]]
> Plot[Table[xt[t][[i]], {i, 1, NPart}], {t, 0, tMax}]
> ***************************************************
>
>
Another, better approach is to avoid Outer and compile the function to C
and use that.
(* Please double check my arithmetic *)
cf = Compile[{{x, _Real, 1}, {yIn, _Real, 1}},
Total[(Function[y, -2 E^-(x + y)^2 (x + y)] /@ -yIn)],
CompilationTarget -> C, RuntimeOptions -> "Speed"]
FFF[x_List, y_List] := cf[x, y]
Equations = {m v'[t] == FFF[x[t], x[t]] + FConf[x[t]], x'[t] == v[t]}
For 10 particles this runs in about 4.5 sec on my laptop.
Depending on what you need, you could also play with
, AccuracyGoal -> 4, PrecisionGoal -> 4
To give a further speedup.
What, I'd try, is the following. Perhaps you can reformulate your function
Total[Outer[FInt, xx, xx]]
in terms of dot products. I could imagine that to be efficient too, Also,
note that the
Total[Outer[FInt, xx, xx]]
is almost symmetric, perhaps you could exploit that too.
hth,
Oliver