Octopus & phonopy calculation

Octopus & phonopy calculation#

Octopus is a real-space (grid based) TDDFT / DFT code using pseudopotentials. This page explains how to calculate phonons with phonopy using Octopus as the force calculator, i.e., using the finite displacement and supercell approach.

How the interface works#

The Octopus interface differs from the other calculator interfaces in the way the crystal structure is exchanged with phonopy:

  • Octopus users normally describe the crystal structure in an Octopus input file (inp), using %LatticeParameters, %LatticeVectors and %ReducedCoordinates blocks. phonopy, however, reads the unit cell from a VASP-style POSCAR file (the default cell file name for --octopus). Octopus can write exactly this file, so the workflow starts by exporting the unit cell from Octopus as a POSCAR (see below). Lattice vectors in POSCAR are given in Angstrom, as usual; phonopy converts them internally to atomic units. As an alternative, phonopy also accepts the unit cell directly as an Octopus geometry include file (auto-detected, lattice in bohr), for users who prefer to keep the structure in Octopus form.

  • The supercells are written by phonopy as Octopus geometry include files named geometry-000, geometry-001, geometry-002, … These contain %LatticeParameters, %LatticeVectors and %ReducedCoordinates blocks in atomic units (bohr), ready to be pulled into an Octopus input file with the include directive.

Octopus uses atomic units, so within this interface the physical units are:

          | Distance   Atomic mass   Force         Force constants
-----------------------------------------------------------------
Octopus   | au (bohr)  AMU           hartree/au    hartree/au^2

The default displacement distance is 0.02 bohr. Phonon frequencies are reported in THz, as for the other interfaces.

Pre-process#

A bulk silicon primitive cell is used as the example throughout this page.

Exporting the unit cell from Octopus#

An Octopus calculation normally starts from an input file inp that already contains the crystal structure, for instance the cell obtained from a geometry optimization with Octopus. phonopy needs this unit cell as a VASP-style POSCAR file, which Octopus can write directly: add the geometry output in POSCAR format to the input file and run a ground-state calculation (CalculationMode = gs) on the (relaxed) structure.

A complete inp for the silicon primitive cell, using the parameters from the standard Octopus bulk-silicon tutorial, is:

# Bulk silicon in the diamond structure (2-atom primitive cell).
# Lengths given without a unit are in atomic units (bohr), the Octopus default.
CalculationMode = gs
PeriodicDimensions = 3

# a is the conventional cubic lattice constant of silicon.
a = 5.43*angstrom
%LatticeParameters
 a | a | a
%
%LatticeVectors
 0.0 | 0.5 | 0.5
 0.5 | 0.0 | 0.5
 0.5 | 0.5 | 0.0
%
%ReducedCoordinates
 "Si" | 0.0 | 0.0 | 0.0
 "Si" | 0.25 | 0.25 | 0.25
%

Spacing = 0.5
%KPointsGrid
 4 | 4 | 4
%

# Write the structure as a VASP-style POSCAR for phonopy.
%Output
 geometry
%
OutputFormat = poscar

This uses Octopus’ default pseudopotential set and LDA exchange-correlation functional, so no further settings are required to run it. The Spacing and %KPointsGrid values are the tutorial defaults; for a production calculation they (and the pseudopotentials/functional) should be converged for your system.

Run Octopus:

% octopus | tee out.log

Octopus writes the cell to static/POSCAR, with lattice vectors in Angstrom and fractional atomic positions. Copy it to the working directory as POSCAR and use it as the phonopy unit cell:

% cp static/POSCAR POSCAR

For the silicon input above, the exported POSCAR is:

Si
1.0
  0.0000000000   2.7150000000   2.7150000000
  2.7150000000   0.0000000000   2.7150000000
  2.7150000000   2.7150000000   0.0000000000
Si
2
Direct
  0.0000000000   0.0000000000   0.0000000000
  0.2500000000   0.2500000000   0.2500000000

Creating supercells with displacements#

In the pre-process, supercell structures with (or without) displacements are created from the unit cell, fully considering crystal symmetry.

To obtain supercells (\(2\times 2\times 2\)) with displacements, run phonopy-init with the –octopus option:

% phonopy-init --octopus -d --dim 2 2 2 --pa auto

You should find the files geometry-000, geometry-{number} and phonopy_disp.yaml:

% ls
geometry-000  geometry-001  phonopy_disp.yaml  POSCAR

