# 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 (MonkhorstPackGrid | RegularKpointGrid) – The q-points for which the phonon density of states should be calculated. Default: MonkhorstPackGrid(nx, ny, nz), where nx, ny, nz is the sampling 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. list of PeriodicTableElement
energies()

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

Returns: The energies used for generating the last spectrum. PhysicalQuantity of type energy.
energyMax()

Query method for the maximum energy eigenvalue.

Returns: The maximum energy eigenvalue. PhysicalQuantity of type energy
energyMin()

Query method for the minimum energy eigenvalue.

Returns: The minimum energy eigenvalue. 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 The entropy. PhysicalQuantity with the unit meV / Kelvin
evaluate()

Return the default phonon density of states spectrum.

Returns: The phonon density of states. 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 The free energy of the molecule / lattice (without internal energy). 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 The density of states spectrum. PhysicalQuantity with the unit 1 / meV
metatext()
Returns: The metatext of the object or None if no metatext is present. str | unicode | 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 | unicode | 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 The density of states spectrum. PhysicalQuantity with the unit 1 / meV
zeroPointEnergy()

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

Returns: The zero point energy of the molecule / lattice 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 [kSW85]. 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

# 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. [kBJA94]

 [kBJA94] 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.
 [kSW85] 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.