Tools#

A number of utility functions, for example for dealing with autocorrelation functions, Fourier transforms, and smoothing.

dynasor.tools.acfs.compute_acf(Z, delta_t=1.0, method='scipy')[source]#

Computes the autocorrelation function (ACF) for a one-dimensional signal \(Z\) in time as

\[\text{ACF}(\tau) = \frac{\left < Z(t) Z^*(t+\tau) \right >}{\left < Z(t) Z^*(t) \right >}\]

Here, only the real part of the ACF is returned since if \(Z\) is complex the imaginary part should average out to zero for any stationary signal.

Parameters:

Z (ndarray[tuple[Any, ...], dtype[float]]) – Complex time signal.
delta_t (Optional[float]) – Spacing in time between two consecutive values in \(Z\).
method (Optional[str]) – Implementation to use; possible values: numpy and scipy (default and usually faster).

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.acfs.fermi_dirac(t, t_0, t_width)[source]#

Evaluates a Fermi-Dirac-like function in time \(f(t)\), which can be applied to an auto-correlation function (ACF) in time to artificially dampen it, i.e., forcing it to go to zero for long times without affecting the short-time correlations too much.

\[f(t) = \frac{1}{\exp{[(t-t_0)/t_\mathrm{width}}] + 1}\]

Parameters:

t (ndarray[tuple[Any, ...], dtype[float]]) – Time array.
t_0 (float) – Starting time for decay.
t_width (float) – Width of the decay.

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.acfs.gaussian_decay(t, t_sigma)[source]#

Evaluates a gaussian distribution in time \(f(t)\), which can be applied to an ACF in time to artificially damp it, i.e., forcing it to go to zero for long times.

\[f(t) = \exp{\left [-\frac{1}{2} \left (\frac{t}{t_\mathrm{sigma}}\right )^2 \right ] }\]

Parameters:

t (ndarray[tuple[Any, ...], dtype[float]]) – Time array.
t_sigma (float) – Width (standard deviation of the gaussian) of the decay.

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.acfs.psd_from_acf(acf, dt=1, even=True)[source]#

Computes the power spectral density (PSD) from an auto-correlation function (ACF).

Let x(t) be a time signal and define its auto-correlation function

C(τ) = ⟨ x(t) x(t+τ) ⟩

where ⟨·⟩ denotes a time or ensemble average. According to the Wiener–Khinchin theorem, the power spectral density (PSD) of x(t) is the Fourier transform of its auto-correlation function:

S(ω) = ∫ C(τ) e^{-i ω τ} dτ

This function computes the PSD by performing a discrete Fourier transform of the provided ACF.

Parameters:

acf (ndarray[tuple[Any, ...], dtype[float]]) – The auto-correlation function C(τ) evaluated at discrete time lags.
dt (Optional[float]) – Time spacing between consecutive samples in the ACF.
even (Optional[bool]) – Whether the ACF is assumed to be even in time. If True, the ACF is mirrored to construct a two-sided correlation function before computing the Fourier transform.

Returns:

Angular frequencies ω and corresponding PSD values.

Return type:

freqs, psd

dynasor.tools.acfs.psd_from_acf_2d(acf, dt=1, even=True)[source]#

Computes the power spectral density (PSD) from an auto-correlation function (ACF).

This is the 2D analogue of psd_from_acf, where the ACF is provided for multiple signals, e.g. F(q, t). The Fourier transform is performed along the last axis.

Parameters:

acf (ndarray[tuple[Any, ...], dtype[float]]) – The ACF as a 2D array with shape (N, Nt), where the last axis corresponds to time lags.
dt (Optional[float]) – The time step between samples.
even (Optional[bool]) – Whether the ACF is even in time and should be mirrored before computing the Fourier transform.

Returns:

Angular frequencies ω and corresponding PSD values with shape (N, Nω).

Return type:

freqs, psd

dynasor.tools.acfs.psd_from_time_signal(x, dt=1, complex=False)[source]#

Computes the power spectral density (PSD) directly from a time signal.

Let x(t) be a time signal. The power spectral density (PSD) can be computed directly from the Fourier transform of the signal:

S(ω) = |F[x(t)]|²

where F denotes the Fourier transform.

This function computes the PSD by performing a discrete Fourier transform of the time signal and taking the squared magnitude of the transform.

This is equivalent to computing the PSD from the Fourier transform of the auto-correlation function.

Parameters:

x (ndarray[tuple[Any, ...], dtype[float]]) – Time signal as a 1D array.
dt (Optional[float]) – Time spacing between consecutive samples.
complex (Optional[bool]) – Whether the time signal is complex. If False (default), the signal is assumed to be real and the PSD is computed using a real FFT (rfft), returning only non-negative frequencies. If True, the full FFT (fft) is used and both positive and negative frequencies are returned.

Returns:

Angular frequencies ω and corresponding PSD values.

Return type:

freqs, psd

dynasor.tools.acfs.smoothing_function(data, window_size, window_type='hamming')[source]#

Smoothing function for 1D arrays. This functions employs the pandas rolling window average function.

Parameters:

data (ndarray[tuple[Any, ...], dtype[float]]) – 1D data array.
window_size (int) – The size of smoothing/smearing window.
window_type (Optional[str]) – What type of window-shape to use, e.g. 'blackman', 'hamming', 'boxcar' (see pandas and scipy documentaiton for more details).

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.damped_harmonic_oscillator.acf_position_dho(t, w0, gamma, A=1.0)[source]#

Calculate the damped damped harmonic oscillator (DHO) autocorrelation function for the position. The definition of this function can be found in the dynasor documentation. # noqa

Parameters:

t (Union[float, ndarray[tuple[Any, ...], dtype[float]]]) – Time, usually as an array.
w0 (float) – Natural angular frequency of the DHO.
gamma (float) – Damping of DHO.
A (Optional[float]) – Amplitude of the DHO.

