Wolfram Function Repository

Instant-use add-on functions for the Wolfram Language

Function Repository Resource:

 RadialDistributionFunction (1.0.0) current version: 1.0.1 »

Source Notebook

Compute the radial distribution function of a set of n-dimensional points

Contributed by: Robert Nachbar

ResourceFunction["RadialDistributionFunction"][points,region,width]

gives an InterpolatingFunction of the radial distribution function of the input points in n-dimensional space in the specified region with a binning resolution of width.

Details

The radial distribuiton function g(r) describes how the density of a system of particles varies as a function of distance from a reference particle.
The region is a Cuboid or a Ball.
The dimension of points and region must be the same.
All elements of points must lie within region.
To avoid edge effects, the distances are computed only for points in an inner region with physical dimensions 1/3 of region.

The histogramming of distances is done over two sets of overlapping bins, effectively giving a moving average.

Examples

Basic Examples (2) 

Define random points within a ball and visualize them:

In[1]:=
pts = RandomPoint[Ball[{0, 0}, 5], 5000]; Graphics[{PointSize[Tiny], Point[pts]}, ImageSize -> Small]
Out[2]=

Compute the RDF for the points:

In[3]:=
\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][pts, Ball[{0, 0}, 5], 0.01]
Out[3]=

Plot it:

In[4]:=
Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 5/3}, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r", "g(r)"}]
Out[4]=

Detect the points in a crystal diffraction image:

In[5]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/a3fb542d-58e1-465f-9052-2f39423875d3"]
Out[7]=
In[8]:=
region = MinMax /@ Transpose@imagePoints // Transpose // MapThread[Map, {{Floor, Ceiling}, #}] & // Apply[Cuboid]
Out[8]=

Compute the radial distribution function:

In[9]:=
\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][imagePoints, region, 1]
Out[9]=

Plot the radial distribution function:

In[10]:=
Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 150}, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r", "g(r)"}]
Out[10]=

Scope (2) 

Detect the points in a crystal diffraction image:

Generate pseudo 1-dimensional data from the image by extracting a narrow slice:

In[11]:=
pointsLine = Cases[imagePoints, {x_, y_} /; 504 <= y < 518 :> {x}];
In[12]:=
ListPlot[Append[#, 0] & /@ pointsLine, AspectRatio -> 1/20, ImageSize -> Large, PlotStyle -> AbsolutePointSize[2]]
Out[12]=

Define a 1-dimensional cuboid region:

In[13]:=
region = Cuboid[{0}, {1024}];

Compute the radial distribution function:

In[14]:=
\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][pointsLine, region, 2]
Out[14]=
In[15]:=
Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 150}, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r", "g(r)"}]
Out[15]=

Use a circular region to generate the radial distribution function for sunflower seeds:

In[16]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/e92a5657-3d09-4f22-b737-5d590e5e5fe4"]
Out[16]=

Use a mask to remove the petals:

In[17]:=
mask = Table[ If[EuclideanDistance[{128, 103}, {x, 196 - y + 1}] < 70, 1, 0], {y, 196}, {x, 257}] // Image; maskedImage = ImageApply[1 &, image, Masking -> ColorNegate@mask]
Out[13]=

Detect the points representing the seeds:

In[18]:=
seedPoints = Values@ComponentMeasurements[ Binarize@ ImageSubtract[maskedImage, BilateralFilter[image, 4, 1]], "Centroid"] // Select[Last[#] > 36 &] // Echo[#, "", Length] &; ListPlot[seedPoints, AspectRatio -> Automatic, PlotRange -> MinMax /@ Transpose@seedPoints, PlotRangePadding -> Scaled[.05], ImageSize -> Small]
Out[15]=

Define the region (essentially the mask):

In[19]:=
region = Ball[{128, 103}, 70];

Compute the radial distribution function:

In[20]:=
\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][seedPoints, region, 0.5]
Out[20]=
In[21]:=
Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 20}, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r", "g(r)"}]
Out[21]=

The integral of the first peak gives the approximate number of nearest nearest neighbors for a seed:

In[22]:=
NIntegrate[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 2, 6}]
Out[22]=

Applications (8) 

Here is a ball of 1118 water molecules approximately 40 Å in diameter:

In[23]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/dae5788c-0980-40ea-8ef6-15c51fdcc53a"]
Out[23]=

Extract the coordinates of the oxygen atoms:

