Radial Distribution Functions ¶

Asap has a module for calculating Radial Distribution Functions, or RDFs. An RDF measures the probability of finding an atom at distance r given that there is an atom at position 0; it is thus essentially a histogram of interatomic distances - and is calculated as such.

When creating the RadialDistributionFunction object, at least three parameters must be given

atoms: The atoms object being analyzed.
rMax: The maximal distance up to which the RDF is calculated, in Ångström. Must be less than the minimal dimension of the simulation, otherwise nonsensical results are returned.
nBins: The number of bins used for the histogram, i.e. the resolution of the RDF.

The RadialDistributionFunction object works by first collecting data (by calling the update() method) and then making them available (e.g. by the get_rdf() method).

Calculating the RDF of a single configuration ¶

Calculating the RDF on a single configuration is easy:

from asap3.analysis.rdf import RadialDistributionFunction
[ ... ]
rdf = RadialDistributionFunction(atoms, rMax, nBins).get_rdf()

When get_rdf() is called on a previously unused RadialDistributionFunction object, data is collected and the RDF data are returned immediately. See the section Interpreting the returned RDF for how the returned data is interpreted.

Calculating the RDF of a previously stored trajectory ¶

The RadialDistributionFunction object is intended to collect data from a changing atoms object during a simulation. However, a Trajectory object returns a new atoms object for each configuration. The RadialDistributionFunction object must be initialized with the first of these, and then subsequently fooled into using the following atoms object. This is however not difficult:

traj = PickleTrajectory(filename)
RDFobj = None
for atoms in traj:
    if RDFobj is None:
        RDFobj = RadialDistributionFunction(atoms, rMax, nBins)
    else:
        RDFobj.atoms = atoms  # Fool RDFobj to use the new atoms
    RDFobj.update()           # Collect data
rdf = RDFobj.get_rdf()

TO DO: Add an optional argument to update(), so the new atoms can be specified here.

Calculating RDFs on the fly ¶

The RadialDistributionFunction object (RDFobj) is intended to collect data as a simulation is being made. To do this, the RDFobj is usually attached to the dynamics, and is called after each N timesteps (typically, N=5). The RDFobj can be instructed to save its data to a file after each M updates has occurred, and/or to notify other objects (observers of the RDFobj observer). These observers can then read the RDF with get_rdf(), and process it.

To do this, extra arguments should be passed to RadialDistributionFunction(). In the following, atoms is supposed to be an Atoms object, and dyn is supposed to be a dynamics object acting on atoms (for example a VelocityVerlet object or similar):

RDFobj = RadialDistributionFunction(atoms, rMax, nBins, average=M,
                                    autoclear=True)
dyn.attach(RDFobj.update, interval=N)
RDFobj.output_file("my_rdf")

This will attach the RDFobj to the dynamics, so data is taken after each N timesteps. The last line enables saving the RDF data to a file, and specifies the file names used: my_rdf0000.rdf, my_rdf0001.rdf etc. The data are saved after each M times RDF data are collected (i.e. after each N*M timesteps). If autoclear=True is omitted, the RDF data will be cumulative. If autoclear=True is specified, the data are cleared after saving.

If you wish to process the RDF data on the fly, this is also possible. A function attached to RDFobj with RDFobj.attach(func) will be called at the same time the data is saved (the function is, of course, also called if saving is not enabled).

Interpreting the returned RDF ¶

The data returned by get_rdf() is a one-dimensional array of length nBins. It contains the RDF. The RDF is defined as the average density of atoms at a given distance from another atom. The RDF is normalized by the average density of the simulation, defined as the number of atoms divided by the volume of the simulation. In a homogeneous system it will therefore go to 1.0 at large distances.

In simulations without periodic boundary conditions, the volume of the simulation is ill-defined. This module uses the volume defined by the computational cell returned by atoms.get_cell() for purposes of normalization.

get_rdf() only returns the values of the RDF, not the corresponding distances. The distances are easily calculated from the number of bins and the cutoff distance:

x = (np.arange(nBins) + 0.5) * rMax / nBins

TO DO: Add a trivial function returning x.

RDFs and the coordination number ¶

Average coordination numbers are often extracted from RDFs by integrating the RDF out to the first minimum. Note that in this case, you need to integrate 4 * pi * x^2 * RDF(x) and multiply with the average density rho = len(atoms) / atoms.get_volume(). This is illustrated in the second example below.

Saving and loading RDFs ¶

