Function Repository Resource:

RadialDistributionFunctionList

Source Notebook

Compute the relative probability of finding a point at a given distance from another point

Contributed by: Nicholas E. Brunk, Wolfram|Alpha Math Team

ResourceFunction["RadialDistributionFunctionList"][{p1,p2,},L]

returns the radial distribution function resulting from each pair of points pi1,pi2 within a periodic box of uniform dimensions L.

ResourceFunction["RadialDistributionFunctionList"][{p1,p2,},L,w]

returns the radial distribution function using separation bin width w.

Details

The radial distribution function is also known as the pair correlation function.
The default value for the bin width (w = 0.005), is a common choice in reduced units quantified relative to the particle diameter.
ResourceFunction["RadialDistributionFunctionList"] will work in up to three dimensions, applying periodicity uniformly in each.
ResourceFunction["RadialDistributionFunctionList"] normalizes the y axis to generate a density relative to that of an ideal gas at the same density in the same space (determined by box size L and the number and dimension of the specified coordinates).
ResourceFunction["RadialDistributionFunctionList"] takes the option "ReferenceCoordinates" {pj1,pj2,} to compute the RDF gij=gji between two pairs of coordinates {pi1,pi2,} and {pj1,pj2,}.
ResourceFunction["RadialDistributionFunctionList"] utilizes Compile for performance enhancement, and takes the same options.

Examples

Basic Examples (2) 

Compute the radial distribution function from two points:

In[1]:=
ResourceFunction[
  "RadialDistributionFunctionList"][{{0, 0, 0}, {0.5, 0, 0}}, 1.] // Short
Out[1]=

Compute the radial distribution function of a set of coordinates enclosed in a periodic 1D "box" of size (L = 4):

In[2]:=
coords = {{-1}, {0}, {1}};
rdfData = ResourceFunction["RadialDistributionFunctionList"][coords, 4];

Display the coordinates in number line format along the axis for which the positions differ:

In[3]:=
NumberLinePlot[First /@ coords]
Out[3]=

The radial distribution function has non-zero peaks at positions near those of the exact distances apart:

In[4]:=
Cases[rdfData, elem_ /; elem[[2]] != 0]
Out[4]=

Plot the radial distribution function characteristic of this system:

In[5]:=
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[5]=

Options (2) 

Cross-correlations can be computed using the option "ReferenceCoordinates":

In[6]:=
(* Evaluate this cell to get the example input *) CloudGet["https://www.wolframcloud.com/obj/9ab5f2b3-4920-4b6a-84bf-dc8612b13d40"]
Out[7]=

Radial distribution function takes the same options as Compile, and benefits from automated compilation to C:

In[8]:=
ResourceFunction["RadialDistributionFunctionList"][
   Table[RandomReal[{0, 1.}, 3], {10000}], 1., 0.005] // AbsoluteTiming // First
Out[8]=
In[9]:=
ResourceFunction["RadialDistributionFunctionList"][
   Table[RandomReal[{0, 1.}, 3], {10000}], 1., 0.005, CompilationTarget -> "C"] // AbsoluteTiming // First
Out[9]=

Applications (2) 

Compute the radial distribution function characteristic of a cube with unit length, enclosed in a sufficiently larger periodic box:

