Compute Irreducible Derivative

Within this package we can compute irreducible derivatives of phonons and phonon interactions at a given order with either LID or BID approach. The theory behind these approaches are described in the following paper:

• L. Fu, M. Kornbluth, Z. Cheng, and C. A. Marianetti, Phys. Rev. B 100, 014303 (2019).

Specifically, we can execute the LID, SS-BID and HS-BID approaches from the command line using the following commands:

1. pm-lid-mesh: compute irreducible derivatives of a given FTG and order using LID method;

2. pm-bid: compute irreducible derivatives of a given FTG and order using SS-BID method;

3. pm-hsbid: compute irreducible derivatives of a given FTG and order using HD-BID method.

These commands take in generally the same command line options, therefore one can switch between approaches by changing the command while keeping the options the same. For example, below is guide on using the LID approach to compute irreducible derivatives of NaCl at second order but one can replace pm-lid-mesh with pm-bid or pm-hsbid to execute SS-BID or HS-BID approach.

LID Approach on NaCl at Second Order

There are many ways to use this package to execute calculations of the irreducible derivatives with LID approach. Here we demonstrate a recommended workflow to execute our LID approach at second order on NaCl from forces computed by first principle method.

Symmetry Analysis and Find the Measurements to Compute

First we create a file xtal.yml describing the structure of NaCl:

# xtal.yml
vec:
  [[ 0.00000000,  2.83000729,  2.83000729],
   [ 2.83000729,  0.00000000,  2.83000729],
   [ 2.83000729,  2.83000729,  0.00000000]]
atoms:
  - Na: 
      [[ 0.00000000,  0.00000000,  0.00000000]]
  - Cl: 
      [[ 0.50000000,  0.50000000,  0.50000000]]

Then, with this structure and the point group \(O_h\), we can write the minimal configuration file, named config.yml, to execute the LID approach with a supercell \(\hat{S}_{BZ}=2\hat1\):

# config.yml
structure: 'xtal.yml'
supa: |
  2 0 0
  0 2 0
  0 0 2  
pg: 'Oh'
order:  2
full_symmetry: True
delta: [ 0.02, 0.03, 0.04]

There are numerous parameters which have default settings that can be changed, and these can be added to the config.yml file listed above:

fdtype: 'c'                  
root_directory: '.'
database: 'database.db'
db_type: 'sqlite'
db_table: 'lid_phonon'

With the configuration file we can initiate the pre-processing step which will execute the symmetry analysis and create the jobs that need to be computed with first principles:

pm-lid-mesh --config config.yml --create 

The above step will create two files:

• lid_mesh.hdf5: A HDF5 file containing the computed symmetry analysis results, which will be used in the post-processing step to compute the irreducible derivatives after first principle calculations are finished.

• database.db: A SQLite database file containing the jobs need to be computed with first principle calculations (db_type is specified to be ‘sqlite’ in the config.yml file, other database can be used if a working interface is provided).

      pm-lid-mesh --save-config --create  --structure xtal.yml \
        --supa '2,0,0;0,2,0;0,0,2' --pg Oh --order 2 --full-symmetry \
        --delta 0.02,0.03,0.04 --fdtype c --database database.db \
        --db-type sqlite --db-table lid_phonon

Prepare and Execute First Principle Calculations

The next step is to create the first principle jobs to be computed, the package generically supports all types of first principle compute engines provided a working interface, and here we demonstrate using VASP. We begin by create the VASP jobs to be executed:

  pm-jobs --config config.yml --create --engine vasp --save-config jobs_config.yml \
    --job-config-path path-to-vasp-input-templates --root-directory jobs-ro-run

Inside the path-to-vasp-input-templates directory, one usually needs to provide “INCAR”, “KPOINTS” and “POTCAR” files and they will be used as templates when creating the vasp jobs, each vasp job will be created in a separate directory for convenience. Once these directories are created, VASP can be executed within each directory. Afterwards, the calculation results can be collected and stored in the database:

  pm-jobs --config jobs_config.yml --retrieve

Compute Irreducible Derivatives from the Results of First Principle Calculations

