# Tag Archive: Metal-Organic Frameworks

Jul 27

## Book chapter: Computational Chemistry Experiment Possibilities

 Authors: Bartłomiej M. Szyja and Danny Vanpoucke Book: Zeolites and Metal-Organic Frameworks, (2018) Chapter Ch 9, p 235-264 Title Computational Chemistry Experiment Possibilities ISBN: 978-94-629-8556-8 export: bibtex pdf:

 Zeolites and Metal-Organic Frameworks (the hard-copy)

## Abstract

Thanks to a rapid increase in the computational power of modern CPUs, computational methods have become a standard tool for the investigation of physico-chemical phenomena in many areas of chemistry and technology. The area of porous frameworks, such as zeolites, metal-organic frameworks (MOFs) and covalent-organic frameworks (COFs), is not different. Computer simulations make it possible, not only to verify the results of the experiments, but even to predict previously inexistent materials that will present the desired experimental properties. Furthermore, computational research of materials provides the tools necessary to obtain fundamental insight into details that are often not accessible to physical experiments.

The methodology used in these simulations is quite specific because of the special character of the materials themselves. However, within the field of porous frameworks, density functional theory (DFT) and force fields (FF)
are the main actors. These methods form the basis of most computational studies, since they allow the evaluation of the potential energy surface (PES) of the system.

Newsflash: here

Jan 19

## Newsflash: Book-chapter on MOFs and Zeolites en route to bookstores near you.

It is almost a year ago that I wrote a book-chapter, together with Bartek Szyja, on MOFs and Zeolites. Coming March 2018, the book will be available through University press. It is interesting to note that in a 13 chapter book, ours was the only chapter dealing with the computational study and simulation of these materials…so there is a lot more that can be done by those who are interested and have the patience to perform these delicate and often difficult but extremely rewarding studies. From my time as a MOF researcher I have learned two important things:

1. Any kind of interesting/extreme/silly physics you can imagine will be present in some MOFs. In this regard, the current state of the MOF/COF field is still in its infancy as most experimental work focuses on  simple applications such as catalysis and gas storage, for which other materials may be better suited. These porous materials may be theoretically interesting for direct industrial application, but the synthesis cost generally will be a bottleneck. Instead, looking toward the fundamental physics applications: Low dimensional magnetism, low dimensional conduction, spin-filters, multiferroics, electron-phonon interactions, interactions between spin and mechanical properties,…. MOFs are a true playground for the theoretician.
2. MOFs are very hard to simulate correctly, so be wary of all (published) results that come computationally cheap and easy. Although the unit-cell of any MOF is huge, with regard to standard solid state materials, the electron interactions are also quite long range, so the first Brillouin zone needs very accurate sampling (something often neglected). Also spin-configurations can have a huge influence, especially in systems with a rather flat potential energy surface.

In the book-chapter, we discuss some basic techniques used in the computational study of MOFs, COFs, and Zeolites, which will be of interest to researchers starting in the field. We discuss molecular dynamics and Monte Carlo, as well as Density Functional Theory and all its benefits and limitations.

Jun 09

## Bachelor Projects Completed: 2 new computational materials scientists initialised

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 07

## VSC-user day 2017: The Poster Edition

Last Friday, the HPC infrastructure in Flanders got celebrated by the VSC user day. Being one of the Tier-1 supercomputer users at UHasselt, I was asked if I could present a poster at the meeting, showcasing the things I do here. Although I was very interested in this event, educational obligations (the presentations of the bachelor projects, on which I will post later) prevented me from attending the meeting.

As means of a compromise, I created a poster for the meeting which Geert Jan Bex, our local VSC/HPC support team, would be so nice to put up at the event. The poster session was preceded by a set of 1-minute presentations of the posters, for which a slide had to be made. As I could not be physically present, I provided the organizers a slide which contained a short description that could be used as the 1-minute presentation. Unfortunately, things got a little mixed up, as Geert Jan accidentally printed this slide as the poster (which gave rise to some difficulties in the printing process 🙄 ). So for those who might have had an interest in the actual poster, let me put it up here:

This poster presents my work on linker functionalisation of the MIL-47, which got recently published in the Journal of physical chemistry C, and the diamond work on the C-vacancy, which is currently submitted. Clicking on the poster above will provide you the full size image. The 1-minute slide presentation, which erroneously got printed as poster:

Mar 30

