pwtools.crys.rpdf

pwtools.crys.rpdf(trajs, dr=0.05, rmax='auto', amask=None, tmask=None, dmask=None, pbc=True, norm_vmd=False, maxmem=2.0)[source]

Radial pair distribution (pair correlation) function for Structures and Trajectories. In case of trajectories, the time-averaged RPDF is returned. Can also handle non-orthorhombic unit cells (simulation boxes). Only fixed-cell MD at the moment.

Parameters:

  • trajs (Structure or Trajectory or list of one or two such objects) – The case len(trajs)==1 is the same as providing the object directly (most common case). Internally we expand the input to [trajs, trajs], i.e. the RPDF of the 2nd coord set w.r.t. to the first is calculated – the order matters! This is like selection 1 and 2 in VMD, but nornmally you would use amask instead. The option to provide a list of two Trajectory objects exists for cases where you don’t want to use amask, but create two different Trajectory objects outside.

  • dr (float, optional) – Radius spacing. Must have the same unit as cell, e.g. Angstrom.

  • rmax ({'auto', float}, optional) –

    Max. radius up to which minimum image nearest neighbors are counted. For cubic boxes of side length L, this is L/2 [AT,MD].

    ’auto’ : the method of [Smith] is used to calculate the max. sphere
    raduis for any cell shape
    float : set value yourself

  • amask (None, list of one or two bool 1d arrays, list of one or two strings) – Optional atom mask. This is the complementary functionality to sel in vmd_measure_gofr(). If len(amask)==1, then we expand to [amask, amask] internally, which would calculate the RPDF between the same atom selection. If two masks are given, then the first is applied to trajs[0] and the second to trajs[1]. Use this to select only certain atoms in each Trajectory. The default is to provide bool arrays. If you provide strings, they are assumed to be atom names and we create a bool array np.array(symbols) == amask[i].

  • tmask (None or slice object, optional) – Time mask. Slice for the time axis, e.g. to use only every 100th step, starting from step 2000 to the end, use tmask=slice(2000,None,100), which is the same as np.s_[2000::100].

  • dmask (None or string, optional) –

    Distance mask. Restrict to certain distances using numpy syntax for creating bool arrays:

    '>=1.0'
'{d} >=1.0' # the same
'({d} > 1.0) & ({d} < 3.0)'

    where {d} is a placeholder for the distance array (you really have to use {d}). The placeholder is optional in some pattern. This is similar to VMD’s “within” (pbc=False) or “pbwithin” (pbc=True) syntax.

  • pbc (bool, optional) – apply minimum image convention to distances

  • norm_vmd (bool, optional) – Normalize g(r) like in VMD by counting duplicate atoms and normalize to natoms0 * natoms1 - duplicates instead of natoms0*natoms1. Affects all-all correlations only. num_int is not affected. Use this only for testing.

  • maxmem (float, optional) – Maximal allowed memory to use, in GB.

Returns:

  • array (len(rad), 3), the columns are

  • rad (1d array) – radius (x-axis) with spacing dr, each value r[i] is the middle of a histogram bin

  • hist (1d array, (len(rad),)) – the function values g(r)

  • num_int (1d array, (len(rad),)) – the (averaged) number integral number_density*hist*4*pi*r**2.0*dr

Notes

rmax : The maximal rmax for which g(r) is correctly normalized is the result of rmax_smith(), i.e. the radius if the biggest sphere which fits entirely into the cell. This is simply L/2 for cubic boxes of side length L and volume L**3, for instance. We do explicitely allow rmax > rmax_smith for testing, but be aware that g(r) and the number integral are wrong for rmax > rmax_smith.

Even though the number integral will always converge to the number of all neighbors for r -> infinity, the integral value (the number of neigbors) is correct only up to rmax_smith.

See examples/rpdf/ for educational evidence. For notes on how VMD does this, see comments in the code below.

selection : The selection mechanism with amask is in principle as capable as VMD’s, but relies completely on the user’s ability to create bool arrays to filter the atoms. In practice, anything more complicated than array(symbols)=='O' (“name O” in VMD) is much more difficult than VMD’s powerful selection syntax.

Curently, the atom distances are calculated by using numpy fancy indexing. That creates (big) arrays in memory. For data from long MDs, you may run into trouble here. For a 20000 step MD, start by using every 200th step or so (use tmask=slice(None,None,200)) and look at the histogram, as you take more and more points into account (every 100th, 50th step, …). Especially for Car Parrinello, where time steps are small and the structure doesn’t change much, there is no need to use every step. See also maxmem.

Examples

>>> # simple all-all RPDF, time-averaged over all MD steps
>>> d = rpdf(traj)

>>> # the same as rpdf(traj,...)
>>> d = rpdf([traj], ...)
>>> d = rpdf([traj, traj], ...)

>>> # 2 selections: RPDF of all H's around all O's, average time step 3000 to
>>> # end, take every 50th step
>>> traj = io.read_cp2k_md('cp2k.out')
>>> d = rpdf(traj, dr=0.1, amask=['O', 'H'],tmask=np.s_[3000::50])
>>> plot(d[:,0], d[:,1], label='g(r)')
>>> twinx()
>>> plot(d[:,0], d[:,2], label='number integral')

>>> # use bool arrays for `amask`, need this for more complicated pattern
>>> sy = np.array(traj.symbols)
>>> # VMD: sel1='name O', sel2='name H', same as amask=['O', 'H']
>>> d = rpdf(traj, dr=0.1, amask=[sy=='O', sy=='H'],tmask=np.s_[3000::50])
>>> # VMD: sel1='name O', sel2='name H Cl', note that the bool arrays must
>>> # be logically OR'ed (| operator) to get the ffect of "H and Cl"
>>> d = rpdf(traj, dr=0.1, amask=[sy=='O', (sy=='H') | (sy=='Cl')],tmask=np.s_[3000::50])

>>> # skip distances >1 Ang
>>> d = rpdf(traj, dr=0.1, amask=['O', 'H'],tmask=np.s_[3000::50]
...          dmask='{d}>1.0')

References

[AT] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids,

Clarendon Press, 1989

[MD] R. Haberlandt, S. Fritzsche, G. Peinel, K. Heinzinger,

Molekulardynamik - Grundlagen und Anwendungen, Friedrich Vieweg & Sohn Verlagsgesellschaft 1995

[Smith] W. Smith, The Minimum Image Convention in Non-Cubic MD Cells,

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.57.1696 1989