Re: Vectorized molecular dynamics with Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg117344] Re: Vectorized molecular dynamics with Mathematica?
- From: DmitryG <einschlag at gmail.com>
- Date: Tue, 15 Mar 2011 06:08:08 -0500 (EST)
- References: <ilksiu$6f7$1@smc.vnet.net>
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}]
***************************************************