Re: Vectorized molecular dynamics with Mathematica?
- To: mathgroup at smc.vnet.net
- Subject: [mg117296] Re: Vectorized molecular dynamics with Mathematica?
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Mon, 14 Mar 2011 06:00:11 -0500 (EST)
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 <einschlag 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 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
>
--
DrMajorBob at yahoo.com