Finally, we can execute the post-processing step of the LID approach to compute the irreducible derivatives and, for example, compute the force constants tensors for Fourier interpolation:

  pm-lid-mesh --load-lid-mesh lid_mesh.hdf5 --config config.yml \
    --post --save-fi --save-ids --save-dts

If the material of interest is an insulator and the dipole-dipole interaction needs to be considered, the dataset for computing the dipole-dipole interaction needs to be provided as an input as well:

  pm-lid-mesh --load-lid-mesh lid_mesh.hdf5 --config config.yml --epsilon path-to-epsilon.yml \
    --post --save-fi --save-ids --save-dts

The dataset for computing the dipole-dipole interaction can be computed from first principle methods. Continuing with the VASP example, we can obtain the dataset by executing a static calculation adding LEPSILON=.TRUE. to the “ICNAR” file. Once the calculation is finished, use the following command inside the calculation directory to parse the output from VASP and create the epsilon.yml file:

  pm-epsilon --vasp -o epsilon.yml .

then replace path-to-epsilon.yml with the path to the generated epsilon.yml file.

The following files are created in this step:

• irreducible_derivatives.hdf5: A HDF5 files containing the names of the irreducible derivatives and their values.

• fourier_interpolation.hdf5: A HDF5 files containing all the necessary information (structure, force constants tensors, FTG and dipole-dipole contributions if at second order) to compute Fourier interpolation.

• dynamic_tensors.hdf5: A HDF5 files containing all the dynamic tensors in the naive basis within irreducible Q.

Computing and analyzing error tails

The results of any finite displacement approach are only reliable if the discretization error has been properly controlled. Forward and backward finite difference guarantee a linear error tail, while central finite difference guarantees a quadratic error tail. Resolving the error tail is a sure sign that the discretization error has been properly extrapolated to zero. For phonons, we will be taking either the first derivatives of the forces or the second derivatives of the energy, which are typically not sensitive. Therefore, it is common practice to use a single finite displacement (i.e. a single delta), and then there is nothing to analyze. However, there are plenty of situations where one does find sensitive results, such as metallic systems with Fermi surface nesting or even band insulators when using the SCAN functional. If one computes three or more deltas, then the error tail can be scrutinized to judge the quality of the result.

To examine the error tails, the corresponding data must be output during the post-processing of results. The above lid-mesh commands needs the –err-output flag as follows:

  pm-lid-mesh --load-lid-mesh lid_mesh.hdf5 --config config.yml \
    --post --save-fi --save-ids --save-dts --err-output error_tails.hdf5

The next step is extract the results from the hdf5 file and output them to individual yaml files:

pm-convert-errortail error_tails.hdf5

The yaml files can then be plotted to PDF files. For example:

pm-plot-errortail 1_0_errortail.yml --iscomplex --output 1_0_errortail.pdf

At present, one can execute a loop in bash over all the files:

for i in *_errortail.yml; do
  echo $i
  pm-plot-errortail $i --iscomplex --output $(basename $i .yml).pdf
done

Computing phonon dispersion and DOS

With the fourier_interpolation.hdf5 file obtained at second order, we can compute phonon dispersion and density of states (DOS).

Let’s continue the NaCl example and first prepare the configuration file phonon_config.yml for computing the phonon dispersion and DOS:

# phonon_config.yml

band:
  - ['0 0 0', '1/2 1/2 0']
  - ['-1/2 1/2 0', '0 0 0']
  - ['0 0 0', '1/2 1/2 1/2']
band_label:
  - $\Gamma$
  - X
  - $\Gamma$
  - L
mesh: 20 20 20

Then we can execute the following command to compute the band structure and DOS:

  pm-phonon-band-dos --load-fi fourier_interpolation.hdf5 --config phonon_config.yml \
    --save-band --save-grid --save-dos --save-plot

The following files are created in this step:

• band_plot_config.yml: A YAML file with the configuration the plot the band structure and DOS;

• bands.dat: A text file with the phonon band structure data computed with Fourier interpolation;

• grid_phonons.dat: A text file with the data of the phonons that are measured at the FTG that fall on the path in the plot;

• dos.dat: A text file with the data of phonon DOS.

We can plot the phonon band structure and DOS with the following command:

  pm-plot-band -s band_plot_config.yml