OK.
Welcome back, everyone.
Today, we have our second part of our fireside chat on density functional theory
and, in particular, it's the one that I call density functional practice.
Just to remind everyone, what is the plan for this three webinar series.
Yesterday, we had a more conceptual part
on electronic structure methods and density functional theory, in particular.
Today will actually be dedicated to the practice.
What do we need to do this calculation using modern plane-wave pseudopotential code.
And tomorrow it will be, I think, again, conceptual understanding.
What are the properties that we can calculate,
those that we cannot calculate and for those that we can, how accurate we are.
Yesterday, I think we reached very quickly the limit of 500 participants.
So I am very sorry that some of you couldn't attend them, but we have already put the slides
and the talk and also just the audio on the Materials Cloud learn section.
And it's also actually on YouTube.
And we'll do the same today and
later on
for tomorrow's webinar.
Let me actually also monitor the chat.
Yeah. So
we have as yesterday
four panellists, Francisco Ramirez, Norma Rivano, Luca Binci and Davide Grassano.
They will monitor on my behalf the chat.
We don't have the time to answer questions in real time, but they will be collecting if something
very simple they'll answer in real time in the chat.
But by next week, you'll also find all the questions that we deem relevant, we collected 40 yesterday,
and answer in a PDF file, also on the Materials Cloud.
Today webinar is also special in the sense that there is homework
if you want to practice some of this calculation,
there are instruction on the Materials Cloud. . There is a GitHub section
where you can download the handout
with instruction on how to either install Quantum Espresso or just install
Virtual Box and run a virtual machine with Linux environment on whatever machine you have,
either Windows, Macs, Mac or Linux and to go through actual calculations.
And I will focus here on teaching you how to do these calculations.
One thing. One slide that I left out from yesterday.
I was planning to present it tomorrow,
but maybe this is relevant already, is a set of pointers on educational material.
Here is, let's say, a long list.
There are some online resources.
I've just pointed out to some girls, say, that I myself ran a class
on atomistic simulation for many, many years, at MIT actually.
Some of the material here is based on that class
that we kept teaching here in Lausanne or in Berkeley.
And so you can use that there are videos or great lecture notes.
We accumulate materials on the Learn section of the Materials Cloud.
There are plenty of other online resources.
Stefaan Cottenier runs an online stream class on the field and you can just google around.
In terms of books, these are some of my favorite books on fundamentals of quantum mechanics.
The two books by Bransden and Joachain are I think particularly good.
If you want to learn more about the fundamental
of density functional theory the Parr and Yang book is very good
or if you want to learn about advanced
electronic structure methods the two books by Richard Martin and Lucia Reining
and David Ceperley are very good.
And if you have more an eye on materials simulation the books by
Efthim Kaxiras, Jorge Kohanoff and Feliciano Giustino.
And
if you needed just one book to bring with you on a desert island,
it would be either the Parr and Yang if you want to, you know, learn about the fundamentals.
Or it would be the Giustino if you wanted to learn about the practice.
And if it's a really nice desert island I can come and actually teach you on the fly.
Now let me sort of get to where we left it yesterday.
So the Hohenberg and Kohn theorem introduce
the fundamentally exact concept that we can deal with an energy that is a functional of the
charge density and we have even a variational principle.
The only catch is that we don't know what that universal function al of the charge density is.
Kohn-Sham
do a very cunning attempt to identify how we could approximate
that functional breaking apart the universal functional F into three terms
T s - kinetic energy of the non-interacting Kohn-Sham electrons, a Hartree term and then they
push everything that is not known into the third term - the exchange-correlation functional.
Now in doing this they reintroduce orbitals. You see it's not any more
a functional theory of the density.
Now there are just again, single particle orbitals.
It's still exact.
Everything is being compressed back
into the exchange-correlation functional and for the exchange-correlation functional,
we need to find approximations. And, as always, we have a functional so we have also
the corresponding Euler-Lagrange equation that today we'll actually
solve in real life, and so the equations that correspond
to the minimum of that functional, we call them in this context, the Kohn-Sham equation.
They look very familiar.
We have a quantum kinetic energy term.
We have a Hartree potential.
We have an external potential and everything that we had pushed back
into this unknown exchange-correlation energy now becomes in the Euler-Lagrange equation,
a functional derivative with respect to the
charge density and so we call that the exchange-correlation potential.
So this is just exchange-correlation.
And as I said, you know,
Kohn and Sham were both bright and lucky.
They followed Fermi's idea of using the local
density approximation, but not for the kinetic energy: for the exchange-correlation.
So for anything that is
non-homegeneous we construct the exchange-correlation energy
for that system by going point by point and at every point assigning a density
of exchange-correlation energy that is the same density of the homogeneous non-interacting gases.
So we build with all these infinitesimal Lego bricks,
each one at its density, the total exchange-correlation energy.
So, why they were brighter and lucky at the same time?
Well, because the local density approximation
for the exchange-correlation satisfy first and very beautifully a fundamental constraint
on the exchange and correlation that, you know, they didn't know at the time.
But basically and you can imagine that in fairly
intuitive physical terms, what happens if you are an electron sitting somewhere.
Basically, you are creating around you
a lack of probability of finding another electron, because, you know, electrons repel each other.
I mean, if you are, you know,
squeezed in an underground train, it's very unlikely that you will find
someone, you know, 1 mm close to you unless, you know, you want to be there.
So, in general, because electron repel each other, they create around themselves
a hole that is called the exchange-correlation hole.
So first the LDA satisfy this exact -1 constraint, that the exchange-correlation hole
should integrate to -1.
And that's the first thing.
But it's not the only thing.
And actually there are so many beautiful calculations done by the group of Richard Needs
in Cambridge looking at this exchange-correlation hole with accurate stochastic methods.
With diffusion Monte Carlo, and here this exchange-correlation hole
is done for a, you know, number of interacting electrons in a very simple potential.
In this sinusoidal potential that is here.
And, you know,
what we are asking ourselves, if I'm on the ground state and I put myself here,
what is the decrease of probability that surrounds me, this exchange-correlation hole?
And I can calculate it in Quantum Monte Carlo, and that's what you have on the left,
the Monte Carlo hole, and they can calculate it with LDA.
So first of all, LDA integrates to the right quantity, but it also starts to look seriously good.
Now, the catch is that if I now look not at the exchange-correlation hole where I'm sitting
at the top of this potential, but I move around, you'll see a little movie of these holes.
As I move to the bottom and I go back to the top.
On the left, you see the exact results and on the right you see the LDA results.
And as you can see, I mean, it's totally so beautiful that I want to run it twice.
You see here, it's a desperation
so and at the bottom of the potential.
the two exchange-correlation hole looks very different. And actually I won't do it
but through thermodynamic integration, you can see that you can derive the exact
exchange-correlation energy from the exact exchange-correlation hole so the two are related.
But as you saw, LDA here was doing a very good job
but at the bottom of that potential, does a very poor job.
With these it's a saving grace.
Well, it's the fact that actually in the exact expression for the exchange-correlation energy
coming from the exchange-correlation hole, the only thing that matter is the spherical average.
So it doesn't matter the explicit form of its hole, but what matters is the spherical average of this hole
and the spherical average of this hole is much better behaved.
These are again Quantum Monte Carlo results and the LDA does quite a good job.
So the idea from 1965 of Kohn-Sham of using LDA turned out to be very successful.
We have already seen
Marvin Cohen calculations of the phase diagram of silicon in its different phases.
And, you know, these are just examples from 20 years ago from
Chris Pickard looking at the very simple mechanical properties.
In this case is just the equilibrium
lattice parameter across a few different classes of compounds.
I show tomorrow much more extensive examples.
But you see, what we have is that the lattice parameter at the LDA level is actually very good.
I mean, the errors that we do and there is no
experimental input so this is completely from first principle
is of the order of, you know, 1% worse in the lattice parameter.
That would be 3% on the volume.
Elastic constant that, you know, come from the derivatives of the energy.
The second derivative around this lattice parameters, our model will say
starts to be a little bit trickier but we'll see that tomorrow.
But the bottom line here is that the LDA was already incredibly good.
So everything that we have done ever since,
I'll briefly mention it in a moment, is a tiny improvement on top of this.
We might think that we have had major improvements and it was very difficult to get, you know, those
improvements that are very important but LDA was very, very good.
And just historically, you know, of course,
if you have built an exchange-correlation from a local density approximation and you want
to make things better, it might be very, you know, natural to do something, perturbative and say,
let's add maybe the gradients of the charge density to make things better and that was done.
But at the beginning it was actually
decreasing the quality of the results and the reason for that was that actually
the exact integral rule on exchange-correlation all integrated to minus one was not satisfied.
And that's why this concept of generalised
gradient approximation was introduced by John Perdew and others that is
gradient approximation that preserved this fundamental rule and that a lot of the effort
in constructing the new exchange-correlation functional has been driven by the desire
to satisfy some of these constraints and that can be very many.
This is just one example.
But once we had the say correct gradient expansions that satisfy the sum rule,
it was similar that we could actually do quite a lot better,
especially when it was a matter of comparing, say, maybe energies or like here, a question of states
under very different chemical regime or coordination regime.
This is an example of the phase transition from alpha-quartz to stishovite.
Stishovite is a high pressure form of alpha-quartz.
Again, it's found if you look at rocks under meteorites, again, if you are lucky
and you see these two equation of states, you remember how well the transition pressure was
coming for silicon from the diamond phase to beta-tin.
Well, here it turns out that if you use the LDA,
the two equation of states have the minimum at the same energy so you don't really need any
pressure to go from one equation of state to the other.
And that's really in disagreement with experiment
and Don Hamann show, you see, more than 25 years ago that, actually, if you did a GGA calculation of
these two phases, the common tangent was giving
a transition pressure that was in much better agreement with experiments.
And from then onwards has been a proliferation of exchange-correlation functional. There is a
LIBXC library in the solid-state community,
developed by Miguel Marques, that I think has reached almost 500 functionals.
I will not try to describe in any sense all of them.
We'll have tomorrow a better discussion.
You have understood very well what LDA is.
GGAs are exchange-correlation functionals
that introduce this generalized gradient.
The solid-state community, the materials community,
is actually quite fond of the PBE functional from 1996
and if you do a lot of solid-state calculation, a very tiny variation on that,
that add a couple of parameters fit not on molecular energies but on solid state
structures is actually quite a bit better, especially the number of materials like
perovskites, oxides that would come out better in PBEsol.
So this would be my goto functional as a standard.
I will discuss the van der Waals corrections of this functional that are again,
a very successful story of the last, I would say, 10, 20 years.
Tomorrow, I will mention some fancier, more complex, extension.
One that I'm fond of and no-one looks at,
that is the weighted density approximation.
Something that, you know,
when it was introduced was quite cumbersome so it never took off but it has very elegant properties
of the response function so it will be still a very interesting field of study.
To go beyond the gradients, of course,
the Laplacian could be, again, a natural step.
We have to be careful on how to introduce them.
Those were introduced only maybe 15 years ago, and around 5 years ago,
meta-GGA functional that is called SCAN (strongly constrained and appropriately-normed)
has become very popular.
And, you know, it improves in a number of cases and a number of important average onto PBEsol.
There are some subtleties that I'll mention later
on that have to do with pseudopotential and pseudopotential generation, because
as you build those, you want to make sure that your
exchange-correlation for those calculation is consistent with what you do later.
And something that has become very popular,
starting from the chemistry community where not only Hartree-Fock calculations are very
common, but tend to be also very inexpensive due to the use of the Gaussian basis set
are hybrid functionals. In the solid-state community,
the Heyd-Scuseria-Ernzerhof - the HSE functional is very common.
Those mix a certain amount
of exchange and correlation coming from something like
PBE or revised-PBE with a certain amount of exact Fock exchange or often adjust the short-range part
and not to inherit from Hartree-Fock this tendency to make...
to turn metal into insulators.
I'll discuss more of this tomorrow so right now I want just to leave you with
zoology of exchange-correlation functional; the comment that LDA was great to begin with,
PBE and PBEsol make some subtle but significant improvements.
SCAN makes some improvement on top that are still a qualitative feature,
so that we'll discuss tomorrow, that are related to self-interaction.
And that is addressed in part by hybrids, by Hubbard functional.
We discuss them tomorrow . And with the meta-GGA keep in mind also this caveat about
the need, in principle, to use a consistent pseudopotential that have been
generated with the same exchange-correlation functional.
But here, if you want, is all the physical quality of your model,
that is, if your exchange-correlation functional were exact you would have the perfect energy.
So anything that is not a numerical error,
that is what mostly we'll discuss now in the rest of the presentation,
comes from either a physical shortcoming of your exchange-correlation functional or, actually,
from the physical quality of your exchange-correlation functional.
But sadly, there is not a systematic way to make these exchange-correlation
functionals better, although there are some very
interesting and nice ideas based on what is called the adiabatic correction.
And so in some ways here, is all the intuition, you know,
there is the perspective of John Perdew that you have to impose a lot of
known constraints on your exchange-correlation
functional and that is absolutely a very good idea.
You see, there are still qualitative failures
that we'll discuss tomorrow and that are important. Some are intrinsic and some are just
related to our imperfect knowledge of what should be the exact exchange-correlation functional.
OK, with this, we take for a moment a pause
from the concepts and the tale of density functional theory and we delve in,
into this general issue, on the fact that we have now created our Kohn-Sham equations.
These do not look very different from Hartree or Hartree-Fock so the only difference
between these three approaches is having either the LDA or GGA or whatever
exchange-correlation functional here or the Hartree in the exchange term here.
Sorry, the Hartree is always there also in density functional theory or just the Hartree equation.
But in general, we have to solve these differential equations.
In addition, you know, all of these contain the solution so the orbitals that solve
this differential equation, you know, build up the charge density,
that is say, a constituent of the Hartree operator in the Hartree energy,
and so we need from the solutions to build the operators, find the solution,
build the operators and self-consistently converge this cycle.
How do we solve a differential equation?
Well, if this was the 19th century,
you would probably learn what the ideal geometric function is.
You would learn about Laguerre polynomials
and, you know, all the great accomplishment of analysis.
We are not going to do that also because the system we are going to solve are very complex.
So we are going to use a computer.
So we are going to turn a differential equation into something that is computable.
How do we do that?
By actually transforming
our unknown object. You know, in the differential equation we need to find a solution,
we need to find a wavefunction into algebra or linear algebra problem.
That is, rather than finding what is
the unknown psi, we are going to find what are the coefficients.
These are real or complex numbers of this
cm that represent the expansion of the unknown
wavefunction psi into a basis of well-known typically analytical function phi.
So give me ... let me give you an example.
Suppose, you know, our ultimate goal was to solve a differential equation
that was giving me as a solution, something that looked like this.
You know, maybe this looks like a Gaussian
or maybe looks almost like a Gaussian, but not exactly.
Maybe, let's say for simplicity, we are in periodic boundary conditions so here my
periodic boundary drawing skill are quite pathetic, but this thing is periodic.
How can we, you know, transform this differential problem?
Well, we could expand with Fourier analysis this unknown red function psi
into a linear combination of sines and cosines or actually probably because this looks
as if it were symmetric around zero, probably, we just need only cosines.
And here is an example.
You see a...
plotted cosines of different wavelength.
This is the one with the
largest wavelength and then we go to finer and finer wavelength and each of these...
cosines... so sorry, I'm not drawing very well here has a certain coefficient c1, c2, c3,
and so maybe if 10 cosines are good enough to describe psi very well,
what I need to solve my problem is actually finding those 10 coefficients.
And you see, those 10 coefficients need to be such that
around zero they, all the cosines,
sum up constructively and the zone boundary is where -8 is, they sum up destructively.
So you have a maximum, then flat and then a maximum again.
So let me actually show how formally this allows
us to solve the Schrodinger equation in our case, but it could be any other differential equation
or it could be the Kohn-Sham equation or it could be the Hartree-Fock equation
whatever this Hamiltonian operator is.
And I'm using again the Dirac notation so this
ket here in real space would just be a wavefunction psi of r.
And so the core
answers, the core hypothesis that I make is that I expand these psi
as a linear combination of basis elements
and these basis elements I choose them to be orthogonal to make my life easier.
So what I say is that the scalar product with..
between, say, the basis element n and the basis element m
is going to be equal to 1 or 0 depending if n is equal or different from m.
In, you know, our standard first quantization picture
this would be the integral of fm star times fn of r.
And of course we understand right away that if n is equal to m,
this is equal to one from the normalisation condition and if they are different...
well I mean posing that the scalar product is different then you can see that it does say if you
take sine of 3 x and cosine of 3 x, you see right away that they are orthogonal by symmetry.
We don't have to have this orthogonality. We can work with a basis set that is non-orthogonal.
It's just to simplify the algebra here.
So how do I make progress here?
If you for a moment, are familiar also with the concept of Hilbert space
and the psi being a vector in Hilbert space
you understand that H psi is also a vector in Hilbert space
and what I'm going to do is I'm going to project the left-hand side of this equation H psi and
the right-hand side of this equation, a constant times psi onto phi 1, phi 2, phi 3
so if H psi was a vector in cartesian space, it would have three coordinates,
say the vector could be one to seven, and these three coordinates are the projection
of the vector on to the three different versos 100, 010, 001
that constitute the basis of my three-dimensional cartesian space.
Here, my Hilbert space has many dimension, in principle, it could have infinite dimension.
I taking a basis set of orthogonal normalised
versos so there's phi n of the equivalent of the versos of Cartesian space and I project
my Schrodinger or my Hartree-Fock or my Kohn-Sham equation onto each of these verso.
So if for simplicity, I suppose that I have ten elements in my basis set - it could be one
thousand, it could be one million, but let's suppose there are only 10,
I have to project my differential equation onto those 10 versos so m goes from 1 to 10.
So now I don't have any more
vectorial identity as it was here, but here I have scalars - I have projected this.
So you see I have scalars on the left and scalars on the right and I have
right now out of this projection 10 different condition for phi 1, phi 2, phi 3 ... phi 10.
But I have still to exploit...
let me change colour just for the sake of it.
I still have to exploit this ansatz that is I have expanded psi as a linear
combination of my basis set with unknown coefficients.
So I'm going to put for psi this linear sum and so the net result is here.
You see, when I calculate the scalar product
phi m psi this is going to be equal to sum over n c n phi m phi n,
but phi m phi n because of orthogonality, that's the simplification, is going to be equal to
Kronecker's Delta and so this is going to be equal to cm and there is my cn.
I can't do the same trick on the left-hand side,
so you have to keep all my complexity so I have the sum for n that goes from 1 to 10
of these unknown coefficients cn times...
this is, again, a scalar - is the expectation value of the Hamiltonian, could be Schrodinger,
could be Kohn-Sham, could be Hartree-Fock between these two basis elements.
And this is actually a matrix element that we could write as Hmn.
So in shorthand form,
my 10 equations, so remember, m goes from 1 to 10,
my 10 equation are represented here
where the matrix Hmn has 100 elements because m can go from 1 to 10 and n can go from 1 to 10,
and this is the expectation value of phi m H phi n.
So we have 10 equations that are a matrix. The matrix is known, Hmn,
and that's why it's good to have a simple analytical function in those phi m and phi n.
And because when I apply H to phi n,
I want to be able to calculate the Laplacian, the second derivatives.
I want something analytical.
So, if I've chosen my basis set I can calculate
my matrix elements and now I have to solve this 10 matrix..
these 10...
these 10 linear equations that are really written
in a much simpler form, as a matrix vector multiplication.
You see, if I look at the multiplication of the first row
H11 to H1,10
times c1, c10 being equal to E c1. This is exactly the first equation for m equal to one
and then for m equal to two,
i have this equation that is the product of this
row time the same column equal to E c2 and so on and so forth.
So you see how I've taken a differential equation, my only hypothesis was choosing a basis set,
the unknown have become the coefficients of the expansion of my unknown wavefunction
into this basis set, and it's become a linear algebra problem.
This is something that a computer understands.
You have to do matrix times vector equal to a constant times a vector.
I can actually bring...
sorry, let me go back.
I can bring the vector
on the left-hand side so I can make this
equal to zero, and what I have here is that on the diagonal, I have H22 minus epsilon
H33, sorry minus E, down to Hkk, minus E. .
So now I have even a simpler equation in linear algebra.
I have matrix time vector equal to zero.
And you know from linear algebra that the only case where this has a non-trivial solution
for the coefficients - you know, if the 10 coefficients are zero,
OK, so the problem that is not interesting. The only time when I have a non-trivial solution
is when the determinant of the matrix is equal to zero - I wrote the determinant as this.
And, you know, a determinant here, the unknown is E. This is the eigenvalue of my equation.
I know the 10x10 Hamiltonian matrix represented on that basis and this is going to be a polynomial
of tenth degree for this specific case and some degree, in general.
So it has n solutions, let's say 10 solution.
Those 10 solutions will be the 10 eigenvalues
of my problem and then for each eigenvalue I can go back and I can solve...
find the non-trivial set of coefficients c1, ck
that give me the eigenvector that I'm looking for.
So this is how we do differential equation on
a computer and you don't have to be interested in quantum mechanics to know these things.
You know, differential equation are ubiquitous in science and engineering.
If you want to make your plane fly,
you probably have to solve to begin with some differential equations.
So this is a very general thing.
It's not the only approach that we have to solve
the problem, because we could have gone back and thought not at the
Euler-Lagrange equation, that were related to our variational functional
but actually we could have attacked the functional directly.
Either rewritten here for you, the Kohn-Sham energy functional.
You see it's a functional of the charge density.
It's a function of psi i and I'm going now to expand the psi i in this basis set.
So my psi i of r is written as a linear
combination with certain coefficients of my basis set phi n.
So you have the n index that sum over all
the possible basis set and the i index that index the particular electronic states that I want
to consider and if I have 10 electrons, you see I have 10 electronic states.
Well, now, again, I can put this expression for each of my
psi i into this energy functional and I can work out all those integrals
and so the unknown remain my coefficients so that I want to find for every possible psi i.
So I will have here c11, c12 up to c1k and then
depending on how many electrons I have cn1
up to cnk.
So a lot of coefficients that our computers have to deal with,
but finding the solution to my problem can also be stated in this equivalent form in which I have to
minimise the functional, the energy functional with respect to
all my coefficients,
so I can see the problem of getting to the ground state, as a non-linear minimisation problem.
We see examples at the end of today.
For the time being, we actually think that solving the differential
equation as a linear algebra problem, as inversion of that matrix.
Now, we have spoken about basic sets, in general,
but let's actually focus in our actual calculation what would be an ideal basis set.
And if you are a quantum chemist,
often you are dealing with molecules so it's convenient because you don't want to describe
wavefunction in the empty space around the molecule.
You want to describe it only where the molecule is, to have maybe localised basis set.
It could be atomic orbitals.
You might want to decide if you want your basis set to be orthogonal or not orthogonal,
and then you have to keep track of the non-orthogonality in terms of overlap matrix.
It's very common to use a basis set based on Gaussians because then
the Hartree-Fock exchange integrals becomes very inexpensive.
We are going to deal with extended system.
Those could be perfectly crystalline systems like a piece of silicon that can be described very well
by two atoms periodically depleted, but it could be also, say, an amorphous system
that is solid but is not periodic or a liquid system that again, it's extended but not periodic.
And even in all those cases,
suppose we want to study, say, water, liquid water as it exists in my water glass [drinks water]
it's much better rather than to describe a droplet of water with all the water molecules
in and then taking the limit of the water droplet being very large
so that the surface effect that go as the square of the radius become negligible with respect
to the volume effect that go as the cube of the radius.
It's much better to actually get to the limit
of the infinite droplet, not by studying a finite droplet that becomes
bigger and bigger, but to study a system that is periodic.
A periodic system has the beauty that is both infinite and infinite.
So you see, if I put 5 water molecules
in my simulation cell, that maybe it's cubic and I repeat them periodically,
what I have is that my system is infinite
in the sense that there are no surfaces so I don't have to worry about the surface effects.
But it's also finite in the sense that I have only, say, 15 independent atoms...
Sorry, yeah, 15 independent atoms, 45 degrees of freedom in the ionic coordinates.
But what I can do is then rather than study five water molecules, I can study ten water molecules
periodically repeated, and then I can study 20, 30 and so on.
And what happens is that because I've never have surfaces in this,
I get to the limit where the properties of my system don't change any more.
I reach what is called the thermodynamic limit
much, much faster if I were to study the increasing droplet.
So studying, you know, systems in periodic boundary conditions
is very convenient to study anything that is extended. A solid for which is natural,
crystalline solid, but also an amorphous solid glass or even a liquid.
If you study molecules, it's not as good,
but you can study still it very easily by, you know, putting the molecules in a unit cell
that is large enough so that it's periodic image doesn't feel
the effects of the original one.
So if these 2 benzene molecules are far away so that their wavefunctions don't overlap
and they don't have, you know, relevant dipoles or quadrupoles with which they interact
electrostatically, molecule A
will think that it's the only molecule in the universe because it
doesn't feel molecule B, and so even if it's studied in periodic boundary
condition, that means every time we move this hydrogen atom, also this hydrogen atom moves.
It doesn't feel the effect so it thinks that it's alone in the universe.
So we'll use periodic boundary conditions to study
solids and for that, I want to give you a quick refresher of periodicity and the translational
operations in three dimensions because this will be actually the language with which we study
periodic systems like silica, or we set up our calculation even for water or for a glass.
And you are probably familiar with this concept
of a Bravais lattice.
A Bravais lattice is a mathematical concept.
If you want, it's embodied in having
three primitive lattice vectors that we call a1, a2 and a3. That they do not lie in the same plane,
and we look at the infinite set of vectors r that are spanned by the linear combination of
a1, a2 and a3 with the integral coefficients l, m and n. l can be...
They belong to Z, so they can be positive integers, negative integers or zero.
And Bravais,
you know, now we remember his name because he was the first in 1848 to classify
correctly the fact that in three dimensions there are 14 qualitatively different Bravais lattices.
Someone before him thought that there were 15.
Everyone has forgotten his name, poor thing. I will not mention it today.
So these are the 14 Bravais lattices.
I will not go into the details.
We divided them in seven crystal classes, divided them in simple volume, body or face-centered.
Maybe in the case of a cubic system it's worth
remembering we have the simple cubic, we have the body-centred cubic
and we have the face-centred cubic.
And if we wanted to represent, say, an infinite crystal of silicon.
Silicon, the crystal, is a diamond lattice that means it's made by
two silicon atoms repeated periodically with the translation operation of the
Bravais cubic face-centred lattice.
So a real physical object like a silicon crystal is made by a physical group of atoms.
We also call that a basis - in the case of silicon, two atoms
periodically repeated with the mathematical translation of the Bravais lattice.
If you take a piece of ferromagnetic iron
that can be represented by one iron atom periodically repeated
using the volume centred, the BCC cubic lattice.
Sadly, you know, this definition is not unique.
I can represent silicon not by taking, say, two atoms and repeating them periodically
with the FCC operation, but I could take these four atoms that I am painting in green right
now and I could repeat them with the simple cubic operations.
So in some ways I've created a super cell that contains four atoms.
I've changed my translation operation.
The net result is an infinite system that has been described equally well.
OK, so with this, let me also introduce the concept of the reciprocal lattice.
Here we have a direct lattice that is a simple cubic - apologies, the symbols are not coming out
well-described, in general, by the linear combination of a1, a2 and a3.
These solid spheres could be the Bravais lattice points
or they could be, you know, real physical atoms
but what I am going to define is the reciprocal lattice
and for that I'll introduce the concept that may be are familiar with, of a plane-wave.
So,
in real space it's the exponential of the imaginary unit i times a vector G, a wave vector G, times r.
So what is the characteristic of a plane-wave?
Well, it's modulation is going to be along the direction of the G wave vector.
This is because any point that lies in a plane
perpendicular to G will have the same scalar
product between any point - this point, this point, this point and G, so the scalar
product is the same so the value of the amplitude of the plane-wave is the same at every point.
And as we move along this direction,
we are going to have a real and an imaginary part that are both, you know,
cosinusoidal and sinusoidal and the wavelength is going to be very large
if the modulus of G is very small, because if the modulus of G is very small,
I have to move a long way with r before G dot r has changed by amount...
by an amount of pi that would be the phase inversion in the plane-wave.
So G is small, wavelength is long.
G is large, wavelength is short.
And we are going to define the ...
a set of G vectors, such that
the plane-wave i G dot r is, we say, compatible with my primitive direct lattice.
So you see, if I have a plane-wave that looks like this,
this is compatible because it has the same amplitude at every point of the Bravais lattice.
That amplitude doesn't have to be zero.
I mean, it could have been this plane-wave shifted.
It still has the same amplitude at all the point of the Bravais lattice.
But if I were to take for a moment this orange plane-wave here,
this is not compatible because it doesn't have the same amplitude.
So if my Bravais lattice was one-dimensional
you understand, right away, that the plane-waves - I'll do the real part
or the imaginary part sines and cosines that are
compatible are these, you know, something like the fundamental
harmonics, and then what I can have is, you know, two wavelengths and then it becomes difficult.
I could have three wavelengths and any higher wavelength.
So how do I define mathematically this condition
if I want the amplitude to be the same at every point?
What I want is this equality here.
That is, I want my plane-wave to have the same amplitude
at an arbitrary point r and I translate that point by a Bravais lattice vector, capital R.
So if I work out the algebra for a generic
capital R, I want this to be true, and you see it right away,
that this is going to be true if G has the form that is in itself a Bravais lattice form,
as a linear combination of three primitive lattice vectors: b1, b2 and b3,
where b1, b2 and b3 have the property that they are dual to the direct lattice vector.
So bi scalar product aj - I have, you know, three b vectors, I have three a vectors
and the scalar product between each one of the three b and each one of the three j
needs to give me 0 or 1, depending if he is equal or different from j.
This is the property of a dual and you can work
out the algebra that if I have a1, a2 and a3, b1, b2 and b3
are defined by this cross product,
while at the denominator here this is nothing else than the volume of the unit cell.
OK, so this was an introduction to Bravais lattice.
Why is it important? Because we need it to construct
our input file for our simulations, because we might want to study silicon and the simplest
way to study silicon is saying that it is given by an FCC Bravais lattice.
So in the language of Quantum Espresso ibrav is
equal to 2 and then given this translational symmetry, you still need to specify that you have
two atoms that would be at a certain distance from each other.
And in particular, if this is the conventional unit cell, a,
you could say put one atom in 000, you could put the other in 1/4, 1/4, 1/4 a
and then you could repeat all of this
with the FCC lattice translations that are given by these three Bravais lattice vectors.
So that's the importance of Bravais lattices for Quantum Espresso.
In general,
Bravais lattices have also relevance on how we classify electronic states in a system.
You know, remember for a moment that the charged density, say,
of density functional theory is going to be given by the sum over all the occupied states
of the square modulus of each one of these Kohn-Sham electrons.
So in general, what we have often to do is to take sum over all
the occupied states to get physical quantities like the charge density or the total energy.
In a solid the things are tricky because
the system is infinite, so we can do infinite sum and we'll do appropriate
coarse graining of this infinite sum using reciprocal space concepts.
Let me start reminding you, you know,
one of the fundamentals theorem of condensed matter physics, the fact that if the Hamiltonian
commutes with any translation given by a lattice...
Bravais lattice vector R, OK, so if applying the Hamiltonian first
and the lattice translation afterwards or the lattice translation first and the Hamiltonian
later gives you the same result if these two operators commute,
what you can conclude, I will not demonstrate this for you is that the eigenstates of the Hamiltonian
can be expressed in a form where they have also very well-defined
symmetry properties with respect to the translation.
Sadly, and interestingly these symmetry properties
are not that the eigenstates of the Hamiltononian are periodic.
The charge density is periodic,
it's a physical quantity but the eigenstates are not periodic the eigenstates of a periodic
Hamiltonian have the Bloch form, that is, they are written in terms of a periodic function
times, again, a plane-wave where the wave vector of the plane-wave belong to the first Brillouin zone.
And I remind you that if we have, say,
a reciprocal lattice so suppose that we are in two dimensions and we had a simple square lattice in
direct space, we have a simple square lattice in reciprocal space.
These are my G vectors.
This is the origin.
Let's call this G1, G2, G3, G4, G5.
These are the... you see the second nearest neighbours, G6, G7, G8.
So these are the reciprocal lattice vectors.
So what is the Brillouin zone is the Wigner-Seitz cell of the reciprocal lattice.
So what I do is I bisect all the lines that connect a lattice point to these nearest neighbours
and I would have also these bisections and then I take the internal convolution,
and the internal convolution is this.
So this is my Brillouin zone, and so my k belongs to this first Brillouin zone.
So Bloch theorem tells me that at variance with a molecule,
where I have molecular orbitals that I classify with integrals,
in an infinite system where I have infinite number of electrons,
I have also some continuum quantum number, and this continuum quantum number is this k
that is called the Bloch quasi-momentum or sometimes the Bloch momentum
and I classify every state in a solid with an integer number that is the band index n
and with that continuum number in the first Brillouin zone, that is the Bloch quasi-momentum.
And so it is very common both in experiments, you can measure those with arcs, or in simulations to
plot what are the possible energy states for any possible integer number n,
for any possible band index n, any possible crystal momentum or
crystal quasi-momentum, and that is the band structure of a solid.
OK, so we are looking at what are all the possible energy levels for these solids,
and those of you that don't need to take this webinar will recognise here a classic
ferroelectric perovskite with a band gap,
with 2P states here and so in orange I have
plotted the occupied bands, maybe in blue I'm plotting the empty bands but these bands
are the values of energy that are allowed by my periodic Hamiltonian.
Could be the true one, it could be the Hartree-Fock, it could be
the Kohn-Sham, for every possible integer n and for every possible Bloch momentum k.
And Bloch's theorem tells us here,
maybe I'm just using x variable one-dimensional, that my wavefunction psi is not periodic
but is the product of a periodic wavefunction u times a plane-wave and you see the periodic
wavefunction u is this, you know, atomic-like function that I'm plotting here.
So this is periodic. So this is unk
and it will vary as a function
of n, as a function of k, but it is periodic and the product of unk in, let's make it, blue
times e to the i k x that is in red is going to give me
the wavefunction psi that you see is this black wavefunction.
It's not periodic but is the product of something periodic times our red modulation.
OK, I think to give everyone some breathing space,
we'll stop for five minutes, so we'll start again at 07.
So I'll stop for five minutes the video and you can all take a breath.
OK, let's... let's go back.
So we have looked at how the wavefunction
looks like and as I move around, say if this is the zone centre
and this is the zone boundary,
what I go is from a modulation that has a very long wavelength
to a modulation that starts to change phase at ... for every time I cross a lattice vector.
So,
the important message here is that Bloch theorem has told me what is the shape of psi.
It's not periodic, but is the product of a periodic part times a phase modulation,
and so it's this periodic part, unk, this blue part, that we need to know
for every band that is for every index n and for every k.
So for every continuum momentum in the Brillouin zone is this
blue periodic term that is going to be expanded in plane-waves.
Because this is periodic the plane-waves,
and that was the important concept of reciprocals space,
the G are chosen such that E to the i G dot r has the same periodicity or higher periodicities
compatible with my direct lattice, and so I can expand unk
as a linear combination of plane-waves.
This is what your code does.
It approximates the periodic part of the Bloch function
as a linear coeffici... combination of plane-waves.
The better you want it to approximate, the more plane-waves you need.
If your unk in your material happens to have a lot of structure, you are not going to be able
to simulate it with just a few plane-waves because you need to go to plane-waves that have
very high periodicity to be able to describe all that fine structures,
and for each of those plane-waves you'll need to know what the coefficient c is.
The coefficient c of the plane-wave G for that band n and for that k-point, k.
So in a typical plane-wave calculation,
you might have a basis set that is of the order of thousands, tens of thousands of plane-waves,
you can have bands that are 10, 20, 100, 1000.
The k can be a lot.
We'll discuss the k in a moment.
So a lot, a lot of coefficients.
So that's why you need computers.
What is good and what is bad about plane-wave
basis set, exactly as in Hi-Fi it's very easy to improve the resolution.
OK, so when in 1981 and 1982, Philips came out with the CD-ROM,
what they were doing is they were sampling a sound so that you could have basically
all frequencies up to around 22 kHz.
Relying on the fact that most rich people that were buying CDs at the time
wouldn't have the sensitivity of, you know, a twelve year old that can go to 25 kHz.
Also, it was good because then you can...
could fit the Ninth Symphony of Beethoven in just one CD.
And that was a major point.
I guess, you know, it's very important over there in Northern Europe.
But so by taking more and more plane-waves, you are able to describe with finer and finer
resolution your wavefunction so you can make it systematically better and better.
You need a lot of them and if you were to, you know,
need to describe, say, in a piece of gold, both the valence electron of golds
that are very smooth, but maybe the 1s electron of gold that are
extremely localised around the nuclei, you have an untenable number of basis elements
and what you will discover, in a moment, is that we are not going to
study explicitly the core electrons in a solid that are very sharply varying in space,
but are also not very active in forming and creating chemical bonds.
We are only studying the valence or so-called
semi-valence electrons and luckily, those are very small.
So
I'll show you an example in a moment that
it's very easy to calculate the gradient of the Laplacian of our plane-wave expression.
I'll comment later on that it's very easy to calculate...
because it's very easy to calculate the kinet ...
the Laplacian is very easy to calculate the kinetic energy operator.
So remember, if we simplify for a moment and we think our Hamiltonian as being a kinetic energies
plus a potential and the kinetic energy is just minus 1/2 Del square...
. the matrix element
of the kinetic energy between two plane-waves is super easy to calculate because basically
if I have a plane-wave e to the i G dot r,
its gradient is equal to i G times e to the i G dot r,
so is here there are derivative automatically and its Laplacian
is going to be equal to minus G square e to the i G dot r.
So twice, you know, I have taken the, if you want, divergence of the gradient.
So very, very easy to deal with this
and it's also very, very easy to go from real to reciprocal space.
So I have a quantity that is written, say, as in this case, the charge density that is
periodic as a linear expansion in plane-waves with coefficient rho of G.
So either I know the coefficients rho of G,
the coefficient of the plane-waves,
and so I have the reciprocal space expression of the rho G but using a fast Fourier transform,
it's very, very fast to go from the rho G
to the rho of r, where r belong to a direct space mesh,
say this could be the Wigner-Seitz cell of my crystal, belong to an equispaced
Wigner-Seitz cell that is nothing else than the dual of my set of G vectors.
So fast Fourier transform, allow us if we have
the rho sampled in real space on a grid to give us the coefficients
rho of G very efficiently.
The cost is just a linear times a logarithm of the linear cost, so it's just above linear.
Or if I have the rho G to go back to the rho
on a mesh, so I can go back from real and reciprocal space very, very easily.
And, you know, one of the fundamental concepts that you'll use
when you do a Quantum Espresso or any kind, you know, it could be a VASP, it could be
a CASTEP, it could be an abinit calculation, "How do you select what is your plane-wave basis set?"
Well, suppose this is my reciprocal space,
so these are my G vectors.
This is gamma.
This is the set of G vectors.
So to every one of these generic G vectors,
this is G1, this is G2, I can give them names, there corresponds a plane-wave.
The larger the G is, so this one has
a periodicity that is shorter than this one, and this one has a periodicity that is longer...
Sorry. The opposite.
So
if I call this A, B, C and D, A has a periodicity that is longer than B, that is...
that has two wavelengths in the space
that A has one wavelength and C has three times, three wavelengths,
and D has four wavelengths and so as I go
to larger and larger G vectors I have plane-waves that become more and more
detailed and so they express better and better
the fine features of my wavefunction.
So how do I choose what is my basis set?
I choose the radius of a sphere centred around the origin and if I give you this radius,
what I am saying is you are going to use all the G vectors that follow inside of this sphere.
And then I can choose a larger radius.
and what I have is a basis set that contains
all the previous basis set elements and more so it has become systematically better.
So choosing a cutoff radius is a very economical and efficient way to express
what is the basis set that I'm going to use and, actually, rather than using the cutoff radius,
we give the energy that corresponds to that radius.
Because if let's say this corresponds to a vector GMAX,
the kinetic energy of that GMAX vector...
remember the Laplacian of e to the i GMAX...
is going to be minus G square, and so this kinetic energy is equal to 1/2 G square.
So this quantity is called, in Quantum Espresso, ECUT.
It could be the wavefunction cutoff or it could be the charge density cutoff.
Let me make a little detour.
If I have that, say, a wavefunction is written as u times e to the i k dot r,
I dropped the index n just for simplicity,
and if I look at the square modulus of psi k, that is the charge density.
So the charge density rho k for these orbitals is equal to the square modulus of psi k
and is equal to the square modulus of uk, because when I take the square modulus,
the i, ik, r modulation goes away.
Well, you see, if my uk was written as an expansion in plane-waves with,
you know, certain G vectors of the coefficient, c, I call them kG e to the i G dot r
where G is a modulus that is less than the G MAX cutoff.
But now if I'm looking not at uk, but I'm looking at the square modulus of uk.
This is uk star times uk, so what I have here is sum over G prime ck G prime star
e to the minus i G prime dot r.
This is uk star times uk, that is sum over G ckG e to the i G dot r.
And lo and behold you discover,
I'll write it here that the charge density, that is the square modulus of psi,
that is the square modulus of u is written also as a plane-wave expansion that is given by
ck star G prime ckG e to the i G minus G prime dot r.
So what do we discover?
You see, we discover that if the cutoff
for the wavefunction was this, that is my G vectors were
inside this cutoff sphere, the charge density is described by
plane-waves that can have a wave vector that span all possible combination
of G minus G prime and of course, G minus G prime, if G and G prime live inside
this blue sphere, so G minus G prime are going
to live in a sphere that has twice the radius,
because if I take G here and I take G prime here, G minus G prime is going to be here.
So the cutoff radius for the wavefunction
is going to give me naturally a cutoff radius for the charge density that is twice as large.
And because we are not talking about radii,
but we are talking with energies associated with radii, it tells me that the energy cutoff
for the wavefunction is related to the energy
cutoff for the charge density by a factor of four.
So if I decided that the energy cutoff for my wave function is 30 Rydberg,
the energy cutoff for my charge density needs to be 120 Rydberg.
So that was one statement.
The second statement that is also important to know - let's do this in black is that,
as I said, Laplacian are easy, and I'll give you an example of how am I going
to solve the Poisson equation using a plane-wave expansion.
So Poisson equation is the Laplacian of V is equal
to 4 pi of rho and suppose that I know that rho is
given by a plane-wave expansion where the rho G coefficient I know.
So,
I know what these are,
my plane-wave expansion.
And V, I don't know, I want to solve this equation and I write V also
as a plane-wave expansion, but the VG are not known.
OK, so I want to solve this differential equation and not in real space, but in reciprocal space.
Given all the set of rho G s what is all the set of VGs?
And that is, you see, super easy because I
insert these equations here and what I get is sum over G VG times minus G square e to the i G dot r,
is equal to 4 pi sum over G rho G e to the i G dot r.
So plane-waves happen to be also orthogonal,
and orthonormal, I won't go into the nomalisation of periodic function here,
but so if I take the integral of this expression times the complex conjugate
of any plane-wave I will have that because of orthogonality
that contribution is zero, and instead if
my scalar product is with that plane-wave that has wave vector G, that is not going to be zero.
And so I get right away the solution that is coefficient by coefficient.
These plane-waves need to have the same coefficient and so if I have the
charge density coefficients,
the coefficients of the potential... it could be, you see,
the Hartree potential that solves the Poisson equation... are just given plane-wave by plane-wave
by the coefficient of the charge density divided by G square.
Now, what happens when you have a very, very large system, a very, very large supercell.
Remember that reciprocal space goes as the inverse of direct space.
So as I make my unit cell larger and larger,
my G vector of my reciprocal lattice becomes, you know, made by vectors that are closer and closer.
If you remember, the expression of the dual said
the volume was at the denominator, so if the direct lattice vectors goes up linearly,
the reciprocal lattice vector, decreases inversely proportional to that.
So in a big system, I will have a lot of G vectors that are very small and if the G can be very
small, it means that a little error of the charge density is going to give me
a large error, because it's divided by G square, on the equivalent component for the potential.
And that's why electrostatics of large system is very tricky.
Things that you're studying by first principle:
a metal slab with vacuum or maybe with some molecule looking here and smiling at you.
Well, and you're solving this in reciprocal space
so you're taking this charge density, you're Fourier transforming,
you're getting the rho G, and then out of that, you're getting the VG are small.
Error on rho G gives you a big error on the VG and so it means that all of a sudden you could put
the charge density here, and this is a large system so in the next
iteration, there will be an electrostatic imbalance,
that was to send back the charge density here and this divergence in the infinite supercell
limit gives rise to what are called sloshing instabilities that make the convergence,
the iterative convergence, of first principle calculation slightly delicate.
OK, so, as I said, we use plane-waves but we didn't have to do that.
A lot of quantum chemistry codes
use Gaussian basis set because the exchange integrals are easy to do.
There are codes like CRYSTAL that do
solids with Gaussian basis set.
There are simplified codes that use linear combination of atomic orbitals.
There are codes that use a linear muffin-tin
approach that has been introduced by O. K. Andersen in the seventies.
I don't know if you know what a muffin tin is.
Muffin tins.
Muffins have become very popular ever since cupcakes have become very popular.
So a muffin tin is something to bake cupcakes in the oven.
And so a muffin-tin approach is an approach where seen from the above
we partition space in atomic spheres around the atoms
and we use one set of basis sets inside and in the interstitial valence region
we use plane-waves and we have to match the basis set at the boundaries very well.
And that, in particular, LMTO was very powerful when it was first
introduced because it could allow to do calculations that before were impossible.
These days we even do better.
We use linearised method
plane-wave techniques that again use atomic like basis sets inside the sphere
and plane-waves in the middle and allow ultimate accuracy.
So, often there are a lot of reference calculation
or calculation on very delicate properties that can be done rather than using plane-waves and
pseudopotential, using all-electron technique that would be impossible to deal with,
with plane-waves, but with LAPW you can use atomic orbitals around
the nucleus and plane-waves in the middle so you can study the entire system.
But again, Quantum Espresso like VASP, like abinit, like CASTEP, like many other codes,
we use plane-waves and that's how we are going to use.
The other technical concept that you need to be
familiar to do an ab initio calculation is the one of Brilouin zone integration.
As I said before, there are a lot of quantities
that require summing over all the states in the system.
Say, let's take silicon, an insulator.
My charge density would be given by the sum over
all the bands and by integral over all the k
in the first Brillouin zone of the square modulus of the unk.
So we can't do continuous integrals in our
calculations, so what we are going to do is we are going to do sampling that means we take
the Brillouin zone and we sample it on a finite discrete mesh of k-points.
So we don't consider all the blue states,
but we consider all the black states and by increasing the fineness
of our k-point sampling mesh, we do better and better.
And the general technique is to use equispaced meshes that are called
Monkhorst-Pack meshes - I'll show them in a moment - so we are doing something, actually,
quite trivial. We are sampling there the Brillouin zone on an equispaced mesh.
There are a couple of subtleties.
One that I'll describe later.
If you have a metallic system you have a Fermi energy that separates the empty states
from the filled states.
And if you want to calculate the charge density, say, in copper, you have only to sum over
the occupied states and so what you have in a metal is actually a discontinuity.
You might be, you know,
following a band and integrate things but as soon as that band crosses the Fermi energy
it doesn't contribute any more to the... to the charge density of the system.
You know, the simplest case that could be a case of,
you know, a very simple, model one-dimensional system that there's a Fermi energy here.
This is the band structure.
So what I have when I do the integral,
I need to count the states but as soon as the states cross the Fermi energy,
their occupation goes to zero, so you have to integrate a discontinuous function.
Discontinuous function are much nastier to integrate carefully rather than smooth functions.
So insulators are easy to sample in k-space than metals and we'll get to metals in a moment.
The other limit that is interesting is the limit
when you study a very large system and maybe the simplest way to think of this is, again,
taking this band structure here in for a one- dimensional system that was looking like this.
You know, this could be, you know,
the band structure of a simple monoatomic chain, and as you know, I'm in real space.
I could study this in the primitive unit cell,
maybe one atom per unit cell, but no one stops me from studying it
in a double unit cell and a triple unit cell and a quadruple unit cell.
And in reciprocal space, let's say if I double the unit cell,
that would correspond to actually having the same system but described in a different way.
Rather than having this Brillouin zone, now you have a Brillouin zone that is one half
and rather than having one band, in black, I have, let's say two bands, in orange.
So the black band and the two orange bands are describing exactly the same system.
But in a different representation,
the black one comes with having in direct space, one atom per unit cell the orange comes
from having two atoms per unit cell. And as I increase my description in real space,
what I get is that my bands get folded and folded and folded.
So, you know, I might end up with having something that looks like this.
So you see, if I study this system in a very,
very big supercell, I have a very small Brillouin zone and I don't have to do a lot of effort
in sampling it because I just sample it at one point, because when I sample it at one point,
I get all these point here in the middle that corresponds to taking all these point here.
So if your unit cell is very large, you need...
you can get away with very coarse sampling.
Maybe just one point.
It could be the gamma point.
It could be the zone boundary.
Or you could ask yourself the question, very smart
of, you know, if you can only afford one point, what is the best point that you can choose?
And it might not be the origin.
It might be something else.
And that's called the Baldereschi point,
and if you're actually dealing with an FCC lattice, that's 1/4 1/4 2 pi over a.
This was especially relevant when, you know,
calculation was so expensive that you could, you know, use...
you could spare any CPU that you...
that you could, so you wanted to sample those as
little as possible, but it's still very relevant these days.
If you are studying your system in a big super
cell, then maybe you can get away with only one point, but then maybe the Baldereschi point is
better than the gamma point because it's actually more accurate.
The gamma point has another advantage,
that is at gamma point the wavefunction are real, so those coefficients don't require
a complex expansion, but you can just get away with half of the informational
content and Quantum Espresso in the gamma only formulation takes advantage, actually, of this.
But in general, this will be what we will do.
We will use equispaced meshes.
Those are called Monkhorst and Pack from the paper in the mid 70s that describe those.
Quantum Espresso generates automatically for you, so you don't have to worry very much.
But just so that you know what they are,
what I have on the right, the (4,4,4) unshifted mesh, you see, is just
sampling the Brillouin zone with
four points in every direction and, of course,
this point here with the cross is equivalent by symmetry with this point here.
And you can have the unshifted mesh or you could
have with the same fineness of sampling what is called the shifted mesh
where this k-point sampling are here in the middle of this division.
And Quantum Espresso would take these
as a 4 4 4 1 1 1, 4 4 4 0 0 0.
One thing that is very important.
Suppose for a moment, we are studying, and this is a real space cell,
a carbon nanotube that, you know, we might describe by something like this
In reciprocal space, you'll have also a Brillouin zone
but what happens in your system is that this y direction is the direction of periodicity
and so is the direction in reciprocal space where there will be band dispersion.
In reciprocal space, in the x and z direction, there will be no dispersion because the system,
the nanotube, is effectively isolated and there's not talking with these periodic images.
So the states don't have a k dependence in the expectation value for the energy.
So the Monkhorst-Pack mesh that you would use in this case is not maybe a sequence of better
and better Monkhort-Pack meshes, be 2 2 2, 3 3 3, 4 4 4
until your quantities are converged.
You don't need to sample in two directions so they will be 1 2 1, 1 3 1, 1 4 1
and, you know, 1 4 1 is going to cost 16 times less than 4 4 4.
You have 64 point in one case and only 4 points in the other.
So that is also very important to remember.
If your system is not truly three-dimensional, but is separated by vacuum.
In this case, two direction, as is the case
for a carbon nanotube or in three direction, if this was a molecule.
For a molecule, you always only use gamma sampling, you know.
For a molecule, there is no point in sampling
anything because the molecular states don't have any k-dispersion.
If they have a k-dispersion,
it means that you don't have a molecule, but the molecule is feeling its periodic image,
so you are actually having something that is an organic crystal.
And if you have, say, a layer of graphene,
you sample in the plane, but you don't sample perpendicularly.
The last thing that makes, actually, our calculations less expensive in solids
is exploiting symmetry operation, the point group symmetry of the system.
You might have mirror planes, you might have a rotation axis.
If you are studying, you know,
a silicon crystal you actually have 48 symmetry operation on the Td point group.
And symmetry is used very effectively in doing these sums,
because rather than, say, in this case, if this was my Brillouin zone,
rather than summing over all the k that belong to the Brillouin zone,
I can use, I will not demonstrate it, but I can use this relation that tells me
that if I know the orbital at a certain k,
I can obtain the Bloch orbital at a certain...
sorry, this is unreadable...
at a certain
S of k that is at a k that is related
to the previous one by a symmetry operation, simply counter rotating my points in real space.
That is to say that if this point here is k
and this point here is S of k, and they are related by a minor operation, if I have psi k,
the psi Sk
is obtained by exactly folding back
the points according to this operation.
So my sum is broken apart.
The sum over all the k-points it becomes a sum over only the k points
that live in what is called the irreducible wedge of this segment of the Brillouin zone that once we
apply all possible symmetry operation blankets the entire Brillouin zone.
Here I have only one mirror symmetry and so you see the orange wedge plus the middle symmetry is
going to cover also the other half of the Brillouin zone.
So here I do the sum over the irreducible wedge and all other points that are outside
the irreducible wedge are not obtained by doing additional calculation, by simply...
but by simply, using this symmetry relation here.
So you see, at the end of the day,
I have to solve the Kohn-Sham equation only at the irreducible points of the Brillouin Zone.
And, as I said, if I'm studying silicon is a saving of 48.
Of course, you know, a factor of 48 is not, you know, life changing these days so
sometimes we don't bother with this, but in very complex codes that use a lot of
calculation for, say, calculating phonons and electrons
and electron-phonon interaction, you know, these savings build up very quickly.
While symmetry, that is always there if we don't have a magnetic field,
so if we don't have a vector potential in our Hamiltonian, is time reversal symmetry.
So what I written here is the Hamiltonian,
not any more for the Bloch function psi.
but what happens if I apply H to psi,
that is if I apply H to unk, e to the i k r
so if I put explicitly the Bloch form it turns out that the Hamiltonian for unk has this form here,
and you see right away that if I solve this Hamiltonian for unk the solution for un minus k
is immediately given by taking the complex conjugate of unk.
So always if this is my Brillouin zone...
OK, this is the gamma point...
I never need more than half of the Brillouin zone because the points in the other half...
OK, if, you know,
the points on the left here have, say, negative k, the points on the right have positive k,
I get from the one to the other, that is I go from k to minus k, just taking the complex conjugate.
So you never have to do a full sampling
of the Brillouin zone but worst case scenario you just have to sample half of it.
PWSCF, that is the self-consistent code
in Quantum Espresso, is very sophisticated about k-points.
It gives you all these option.
It gives you the option of just choosing
the gamma point for a big super cell and then your calculation becomes
twice as fast and twice as
less memory intensive.
Most of the time you lose automatic k-point meshes,
maybe unshifted with k1, k2, k3
0 0 0 shifted 1 1 1
but you need to test when your calculations are converged with respect to nk1, nk2, nk3.
And the calculation goes linear in the product of these three numbers.
So doing a 2 2 2 mesh is 8 times
cheaper than doing a 4 4 mesh, that is in itself 8 times cheaper than doing
an 8 8 8 mesh, that is in itself 64 times more expensive than a 2 2 2 mesh.
So you want to, you know, proceed carefully but you need to check
that your physical quantities are converged with respect to the k-point sampling.
Remember again, you know, if you have a molecule you just use gamma.
If you have a nanotube, you sample only one direction or you have a wire, a nanowire.
If you have a two-dimensional system or a metal slab with a molecule absorbed, you don't need...
you don't want to sample in the direction perpendicular to this.
OK, so this was k-point sampling, so we are starting to look at all
the computational tools that we need to run a proper calculation.
We have looked at the plane-wave cutoff, the basis set, the k-point sampling.
Let me introduce
in very handwaving form the concept of pseudopotentials, that is the idea that,
actually, when you look at a solid or you look
at a molecule, the electrons deep down in the core,
say, for aluminum, it would be the 1s, the 2s and the 2p electrons,
they are so tightly bound to the nucleus that they don't care if aluminum is in a molecule,
is in a bulk solid, is at the surface or is in molten aluminum.
And so for this reason, what we want to do is only
study the valence electrons and the core electrons that are very low
in energy, we consider frozen in the atomic configuration.
Of course, if you want to study, say, X-ray spectroscopy,
is where you excite electrons from the core, well then you need to do something else.
We will not go into that but that are also
something that you can calculate with standard density functional codes.
So how do we get rid of the need to study core electrons?
You see, core electrons have very large negative energies,
so a small error on a core electron messes up our calculation because, you know,
binding energies are of the order of a fraction of an electron volt.
So doing, you know, a .1% error in the 1s state completely messes up our calculation.
So we don't want to deal with this and the core electrons are very localised around
the nucleus, so they would need a lot of plane-waves.
So to do that, we introduce the concept of the pseudopotential.
What is a pseudopotential?
Well, say for the case of aluminum,
a pseudopotential is a potential that we create that reproduce the effects,
not only of the nucleus, minus 13 over r, the coulombic attraction,
but it reproduces the effect of the nucleus screened by the frozen core electrons.
And so we want to create that pseudopotential
so that when we throw the three valence electron of aluminum in the field of those pseudopotential,
the energy of those three electrons will be these energies here,
-7.8 and -2.7, that are exactly the energy that the true...
that, sorry, that they would have had if they had
felt the full explicit electronic core configuration.
So we build pseudopotentials so that
valence electron thrown in that pseudopotential field have the same energies that
valence electrons thrown in the through field of the nucleus and all the core states would have.
But not only that, we want also that the valence wavefunction is correct
and we are going to care for the valence wavefunctions to be correct outside of the core.
So we want the valence wavefunction to be correct
outside of the core because that's where chemical bonding takes place.
And so here you see what I have as a thin line is
a coulombic attractive potential if, say, I have an atom sitting here at the origin.
In blue, sorry, in solid red, in thick red, I have a pseudopotential that agrees exactly
with the true coulombic potential outside the core and is different inside the core,
and this thick red pseudopotential has been built such that the electron in the field
of the pseudopotential obtains a ground state that is the thick blue curve,
and that thick blue curve has two characteristics.
One is identical to the thin blue curve that was, you know, this could be the
6s, 7s state in an atom.
So outside the core, the pseudo wavefunction and the true wavefunction are identical
and also the energy of the blue state, thick, is equal to the energy of the thin blue state.
How to do that?
How to build this pseudopotential is very complex,
but you gain a lot because all of a sudden you don't have to study all those core electrons
that are very high energy, so very little error are going to be very harmful
and they are very sharply peaked, so they require an impossibly large cutoff.
How do I build these pseudopotentials?
There were a number of techniques developed, but really at the end of the 70s Don Hamann,
together with Michael Schluter and Chiang, developed these ideas of norm-conserving
pseudopotentials and they developed a construction protocol for pseudopotential,
where as I said, the eigenvalues, and they agree for the real states
and for the pseudo states, the wavefunctions agree
beyond core radius.
The concept of norm-conservation told us that the integral of the pseudo wavefunction
and the integral of the pseudo and the real wavefunction both squared, that is the real
and the pseudo charge are going to be identical within the core, the integral,
and then because the pseudo and the
real wavefunction agree outside the core, they're also going to agree outside the core.
And then there is another subtlety
that I'll just mention in passing, there is going to be an agreement for,
and I'll explain what this is in a moment, what the logarithmic derivatives are,
and there are first energy derivatives.
I'll get to that in a moment.
But in order to make this norm-conserving pseudopotential,
it turns out that what you have to do is,
you know, you can't really do it with just a local pseudopotential, that is a pseudopotential
that depends only on position.
You need to act differently on the wavefunction, that is, you need to act differently
on the wavefunction that looks like an s orbital, or on the wavefunction that looks like a p orbital
or a wavefunction that looks like a d orbital.
That is to say, in order to act, your pseudopotential
takes a scalar product on the spherical harmonics and you are going to act differently...
so this would be the matrix element...
you are going to act differently on the component
of the wavefunction that has a non-zero projection onto an s spherical harmonic, onto a p spherical harmonic
onto a d spherical harmonic and so on and so forth.
So, proper norm-conserving pseudopotential will act
and it's, you know,
represented here differently on the s component of the wavefunction, on the p component,
on the d component, on the f component, and at a certain point for higher order
angular momentum components, it will stop doing that because it's not necessary anymore.
So maybe for the case of gallium,
we have an s component, a p component, a d component,
and we say for anything that has a higher order, angular symmetry will just act in the same way
and that will be the local component of the pseudopotential that might look like this.
And again, you see outside a certain core radius, they all look identical.
And if I study in this case, the gallium atom,
you see, I'm going to reproduce the s, the p and the d wavefunction of the atom
exactly outside the core.
So this is inside the core.
This is outside of the core.
And inside the core, they are different.
They are much smoother,
so you see they are much easier.
You see how sharp the s wavefunction was here.
So the 3s,
so they are much smoother but they have still the same energy.
And the last concept was this concept of a logarithmic derivatives.
I won't go into the detail of this
but it comes from scattering theory. The way at the end of the pseudopotential are built
comes from the idea that a real core, that is the nucleus plus the real core electrons
and a pseudopotential, that is a potential that represents all of this in a, you know, a unique potential.
So there are no core electron explicit in the pseudopotential,
so these two object, the real core
with the nucleus and the pseudopotential should scatter an incoming plane-wave in similar ways.
And what does it mean in similar ways?
It means that you can take the plane-wave and after the scattering you can do a Born
expansion in spherical waves and you can represent the effect of the scattering
in terms of the phase shifts of all these coefficients.
And the norm-conserving pseudopotential are built to have, say, the correct phase shift,
that's the logarithmic derivative at one given energy,
but they are also transferable.
That means that they have the correct derivative and having not only the absolute value,
but also the correct derivative means that those phase shifts are going to be actually good
also in the neighborhood and, actually, you see,
you can make this phase shift to be good in a very large neighbourhood of energy.
So that's why the norm-conservation concepts
of Hamann, Schluter and Chiang was so powerful.
Later on in 1985, David Vanderbilt extended this concept and it developed a formalism
by which you could impose this condition and not at only one value of energies,
but at multiple values.
So you could reproduce better and better
the scattering from the real core and the pseudopotential in a wider region in energy.
And these are the extended norm-conserving
pseudopotential that are now becoming more and more used because
you can make them more accurate than standard norm-conserving because they have better
scattering properties but are not so complicated as
what had been developed in between, also by David Vanderbilt,
of the ultrasoft pseudopotentials or by Peter Blochl, the projected augmented wave method.
These two approaches are somehow similar.
They tend to make the calculation much less
expensive, but they also tend to make the entire algorithm more expensive.
I'm looking at my notes, I think, today I'm going long.
So I'll mention I think we'll have 20 more minutes to go.
If I don't faint. If I faint, that's it,
if I don't faint, you know, and you also know that all of this is recorded and will be posted.
So if you have to go, you can also watch it, watch it later.
One thing that is interesting about the ultrasoft pseudopotential, I just draw a graph here.
Suppose that this was the charge density, the radial charge density of a 3d orbital,
well, the ultrasoft and similarly the projected
augmented wave, does very smart, introduces some very smart formalism
to make the wavefunction, the pseudo wavefunctions that we have in our calculation,
even smoother than what norm- conserving pseudopotential can do.
And when they were introduced in the early 90s,
it means that materials that were very challenging for plane-wave calculation, materials with 2p electrons
or 3d electrons... you see 2p electrons or 3d electrons
don't have to be orthogonal to 1p
because there is no 1p, or to 2d because there is no 2d.
So 2p and 3d electrons are very compressed towards the nucleus
so they need a lot of plane-waves.
They are very tough to describe and ultrasoft
pseudopotential manage to describe those with very, very smooth pseudo wavefunction
but in that, what you lost is that the charge density was not any more the square of the
wavefunction, but was actually the square of the wavefunction
plus some so-called augmentation terms because the charge density is a physical observable,
and so if I make my wavefunction smoother than it should be,
I'm not going to have any more that the square
of the wavefunction gives me the true charge density.
I need to have these augmentation terms.
And so because of this,
it means that the relation between the cutoff for the charge density, being 4 times the cutoff
for the wavefunction, has become broken. The cutoff for the charge density cannot change
because of any computational trick that I do,
because that's a physical quantity and that depends on the system I'm studying.
With the ultrasoft pseudopotential, I can make the cutoff on the wavefunction
smaller and so it means that the ratio between these two
is not any more 4, but it could be 6, 7, 8, 10, 12.
It's something that, in principle, you also have to test.
Luckily, if you go to the Materials Cloud,
we have tested for you a standard set of pseudopotential.
We call them standard solid-state pseudopotential, SSSP, and we have prepared them for the case that you want to
use PBE as your exchange-correlation potential or PBEsol as your exchange-correlation potential.
No one has generated, as far as I know, a set of pseudopotential for the SCAN functional.
That would be very good.
There are some, but not an entire periodic table.
So there are going to be inconsistency when you do SCAN calculation between your pseudopotential
that has been generated using density functional theory for the atom, but with an exchange-correlation
that is not the one that's used in the solid and I don't like it and a number of us don't like it.
But this is to say that, you know, given the pseudopotential that you are using,
those can be hard, like if you use
norm-conserving pseudopotential
for oxygen, they might require a cutoff of 80 Rydberg or 100 Rydberg.
If you use ultrasoft pseudopotential or PAW,
you might get away with 30 Rydberg, and there is a huge difference in the cost
but you need to make sure that also your cutoff
of the charge density is correct. There is not going to be a factor of four between the [inaudible].
Let me go to something else as well.
The case I wanted to delve a little bit deeper in this concept of methods.
That is what happens when you have, actually,
a Fermi energy separating the occupied states by the unoccupied states.
And, as I said, that makes our integrals, that here have been discretised as
Monkhorst-Pack meshes, much more tricky
because now we have to say if that state is occupied or not.
OK, so f is equal to one, occupied, or is equal to zero if the state is unoccupied.
So this is really a sharp discontinuity.
I was doing my happy integral here,
but things dropped to zero and this one cuts and this is very difficult.
So during the years there has been the concept of trying to smooth these discontinuities
by introducing a fictitious temperature in the system, you know. Suppose that we had a real temperature,
now, the occupation of my states would be Fermi-Dirac occupations,
and while at zero temperatures, the Fermi-Dirac is a step function that goes
from one or two, if you are, say, spin degenerate to zero at finite temperature.
It's not. It's smoother,
and so it means that this, you know, discontinuity has now become
much smoother.
So if you have a temperature in your system,
your functional is not any more the total energy, but is the Helmholtz free-energy,
that is equal to the total energy minus the temperature
times the entropy.
Well, this would be the... say, if I have
a Fermi-Dirac occupation would be f log f, one minus f, log of one minus f.
So this is the proper extension to a finite temperature formalism and often
the temperature, to remind ourselves that it's not a real physical temperature,
we call it a smearing, sigma.
That is just there to smoothen our discontinuity.
Quantum Espresso, to confuse you, calls this quantity, the total energy also for a metal,
but remember, it's not a total energy it's a
total....
total free energy.
And so, all of a sudden, what happens is that now you have an additional thing to figure out.
That is how much temperature or how much smearing you have to give to your system.
And things are delicate.
Here I am looking at a physical quantity.
I think this is an iron bulk system with a few atoms in the primitive cell.
I am displacing one. I am measuring... calculating the force on an atom.
So you see, what I am plotting here is
the force on an atom as a function of the smearing. You know, .02 Rydberg
is equal to .27 electron volts,
that is roughly equal to... well, let's say 3000 K.
So if we think in physical terms,
you know, already .02 Rydberg of smearing is an enormous temperature.
So what is the right result here? You see, what I'm seeing is that the force that I
get changes if I change the smearing, and changes if I change the k-point sampling.
What I haven't shown you is that if we go to very,
very high smearing, we would always get one curve, because at very high smearing,
my discontinuities become so smooth that I can integrate it very well.
And you see a remnant of this here.
At .02, you have three or four curves
that are almost identical, the ones with 10 or 20, 24, 36 point meshes.
But as I decrease the smearing, you know, all this curve starts to become different
Even the best curve that I have, the best curve is the green diamonds here and the black dots here.
You see these two best curves are identical
in the high temperature, high smearing region, and then they start to separate.
And as you go to very small smearing, things go actually berserk.
So it's not that, you know, getting to zero smearing, gives you a good result.
You are here. You are all over the place.
What you would want ideally is the limit
of infinite k-point sampling and then zero smearing, but you cannot reach that.
So, you know, what you need to figure out roughly... let me remove these very confusing pen marks...
is finding, if you want the plateau around here and, you know, you need to look at what is,
you know, the error that you can tolerate here,
and this is a matter of, you know, around .001, .002 electron volt
so, it's not that big error.
But you need to look at the plateau where,
you know, with the highest sampling, you roughly have the same results.
And that's a function of smearing
things are not changing too much. If you go to larger smearings
it's easier to converge
but you converge to the wrong result because you have a lot of temperature.
If you go to smaller smearing,
you have the right temperature, but your sampling error is killing you.
So here, sampling is killing you.
And here,
temperature is killing you.
It's like being in a place where lava falls into ice water,
You need to find that comfort zone and it is very thin in between these two.
But, you know, to make things a little bit
more reassuring,
I would say that, you know, this interval of smearing tends to be conservatively good.
So never go below, never go above and make sure you have enough k-points.
And the reason why this is actually a good range is because we have developed 20, 25 years ago
what are called generalised entropic formulation that I'll just
briefly go through. That is we have been able to construct
formulation that are not leading to a Fermi-Dirac entropy, are leading to a different entropy.
And they are such that A,
the free energy as a function of the fictitious temperature sigma, is actually quite flat.
That means that if A is flat as a function of temperature, I can do my calculation
at high temperature where my k-point sampling is converged,
but this number that I get is actually very, very close to the number I can get at zero temperature.
How do we do that?
Well, we write the generalised Helmholtz free energy,
and we look at the minimum conditions of these functionals.
So I minimise it.
I have two Lagrange multipliers on the total number of electrons on orthogonality
and I have the derivatives, functional derivatives, with respect to the wavefunction
and I have the derivative with respect to the occupations.
And in particular, these derivatives
gives me a constructive relation between the expectation value of the Hamiltonian,
that's not the eigenvalues, the expectation value of the Hamiltonian on the psi i
minus the Fermi energy, the chemical potential, mu, and this is equal to T times
the derivative of the entropy with respect to the occupation.
So this can be written also
epsilon i minus mu divided by T needs to be equal to the derivative of the
entropy with respect to the occupation.
So from this
we can do something very simple, that is, you know, these ideas of smearing,
this idea of an occupation being a smooth function of the temperature,
And this variable is epsilon minus mu versus T was introduced already in the mid eighties by Fu and Ho
that said, well, you know, let's take a broadening of a Dirac delta, it could be a Gaussian,
and we calculated the occupation as this integral.
You see, for epsilon is very large,
this integral is going to give me one... Sorry, this is...
Yeah, so.
If epsilon is very large, this integral is going to give me zero,
and if epsilon is very negative, this integral is going to give me one.
If epsilon is equal to mu this integral
is going to give me one half, and, you know, they proposed different forms.
Let me not get into this now.
And...
So, remember what we had before,
so we had this constitutive relation between the expectation value of the Hamiltonian and the entropy.
We have defined the occupation as an integral
of a broadening function g that we decide. This is a broadening function.
And there is a very beautiful relation that Alessandro De Vita
and Mike Gillan developed in the early 90s that then tells me, well, if that is the case,
so if f is written as an integral of a broaden
occupation, the entropy is the integral of minus
whatever variables t times the broaden occupation g(t).
So there is if you want a very simple way to go
from broadening g to entropy S so you can choose your broadening.
You can take your Dirac delta and you can broaden it in very many different ways.
And this construction will tell you what is
the form of the entropy that needs to go into your Helmholtz free energy here.
And so we use this concept to build
broadening function
delta tilde that had a number of quality,
that is, of course, is normalised, that for which this integral was giving zero.
That means that the entropy at zero temperature is zero
but in particular, we built
distribution functions, broader delta that add a second moment that is equal
to zero and adding this second moment equal to zero was demonstrated by Stefano de Gironcoli
that is equivalent to having the free energy, not be, say linear or quadratic as a function
of temperature, but actually to have only cubic or quartic terms as a function of temperature.
And this is an example of the broadening function.
This is called cold smearing in Quantum Espresso.
There are other variation on the theme.
There is the Methfessel and Paxton who's actually introduced first and doesn't go through this.
If you want conceptual derivation of a free energy functional,
the Methfessel-Paxton has slight defects that can give rise to negative occupancies.
That is something that we do not like but both
of these broadening distribution have a zero second moment.
That means that the free energy does not depend on temperature up to the cubic order.
And I have a plot here from the case of aluminum in which I plot here E minus T S and I plot here E
and you see how the free energy is quadratic in temperature.
The total energy is actually quadratic with the opposite coefficient.
It was realised very early on that if we were to take the semi-sum, this is what Mike Gillan called
the corrected energy, E minus 1/2 T S, we would have, you see, with very large smearing...
now these are electron volts... a very good estimate of the zero value term.
Remember that this is the no-go zone where we can't do good sampling in there.
But this is just good for the total energy.
The functional is still the free energy,
so the stress and the forces that would come out of would not be corrected,
but the generalised smearing formulation, like the cold smearing, the Methfessel-Paxton,
do both.
They both give you a free energy that in the case of the cold smearing, say, you see, it goes like this.
You see how much flatter it is than the free energy.
The free energy, you see it's exploding,
and remember, this is the region where we are doing our calculations.
But because it's a variational functional,
we can calculate also forces and stresses and all derived quantities.
So let me summarise, what are, you know,
all the ingredients that you'll find in your homework to do a self-consistent calculation.
You need to figure out what is the translational symmetry of your system?
What is your Bravais lattice?
What is the group of atoms that we call the basis for the crystal structure?
That is, the atoms that are periodically repeated
by the Bravais translation, and they will give rise to the entire real crystal.
So for silicon, you would have a basis of two atoms.
For gallium arsenide, also a basis of two atoms and an FCC Bravais lattice.
For aluminum, we would have one atom in the FCC Bravais lattice.
For BCC iron, you would have one atom and then BCC Bravais lattice,
but, of course, if you want to study a defect
in silicon, then you might need to construct a much bigger primitive unit cell with maybe
sixty four atoms where then you remove one of these atoms.
So you need to construct your crystals.
You need to get your pseudopotentials from somewhere.
From the SSSP, from the Quantum Espresso website, from the PseudoDojo website.
You need to find out what is the correct cutoff energy for your wavefunction.
The SSSP website gives you suggestions.
What is the correct cutoff for the charge density?
Four times larger for now, concerning pseudopotential
but may be more and more for ultrasoft pseudopotentials.
What is the appropriate Monkhorst-Pack mesh?
With this, the fictitious temperature that you can use.
And then you need to make sure that you have
a good technique to get to self-consistency.
And I'll discuss that in a moment.
If you are desperate, you can go to the Materials Cloud.
In the the work tool, there is something called the Quantum Espresso input generator
that allows you to very simple...
very simply upload a structure, choose a table of pseudopotential for PBE or PBEsol,
either efficiency, that is less accurate but less expensive, or accuracy, more accurate.
And then you just have to decide if your system is
insulating or metallic and if it's magnetic or non-magnetic.
And you decide one of the standard three meshes of k-points and Quantum Espresso
generates... the Quantum Espresso input generator generates for you the...
the input file.
And if you want to plot band structure, there is another tool called SeeK-path,
that gives you also the appropriate crystallographic approved path in the Brillouin zone.
So let me conclude for today... I think I'm going a little bit longer.
let's try to do this in eight minutes...
with how to get to self-consistence and exactly as in Hartree and Hartree-Fock, in Kohn-Sham.
We need to get to the self-consistent solution of these equations.
So, you know, once we set up our problem with all these computational parameters,
we need to construct, you know, the Hartree operator, the exchange-correlation operator.
So we need to, you know, have actually a starting electronic density to then iterate,
and so this is typically given by sums of atomic densities.
At these points we can construct the Hartree operator, the exchange-correlation operator.
So we have the Hamiltonian. For that Hamiltonian,
we solve the Kohn-Sham equation.
This Kohn-Sham equation gives us a set of, you know, phi 1,
phi N for a molecule or phi 1k,
phi Nk for a solid with k belonging to the Monkhorst-Pack mesh.
Now with these
phis, we calculate the new charge density.
This was maybe rho 0.
Now we have rho 1 and with rho 1
we construct the Hartree and exchange-correlation operator.
We solve, we get a rho 2. We iterate, iterate, iterate
until things stop changing.
And, you know, this process can get to self consistency.
More often than not, this process would actually not converge
but you can get a little bit smarter and introduce a damping term.
You can say that
the new charge density that enters into this loop is not equal to the old
but is actually just equal to a fraction of the old
plus the remaining of the previous one.
So you sort of mix the charge density in an iterative way by making beta very small.
You slow down the process of
updating what was your initial starting condition, so you much more smoothly move around
and if you're lucky, it becomes easier to converge,
although you might also be slower.
And this is actually what is taken from the review of modern physics
of Mike Payne on the total energy pseudopotential plane-wave method from 1992.
This is exactly what I said, just in graphical form.
So but, you know, if you don't converge
and you decrease your beta and nothing happens, you just might be unlucky in your system.
You know, it's very tough to converge.
Maybe it's very large, maybe it has sloshing instability.
Maybe it has two levels, one empty and one full that, you know,
keep oscillating back and forth and you know which one to fill.
And so I wanted to introduce the complementary approach that maybe might be less efficient,
but is much more robust to converge total energy plane-wave calculation
that is based on variational recipes.
That is the idea that I've given before that rather than
iterating to self-consistency, the Hamiltonian and the Kohn-Sham equation,
you could do a direct nonlinear minimisation of the energy functional.
And to give you an example, I've written here the kinetic energy for a simple case.
I haven't put k-points.
I've just put, say, orbitals with just a discrete index n expanded in plane-waves.
And, you know, in the energy, as I already told you, there is,
you know, for the every plane-wave, we have a one half G square component.
And so the expression of the kinetic energy for a
given orbital and then summed over all
the orbital is, you see, sum over all the orbitals, one half sum over all
the plane-waves, the square modulus of the coefficient of the plane-waves.
So the kinetic energy is quadratic in the coefficient of the plane-waves.
But there is not only the kinetic energy, there is also the potential energies in the system.
And now the potential energy is not diagonal like the kinetic energy in the basis of plane-waves.
So the matrix element between G and G prime is not zero as in the kinetic energy.
So the Hamiltonian is... has kinetic energy on the diagonal and potential energy everywhere.
This is the Hamiltonian in a basis of plane-waves
and the matrix element of the potential are nothing else.
You see, this is just the coefficient of the Fourier transform of the potential in real space
with respect to a plane-wave exponential of i G prime minus G.
So if we put the two together and if for a moment
the Hamiltonian and the total energy of the system was just given by a kinetic energy term
and a potential term without considering that in Kohn-Sham we have
the external potential, the Hartree potential, the exchange-correlation potential.
But you see, there is a little star here that's difficult to see.
This would be the expression of the energy and we
can minimise it with respect to all these coefficients.
So finding the ground state
can be seen not as an iterative process, but as a minimisation of this form that is ...
here it looks like a quadratic form on the coefficients, cG, but whilst we have
not only the, say, external potential, but we have the external, the Hartree
and the exchange-correlation in this potential, the charge density enters
so this coefficient enters inside the Fourier
transform of the potential, making this a non-linear minimisation problem.
But, you know, it's still a problem that, you know, people have studied and so we can try to minimise this stuff
in the space of all the possible coefficients.
So I've written it here again, my approximate, here, energy function because I
don't have Hartree and exchange correlation so I can minimise it in the space
of all these coefficients, so it's a a multidimensional minimisation.
Remember that you might have
thousands or tens of thousands or hundreds of thousands of coefficients.
How do I minimise in a space of coefficients?
Well, maybe I'm somewhere here in the space of coefficient,
that is, I have certain values and I need
to evolve to get to the minimum of the energy in this multidimensional space of, you know,
10,000 dimensions and each coefficient has its own dimension.
And I need to get to the global minimum.
But the process is similar to what you would see in, say, this case,
one-dimensional, two-dimensional system, and maybe we have a two-dimensional space
of coefficients and we want to get to the minimum of this energy in the space of c1 and c2.
There are a lot of techniques
but one thing that you want to know is what is the derivative of the energy with respect to c1 or c2
because that gives you the right direction, the right gradient, where to move.
So in general, this is written in a compact form.
We need to know what is the derivative
of the energy. In this case, a functional derivative with respect...
maybe these are just real coefficients and, actually, it turns out these derivatives is just
given by the application of the Hamiltonian to psi i.
So if I have a wavefunction
the application of the Hamiltonian gives me the direction where to move minus.
And so I can start to find ways with these derivatives to go to the bottom.
And there are sort of two main ways to go to the bottom.
One is to say, well, if I'm here and I feel
a force like this, I can use that force as an acceleration.
So, you know, I'll start moving in the space of coefficients towards the bottom,
but because I'm feeling the force of
some acceleration, by the time I'm here, I have zero acceleration, but I have still a velocity.
And so, you know, without any damping, I would oscillate forever.
And so I need to put some damping here to make sure that at the end I converge to this ground state.
The other possibility would be to say that
I instead move with a velocity,
this is, you see, acceleration versus velocity. One dot or two dots.
So if I have a force that is proportional to the velocity, again I'll start moving,
but a force that is proportional to the velocity
means that I have critical slowing down. As I get closer and closer to the minimum.
I never reach it because my driving force goes to zero.
So the first recipe that we call Car-Parrinello
and the second recipe that we call conjugate gradients have the defect, one,
that Car-Parrinello overshoots the minimum. That has to be damped down.
Conjugate gradient has a critical slowing down and it has to be pushed to get to the minimum.
But these are the two main techniques:
conjugate gradient comes again from some mathematical analysis in which if you have, say,
a potential energy surface,
that is very asymmetric and you use the gradients to get to the bottom.
So you see here, I'm moving along the gradient until the derivative along that direction is zero.
Then I move along the gradient and move along the gradient. I zigzag to the bottom
and it's much better once I've taken the first step not to move in the
steepest descent direction, but to move in what is called a conjugate
direction that takes the steepest direction, but tries to make sure that as you move
in the direction, you don't lose the minimisation that you have done before.
You see fine, we're here.
And moving down the steepest direction,
it's true that I would improve around that,
but I would do more so with respect to what I done on the previous iteration.
So conjugate gradients...
Thus the conjugate gradients techniques
have what have been used in a lot of codes, CASTEP, especially, to get in a very robust way
to the ground state and thus give that when you have to the ground state,
we'll see that tomorrow, you can calculate forces,
and then maybe you let atoms move using molecular dynamics techniques.
But
exactly
if I'm
looking at a system during a time evolution.
So if the atoms move, every time an atom move I have a new potential energy surface,
and so you have to minimize. I have to minimize.
And so getting to the minimum
at every step is an expensive operation.
And that's the idea that came
to Car and Parrinello in the mid 1980s,
that is to say, well, let's try to find a way by which we don't have to reach the minimum
of this, you know, total energy as the function of the coefficient space at every step,
but with each of the coefficient, how to evolve themselves so that they follow,
maybe not exactly, but they follow very closely
what would be their minimum at every instant, at every time step.
So this could be time zero, this could be time zero plus delta T,
this could be time zero plus two delta T and so on and so forth.
And Car and Parrinello did that through the idea of extended Lagrangian, that is,
rather than having a standard Lagrangian with the kinetic energy minus
the potential energy of our system, that is the electronic energy, the Kohn-Sham energy.
What they did is that they introduced
dynamics, also for the coefficients
of, sorry, of the orbitals. You see here there are the velocity of the orbital.
So this is like a physical kinetic energy for the atoms
and this is a fictitious kinetic energies for the plane-wave coefficients.
But by making the mass of these fictitious plane-wave coefficients
that are now thrown in into my dynamical algorithm by making this mass very, very small.
It's like making them the flies in the cow. The cow that are the atoms can move,
and these coefficients don't affect the dynamics,
but just follows the atoms and very quickly rearrange themselves in that blue
line at the bottom of the evolving potential surface.
So this would be the extended Car-Parrinello Lagrangian. It gives rise to equation of motion,
not only for the atoms, but also for the coefficients that learn how to follow.
And once you look at the,
actually, constant of motions, you discover that because you are
integrating the extended Lagrangian, you have a constant of motion where the
kinetic energy of the atom sees the potential energy of the atom is.
But there is also these
kinetic energies of the coefficient of the plane-waves,
and that's, you know, it's just a computational trick,
but by making this effective mass very small, it can be negligible.
And also the blue curve that would be this
expression becomes closer and closer to the red curve
that would be the right expression.
The only cost that you have to pay if you make a mu too small
is that the coefficient of the plane-waves change so fast that you're time step
for the Car-Parrinello integration needs to become very small.
And that's a bit of a technicality.
But if you're able to do this calculation well, things can become very efficient.
You have a very flat physical constant of motion,
because the total constant of motion is very well preserved.
And this, you know,
unphysical quantity is very small and very flat. These kinetic energies,
you know, we have hundreds of thousands of coefficients
of plane-waves, but their total kinetic energy can be a tiny fraction,
one tenth, one twentieth of the true physical kinetic energy of the ions.
So the coefficient of the plane-waves follow, but they don't give back energies to the ions.
So in a nutshell, you can do
Born-Oppenheimer
molecular dynamic simulation, where you move the ions. You converge the electrons back to...
to the ground state with, say, iterative procedure or conjugate gradient procedures,
but you're always, sort of, you know, you have always a little bit of a systematic error,
so your constant of motion will not be very constant.
Instead, with Car-Parrinello techniques,
your constant of motion is exact but you need to make sure that you have integrated them.
Things can go well in a semi-conductor or insulator.
The issue, and then I'll conclude, is what happens
if you have a metal. And if you happen, what if ... what happens if you have a metal
in a variational algorithm is that you start to have three degrees of freedom,
that is, you need to evolve your wavefunction but you need also to evolve your occupations.
And you need that to make sure that your occupations are calculated on the eigenvalues of
the Hamiltonian, so, you know, what is the representation that diagonalises
the Hamiltonian, so that you can get the actual eigenvalues of the system.
So in a metal, as opposed to a semiconductor
or an insulator, you don't have to evolve only the wavefunctions.
In a semiconductor or insulator, the total energy is invariant for unitary transformation.
So you don't care about that.
But instead in a metal, you need to make sure that you evolve
also the occupations that can change from zero to one, or they can be fractionally.
And you always need to know what is
the eigenvalue of the Hamiltonian that is unique to being the representation where you can
calculate eigenvalues because the problem is not invariant anymore.
So you have, if you want, two additional degrees of freedom,
occupation numbers and unitary rotations.
And these occupation numbers can be very effectively combined with the unitary
rotation to, say, in a metal, the additional degree of freedom is actually a matrix
that is the representation.
Say, if this was a standard finite temperature problem,
it would be the representation of the Fermi-Dirac, operator on the basis of wavefunction psi i.
And, in general, this matrix is non-diagonal.
And so to deal with this non-diagonal matrix, you need to write the charge density in a covariant
form like this, and you need to write the total energy
and the external potential in a covariant form, as a double form.
And of course, you could always
choose a special representation in which this f is diagonal, but then you would need also
to evolve the unitary matrices that connect the psi H psi
with its eigenvalues, epsilon i.
So there is no winning.
I mean, you have all this complexity.
And the nice thing about this covariant
representation is that it allows you to define this functional G
that is defined as the minimum of the free-energy functional with respect
to this matrix representation of the occupations.
And you see this functional G is a functional where occupations have disappeared.
Not only that but is a functional that is
invariant for a unitary transformation of the orbital.
So if I take the psi and I send them into psi U, U is a unitary matrix,
I have done a unitary transformation.
What happens is that the values of this functional
doesn't change because if I do a transformation of the psi into psi U,
what I can do is a transformation of the U into U dagger f U and the value has not changed.
So if I rotated psi, the f counter rotates and the minimum doesn't change.
So this trick gives me a functional G that looks like the functional of a semiconductor or insulator.
It doesn't have explicit occupation numbers and it's invariant for unitary transformation,
so I can do conjugate gradient iterative
techniques and so on in this and I can get to very stable molecular dynamics again.
This is called... we called it ensemble DFT.
It's coded in VASP.
This was actually done by Freysoldt and collaborators. It's coded in CASTEP.
It's not coded in Quantum Espresso.
That's somehow very sad
but I think in a few months we'll be able to make a release of this together with the GPU library
serials that Anton Kozhevnikov and collaborator are developing at CSCS.
So, I think this, and that's the second lesson, I've been only 45 minutes longer.
I think there were a lot of technicalities
so I figure out that only people that want to do actual calculations
were going to follow this and be interested in this.
Tomorrow we'll go back to 2 hrs
and a much more conceptual discussion. This will be recorded, will be put on the Materials Cloud.
If you want to do the calculation, go to the Materials Cloud.
In the fireside section, there is a link to GitHub.
There is a complete handout that tells you how to download and install Quantum Espresso
or how to download and install the Quantum Mobile, the virtual machine with Quantum Espresso
preloaded into it and one by one there are homeworks that tells you and, you know,
teach you how to set up the cutoff for your calculation,
how to check on the k-point sampling,
how to do a simple calculation, an equation of state for sodium chloride.
If you are lost, the Quantum Espresso input generator always gives you a good start
to, you know, have a good chance for an input file for Quantum Espresso.
The SSSP library gives you tested pseudopotential and tested cutoff.
And there is always a mailing list where you can ask questions for the Quantum Espresso users.
And if you're interested also in more on this, we have put on the lower section, also,
you can enrol in a mailing list for
more fireside chats. With this I'll conclude. The panellists have collected all your questions.
We'll put them in a PDF in the next few days.
We'll find them and we'll post them next week.
Other than this, I look forward to seeing you tomorrow.
And this we'll also post on the lower section.
Thanks a lot and thanks especially exactly
to Davide, Francisco, Luca and Norma that have prepared the handout
and then been sitting here all of these 2 hours and 45 minutes listening to this,
and, you know, see you indirectly tomorrow.