index | modules | gitlab | page source

# Vibrational modes

You can calculate the vibrational modes of an `Atoms` object in the harmonic approximation using the `Vibrations`.

class ase.vibrations.Vibrations(atoms, indices=None, name='vib', delta=0.01, nfree=2)

  Class for calculating vibrational modes using finite difference.

  The vibrational modes are calculated from a finite difference approximation of the Hessian matrix.

  The summary(), get_energies() and get_frequencies() methods all take an optional method keyword. Use method=’Frederiksen’ to use the method described in:

  > T. Frederiksen, M. Paulsson, M. Brandbyge, A. P. Jauho: “Inelastic transport theory from first-principles: methodology and applications for nanoscale devices”, Phys. Rev. B 75, 205413 (2007)

  atoms: Atoms object

    The atoms to work on.

  indices: list of int

    List of indices of atoms to vibrate. Default behavior is to vibrate all atoms.

  name: str

    Name to use for files.

  delta: float

    Magnitude of displacements.

  nfree: int

    Number of displacements per atom and cartesian coordinate, 2 and 4 are supported. Default is 2 which will displace each atom +delta and -delta for each cartesian coordinate.

  Example:

      >>> from ase import Atoms
      >>> from ase.calculators.emt import EMT
      >>> from ase.optimize import BFGS
      >>> from ase.vibrations import Vibrations
      >>> n2 = Atoms('N2', [(0, 0, 0), (0, 0, 1.1)],
      ...            calculator=EMT())
      >>> BFGS(n2).run(fmax=0.01)
      BFGS:   0  16:01:21        0.440339       3.2518
      BFGS:   1  16:01:21        0.271928       0.8211
      BFGS:   2  16:01:21        0.263278       0.1994
      BFGS:   3  16:01:21        0.262777       0.0088
      >>> vib = Vibrations(n2)
      >>> vib.run()
      >>> vib.summary()
      ---------------------
      #    meV     cm^-1
      ---------------------
      0    0.0       0.0
      1    0.0       0.0
      2    0.0       0.0
      3    1.4      11.5
      4    1.4      11.5
      5  152.7    1231.3
      ---------------------
      Zero-point energy: 0.078 eV
      >>> vib.write_mode(-1)  # write last mode to trajectory file

  clean(empty_files=False, combined=True)

    Remove json-files.

    Use empty_files=True to remove only empty files and combined=False to not remove the combined file.

  combine()

    Combine json-files to one file ending with ‘.all.json’.

    The other json-files will be removed in order to have only one sort of data structure at a time.

  fold(frequencies, intensities, start=800.0, end=4000.0, npts=None, width=4.0, type='Gaussian', normalize=False)

    Fold frequencies and intensities within the given range and folding method (Gaussian/Lorentzian). The energy unit is cm^-1. normalize=True ensures the integral over the peaks to give the intensity.

  get_energies(method='standard', direction='central', **kw)

    Get vibration energies in eV.

  get_frequencies(method='standard', direction='central')

    Get vibration frequencies in cm^-1.

  get_mode(n)

    Get mode number .

  get_vibrations(method='standard', direction='central', read_cache=True, **kw)

    Get vibrations as VibrationsData object

    If read() has not yet been called, this will be called to assemble data from the outputs of run(). Most of the arguments to this function are options to be passed to read() in this case.

    Parameters:

      - method (str) – Calculation method passed to read()

      - direction (str) – Finite-difference scheme passed to read()

      - read_cache (bool) – The VibrationsData object will be cached for quick access. Set False to force regeneration of the cache with the current atoms/Hessian/indices data.

      - **kw – Any remaining keyword arguments are passed to read()

    Returns:

      VibrationsData

  iterdisplace(inplace=False) → Iterator[Tuple[Displacement, Atoms]]

    Iterate over initial and displaced structures.

    Use this to export the structures for each single-point calculation to an external program instead of using `run()`. Then save the calculated gradients to <name>.json and continue using this instance.

    ## Parameters:

    inplace: bool

      If True, the atoms object will be modified in-place. Otherwise a copy will be made.

    ## Yields:

    disp: Displacement

      Displacement object with information about the displacement.

    atoms: Atoms

      Atoms object with the displaced structure.

  iterimages()

    Yield initial and displaced structures.

  run()

    Run the vibration calculations.

    This will calculate the forces for 6 displacements per atom +/-x, +/-y, +/-z. Only those calculations that are not already done will be started. Be aware that an interrupted calculation may produce an empty file (ending with .json), which must be deleted before restarting the job. Otherwise the forces will not be calculated for that displacement.

    Note that the calculations for the different displacements can be done simultaneously by several independent processes. This feature relies on the existence of files and the subsequent creation of the file in case it is not found.

    If the program you want to use does not have a calculator in ASE, use `iterdisplace` to get all displaced structures and calculate the forces on your own.

  split()

    Split combined json-file.

    The combined json-file will be removed in order to have only one sort of data structure at a time.

  write_dos(out='vib-dos.dat', start=800, end=4000, npts=None, width=10, type='Gaussian', method='standard', direction='central')

    Write out the vibrational density of states to file.

    First column is the wavenumber in cm^-1, the second column the folded vibrational density of states. Start and end points, and width of the Gaussian/Lorentzian should be given in cm^-1.

  write_jmol()

    Writes file for viewing of the modes with jmol.

  write_mode(n=None, kT=0.02585199101165164, nimages=30)

    Write mode number n to trajectory file. If n is not specified, writes all non-zero modes.

