1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
.. _workflow:
.. index:: workflow
Workflow
********
.. figure:: _static/workflow.svg
:align: center
:scale: 6%
Schematic illustration of the workflow in :program:`dynasor`.
Text in red marks user inputs.
The time window, indicated in blue, is moved through the trajectory until the specified maximum number of frames is reached.
Here, some additional information regarding the input parameters is provided, specifically for calculating dynamic and static structure factors.
Additional details concerning these input parameters can be found in the :ref:`documentation of the python interface <interface_python>`.
Internal units
--------------
Internally :program:`dynasor` uses units of Ångström (Å) for length, femtosecond (fs) for time, and electronvolt (eV) for energy.
This means that all calculated correlation functions will be in terms of these units.
Note that :math:`\boldsymbol{q}`-points are defined to include the factor :math:`2\pi` as is commonly done in physics `wavenumber wikipedia <https://en.wikipedia.org/wiki/Wavenumber#Definition>`_. This means that :math:`\boldsymbol{q}`-points technically have units of rad/Å, but note that this is usually simply written as 1/Å.
The trajectory
--------------
The MD trajectory must have a constant (not varying with time) simulation cell.
Furthermore, when analyzing dynamic properties of systems the microcanonical (NVE) ensemble is usually preferred over the canonical ensemble (NVT), as thermostats influence the dynamics of the system.
In order to resolve certain frequencies, e.g., in :math:`S(q, \omega)`, you need to make sure the time between two consecutive snapshots complies with the `Nyquist-Shannon sampling theorem <https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem>`_.
It is important to specify the correct units in your trajectory to the trajectory reader of :program:`dynasor`, such that they can be correctly converted to internal :program:`dynasor` units.
Trajectory formats
^^^^^^^^^^^^^^^^^^
:program:`dynasor` can parse multiple different trajectory formats.
It has native support for parsing trajectories in :program:`lammps` dump format :cite:`LAMMPSweb`.
Trajectories in `extended xyz format <https://github.com/libAtoms/extxyz>`_ can be read with the help of :program:`ase` :cite:`HjorthLarsen2017`.
Note that velocities are stored as arrays with the name ``vel`` as written by :program:`gpumd` (https://gpumd.org/).
The reader is multithreaded.
:program:`dynasor` also interfaces with python package `MDAnalysis <https://www.mdanalysis.org/>`_, which allows users to parse a vast set of different trajectory formats.
Atomic types
------------
If you want to calculate the correlation functions for each atom type (or different groupings of atoms) separately you need to provide the :attr:`atomic_indices` parameter in the form of a dictionary, for example, ::
atomic_indices['Cs'] = [0, 3, 6, 9]
atomic_indices['Br'] = [1, 2, 4, 5, 7, 8, 10, 11]
Alternativley, you can provide an index file (`Gromacs ndx-style <https://manual.gromacs.org/current/reference-manual/file-formats.html#ndx>`_) to group atoms. For some trajectory formats, it is also possible to read the :attr:`atomic_indices` directly from the trajectory by setting it to ``'read_from_trajectory'``.
If :attr:`atomic_indices` is left unspecified all atoms are assumed to be of the same species.
For example, in the case of a small system containing eight water molecules one could employ the following index file ::
[ Hydrogen ]
2 5 8 11 14 17 20 23
3 6 9 12 15 18 21 24
[ Oxygen ]
1 4 7 10 13 16 19 22
Here, the numbers following each tag (``[...]``) specify the atomic indices that belong to the respective species.
.. index:: Parameters; q-space sampling
q-space sampling parameters
---------------------------
:program:`dynasor` can either directly evaluate the correlation function at certain :math:`\boldsymbol{q}`-points (typically desirable for solid or crystalline materials) or average the correlation functions over :math:`\boldsymbol{q}`-points with similar length.
The latter means constructing spherical average over :math:`\boldsymbol{q}`-points which (often) requires thousands of :math:`\boldsymbol{q}`-points to be sampled to get good statistics.
Note that :math:`\boldsymbol{q}`-points in :program:`dynasor` are not in reduced coordinates but given as Cartesian coordinates (in inverse Å) and should include the factor :math:`2\pi`.
Spherical q-point sampling and averaging
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
:program:`dynasor` provides functionality to generate :math:`\boldsymbol{q}`-points on a 3D grid with the spacing set by the simulation cell up to to the maximum radius set by :attr:`q_max` in units of rad/Å.
The resulting correlation functions can then for example be averaged over :math:`\boldsymbol{q}`-points that fall into the same radial bin, where the parameter :attr:`q_bins` determines the number of radial bins.
The number of :math:`\boldsymbol{q}`-points can grow very large using a 3D grid, and hence slow down calculations.
The maximum number of :math:`\boldsymbol{q}`-points to scan is specified by :attr:`max_points`, which allows one to set an approximate upper limit for the number of :math:`\boldsymbol{q}`-points.
If the number of :math:`\boldsymbol{q}`-points on the grid exceeds :attr:`max_points` :math:`\boldsymbol{q}`-points will be randomly removed in such a way that the number of :math:`\boldsymbol{q}`-points inside each bin is roughly uniformly distributed.
.. index:: Parameters; time sampling
Time sampling parameters
------------------------
The two main time sampling parameters are the number of time lags :attr:`window_size` and maximum number of frames :attr:`frame_stop`.
The window size is only relevant for dynamic correlation functions but not for static ones.
The maximum number of frames sets a limit on how many frames are read from the trajectory file.
The window size is the longest time for which correlations are computed and is given in terms of the number of frames.
Additionally, you can specify :attr:`frame_step` to only use every :attr:`frame_step`-th snapshot in the trajectory.
For example if :attr:`frame_step=1` every snapshot will be used while for :attr:`frame_step=2` every other snapshot is used.
The :attr:`window_step` parameter controls how far the window is moved each time, with a default value of 1.
A larger value may lead to the calculation running faster but also possibly less converged correlation functions.
Note that if :attr:`window_step` is larger than the time window some snapshots in the trajectory will be skipped entirely.
To demonstrate how :attr:`frame_step` and :attr:`window_step` interact consider a trajectory with snapshots indexed from 1 to 100, and with :attr:`window_size=5`.
The following snapshots will enter in the two first windows for which the correlation function is computed::
frame_step = 1
window_step = 1
window 0 : snapshots [0, 1, 2, 3, 4]
window 1 : snapshots [1, 2, 3, 4, 5]
frame_step = 2
window_step = 1
window 0 : snapshots [0, 2, 4, 6, 8]
window 1 : snapshots [2, 4, 6, 8, 10]
frame_step = 1
window_step = 2
window 0 : snapshots [0, 1, 2, 3, 4]
window 1 : snapshots [2, 3, 4, 5, 6]
frame_step = 2
window_step = 2
window 0 : snapshots [0, 2, 4, 6, 8]
window 1 : snapshots [4, 6, 8, 10, 12]
Lastly, it is required to set the time difference, :attr:`dt`, between two consecutive snapshots in the trajectory in femto seconds in order to get everything in the correct units.
Note that you should not change :attr:`dt` when you change :attr:`frame_step`.
.. figure:: _static/window_sampling.png
:scale: 16 %
:align: center
Schematic illustration of sliding time-window averaging over a trajectory.
Incoherent intermediate scattering function
--------------------------------------------------
The calculation of the incoherent part of the intermediate scattering function is turned off by default.
It can be activated by setting :attr:`calculate_incoherent` to ``True`` when using the python interface and ``--calculate-incoherent`` when using the command line interface.
Post-processing
---------------
After the correlation functions have been calculated the results are stored as `numpy <https://numpy.org/>`_ array in the :class:`Sample <sample>` object.
There are multiple post-processing functions available, such as conducting spherical averages over :math:`\boldsymbol{q}`-points and weighting with atomic form factors.
Limitations
-----------
Sometimes problems can occur due to known FFT features such as folding and aliasing.
We have found the `filon` integration to work well for most applications, but sometimes it might be useful to apply FFT filters/windows to get a more smooth spectrum in the frequency domain.
The concept of :math:`\boldsymbol{q}`-points becomes muddy if the simulation box is changing (which is the case, e.g., for NPT MD simulations).
Differences in positions between two frames with different boxes are not uniquely defined.
Therefore only constant volume MD trajectories are supported in :program:`dynasor`.
|