The RDF data can be saved automatically as shown in the section Calculating RDFs on the fly above. This mechanism can also be tricked into saving manually. After collecting the data save it like this:

RDFobj.output_file("my_rdf")
RDFobj.save()

TO DO: Add an optional file name to .save(), so RDFs can be saved manually without interfering with the autosave mechanism.

The saved RDFs can later be loaded using the class method RadialDistributionFunction.load(), the data can then be fetched with .get_rdf(). Note that load is a class method, you do not need to create a RadialDistributionFunction instance first, but can call it directly with the class name. It returns an instance.

RDFobj = RadialDistributionFunction.load("my_rdf0000.rdf")
rdf = RDFobj.get_rdf()

Partial RDFs: Looking at specific elements or subsets of atoms ¶

It is often useful to look at partial RDFs, for example RDFs only taking some elements into account, for example to get the distribution of atoms of element B around an atom of element A. Do do this, call get_rdf() with the optional argument, elements. It must be a tuple of two atomic numbers (a, b), the returned RDF then tells how many b neighbors an a atom has.

It is also possible to group atoms according to other criteria, for example to calculate RDFs in different parts of space. In this case, the atoms must be divided into groups when the RadialDistributionFunction object is created. Pass the extra argument groups when creating the object, it must be an array containing a non-negative integer for each atom specifying its group. When calling get_rdf() use the argument group to specify for which group of atoms you want the RDF.

IMPORTANT: The partial RDFs are normalized such that they sum up to the global RDF. This means that integrating the first peak of a partial RDF obtained with elements=(a,b) does not give you the number of B atoms in the first shell around the A atom. Instead it gives this coordination number multiplied with the concentration of A atoms.

Documentation of the RadialDistributionFunction class ¶

Creating RadialDistributionFunction objects ¶

The constructor takes the following arguments:

atoms: The atoms being analyzed.
rMax: (A float). The maximal distance out to which the RDF is calculated. Should be smaller than the system size. Both memory consumption and calculation time scale with rMax cubed!
nBins: (A positive integer). The resolution of the RDF, i.e. how many points are calculated on the curve. Increasing this has only a small effect on memory consumption and calculation time (unless it is increased to absurd values).
groups: (Optional. An array of non-negative integers, one for each atom). Defines which groups the atoms belong to.
average: (Optional integer, defaults to 1). Specifies how often observers to the RadialDistributionFunction object are called, and how often the data are saved (if enabled). This is thus the number of snapshots over which the RDF is averaged.
autoclear: (Optional boolean, defaults to False). Specifies if the data is cleared after observers are called and data is saved. (The actually clearing is delayed until the next time update() is called).
verbose: (Optional boolean, defaults to False). If true, prints a line to standard output when a calculation starts and stops.

The methods ¶

update(): Collect data from the atoms. Attach this function to the dynamics, or call it manually when you want to collect data.
clear(): Clear the accumulated data. Normally called automatically, but can also be called explicitly.
get_rdf(group=None, elements=None): Get the RDF. If group and/or elements are specified, a partial RDF is returned. If present, group must be an integer found in the groups array specified when the object was created. If present, elements must be a tuple of two atomic numbers. See the discussion of partial RDFs above.
output_file(prefix): Turn on the autosaving feature, and specify the prefix of the file name. The prefix should be a string, and should not contain the “%” character unless you know what you are doing.
save(): Trigger a save into a file. Not intended to be called manually, but you can do it as shown in the section Saving and loading RDFs.
load(): A class method creating a new object from a saved file.
attach(callable, ...): Attach an observer. The observer is called each average time update() is called. THIS IS APPARENTLY BROKEN.

Example: Plotting the RDF ¶

RDF.py

from asap3 import *
from asap3.md.velocitydistribution import MaxwellBoltzmannDistribution
from asap3.analysis.rdf import RadialDistributionFunction
from ase.lattice.cubic import *
import matplotlib.pyplot as plt

# We want RDF out to a distance of 15 Angstrom, with 200 bins
rng=15.0
bins = 200

# Create a block of Cu atoms, expand the crystal 10% to leave space
# for thermal expansion and to facilitate melting.  Describe
# interactions with EMT, and set a very high initial temperature so it
# melts.
atoms = FaceCenteredCubic(size=(20,20,20), symbol="Cu", pbc=True)
atoms.set_cell(atoms.get_cell() * 1.1, scale_atoms=True)

atoms.set_calculator(EMT())
MaxwellBoltzmannDistribution(atoms, 3500*units.kB)

