10.3. Dynamical matrix and Born effective charges

In this section we consider case of atomic displacements as perturbations. For ground state calculation, it is assumed that atoms are fixed at their equilibrium position. However, in reality they often vibrate around their equilibrium positions, which is the origin of important physical phenomena related to lattice dynamics. The dynamical matrix characterizes the vibration of atoms, known as phonons, and is related to “2nd” order derivatives of total energy with respect to atomic displacements. DFPT provides an efficient framework to calculate the dynamical matrix and related physical properties [BGCG01].

Phonons, like any other wavelike quantities, are characterized by a wave-vector \(\mathbf{q}\). Some physical quantities depend on phonons with general wave-vectors, such as phonon band-structure and density of states, and others are only coupled to phonons at \(\Gamma\), i.e. \(\mathbf{q}=0\) point, such as optical properties at infrared region. This determines reciprocal-space sampling of phonons in a DFPT calculation.

Note

Number of independent calculations to obtain all possible phonon modes for \(N_{q}\) wave-vectors in a system with \(N_{atom}\) is \(3\times N_{atom}\times N_{q}\). Available symmetries can be used to reduce the total number of calculations.

Born effective charges (BECs) quantify the coupling between the optical phonons in long wavelength limit and electric fields, in other words they are the coefficient of proportionality between the polarization response in one direction caused by an atomic displacement in another direction [GL97]. In RESCU BECs are obtained during phonon calculations at \(\Gamma\) wave-vector.

In the following we perform a phonon DFPT calculation for boron nitride example by using the ground state data obtained in Preparing the ground state data. A typical input file for phonon calculation is given below:

info.calculationType = 'dfpt-phonon'
info.savepath = './results/bn_real_phonon'
rho.in = './results/bn_real_scf.mat'
psi = './results/bn_real_scf.h5'
dfpt.qpointGridn = [3,3,3]
mixing.tol = 1e-5*[1 1]
symmetry.spacesymmetry = 1    % MUST be set to false if perturbed wavefunctions and potentials are
                              % supposed to be used in subsequent nonself-consistent calculation
symmetry.pointsymmetry = 1    % MUST be set to false if perturbed wavefunctions and potentials are
                              % supposed to be used in subsequent nonself-consistent calculation
symmetry.timereversal = 1     % use to modify the default (true) value
option.saveWavefunction = 1   % only required if perturbed wavefunctions are supposed
                              % to be used for subsequent nonself-consistent calculation

The keywords are explained as follows:

  • info.calculationType determines the type of the calculation;

  • info.savepath points to the file which stores the output data of the phonon calculation;

  • rho.in this DFPT calculation loads the ground state electronic density at the beginning of the computation, and this keyword points to the file which contains the output of the ground state calculation, so that RESCU can load corresponding electronic density data;

  • psi this DFPT calculation also loads the ground state wavefunctions at the beginning of the computation, and this keyword points to the file which contains the stored ground state wavefunctions;

  • dfpt.qpointGridn determines the sampling of the reciprocal space for phonon wave-vectors. If the data is to be used to obtain a phonon BS or DOS, it must be identical to kpoint.gridn of the ground state data, or can be simply removed from the input file. To perform the calculation at only \(\Gamma\)-point it should be set to [1,1,1];

Important

Current DFPT implementation in RESCU requires a \(\Gamma\)-centered and odd reciprocal-space sampling for both ground state and phonons. Furthermore, kpoint.gridn and dfpt.qpointGridn must be identical for case of phonon with general q-sampling. If one wishes to only calculate phonons at \(\Gamma\) wave-vector, dfpt.qpointGridn can be set to [1,1,1].

  • mixing.tol first entry sets the convergence criterion for the perturbed self-consistent density and the second entry sets the criterion for the perturbed self-consistent potential;

  • option.saveWavefunction determines whether first-order perturbed wavefunctions should be saved to disk. This is only necessary if one wishes to perform a subsequent DFPT calculation based “2n+1” theorem to compute non-linear response functions.

Warning

Storing the perturbed wavefunctions requires a large amount of disk space, and one must make sure that they are only saved if are to be used later

  • symmetry.spacesymmetry determine whether or not crystal space-group symmetry should be used in calculating the dynamical matrix. The default value is true;

  • symmetry.pointsymmetry determine whether or not crystal point-group symmetry should be used in calculating the dynamical matrix. The default value is true;

  • symmetry.timereversal determine whether or not time-reversal symmetry should be used in calculating the dynamical matrix. The default value is true.

