class VibrationalDensityOfStates(md, start_time=None, end_time=None, atom_selection=None, time_resolution=None, info_panel=None)

Class for calculating the vibrational density of states of an MD simulation.

  • md_trajectory (MDTrajectory) – The MDtrajectory to calculate the vibrational density of states for.
  • start_time (PhysicalQuantity of type time) – The start time.
    Default: 0.0 * fs
  • end_time (PhysicalQuantity of type time) – The end time.
    Default: The last time frame
  • atom_selection (PeriodicTableElement | str | list of ints) – Only include contributions from this selection. The atoms can be selected by element i.e. PeriodicTableElement, tag or a list of atomic indices.
    Default: All atoms.
  • time_resolution (PhysicalQuantity of type time) – The time interval between snapshots in the MD trajectory that are included in the analysis.
  • info_panel (InfoPanel (Plot2D)) – Info panel to show the calculation progress.
    Default: No info panel
Returns:The frequencies associated with the vibrational density of states.
Return type:PhysicalQuantity of type 1/time
Returns:The vibrational density-of-states, calculated from the Fourier transformation of the velocities.
Return type:PhysicalQuantity of type time

Usage Examples

Load an MDTrajectory and calculate the vibrational density of states from the simulation:

md_trajectory = nlread('alumina_trajectory.hdf5', class_type=MDTrajectory)[-1]

# Initialize a velocity autocorrelation object.
vacf = VelocityAutocorrelation(md_trajectory, atom_selection=Oxygen)

# Get the vibrational density of states.
vdos = vacf.vibrationalDensityOfStates()

# Get the frequencies associated with the DOS.
frequencies = vacf.frequencies()

# Convert the frequencies to energies in meV.
energies = (frequencies*planck_constant).inUnitsOf(meV)

# Plot the data using pylab.
import pylab

pylab.plot(energies, vdos, label='Vibrational density of states')
pylab.xlabel('E (meV)')


The vibrational density of states is calculated from the Fourier transform of the atomic velocities (cf. [LBGI03]):

\[S(\nu) = \frac{2}{k_B T} \sum_{i=1}^{N} \sum_{\alpha=x,y,z} m_i s_i^{\alpha}(\nu),\]


\[s_i^{\alpha}(\nu) = \frac{ \left \lvert \int_0^{t_{max}} dt' v_i^{\alpha}(t')e^{-i2\pi\nu t} \right \rvert^2 } { t_{max} } \, .\]

The temperature entering the first expression is calculated directly from the velocities during the simulation.

The vibrational DOS are normalized in such a way that the integration over the frequencies (if given in \(fs^{-1}\) ) gives the total number of degrees of freedom in the system

\[\int_0^{\infty}d\nu S(\nu) = 3N \, .\]

Apart from the normalization, the vibrational density of states should give similar results as the PhononDensityOfStates with the important distinction, that in the former, anharmonic effects are seamlessly included. You also have the possibility to extract the contribution of only a part of the system, specified via the atom_selection parameter. Moreover, in contrast to the PhononDensityOfStates analysis object, the vibrational density of states can be calculated for soft systems or even for liquid systems.

In order to sample all relevant vibrations, the system size must be chosen large enough to accommodate vibrational modes with long wavelengths. The log interval for saving snapshots to the MDTrajectory must be chosen such that the sampling rate (the reciprocal of the time step times the logging interval) is twice the highest vibrational frequency that should be included in the density of states. For example, if the time step is 1 femtosecond and the logging interval is 10 steps, then the sampling rate is 100 THz and the maximum frequency in the density of states is 50 THz.

[LBGI03]S.-T. Lin, M. Blanco, and W. A Goddard III. The two-phase model for calculating thermodynamic properties of liquids from molecular dynamics: Validation for the phase diagram of Lennard-Jones fluids. J. Chem. Phys., 2003. doi:10.1063/1.1624057.