## Linker Functionalization in MIL-47(V)-R Metal–Organic Frameworks: Understanding the Electronic Structure

 Authors: Danny E. P. Vanpoucke Journal: J. Phys. Chem. C 121(14), 8014-8022 (2017) doi: 10.1021/acs.jpcc.7b01491 IF(2017): 4.484 export: bibtex pdf:
 Graphical Abstract: Evolution of the electronic band structure of MIL-47(V) upon OH-functionalization of the BDC linker. The π-orbital of the BDC linker splits upon functionalisation, and the split-off π-band moves up into the band gap, effectively reducing the latter.

## Abstract

Metal–organic frameworks (MOFs) have gained much interest due to their intrinsic tunable nature. In this work, we study how linker functionalization modifies the electronic structure of the host MOF, more specifically, the MIL-47(V)-R (R = −F, −Cl, −Br, −OH, −CH3, −CF3, and −OCH3). It is shown that the presence of a functional group leads to a splitting of the π orbital on the linker. Moreover, the upward shift of the split-off π-band correlates well with the electron-withdrawing/donating nature of the functional groups. For halide functional groups the presence of lone-pair back-donation is corroborated by calculated Hirshfeld-I charges. In the case of the ferromagnetic configuration of the host MIL-47(V+IV) material a half-metal to insulator transition is noted for the −Br, −OCH3, and −OH functional groups, while for the antiferromagnetic configuration only the hydroxy group results in an effective reduction of the band gap.

Jan 06

## VASP tutor: Structure optimization through Equation-of-State fitting

Materials properties, such as the electronic structure, depend on the atomic structure of a material. For this reason it is important to optimize the atomic structure of the material you are investigating. Generally you want your system to be in the global ground state, which, for some systems, can be very hard to find. This can be due to large barriers between different conformers, making it easy to get stuck in a local minimum. However, a very shallow energy surface will be problematic as well, since optimization algorithms can get stuck wandering these plains forever, hopping between different local minima (Metal-Organic Frameworks (MOFs) and other porous materials like Covalent-Organic Frameworks and Zeolites are nice examples).

VASP, as well as other ab initio software, provides multiple settings and possibilities to perform structure optimization. Let’s give a small overview, which I also present in my general VASP introductory tutorial, in order of increasing workload on the user:

1. Experimental Structure: This the most lazy option, as it entails just taking an experimentally obtained structure and not optimizing it at all. This should be avoided unless you have a very specific reason why you want to use specifically this geometry. (In this regard, Force-Field optimized structures fall into the same category.)
2. Simple VASP Optimization: You can let VASP do the heavy lifting. There are several parameters which help with this task.
1. IBRION = 1 (RMM-DIIS, good close to a minimum), 2 (conjugate gradient, safe for difficult problems, should always work), 3 (damped molecular dynamics, useful if you start from a bad initial guess) The IBRION tag determines how ions are moved during relaxation.
2. ISIF = 2 (Ions only, fixed shape and volume), 4 (Ions and cell shape, fixed volume), 3 (ions, shape and volume relaxed) The ISIF tag determines how the stress tensor is calculated, and which degrees of freedom can change during a relaxation.
3. ENCUT = max(ENMAX)x1.3  To reduce Pulay stresses, it is advised to increase the basis set to 1.3x the default value, which is the largest ENMAX value for the atoms used in your system.
1. Volume Scan (Quick and dirty): For many systems, especially simple systems, the internal coordinates of the ions are often well represented in available structure files. The main parameter which needs optimization is the lattice parameter. This is also often the main change if different functional are used. In a quick and dirty volume scan, one performs a set of static calculations, only the volume of the cell is changed. The shape of the cell and the internal atom coordinates are kept fixed. Fitting a polynomial to the resulting Energy-Volume data can then be used to obtain the optimum volume. This option is mainly useful as an initial guess and should either be followed by option 2, or improved to option 4.
2. Equation of state fitting to fixed volume optimized structures: This approach is the most accurate (and expensive) method. Because you make use of fixed volume optimizations (ISIF = 4), the errors due to Pulay stresses are removed. They are still present for each separate fixed volume calculation, but the equation of state fit will average out the basis-set incompleteness, as long as you take a large enough volume range: 5-10%. Note that the 5-10% volume range is generally true for small systems. In case of porous materials, like MOFs, ±4% can cover a large volume range of over 100 Å3. Below you can see a pseudo-code algorithm for this setup. Note that the relaxation part is split up in several consecutive relaxations. This is done to further reduce basis-set incompleteness errors. Although the cell volume does not change, the shape does, and the original sphere of G-vectors is transformed into an ellipse. At each restart this is corrected to again give a sphere of G-vectors. For many systems the effect may be very small, but this is not always the case, and it can be recognized as jumps in the energy going from one relaxation calculation to the next. The convergence is set the usual way for a relaxation in VASP (EDIFF and EDIFFG parameters) and a threshold in the number of ionic steps should be set as well (5-10 for normal systems is reasonable, while for porous/flexible materials you may prefer a higher value). There exist several possible equations-of-state which can be used for the fit of the E(V) data. The EOSfit option of HIVE-4 implements 3:
1. Birch-Murnaghan third order isothermal equation of state
2. Murnaghan equation of state
3. Rose-Vinet equation of state (very well suited for (flexible) MOFs)

