77 lines
3.5 KiB
ReStructuredText
77 lines
3.5 KiB
ReStructuredText
|
Computing the Overlap Function
|
||
|
|
==============================
|
||
|
|
|
||
|
|
The overlap function is defined as the portion of particles of a given set,
|
||
|
|
whose positions *overlap* after a given time :math:`t` with the reference configuration at :math:`t=0`.
|
||
|
|
This is calculated as follows:
|
||
|
|
The Initial positions define spheres of a given radius :math:`r` which then are used
|
||
|
|
to test how many of the particles at a later time are found within those spheres.
|
||
|
|
Normalized by the number of spheres this gives the correlation of the configurational overlap.
|
||
|
|
|
||
|
|
.. math::
|
||
|
|
|
||
|
|
Q(t) = \frac{1}{N} \left\langle \sum\limits_{i=1}^N n_i(t) \right\rangle
|
||
|
|
|
||
|
|
Where :math:`n_i(t)` defines the :math:`N` spheres, with :math:`n_i(t)=1` if a particle
|
||
|
|
is found within this sphere at time :math:`t` and :math:`n_i(0) = 1` for :math:`1\leq i \leq N`.
|
||
|
|
|
||
|
|
Evaluation with mdevaluate
|
||
|
|
--------------------------
|
||
|
|
|
||
|
|
Computation of the overlap requires the relatively expensive computation of next neighbor distances,
|
||
|
|
which scales with the order of :math:`\mathcal{O}(N^2)`.
|
||
|
|
There are more efficient ways for the solution of this problem, the one used here is
|
||
|
|
the so called :class:`~scipy.spatial.cKDTree`.
|
||
|
|
This is much more efficient and allows to compute the overlap relatively fast::
|
||
|
|
|
||
|
|
OW = md.open('/path/to/sim').subset(atom_name='OW')
|
||
|
|
tree = md.coordinates.CoordinatesKDTree(OW)
|
||
|
|
Qol = md.correlation.shifted_correlation(
|
||
|
|
partial(md.correlation.overlap, crds_tree=tree, radius=0.11),
|
||
|
|
OW
|
||
|
|
)
|
||
|
|
|
||
|
|
As seen above, mdevaluate provides the function :func:`~mdevaluate.correlation.overlap`
|
||
|
|
for this evaluation, which uses a special object of type :class:`~mdevaluate.coordinates.CoordinatesKDTree`
|
||
|
|
for the neighbor search.
|
||
|
|
The latter provides two features, necessary for the computation:
|
||
|
|
First it computes a :class:`~scipy.spatial.cKDTree` for each necessary frame of the trajectory;
|
||
|
|
second it caches those trees, since assembly of KDTrees is expensive.
|
||
|
|
The size of the cache can be controlled with the keyword argument ``maxsize`` of the CoordinatesKDTree initialization.
|
||
|
|
|
||
|
|
Note that this class uses the C version (hence the lowercase C) rather than
|
||
|
|
the pure Python version :class:`~scipy.spatial.KDTree` since the latter is significantly slower.
|
||
|
|
The only downside is, that the C version had a memory leak before SciPy 0.17,
|
||
|
|
but as long as a recent version of SciPy is used, this shouldn't be a problem.
|
||
|
|
|
||
|
|
Overlap of a Subsystem
|
||
|
|
----------------------
|
||
|
|
|
||
|
|
In many cases the overlap of a subsystem, e.g. a spatial region, should be computed.
|
||
|
|
This is done by selecting a subset of the initial configuration before defining the spheres.
|
||
|
|
The overlap is then probed with the whole system.
|
||
|
|
This has two benefits:
|
||
|
|
|
||
|
|
1. It yields the correct results
|
||
|
|
2. The KDTree structures are smaller and thereby less computation and memory expensive
|
||
|
|
|
||
|
|
An example of a spatial resolved analysis, where ``OW`` is loaded as above::
|
||
|
|
|
||
|
|
selector = partial(
|
||
|
|
md.coordinates.spatial_selector,
|
||
|
|
transform=md.coordinates.spherical_radius,
|
||
|
|
rmin=1.0,
|
||
|
|
rmax=1.5
|
||
|
|
)
|
||
|
|
tree = md.coordinates.CoordinatesKDTree(OW, selector=selector)
|
||
|
|
Qol = md.correlation.shifted_correlation(
|
||
|
|
partial(md.correlation.overlap, crds_tree=tree, radius=0.11),
|
||
|
|
OW
|
||
|
|
)
|
||
|
|
|
||
|
|
This computes the overlap of OW atoms in the region :math:`1.0 \leq r \leq 1.5`.
|
||
|
|
This method can of course be used to probe the overlap of any subsystem, which is selected by the given selector function.
|
||
|
|
It should return a viable index for a (m, 3) sized NumPy array when called with original frame of size (N, 3)::
|
||
|
|
|
||
|
|
subset = frame[selector(frame)]
|