Re: Vectorized molecular dynamics with Mathematica?

*To*: mathgroup at smc.vnet.net*Subject*: [mg117463] Re: Vectorized molecular dynamics with Mathematica?*From*: DmitryG <einschlag at gmail.com>*Date*: Sat, 19 Mar 2011 05:20:05 -0500 (EST)*References*: <ilq74c$bqd$1@smc.vnet.net>

On 16 Mrz., 07:32, Daniel Lichtblau <d... at wolfram.com> wrote: > DmitryG wrote: > > 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 part= ic= > > les > >>> (screened attraction) *) > >>> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces betwee= n = > > all > >>> particles *) > >>> (* Demonstrating pair interaction *) > >>> xx = {x1, x2, x3, x4}; (* Position vector for a system of 4 part= ic= > > 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 tha= t > >>> 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. Thi= s > >>> 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 b= et= > > ween > >>> two particles (screened > >>> attraction) *) > >>> F[x_List, y_List] := Total[Outer[FInt, x, y]]; (* Forces betwee= n = > > 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}, Max= St= > > 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 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 *) > > 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}, MaxSt= eps > > -> Infinity];] > > xt[t_] := x[t] /. Sol[[1]] > > Plot[Table[xt[t][[i]], {i, 1, NPart}], {t, 0, tMax}] > > *************************************************** > > You are using all pairs in your force computations. That scales > quadratically. To do better than this you might want some m > method-of-moments type of approach. At the least you'd need a way to > determine, quickly, which few points are close to a given point at any > step, and only use those for the updates to that one point. > > My inclination here would be to compute Nearest[pts] at each step, and > use that to find who is within who's sphere of attraction, so to speak. > I am not sure how well this will play with Compile. It will certainly > require external evaluations. But it could reduce the complexity of > update steps to perhaps O(n log(n)). > > Daniel Lichtblau > Wolfram Research Hi Daniel, Thank you for the thoughts! What you suggest is quite reasonable and is actually used in professional-grade MD simulations. My purpose here is only to write a simple MD program for educational purposes, and I am not doing MD in my research. I only wanted to illustrate for students how Mathematica solves systems of many equations. My existing 2d program can handle up to 50-70 particles on a regular computer, and I thought with vectorizations I could do several hundreds, maybe up to a thousand. I am going to the APS Meeting in Dallas tomorrow and upon returning I will look into Oliver's suggestions. But it is clear that the problem of vectorization is not a easy one for MD. Best, Dmitry

**Re: Does MemoryConstrained count stack space?**

**Re: How to kill slave kernel securely?**

**Re: Vectorized molecular dynamics with Mathematica?**

**Export strings to .pdf**