Group Abstract Group Abstract

Message Boards Message Boards

1
|
6.6K Views
|
8 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Using NBodySimulation with a custom pairwise function?

Madeleine Sutherland, MIT
Posted 4 years ago

Background:

Popular solid-state Nuclear Magnetic Resonance (NMR) techniques measure inter-atomic distances by observing the rate of transfer of magnetization between two spins by a spin diffusion mechanism. Plots of normalized magnetization vs spin diffusion mixing time are matched with simulated spin diffusion buildup curves for various inter-atomic distances (r) to find the best fit.

For some interesting biomolecular problems, such as the distance from membrane lipids chains to a site on a protein, the following 1D lattice simulation procedure is implemented:

  1. Arrange (r+1) points i along a 1D grid with width equal to your "guess" distance r.
  2. Set initial magnetization: M1,t=0 = 1 and Mi=/=1, t=0 = 0
  3. Evolve magnetization for all points by ∆Mi/∆t = D/a2[(Mi+1 - Mi) - (Mi - Mi-1)] (where D and a2 are constants) for some # of steps.
  4. Repeat for different simulation times and make a buildup curve; determine fit with data.

Question How do I implement Step 3 above in the law argument of the NBodySimulation function? Follow-up question if this is impossible: Should I be going a SystemModelParametricSimulate direction instead?

Note: this Mathematica Stack Exchange thread ends with a similar unanswered query: https://mathematica.stackexchange.com/questions/105342/how-best-to-simulate-n-body-systems-in-a-functional-way
POSTED BY: Madeleine Sutherland
8 Replies
Sort By:
POSTED BY: Madeleine Sutherland

Hello Madeleine,

thank you for you thank you. That is kind.

Concerning the diffusion in two phases : Long time ago I wrote a code dealing with that problem. Perhaps it is helpful.

a = 2.; b = 6; du = .1; dv = .8; u0 = 5.; v0 = 1.5;ux = 1/a; vx = 1/(a - b);

Clear[sol]
sol[t1_] := Module[{},
  {g, h} /. NDSolve[
     {
      D[ g[x, t], t] == ux^2 du  D[g[x, t], x, x],
      D[ h[x, t], t] == vx^2 dv D[ h[x, t], x, x],
      g[x, 0] == u0,
      h[x, 0] == v0,
      (D[ g[x, t], x] == 0) /. x -> 0,
      (D[ h[x, t], x] == 0) /. x -> 0,
      (  ux  du D[ g[x, t], x] == vx  dv D[h[x, t], x]) /. x -> 1,
      h[1, t] == g[1, t] (1 - Exp[-30 t]) + v0 Exp[-30 t]
      },
     {g, h},
     {x, 0, 1}, {t, 0, t1}][[1]]
  ]

lsg = sol[20];
f1 = lsg[[1]];
f2 = lsg[[2]];

fp[x_, t_] := If[x < a, f1[x/a, t], f2[(x - b)/(a - b), t]]

Manipulate[
 Plot[fp[x, tt], {x, 0, b}, PlotRange -> {0, 6}],
 {tt, 0, 20}]

And good luck for your thesis.
POSTED BY: Hans Dolhaine
POSTED BY: Henrik Schachner

Hello Joshua,

that looks intriguing, but: the equation given by Madeline (point 3 above) is a diffusion equation, so Total/@results should always give one. This at least is the case with mm[i] calculated in my code above. What is going wrong?
POSTED BY: Hans Dolhaine

I suspect it might be related to the boundary conditions...there is a non-conservation at the boundaries because of the interpretation of how to treat the (Mi-1 - Mi) term when i=1. (this can be seen by doing a MatrixForm@update[10, 0.1]). My reading of @MadeleineSutherland's post was that the i-1 term can be discarded, while retaining the i term...but maybe that's not such a physically appropriate idea. Modifying the ends of the matrix to remove these terms should suffice.
POSTED BY: Joshua Schrier
POSTED BY: Hans Dolhaine

To answer the original question: the documentation for NBodySimulation seems to suggest that only position changes are allowed, which would preclude using it for a problem like this.

You could also write this in a functional style (avoiding the use of Do loops) like the following. This is equivalent to @HansDolhaine's first solution above:

    initialState[r_Integer] := {1}~Join~ConstantArray[0, r] 

(*  define a matrix to perform the update step; for large systems, consider using SparseArray*)
(* lump all constants together*)
    update[r_Integer, constants_] := constants*Table[
       (KroneckerDelta[i + 1, j] - KroneckerDelta[i, j]) - (KroneckerDelta[i, j] -KroneckerDelta[i - 1, j]),
     {i, r + 1}, {j, r + 1}]

(*compute a lists of states at each tilmestep*)
results = With[
   {start = initialState[10],
    step = update[10, 0.1],
    timesteps = 20},
   NestList[ (step . # + #) &, start, timesteps]];

ArrayPlot[results] (*watch the magnitization*)

ListLinePlot[Mean /@ results] (*compute properties at each timestep and plot*)
POSTED BY: Joshua Schrier
POSTED BY: Hans Dolhaine
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard