Correlator.hpp File Reference

This module computes correlations (and other two time averages) on the fly and from files. More...

#include "AccumulatorBase.hpp"
#include "observables/Observable.hpp"
#include "system/System.hpp"
#include <utils/Vector.hpp>
#include <boost/multi_array.hpp>
#include <boost/serialization/access.hpp>
#include <cstddef>
#include <memory>
#include <string>
#include <utility>
#include <vector>

Include dependency graph for core/accumulators/Correlator.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
class	Accumulators::Correlator
	The main correlator class. More...

Namespaces
namespace	Accumulators

Detailed Description

This module computes correlations (and other two time averages) on the fly and from files.

The basic idea is that the user can write arbitrary functions A and B that can depend on e.g. particle coordinates or whatever state of the MD box. These functions can be vector-valued, always indicated by dim_A and dim_B.

The way they can be correlated is can be formulated as a (vector-valued) operation on A and B. One example would be to calculate the product component by component, and if A and B are both the particle velocities, then one would obtain

\[ <v_1^x v_1^x(t)>, <v_1^y v_1^y(t)>, <v_1^z v_1^z(t)>, <v_2^x v_2^x(t)>, <v_2^y v_2^y(t)>, <v_2^z v_2^z(t)>, \ldots \]

The idea of the algorithm is to not keep all As and Bs in memory, but feed the algorithm (or the container object) successively by new As and Bs and new correlation estimates are calculated on the fly. Those As and Bs which have exceeded a certain age are first compressed (see below) and then discarded.

To save memory, increase statistics and make the calculation possible for many orders of magnitude in time, the order-n blocking algorithm in [19] (algorithm 9, section 4.4.2 p. 95) is used. Thus not all As and Bs of the whole "past" are stored but some of them are blocked. In this implementation, a blocking based on 2 is always applied: all As and Bs not older than a certain tau_lin are stored as they were, those which are older are stored as average values. All As and Bs older than 2*tau_lin are compressed in blocks of four, etc.

This leads to a hierarchical "history": on level 0 the last tau_lin values are stored. This is done in a cyclic array: the newest is always appended at the end of the array, and if the array length is reached, values are appended at the beginning, overwriting older values. We therefore have to carry the index newest[0] which tells, what is the index of the newest value of A and B.

As soon as the first row of As and Bs is full, two of them are compressed by calculating the arithmetic mean, and stored in the first element on level 1. Then we can overwrite these two values on level 0, because we have saved them. If necessary, we produce space on level 0 by compressing to level 1. On level 1 an array with tau_lin entries is available to store the level-1-compressed values. It is successively filled and also cyclic. When it is filled, the values are stored on level 2 and so on.

This allows to have a "history" over many orders of magnitude in time, without the full memory effort.

Correlations are only calculated on each level. For tau=1,2,..,tau_lin the values are taken from level 1. For tau=tau_lin, tau_lin+2, .., 2*tau_lin we take the values from level 2. On level 2 we also have values for 0,..tau_lin-2, but these are discarded as we have "better" estimates on level 1.

The functions A and B can get extra arguments. This can e.g. be an integer describing the "type" of the particles to be considered, or in case of files, it is a file_data_source pointer, which tells the function from where to read. These arguments have to be a single pointer to an object that carries all relevant information to obtain A and B (from the MD Box or somewhere else).

The correlation has to be initialized with all necessary information, i.e. all function pointers, the dimensions of A and B and their dimensions, etc.

TODO: There is a lot of stuff to do:

• Expand the file_data_source so that one can specify which columns of the file are to be processed
• calculate an estimate of average values. This might be even necessary to calculate \(\left<\left(A(\tau)-\left<A\right>\right)\left(B(\tau)-\left<B\right>\right)\right>\), which is often probably what people want
• Use the A_args to calculate As and Bs only for particular particle types (especially and example, so that other people can follow)
• Use the A_args to calculate molecular stuff in combination with the topology concept
• Write a destructor

Definition in file core/accumulators/Correlator.hpp.