In[10]:=
boxLength = 5.0;
coords = N@PolyhedronData["Cube", "Vertices"];
rdfData = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength, 0.005];
Graphics3D[{{Opacity[0.25], Cube[boxLength]}, {Sphere[#, 0.5] & /@ coords}}, Sequence[
 Boxed -> False, Axes -> True]]
Out[11]=

Plot the resulting radial distribution function, where you can see nearest neighbor, polygonal face hypotenuse, and internal hypotenuse distances:

In[12]:=
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[12]=

Properties and Relations (4) 

By definition, the RDF of an ideal gas converges to 1 at all distances when provided sufficient data (here, 1000 time steps on 100 positions of 100 particles):

In[13]:=
bulkCoordsData = Table[RandomReal[{0, 1.}, 3], {1000}, {100}];

Compute both the radial distribution function of just the first time step, and then that of the ensemble average over all time steps:

In[14]:=
singleTimestepRDF = ResourceFunction["RadialDistributionFunctionList"][
   First@bulkCoordsData, 1.0, 0.005];
ensAvgRDF = Mean[Map[
    ResourceFunction["RadialDistributionFunctionList"][#, 1.0, 0.005] &, bulkCoordsData]];

As expected, the radial distribution function converges to nearly 1 as sufficient data is introduced:

In[15]:=
ListLinePlot[{singleTimestepRDF, ensAvgRDF}, Sequence[
 PlotRange -> All, PlotLegends -> {"One Timestep", "Ensemble Average"},
   Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[15]=

Normalization by the expected number of particles in that vicinity works in 2D cases as well:

In[16]:=
bulkCoordsData = Table[RandomReal[{0, 1.}, 2], {1000}, {100}];
singleTimestepRDF = ResourceFunction["RadialDistributionFunctionList"][
   First@bulkCoordsData, 1.0, 0.005];
ensAvgRDF = Mean[Map[
    ResourceFunction["RadialDistributionFunctionList"][#, 1.0, 0.005] &, bulkCoordsData]];
In[17]:=
ListLinePlot[{singleTimestepRDF, ensAvgRDF}, Sequence[
 PlotRange -> All, PlotLegends -> {"One Timestep", "Ensemble Average"},
   Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[17]=

The same is true for 1D cases:

In[18]:=
bulkCoordsData = Table[{RandomReal[{0, 1.0}]}, {1000}, {100}];
singleTimestepRDF = ResourceFunction["RadialDistributionFunctionList"][
   First@bulkCoordsData, 1.0, 0.005];
ensAvgRDF = Mean[Map[
    ResourceFunction["RadialDistributionFunctionList"][#, 1.0, 0.005] &, bulkCoordsData]];
In[19]:=
ListLinePlot[{singleTimestepRDF, ensAvgRDF}, Sequence[
 PlotRange -> All, PlotLegends -> {"One Timestep", "Ensemble Average"},
   Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[19]=

Ensemble average cross-correlations between two distinct populations (lists of coordinates) can also be computed:

In[20]:=
bulkCoordsData1 = Table[RandomReal[{0, 1.0}, 3], {1000}, {100}];
bulkCoordsData2 = Table[RandomReal[{0, 1.0}, 3], {1000}, {100}];
rdfData12 = Mean@Table[
    ResourceFunction["RadialDistributionFunctionList"][
     bulkCoordsData1[[i]], 1.0, .005, "ReferenceCoordinates" -> bulkCoordsData2[[i]]], {i, Length@bulkCoordsData1}];

ListLinePlot[rdfData12, PlotRange -> All]
Out[21]=

Cross-correlations are symmetric, that is, the computation is the same for swapped populations g{XMLElement[i, {}, {XMLElement[span, {class -> stylebox}, {ij}]}]}(r)=g{XMLElement[i, {}, {XMLElement[span, {class -> stylebox}, {ji}]}]}(r):

In[22]:=
rdfData21 = Mean@Table[
    ResourceFunction["RadialDistributionFunctionList"][
     bulkCoordsData2[[i]], 1.0, .005, "ReferenceCoordinates" -> bulkCoordsData1[[i]]], {i, Length@bulkCoordsData1}];

ListLinePlot[{rdfData12, rdfData21}, PlotRange -> All, PlotStyle -> {Automatic, {Dashed, Thick}}]
Out[8]=
In[23]:=
rdfData12 == rdfData21
Out[23]=

Possible Issues (4) 

For finite-sized systems, the radial distribution function is accurate out to and is thus truncated at this point, beyond which it would artificially decay:

In[24]:=
boxLength = 5;
coords = {{0, 0, 0}, {1, 0, 0}};
rdfData = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength];
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[3]=

If the box is sufficiently larger than the domain of the particle coordinates, the calculation is as expected (using real particle positions, rather than the minimum image convention under periodicity):

In[25]:=
boxLength = 5;
coords = N@PolyhedronData["Cube", "Vertices"];
rdfData = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength, 0.005];
Graphics3D[{{Opacity[0.25], Cube[boxLength]}, {Sphere[#, 0.5] & /@ coords}}, Sequence[
 Boxed -> False, Axes -> True]]
Out[26]=

The peaks are at the expected unit distance for the touching particles, at the polygonal face hypotenuse distance, and the interior hypotenuse:

In[27]:=
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, PlotLegends -> {"One Timestep", "Ensemble Average"},
   Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[27]=

If the box size (L) is too small, the minimum image convention is employed, applying periodicity in each dimension:

In[28]:=
boxLength = 1;
coords = N@PolyhedronData["Cube", "Vertices"];
rdfData = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength, 0.005];
Graphics3D[{{Opacity[0.25], Cube[boxLength]}, {Sphere[#, 0.25] & /@ coords}}, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[31]=

In this case, particles are perfectly overlapping and the separation distance between all pairs is r{XMLElement[i, {}, {XMLElement[span, {class -> stylebox}, {ij}]}]}=0, however the convention is to put them in the first non-zero bin:

In[32]:=
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[32]=

Modifying the concentric shell bin width can affect normalization and positional accuracy, similar to the behavior in Riemann sums with too few discretizing rectangles:

In[33]:=
boxLength = 10;
coords = {{-1}, {0}, {1}};
NumberLinePlot@Flatten@coords
Out[35]=

Compare the difference when using two different bin widths:

In[36]:=
rdfDataLowBinWidth = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength, 0.005];
rdfDataHighBinWidth = ResourceFunction["RadialDistributionFunctionList"][coords, boxLength, 0.5];

ListLinePlot[{rdfDataLowBinWidth, rdfDataHighBinWidth}, PlotRange -> {All, {0, 20}}, PlotLegends -> {"Low Bin Width (w = 0.005)", "High Bin Width (w = 0.5)"}, Frame -> True, FrameLabel -> (Style[#1, 16] &) /@ {"Distance", "Relative Density"}]
Out[38]=

Neat Examples (3) 

Compute the radial distribution function characteristic of an icosahedral distribution of particles and compare it with that of a cube, each with a minimum separation distance of 1:

In[39]:=
boxLength = 5;
icosaCoords = N@PolyhedronData["Icosahedron", "Vertices"];
cubeCoords = N@PolyhedronData["Cube", "Vertices"];
rdfData = ResourceFunction["RadialDistributionFunctionList"][#, boxLength, 0.005] & /@ {icosaCoords, cubeCoords};
Row@(Graphics3D[{{Opacity[0.25], Cube[boxLength]}, {Sphere[#, 0.5] & /@ #}}, Boxed -> False, Axes -> True, ViewPoint -> {1.35, -3, 0.8}, ImageSize -> 350] & /@ {icosaCoords, cubeCoords})
Out[43]=

Plot the resulting radial distribution functions to compare:

In[44]:=
ListLinePlot[rdfData, Sequence[
 PlotRange -> All, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[44]=

Compare the radial distribution function of both a simple cubic and face-centered cubic (FCC) crystalline lattice (neglecting periodicity, each with nearest neighbors touching at unit lengths):

In[45]:=
cubicCoords = N@Flatten[
    Table[{i, j, k}, {i, -2, 2, 1}, {j, -2, 2, 1}, {k, -2, 2, 1}], 2];

\[ScriptCapitalB] = Normal@LatticeData["FaceCenteredCubic", "Basis"];
fccCoords = N@Flatten[
    Table[i \[ScriptCapitalB][[1]] + j \[ScriptCapitalB][[2]] + k \[ScriptCapitalB][[3]], {i, 1, 5}, {j, 1, 5}, {k, 1, 5}], 2];
fccCoords = fccCoords/
   Min@DeleteCases[Flatten[DistanceMatrix@fccCoords], x_Real /; x == 0.];
Row@(Graphics3D[{{Sphere[#, 0.5] & /@ #}}, Boxed -> False, Axes -> True, ImageSize -> 300] & /@ {cubicCoords, fccCoords})
Out[49]=

Compute and plot the radial distribution functions to compare:

In[50]:=
rdfData = ResourceFunction["RadialDistributionFunctionList"][#, 20, 0.005] & /@ {cubicCoords, fccCoords};
ListLinePlot[rdfData, Sequence[
 PlotRange -> {{0, 5}, Automatic}, Frame -> True, FrameLabel -> Map[Style[#, 16]& , {"Distance", "Relative Density"}]]
 ]
Out[51]=

Compute the Potential of Mean Force that has given rise to the interactions, in this example showing that the 3D ideal gas has no potential of interaction:

In[52]:=
bulkCoordsData = Table[RandomReal[{0, 1.}, 3], {1000}, {100}];
ensAvgRDF = Mean[Map[
    ResourceFunction["RadialDistributionFunctionList"][#, 1.0, 0.005] &, bulkCoordsData]];
In[53]:=
ListLinePlot[ensAvgRDF, PlotRange -> All, Frame -> True, FrameLabel -> (Style[#1, 16] &) /@ {"Distance", "Relative Density"}]
Out[53]=

Compute the corresponding Potential of Mean Force (in units where kBT=1):

In[54]:=
pmf = ensAvgRDF /. {r_, gr_} :> {r, -Log[gr]};
In[55]:=
ListLinePlot[ensAvgRDF /. {r_, gr_} :> {r, -Log[gr]}, PlotRange -> All, Frame -> True, FrameLabel -> (Style[#1, 16] &) /@ {"Distance", "Potential of Mean Force"}]
Out[55]=

It can be seen that despite the noise associated finite data, the Potential of Mean Force is fairly close to zero when averaged over the entire domain:

In[56]:=
Mean[Last /@ pmf]
Out[56]=

Publisher

Wolfram|Alpha Math Team

Version History

  • 2.0.0 – 23 March 2023
  • 1.1.0 – 27 September 2022
  • 1.0.0 – 16 September 2020

Source Metadata

Related Resources

Author Notes

This function's source code could be generalized to periodic boxes of unequal dimensions (Lx LyLz).

FileNameJoin[ReplacePart[FileNameSplit[FindFile["ResourceFunctionHelpers`"]], -1 "RadialDistributionFunctionList.wl"]] // SystemOpen

License Information