## Example

The example of a water molecule in the EAM potential

    from ase.build import molecule
    from ase.calculators.emt import EMT
    from ase.optimize import BFGS
    from ase.vibrations import Vibrations

    h2o = molecule('H2O')
    h2o.calc = EMT()
    BFGS(h2o).run(fmax=0.01)

    vib = Vibrations(h2o)
    vib.run()
    vib.summary(log='H2O_EMT_summary.txt')
    vib.write_mode(-1)

where the output

    ---------------------
      #    meV     cm^-1
    ---------------------
      0    6.3i     51.0i
      1    5.9i     47.6i
      2    0.0i      0.3i
      3    0.0i      0.1i
      4    0.1       0.4
      5    5.4      43.8
      6   32.1     258.9
      7  296.7    2392.7
      8  387.4    3124.9
    ---------------------
    Zero-point energy: 0.361 eV

shows 3 meaningful vibrations (the last 3 with highest energies.)

These vibrations can be viewed in `ase gui` either by writing them out as a “movie”:

    vib.write_mode(-1)

which writes out the file `vib.8.traj` The vibrations can also be encoded as forces:

    vib.show_as_force(8)

which opens `ase gui` automatically and the forces point into directions of the movement of the atoms.

## Old calculations

The output format of vibrational calculations was changed from `pickle` to `json`. There is a tool to convert old `pickle`-files:

    > python3 -m ase.vibrations.pickle2json mydirectory/vib.*.pckl

## Vibrational Data

Vibrational data is stored inside the `VibrationsData` class.