Using the obtained equilibrium volume a final round of fixed volume relaxations should be done to get the fully optimized structure.

 Pseudo-code
For (set of Volumes: equilibrium volume ±5%){
Step 1          : Fixed Volume relaxation
(IBRION = 2, ISIF=4, ENCUT = 1.3x ENMAX, LCHARG=.TRUE., NSW=100)
Step 2→n-1: Second and following fixed Volume relaxation (until a threshold is crossed and the structure is relaxed in fewer than N ionic steps) (IBRION = 2, ISIF=4, ENCUT = 1.3x ENMAX, ICHARG=1, LCHARG=.TRUE., NSW=100)
Step n : Static calculation (IBRION = -1, no ISIF parameter, ICHARG=1, ENCUT = 1.3x ENMAX, ICHARG=1, LCHARG=.TRUE., NSW=0)
}
Fit Volume-Energy to Equation of State.
Fixed volume relaxation at equilibrium volume. (With continuations if too many ionic steps are required.)
Static calculation at equilibrium volume

Top-left: Volume scan of Diamond. Top-right: comparison of volume scan and equation of state fitting to fixed volume optimizations, showing the role of van der Waals interactions. Bottom: Inter-layer binding in graphite for different functionals.

## Some examples

Let us start with a simple and well behaved system: Diamond. This material has a very simple internal structure. As a result, the internal coordinates should not be expected to change with reasonable volume variations. As such, a simple volume scan (option 3), will allow for a good estimate of the equilibrium volume. The obtained bulk modulus is off by about 2% which is very good.

Switching to graphite, makes things a lot more interesting. A simple volume scan gives an equilibrium volume which is a serious overestimation of the experimental volume (which is about 35 Å3), mainly due to the overestimation of the c-axis. The bulk modulus is calculated to be 233 GPa a factor 7 too large. Allowing the structure to relax at fixed volume changes the picture dramatically. The bulk modulus drops by 2 orders of magnitude (now it is about 24x too small) and the equilibrium volume becomes even larger. We are facing a serious problem for this system. The origin lies in the van der Waals interactions. These weak forces are not included in standard DFT, as a result, the distance between the graphene sheets in graphite is gravely overestimated. Luckily several schemes exist to include these van der Waals forces, the Grimme D3 corrections are one of them. Including these the correct behavior of graphite can be predicted using an equation of state fit to fixed volume optimizations.(Note that the energy curve was shifted upward to make the data-point at 41 Å3 coincide with that of the other calculations.) In this case the equilibrium volume is correctly estimated to be about 35 Å3 and the bulk modulus is 28.9 GPa, a mere 15% off from the experimental one, which is near perfect compared to the standard DFT values we had before.

In case of graphite, the simple volume scan approach can also be used for something else. As this approach is well suited to check the behaviour of 1 single internal parameter, we use it to investigate the inter-layer interaction. Keeping the a and b lattice vectors fixed, the c-lattice vector is scanned. Interestingly the LDA functional, which is known to overbind, finds the experimental lattice spacing, while both PBE and HSE06 overestimate it significantly. Introducing D3 corrections for these functionals fixes the problem, and give a stronger binding than LDA.

Comparison of a volume scan and an EOS-fit to fixed volume optimizations for a Metal-Organic Framework with MIL53/47 topology.

We just saw that for simple systems, the simple volume scan can already be too simple. For more complex systems like MOFs, similar problems can be seen. The simple volume scan, as for graphite gives a too sharp potential (with a very large bulk modulus). In addition, internal reordering of the atoms gives rise to very large changes in the energy, and the equilibrium volume can move quite a lot. It even depends on the spin-configuration.

In conclusion: the safest way to get a good equilibrium volume is unfortunately also the most expensive way. By means of an equation of state fit to a set of fixed volume structure optimizations the ground state (experimental) equilibrium volume can be found. As a bonus, the bulk modulus is obtained as well.

Sep 08

## Colloquium on Porous Frameworks: Day 2