geometry-000 is the perfect supercell structure, phonopy_disp.yaml contains the information on displacements, and geometry-{number} are the supercells with atomic displacements. Each geometry-{number} corresponds to one of the displacements written in phonopy_disp.yaml. Because of the high symmetry of the diamond structure, only a single displacement (geometry-001) is generated in this example.

A generated geometry-001 looks like:

%LatticeParameters
 14.511546484 | 14.511546484 | 14.511546484
%

%LatticeVectors
  0.000000000 | 0.707106781 | 0.707106781
  0.707106781 | 0.000000000 | 0.707106781
  0.707106781 | 0.707106781 | 0.000000000
%

%ReducedCoordinates
  "Si" | 0.000689106 | 0.000000000 | -0.000000000
  "Si" | 0.500000000 | 0.000000000 | -0.000000000
  ...
%

Calculation of sets of forces#

For each geometry-{number} file an Octopus ground-state calculation (CalculationMode = gs) is run to obtain the forces on the atoms of the displaced supercell. The geometry-{number} file is pulled into the Octopus input file inp with the include directive, together with the calculation settings that are appropriate for your system. It is convenient to run each displacement in its own directory:

% mkdir disp-001
% cd disp-001

An example inp for the silicon supercell (disp-001/inp) is:

CalculationMode = gs
PeriodicDimensions = 3
Spacing = 0.5*angstrom
BoxShape = parallelepiped

include ../geometry-001

%KPointsGrid
 2 | 2 | 2
%

Then run Octopus:

% octopus | tee out.log

Octopus writes the forces into static/info, in the “Forces on the ions” block (in Hartree/bohr by default):

Forces on the ions [H/b]
 Ion                        x              y              z
   1        Si  -1.00521691E-05  -1.11410816E-03  -1.11410816E-03
   2        Si   1.34774669E-03   4.66640181E-05   4.66640168E-05
   ...

Note

Be careful not to relax the structures. The atomic forces induced by the small displacement written in geometry-{number} are exactly what is needed for the phonon calculation, so the supercells with displacements must not be relaxed. Use a ground-state calculation (CalculationMode = gs, no geometry optimization).

Since the calculation is a supercell calculation, the convergence parameters (grid Spacing, %KPointsGrid, exchange-correlation functional, pseudopotentials, …) have to be chosen for your system. The settings above are only a minimal, fast example.

After the Octopus calculations of all displacements have finished, create the FORCE_SETS file with the -f option, passing the static/info files of the displacement calculations:

% phonopy-init -f disp-001/static/info

or, for several displacements,

% phonopy-init -f disp-{001..003}/static/info

The calculator (octopus) is read from phonopy_disp.yaml, so the --octopus option is not needed again here.

Post-process#

The post-processing is identical to the other calculators: it reads phonopy_disp.yaml and FORCE_SETS, so the --octopus option is not required (the calculator is stored in phonopy_disp.yaml).

In the post-process,

  1. Force constants are calculated from the sets of forces,

  2. A part of the dynamical matrix is built from the force constants,

  3. Phonon frequencies and eigenvectors are calculated from the dynamical matrices at the specified q-points.

The density of states (DOS) is plotted by

% phonopy --mesh 20 20 20 -p

Thermal properties are calculated with the sampling mesh by

% phonopy --mesh 20 20 20 -t

You should check the convergence with respect to the mesh numbers. Thermal properties can be plotted by

% phonopy --mesh 20 20 20 -t -p

Projected DOS is calculated and plotted by

% phonopy --mesh 20 20 20 --pdos "1 2, 3 4 5 6" -p

Band structure is plotted by

% phonopy --band "0.5 0.5 0.5  0.0 0.0 0.0  0.5 0.5 0.0  0.0 0.5 0.0" -p

In either case, by setting the -s option, the plot is going to be saved in the PDF format. If you don’t need to plot the DOS, the (partial) DOS is just calculated using the --dos option.

Non-analytical term correction (Optional)#

To activate the non-analytical term correction, a BORN file is required. It contains the macroscopic dielectric constant and the Born effective charges of the atoms in the primitive cell. Both quantities can be obtained from Octopus by a linear-response (Sternheimer) calculation of the electromagnetic response.

The BORN file can be generated automatically from the Octopus response calculation with the phonopy-octopus-born helper (analogous to phonopy-vasp-born, phonopy-qe-born and phonopy-crystal-born); it can also be assembled by hand, as described below. The Octopus output already uses the units expected by phonopy (Born charges in units of the elementary charge, dielectric tensor dimensionless), so the values are used as-is.