Return type:

Union[float, ndarray[tuple[Any, ...], dtype[float]]]

dynasor.tools.damped_harmonic_oscillator.acf_velocity_dho(t, w0, gamma, A=1.0)[source]#

Calculate the damped damped harmonic oscillator (DHO) autocorrelation function for the velocity. The definition of this function can be found in the dynasor documentation. # noqa

Parameters:

t (Union[float, ndarray[tuple[Any, ...], dtype[float]]]) – Time, usually as an array.
w0 (float) – Natural angular frequency of the DHO.
gamma (float) – Damping of DHO.
A (Optional[float]) – Amplitude of the DHO.

Return type:

Union[float, ndarray[tuple[Any, ...], dtype[float]]]

dynasor.tools.damped_harmonic_oscillator.psd_position_dho(w, w0, gamma, A=1.0)[source]#

Calculate the power spectral density (PSD) function for the damped harmonic oscillator (DHO), i.e., the Fourier transform of the autocorrelation function) for the position.

The definition of this function can be found in the dynasor documentation. # noqa

Parameters:

w (Union[float, ndarray[tuple[Any, ...], dtype[float]]]) – Angular frequency, usually as an array.
w0 (float) – Natural angular frequency of the DHO.
gamma (float) – Damping of DHO.
A (Optional[float]) – Amplitude of the DHO.

Return type:

Union[float, ndarray[tuple[Any, ...], dtype[float]]]

dynasor.tools.damped_harmonic_oscillator.psd_velocity_dho(w, w0, gamma, A=1.0)[source]#

Calculate the power spectral density (PSD) function for the damped harmonic oscillator (DHO), i.e., the Fourier transform of the autocorrelation function) for the position.

The definition of this function can be found in the dynasor documentation. # noqa

Parameters:

w (Union[float, ndarray[tuple[Any, ...], dtype[float]]]) – Angular frequency, usually as an array.
w0 (float) – Natural angular frequency of the DHO.
gamma (float) – Damping of DHO.
A (Optional[float]) – Amplitude of the DHO.

Return type:

Union[float, ndarray[tuple[Any, ...], dtype[float]]]

dynasor.tools.structures.align_structure(atoms, atol=1e-05)[source]#

Rotates and realigns a structure such that * the first cell vector points along the x-directon * the second cell vector lies in the xy-plane

Note that this function modifies the atoms object in place.

Parameters:

atoms (Atoms) – Input structure to be rotated aligned with the x,y,z coordinte system.
atol (Optional[float]) – Absolute tolerance used for sanity checking the cell.

Return type:

None

dynasor.tools.structures.find_permutation(atoms, atoms_ref)[source]#

Returns the best permutation of atoms for mapping one configuration onto another.

Parameters:

atoms (Atoms) – Configuration to be permuted.
atoms_ref (Atoms) – Configuration onto which to map.

Return type:

list[int]

Example

After obtaining the permutation via p = find_permutation(atoms1, atoms2) the reordered structure atoms1[p] will give the closest match to atoms2.

dynasor.tools.structures.get_P_matrix(c, S)[source]#

Returns the P matrix, i.e., the 3x3 integer matrix \(P\) that satisfies

\[P c = S\]

Here, \(c\) is the primitive cell metric and \(S\) is the supercell metric as row vectors. Note that the above condition is equivalent to:

\[c^T P^T = S^T\]

Parameters:

c (ndarray[tuple[Any, ...], dtype[float]]) – Cell metric of the primitive structure.
S (ndarray[tuple[Any, ...], dtype[float]]) – Cell metric of the supercell.

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.structures.get_displacements(atoms, atoms_ideal, check_mic=True, cell_tol=0.0001)[source]#

Returns the the smallest possible displacements between a displaced configuration relative to an ideal (reference) configuration.

Parameters:

atoms (Atoms) – Structure with displaced atoms.
ideal – Ideal configuration relative to which displacements are computed.
check_mic (Optional[bool]) – Whether to check minimum image convention.
cell_tol (Optional[float]) – Cell tolerance; if cell missmatch more than tol value error is raised.

Return type:

ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.structures.get_displacements_from_u(u, cell, check_mic=True)[source]#

wraps displacements using mic

Return type:: ndarray[tuple[Any, ...], dtype[float]]

dynasor.tools.structures.get_offset_index(primitive, supercell, tol=0.01, wrap=True)[source]#

Returns the basis index and primitive cell offsets for a supercell.

This implementation uses a simple iteration procedure that should be fairly quick. If more stability is needed consider the following approach:

find the P-matrix: P = ideal.cell @ prim.cell_inv.T
compensate for strain: P *= len(ideal)/len(prim)/det(P)
generate the reference structure: ref_atoms = make_supercell(round(P), prim)
find the assignment using ref_atoms via the Hungarian algorithm using the mic distances

Parameters:

primitive (Atoms) – Primitive cell.
supercell (Atoms) – Some ideal repetition of the primitive cell.
tol (Optional[float]) – Tolerance length parameter. Increase to allow for slgihtly rattled or strained cells.
wrap (Optional[bool]) – It might happen that the ideal cell boundary cuts through a unit cell whose lattice points lie inside the ideal cell. If there is a basis, an atom belonging to this unit cell might get wrapped while another is not. Then the wrapped atom now belongs to a lattice point outside the P matrix so to say. This would result in more lattice points than expected from N_unit = len(ideal)/len(prim).

Return type:

tuple[ndarray[tuple[Any, ...], dtype[float]], ndarray[tuple[Any, ...], dtype[float]]]

Returns:

offsets – The lattice points as integers in (N, 3)-array.
index – The basis indices as integers in (N,)-array.