On Monday, we had the second day of our colloquium on Porous Frameworks, containing no less than 4 full sessions, covering all types of frameworks. We started the day with the invited presentation of Prof. Dirk De Vos of the KU Leuven, who discussed the breathing behavior in Zr and Ti containing MOFs, including the work on the COK-69 in which I was involved myself. In the MOFs presented, the breathing behavior was shown to originate from the folding of the linkers, in contrast to breathing due to the hinging motion of the chains in MIL-47/53 MOFs.

After the transition metals, things were stepped up even further by Dr. Stefania Tanase who talked about the use of lanthanide ions in MOFs. These lanthanides give rise to coordinated water molecules which appear to be crucial to their luminescence. Prof. Donglin Jiang, of JAIST in Japan, changed the subject to the realm of COFs, consisting of 2D porous sheets which, through Van Der Waals interactions form 3D structures (similar to graphite). The tunability of these materials would make them well suited for photoconductors and photoenergy conversion (i.e. solar cells).

With Prof. Rochus Schmid of the University of Bochum we delved into the nitty-gritty details of developing Force-Fields for MOFs. He noted that such force-fields can provide good first approximations for structure determination of new MOFs, and if structure related terms are missing in the force-field these will pop up as missing phonon-frequencies.

Prof. Monique Van der Veen showed us how non-polar guest molecules can make a MOF polar, while Agnes Szecsenyi bravely tackled the activity in Iron based MIL-53 MOFs from the DFT point of view. The row of 3 TU Delft contributions was closed by the invited presentation of Prof. Jorge Gascon who provided an overview of the work in his group and discussed how the active sites in MOFs can be improved through cooperative effects.

Prof. Jaroslaw Handzlik provided the last invited contribution, with a comparative theoretical study of Cr-adsorption on various silicate based materials (from amorphous silicate to zeolites). The final session was then closed by the presentations of Dr. Katrine Svane (Bath University) who discussed the effect of defects in UiO-66 MOFs in further detail and Marcus Rose presenting his findings on hyper-crosslinked Polymers, a type of COFs with an amorphous structure and a wide distribution in different pore sizes.

This brought us to a happy end of a successful colloquium, which was celebrated with a drink in the city center of Groningen. Tuesday we traveled back home, such that Wednesday Sylvia could start at the third part of the conference-holiday roller coaster by leaving for Saltzburg.

Sep 04

## Colloquium on Porous Frameworks: Day 1

Today the CMD26 conference started in Groningen, and with its kick-off also our own 2-day colloquium on porous frameworks (aka MOFs, COFs and Zeolites) was launched. During the two sessions of the day, the focus mainly went out to the Zeolites, with Prof. Emiel Hensen of the Technical university of Eindhoven introducing us to the subject and discussing how new zeolites could be designed in a more rational way. He showed us how the template used during synthesis plays a crucial role in the final growth and structure. Dr. Nakato explained how alkali-metal nanoclusters can undergo insulator to metal transitions when incorporated in zeolites (it is due to the competition between electron-electron repulsion and electron-phonon coupling), while Dr. De Wijs informed us on how Al T-sites need to be ordered and assigned in zeolites to allow for the prediction of NMR parameters.

After the coffee break Dr. Palcic, from the Rudjer Boskovic Institute in Croatia, taught us about the role of heteroatoms in zeolites. She told us that even though more than 2 million theoretical structures exist, only 231 have officially been recognized as having been synthesized, so there is a lot more work to be done. She also showed that to get stable zeolites with pores larger than 7-8 Angstrom one needs to have 3 and 4-membered rings in the structure, since these lead to more rigid configurations. Unfortunately these rings are themselves less stable, and need to be stabilized by different atoms at the T-sites.

Dr. Vandichel, still blushing from his tight traveling scheme, changed the subject from zeolites to MOFs, in providing new understanding in the role of defects in MOFs on their catalytic performance. Dr. Liu changed the subject even further with the introduction of COFs and showing us how Hydrogen atoms migrate through these materials. Using the wisdom of Bruce Lee :

You must be shapeless, formless, like water. When you pour water in a cup, it becomes the cup. When you pour water in a bottle, it becomes the bottle. When you pour water in a teapot, it becomes the teapot.

he clarified how water behaves inside these porous materials. Our first colloquium day was closed by Ir. Rohling, who took us back to the zeolite scene (although he was comparing the zeolites to enzymes). He discussed how reactivity in zeolites can be tweaked by the confinement of the reacting agents, and how this can be used for molecule identification. More importantly he showed how multiple active site collaborate, making chemical reactions much easier than one would expect from single active site models.

After all was said and done, it was time to relax a little during the conference welcome reception. And now time to prepare for tomorrow, day 2 of our colloquium on porous frameworks.

Jun 30

