MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Joining points of ListPlot
  • Next by Date: Re: Joining points of ListPlot
  • Previous by thread: Vectorized molecular dynamics with Mathematica?
  • Next by thread: Re: Vectorized molecular dynamics with Mathematica?