Phonons and vibration modes

Phonons from finite differences

The finite difference method can be used to evaluate the response of an electronic system to perturbation without relying explicitly on density functional perturbation theory. It is for instance possible to use finite differences to compute the second-order force constants of a system which is defined as

$$ \Psi_{I\alpha J\beta} = \frac{\partial^2 E_{\rm tot}}{\partial u_{I\alpha}u_{J\beta}} = -\frac{\partial F_{I\alpha}}{\partial u_{J\beta}} .\,\! $$ where $I,J$ are atoms of the systems, $\alpha,\beta=x,y,z$, $E_{\rm tot}$ is the total energy of the system, and $\mathbf{F}_I[\psi]$ is the force acting on ion $I$ given a set of Kohn-Sham orbitals $\{\Psi_{n,\mathbf{k}}\}$.

In the finite differences approach, this is evaluate using $$ \Psi_{I\alpha J\beta} \approx - \frac{(F_{I\alpha}[\{\Psi_{u_{J\beta}\lambda}\}]-F_{I\alpha}[\{\Psi_{u_{J\beta}-\lambda}\}])}{2\lambda}\,\! $$

This is what is implemented in Octopus when using the following input option:


  CalculationMode = vib_modes
  ResponseMethod = finite_differences 

In particular, in the finite difference approach, ions get displaced with a magnitude of $\lambda$, which is controlled in the code by the variable Displacement. The displacements are refered in the above formula as $u_{I\alpha}$.

Phonons from linear response

Alternatively to the finite differences method, it is possible to compute the second-order force constants based on perturbation theory. For this, we use the fact that

$$ \frac{dE_{\rm tot}}{d\lambda^2} = \Big\langle \frac{d\psi}{d\lambda} \Big| \frac{dV}{d\lambda} \Big|\psi \big\rangle + \Big\langle \psi \Big| \frac{dV}{d\lambda} \Big|\frac{d\psi}{d\lambda} \big\rangle - \Big\langle \psi \Big| \frac{d^2V}{d^2\lambda} \Big|\psi \big\rangle .\,\! $$

This expression requires to know the first order change of the wavefunctions with respect to an ionic displacement. In Octopus, this is obtained solving the Sternheimer equation.

In order to use the Sterheimer approach, one can use the following options:


  CalculationMode = vib_modes
  ResponseMethod = sternheimer

  RestartFixedOccupations = no

The line containing ResponseMethod is not mandatory here, as the “sternheimer” option is the default value. The Sternheimer version offers several related options: CalcNormalModeWfs, CalcInfrared, and SymmetrizeDynamicalMatrix, see the documetation of these variables for more details. Note that the last line is required due to internal constrains.