28
$\begingroup$

I have the following function

$$V(r) = \sum_{i=1}^N 4 \epsilon_i \left(\frac{\sigma_i^{12}}{\|r-r_{0i}\|^{12}}-\frac{\sigma_i^6}{\|r-r_{0i}\|^6}\right)$$

which -for those interested- corresponds to a sum of Lennard-Jones potentials, with the following real life set of parameters

sig = {0.329633, 0.0400014, 0.405359, 0.235197, 0.387541, 0.235197,
       0.235197, 0.387541, 0.235197, 0.235197, 0.387541, 0.235197,
       0.235197, 0.387541, 0.235197, 0.235197, 0.329633, 0.0400014,
       0.0400014, 0.0400014, 0.356359, 0.302906, 0.329633, 0.0400014,
       0.387541, 0.235197, 0.235197, 0.356359, 0.302906, 0.329633,
       0.0400014, 0.405359, 0.235197, 0.387541, 0.235197, 0.235197,
       0.405359, 0.235197, 0.36705, 0.235197, 0.235197, 0.235197, 0.36705,
       0.235197, 0.235197, 0.235197, 0.356359, 0.302906, 0.329633,
       0.0400014, 0.405359, 0.235197, 0.387541, 0.235197, 0.235197,
       0.387541, 0.235197, 0.235197, 0.356359, 0.302906, 0.329633,
       0.0400014, 0.0400014, 0.356359, 0.302906, 0.329633, 0.0400014,
       0.405359, 0.235197, 0.405359, 0.235197, 0.36705, 0.235197, 0.235197,
       0.235197, 0.387541, 0.235197, 0.235197, 0.36705, 0.235197, 0.235197,
       0.235197, 0.356359, 0.302906}

eps = {0.8368, 0.192464, 0.08368, 0.092048, 0.23012, 0.092048, 0.092048,
       0.23012, 0.092048, 0.092048, 0.23012, 0.092048, 0.092048, 0.23012,
       0.092048, 0.092048, 0.8368, 0.192464, 0.192464, 0.192464, 0.46024,
       0.50208, 0.8368, 0.192464, 0.23012, 0.092048, 0.092048, 0.46024,
       0.50208, 0.8368, 0.192464, 0.08368, 0.092048, 0.23012, 0.092048,
       0.092048, 0.08368, 0.092048, 0.33472, 0.092048, 0.092048, 0.092048,
       0.33472, 0.092048, 0.092048, 0.092048, 0.46024, 0.50208, 0.8368,
       0.192464, 0.08368, 0.092048, 0.23012, 0.092048, 0.092048, 0.23012,
       0.092048, 0.092048, 0.29288, 0.50208, 0.8368, 0.192464, 0.192464,
       0.46024, 0.50208, 0.8368, 0.192464, 0.08368, 0.092048, 0.08368,
       0.092048, 0.33472, 0.092048, 0.092048, 0.092048, 0.23012, 0.092048,
       0.092048, 0.33472, 0.092048, 0.092048, 0.092048, 0.46024, 0.50208}