Important

Currently symmetry can be used in RESCU to reduce the total number of required calculations to compute the dynamical matrix. However, if perturbed wavefunctions and potentials are also required for subsequent nonself-consistent DFPT calculations (e.g. see The Raman spectrum and non-linear optical susceptibility and Electron-phonon interaction), then only time-reversal symmetry can be applied to phonon calculations. In such cases, the following setup for symmetry keywords must be used

symmetry.timereversal   = 1
symmetry.spacesymmetry  = 0
symmetry.pointsymmetry  = 0

Once the input file is ready, RESCU can be executed by

rescu -i bn_real_phonon

Note

Examining the output log shows that symmetry has drastically reduces the number of independent atomic displacements perturbations from 162 to 16.

Symmetry status for dynamical matrix calculation
(time-reversal,space-group,qpoint-group)     = (1,1,1)
Number of qpoints (total)                    = 27
Number of qpoints (to consider)              = 4
Number of atomic displacements (total)       = 162
Number of atomic displacements (to consider) = 16

Upon completion, the output data can be found in dfpt structure in output file ./results/bn_real_phonon.mat.

>> load results/bn_real_phonon.mat
>> dfpt

dfpt =

struct with fields:

         BornEC: [3x3x2 double]
            desc: [1x1 struct]
         dynMat: [6x6x27 double]
   qpointDirect: [27x3 double]
   qpointGridn: [3 3 3]
            type: 'phonon'

Two main output properties are dynamical matrix which is stored as a \(3N_{atoms}\times 3N_{atoms}\times N_{q}\) array in dfpt.dynMat and Born effective charges (BECs) which is stored as \(3\times 3\times N_{atoms}\) in dfpt.BornEC

>> dfpt.dynMat(:,:,1)

ans =

Columns 1 through 4

   0.2409 + 0.0000i   0.0000 + 0.0000i   0.0000 - 0.0000i  -0.2403 - 0.0000i
   0.0000 - 0.0000i   0.2409 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 - 0.0000i   0.0000 - 0.0000i   0.2409 + 0.0000i   0.0000 + 0.0000i
  -0.2403 - 0.0000i  -0.0000 - 0.0000i  -0.0000 - 0.0000i   0.2419 - 0.0000i
  -0.0000 - 0.0000i  -0.2403 - 0.0000i  -0.0000 - 0.0000i   0.0002 + 0.0000i
  -0.0000 - 0.0000i  -0.0000 + 0.0000i  -0.2403 - 0.0000i   0.0002 + 0.0000i

Columns 5 through 6

   0.0000 + 0.0000i   0.0000 + 0.0000i
  -0.2403 + 0.0000i   0.0000 - 0.0000i
   0.0000 + 0.0000i  -0.2403 + 0.0000i
   0.0002 - 0.0000i   0.0002 + 0.0000i
   0.2419 - 0.0000i   0.0002 - 0.0000i
   0.0002 + 0.0000i   0.2419 - 0.0000i

>> dfpt.BornEC

ans(:,:,1) =

    1.7356    0.0000    0.0000
   -0.0000    1.7356   -0.0000
   -0.0000   -0.0000    1.7356


ans(:,:,2) =

   -2.8601   -0.0000   -0.0000
   -0.0000   -2.8601   -0.0000
   -0.0000   -0.0000   -2.8601

Note

Born effective charges are set to zero for metallic calculations, and also for non-metallic systems with only one species type (non-polar materials).

DFPT dielectric and phonon calculations provide three important response functions, i.e., ion-clamped dielectric constant, dynamical matrix, and Born effective charges which are used in variety of nonself-consistent DFPT calculations to extract further information regarding the response of materials to static electric field and atomic displacements perturbations. In particular, one can continue to compute full phonons band-structure and density of states which are covered in Phonon band structure and Phonon density of states. For other available nonself-consistent DFPT features please refer to DFPT.

[BGCG01]

Stefano Baroni, Stefano De Gironcoli, Andrea Dal Corso, and Paolo Giannozzi. Phonons and related crystal properties from density-functional perturbation theory. Reviews of Modern Physics 73.2 (2001), p. 515.