# First, run 25 time steps before taking data.
dyn = VelocityVerlet(atoms, 4*units.fs, logfile='-')
dyn.run(25)

# Make the RDF object, attach it to the dynamics
RDFobj = RadialDistributionFunction(atoms, rng, bins, verbose=True)
dyn.attach(RDFobj.update, interval=5)

# Run the simulation
dyn.run(100)

# Get the RDF and plot it.
rdf = RDFobj.get_rdf()
x = np.arange(bins) * rng / bins
plt.plot(x, rdf)
plt.show()

The resulting RDF will look like this:

Example: Element-specific RDFs and coordination numbers ¶

RDF2.py

from asap3 import *
from asap3.md.velocitydistribution import MaxwellBoltzmannDistribution
from asap3.analysis.rdf import RadialDistributionFunction
from ase.lattice.compounds import *
from ase import data
import matplotlib.pyplot as plt
import numpy as np

# We want RDF out to a distance of 15 Angstrom, with 200 bins
rng=10.0
bins = 200

# Helper function calculating coordination number from RDF
def getcoord(rdf, r, dens):
    "Calculate coordination number from RDF."
    # Find first peak by assuming it is the largest
    imax = rdf.argmax()
    # Find first minimum, assuming it is the deepest after the maximum.
    imin = imax + rdf[imax:].argmin()
    print("Position of first peak and following minimum:", imax, imin)
    dr = r[1] - r[0]
    integral = 0.0
    for i in range(imin):
        integral += 4 * np.pi * r[i] * r[i] * dr * rdf[i]
    return integral * dens

# Create a block of Cu atoms, expand the crystal 10% to leave space
# for thermal expansion and to facilitate melting.  Describe
# interactions with EMT, and set a very high initial temperature so it
# melts.
atoms = L1_2(size=(20,20,20), symbol=("Au", "Cu"), pbc=True,
             latticeconstant=3.75)

atoms.set_calculator(EMT())
MaxwellBoltzmannDistribution(atoms, 600*units.kB)

# First, run 25 time steps before taking data.
dyn = VelocityVerlet(atoms, 5*units.fs, logfile='-')
dyn.run(25)

# Make the RDF object, attach it to the dynamics
RDFobj = RadialDistributionFunction(atoms, rng, bins, verbose=True)
dyn.attach(RDFobj.update, interval=5)

# Run the simulation
dyn.run(100)

# Get the RDF and plot it.
rdf = RDFobj.get_rdf()
x = np.arange(bins) * rng / bins
plt.plot(x, rdf, 'k')  # Black

# Get the partial RDFs and plot them
Au = data.atomic_numbers['Au']
Cu = data.atomic_numbers['Cu']
rdfAuAu = RDFobj.get_rdf(elements=(Au, Au))
plt.plot(x, rdfAuAu, 'r')
rdfCuAu = RDFobj.get_rdf(elements=(Cu, Au))
plt.plot(x, rdfCuAu, 'b--')
rdfAuCu = RDFobj.get_rdf(elements=(Au, Cu))
plt.plot(x, rdfAuCu, 'c')
rdfCuCu = RDFobj.get_rdf(elements=(Cu, Cu))
plt.plot(x, rdfCuCu, 'g')

rho = len(atoms)/atoms.get_volume()
c = getcoord(rdf, x, rho)
cAuAu = getcoord(rdfAuAu, x, rho)
cAuAu = 0.0  # RDF found next-nearest neighbors as no nearest.
cCuAu = getcoord(rdfCuAu, x, rho) / 0.75
cAuCu = getcoord(rdfAuCu, x, rho) / 0.25
cCuCu = getcoord(rdfCuCu, x, rho) / 0.75

print("Total coodination number:", c)
print("Cu atoms around an Cu atom:", cCuCu)
print("Au atoms around an Au atom:", cAuAu)
print("Cu atoms around an Au atom:", cAuCu)
print("Au atoms around an Cu atom:", cCuAu)

plt.show()

Note that the coordination numbers obtained from the RDF are not 12, 8 and 4 as one would expect, but slightly below. Presumably this is due to the numerical integration of the RDF, but it should be investigated.

Also note that the algorithm for finding the nearest-neighbor peak is quite primitive, it assumes that the first peak is the highest. It fails in the case of the Au-Au RDF, since the first peak is zero here (in the ordered AuCu3 lattice, gold atoms are never nearest neighbors).