mdse.md.resultMD

Post-processing and analysis of molecular dynamics simulation data.

This module defines the ResultMD class, which serves as a primary tool for analyzing the output of molecular dynamics simulations. It acts as a container for a sequence of simulation frames (ASE Atoms objects) and provides a rich set of methods to compute various physical properties from the trajectory data.

The ResultMD class can be instantiated directly with a list of frames or loaded from a .traj file using the from_file classmethod.

Key functionalities include the calculation of:

Structural and Transport Properties:
- Mean Squared Displacement (MSD)
- Self-diffusion coefficient
- Lindemann index
- Nearest-neighbor distance
Thermodynamic Properties:
- Isochoric heat capacity
- Isobaric specific heat and enthalpy
- Cohesive energy
- Temperature, pressure, and energy evolution over time
Vibrational Properties:
- Velocity Autocorrelation Function (VACF)
- Density of States (DOS)
- Debye temperature
Mechanical Properties:
- Elastic constants (C11, C12, C44)
- Bulk, Shear, and Young’s moduli

Classes

ResultMD(data[, conv_crystal, calc_params])

Class representing the result of a molecular dynamics (MD) simulation.

class mdse.md.resultMD.ResultMD(data, conv_crystal=None, calc_params=None)

Bases: object

Class representing the result of a molecular dynamics (MD) simulation.

This class stores frames from a simulation and provides methods to calculate and visualize the mean squared displacement (MSD).

calc_C11(crystal_equil, epsilon=0.01)

Calculate the C11 elastic constant.

C11 is calculated by applying a uniaxial strain \(\epsilon_{xx}\) and using a finite difference approximation for the second derivative of energy.

Parameters:

crystal_equil (ase.Atoms) – A relaxed equilibrium structure. A calculator will be attached inside the function.
epsilon (float, optional) – Magnitude of the applied strain. Defaults to 0.01.

Returns:

The C11 elastic constant in Pascals.

Return type:

float

calc_C12(crystal_equil, epsilon=0.01)

Calculate the C12 elastic constant.

C12 is calculated by applying a biaxial strain (\(\epsilon_{xx} = \epsilon\), \(\epsilon_{yy} = -\epsilon\)) and using a finite difference approximation.

Parameters:

crystal_equil (ase.Atoms) – A relaxed equilibrium structure. A calculator will be attached inside the function.
epsilon (float, optional) – Magnitude of the applied strain. Defaults to 0.01.

Returns:

The C12 elastic constant in Pascals.

Return type:

float

calc_C44(crystal_equil, gamma)

Calculate the C44 elastic constant.

C44 is calculated by applying a shear strain (\(\epsilon_{xy}\)) and using a finite difference approximation.

Parameters:

crystal_equil (ase.Atoms) – A relaxed equilibrium structure. A calculator will be attached inside the function.
gamma (float) – Magnitude of the applied shear strain.

Returns:

The C44 elastic constant in Pascals.

Return type:

float

calc_bulk_modulus(strain=0.01)

Calculate the bulk modulus (B) for a cubic crystal.

The bulk modulus is calculated from the elastic constants C11 and C12 using the relation:

\[B = \frac{C_{11} + 2C_{12}}{3}\]

Parameters:: strain (float, optional) – Magnitude of the strain applied for calculating the underlying elastic constants. Defaults to 0.01.
Returns:: The bulk modulus in Pascals.
Return type:: float

calc_debye_temperature(frame_skip=0.5)

Calculates the Debye temperature.

This method estimates the Debye temperature by first finding the Debye frequency from the (vibrational) Density of States (DOS).

The Debye frequency is defined as the frequency cutoff required for the total number of vibrational modes to equal the system’s total degrees of freedom. It is found by solving the following equation:

\[\int_{0}^{\omega_D} DOS(\omega) d\omega = 3N\]

The calculation performs the following steps:

Calls calc_density_of_states() to get the dos and angular frequency \(\omega\) arrays.
Calculates the cumulative integral of the DOS with respect to \(\omega\) (i.e., the total number of modes up to a given frequency).
Finds the index where this cumulative integral first matches or exceeds the target degrees of freedom (3N).
The frequency at this index is taken as the Debye frequency.
Converts the Debye frequency to the Debye temperature using the relation

\[\Theta_D = \frac{\hbar \omega_D}{k_B}.\]

