Tag Archive: phonons

Jun 09

Bachelor Projects Completed: 2 new computational materials scientists initialised

The black arts of computational materials science.

Black arts of computational materials science.

Just over half a year ago, I mentioned that I presented two computational materials science related projects for the third bachelor physics students at the UHasselt. Both projects ended up being chosen by a bachelor student, so I had the pleasure of guiding two eager young minds in their first steps into the world of computational materials science. They worked very hard, cursed their machine or code (as any good computational scientist should do once in a while, just to make sure that he/she is still at the forefront of science) and survived. They actually did quite a bit more than “just surviving”, they grew as scientists and they grew in self-confidence…given time I believe they may even thrive within this field of research.

One week ago, they presented their results in a final presentation for their classmates and supervisors. The self-confidence of Giel, and the clarity of his story was impressive. Giel has a knack for storytelling in (a true Pan Narrans as Terry Pratchett would praise him). His report included an introduction to various topics of solid state physics and computational materials science in which you never notice how complicated the topic actually is. He just takes you along for the ride, and the story unfolds in a very natural fashion. This shows how well he understands what he is writing about.

This, in no way means his project was simple or easy. Quite soon, at the start of his project Giel actually ran into a previously unknown VASP bug. He had to play with spin-configurations of defects and of course bumped into a hand full of rookie mistakes which he only made once *thumbs-up*. (I could have warned him for them, but I believe people learn more if they bump their heads themselves. This project provided the perfect opportunity to do so in a safe environment. ūüėé ) ¬†His end report was impressive and his results on the Ge-defect in diamond are of very good quality.

The second project was brought to a successful completion by Asja. This very eager student actually had to learn how to program in fortran before he could even start. He had to implement code to calculate partial phonon densities with the existing HIVE code. Along the way he also discovered some minor bugs (Thank you very much ūüôā ¬†) and crashed into a rather unexpected hard one near the end of the project. For some time, things looked very bleak indeed: the partial density of equivalent atoms was different, and the sum of all partial densities did not sum to the total density. As a result there grew some doubts if it would be possible to even fulfill the goal of the project. Luckily, Asja never gave up and stayed positive, and after half a day of debugging on my part the culprit was found (in my part of the code as well). Fixing this he quickly started torturing his own laptop calculating partial phonon densities of state for Metal-organic frameworks and later-on also the Ge-defect in diamond, with data provided by Giel. Also these results are very promising and will require some further digging, but they will definitely be very interesting.

For me, it has been an interesting experience, and I count myself lucky with these two brave and very committed students. I wish them all the best of luck for the future, and maybe we meet again.

Jun 27

Folding Phonons

Game-diamondsAbout a year ago, I discussed the possibility of calculating phonons (the collective vibration of atoms) in the entire Brillouin zone for Metal-Organic Frameworks. Now, one year later, I return to this topic, but this time the subject matter is diamond. In contrast to Metal-Organic Frameworks, the unit-cell of diamond is very small (only 2 atoms). Because a phonon spectrum is calculated through the gradients of forces felt by one atom due to all other atoms, it is clear that within one diamond unit-cell these forces will not be converged. As such, a supercell will be needed to make sure the contribution, due to the most distant atoms, to the experienced forces, are negligible.

Using such a supercell has the unfortunate drawback that the dynamical matrix (which is 3N \times 3N, for N atoms) explodes in size, and, more importantly, that the number of eigenvalues, or phonon-frequencies also increases (3N) where we only want to have 6 frequencies ( 3 \times 2 atoms) for diamond. For an M \times M \times M supercell we end up with 24M^3 -6  additional phonon bands which are the result of band-folding. Or put differently, 24M^3 -6 phonon bands coming from the other unit-cells in the supercell. This is not a problem when calculating the phonon density of states. It is, however, a problem when one is interested in the phonon band structure.

The phonon spectrum at a specific q-point in the first Brillouin zone is given by the square root of the eigenvalues of the dynamical matrix of the system. For simplicity, we first assume a finite system of n atoms (a molecule). In that case, the first Brillouin zone is reduced to a single point q=(0,0,0) and the dynamical matrix looks more or less like the hessian:

With \varphi (N_a,N_b) = [\varphi_{i,j}(N_a , N_b)] 3 \times 3 matrices \varphi_{i,j}(N_a,N_b)=\frac{\partial^2\varphi}{\partial x_i(N_a) \partial x_j(N_b)} = - \frac{\partial F_i (N_a)}{\partial x_j (N_b)}  with i, j = x, y, z. Or in words, \varphi_{i,j}(N_a , N_b) represents the derivative of the force felt by atom N_a due to the displacement of atom N_b. Due to Newton’s second law, the dynamical matrix is expected to be symmetric.

When the system under study is no longer a molecule or a finite cluster, but an infinite solid, things get a bit more complicated. For such a solid, we only consider the symmetry in-equivalent atoms (in practice this is often a unit-cell). Because the first Brillouin zone is no longer a single point, one needs to sample multiple different points to get the phonon density-of-states. The role of the q-point is introduced in the dynamical matrix through a factor e^{iq \cdot (r_{N_a} - r_{N_b}) }, creating a dynamical matrix for a single unit-cell containing n atoms:

Because a real solid contains more than a single unit-cell, one should also take into account the interactions of the atoms of one unit-cell with those of all other unit-cells in the system, and as such the dynamical matrix becomes a sum of matrices like the one above:

Where the sum runs over all unit-cells in the system, and Ni indicates an atom in a specific reference unit-cell, and MRi  an atom in the Rth unit-cells, for which we give index 1 to the reference unit-cell. As the forces decay with the distance between the atoms, the infinite sum can be truncated. For a Metal-Organic Framework a unit-cell will quite often suffice. For diamond, however, a larger cell is needed.

An interesting aspect to the dynamical matrix above is that all matrix-elements for a sum over n unit-cells are also present in a single dynamical matrix for a supercell containing these n unit-cells. It becomes even more interesting if one notices that due to translational symmetry one does not need to calculate all elements of the entire supercell dynamical matrix to construct the full supercell dynamical matrix.

Assume a 2D 2×2 supercell with only a single atom present, which we represent as in the figure on the right. A single periodic copy of the supercell is added in each direction. The dynamical matrix for the supercell can now be constructed as follows: Calculate the elements of the first column (i.e. the gradient of the force felt by the atom in the reference unit-cell, in black, due to the atoms in each of the unit-cells in the supercell). Due to Newton’s third law (action = reaction), this first column and row will have the same elements (middle panel).

Translational symmetry on the other hand will allow us to determine all other elements. The most simple are the diagonal elements, which represent the self-interaction (so all are black squares). The other you can just as easily determine by looking at the schematic representation of the supercell under periodic boundary conditions. For example, to find the derivative of the force on the second cell (=second column, green square in supercell) due to the third cell (third row, blue square in supercell), we look at the square in the same relative position of the blue square to the green square, when starting from the black square: which is the red square (If you read this a couple of times it will start to make sense). Like this, the dynamical matrix of the entire supercell can be constructed.

This final supercell dynamical matrix can, with the same ease, be folded back into the sum of unit-cell dynamical matrices (it becomes an extended lookup-table). The resulting unit-cell dynamical matrix can then be used to create a band structure, which in my case was nicely converged for a 4x4x4 supercell. The bandstructure along high symmetry lines is shown below, but remember that these are actually 3D surfaces. A nice video of the evolution of the first acoustic band (i.e. lowest band) as function of its energy can be found here.

The phonon density of states can also be obtained in two ways, which should, in contrast to the bandstructure, give the exact same result: (for an M \times M \times M supercell with n atoms per unit-cell)

  1. Generate the density of states for the supercell and corresponding Brillouin zone. This has the advantage that the smaller Brillouin zone can be sampled with fewer q-points, as each q-point acts as M3 q-points in a unit-cell-approach. The drawback here is the fact that for each q-point a (3nM3)x(3nM3) dynamical matrix needs to be solved. This solution scale approximately as O(N3) ~ (3nM3)3 =(3n)3M9. Using linear algebra packages such as LAPACK, this may be done slightly more efficient (but you will not get O(N2) for example)
  2. Generate the density of states for the unit-cell and corresponding Brillouin zone. In this approach, the dynamical matrix to solve is more complex to construct (due to the sum which needs to be taken) but much smaller: 3nx3n. However to get the same q-point density, you will need to calculate M3 times as many q-points as for the supercell.

In the end, the choice will be based on whether you are limited by the accessible memory (when running a 32-bit application, the number of q-point will be detrimental) or CPU-time (solving the dynamical matrix quickly becomes very expensive).

 

Jun 19

Phonons: shake those atoms

In physics, a phonon is a collective excitation in a periodic, elastic arrangement of atoms or molecules in condensed matter, like solids and some liquids. Often designated a quasi-particle, it represents an excited state in the quantum mechanical quantization of the modes of vibrations of elastic structures of interacting particles.

‚ÄĒ source:¬†wikipedia

Or for simplicity: sound waves; the ordered shaking of atoms or molecules. When you hit a metal bell with a (small) hammer, or make a wineglass sing by rubbing the edges, you experience the vibrations of the object as sound. All objects in nature (going from atoms to stars) can be made to vibrate, and they do this at one or more specific frequencies : their eigenfrequencies or normal frequencies.

Also single molecules,¬†if they are hit (for example by another molecule bumping into them) or receive extra energy in another way, start to vibrate. These vibrations can take many forms (elongating and shortening of bonds, rotating of parts of the molecule with respect to other parts, flip-flopping of loose ends, and so forth) and give a unique signature to the molecule since each of these vibrations (so-called eigen-modes) corresponds with a certain energy given to the molecule. As a result, if you know all the eigen-modes of a molecule, you also know which frequencies of infrared light they should absorb, which is very useful, since in experiment we do not “see” molecules (if we see them at all) as nice ball-and-stick objects.

From the computational point of view, this is not the only reason why in molecular modeling the vibrational frequencies of a system (i.e. the above eigen-modes) are calculated. In addition, they also tell if a system is in its ground state (which is what one is looking for most of the time) or not. Although this tool has¬†wide-spread usage in molecular modeling, it is seldom used in ab initio solid state physics because of the associated computational cost. In addition, because of the finite size of the unit cell, the reciprocal space in which phonons¬†live also has a finite size, in contrast to the single point for a molecule…making life complex. ūüėé

Read the rest of this entry »