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.

1. The only supported Bravais lattice option is currently ibrav=0, which means that the CELL_PARAMETERS have to be input manually
2. 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