r0 = {{0.681, -2.673}, {0.605, -2.736}, {0.715, -2.578}, {0.812, -2.583},
      {0.628, -2.607}, {0.654, -2.698}, {0.533, -2.609}, {0.63, -2.515},
      {0.559, -2.545}, {0.609, -2.423}, {0.763, -2.509}, {0.825, -2.446},
      {0.804, -2.6}, {0.742, -2.461}, {0.709, -2.367}, {0.829, -2.465},
      {0.642, -2.547}, {0.629, -2.515}, {0.675, -2.642}, {0.555, -2.543},
      {0.693, -2.445}, {0.778, -2.359}, {0.585, -2.422}, {0.526, -2.499},
      {0.55, -2.291}, {0.543, -2.236}, {0.459, -2.299}, {0.638, -2.219},
      {0.661, -2.105}, {0.715, -2.293}, {0.688, -2.387}, {0.829, -2.244},
      {0.791, -2.156}, {0.867, -2.346}, {0.893, -2.431}, {0.946, -2.309},
      {0.754, -2.38}, {0.674, -2.422}, {0.697, -2.255}, {0.627, -2.281},
      {0.77, -2.205}, {0.657, -2.197}, {0.802, -2.482}, {0.729, -2.503},
      {0.825, -2.566}, {0.883, -2.448}, {0.956, -2.216}, {1.034, -2.126},
      {0.986, -2.287}, {0.921, -2.356}, {1.104, -2.268}, {1.178, -2.271},
      {1.13, -2.381}, {1.056, -2.379}, {1.217, -2.362}, {1.135, -2.523},
      {1.218, -2.53}, {1.055, -2.534}, {1.137, -2.64}, {1.171, -2.628},
      {1.083, -2.751}, {1.043, -2.753}, {1.08, -2.832}, {1.099, -2.133},
      {1.196, -2.059}, {0.987, -2.095}, {0.911, -2.161}, {0.964, -1.963},
      {1.042, -1.957}, {0.833, -1.955}, {0.818, -1.859}, {0.844, -2.053},
      {0.759, -2.051}, {0.92, -2.024}, {0.861, -2.146}, {0.714, -1.996},
      {0.729, -2.089}, {0.708, -1.934}, {0.581, -1.991}, {0.505, -2.018},
      {0.565, -1.898}, {0.586, -2.053}, {0.98, -1.845}, {1.047, -1.739}}

which -again, for those interested- correspond to a 2D projection of the first five residues of a S4S5 alpha helix for the Kv1.2 Ion Channel.

Defining $V(r)$ as

v[r_, r0_, s_, ep_] := 4 ep (s^12/EuclideanDistance[r, r0]^12 -
                             s^6/EuclideanDistance[r, r0]^6);

I can plot the potential with

Plot3D[Sum[v[{x, y}, r0[[i]], sig[[i]], eps[[i]]], {i, 1, 84}], 
       {x, -0.5, 2}, {y, -3.5, -1}, 
       PlotStyle -> Directive[Opacity[0.35], Blue], 
       AxesLabel -> {x, y}, PlotRange -> {-5, 1}
]

obtaining the following output

plot of sum of potentials

or, using PlotPoints -> 50

refined plot of sum of potentials

where the minima/maxima can be seen really well.

The thing is, I have a lot of these objects, with a lot more elements, and a lot of minima/maxima, in a way that is very expensive for my (old) computer to simply increase PlotPoints for smoother graphics, and I was wondering, due the fact that $V(r)$ rapidly decreases, if there is a way to ask MMA to increase resolution near the minima/maxima, and reduce it far from the data set r0.

Hope my question is clear and interesting.

--FINAL EDIT--

First of all, I want to thank you all for your comments and answers. If I could, I'd accept all of them, since each one gave me the insight needed to solve my problem. For obvious reasons, Silvia's answer deserves maximum recognition, but readers should also check PlatoManiac and Sjoerd C. de Vries responses, as they will provide a full picture.

Now, to honour the work of all the people involved, I'll show you the beautiful application of the code they've worked out.

Here is a caricaturization of the S4S5 linker helix, believed to play a major role in the opening and closing of voltage gated potassium channels.

S4S5 linker helix

This helix generates all sort of van der Waals interactions, that can be modelled by the Lennard-Jones potential.

For specific reasons, I need to see how this potential looks on a given plane, namely the XY plane, and it is very important to capture all the maxima and minima of it, as it will provide a full picture of the dynamics in that plane.

Thanks to Silvia's code, one can see:

a view from below of the potential generated by the helix, View from below

a sideways view, where the peaks are actually minima, Sideways veiw, the peaks are actually minima

and a view form above,

View from above

What you are seeing is the van der Waals interactions generated by the helix over that plane, and if you're wondering what does that sharp barrier around the helix is, it's the macroscopic consequence of Pauli's Exclusion Principle!

Thank you all for your help, you've put a big smile on my face!

