PhononDensityOfStates

class PhononDensityOfStates(dynamical_matrix, configuration=None, qpoints=None, number_of_bands=None)

Constructor for the phonon density of states object.

Parameters:

• dynamical_matrix (DynamicalMatrix) – The dynamical matrix to calculate the phonon density of states from.

• configuration (BulkConfiguration | MoleculeConfiguration) – The configuration for which the vibrational modes should be calculated. If no calculator is attached to the configuration, the calculator from dynamical_matrix is used. Currently, configuration must be the same as the configuration from dynamical_matrix and is therefore not needed.
  Default: The configuration from dynamical_matrix

• qpoints (sequence (size 3) of int | MonkhorstPackGrid | KpointDensity | RegularKpointGrid) – The q-points for which the phonon density of states should be calculated.
  Default: The Monkhorst-Pack grid used for the self-consistent calculation.

• number_of_bands (int) – The number of bands.
  Default: All bands

elements()

Query method for the elements in the configuration used for calculating the density of states.

Returns:

A list of the elements.

Return type:

list of PeriodicTableElement

energies()

Query method for the energies used for generating the last spectrum.

Returns:

The energies used for generating the last spectrum.

Return type:

PhysicalQuantity of type energy.

energyMax()

Query method for the maximum energy eigenvalue.

Returns:

The maximum energy eigenvalue.

Return type:

PhysicalQuantity of type energy

energyMin()

Query method for the minimum energy eigenvalue.

Returns:

The minimum energy eigenvalue.

Return type:

PhysicalQuantity of type energy

entropy(temperature=PhysicalQuantity(300.0, K))

Calculate the entropy of the phonons.

The entropy is evaluated using the quasi-harmonic approximation (see Eq.(5.5) in Dove, Martin T. (1993), Introduction to lattice dynamics, Cambridge university press).

Parameters:

temperature (PhysicalQuantity of type temperature.) – The temperature to evaluate the entropy at.
Default: 300 * Kelvin

Returns:

The entropy.

Return type:

PhysicalQuantity with the unit meV / Kelvin

evaluate()

Return the default phonon density of states spectrum.

Returns:

The phonon density of states.

Return type:

PhysicalQuantity with the unit 1 / meV

freeEnergy(temperature=PhysicalQuantity(300.0, K))

Calculate the Helmholtz free energy of the molecule / lattice.

Parameters:

temperature (PhysicalQuantity of type temperature.) – The temperature to evaluate the entropy at.
Default: 300 * Kelvin

Returns:

The free energy of the molecule / lattice (without internal energy).

Return type:

PhysicalQuantity with the unit eV

gaussianSpectrum(energies=None, broadening=None)

Return the phonon density of states spectrum using Gaussian broadening.

Parameters:

• energies (PhysicalQuantity of type energy) – List of energies for which the spectrum should be calculated.
  Default: Energies from emin to emax

• broadening (PhysicalQuantity of type energy) – Broadening of the Gaussian.
  Default: 0.01 * eV

Returns:

The density of states spectrum.

Return type:

PhysicalQuantity with the unit 1 / meV

internalEnergy(temperature=PhysicalQuantity(300.0, K))

Calculate the internal harmonic energy of the molecule / lattice at the given temperature.

Parameters:

temperature (PhysicalQuantity of type temperature.) – The temperature to evaluate the internal energy at.
Default: 300 * Kelvin

Returns:

The internal energy of the molecule / lattice.

Return type:

PhysicalQuantity with the unit eV

metatext()

Returns:

The metatext of the object or None if no metatext is present.

Return type:

str | None

nlprint(stream=None)

Print a string containing an ASCII table useful for plotting the AnalysisSpin object.

Parameters:

stream (python stream) – The stream the table should be written to.
Default: NLPrintLogger()

setMetatext(metatext)

Set a given metatext string on the object.

Parameters:

metatext (str | None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

tetrahedronSpectrum(energies=None)

Return the phonon density of states spectrum using the tetrahedron method.

Parameters:

energies (PhysicalQuantity of type energy) – List of energies for which the spectrum should be calculated.
Default: Energies from emin to emax

Returns:

The density of states spectrum.

Return type:

PhysicalQuantity with the unit 1 / meV

uniqueString()

Return a unique string representing the state of the object.

zeroPointEnergy()

Calculate the zero point energy in the quasi-harmonic approximation.

Returns:

The zero point energy of the molecule / lattice

Return type:

PhysicalQuantity with the unit eV

Usage Examples

Calculate the phonon density of states of a silicon crystal using the TremoloXCalculator calculator with the Stillinger-Weber potential [1]. The example also shows how to calculate the vibrational entropy of the crystal unit cell.

# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

# Set up lattice
lattice = FaceCenteredCubic(5.4306*Angstrom)

# Define elements
elements = [Silicon, Silicon]

# Define coordinates
fractional_coordinates = [[0.0 , 0.0 , 0.0 ],
                          [0.25, 0.25, 0.25]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
potentialSet = StillingerWeber_Si_1985()
calculator = TremoloXCalculator(parameters=potentialSet)

bulk_configuration.setCalculator(calculator)

# -------------------------------------------------------------
# Dynamical matrix
# -------------------------------------------------------------
dynamical_matrix = DynamicalMatrix(
    configuration=bulk_configuration,
    repetitions=Automatic,
    atomic_displacement=0.01*Angstrom,
    acoustic_sum_rule=False,
    symmetrize=True,
    finite_difference_method=Central,
    )

# -------------------------------------------------------------
# Phonon density of states
# -------------------------------------------------------------
phonon_density_of_states = PhononDensityOfStates(
    configuration=bulk_configuration,
    dynamical_matrix=dynamical_matrix,
    qpoints=MonkhorstPackGrid(10,10,10),
    )

nlsave('si_phonon_dos.nc', phonon_density_of_states)

# Print the results ot the log
nlprint(phonon_density_of_states)

# Calculate the vibrational entropy at 300 K
entropy = phonon_density_of_states.entropy(300.0*Kelvin)

si_phonon_dos.py

Note

The spectrum can be shown by selecting the file si_phonon_dos.nc in the Lab Floor of QuantumATK and opening it using the 2D Plot tool.

Plot the Gaussian and tetrahedron spectrum using pylab:

# read phonon DOS object from file
dos = nlread('si_phonon_dos.nc', PhononDensityOfStates)[0]

# make list of energies
energies = numpy.arange(0.0, 80.0, 0.2)*meV

# calculate the DOS spectrum with two different methods
dos_t = dos.tetrahedronSpectrum(energies)
dos_g = dos.gaussianSpectrum(energies, broadening = 1.0*meV)

#plot the spectra using pylab
import pylab
pylab.figure()
pylab.plot(energies.inUnitsOf(meV), dos_t.inUnitsOf(meV**-1))
pylab.plot(energies.inUnitsOf(meV), dos_g.inUnitsOf(meV**-1))
pylab.xlabel("Energy (meV)")
pylab.ylabel("Phonon DOS (1/meV)")
pylab.show()

si_phonon_dos_plot.py

Notes

• The details of the force constant calculation to obtain the dynamical matrix can be found under DynamicalMatrix.

• The routine utilizes the symmetries of the crystal to reduce the computational load.

• The default spectrum is a Gaussian spectrum for a MoleculeConfiguration and a BulkConfiguration with less than 10 k-points. For a BulkConfiguration with more than 10 k-points, the tetrahedron method is used.

• The implementation of the tetrahedron method follows Ref. [2]

The PhononDensityOfStates class can also be used to calculate thermodynamic properties in the harmonic approximation. The following properties can be calculated:

• The free energy at a given temperature can be calculated via the freeEnergy() method:

\[F = -\frac{1}{N_q} \sum_{q,s} \frac{E(q, s)}{2} + \frac{1}{\beta N_q} \sum_{q,s} \ln (1-\exp(-\beta E(q, s)))\]

where \(\beta=1/k_B T\) and \(N_q\) denotes the number of q-points to sum over.

• The internal energy at a given temperature, including zero-point energy, can be calculated via the internalEnergy() method:

\[U = -\frac{1}{N_q} \sum_{q,s} \frac{E(q, s)}{2} + \frac{1}{N_q} \sum_{q,s} \frac{ E(s,q) }{ \exp(\beta E(s,q)) - 1 }\]

In this formula, the first part describes the zero-point energy and the second part stands for the finite temperature contribution of the phonons to the internal energy.

• The entropy at a given temperature can be calculated via the entropy() method:

\[S = \frac{1}{N_q} \sum_{q,s} \frac{E(q,s)}{T} \frac{\exp(-\beta E(q,s))}{1 - \exp(-\beta E(q,s))} - \frac{k_B}{N_q} \sum_{q,s} \ln(1 - \exp(-\beta E(q,s)))\]

Note that all these quantities do not include any potential energy contributions, as calculated via TotalEnergy, so in order to get the absolute free energy or internal energy values, one would have to add these contributions on top.

[1]

F. H. Stillinger and T. A. Weber. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B, 31:5262–5271, Apr 1985. doi:10.1103/PhysRevB.31.5262.

[2]

P. E. Blöchl, O. Jepsen, and O. K. Andersen. Improved tetrahedron method for brillouin-zone integrations. Phys. Rev. B, 49:16223–16233, Jun 1994. doi:10.1103/PhysRevB.49.16223.