In[24]:=
oxygenAtoms = Thread[{waterBall["AtomList"], waterBall["AtomCoordinates"]}] // Cases[#, {Atom["O"], xyz_} :> QuantityMagnitude@xyz] &; oxygenAtoms // Short
Out[25]=

Compute the region by first finding the convex hull of the oxygen atoms:

In[26]:=
(hull = ConvexHullRegion@oxygenAtoms) // Graphics3D[#, ImageSize -> Small] &
Out[26]=

Then find its centroid, find the distance of the hull point furthest from the centroid, and pad for the approximate diameter of a water molecule of 1.93 Å:

In[27]:=
region = With[{points = First@hull}, Ball[Mean@points, Max[DistanceMatrix[{Mean@points}, points]] + 1.93]]
Out[27]=

Compute the radial distribution function:

In[28]:=
\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][oxygenAtoms, region, 1/10]
Out[28]=

Plot it:

In[29]:=
Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 8}, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r (\[Angstrom])", "g(r)"}]
Out[29]=

The approximate number of nearest oxygen atoms (water molecules) is:

In[30]:=
NIntegrate[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 2.5, 3.1}]
Out[30]=

The average distance is:

In[31]:=
First@Quiet@ FindArgMax[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][ r], {r, 2.8}]
Out[31]=

Neat Examples (2) 

Interactive explorer:

In[32]:=
DynamicModule[{nmax = 3}, Manipulate[ Module[{M = Floor[10^nmax], pts, \[ScriptR]\[ScriptD]\[ScriptF]}, pts = Table[ Sqrt[k/M] {Cos[\[Omega] k], Sin[\[Omega] k]}, {k, 0., M}]; \[ScriptR]\[ScriptD]\[ScriptF] = ResourceFunction["RadialDistributionFunction"][pts, Ball[{0, 0}, 1], binWidth]; GraphicsRow[{Graphics[{GrayLevel[0.95], Disk[], GrayLevel[0.9], Disk[{0, 0}, 1/2], Darker[Blue], AbsolutePointSize[1], Point[pts]}, ImageSize -> Small], Plot[\[ScriptR]\[ScriptD]\[ScriptF][r], {r, 0, 1/2}, PlotRange -> {0, 5}, PlotRangePadding -> Scaled[.05], Axes -> False, Frame -> True, ImageSize -> Small(*,AspectRatio-> 4/5*)]}]], Row[{Control[{{nmax, 3}, 1, 4.5}], " ", Dynamic[Floor[10^nmax]]}], {{\[Omega], 1}, 0, 10, Appearance -> "Labeled"}, Delimiter, {{binWidth, 0.01}, 0.005, 0.1, Appearance -> "Labeled"}] ]
Out[32]=

Generate a large patch of the Hat tiling:

In[33]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/369cf84d-19e7-46e1-9638-db13e7a3e1bb"]
Out[17]=

Compute the mean of the coordinates of each polygon, and then group by color:

In[34]:=
points = GroupBy[First@ tiling, (Replace[#[[2]], {LightBlue -> Blue, Green -> Darker@Green}] &) -> (Apply[Mean, #[[3]]] &)] // N;
In[35]:=
KeyValueMap[{color, coords} |-> Graphics[{color, AbsolutePointSize[2], Point[coords]}, ImageSize -> Small], points] // Row
Out[35]=

Select the points within the inscribed circle of the hexagonal patch:

In[36]:=
points = Select[#, With[{center = Mean[#]}, RegionMember[Ball[center, 94]]]] & /@ points;

Compute the radial distribution function for each group of points and plot it along with the points:

In[37]:=
KeyValueMap[{color, coords} |-> Column[{Graphics[{{color, AbsolutePointSize[2], Point[coords]}}, ImageSize -> Small], \[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF] = ResourceFunction["RadialDistributionFunction"][coords, Ball[Mean[coords], 95], 1]; Plot[\[ScriptCapitalR]\[ScriptCapitalD]\[ScriptCapitalF][r], {r, 0, 60}, PlotStyle -> color, ImageSize -> Small, PlotRange -> All, Axes -> False, Frame -> True, FrameLabel -> {"r", "g(r)"}]}], points] // Row[#, Spacer[9]] &
Out[37]=

Publisher

Robert Nachbar

Version History

  • 1.0.1 – 06 September 2024
  • 1.0.0 – 30 August 2024

Related Resources

License Information