class ase.vibrations.VibrationsData(atoms: Atoms, hessian: Sequence[Sequence[Sequence[Sequence[Real]]]] | ndarray, indices: Sequence[int] | ndarray = None)

  Class for storing and analyzing vibrational data (i.e. Atoms + Hessian)

  This class is not responsible for calculating Hessians; the Hessian should be computed by a Calculator or some other algorithm. Once the VibrationsData has been constructed, this class provides some common processing options; frequency calculation, mode animation, DOS etc.

  If the Atoms object is a periodic supercell, VibrationsData may be converted to a PhononData using the VibrationsData.to_phonondata() method. This provides access to q-point-dependent analyses such as phonon dispersion plotting.

  If the Atoms object has FixedAtoms or FixedCartesian constraints, these will be respected and the Hessian should be sized accordingly.

  Parameters:

    - atoms – Equilibrium geometry of vibrating system. This will be stored as a full copy.

    - hessian – Second-derivative in energy with respect to Cartesian nuclear movements as an (N, 3, N, 3) array.

    - indices – indices of atoms which are included in Hessian. Default value (None) includes all freely moving atoms (i.e. not fixed ones). Leave at None if constraints should be determined automatically from the atoms object.

  classmethod from_2d(atoms: Atoms, hessian_2d: Sequence[Sequence[Real]] | ndarray, indices: Sequence[int] = None) → VibrationsData

    Instantiate VibrationsData when the Hessian is in a 3Nx3N format

    Parameters:

      - atoms – Equilibrium geometry of vibrating system

      - hessian – Second-derivative in energy with respect to Cartesian nuclear movements as a (3N, 3N) array.

      - indices – Indices of (non-frozen) atoms included in Hessian

  get_dos() → RawDOSData

    Total phonon DOS

  get_energies() → ndarray

    Diagonalise the Hessian to obtain eigenvalues

    Results are cached so diagonalization will only be performed once for this object instance.

    Returns:

      Harmonic mode energies in units of eV

  get_energies_and_modes(all_atoms: bool = False) → Tuple[ndarray, ndarray]

    Diagonalise the Hessian to obtain harmonic modes

    Results are cached so diagonalization will only be performed once for this object instance.

    Parameters:

      all_atoms – If True, return modes as (3N, [N + N_frozen], 3) array where the second axis corresponds to the full list of atoms in the attached atoms object. Atoms that were not included in the Hessian will have displacement vectors of (0, 0, 0).

    Returns:

      tuple (energies, modes)

      Energies are given in units of eV. (To convert these to frequencies in cm-1, divide by ase.units.invcm.)

      Modes are given in Cartesian coordinates as a (3N, N, 3) array where indices correspond to the (mode_index, atom, direction).

  get_frequencies() → ndarray

    Diagonalise the Hessian to obtain frequencies in cm^-1

    Results are cached so diagonalization will only be performed once for this object instance.

    Returns:

      Harmonic mode frequencies in units of cm^-1

  get_hessian() → ndarray

    The Hessian; second derivative of energy wrt positions

    This format is preferred for iteration over atoms and when addressing specific elements of the Hessian.

    Returns:

      array with shape (n_atoms, 3, n_atoms, 3) where

      - the first and third indices identify atoms in self.get_atoms()

      - the second and fourth indices cover the corresponding Cartesian movements in x, y, z

      e.g. the element h[0, 2, 1, 0] gives a harmonic force exerted on atoms[1] in the x-direction in response to a movement in the z-direction of atoms[0]

  get_hessian_2d() → ndarray

    Get the Hessian as a 2-D array

    This format may be preferred for use with standard linear algebra functions

    Returns:

      array with shape (n_atoms * 3, n_atoms * 3) where the elements are ordered by atom and Cartesian direction:

          >> [[at1x_at1x, at1x_at1y, at1x_at1z, at1x_at2x, ...],
          >> [at1y_at1x, at1y_at1y, at1y_at1z, at1y_at2x, ...],
          >> [at1z_at1x, at1z_at1y, at1z_at1z, at1z_at2x, ...],
          >> [at2x_at1x, at2x_at1y, at2x_at1z, at2x_at2x, ...],
          >> ...]

      e.g. the element h[2, 3] gives a harmonic force exerted on atoms[1] in the x-direction in response to a movement in the z-direction of atoms[0]

  get_mask() → ndarray

    Boolean mask of atoms selected by indices

  get_modes(all_atoms: bool = False) → ndarray

    Diagonalise the Hessian to obtain harmonic modes

    Results are cached so diagonalization will only be performed once for this object instance.

    all_atoms:

      If True, return modes as (3N, [N + N_frozen], 3) array where the second axis corresponds to the full list of atoms in the attached atoms object. Atoms that were not included in the Hessian will have displacement vectors of (0, 0, 0).

    Returns:

      Modes in Cartesian coordinates as a (3N, N, 3) array where indices correspond to the (mode_index, atom, direction).

  get_pdos() → DOSCollection

    Phonon DOS, including atomic contributions

  get_zero_point_energy() → float

    Diagonalise the Hessian and sum hw/2 to obtain zero-point energy

    Parameters:

      energies – Pre-computed energy eigenvalues. Use if available to avoid re-calculating these from the Hessian.

    Returns:

      zero-point energy in eV

  static indices_from_constraints(atoms: Atoms) → List[int]

    Indices corresponding to Atoms Constraints

    Deduces the freely moving atoms from the constraints set on the atoms object. VibrationsData only supports FixCartesian and FixAtoms. All others are neglected.

    Parameters:

      atoms – Atoms object.

    Retruns:

      indices of free atoms.

  static indices_from_mask(mask: Sequence[bool] | ndarray) → List[int]

    Indices corresponding to boolean mask

    This is provided as a convenience for instantiating VibrationsData with a boolean mask. For example, if the Hessian data includes only the H atoms in a structure:

        h_mask = atoms.get_chemical_symbols() == 'H'
        vib_data = VibrationsData(atoms, hessian,
                                  VibrationsData.indices_from_mask(h_mask))

    Take care to ensure that the length of the mask corresponds to the full number of atoms; this function is only aware of the mask it has been given.

    Parameters:

      mask – a sequence of True, False values

    Returns:

      indices of True elements

  iter_animated_mode(mode_index: int, temperature: float = 0.02585199101165164, frames: int = 30) → Iterator[Atoms]

    Obtain animated mode as a series of Atoms

    Parameters:

      - mode_index – Selection of mode to animate

      - temperature – In energy units - use units.kB * T_IN_KELVIN

      - frames – number of image frames in animation

    Yields:

      Displaced atoms following vibrational mode

  classmethod read(fd)

    Read new instance from JSON file.

  show_as_force(mode: int, scale: float = 0.2, show: bool = True) → Atoms

    Illustrate mode as “forces” on atoms

    Parameters:

      - mode – mode index

      - scale – scale factor

      - show – if True, open the ASE GUI and show atoms

    Returns:

      Atoms with scaled forces corresponding to mode eigenvectors (using attached SinglePointCalculator).

  tabulate(im_tol: float = 1e-08) → str

    Get a summary of the vibrational frequencies.

    Parameters:

      im_tol – Tolerance for imaginary frequency in eV. If frequency has a larger imaginary component than im_tol, the imaginary component is shown in the summary table.

    Returns:

      Summary table as formatted text

  with_new_masses(masses: Sequence[float] | ndarray) → VD

    Get a copy of vibrations with modified masses and the same Hessian

    Parameters:

      masses – New sequence of masses corresponding to the atom order in self.get_atoms()

    Returns:

      A copy of the data with new masses for the same Hessian

  write(fd)

    Write to JSON file.

  write_jmol(filename: str = 'vib.xyz', ir_intensities: Sequence[float] | ndarray = None) → None

    Writes file for viewing of the modes with jmol.

    This is an extended XYZ file with eigenvectors given as extra columns and metadata given in the label/comment line for each image. The format is not quite human-friendly, but has the advantage that it can be imported back into ASE with ase.io.read.

    Parameters:

      - filename – Path for output file

      - ir_intensities – If available, IR intensities can be included in the header lines. This does not affect the visualisation, but may be convenient when comparing to experimental data.