Note

The non-analytical term correction is only relevant for polar (ionic) crystals; for a non-polar crystal such as the silicon used in the rest of this page the Born effective charges vanish by symmetry and the correction has no effect. The example below therefore uses rock-salt NaCl, independently of the silicon force-calculation workflow above.

Computing Born charges and the dielectric tensor with Octopus#

The Born charges and the dielectric tensor are properties of the unit cell (not the supercell), and this calculation is run entirely with Octopus, so the structure is provided in the usual Octopus form. Here we use the rock-salt NaCl primitive cell, written as a geometry include file geometry-unitcell:

a = 5.64*angstrom
%LatticeParameters
 a | a | a
%
%LatticeVectors
 0.0 | 0.5 | 0.5
 0.5 | 0.0 | 0.5
 0.5 | 0.5 | 0.0
%
%ReducedCoordinates
 "Na" | 0.0 | 0.0 | 0.0
 "Cl" | 0.5 | 0.5 | 0.5
%

(If you only have the unit cell as a POSCAR, the same include file can be produced with phonopy-calc-convert -i POSCAR -o geometry-unitcell --calcin vasp --calcout octopus.)

For a periodic system the electric response is computed with the \(\vec{k}\cdot\vec{p}\) perturbation, so three runs are performed in the same directory: a ground state, a kdotp calculation, and finally the electromagnetic response. Each run restarts from the previous one.

The three runs share the same computational settings, which we collect in a file common.inp. The geometry is not included here; instead each run’s inp includes geometry-unitcell directly, alongside common.inp. This keeps the structure include at the top level — the same pattern as the force calculation above, where each disp-*/inp includes its own geometry-{number} — rather than nesting it inside common.inp:

PeriodicDimensions = 3
Spacing = 0.3*angstrom
BoxShape = parallelepiped
PseudopotentialSet = hgh_lda
%KPointsGrid
 4 | 4 | 4
%
KPointsUseSymmetries = no
ExperimentalFeatures = yes

  1. Ground state:

    CalculationMode = gs
include geometry-unitcell
include common.inp

  2. \(\vec{k}\cdot\vec{p}\) perturbation (required for the periodic electric response):

    CalculationMode = kdotp
KdotPCalcSecondOrder = yes
include geometry-unitcell
include common.inp

  3. Electromagnetic response with Born charges:

    CalculationMode = em_resp
RestartFixedOccupations = no
include geometry-unitcell
include common.inp
# Static (zero-frequency) response
%EMFreqs
 1 | 0.0
%
EMCalcBornCharges = yes

Run Octopus once per step (replacing the inp file between runs).

Note

Several constraints apply to this (experimental) calculation, hence ExperimentalFeatures = yes:

  • kdotp first. For a periodic system the electric response reads the \(\vec{k}\cdot\vec{p}\) wavefunctions, so the kdotp run must precede em_resp.

  • RestartFixedOccupations = no in the em_resp run, so that occupations are recomputed with semiconducting smearing (a gap is required for the \(\vec{k}\cdot\vec{p}\) electric response).

  • LDA only. The Sternheimer linear response used by em_resp evaluates the XC kernel without its gradient terms, so it currently supports only LDA functionals; a GGA functional stops with “GGA functionals are not allowed for now in XCKernel”. Use an LDA functional and matching pseudopotentials. (This restricts the em_resp path only; Octopus’ Casida TDDFT does handle GGA kernels.)

  • No nonlinear core corrections. The Born-charge force derivatives are not implemented for pseudopotentials with NLCC.

    The hgh_lda (Hartwigsen–Goedecker–Hutter LDA) set used above satisfies both the LDA and no-NLCC requirements. As with any response calculation, the grid spacing, k-point mesh and solver convergence have to be checked.

After the em_resp run, the results are written under em_resp/freq_0.0000/.

The macroscopic dielectric tensor \(\epsilon\) is in em_resp/freq_0.0000/epsilon:

# Real part of dielectric constant
        2.588022        0.000000       -0.000000
        0.000000        2.588022       -0.000000
       -0.000000       -0.000000        2.588022
Isotropic average        2.588022
...

The Born effective-charge tensors \(Z^*\) are in em_resp/freq_0.0000/born_charges, one \(3\times3\) tensor per atom:

# (Frequency-dependent) Born effective charge tensors
Index:     1   Label:    Na   Ionic charge:     1.0000
        1.141794        0.000000       -0.000000
        0.000000        1.141794        0.000000
        0.000000       -0.000000        1.141794
