:sequential_nav: next .. _tutorial-advanced-phonons: Advanced topics in phonons ========================== :author: Andrei Postnikov (Univ. Lorraine) .. note:: For background, see the excellent `slide presentation `_ for a lecture by Andrei Postnikov in a previous Siesta school. .. note:: There is some extra material for this tutorial `here `_. .. note:: To fully appreciate the issues in this tutorial you need to have understood the concepts presented in the :ref:`basic vibrations tutorial`. .. This intermediate/advanced tutorial will focus on topics such as: * Post-processing with Andrei Postnikov's lattice-dynamics tools. * Phonon (projected)DOS * Born-effective charges and infrared activity. .. Exercises: Use the newer examples by Andrei on a BN monolayer and a Au chain. .. For LO-TO and infrared, we need examples. Actually, for LO-TO splitting we need to fish out Thomas Archer's vibra code. The supercell explosion ----------------------- We are going to consider a monolayer of BN, with two atoms per unit cell. In directory *Phonon_Q* we have already generated a 5x5x1 supercell with *fcbuild* (see the :ref:`basic vibrations tutorial`), by specifying:: Supercell_1 2 Supercell_2 2 Supercell_3 0 The supercell contains 50 atoms, and we should do in principle 2x6=12 calculations displacing the two atoms in the unit cell in the positive and negative cartesian directions to compute the forces on all the atoms in the supercell. That is a lot of work, so we have prepared a FC file for the rest of the exercise. You need to copy this pre-generated FC file: ``cp BN.FC_save BN.FC``, and then run vibra: ``vibra < BN_vibra.fdf``. You will get the *BN.bands* and *BN.vectors* files. The *BN.bands* file holds the phonon dispersion information, that you can plot with:: gnubands -G -o bnq BN.bands gnuplot --persist bnq.gplot You will get something like .. figure:: images/bnq-dispersion.png :align: center which looks quite bad in a neighborhood of Gamma. Obviously, to improve on this we should increase the size of the supercell, maybe going to 7x7x1, 9x9x1, or even larger. The prospect is not pleasing. But there is another way. Instead of generating a "brute force" supercell, we can generate a supercell that will describe correctly the dispersion along the high-symmetry directions that we care about. In this particular case, we can construct, following the indications in the slide presentation, a supercell that will exactly map the K and M points to Gamma. This is not the supercell for bookkeeping of the force-constants, but a new "unit cell" that serves our purposes. The process starts in the BZ, chosing new unit vectors to define a lattice that will include M and K: .. figure:: images/new-bz.png :align: center The reciprocal cell of this new unit cell is the real-space supercell we seek: .. figure:: images/new-cell.png :align: center In directory *Phonon_GMK* we have the necessary files to compute the phonons **at Gamma** in this new "enlarged unit cell". Now, with 24 atoms, we would need 24*6=144 calculations. It is a bit more CPU than the original case above (12 calculations for 50 atoms, each about 8 times more expensive than one for 24 atoms). Instead of running the calculations, you can use the provided FC file. (Copy the *BN_save.* files to use simply the prefix *BN*) Run ``vibq0`` and use the *BN* SystemLabel when asked. It will read the `BN.FC` and `BN.XV` files and produce *.bands* and *.vectors* files, just like ``vibra``, but with a much simpler interface when only the Gamma point is handled. In this case, the plotting is different because there is a single q-point:: gnubands -G -o gamma BN.bands gnuplot gnuplot> plot "gamma" using 1:2 with points and you will get a column of points marking the frequencies of the zone-center phonons in this enlarged unit cell. Not very useful compared to the original dispersion. All the modes are piled up at Gamma, and we need to **unfold** them. .. note:: There are fundamental similarities between this problem and that of the unfolding of electronic bands covered in :ref:`this tutorial`, but the techniques are different. Unfolding the Gamma "phonon structure" in the supercell ....................................................... The character of these different modes can be revealed by projecting their eigenvectors onto waves with different Q-vector (see the lecture liked-to above): .. figure:: images/spectral.png :align: center The projection is done by the ``phdos`` program which is part of the lattice dynamics analysis package. In order to probe projections onto different symmetric points of the BZ, the following coordinates of Q-points can be specified (in Cartesian coordinates and units of 1/Bohr, omitting the 2*pi factor). For A=2.495 Bohr:: M [0, 2*pi/(a*sqrt(30), 0] -> [0, 0.2314029, 0], or M [pi/a, pi/(a*sqrt(3)), 0] -> [0.2004008, 0.1157014, 0]; K [2*pi/(3*a), 2*pi/(a*sqrt(3)), 0] -> [0.1336005, 0.2314029, 0]; K [4*pi/(3*a), 0, 0] -> [0.2672011, 0, 0] according to the picture: .. figure:: images/paths-bn.png :align: center We have prepared the *G.inp*, *K.inp*, *M.inp* files to run the program. For example:: phdos < G.inp mv BN.WSQ BN.G.WSQ mv BN.VSQ BN.G.VSQ will produce a *.VSQ* and a *.WSQ* file, which we rename for convenience and clarity. To plot the contents of the second one, which is more useful:: gnuplot gnuplot> plot "BN.G.WSQ" u 1:3 w l gnuplot> plot "BN.G.WSQ" u 1:2 w l We will get sharp peaks at the positions of the "true" Gamma phonons (three acoustic modes with zero frequency, an intermediate one, and a degenerate optical mode), as in this figure: .. figure:: images/BN-GMK-Gamma-WSQ.png :align: center We can do the same for the K and M points. For K: .. figure:: images/BN-GMK-K-WSQ.png :align: center For M: .. figure:: images/BN-GMK-M-WSQ.png :align: center The two different colors correspond to the B and N atoms' role in each mode. Being lighter, B will feature more prominently in the higher-frequency modes. Can you check whether the phonons at K and M match the values in the general dispersion in the first section? Can you think of which q-vectors to use with ``phdos`` to see whether we get better results in the lines Gamma-K and Gamma-M? .. note:: One would be tempted to get a "dispersion" by selecting arbitrary q-vectors along a line and feeding them to ``phdos``. However, it should be noted that we have really computed only 3N modes, where N is the number of atoms in the supercell. For some *q*'s we will get sharp peaks, but for others the signal will degrade and it would be difficult to make out "frequencies". .. hint:: The *.VSQ* files might be useful to get the frequencies of the peaks, as they have no broadening. .. note:: There is an alternative method to compute the phonon dispersion at arbitrary *q*s in the BZ, using only *unit cell* calculations. It is based on Density-Functional Perturbation Theory. While convenient, this route removes half the fun of understanding supercells and folding issues. LO-TO splitting and Born dynamical charges -------------------------------------------- Here is a picture of the phonon dispersion in hexagonal BN: .. figure:: images/h-bn-dispersion.png :align: center The LO-TO splittings at Gamma, in a *3D* system, are caused by the effect of the electric field that the optical modes produce. There is an extra term in the force constants: .. figure:: images/lo-to-term.png :align: center The second term is non-analytic, in the sense that it has different limiting values when q-->0 depending on the direction in which we approach it. The *Z*'s in the above equation are the Born effective or dynamical charges, and can be computed as the derivative of the polarization in the crystal with respect to the atomic displacement, as shown in :ref:`this tutorial`, which you might enjoy following. In fact, the Born charges for hexagonal BN are explicitly calculated there, so you could go ahead and estimate roughly the size of the extra terms in the dynamical matrix for specific directions. There is one extra ingredient: the :math:`\varepsilon_\infty` in the non-analytic term is the clamped-ion dielectric constant, which can be computed using the methods outlined in the :ref:`tutorial on optical properties`. .. Miguel Pruneda says: In a 2D system there is no LO-TO splitting. The LO-TO comes from planes of atoms generating dipoles, which give constant electric fields in space, to which the ions are coupled, hence modifying the LO branch. In 2D systems the dipoles form 1D lines, and the field has a different long-range behaviour. I learnt this from an old paper by Eduardo Hernandez and Daniel Sanchez-Portal (10.1103/PhysRevB.66.235415), and more recently Nicola Marzari & co published another work (10.1021/acs.nanolett.7b01090) We will not carry out any more specific calculations here, but it is worth noting the subtle links of the topic of phonons with other domains. Also, it is apparent that to fully analyze a given system you might have to do several kinds of simulations: force-constants, Born charges, dielectric constant, etc. Ideally, you would set up an automated **workflow** that can help you manage everything. See the :ref:`AiiDA framework tutorial` for more information. .. note:: Another use of the Born charges is in computing the infrared intensities of the modes. Only a few modes (those which involve changes in the polarization) would be infrared-active. The ``vibra`` program can compute the infrared intensities if both the FC and the BC (for Born charges) files are available. See the Siesta and Vibra manuals for more information.