Re: Vectorized molecular dynamics with Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg117463] Re: Vectorized molecular dynamics with Mathematica?
- From: DmitryG <einschlag at gmail.com>
- Date: Sat, 19 Mar 2011 05:20:05 -0500 (EST)
- References: <ilq74c$bqd$1@smc.vnet.net>
On 16 Mrz., 07:32, Daniel Lichtblau <d... at wolfram.com> wrote:
> 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 part=
ic=
> > les
> >>> (screened attraction) *)
> >>> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces betwee=
n =
> > all
> >>> particles *)
> >>> (* Demonstrating pair interaction *)
> >>> xx = {x1, x2, x3, x4}; (* Position vector for a system of 4 part=
ic=
> > 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 tha=
t
> >>> 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. Thi=
s
> >>> 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 b=
et=
> > ween
> >>> two particles (screened
> >>> attraction) *)
> >>> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces betwee=
n =
> > 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}, Max=
St=
> > 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 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 *)
> > 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}, MaxSt=
eps
> > -> 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
Hi Daniel,
Thank you for the thoughts! What you suggest is quite reasonable and
is actually used in professional-grade MD simulations.
My purpose here is only to write a simple MD program for educational
purposes, and I am not doing MD in my research. I only wanted to
illustrate for students how Mathematica solves systems of many
equations. My existing 2d program can handle up to 50-70 particles on
a regular computer, and I thought with vectorizations I could do
several hundreds, maybe up to a thousand.
I am going to the APS Meeting in Dallas tomorrow and upon returning I
will look into Oliver's suggestions. But it is clear that the problem
of vectorization is not a easy one for MD.
Best,
Dmitry