Parameters:: frame_skip (float, optional) – Fraction of the total initial frames to skip for equilibration. This value is passed directly to calc_density_of_states. Defaults to 0.5.
Returns:: The calculated Debye temperature (\(\Theta_D\)) in Kelvin.
Return type:: float

calc_density_of_states(frame_skip=0.5)

Calculates the (vibrational) Density of States (DOS).

This method computes the DOS by taking the Fourier transform of the Velocity Autocorrelation Function (VACF), a relationship described by the Wiener-Khinchin theorem.

The resulting DOS is cached in self.dos. If this method is called again, it will return the cached values without recalculation.

Parameters:

frame_skip (float, optional) – Fraction of the total initial frames to skip for equilibration. This value is passed directly to _calc_vacf(). Defaults to 0.5.

Returns:

A tuple of (dos, omega):

dos (numpy.ndarray): The 1D array of the normalized density of states.
omega (numpy.ndarray): The 1D array of the corresponding angular frequencies (in rad/s).

Return type:

tuple

calc_isobaric_enthalpy()

Calculates isobaric enthalpy after a NPT ensemble.

The enthalpy (H) is calculated as:

\[H = E + PV\]

where E is the total energy (potential + kinetic), P is the pressure, and V is the volume.

Returns:: An array of enthalpy values for each frame in Joules.
Return type:: numpy.ndarray

Note

Developers, you need a calc to get the total energy. Btw, if you just began reading theese docs, calc stands for calculator.

calc_isobaric_specific_heat()

Calculate the isobaric specific heat (Cp) per unit mass.

This method is intended for results from an NPT ensemble. It calculates Cp from the fluctuations in enthalpy (H) using the formula:

\[C_p = \frac{\langle (H - \langle H \rangle)^2 \rangle}{k_B T^2}\]

Returns:: The isobaric specific heat in units of J / (kg * K).
Return type:: float

calc_isochoric_heat_capacity_per_atom()

Calculate the isochoric heat capacity (Cv) per atom.

This method is intended for results from an NVT ensemble. It calculates Cv from the fluctuations in total energy (E) using the formula:

\[C_v = \frac{\langle (E - \langle E \rangle)^2 \rangle}{k_B T^2}\]

Returns:: The isochoric heat capacity per atom in units of J/K per atom.
Return type:: float

calc_lattice()

Determine the equilibrium lattice constants for the conventional cell.

This method calls _estimate_lattice() to get energy-vs-lattice data and then fits this data to find the lattice parameters that minimize the potential energy. The fitting method depends on the crystal structure:

For cubic or trigonal structures, it performs an Equation of State (EOS) fit.
For hexagonal or tetragonal structures, it performs a 2D quadratic fit to the energy surface.
For orthorhombic, monoclinic, or triclinic structures, it returns the constants found by the sequential optimization in _estimate_lattice().

Returns:

A tuple containing:

cov_structure (str): The name of the crystal structure type.
list: The optimal lattice constants [a, b, c] in Angstroms.

Return type:

tuple

calc_lindemann(a=None)

Compute the global Lindemann parameter.

The Lindemann index is a measure of structural disorder, defined as the root-mean-square displacement divided by the nearest-neighbor distance.

Parameters:: a (float, optional) – Average nearest-neighbor distance (Å). If not provided, it will be calculated using estimate_average_a(). Defaults to None.
Returns:: The dimensionless Lindemann parameter delta_L.
Return type:: float

calc_msd(frame_skip=0.2)

Compute the overall mean squared displacement (MSD).

This method calculates the total MSD, averaged over all three spatial directions.

Parameters:: frame_skip (float, optional) – Fraction of initial frames to skip. Passed to _calc_msd_list(). Defaults to 0.2.
Returns:: The average MSD value across all directions (Å²).
Return type:: float

calc_self_diff()

Calculate the self-diffusion coefficient from the MSD.

The diffusion coefficient (D) is determined from the slope of the mean squared displacement (MSD) vs. time lag (\(\tau\)) plot, based on the Einstein relation in 3D:

\[MSD(\tau) = 6D\tau\]

A linear regression is performed on a central portion of the MSD data to find the slope.

Returns:: The total self-diffusion coefficient in units of Å² per unit of time from the simulation (e.g., Å²/fs).
Return type:: float

calc_shear_modulus(strain=0.01)

Calculate the shear modulus (G) for a cubic crystal.

This method computes the Voigt, Reuss, and Hill averages for the shear modulus based on the calculated elastic constants C11, C12, and C44. It returns the Hill average, which is the arithmetic mean of the Voigt (upper bound) and Reuss (lower bound) moduli.

