Tools#

Data analysis#

calorine.tools.analyze_data(data, max_lag=None)[source]#

Carries out an extensive analysis of the data series.

Parameters:

data (ndarray) – data series to compute autocorrelation function for
max_lag (int) – maximum lag between two data points, used for computing autocorrelation

Returns:

calculated properties of the data including, mean, standard deviation, correlation length and a 95% error estimate.

Return type:

dict

calorine.tools.get_autocorrelation_function(data, max_lag=None)[source]#

Returns autocorrelation function.

The autocorrelation function is computed using pandas.Series.autocorr.

Parameters:

data (ndarray) – data series to compute autocorrelation function for
max_lag (int) – maximum lag between two data points

Return type:

ndarray

Returns:

calculated autocorrelation function

calorine.tools.get_correlation_length(data)[source]#

Returns estimate of the correlation length of data.

The correlation length is taken as the first point where the autocorrelation functions is less than \(\exp(-2)\). If the correlation function never drops below \(\exp(-2)\) np.nan is returned.

If the correlation length cannot be computed since the auto-correlation function is unconverged the function returns None.

Parameters:: data (ndarray) – data series for which to the compute autocorrelation function
Return type:: Optional[int]
Returns:: correlation length

calorine.tools.get_error_estimate(data, confidence=0.95)[source]#

Returns estimate of standard error \(\mathrm{error}\) with confidence interval.

\[\mathrm{error} = t_\mathrm{factor} * \mathrm{std}(\mathrm{data}) / \sqrt{N_s}\]

where \(t_{factor}\) is the factor corresponding to the confidence interval and \(N_s\) is the number of independent measurements (with correlation taken into account).

If the correlation length cannot be computed since the auto-correlation function is unconverged the function returns None.

Parameters:: data (ndarray) – data series for which to estimate the error
Return type:: Optional[float]
Returns:: error estimate

Property calculation#

calorine.tools.batch_predict_properties(structures, model, command=None, directory=None)[source]#

Evaluates NEP model properties for a list of structures in a single pass, using the prediction mode of the nep executable (see here). This is substantially faster than evaluating structures one at a time with CPUNEP or GPUNEP, since all structures are transferred to the GPU in a single pass.

Parameters:

structures (List[Atoms]) – Structures for which to evaluate properties.
model (Union[str, Model]) – Either a path to a NEP model in nep.txt format, or a Model object.
command (Optional[str]) – Command used to invoke the nep executable. Default: nep, or the value of the CALORINE_NEP_COMMAND environment variable if set.
directory (Optional[str]) – Directory in which to run nep. If None, a temporary directory is created and removed once the calculation is finished. If specified, the directory is created if needed and is not deleted afterward, which is useful for debugging or for further analysis of the raw nep output files (e.g. via read_structures).

Returns:

A new list of Atoms objects, in the same order as structures, each with a SinglePointCalculator attached exposing the predicted properties in the standard way (energy, forces, and stress, plus charges and born_effective_charges for qNEP models, or dipole/polarizability for TNEP models). The input structures are not modified.

Return type:

List[Atoms]

calorine.tools.get_elastic_stiffness_tensor(structure, clamped=False, epsilon=0.001, **kwargs)[source]#

Calculate and return the elastic stiffness tensor in units of GPa for the given structure in Voigt form.

Parameters:

structure (Atoms) – input structure; should be fully relaxed
clamped (bool) – if False (default) return the relaxed elastic stiffness tensor; if True return the clamped ion elastic stiffness tensor
epsilon (float) – magnitude of the applied strain
kwargs – keyword arguments forwarded to the relax_structure function used for relaxing the structure when computing the relaxed stiffness tensor; it should not be necessary to change the default for the vast majority of use cases; use with care

Return type:

ndarray

calorine.tools.get_entropy(descriptors, width, use_tqdm=True, device='cuda', block=1024, dtype=torch.float32, eps=1e-12)[source]#

