Re: Vectorized molecular dynamics with Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg117384] Re: Vectorized molecular dynamics with Mathematica?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 16 Mar 2011 06:31:18 -0500 (EST)
DmitryG wrote:
> On 14 Mrz., 07:01, DrMajorBob <btre... at austin.rr.com> wrote:
>> vx is neither defined outside NDSolve nor used as a variable in NDSolve,
>> so the equation
>>
>> Derivative[1][x][t] == vx[t]
>>
>> is a problem.
>>
>> To make it clearer that F has vector arguments and can't be computed for
>> symbolic t, you might try:
>>
>> Clear[F1, F]
>> F1[{a_, b_, c_}, {d_, e_, f_}] =
>> Total[Outer[FInt, {a, b, c}, {d, e, f}]];
>> F[x_List, y_List] := F1[x, y]
>>
>> Bobby
>>
>>
>> On Sun, 13 Mar 2011 05:27:14 -0500, DmitryG <einsch... at gmail.com> wrote:
>>> Dear All,
>>> Recently with the help of Oliver Ruebenkoenig and Daniel Lichtblau I
>>> have discovered the advantage of writing systems of ODEs in a
>>> vectorized form that can be solved by NDSolve much faster and with a
>>> smaller usage of memory. It seems that probably every system of ODEs
>>> with individual equations following from the same general formula can
>>> be vectorized.
>>> To check this conjecture, I have tried to vectorize the molecular
>>> dynamics problem of many particles moving under the influence of one-
>>> and two-particle forces. My idea was to represent two-particle forces
>>> with the help of the Outer[..] comand. Below is my code that,
>>> unfortunately, does not work. Although the method can be formulated
>>> for motion in 2d and 3d, I illustrate the problem for a more
>>> transparent one-dimentional model:
>>> ***********************************************************************
>>> NPart = 3; (* 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 partic=
> les
>>> (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 partic=
> les
>>> *)
>>> 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 - vector form *)
>>> Equations = {
>>> m v'[t] == F[x[t], x[t]] + FConf[x[t]],
>>> x'[t] == vx[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];]
>>> *********************************************************
>>> Probably, in spite of my efforts, Mathematica does not understand that
>>> the argument x in Outer[FInt, x, x] is a list. What can be done to
>>> make this code work?
>>> Below is a non-vectorized version of the same problem. It is very
>>> similar to the above except of making the lists x and y explicit. This
>>> kills vectorization but the program works.
>>> **********************************************************
>>> NPart = 3; (* 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 (* Interaction force bet=
> ween
>>> two particles (screened
>>> attraction) *)
>>> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces between =
> all
>>> particles *)
>>> (* Confinement *)
>>> PowConf = 20; xConf = 30.;
>>> UConf[x_] := (x/xConf)^PowConf; (* Confinement potential *)
>>> FConf[x_] = -UConf'[x]; (* Confiment force *)
>>> (* Setting up and solving equations *)
>>> x[t_] := Table[xx[i][t], {i, 1, NPart}];
>>> v[t_] := Table[vv[i][t], {i, 1, NPart}];
>>> Equations = Join[
>>> Thread[m v'[t] == F[x[t], x[t]] + FConf[x[t]]],
>>> Thread[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 = Join[Thread[x[0] == x0], Thread[v[0] == v0]];
>>> Vars = Join[Table[xx[i], {i, 1, NPart}], Table[vv[i], {i, 1, NPart}]]=
> ;
>>> Timing[
>>> Sol = NDSolve[Join[Equations, IniConds], Vars, {t, 0, tMax}, MaxSt=
> eps
>>> -> Infinity];]
>>> xt[i_][t_] := xx[i][t] /. Sol[[1]]
>>> Plot[Evaluate[Table[xt[i][t], {i, 1, NPart}]], {t, 0, tMax}]
>>> ********************************************************
>>> Thank you for your help, as ever,
>>> Dmitry
>> --
>> DrMajor... at yahoo.com
>
> 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.
>
> 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];]
> xt[t_] := x[t] /. Sol[[1]]
> Plot[Table[xt[t][[i]], {i, 1, NPart}], {t, 0, tMax}]
> ***************************************************
>
You are using all pairs in your force computations. That scales
quadratically. To do better than this you might want some m
method-of-moments type of approach. At the least you'd need a way to
determine, quickly, which few points are close to a given point at any
step, and only use those for the updates to that one point.
My inclination here would be to compute Nearest[pts] at each step, and
use that to find who is within who's sphere of attraction, so to speak.
I am not sure how well this will play with Compile. It will certainly
require external evaluations. But it could reduce the complexity of
update steps to perhaps O(n log(n)).
Daniel Lichtblau
Wolfram Research