Parameters:: strain (float, optional) – Magnitude of the strain applied for calculating the underlying elastic constants. Defaults to 0.01.
Returns:: The Hill average shear modulus in Pascals.
Return type:: float

calc_youngs_modulus(strain=0.01)

Calculate Young’s modulus (E).

Young’s modulus is calculated from the bulk (B) and shear (G) moduli using the relation:

\[E = \frac{9BG}{3B + G}\]

Parameters:: strain (float, optional) – Magnitude of the strain applied for calculating the underlying bulk and shear moduli. Defaults to 0.01.
Returns:: Young’s modulus in Pascals.
Return type:: float

check_equilibrium()

Check if the simulation has reached equilibrium.

This method attempts to identify the frame at which the system equilibrates by checking for three conditions in order:

Total energy becomes constant.
Temperature becomes constant.
Total energy begins to oscillate around a stable mean.

It sets self.reached_equilibrium to True if any condition is met.

Returns:: The index of the first frame considered to be in equilibrium. Returns 0 if no equilibrium is detected.
Return type:: int

estimate_average_a()

Estimate the average nearest-neighbor distance over all frames.

Returns:: The average nearest-neighbor distance across all frames in the trajectory (Å).
Return type:: float

estimate_nearest_neighbor_distance(positions)

Estimate average nearest-neighbor distance for one frame.

Parameters:: positions (numpy.ndarray) – A (N, 3) array of atomic positions for N atoms.
Returns:: The average nearest-neighbor distance for the given frame (Å).
Return type:: float

classmethod from_file(filepath)

Create a ResultMD object from a trajectory file.

Parameters:: filepath (str) – Path to the trajectory file (e.g., a .traj file).
Returns:: An instance of the class containing the trajectory frames.
Return type:: ResultMD

get_cohesive_energy()

Calculate the cohesive energy per atom.

Cohesive energy is the energy required to separate a material into isolated, neutral atoms. It is calculated as:

\[E_{\text{cohesive}} = E_{\text{isolated}} - E_{\text{bulk}}\]

where \(E_{\text{isolated}}\) is the average energy of an isolated atom and \(E_{\text{bulk}}\) is the average potential energy per atom in the bulk material.

Returns:: The cohesive energy per atom (eV).
Return type:: float

get_kin_energies()

Get the kinetic energy for each frame.

Returns:: A list of kinetic energy values (eV) for each frame.
Return type:: list

get_pot_energies()

Get the potential energy for each frame.

If potential energy was not saved during the simulation, this method will calculate it for all frames.

Returns:: A list of potential energy values (eV) for each frame.
Return type:: list

get_pressure(frame)

Get pressure for a single frame.

Pressure is calculated as the negative average of the diagonal elements of the stress tensor:

\[P = - \frac{\sigma_{xx} + \sigma_{yy} + \sigma_{zz}}{3}\]

Parameters:: frame (ase.Atoms) – The Atoms object for which to calculate the pressure.
Returns:: The pressure of the frame.
Return type:: float

get_pressures()

Get the pressure for each frame.

Returns:: A list of pressure values for each frame, calculated from the stress tensor.
Return type:: list

get_temperatures()

Get the temperature for each frame.

Returns:: A list of temperature values (K) for each frame.
Return type:: list

get_time_axis()

Generate a time axis for the simulation frames.

Returns:: An array of time values corresponding to each frame, in the same time units as the simulation timestep (e.g., fs).
Return type:: numpy.ndarray

get_tot_energies()

Get the total energy (potential + kinetic) for each frame.

Returns:: A list of total energy values (eV) for each frame.
Return type:: list

plot_density_of_states()

Visualize the density of states (DOS).

This method plots the cached DOS as a function of angular frequency. calc_density_of_states() must be called first to compute and cache the DOS data.

single_atom_energy()

Calculate the average potential energy of isolated atoms.

This method calculates the total potential energy of the individual, non-interacting atoms that make up the material’s formula unit. The result is then averaged per atom in the formula unit. This value is a key component for calculating the cohesive energy.

Returns:: The average potential energy per isolated atom (eV).
Return type:: float

visualize_msd()

Visualize the mean squared displacement (MSD) vs. time lag.

This method generates a plot showing the MSD for the x, y, and z directions as a function of the time lag, tau. It is a convenience wrapper around _calc_msd_list.