First of all happy new year to everyone! I have received a lot of positive feedback on my posts and I apologize when it takes me a long time to get back to you – my new year’s resolution for 2020 is to post and reply more frequently 😉

Today I want to show two different ways of calculating phonons for molecules. The first one is by density functional perturbation theory (DFPT) as implemented in PWSCF using ph.x, the second one is by finite differences and using phonopy. Generally DFPT is fast and reliable but in the current version of PWSCF (6.5) it is still limited and can’t do hybrid functionals for example. This is where phonopy comes in. In the finite difference method the the equilibrium geometry is distorted, then a self-consistent calculation (pw.x) is performed and in the end all the forces are collected. This means that this method can be used whenever pw.x can be used which includes hybrid functionals (keyword *input_dft=’HSE’*). However finite differences cannot calculate the dielectric tensor or effective charges for which we are still limited to ph.x.

To calculate a molecule in a plane-wave code one uses the “molecule in a vacuum box” setup, where the molecule is sufficiently far away from its periodic images and this “isolated”. The size of the vacuum box needs to be converged but I find separations of about 15 Angstrom to be usually sufficient. This is shown on the right hand side where one periodic image in z-direction is included to emphasize the intrinsic periodic nature.

There are a few limitations of how the phonopy code reads the PWSCF in file and generates the displacement patterns.

- The only supported Bravais lattice option is currently
*ibrav=0*, which means that the*CELL_PARAMETERS*have to be input manually *CELL_PARAMETERS*have to be in Bohr units but can be converted from angstrom. However, be aware, phonopy currently only detects angstrom and not {angstrom} or (angstrom) – which are fine for PWSCF so be careful.

The workflow is as follows:

A) on GGA level perform a pw.x SCF calculation using *occupations=’fixed’* followed by a ph.x DFPT calculation where *epsil = .true.* . This will give the dielectric tensor as well as the effective charges which will serve as an input for phonopy calculations. The have to be copied in a file called BORN in a certain format (which will be done automatically by a primitive script if you run the attached example but you will have to adapt this for any other calculation). An explanation of the format is below. If phonopy detects symmetries in your system you don’t have to add all atoms in the BORN file – this can be tested by running phonopy with the -v key and check which atom has an asterisk. In this particular case all (12) atoms have to be added.

B) moving on to phonopy we need a PWSCF scf input file (which I named scf.in) and ask phonopy to create all displacements (which phonopy calls supercell-XXX.in) This is done by invoking

** phonopy –qe -d –dim=”1 1 1″ -c scf.in -v**

The options used here are –qe (using PWSCF), -d (create displacements) –dim=”1 1 1″ (the supercell has the same dimension as the original one since we are working with a vacuum box), -c scf.in (read the crystal structure from scf.in) and -v (verbose). The output will look something like

-------------------------------- super cell -------------------------------- Lattice vectors: a 19.314090000000000 0.000000000000000 0.000000000000000 b 0.000000000000000 19.963010000000001 0.000000000000000 c 0.000000000000000 0.000000000000000 15.000000000000000 Atomic positions (fractional): *1 C 0.43756745800000 0.53539552000000 0.50000000000000 12.011 > 1 *2 C 0.43711714400000 0.46541472000000 0.49999999900000 12.011 > 2 *3 C 0.50045009500000 0.56998119400000 0.50000000200000 12.011 > 3 *4 C 0.56285932100000 0.53460011500000 0.50000000300000 12.011 > 4 *5 C 0.56240056200000 0.46462296300000 0.50000000200000 12.011 > 5 *6 C 0.49952532300000 0.43003191600000 0.49999999900000 12.011 > 6 *7 H 0.38890741900000 0.56298252000000 0.50000000200000 1.008 > 7 *8 H 0.50083046500000 0.62454308600000 0.50000000200000 1.008 > 8 *9 H 0.38810367000000 0.43842476300000 0.49999999900000 1.008 > 9 *10 H 0.49916408600000 0.37546922500000 0.50000000000000 1.008 > 10 *11 H 0.61188283800000 0.56157435100000 0.50000000400000 1.008 > 11 *12 H 0.61105911200000 0.43703670400000 0.50000000300000 1.008 > 12 ---------------------------------------------------------------------------- "phonopy_disp.yaml" and supercells have been created.

and as mentioned above all atoms with an asterisk need to be included in the BORN file. At the same time phonopy will have created a series of displacements (48 in the case of benzene) which all contain only the displaced geometry but not the other input cards (&control, .. etc). The easiest way to use these is to append them to a “PWSCF template” of sorts and just add the coordinates. Then one has to run single-point force calculations using pw.x for each file which can be done at any hybrid exchange-correlation functional that is currently available – neat! Once those are done phonopy can collect the forces from all the separate out-files by using

phonopy –qe -f supercell-{001..048}.out –nac

This will create a file called FORCE_SETS and it will read the BORN file to add the “non-analytical correction” (–nac).

To check the phonon frequencies one can easily grep “frequency” from the mesh.yaml file which in this particula case will look like

frequency: -0.8565677184 frequency: -0.4716156352 frequency: -0.3672091881 frequency: -0.0000013963 frequency: -0.0000008147 frequency: -0.0000002343 frequency: 16.2961123012 frequency: 16.3088957474 frequency: 24.8072360316 frequency: 24.8201310171 frequency: 27.2807808008 frequency: 29.0709809006 frequency: 34.4192607530 frequency: 34.5190166527 frequency: 39.4455527017 frequency: 39.5271778453 frequency: 40.6711523243 frequency: 40.9449176736 frequency: 41.1184463866 frequency: 42.6869160485 frequency: 42.7228524891 frequency: 47.0752882284 frequency: 48.0036987657 frequency: 48.0300599224 frequency: 55.1675013606 frequency: 55.3599357594 frequency: 60.4965619592 frequency: 60.5350132184 frequency: 65.5797175182 frequency: 65.5925203426 frequency: 127.3836904483 frequency: 127.7387235024 frequency: 127.8205031638 frequency: 128.3666018221 frequency: 128.4551561033 frequency: 128.8032547943

negative frequencies are imaginary modes. If they are small enough (here all less than 1 THz ~ 4 meV) one can ignore them, if they are large the structure was probably not sufficiently relaxed or the calculation was not converged in terms of the cutoff, vacuum box size, ..). With the phonopy-spectroscopy code one can also easily calculate nice (broadened) infrared spectra with one simple line:

phonopy-ir –linewidth=16.5 –spectrum_range=”0.0 4000″

For visual comparison of many spectra this is very convenient so check it out!

Finally here it is: the spectra of benzene using different exchange-correlation functionals compared to experiments (from the NIST database) and the DFPT results.

Was this post helpful or am I missing something? Please let me know in the comments or via mail!

To reproduce the results just run the attached run.sh file. Have fun! run_phonopy_qe.zip