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]]],
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}, 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