Computes an estimate for the information entropy \(H(\mathbf{X})\) of a set of descriptors \(\mathbf{X}\). The estimate is described in [Nat. Comm. 16, 4014 (2025)](https://doi.org/10.1038/s41467-025-59232-0) and given by

\[\mathcal{H}(\{\mathbf{X}\}) = -\frac{1}{n} \sum_{i=1}^{n} p_i\]

where

\[p_i = \log \left[ \frac{1}{n} \sum_{j=1}^{n} K_h(\mathbf{X}_i, \mathbf{X}_j) \right]\]

with a Gaussian kernel

\[K_h(\mathbf{X}_i, \mathbf{X}_j) = \exp\!\left( -\frac{\lVert \mathbf{X}_i - \mathbf{X}_j \rVert^2}{2h^2}. \right)\]

The calculation is done via torch if the library has been installed, and numpy otherwise. When using torch the calculation is run via CUDA. The latter behavior can be controlled using the device argument.

Parameters:

descriptors (ndarray | Tensor) – The set of descriptors \(\mathbf{X}\) for which to evaluate to the entropy. Typically each row corresponds to one atom and the columns correspond to the different descriptor components.
width (float) – Width \(h\) of the Gaussian kernel.
use_tqdm (bool) – Use tqdm to show a progress bar. Note that this requires tqdm to be installed.
block (int) – In order to limit the memory needs, the kernel density estimate matrix is handled in blocks. This parameter controls the size of each block. Smaller numbers imply a smaller memory footprint.
eps (float) – Smallest (absolute) permissible value.
device (str | device) – Device to use for calculation. The documentation of [torch.device](https://docs.pytorch.org/docs/stable/tensor_attributes.html#torch.device) provides more information. Only used when pytorch is available.
dtype (dtype) – Floating point precision used for the computation. The documentation of [torch.dtype](https://docs.pytorch.org/docs/stable/tensor_attributes.html) provides an overview. Only used when pytorch is available.

Return type:

tuple[float, ndarray]

Returns:

A tuple comprising the total entropy \(H(\mathbf{X})\) and the entropy contributions
\(p_i\) from each row in the input descriptor matrix.

calorine.tools.get_force_constants(structure, calculator, supercell_matrix, kwargs_phonopy={}, kwargs_generate_displacements={})[source]#

Calculates the force constants for a given structure using phonopy, which needs to be cited if this function is used for generating data for publication. The function returns a Phonopy object that can be used to calculate, e.g., the phonon dispersion, the phonon density of states as well as related quantities such as the thermal displacements and the free energy.

Parameters:

structure (Atoms) – structure for which to compute the phonon dispersion; usually this is a primitive cell
calculator (SinglePointCalculator) – ASE calculator to use for the calculation of forces
supercell_matrix (ndarray) – specification of supercell size handed over to phonopy; should be a tuple of three values or a matrix
kwargs_phonopy (Dict[str, Any]) – Expert option: keyword arguments used when initializing the Phonopy object; this includes, e.g., the tolerance used when determining the symmetry (symprec) and parameters for the non-analytical corrections (nac_params)
kwargs_generate_displacements (Dict[str, Any]) – Expert option: keyword arguments to be handed over to the generate_displacements method; this includes in particular the distance keyword, which specifies the magnitude of the atomic displacement imposed when calculating the force constant matrix

Return type:

Phonopy

calorine.tools.get_rtc_from_hac(hac, V, T, dt)[source]#

Returns the running thermal conductivity (RTC) in W/m/K using a heat auto-correlation (HAC) as input.

Parameters:

hac (ndarray) – The HAC \(\langle j(t)j(0)\rangle\) in units of eV³/amu as given by GPUMD
V (float) – Volume of cell in Å³
T (float) – Temperature in Kelvin
dt (float) – Time step in the HAC in ps running thermal conductivity in W/m/K

Return type:

ndarray

Spectra#

calorine.tools.apply_kramers_kronig(df, method='vectorized')[source]#

Apply the Kramers-Kronig relation to add the real part of the dielectric function.

Takes the output of get_dielectric_function(), computes \(\epsilon_1(\omega)\) from \(\epsilon_2(\omega)\) via the Kramers-Kronig relation, and returns the DataFrame with new epsilon_real* columns appended.

Each epsilon_imag{suffix} column in df produces a corresponding epsilon_real{suffix} column (e.g. epsilon_imag_xx → epsilon_real_xx).

Parameters:

df (DataFrame) – DataFrame as returned by get_dielectric_function() or get_ir_spectrum(); must contain angular_frequency and at least one column whose name starts with epsilon_imag.
method (str) – Integration method. 'vectorized' (default) uses an exact trapezoid rule over an \(n \times n\) matrix (\(\mathcal{O}(n^2)\) time and memory); 'fft' uses an \(\mathcal{O}(n \log n)\) Hilbert transform and is faster for large arrays but approximates the principal-value integral. The two methods agree to within a few percent on typical grids.

Returns:

Input DataFrame with additional epsilon_real* columns.

Return type:

DataFrame

calorine.tools.apply_quantum_correction(df, temperature, column, order='first', force=False)[source]#

Apply a quantum correction to a classically computed spectrum.

Classical MD underestimates spectral intensities because it samples the classical Boltzmann distribution rather than the Bose-Einstein distribution. The correction factor depends on the scattering order and mode type [Cardona1982] [Rosander2025].

Parameters:

df (DataFrame) – DataFrame with an angular_frequency column (THz) and the column to correct.
temperature (float) – Temperature in K.
column (str) – Name of the spectral column to correct (e.g. 'raman_isotropic', 'ir_intensity', 'epsilon_imag').
order (str) –
Correction type:

'first'
First-order scattering (IR absorption and first-order Raman): \(f = \beta\hbar\omega / (1 - e^{-\beta\hbar\omega})\).

'overtone'
Second-order overtone (same mode, \(\omega_1 = \omega/2\)): \(f = [y/(1 - e^{-y})]^2 (2 - e^{-y})\) with \(y = \beta\hbar\omega/2\).

'combination'
Second-order combination band upper bound (\(\omega_1 = \omega_2 = \omega/2\)): \(f = [y/(1 - e^{-y})]^2\) with \(y = \beta\hbar\omega/2\).
force (bool) – If True, overwrite an existing column + '_qm' column. Default False raises ValueError to prevent accidental double-correction.

Returns:

Copy of df with a new column column + '_qm' containing the quantum-corrected intensities. The input is not mutated.

Return type:

DataFrame

References

[Cardona1982]

M. Cardona, Resonance phenomena, in Topics in Applied Physics, Vol. 50, edited by M. Cardona and G. Güntherodt (Springer, Berlin, 1982).

[Rosander2025]

Rosander et al., Phys. Rev. B 111, 064107 (2025), Supp. Eqs. S5, S7, S10.

calorine.tools.get_dielectric_function(signal, volume, temperature, column='dP', dt=None, t_sigma=None, max_lag=0.25, return_real_part=True, kk_method='vectorized')[source]#

Compute the dielectric function from a three-component polarization time series.

Computes all six unique Voigt components of \(\epsilon_2^{\alpha\beta}(\omega)\) (xx, yy, zz, yz, xz, xy) from the three-component time series written by GPUMD. The real part \(\epsilon_1^{\alpha\beta}(\omega)\) is obtained via the Kramers-Kronig relation when return_real_part is True (the default).

The imaginary part is related to the symmetrized cross-correlation via

\[\epsilon_2^{\alpha\beta}(\omega) = \frac{\beta}{\epsilon_0 \omega V} \,\mathrm{Re}\!\left[\int_0^\infty \langle \dot{P}_\alpha(0)\,\dot{P}_\beta(t) \rangle e^{-i\omega t}\,dt\right],\]

where \(\beta = 1/(k_\mathrm{B}T)\) and \(V\) is the cell volume.

Parameters:

signal (DataFrame) – DataFrame as returned by read_dpdt() or read_dipole(). The columns to use are selected via column.
volume (float) – Simulation cell volume in Å³.
temperature (float) – Temperature in K.
column (str) –
Selects which columns of signal to use and implies whether the input is a time derivative:

'dP' (default)
Uses dPx, dPy, dPz from read_dpdt() (qNEP); input is \(\dot{\mathbf{P}}(t)\) in e·Å/fs.

'P'
Uses Px, Py, Pz from read_dpdt() (qNEP); input is the polarization \(\mathbf{P}(t)\) in e·Å.

'mu'
Uses mu_x, mu_y, mu_z from read_dipole() (TNEP); input is the dipole moment \(\boldsymbol{\mu}(t)\) in e·Å.
dt (float) – Time between consecutive frames in ps. Auto-extracted from the time column when present (i.e. for read_dpdt() output). Must be supplied explicitly for step-indexed inputs such as read_dipole() output.
t_sigma (float) – Width of the Gaussian window applied to the ACF in ps. None uses no windowing.
max_lag (float) – Maximum ACF length as a fraction of the trajectory length; must be in (0, 1]. Defaults to 0.25 (25 %).
return_real_part (bool) – When True (default), the real part \(\epsilon_1^{\alpha\beta}(\omega)\) is computed via the Kramers-Kronig relation and appended as epsilon_real_{xx,yy,zz,yz,xz,xy} columns.
kk_method (str) – Integration method passed to apply_kramers_kronig() when return_real_part is True. 'vectorized' (default) uses an exact trapezoid rule; 'fft' uses a faster Hilbert transform approximation. Ignored when return_real_part is False.

Returns:

Contains angular_frequency (THz), wavenumber_invcm (\(\mathrm{cm}^{-1}\)), epsilon_imag_{xx,yy,zz,yz,xz,xy} (dimensionless, imaginary dielectric tensor components), and conductivity_{xx,yy,zz,yz,xz,xy} (S/m). When return_real_part is True, epsilon_real_{xx,yy,zz,yz,xz,xy} columns are appended.

Return type:

DataFrame

calorine.tools.get_ir_spectrum(signal, volume=None, temperature=None, column='dP', dt=None, polarization=None, t_sigma=None, max_lag=0.25)[source]#

Compute the IR spectrum from a time series of dipole moments or polarizations.

Parameters:

signal (DataFrame) – DataFrame as returned by read_dpdt() or read_dipole(). The columns to use are selected via column.
volume (float) – Simulation cell volume in Å³. Required for extended systems; when provided the function returns the imaginary dielectric function and conductivity. Pass None for molecules to obtain an unnormalized line shape in arbitrary units.
temperature (float) – Temperature in K. Required when volume is not None.
column (str) –
Selects which columns of signal to use and implies whether the input is a time derivative:

'dP' (default)
Uses dPx, dPy, dPz from read_dpdt() (qNEP); signal is \(\mathrm{d}P/\mathrm{d}t\) in e·Å/fs.

'P'
Uses Px, Py, Pz from read_dpdt() (qNEP); signal is the polarization \(P(t)\) in e·Å.

'mu'
Uses mu_x, mu_y, mu_z from read_dipole() (TNEP); signal is the dipole moment \(\mu(t)\) in e·Å.
dt (float) – Time between consecutive frames in ps. Auto-extracted from the time column when present (i.e. for read_dpdt() output). Must be supplied explicitly for step-indexed inputs such as read_dipole() output.
polarization (ndarray) – Unit vector (3,) defining the electric-field polarization direction. When given, the spectrum is computed from the projected signal \(s(t) = \hat{n} \cdot \mathrm{signal}(t)\) instead of the isotropic average.
t_sigma (float) – Width of the Gaussian window applied to the ACF in ps. None uses no windowing.
max_lag (float) – Maximum ACF length as a fraction of the trajectory length; must be in (0, 1]. Defaults to 0.25 (25 %).

Returns:

Always contains angular_frequency (THz) and wavenumber_invcm (\(\mathrm{cm}^{-1}\)). When volume and temperature are given: epsilon_imag (dimensionless) and conductivity (S/m). When volume is None: ir_intensity (arbitrary units, proportional to \(\mathrm{PSD}[\dot{\mu}]\)).

Return type:

DataFrame

calorine.tools.get_raman_spectrum(dt, polarizability, polarization_in=None, polarization_out=None, t_sigma=None, max_lag=0.25)[source]#

Compute the Raman spectrum from a time series of polarizability tensors.

Parameters:

dt (float) – Time between consecutive frames in ps.
polarizability (DataFrame) – DataFrame with columns xx, yy, zz, xy, yz, xz containing the polarizability \(\alpha\) (molecules) or susceptibility \(\chi\) (extended systems), as returned by read_polarizability(). Units are those of the TNEP training data (typically bohr³).
polarization_in (ndarray) – Unit vector (3,) for the polarization of the incoming light. If both polarization_in and polarization_out are given, the polarization-resolved intensity \(I(\omega) \propto \mathrm{FT}[\langle s(0)s(t)\rangle]\) with \(s(t) = \hat{n}^\mathrm{out} \cdot \alpha(t) \cdot \hat{n}^\mathrm{in}\) is added as the column raman_polarized.
polarization_out (ndarray) – Unit vector (3,) for the polarization of the outgoing (scattered) light.
t_sigma (float) – Width of the Gaussian window applied to the ACF in ps. None uses no windowing.
max_lag (float) – Maximum ACF length as a fraction of the trajectory length; must be in (0, 1]. Defaults to 0.25 (25 %).

Returns:

Always contains angular_frequency (THz), wavenumber_invcm (\(\mathrm{cm}^{-1}\)), raman_isotropic (proportional to \(\mathrm{FT}[\langle\gamma(0)\gamma(t)\rangle]\)), and raman_anisotropic (proportional to \(\mathrm{FT}[\langle\mathrm{Tr}[\beta(0)\beta(t)]\rangle]\)). If both polarization vectors are given, also contains raman_polarized.

Return type:

DataFrame

Structure manipulation#

calorine.tools.relax_structure(structure, fmax=0.001, steps=500, minimizer='bfgs', constant_cell=False, constant_volume=False, scalar_pressure=0.0, **kwargs)[source]#

Relaxes the given structure.

Parameters:

structure (Atoms) – Atomic configuration to relax.
fmax (float) – Stop relaxation if the absolute force for all atoms falls below this value.
steps (int) – Maximum number of relaxation steps the minimizer is allowed to take.
minimizer (str) – Minimizer to use; possible values: ‘bfgs’, ‘lbfgs’, ‘fire’, ‘gpmin’, ‘bfgs-scipy’.
constant_cell (bool) – If True do not relax the cell or the volume.
constant_volume (bool) – If True relax the cell shape but keep the volume constant.
kwargs – Keyword arguments to be handed over to the minimizer; possible arguments can be found in the ASE documentation
scalar_pressure (float) – External pressure in GPa.

Return type:

None

Symmetry analysis#

calorine.tools.get_primitive_structure(structure, no_idealize=True, to_primitive=True, symprec=1e-05)[source]#

Returns the primitive structure using spglib. This is a convenience interface to the standardize_cell() function of spglib that works directly with ase Atoms objects.

Parameters:

structure (Atoms) – Input atomic structure.
no_idealize (bool) – If True lengths and angles are not idealized.
to_primitive (bool) – If True convert to primitive structure.
symprec (float) – Tolerance imposed when analyzing the symmetry.

Return type:

Atoms

calorine.tools.get_spacegroup(structure, symprec=1e-05, angle_tolerance=-1.0, style='international')[source]#

Returns the space group of a structure using spglib. This is a convenience interface to the get_spacegroup() function of spglib that works directly with ase Atoms objects.

Parameters:

structure (Atoms) – Input atomic structure.
symprec (float) – Tolerance imposed when analyzing the symmetry.
angle_tolerance (float) – Tolerance imposed when analyzing angles.
style (str) – Space group notation to be used. Can be 'international' for the interational tables of crystallography (ITC) style (Hermann-Mauguin and ITC number) or 'Schoenflies' for the Schoenflies notation.

Return type:

str

calorine.tools.get_wyckoff_sites(structure, map_occupations=None, symprec=1e-05, include_representative_atom_index=False)[source]#

Returns the Wyckoff symbols of the input structure. The Wyckoff labels can be conveniently attached as an array to the structure object as demonstrated in the examples section below.

By default the occupation of the sites is part of the symmetry analysis. If a chemically disordered structure is provided this will usually reduce the symmetry substantially. If one is interested in the symmetry of the underlying structure one can control how occupations are handled. To this end, one can provide the map_occupations keyword argument. The latter must be a list, each entry of which is a list of species that should be treated as indistinguishable. As a shortcut, if all species should be treated as indistinguishable one can provide an empty list. Examples that illustrate the usage of the keyword are given below.

Parameters:

structure (Atoms) – Input structure. Note that the occupation of the sites is included in the symmetry analysis.
map_occupations (List[List[str]]) – Each sublist in this list specifies a group of chemical species that shall be treated as indistinguishable for the purpose of the symmetry analysis.
symprec (float) – Tolerance imposed when analyzing the symmetry using spglib.
include_representative_atom_index (bool) – If True the index of the first atom in the structure that is representative of the Wyckoff site is included in the symbol. This is in particular useful in cases when there are multiple Wyckoff sites sites with the same Wyckoff letter.

Return type:

List[str]

Examples

Wyckoff sites of a hexagonal-close packed structure:

>>> from ase.build import bulk
>>> structure = bulk('Ti')
>>> wyckoff_sites = get_wyckoff_sites(structure)
>>> print(wyckoff_sites)
['2d', '2d']

The Wyckoff labels can also be attached as an array to the structure, in which case the information is also included when storing the Atoms object:

>>> from ase.io import write
>>> structure.new_array('wyckoff_sites', wyckoff_sites, str)
>>> write('structure.xyz', structure)

The function can also be applied to supercells:

>>> structure = bulk('GaAs', crystalstructure='zincblende', a=3.0).repeat(2)
>>> wyckoff_sites = get_wyckoff_sites(structure)
>>> print(wyckoff_sites)
['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']

Now assume that one is given a supercell of a (Ga,Al)As alloy. Applying the function directly yields much lower symmetry since the symmetry of the original structure is broken:

>>> structure.set_chemical_symbols(
...        ['Ga', 'As', 'Al', 'As', 'Ga', 'As', 'Al', 'As',
...         'Ga', 'As', 'Ga', 'As', 'Al', 'As', 'Ga', 'As'])
>>> print(get_wyckoff_sites(structure))
['8g', '8i', '4e', '8i', '8g', '8i', '2c', '8i',
 '2d', '8i', '8g', '8i', '4e', '8i', '8g', '8i']

Since Ga and Al occupy the same sublattice, they should, however, be treated as indistinguishable for the purpose of the symmetry analysis, which can be achieved via the map_occupations keyword:

>>> print(get_wyckoff_sites(structure, map_occupations=[['Ga', 'Al'], ['As']]))
['4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c',
 '4a', '4c', '4a', '4c', '4a', '4c', '4a', '4c']

If occupations are to ignored entirely, one can simply provide an empty list. In the present case, this turns the zincblende lattice into a diamond lattice, on which case there is only one Wyckoff site:

>>> print(get_wyckoff_sites(structure, map_occupations=[]))
['8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a',
 '8a', '8a', '8a', '8a', '8a', '8a', '8a', '8a']