Improve this question
$\endgroup$
6
  • $\begingroup$ You might already get some speed improvement by defining expr and then Plot3D[expr,...]. $\endgroup$
    – b.gates.you.know.what
    Sep 10, 2012 at 20:27
  • $\begingroup$ @b.gatessucks You're right. In my code it's done that way. I guess I didn't gave it importance due that the major delay comes from creating the plot and not from the evaluation of the sum. $\endgroup$
    – Pragabhava
    Sep 10, 2012 at 20:31
  • $\begingroup$ I think generating a custom grid and plotting discrete data by ListPlot3D will be good. $\endgroup$
    – Silvia
    Sep 11, 2012 at 17:23
  • 2
    $\begingroup$ You can squeeze out a little more speed by defining v using Dot instead of EuclideanDistance to avoid all the Abs : v[r_,r0_,s_,ep_]:=4ep(s^12/(#.#&[r-r0])^6-s^6/(#.#&[r-r0])^3). It's also worth putting expr=Expand@Sum[...] $\endgroup$
    – Simon Woods
    Sep 11, 2012 at 19:48
  • 1
    $\begingroup$ @Pragabhava, Expand here just reduces the number of multiplications in the expression, for a minor speed increase. In your definition of vTot you have used SetDelayed - this is unwise as it means the Sum will be recomputed every time you evaluate vTot. $\endgroup$
    – Simon Woods
    Sep 12, 2012 at 19:10

3 Answers 3

Reset to default
33
$\begingroup$

In order to emphasize the wanted features, my idea is to manually generate a grid which is fine near the ridge line and sketchy far away:

Mathematica graphics

1:

First we define our functions following Simon's suggestion:

v[r_, r0_, s_, ep_] := 4*ep*(s^12/(#.#&)[r-r0]^(12/2)-s^6/(#.#&)[r-r0]^(6/2))

expr = Sum[v[{x, y}, r0[[i]], sig[[i]], eps[[i]]], {i, 1, 84}]; 

func = Compile[{x, y}, Evaluate[expr // Expand]];

funcWrap[x_?NumericQ, y_?NumericQ] := func[x, y]

2:

The shape of the plot makes me feel that it's more convenient to work in a polar coordinate system in the $x$-$y$ plane, with rC (center of r0) as the origin. In order to generate a grid fit for expr, two key lines need to be known first: one is the boundary innerBoundLine where expr==2 (because we want to cut the PlotRange below $1<2$), and the other is the ridge line minPointsPolar, highlighted blue in the plot above:

rC = Mean[r0];

innerBoundLine = ContourPlot[funcWrap[x, y] == 2,
   {x, -0.5, 2}, {y, -3.5, -1}, PlotPoints -> 50];

innerBoundPoints = Cases[
     Normal[innerBoundLine[[1]] /.Tooltip[expr_, _] :> expr],
     Line[pts_] :> pts, ∞][[1]] // Most;

innerBoundPointsPolar =
  {Arg[{1, I}.#], Norm@#} &[# - rC] & /@ innerBoundPoints // SortBy[#, First] &;

innerBoundFunc = Interpolation[
   ReplacePart[#,{{1, 2}, {-1, 2}} -> Mean[#[[{1, -1}, 2]]]] &[innerBoundPointsPolar],
   PeriodicInterpolation -> True, InterpolationOrder -> 1];

minPointsPolar = Module[
   {ρmin, linefunc, ρinit},
   Table[ρinit = innerBoundFunc[φ];
    ρmin = ρ /. FindMinimum[ funcWrap @@ (ρ {Cos[φ], Sin[φ]} + rC),
        {ρ, ρinit}, Method -> "PrincipalAxis"][[2]];
    {φ, ρmin},
    {φ, 0, 2 π, 1. Degree}]];

minFunc = Interpolation[
   ReplacePart[#, {{1, 2}, {-1, 2}} -> Mean[#[[{1, -1}, 2]]]] &[minPointsPolar],
   PeriodicInterpolation -> True, InterpolationOrder -> 1];

3:

Now we can proceed to the grid generation step. Here, four grids with different fineness are generated. minLineGrid is near the ridge line and is the finest. From fineGrid to transitionalGrid to outlineGrid, the grids are farther and farther from the ridge, and sketchier and sketchier.

minLineGrid = Module[{ρmin, linefunc},
   Table[
     ρmin = minFunc[φ];
     linefunc = Append[#,
         funcWrap @@ #] &[
       ρ {Cos[φ], Sin[φ]} + rC];
     Table[linefunc, {ρ,
       Range[.98, 1.02, .01] ρmin
       }],
     {φ, 0, 2 π, 1. Degree}] // Flatten[#, 1] &
   ];

fineGrid = Module[{ρmin, linefunc, ρinit},
   Table[
     ρinit = innerBoundFunc[φ];
     ρmin = minFunc[φ];
     linefunc = Append[#,
         funcWrap @@ #] &[
       ρ {Cos[φ], Sin[φ]} + rC];
     Table[linefunc, {ρ, Join[
        Rescale[Range[0, 1, .5], {0, 1}, {ρinit, ρmin}],
        Rescale[Range[0, 1, .3], {0, 1}, {ρmin, 1.1 ρmin}]
        ]}],
     {φ, 0, 2 π, 2. Degree}] // Flatten[#, 1] &
   ];

transitionalGrid = Module[{ρmin, linefunc},
   Table[
     ρmin = minFunc[φ];
     linefunc = Append[#,
         funcWrap @@ #] &[
       ρ {Cos[φ], Sin[φ]} + rC];
     Table[linefunc, {ρ,
       Rescale[Range[0, 1, .2], {0, 1}, {1.1 ρmin, 1.5 ρmin}]
       }],
     {φ, 0, 2 π, 5. Degree}] // Flatten[#, 1] &
   ];

outlineGrid = Module[{ρmin, linefunc},
   Table[
     ρmin = minFunc[φ];
     linefunc = Append[#,
         funcWrap @@ #] &[
       ρ {Cos[φ], Sin[φ]} + rC];
     Table[linefunc, {ρ,
       Rescale[Range[0, 1, .5], {0, 1}, {1.5 ρmin, 2}]
       }],
     {φ, 0, 2 π, 20. Degree}] // Flatten[#, 1] &
   ];

dataGrid = Join[outlineGrid, transitionalGrid, fineGrid, minLineGrid];

The grid points would look like this in the $x$-$y$ plane:

Graphics[{PointSize[.001], Black, Point[minLineGrid[[All, 1 ;; 2]]],
  PointSize[.002], Red, Point[fineGrid[[All, 1 ;; 2]]],
  PointSize[.005], Darker@Green, Point[transitionalGrid[[All, 1 ;; 2]]],
  PointSize[.008], Blue, Point[outlineGrid[[All, 1 ;; 2]]]},
 Frame -> True, FrameLabel -> (Style[#, Bold, 20] & /@ {x, y})]

Mathematica graphics

4:

Plot the dataGrid by ListPlot3D:

ListPlot3D[dataGrid,
 PlotRange -> {{-.5, 2}, {-3.5, -1}, {-5, 1}},
 ClippingStyle -> Gray, BoundaryStyle -> Blue, 
 AxesLabel -> (Style[#, Bold, 20] & /@ {x, y})]

Mathematica graphics

($x$,$y$) grid of the plot:

ListPlot3D[dataGrid,
 PlotRange -> {{-.5, 2}, {-3.5, -1}, {-5, 1}},
 Mesh -> All, PlotStyle -> None, 
 AxesLabel -> (Style[#, Bold, 20] & /@ {x, y}), 
 ViewPoint -> {0, 0, ∞}]

Mathematica graphics

Remarks: I believe it can be made more adaptive by finding an appropriate coordinate transformation which converts the ridge-line and the outside borders (-0.5 <= x <= 2 && -3.5 <= y <= -1) to centered concentric circles, and then generate polar grids in this new coordinate system.

Edit: The typical time spent for generating above grids on my computer:

Timing[
 rC = Mean[r0];

 ...

 outlineGrid = ...;
 ]

{1.482, Null}

Improve this answer
$\endgroup$
19
  • $\begingroup$ This looks awesome. I'll be able to test it in a couple of hours and then comment. Thank you. $\endgroup$
    – Pragabhava
    Sep 12, 2012 at 14:43
  • $\begingroup$ @Pragabhava Thanks :) Well I have to say this process highly relies on the special shape of expr. If you need more general method, MaxRecursion would be handy. $\endgroup$
    – Silvia
    Sep 12, 2012 at 14:48
  • 2
    $\begingroup$ @Silvia I get your point but lets try to figure out something in the line of your algorithm so that we have something in our arsenal to even attack general problems for which Plot3D family is too slow. One great article that I found considering dl.acm.org/citation.cfm?doid=383259.383267 $\endgroup$
    – PlatoManiac
    Sep 12, 2012 at 15:25
  • 3
    $\begingroup$ If you got no ACM dgp.toronto.edu/~mooncake/papers/SIGGRAPH2001_Tupper.pdf $\endgroup$
    – PlatoManiac
    Sep 12, 2012 at 15:31
  • 2
    $\begingroup$ @Silvia Sorry for the late acceptance, I wanted to show you how useful your code was. Thanks again! $\endgroup$
    – Pragabhava
    Sep 14, 2012 at 18:08
13
$\begingroup$

Try

exp1 = Sum[v[{x, y}, r0[[i]], sig[[i]], eps[[i]]], {i, 1, 84}];
Plot3D[exp1, {x, -0.5, 2}, {y, -3.5, -1}, PlotPoints -> 60, 
  Mesh -> None, MaxRecursion -> 4, 
  PlotStyle -> 
   Directive[Yellow, Specularity[White, 20], Opacity[0.8]], 
  ExclusionsStyle -> {None, Red}, 
  ColorFunction -> Function[{x, y, z}, Hue[.95 (1. - z)]], 
  AxesLabel -> {x, y}, PlotRange -> {-5, 1}, 
  ClippingStyle -> Directive[Opacity[0.65], Red], 
  BoundaryStyle -> Directive[Red, Thick]] // AbsoluteTiming

enter image description here

4.8882796

This takes about $5$ sec. in my PC. Right combination of PlotPoints and MaxRecursion can help you to have smoother plots. But the more recursion you allow the slower the Plot3D evaluates. Same logic is true for PlotPoints option.

Update

With Compile I could bring this timing to a 2.75 sec.

Clear[exp, v];
v[r_, r0_, s_, ep_] := 4 ep (s^12/Norm[r - r0]^12 - s^6/Norm[r - r0]^6);
exp[s1_, s2_] := Sum[v[{s1, s2}, r0[[i]], sig[[i]], eps[[i]]], {i, 1, 84}];
cf = Compile[{{x, _Real}, {y, _Real}}, 
             Evaluate@(exp[s1, s2] /. s1 -> x /. s2 -> y), 
             CompilationTarget -> "C", RuntimeOptions -> "Speed" ];
Needs["CompiledFunctionTools`"];
CompilePrint[cf] // Short

Returning no error and no call to MainEvaluate!

2 arguments
2 Integer registers
841 Real registers
Underflow checking off
Overflow checking off
Integer overflo...  + R823 + R824 + R825 + R826 + R827 + R828 + R829 + R830 + R732 +
R833 + R834 + R835 + R836 + R837 + R838 + R839
1241    Return

Now call Plot3D

Quiet@Plot3D[cf[x, y], {x, -0.5, 2.}, {y, -3.5, -1.},options]
Improve this answer
$\endgroup$
10
$\begingroup$

You could try to find the optimal combination of PlotPoints and MaxRecursion. High values of both will give you the prettiest plots, but they will take a lot of time. Often a low number of PlotPoints with a higher number of MaxRecursion or vice versa will do. Below you will find your plots with a few combinations of these values and the timings needed (given in the plot title)

Mathematica graphics

Improve this answer
$\endgroup$
2
  • $\begingroup$ So, my best bet would be to separate the plot in two: a plot of the annular region starting from inside the singular zone and extending outside as to contain most of the minima/maxima, and the second region where $V$ is slowly varying (i.e. far from r0). Then use a convenient combination of PlotPoints and MaxRecursion in the first one, and few PlotPoints on the second. Is this a good idea? $\endgroup$
    – Pragabhava
    Sep 11, 2012 at 4:46
  • $\begingroup$ @Pragabhava I believe that could work. Perhaps that the connection of the two will yield visible artifacts. $\endgroup$
    – Sjoerd C. de Vries
    Sep 11, 2012 at 6:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.