## I have a Question: about thermal expansion

“I have a question”(ik heb een vraag). This is the name of a Belgian (Flemisch) website aimed at bringing Flemisch scientists and the general public together through scientific or science related questions. The basic idea is rather simple. Someone has a scientific question and poses it on this website, and a scientist will provide an answer. It is an excellent opportunity for the latter to hone his/her own science communication skills (and do some outreach) and for the former to get an good answer to his/her question.

All questions and answers are collected in a searchable database, which currently contains about fifteen thousand questions answered by a (growing) group of nearly one thousand scientists. This is rather impressive for a region of about 6.5 Million people. I recently joined the group of scientists providing answers.

An interesting materials-related question was posed by Denis (my translation of his question and context):

#### What is the relation between the density of a material and its thermal expansion?

I was wondering if there exists a relation between the density of a material and the thermal expansion (at the same temperature)? In general, gasses expand more than solids, so can I extend this to the following: Materials with a small density will expand more because the particles are separated more and thus experience a small cohesive force. If this statement is true, then this would imply that a volume of alcohol should expand more than the same volume of air, which I think is puzzling. Can you explain this to me?

#### Answer (a bit more expanded than the Dutch one):

Unfortunately there exists no simple relation between the density of a material and its thermal expansion coefficient.

Let us first correct something in the example given: the density of alcohol (or ethanol) is 46.07 g/mol (methanol would be 32.04 g/mol) which is significantly more than the density of air which is 28.96 g/mol. So following the suggested assumption, air should expand more. If we look at liquids, it is better to compare ethanol (0.789 g/cm3) to compare water (1 g/cm3) as liquid air (0.87 g/cm3) needs to be cooled below  -196 °C (77K). The thermal expansion coefficients of wtare and ethanol are 207×10-6/°C and 750×10-6/°C, respectively. So in this case, we see that alcohol will expand more than water (at 20°C). Supporting Denis’ statement.

Unfortunately, these are just two simple materials at a very specific temperature for which this statement is true. In reality, there are many interesting aspects complicating life. A few things to keep in mind are:

• A gas (in contrast to a liquid or solid) has no own boundary. So if you do not put it in any type of a container, then it will just keep expanding. The change in volume observed when a gas is heated is due to an increase in pressure (the higher kinetic energy of the gas molecules makes them bounce harder of the walls of your container, which can make a piston move or a balloon grow). In a liquid or a solid on the other hand, the expansion is rather a stretching of the material itself.
• Furthermore, the density does not play a role at all, in case of the expansion of an ideal gas, since p*V=n*R*T. From this it follows that 1 mole of H2 gas, at 20°C and a pressure of 1 atmosphere, has the exact same volume as 1 mole of O2 gas, at 20°C and a pressure of 1 atmosphere, even though the latter has a density which is 16 times higher.
• There are quite a lot of materials which show a negative thermal expansion in a certain temperature region (i.e. they shrink when you increase the temperature). One well-known example is water. The density of liquid water at 0 °C is lower than that of water at 4 °C. This is the reason why there remains some liquid water at the bottom of a pond when it is frozen over.
• There are also materials which show “breathing” behavior (this are reversible volume changes in solids which made the originators of the term think of human breathing: inhaling expands our lungs and chest, while exhaling contracts it again.) One specific class of these materials are breathing Metal-Organic Frameworks (MOFs). Some of these look like wine-racks (see figure here) which can open and close due to temperature variations. These volume variations can be 50% or more! 😯

The way a material expands due to temperature variations is a rather complex combination of different aspects. It depends on how thermal vibrations (or phonons) propagate through the material, but also on the possible presence of phase-transitions. In some materials there are even phase-transitions between solid phases with a different crystal structure. These, just like solid/liquid phase transitions can lead to very sudden jumps in volume during heating or cooling. These different crystal phases can also have very different physical properties. During the middle-ages, tin pest was a large source of worries for organ-builders. At a temperature below 13°C β-tin is more stable α-tin, which is what was used in organ pipes. However, the high activation energy prevents the phase-transformation from α-tin to β-tin to happen too readily. At temperatures of -30 °C and lower this barrier is more easily overcome.This phase-transition gives rise to a volume reduction of 27%. In addition, β-tin is also a brittle material, which easily disintegrates. During the middle ages this lead to the rapid deterioration and collapse of organ-pipes in church organs during strong winters. It is also said to have caused the buttons of the clothing of Napoleon’s troops to disintegrate during his Russian campaign. As a result, the troops’ clothing fell apart during the cold Russian winter, letting many of them freeze to death.

Jun 27

## Folding Phonons

About 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).