Isotropic average        1.141794

Index:     2   Label:    Cl   Ionic charge:     7.0000
       -1.141794        0.000000        0.000000
        0.000000       -1.141794       -0.000000
       -0.000000        0.000000       -1.141794
Isotropic average       -1.141794

# Discrepancy of Born effective charges from acoustic sum rule before correction, per atom
       -0.025770        0.000000       -0.000000
        0.000000       -0.025770        0.000000
        0.000000       -0.000000       -0.025770
Isotropic average       -0.025770

Each tensor is printed as three rows of the matrix (\(xx\,xy\,xz\) / \(yx\,yy\,yz\) / \(zx\,zy\,zz\)). The Isotropic average lines and the acoustic-sum-rule discrepancy block are not part of the BORN file.

The values are physically sensible: \(Z^*(\mathrm{Na}) \approx +1.14\), \(Z^*(\mathrm{Cl}) \approx -1.14\) and \(\epsilon_\infty \approx 2.59\), close to the experimental NaCl values (\(Z^* \approx \pm 1.1\), \(\epsilon_\infty \approx 2.3\)). The settings above are only a small example; check convergence for production use.

Assembling the BORN file#

The simplest route is the phonopy-octopus-born helper, which reads the dielectric tensor and Born charges from the em_resp output, symmetrizes them, keeps the symmetry-independent atoms, and writes the BORN file. Give it the unit cell and the em_resp results directory:

% phonopy-octopus-born geometry-unitcell em_resp/freq_0.0000 > BORN

The unit cell may be a POSCAR or an Octopus geometry file, but the geometry file must be in numeric form: the geometry-unitcell produced by phonopy-calc-convert above qualifies, whereas an input written with variables such as a = 5.64*angstrom does not (pass a POSCAR in that case). Passing the geometry file that was included in the em_resp run keeps the atom order aligned with born_charges.

Alternatively, the BORN file can be written by hand. Following the BORN format, write:

  1. the unit conversion factor on the first line (use the default by giving a non-numeric placeholder such as default, or the Octopus factor from Default unit conversion factor for non-analytical term correction),

  2. the nine dielectric-tensor components (\(xx\,xy\,xz\,yx\,yy\,yz\,zx\,zy\,zz\)) from epsilon on the second line,

  3. from the third line on, the nine \(Z^*\) components for each symmetry-independent atom of the primitive cell, taken from the corresponding Index: block of born_charges. The symmetry-independent atoms can be identified from the atom_mapping section printed by phonopy-init --octopus --symmetry -c POSCAR; for NaCl the two atoms (Na and Cl) are inequivalent.

For the NaCl output above, the resulting BORN file is:

default
2.588022 0.0 0.0 0.0 2.588022 0.0 0.0 0.0 2.588022
1.141794 0.0 0.0 0.0 1.141794 0.0 0.0 0.0 1.141794
-1.141794 0.0 0.0 0.0 -1.141794 0.0 0.0 0.0 -1.141794

Post-process with the non-analytical term correction#

With a BORN file in the working directory, phonopy activates the correction through the --nac option. This also requires the NaCl force constants: generate phonopy_disp.yaml and FORCE_SETS for NaCl with the same finite-displacement procedure used for silicon above (create the displaced supercells with phonopy-init --octopus -d, run the Octopus gs force calculations, and collect the forces with -f). With phonopy_disp.yaml, FORCE_SETS and BORN all present in the directory, add --nac to any post-process command.

The correction is most visible in the band structure, where it produces the LO–TO splitting of the optical branches at \(\Gamma\):

% phonopy --nac --band "0.5 0.5 0.5  0.0 0.0 0.0  0.5 0.5 0.0  0.0 0.5 0.0" -p

The DOS, thermal properties, and projected DOS (onto the Na and Cl atoms) are obtained exactly as before, with --nac added:

% phonopy --nac --mesh 20 20 20 -p
% phonopy --nac --mesh 20 20 20 -t
% phonopy --nac --mesh 20 20 20 --pdos "1, 2" -p

Without --nac the longitudinal- and transverse-optical modes stay degenerate at \(\Gamma\); the correction lifts that degeneracy.

Note

Octopus can alternatively compute Born effective charges and infrared intensities directly through its own linear-response phonon calculation (CalculationMode = vib_modes with CalcInfrared = yes), which writes them to vib_modes/infrared. That route is independent of the phonopy finite- displacement workflow described here.