Vectorized molecular dynamics with Mathematica?

• To: mathgroup at smc.vnet.net
• Subject: [mg117267] Vectorized molecular dynamics with Mathematica?
• From: DmitryG <einschlag at gmail.com>
• Date: Sun, 13 Mar 2011 05:27:14 -0500 (EST)

```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 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 - 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 between
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]]],
]

aIni = 0.5 xConf;   vv0 = 0.1;
x0 = Table[RandomReal[{-aIni, aIni}], {i, 1, NPart}];
v0 = Table[RandomReal[{-vv0, vv0}], {i, 1, NPart}];

Vars = Join[Table[xx[i], {i, 1, NPart}], Table[vv[i], {i, 1, NPart}]];

Timing[
Sol = NDSolve[Join[Equations, IniConds], Vars, {t, 0, tMax}, MaxSteps
-> 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

```

• Prev by Date: Re: Contour line colors from z coord of a 3D plot
• Next by Date: Re: non linear system with 8 equations
• Previous by thread: Re: How to interactively re-scale graphics with mouse?
• Next by thread: Re: Vectorized molecular dynamics with Mathematica?