WEBVTT
1
00:00:04.040 --> 00:00:05.380
OK.
2
00:00:05.480 --> 00:00:08.260
Welcome back, everyone.
3
00:00:08.360 --> 00:00:14.620
Today, we have our second part of our
fireside chat on density functional theory
4
00:00:14.720 --> 00:00:19.720
and, in particular, it's the one that
I call density functional practice.
5
00:00:19.820 --> 00:00:26.800
Just to remind everyone, what is
the plan for this three webinar series.
6
00:00:26.801 --> 00:00:30.386
Yesterday, we had a more conceptual part
7
00:00:30.387 --> 00:00:34.740
on electronic structure methods and
density functional theory, in particular.
8
00:00:34.840 --> 00:00:37.460
Today will actually be dedicated to the practice.
9
00:00:37.560 --> 00:00:45.060
What do we need to do this calculation using
modern plane-wave pseudopotential code.
10
00:00:45.160 --> 00:00:49.860
And tomorrow it will be, I think,
again, conceptual understanding.
11
00:00:49.960 --> 00:00:52.960
What are the properties that we can calculate,
12
00:00:53.720 --> 00:00:59.340
those that we cannot calculate and for those
that we can, how accurate we are.
13
00:00:59.440 --> 00:01:05.480
Yesterday, I think we reached very quickly
the limit of 500 participants.
14
00:01:05.580 --> 00:01:14.860
So I am very sorry that some of you couldn't
attend them, but we have already put the slides
15
00:01:14.960 --> 00:01:21.360
and the talk and also just the audio
on the Materials Cloud learn section.
16
00:01:21.460 --> 00:01:23.340
And it's also actually on YouTube.
17
00:01:23.440 --> 00:01:28.040
And we'll do the same today and
18
00:01:28.640 --> 00:01:30.240
later on
19
00:01:30.400 --> 00:01:33.180
for tomorrow's webinar.
20
00:01:33.280 --> 00:01:37.600
Let me actually also monitor the chat.
21
00:01:39.920 --> 00:01:41.400
Yeah.
So
22
00:01:41.560 --> 00:01:43.520
we have as yesterday
23
00:01:43.640 --> 00:01:49.180
four panellists, Francisco Ramirez,
Norma Rivano, Luca Binci and Davide Grassano.
24
00:01:49.280 --> 00:01:53.300
They will monitor on my behalf the chat.
25
00:01:53.400 --> 00:01:59.700
We don't have the time to answer questions in real
time, but they will be collecting if something
26
00:01:59.800 --> 00:02:02.980
very simple they'll answer in real time in the chat.
27
00:02:03.080 --> 00:02:11.586
But by next week, you'll also find all the questions
that we deem relevant, we collected 40 yesterday,
28
00:02:11.671 --> 00:02:18.240
and answer in a PDF file, also on the Materials Cloud.
29
00:02:18.560 --> 00:02:23.200
Today webinar is also special
in the sense that there is homework
30
00:02:24.240 --> 00:02:27.860
if you want to practice some of this calculation,
31
00:02:27.960 --> 00:02:31.760
there are instruction on the Materials Cloud.
. There is a GitHub section
32
00:02:31.920 --> 00:02:34.220
where you can download the handout
33
00:02:34.320 --> 00:02:40.328
with instruction on how to either install
Quantum Espresso or just install
34
00:02:40.329 --> 00:02:47.471
Virtual Box and run a virtual machine with Linux
environment on whatever machine you have,
35
00:02:47.472 --> 00:02:55.180
either Windows, Macs, Mac or Linux
and to go through actual calculations.
36
00:02:55.280 --> 00:03:02.260
And I will focus here on teaching
you how to do these calculations.
37
00:03:02.360 --> 00:03:06.240
One thing.
One slide that I left out from yesterday.
38
00:03:06.340 --> 00:03:08.960
I was planning to present it tomorrow,
39
00:03:09.060 --> 00:03:16.560
but maybe this is relevant already,
is a set of pointers on educational material.
40
00:03:16.660 --> 00:03:19.140
Here is, let's say, a long list.
41
00:03:19.240 --> 00:03:21.780
There are some online resources.
42
00:03:21.880 --> 00:03:27.680
I've just pointed out to some girls,
say, that I myself ran a class
43
00:03:27.780 --> 00:03:31.180
on atomistic simulation for many,
many years, at MIT actually.
44
00:03:31.280 --> 00:03:33.760
Some of the material here is based on that class
45
00:03:33.860 --> 00:03:38.140
that we kept teaching here in Lausanne or in Berkeley.
46
00:03:38.240 --> 00:03:42.700
And so you can use that there are
videos or great lecture notes.
47
00:03:42.800 --> 00:03:47.340
We accumulate materials on the
Learn section of the Materials Cloud.
48
00:03:47.440 --> 00:03:49.940
There are plenty of other online resources.
49
00:03:50.040 --> 00:03:58.000
Stefaan Cottenier runs an online stream class
on the field and you can just google around.
50
00:03:58.200 --> 00:04:03.000
In terms of books, these are some of my favorite
books on fundamentals of quantum mechanics.
51
00:04:03.100 --> 00:04:08.140
The two books by Bransden and Joachain
are I think particularly good.
52
00:04:08.240 --> 00:04:11.200
If you want to learn more about the fundamental
53
00:04:11.300 --> 00:04:15.220
of density functional theory the
Parr and Yang book is very good
54
00:04:15.320 --> 00:04:18.540
or if you want to learn about advanced
55
00:04:18.640 --> 00:04:21.420
electronic structure methods the two
books by Richard Martin and Lucia Reining
56
00:04:21.520 --> 00:04:23.860
and David Ceperley are very good.
57
00:04:23.960 --> 00:04:29.420
And if you have more an eye on
materials simulation the books by
58
00:04:29.520 --> 00:04:31.540
Efthim Kaxiras, Jorge Kohanoff and Feliciano Giustino.
59
00:04:31.640 --> 00:04:33.000
And
60
00:04:33.120 --> 00:04:39.220
if you needed just one book to bring
with you on a desert island,
61
00:04:39.320 --> 00:04:44.486
it would be either the Parr and Yang if you want
to, you know, learn about the fundamentals.
62
00:04:44.487 --> 00:04:49.460
Or it would be the Giustino if you
wanted to learn about the practice.
63
00:04:49.560 --> 00:04:56.880
And if it's a really nice desert island I
can come and actually teach you on the fly.
64
00:04:57.200 --> 00:05:01.820
Now let me sort of get to where we left it yesterday.
65
00:05:01.920 --> 00:05:06.340
So the Hohenberg and Kohn theorem introduce
66
00:05:06.440 --> 00:05:14.560
the fundamentally exact concept that we can deal
with an energy that is a functional of the
67
00:05:14.660 --> 00:05:18.380
charge density and we have
even a variational principle.
68
00:05:18.480 --> 00:05:23.580
The only catch is that we don't know what that
universal function al of the charge density is.
69
00:05:25.471 --> 00:05:27.186
Kohn-Sham
70
00:05:27.360 --> 00:05:34.060
do a very cunning attempt to
identify how we could approximate
71
00:05:34.160 --> 00:05:42.000
that functional breaking apart
the universal functional F into three terms
72
00:05:42.320 --> 00:05:48.920
T s - kinetic energy of the non-interacting
Kohn-Sham electrons, a Hartree term and then they
73
00:05:49.020 --> 00:05:55.820
push everything that is not known into the third
term - the exchange-correlation functional.
74
00:05:55.920 --> 00:06:01.080
Now in doing this they reintroduce
orbitals. You see it's not any more
75
00:06:03.200 --> 00:06:05.160
a functional theory of the density.
76
00:06:05.260 --> 00:06:09.780
Now there are just again, single particle orbitals.
77
00:06:09.880 --> 00:06:10.960
It's still exact.
78
00:06:11.060 --> 00:06:14.160
Everything is being compressed back
79
00:06:14.260 --> 00:06:19.740
into the exchange-correlation functional
and for the exchange-correlation functional,
80
00:06:19.840 --> 00:06:26.920
we need to find approximations. And, as always,
we have a functional so we have also
81
00:06:27.020 --> 00:06:31.160
the corresponding Euler-Lagrange
equation that today we'll actually
82
00:06:31.260 --> 00:06:36.420
solve in real life, and so
the equations that correspond
83
00:06:36.520 --> 00:06:41.660
to the minimum of that functional, we call
them in this context, the Kohn-Sham equation.
84
00:06:41.760 --> 00:06:42.980
They look very familiar.
85
00:06:43.080 --> 00:06:46.660
We have a quantum kinetic energy term.
86
00:06:46.760 --> 00:06:48.820
We have a Hartree potential.
87
00:06:48.920 --> 00:06:54.340
We have an external potential and
everything that we had pushed back
88
00:06:54.440 --> 00:07:02.520
into this unknown exchange-correlation energy
now becomes in the Euler-Lagrange equation,
89
00:07:02.960 --> 00:07:05.200
a functional derivative with respect to the
90
00:07:05.201 --> 00:07:10.420
charge density and so we call that
the exchange-correlation potential.
91
00:07:10.520 --> 00:07:13.540
So this is just exchange-correlation.
92
00:07:13.640 --> 00:07:16.040
And as I said, you know,
93
00:07:16.280 --> 00:07:19.080
Kohn and Sham were both bright and lucky.
94
00:07:19.180 --> 00:07:23.000
They followed Fermi's idea of using the local
95
00:07:23.100 --> 00:07:28.600
density approximation, but not for the kinetic
energy: for the exchange-correlation.
96
00:07:28.700 --> 00:07:31.500
So for anything that is
97
00:07:31.600 --> 00:07:36.660
non-homegeneous we construct
the exchange-correlation energy
98
00:07:36.760 --> 00:07:43.700
for that system by going point by point
and at every point assigning a density
99
00:07:43.800 --> 00:07:51.320
of exchange-correlation energy that is the same
density of the homogeneous non-interacting gases.
100
00:07:51.420 --> 00:07:54.957
So we build with all these infinitesimal Lego bricks,
101
00:07:54.971 --> 00:08:00.520
each one at its density, the total
exchange-correlation energy.
102
00:08:01.120 --> 00:08:05.180
So, why they were brighter
and lucky at the same time?
103
00:08:05.280 --> 00:08:09.240
Well, because the local density approximation
104
00:08:09.340 --> 00:08:19.260
for the exchange-correlation satisfy first
and very beautifully a fundamental constraint
105
00:08:19.360 --> 00:08:24.420
on the exchange and correlation that,
you know, they didn't know at the time.
106
00:08:24.520 --> 00:08:27.000
But basically and you can imagine that in fairly
107
00:08:27.100 --> 00:08:32.940
intuitive physical terms, what happens if
you are an electron sitting somewhere.
108
00:08:33.040 --> 00:08:37.180
Basically, you are creating around you
109
00:08:37.280 --> 00:08:42.280
a lack of probability of finding another electron,
because, you know, electrons repel each other.
110
00:08:42.560 --> 00:08:44.160
I mean, if you are, you know,
111
00:08:44.600 --> 00:08:51.420
squeezed in an underground train,
it's very unlikely that you will find
112
00:08:51.520 --> 00:08:55.820
someone, you know, 1 mm close to you
unless, you know, you want to be there.
113
00:08:55.920 --> 00:09:01.900
So, in general, because electron repel each other,
they create around themselves
114
00:09:02.000 --> 00:09:05.180
a hole that is called the exchange-correlation hole.
115
00:09:05.280 --> 00:09:13.800
So first the LDA satisfy this exact -1 constraint,
that the exchange-correlation hole
116
00:09:13.960 --> 00:09:15.840
should integrate to -1.
117
00:09:16.080 --> 00:09:17.240
And that's the first thing.
118
00:09:17.240 --> 00:09:18.520
But it's not the only thing.
119
00:09:18.620 --> 00:09:24.660
And actually there are so many beautiful
calculations done by the group of Richard Needs
120
00:09:24.760 --> 00:09:31.400
in Cambridge looking at this exchange-correlation
hole with accurate stochastic methods.
121
00:09:31.500 --> 00:09:35.940
With diffusion Monte Carlo,
and here this exchange-correlation hole
122
00:09:36.040 --> 00:09:41.999
is done for a, you know, number of interacting
electrons in a very simple potential.
123
00:09:42.000 --> 00:09:45.460
In this sinusoidal potential that is here.
124
00:09:45.560 --> 00:09:46.660
And, you know,
125
00:09:46.760 --> 00:09:51.660
what we are asking ourselves, if I'm on
the ground state and I put myself here,
126
00:09:51.760 --> 00:09:58.420
what is the decrease of probability that surrounds
me, this exchange-correlation hole?
127
00:09:58.520 --> 00:10:03.535
And I can calculate it in Quantum Monte Carlo,
and that's what you have on the left,
128
00:10:03.735 --> 00:10:07.100
the Monte Carlo hole, and they
can calculate it with LDA.
129
00:10:07.200 --> 00:10:14.660
So first of all, LDA integrates to the right quantity,
but it also starts to look seriously good.
130
00:10:14.760 --> 00:10:21.540
Now, the catch is that if I now look not at the
exchange-correlation hole where I'm sitting
131
00:10:21.640 --> 00:10:26.771
at the top of this potential, but I move around,
you'll see a little movie of these holes.
132
00:10:26.772 --> 00:10:30.300
As I move to the bottom and I go back to the top.
133
00:10:30.400 --> 00:10:37.800
On the left, you see the exact results
and on the right you see the LDA results.
134
00:10:39.040 --> 00:10:44.820
And as you can see, I mean, it's totally so
beautiful that I want to run it twice.
135
00:10:44.920 --> 00:10:48.160
You see here, it's a desperation
136
00:10:48.480 --> 00:10:50.420
so and at the bottom of the potential.
137
00:10:50.520 --> 00:10:56.300
the two exchange-correlation hole
looks very different. And actually I won't do it
138
00:10:56.400 --> 00:11:00.880
but through thermodynamic integration,
you can see that you can derive the exact
139
00:11:00.980 --> 00:11:07.700
exchange-correlation energy from the exact
exchange-correlation hole so the two are related.
140
00:11:07.800 --> 00:11:12.960
But as you saw, LDA here was doing a very good job
141
00:11:13.060 --> 00:11:18.180
but at the bottom of that potential, does a very poor job.
142
00:11:18.280 --> 00:11:21.620
With these it's a saving grace.
143
00:11:21.720 --> 00:11:28.180
Well, it's the fact that actually in the exact
expression for the exchange-correlation energy
144
00:11:28.280 --> 00:11:34.460
coming from the exchange-correlation hole, the
only thing that matter is the spherical average.
145
00:11:34.560 --> 00:11:42.207
So it doesn't matter the explicit form of its hole,
but what matters is the spherical average of this hole
146
00:11:42.407 --> 00:11:46.680
and the spherical average
of this hole is much better behaved.
147
00:11:46.680 --> 00:11:53.780
These are again Quantum Monte Carlo results
and the LDA does quite a good job.
148
00:11:53.880 --> 00:12:01.960
So the idea from 1965 of Kohn-Sham of using
LDA turned out to be very successful.
149
00:12:02.060 --> 00:12:03.540
We have already seen
150
00:12:03.640 --> 00:12:10.540
Marvin Cohen calculations of the phase
diagram of silicon in its different phases.
151
00:12:10.640 --> 00:12:17.180
And, you know,
these are just examples from 20 years ago from
152
00:12:17.280 --> 00:12:22.760
Chris Pickard looking at the very simple
mechanical properties.
153
00:12:22.860 --> 00:12:25.820
In this case is just the equilibrium
154
00:12:25.920 --> 00:12:29.520
lattice parameter across a few
different classes of compounds.
155
00:12:29.620 --> 00:12:34.500
I show tomorrow much more extensive examples.
156
00:12:34.600 --> 00:12:41.120
But you see, what we have is that the lattice
parameter at the LDA level is actually very good.
157
00:12:41.220 --> 00:12:43.920
I mean, the errors that we do and there is no
158
00:12:44.020 --> 00:12:48.700
experimental input so this is
completely from first principle
159
00:12:48.800 --> 00:12:53.880
is of the order of, you know,
1% worse in the lattice parameter.
160
00:12:53.980 --> 00:12:58.460
That would be 3% on the volume.
161
00:12:58.560 --> 00:13:04.100
Elastic constant that, you know,
come from the derivatives of the energy.
162
00:13:04.200 --> 00:13:10.500
The second derivative around
this lattice parameters, our model will say
163
00:13:10.600 --> 00:13:15.660
starts to be a little bit trickier
but we'll see that tomorrow.
164
00:13:15.760 --> 00:13:21.140
But the bottom line here is that the LDA
was already incredibly good.
165
00:13:21.240 --> 00:13:23.960
So everything that we have done ever since,
166
00:13:24.060 --> 00:13:30.920
I'll briefly mention it in a moment,
is a tiny improvement on top of this.
167
00:13:31.020 --> 00:13:36.060
We might think that we have had major improvements
and it was very difficult to get, you know, those
168
00:13:36.160 --> 00:13:42.020
improvements that are very important
but LDA was very, very good.
169
00:13:42.120 --> 00:13:45.220
And just historically, you know, of course,
170
00:13:45.320 --> 00:13:51.840
if you have built an exchange-correlation
from a local density approximation and you want
171
00:13:51.940 --> 00:13:59.760
to make things better, it might be very, you know,
natural to do something, perturbative and say,
172
00:13:59.860 --> 00:14:05.660
let's add maybe the gradients of the charge
density to make things better and that was done.
173
00:14:05.760 --> 00:14:08.220
But at the beginning it was actually
174
00:14:08.320 --> 00:14:15.700
decreasing the quality of the results
and the reason for that was that actually
175
00:14:15.800 --> 00:14:23.240
the exact integral rule on exchange-correlation
all integrated to minus one was not satisfied.
176
00:14:23.340 --> 00:14:26.420
And that's why this concept of generalised
177
00:14:26.520 --> 00:14:31.620
gradient approximation was introduced
by John Perdew and others that is
178
00:14:31.720 --> 00:14:38.860
gradient approximation that preserved this
fundamental rule and that a lot of the effort
179
00:14:38.960 --> 00:14:45.360
in constructing the new exchange-correlation
functional has been driven by the desire
180
00:14:45.460 --> 00:14:49.800
to satisfy some of these constraints
and that can be very many.
181
00:14:49.900 --> 00:14:51.940
This is just one example.
182
00:14:52.040 --> 00:14:58.340
But once we had the say correct gradient
expansions that satisfy the sum rule,
183
00:14:58.440 --> 00:15:02.780
it was similar that we could actually
do quite a lot better,
184
00:15:02.880 --> 00:15:12.020
especially when it was a matter of comparing, say,
maybe energies or like here, a question of states
185
00:15:12.120 --> 00:15:17.200
under very different chemical regime or coordination regime.
186
00:15:17.300 --> 00:15:24.780
This is an example of the phase transition
from alpha-quartz to stishovite.
187
00:15:24.880 --> 00:15:29.060
Stishovite is a high pressure form of alpha-quartz.
188
00:15:29.160 --> 00:15:34.340
Again, it's found if you look at rocks under
meteorites, again, if you are lucky
189
00:15:34.440 --> 00:15:39.080
and you see these two equation of states, you
remember how well the transition pressure was
190
00:15:39.180 --> 00:15:43.140
coming for silicon from the
diamond phase to beta-tin.
191
00:15:43.240 --> 00:15:46.700
Well, here it turns out that if you use the LDA,
192
00:15:46.800 --> 00:15:53.200
the two equation of states have the minimum
at the same energy so you don't really need any
193
00:15:53.300 --> 00:15:58.820
pressure to go from one equation of state to the other.
194
00:15:58.920 --> 00:16:01.520
And that's really in disagreement with experiment
195
00:16:01.529 --> 00:16:08.914
and Don Hamann show, you see, more than 25 years
ago that, actually, if you did a GGA calculation of
196
00:16:08.915 --> 00:16:11.786
these two phases, the common tangent was giving
197
00:16:11.800 --> 00:16:18.560
a transition pressure that was in much
better agreement with experiments.
198
00:16:18.840 --> 00:16:28.280
And from then onwards has been a proliferation
of exchange-correlation functional. There is a
199
00:16:28.800 --> 00:16:32.600
LIBXC library in the solid-state community,
200
00:16:32.700 --> 00:16:39.800
developed by Miguel Marques, that I think
has reached almost 500 functionals.
201
00:16:40.400 --> 00:16:45.400
I will not try to describe in any sense all of them.
202
00:16:45.500 --> 00:16:48.020
We'll have tomorrow a better discussion.
203
00:16:48.120 --> 00:16:51.660
You have understood very well what LDA is.
204
00:16:51.760 --> 00:16:54.920
GGAs are exchange-correlation functionals
205
00:16:55.020 --> 00:16:57.894
that introduce this generalized gradient.
206
00:16:58.094 --> 00:17:01.540
The solid-state community,
the materials community,
207
00:17:01.640 --> 00:17:08.300
is actually quite fond of the PBE functional from 1996
208
00:17:08.400 --> 00:17:14.540
and if you do a lot of solid-state calculation,
a very tiny variation on that,
209
00:17:14.640 --> 00:17:20.560
that add a couple of parameters fit
not on molecular energies but on solid state
210
00:17:20.660 --> 00:17:27.360
structures is actually quite a bit better,
especially the number of materials like
211
00:17:27.480 --> 00:17:30.680
perovskites, oxides that would
come out better in PBEsol.
212
00:17:30.780 --> 00:17:36.700
So this would be my goto functional as a standard.
213
00:17:36.800 --> 00:17:41.120
I will discuss the van der Waals corrections
of this functional that are again,
214
00:17:41.220 --> 00:17:46.500
a very successful story of the last,
I would say, 10, 20 years.
215
00:17:46.600 --> 00:17:52.620
Tomorrow, I will mention some
fancier, more complex, extension.
216
00:17:52.720 --> 00:17:55.740
One that I'm fond of and no-one looks at,
217
00:17:55.840 --> 00:17:58.620
that is the weighted density approximation.
218
00:17:58.720 --> 00:17:59.860
Something that, you know,
219
00:17:59.960 --> 00:18:06.940
when it was introduced was quite cumbersome so it
never took off but it has very elegant properties
220
00:18:07.040 --> 00:18:14.180
of the response function so it will be still
a very interesting field of study.
221
00:18:14.280 --> 00:18:18.520
To go beyond the gradients, of course,
222
00:18:18.640 --> 00:18:21.000
the Laplacian could be, again, a natural step.
223
00:18:21.100 --> 00:18:24.180
We have to be careful on how to introduce them.
224
00:18:24.280 --> 00:18:30.220
Those were introduced only maybe
15 years ago, and around 5 years ago,
225
00:18:30.320 --> 00:18:37.480
meta-GGA functional that is called SCAN
(strongly constrained and appropriately-normed)
226
00:18:37.680 --> 00:18:39.880
has become very popular.
227
00:18:39.980 --> 00:18:48.860
And, you know, it improves in a number of cases
and a number of important average onto PBEsol.
228
00:18:48.960 --> 00:18:51.360
There are some subtleties that I'll mention later
229
00:18:51.460 --> 00:18:58.100
on that have to do with pseudopotential
and pseudopotential generation, because
230
00:18:58.200 --> 00:19:00.920
as you build those, you want to make sure that your
231
00:19:01.020 --> 00:19:06.940
exchange-correlation for those calculation
is consistent with what you do later.
232
00:19:07.040 --> 00:19:10.420
And something that has become very popular,
233
00:19:10.520 --> 00:19:17.000
starting from the chemistry community where
not only Hartree-Fock calculations are very
234
00:19:17.100 --> 00:19:24.640
common, but tend to be also very inexpensive
due to the use of the Gaussian basis set
235
00:19:25.440 --> 00:19:30.200
are hybrid functionals.
In the solid-state community,
236
00:19:31.080 --> 00:19:35.220
the Heyd-Scuseria-Ernzerhof
- the HSE functional is very common.
237
00:19:35.320 --> 00:19:38.460
Those mix a certain amount
238
00:19:38.560 --> 00:19:44.520
of exchange and correlation
coming from something like
239
00:19:45.240 --> 00:19:52.600
PBE or revised-PBE with a certain amount of exact
Fock exchange or often adjust the short-range part
240
00:19:52.700 --> 00:19:58.300
and not to inherit from Hartree-Fock
this tendency to make...
241
00:19:58.400 --> 00:20:01.940
to turn metal into insulators.
242
00:20:02.040 --> 00:20:07.460
I'll discuss more of this tomorrow so right
now I want just to leave you with
243
00:20:07.560 --> 00:20:15.540
zoology of exchange-correlation functional;
the comment that LDA was great to begin with,
244
00:20:15.640 --> 00:20:22.600
PBE and PBEsol make some subtle
but significant improvements.
245
00:20:22.760 --> 00:20:29.300
SCAN makes some improvement on top
that are still a qualitative feature,
246
00:20:29.400 --> 00:20:32.940
so that we'll discuss tomorrow, that
are related to self-interaction.
247
00:20:33.040 --> 00:20:39.680
And that is addressed in part by hybrids,
by Hubbard functional.
248
00:20:39.780 --> 00:20:45.940
We discuss them tomorrow . And with the meta-GGA
keep in mind also this caveat about
249
00:20:46.040 --> 00:20:50.440
the need, in principle, to use a consistent
pseudopotential that have been
250
00:20:50.540 --> 00:20:54.920
generated with the same exchange-correlation functional.
251
00:20:55.200 --> 00:21:01.780
But here, if you want, is all
the physical quality of your model,
252
00:21:01.880 --> 00:21:07.500
that is, if your exchange-correlation functional
were exact you would have the perfect energy.
253
00:21:07.600 --> 00:21:11.560
So anything that is not a numerical error,
254
00:21:11.660 --> 00:21:17.840
that is what mostly we'll discuss now
in the rest of the presentation,
255
00:21:18.080 --> 00:21:24.400
comes from either a physical shortcoming of your
exchange-correlation functional or, actually,
256
00:21:24.500 --> 00:21:28.460
from the physical quality of your
exchange-correlation functional.
257
00:21:28.560 --> 00:21:34.620
But sadly, there is not a systematic way
to make these exchange-correlation
258
00:21:34.720 --> 00:21:36.980
functionals better, although there are some very
259
00:21:37.080 --> 00:21:42.340
interesting and nice ideas based on what
is called the adiabatic correction.
260
00:21:42.440 --> 00:21:46.380
And so in some ways here,
is all the intuition, you know,
261
00:21:46.480 --> 00:21:52.320
there is the perspective of John Perdew
that you have to impose a lot of
262
00:21:52.760 --> 00:21:56.160
known constraints on your exchange-correlation
263
00:21:56.260 --> 00:21:59.340
functional and that is absolutely a very good idea.
264
00:21:59.440 --> 00:22:03.260
You see, there are still qualitative failures
265
00:22:03.360 --> 00:22:09.660
that we'll discuss tomorrow and that are
important. Some are intrinsic and some are just
266
00:22:09.760 --> 00:22:16.920
related to our imperfect knowledge of what should
be the exact exchange-correlation functional.
267
00:22:17.760 --> 00:22:21.980
OK, with this, we take for a moment a pause
268
00:22:22.080 --> 00:22:28.740
from the concepts and the tale of density
functional theory and we delve in,
269
00:22:28.840 --> 00:22:37.960
into this general issue, on the fact that we
have now created our Kohn-Sham equations.
270
00:22:38.240 --> 00:22:43.960
These do not look very different from
Hartree or Hartree-Fock so the only difference
271
00:22:44.060 --> 00:22:50.840
between these three approaches is having
either the LDA or GGA or whatever
272
00:22:50.940 --> 00:22:56.040
exchange-correlation functional here
or the Hartree in the exchange term here.
273
00:22:56.140 --> 00:23:02.580
Sorry, the Hartree is always there also in density
functional theory or just the Hartree equation.
274
00:23:02.680 --> 00:23:08.200
But in general, we have to solve
these differential equations.
275
00:23:08.680 --> 00:23:17.160
In addition, you know, all of these
contain the solution so the orbitals that solve
276
00:23:17.260 --> 00:23:21.140
this differential equation,
you know, build up the charge density,
277
00:23:21.240 --> 00:23:27.140
that is say, a constituent of the Hartree
operator in the Hartree energy,
278
00:23:27.240 --> 00:23:32.240
and so we need from the solutions to build
the operators, find the solution,
279
00:23:32.340 --> 00:23:38.760
build the operators and self-consistently
converge this cycle.
280
00:23:39.760 --> 00:23:43.340
How do we solve a differential equation?
281
00:23:43.440 --> 00:23:46.180
Well, if this was the 19th century,
282
00:23:46.280 --> 00:23:50.120
you would probably learn what
the ideal geometric function is.
283
00:23:50.220 --> 00:23:53.580
You would learn about Laguerre polynomials
284
00:23:53.680 --> 00:23:58.240
and, you know, all the great
accomplishment of analysis.
285
00:23:58.480 --> 00:24:03.300
We are not going to do that also because the
system we are going to solve are very complex.
286
00:24:03.400 --> 00:24:05.200
So we are going to use a computer.
287
00:24:05.300 --> 00:24:12.220
So we are going to turn a differential
equation into something that is computable.
288
00:24:12.320 --> 00:24:14.420
How do we do that?
289
00:24:14.520 --> 00:24:18.380
By actually transforming
290
00:24:18.480 --> 00:24:24.920
our unknown object. You know, in the differential
equation we need to find a solution,
291
00:24:25.020 --> 00:24:31.780
we need to find a wavefunction
into algebra or linear algebra problem.
292
00:24:31.880 --> 00:24:34.820
That is, rather than finding what is
293
00:24:34.920 --> 00:24:40.460
the unknown psi, we are going
to find what are the coefficients.
294
00:24:40.560 --> 00:24:45.280
These are real or complex numbers of this
295
00:24:45.640 --> 00:24:52.920
cm that represent the expansion of the unknown
296
00:24:53.560 --> 00:25:00.700
wavefunction psi into a basis of well-known
typically analytical function phi.
297
00:25:00.800 --> 00:25:03.280
So give me ... let me give you an example.
298
00:25:03.380 --> 00:25:08.399
Suppose, you know, our ultimate goal
was to solve a differential equation
299
00:25:08.500 --> 00:25:13.060
that was giving me as a solution,
something that looked like this.
300
00:25:13.160 --> 00:25:15.080
You know, maybe this looks like a Gaussian
301
00:25:15.180 --> 00:25:18.820
or maybe looks almost like
a Gaussian, but not exactly.
302
00:25:18.920 --> 00:25:24.380
Maybe, let's say for simplicity, we are in
periodic boundary conditions so here my
303
00:25:24.480 --> 00:25:29.900
periodic boundary drawing skill are quite
pathetic, but this thing is periodic.
304
00:25:30.000 --> 00:25:36.260
How can we, you know, transform
this differential problem?
305
00:25:36.360 --> 00:25:45.460
Well, we could expand with Fourier analysis
this unknown red function psi
306
00:25:45.560 --> 00:25:52.140
into a linear combination of sines and cosines
or actually probably because this looks
307
00:25:52.240 --> 00:25:57.700
as if it were symmetric around zero,
probably, we just need only cosines.
308
00:25:57.800 --> 00:25:59.160
And here is an example.
309
00:25:59.160 --> 00:26:00.480
You see a...
310
00:26:00.760 --> 00:26:05.120
plotted cosines of different wavelength.
311
00:26:05.220 --> 00:26:07.700
This is the one with the
312
00:26:07.800 --> 00:26:17.320
largest wavelength and then we go to finer
and finer wavelength and each of these...
313
00:26:17.960 --> 00:26:23.540
cosines... so sorry, I'm not drawing very well here
has a certain coefficient c1, c2, c3,
314
00:26:23.640 --> 00:26:30.340
and so maybe if 10 cosines are good
enough to describe psi very well,
315
00:26:30.440 --> 00:26:38.660
what I need to solve my problem is
actually finding those 10 coefficients.
316
00:26:38.760 --> 00:26:43.460
And you see, those 10 coefficients need to be such that
317
00:26:43.560 --> 00:26:46.700
around zero they, all the cosines,
318
00:26:46.800 --> 00:26:53.840
sum up constructively and the zone boundary is
where -8 is, they sum up destructively.
319
00:26:53.940 --> 00:26:59.320
So you have a maximum,
then flat and then a maximum again.
320
00:27:00.160 --> 00:27:05.680
So let me actually show how formally this allows
321
00:27:05.780 --> 00:27:12.940
us to solve the Schrodinger equation in our case,
but it could be any other differential equation
322
00:27:13.040 --> 00:27:17.880
or it could be the Kohn-Sham equation
or it could be the Hartree-Fock equation
323
00:27:18.000 --> 00:27:21.300
whatever this Hamiltonian operator is.
324
00:27:21.400 --> 00:27:26.100
And I'm using again the Dirac notation so this
325
00:27:26.200 --> 00:27:34.100
ket here in real space would just be
a wavefunction psi of r.
326
00:27:34.200 --> 00:27:37.660
And so the core
327
00:27:37.760 --> 00:27:46.100
answers, the core hypothesis that
I make is that I expand these psi
328
00:27:46.329 --> 00:27:51.879
as a linear combination of basis elements
329
00:27:51.880 --> 00:27:59.780
and these basis elements I choose them
to be orthogonal to make my life easier.
330
00:27:59.880 --> 00:28:06.040
So what I say is that the scalar product with..
331
00:28:06.720 --> 00:28:11.180
between, say, the basis element
n and the basis element m
332
00:28:11.280 --> 00:28:16.743
is going to be equal to 1 or 0 depending
if n is equal or different from m.
333
00:28:16.744 --> 00:28:20.200
In, you know, our standard first quantization picture
334
00:28:20.300 --> 00:28:27.740
this would be the integral of fm star times fn of r.
335
00:28:27.840 --> 00:28:31.620
And of course we understand right away
that if n is equal to m,
336
00:28:31.720 --> 00:28:37.520
this is equal to one from the normalisation
condition and if they are different...
337
00:28:37.620 --> 00:28:44.440
well I mean posing that the scalar product is
different then you can see that it does say if you
338
00:28:44.540 --> 00:28:52.460
take sine of 3 x and cosine of 3 x, you see
right away that they are orthogonal by symmetry.
339
00:28:52.560 --> 00:28:58.100
We don't have to have this orthogonality. We can
work with a basis set that is non-orthogonal.
340
00:28:58.200 --> 00:29:02.160
It's just to simplify the algebra here.
341
00:29:03.400 --> 00:29:07.140
So how do I make progress here?
342
00:29:07.240 --> 00:29:12.860
If you for a moment, are familiar also
with the concept of Hilbert space
343
00:29:12.960 --> 00:29:17.140
and the psi being a vector in Hilbert space
344
00:29:17.240 --> 00:29:22.420
you understand that H psi
is also a vector in Hilbert space
345
00:29:22.520 --> 00:29:30.360
and what I'm going to do is I'm going to project
the left-hand side of this equation H psi and
346
00:29:30.460 --> 00:29:39.640
the right-hand side of this equation, a
constant times psi onto phi 1, phi 2, phi 3
347
00:29:39.760 --> 00:29:47.140
so if H psi was a vector in cartesian space,
it would have three coordinates,
348
00:29:47.240 --> 00:29:53.220
say the vector could be one to seven,
and these three coordinates are the projection
349
00:29:53.320 --> 00:30:01.520
of the vector on to the three different
versos 100, 010, 001
350
00:30:01.840 --> 00:30:06.920
that constitute the basis of my
three-dimensional cartesian space.
351
00:30:07.020 --> 00:30:12.500
Here, my Hilbert space has many dimension,
in principle, it could have infinite dimension.
352
00:30:12.600 --> 00:30:17.180
I taking a basis set of orthogonal normalised
353
00:30:17.280 --> 00:30:25.740
versos so there's phi n of the equivalent
of the versos of Cartesian space and I project
354
00:30:25.840 --> 00:30:32.140
my Schrodinger or my Hartree-Fock or my
Kohn-Sham equation onto each of these verso.
355
00:30:32.240 --> 00:30:39.680
So if for simplicity, I suppose that I have
ten elements in my basis set - it could be one
356
00:30:39.780 --> 00:30:43.780
thousand, it could be one million,
but let's suppose there are only 10,
357
00:30:43.880 --> 00:30:52.660
I have to project my differential equation onto
those 10 versos so m goes from 1 to 10.
358
00:30:52.760 --> 00:30:55.980
So now I don't have any more
359
00:30:56.080 --> 00:31:02.540
vectorial identity as it was here, but here
I have scalars - I have projected this.
360
00:31:02.640 --> 00:31:09.860
So you see I have scalars on the left and
scalars on the right and I have
361
00:31:09.960 --> 00:31:17.720
right now out of this projection 10 different
condition for phi 1, phi 2, phi 3 ... phi 10.
362
00:31:18.120 --> 00:31:21.020
But I have still to exploit...
363
00:31:21.120 --> 00:31:24.820
let me change colour just for the sake of it.
364
00:31:24.920 --> 00:31:32.600
I still have to exploit this
ansatz that is I have expanded psi as a linear
365
00:31:32.700 --> 00:31:36.140
combination of my basis set
with unknown coefficients.
366
00:31:36.240 --> 00:31:44.220
So I'm going to put for psi this linear
sum and so the net result is here.
367
00:31:44.320 --> 00:31:49.760
You see, when I calculate the scalar product
368
00:31:50.720 --> 00:31:58.060
phi m psi this is going to be equal to
sum over n c n phi m phi n,
369
00:31:58.160 --> 00:32:04.620
but phi m phi n because of orthogonality,
that's the simplification, is going to be equal to
370
00:32:04.720 --> 00:32:14.160
Kronecker's Delta and so this is going
to be equal to cm and there is my cn.
371
00:32:14.760 --> 00:32:18.214
I can't do the same trick on the left-hand side,
372
00:32:18.215 --> 00:32:24.971
so you have to keep all my complexity
so I have the sum for n that goes from 1 to 10
373
00:32:24.972 --> 00:32:29.020
of these unknown coefficients cn times...
374
00:32:29.120 --> 00:32:34.680
this is, again, a scalar - is the expectation
value of the Hamiltonian, could be Schrodinger,
375
00:32:34.780 --> 00:32:40.260
could be Kohn-Sham, could be Hartree-Fock
between these two basis elements.
376
00:32:40.360 --> 00:32:48.360
And this is actually a matrix
element that we could write as Hmn.
377
00:32:48.560 --> 00:32:51.220
So in shorthand form,
378
00:32:51.320 --> 00:32:58.686
my 10 equations, so remember, m goes from 1 to 10,
379
00:32:58.743 --> 00:33:02.043
my 10 equation are represented here
380
00:33:02.100 --> 00:33:12.614
where the matrix Hmn has 100 elements because
m can go from 1 to 10 and n can go from 1 to 10,
381
00:33:12.914 --> 00:33:18.357
and this is the expectation value of phi m H phi n.
382
00:33:19.560 --> 00:33:25.860
So we have 10 equations that
are a matrix. The matrix is known, Hmn,
383
00:33:25.960 --> 00:33:31.760
and that's why it's good to have a simple
analytical function in those phi m and phi n.
384
00:33:31.860 --> 00:33:35.500
And because when I apply H to phi n,
385
00:33:35.600 --> 00:33:39.220
I want to be able to calculate
the Laplacian, the second derivatives.
386
00:33:39.320 --> 00:33:40.700
I want something analytical.
387
00:33:40.800 --> 00:33:46.180
So, if I've chosen my basis set I can calculate
388
00:33:46.280 --> 00:33:52.700
my matrix elements and now I
have to solve this 10 matrix..
389
00:33:52.800 --> 00:33:54.680
these 10...
390
00:33:55.320 --> 00:34:00.220
these 10 linear equations that are really written
391
00:34:00.320 --> 00:34:06.000
in a much simpler form,
as a matrix vector multiplication.
392
00:34:06.100 --> 00:34:11.520
You see, if I look at the
multiplication of the first row
393
00:34:12.240 --> 00:34:13.440
H11 to H1,10
394
00:34:13.720 --> 00:34:23.440
times c1, c10 being equal to E c1. This is
exactly the first equation for m equal to one
395
00:34:23.680 --> 00:34:25.280
and then for m equal to two,
396
00:34:25.380 --> 00:34:28.540
i have this equation that is the product of this
397
00:34:28.640 --> 00:34:33.300
row time the same column equal to E c2
and so on and so forth.
398
00:34:33.400 --> 00:34:41.500
So you see how I've taken a differential equation,
my only hypothesis was choosing a basis set,
399
00:34:41.600 --> 00:34:46.640
the unknown have become the coefficients
of the expansion of my unknown wavefunction
400
00:34:46.740 --> 00:34:50.860
into this basis set, and it's become
a linear algebra problem.
401
00:34:50.960 --> 00:34:53.780
This is something that a computer understands.
402
00:34:53.880 --> 00:34:58.900
You have to do matrix times vector
equal to a constant times a vector.
403
00:34:59.000 --> 00:35:00.320
I can actually bring...
404
00:35:00.420 --> 00:35:01.700
sorry, let me go back.
405
00:35:01.800 --> 00:35:04.820
I can bring the vector
406
00:35:04.920 --> 00:35:08.900
on the left-hand side so I can make this
407
00:35:09.000 --> 00:35:17.860
equal to zero, and what I have here is
that on the diagonal, I have H22 minus epsilon
408
00:35:17.960 --> 00:35:22.660
H33, sorry minus E, down to Hkk, minus E.
.
409
00:35:22.760 --> 00:35:27.520
So now I have even a simpler
equation in linear algebra.
410
00:35:27.620 --> 00:35:31.420
I have matrix time vector equal to zero.
411
00:35:31.520 --> 00:35:37.840
And you know from linear algebra that the only
case where this has a non-trivial solution
412
00:35:37.940 --> 00:35:41.560
for the coefficients - you know,
if the 10 coefficients are zero,
413
00:35:41.660 --> 00:35:46.820
OK, so the problem that is not interesting.
The only time when I have a non-trivial solution
414
00:35:46.920 --> 00:35:55.520
is when the determinant of the matrix is
equal to zero - I wrote the determinant as this.
415
00:35:56.000 --> 00:36:04.420
And, you know, a determinant here, the unknown is E.
This is the eigenvalue of my equation.
416
00:36:04.520 --> 00:36:13.620
I know the 10x10 Hamiltonian matrix represented
on that basis and this is going to be a polynomial
417
00:36:13.720 --> 00:36:19.360
of tenth degree for this specific
case and some degree, in general.
418
00:36:19.460 --> 00:36:23.060
So it has n solutions, let's say 10 solution.
419
00:36:23.160 --> 00:36:27.180
Those 10 solutions will be the 10 eigenvalues
420
00:36:27.280 --> 00:36:33.280
of my problem and then for each eigenvalue
I can go back and I can solve...
421
00:36:33.380 --> 00:36:38.300
find the non-trivial set of coefficients c1, ck
422
00:36:38.400 --> 00:36:44.100
that give me the eigenvector that I'm looking for.
423
00:36:44.200 --> 00:36:47.080
So this is how we do differential equation on
424
00:36:47.180 --> 00:36:51.860
a computer and you don't have to be interested
in quantum mechanics to know these things.
425
00:36:51.960 --> 00:36:55.040
You know, differential equation are
ubiquitous in science and engineering.
426
00:36:55.140 --> 00:36:57.400
If you want to make your plane fly,
427
00:36:57.400 --> 00:37:00.480
you probably have to solve to begin
with some differential equations.
428
00:37:00.580 --> 00:37:04.240
So this is a very general thing.
429
00:37:04.880 --> 00:37:08.420
It's not the only approach that we have to solve
430
00:37:08.520 --> 00:37:15.420
the problem, because we could have
gone back and thought not at the
431
00:37:15.520 --> 00:37:19.620
Euler-Lagrange equation,
that were related to our variational functional
432
00:37:19.720 --> 00:37:25.380
but actually we could have
attacked the functional directly.
433
00:37:25.480 --> 00:37:31.440
Either rewritten here for you, the
Kohn-Sham energy functional.
434
00:37:31.440 --> 00:37:33.740
You see it's a functional of the charge density.
435
00:37:33.840 --> 00:37:42.120
It's a function of psi i and I'm going now
to expand the psi i in this basis set.
436
00:37:42.220 --> 00:37:46.800
So my psi i of r is written as a linear
437
00:37:46.900 --> 00:37:51.860
combination with certain coefficients
of my basis set phi n.
438
00:37:51.960 --> 00:37:54.840
So you have the n index that sum over all
439
00:37:54.940 --> 00:38:03.440
the possible basis set and the i index that
index the particular electronic states that I want
440
00:38:03.540 --> 00:38:09.400
to consider and if I have 10 electrons,
you see I have 10 electronic states.
441
00:38:10.160 --> 00:38:17.160
Well, now, again, I can put this
expression for each of my
442
00:38:17.800 --> 00:38:25.020
psi i into this energy functional
and I can work out all those integrals
443
00:38:25.120 --> 00:38:34.980
and so the unknown remain my coefficients so
that I want to find for every possible psi i.
444
00:38:35.080 --> 00:38:43.260
So I will have here c11, c12 up to c1k and then
445
00:38:43.360 --> 00:38:46.920
depending on how many electrons I have cn1
446
00:38:47.800 --> 00:38:50.220
up to cnk.
447
00:38:50.320 --> 00:38:53.860
So a lot of coefficients that our computers
have to deal with,
448
00:38:53.960 --> 00:39:02.500
but finding the solution to my problem can also be
stated in this equivalent form in which I have to
449
00:39:02.600 --> 00:39:10.040
minimise the functional,
the energy functional with respect to
450
00:39:10.160 --> 00:39:12.300
all my coefficients,
451
00:39:12.400 --> 00:39:19.300
so I can see the problem of getting to the ground
state, as a non-linear minimisation problem.
452
00:39:19.400 --> 00:39:22.900
We see examples at the end of today.
453
00:39:23.000 --> 00:39:27.600
For the time being, we actually
think that solving the differential
454
00:39:27.700 --> 00:39:34.160
equation as a linear algebra problem,
as inversion of that matrix.
455
00:39:34.320 --> 00:39:39.260
Now, we have spoken about basic sets, in general,
456
00:39:39.360 --> 00:39:45.860
but let's actually focus in our actual
calculation what would be an ideal basis set.
457
00:39:45.960 --> 00:39:48.620
And if you are a quantum chemist,
458
00:39:48.720 --> 00:39:55.900
often you are dealing with molecules so it's
convenient because you don't want to describe
459
00:39:56.000 --> 00:39:59.360
wavefunction in the empty space around the molecule.
460
00:39:59.460 --> 00:40:05.160
You want to describe it only where the molecule
is, to have maybe localised basis set.
461
00:40:05.260 --> 00:40:07.580
It could be atomic orbitals.
462
00:40:07.680 --> 00:40:13.280
You might want to decide if you want your
basis set to be orthogonal or not orthogonal,
463
00:40:13.380 --> 00:40:19.580
and then you have to keep track of the
non-orthogonality in terms of overlap matrix.
464
00:40:19.680 --> 00:40:24.514
It's very common to use a basis set
based on Gaussians because then
465
00:40:24.529 --> 00:40:31.360
the Hartree-Fock exchange integrals
becomes very inexpensive.
466
00:40:31.480 --> 00:40:35.660
We are going to deal with extended system.
467
00:40:35.760 --> 00:40:41.400
Those could be perfectly crystalline systems like
a piece of silicon that can be described very well
468
00:40:41.500 --> 00:40:47.300
by two atoms periodically depleted,
but it could be also, say, an amorphous system
469
00:40:47.400 --> 00:40:56.220
that is solid but is not periodic or a liquid
system that again, it's extended but not periodic.
470
00:40:56.320 --> 00:40:58.280
And even in all those cases,
471
00:40:58.380 --> 00:41:08.680
suppose we want to study, say, water, liquid water
as it exists in my water glass [drinks water]
472
00:41:10.120 --> 00:41:17.820
it's much better rather than to describe a
droplet of water with all the water molecules
473
00:41:17.920 --> 00:41:22.820
in and then taking the limit
of the water droplet being very large
474
00:41:22.920 --> 00:41:28.800
so that the surface effect that go as the square
of the radius become negligible with respect
475
00:41:28.900 --> 00:41:33.220
to the volume effect
that go as the cube of the radius.
476
00:41:33.320 --> 00:41:36.640
It's much better to actually get to the limit
477
00:41:36.740 --> 00:41:41.700
of the infinite droplet,
not by studying a finite droplet that becomes
478
00:41:41.800 --> 00:41:48.740
bigger and bigger, but to study
a system that is periodic.
479
00:41:48.840 --> 00:41:52.560
A periodic system has the beauty
that is both infinite and infinite.
480
00:41:52.660 --> 00:41:56.780
So you see, if I put 5 water molecules
481
00:41:56.880 --> 00:42:03.420
in my simulation cell, that maybe it's cubic
and I repeat them periodically,
482
00:42:03.520 --> 00:42:06.880
what I have is that my system is infinite
483
00:42:06.980 --> 00:42:13.580
in the sense that there are no surfaces so I
don't have to worry about the surface effects.
484
00:42:13.680 --> 00:42:20.820
But it's also finite in the sense that I
have only, say, 15 independent atoms...
485
00:42:20.920 --> 00:42:26.900
Sorry, yeah, 15 independent atoms,
45 degrees of freedom in the ionic coordinates.
486
00:42:27.000 --> 00:42:34.460
But what I can do is then rather than study five
water molecules, I can study ten water molecules
487
00:42:34.560 --> 00:42:39.100
periodically repeated, and then
I can study 20, 30 and so on.
488
00:42:39.200 --> 00:42:44.900
And what happens is that because I've never
have surfaces in this,
489
00:42:45.000 --> 00:42:50.840
I get to the limit where the properties
of my system don't change any more.
490
00:42:50.860 --> 00:42:53.160
I reach what is called the thermodynamic limit
491
00:42:53.800 --> 00:42:59.420
much, much faster if I were to study
the increasing droplet.
492
00:42:59.520 --> 00:43:04.380
So studying, you know, systems in periodic
boundary conditions
493
00:43:04.480 --> 00:43:09.780
is very convenient to study anything that is
extended. A solid for which is natural,
494
00:43:09.880 --> 00:43:15.320
crystalline solid, but also an amorphous
solid glass or even a liquid.
495
00:43:15.520 --> 00:43:17.720
If you study molecules, it's not as good,
496
00:43:17.820 --> 00:43:25.900
but you can study still it very easily by,
you know, putting the molecules in a unit cell
497
00:43:26.000 --> 00:43:33.400
that is large enough so that it's
periodic image doesn't feel
498
00:43:33.520 --> 00:43:36.360
the effects of the original one.
499
00:43:36.460 --> 00:43:45.060
So if these 2 benzene molecules are
far away so that their wavefunctions don't overlap
500
00:43:45.160 --> 00:43:50.743
and they don't have, you know, relevant dipoles
or quadrupoles with which they interact
501
00:43:50.757 --> 00:43:54.740
electrostatically, molecule A
502
00:43:54.840 --> 00:43:59.040
will think that it's the only molecule
in the universe because it
503
00:43:59.140 --> 00:44:04.360
doesn't feel molecule B, and so
even if it's studied in periodic boundary
504
00:44:04.460 --> 00:44:10.020
condition, that means every time we move this
hydrogen atom, also this hydrogen atom moves.
505
00:44:10.120 --> 00:44:16.500
It doesn't feel the effect so it
thinks that it's alone in the universe.
506
00:44:16.600 --> 00:44:20.480
So we'll use periodic boundary conditions to study
507
00:44:20.580 --> 00:44:27.600
solids and for that, I want to give you a quick
refresher of periodicity and the translational
508
00:44:27.700 --> 00:44:33.780
operations in three dimensions because this will
be actually the language with which we study
509
00:44:33.880 --> 00:44:41.340
periodic systems like silica, or we set up our
calculation even for water or for a glass.
510
00:44:41.440 --> 00:44:43.420
And you are probably familiar with this concept
511
00:44:43.520 --> 00:44:44.743
of a Bravais lattice.
512
00:44:44.771 --> 00:44:47.900
A Bravais lattice is a mathematical concept.
513
00:44:47.901 --> 00:44:51.740
If you want, it's embodied in having
514
00:44:51.840 --> 00:45:00.980
three primitive lattice vectors that we call a1,
a2 and a3. That they do not lie in the same plane,
515
00:45:01.080 --> 00:45:11.660
and we look at the infinite set of vectors r
that are spanned by the linear combination of
516
00:45:11.760 --> 00:45:18.200
a1, a2 and a3 with the integral
coefficients l, m and n. l can be...
517
00:45:19.920 --> 00:45:28.600
They belong to Z, so they can be
positive integers, negative integers or zero.
518
00:45:28.880 --> 00:45:30.600
And Bravais,
519
00:45:30.643 --> 00:45:36.140
you know, now we remember his name because
he was the first in 1848 to classify
520
00:45:36.240 --> 00:45:42.520
correctly the fact that in three dimensions there
are 14 qualitatively different Bravais lattices.
521
00:45:42.680 --> 00:45:45.700
Someone before him thought that there were 15.
522
00:45:45.800 --> 00:45:51.500
Everyone has forgotten his name, poor thing.
I will not mention it today.
523
00:45:51.600 --> 00:45:54.140
So these are the 14 Bravais lattices.
524
00:45:54.240 --> 00:45:55.840
I will not go into the details.
525
00:45:55.940 --> 00:46:03.260
We divided them in seven crystal classes, divided
them in simple volume, body or face-centered.
526
00:46:03.360 --> 00:46:05.800
Maybe in the case of a cubic system it's worth
527
00:46:05.900 --> 00:46:09.700
remembering we have the simple cubic,
we have the body-centred cubic
528
00:46:09.800 --> 00:46:12.580
and we have the face-centred cubic.
529
00:46:12.680 --> 00:46:19.060
And if we wanted to represent, say,
an infinite crystal of silicon.
530
00:46:19.160 --> 00:46:23.686
Silicon, the crystal, is a diamond lattice
that means it's made by
531
00:46:23.687 --> 00:46:33.080
two silicon atoms repeated periodically
with the translation operation of the
532
00:46:33.240 --> 00:46:33.980
Bravais cubic face-centred lattice.
533
00:46:34.080 --> 00:46:42.500
So a real physical object like a silicon crystal
is made by a physical group of atoms.
534
00:46:42.600 --> 00:46:47.780
We also call that a basis - in the case
of silicon, two atoms
535
00:46:47.880 --> 00:46:52.460
periodically repeated with the mathematical
translation of the Bravais lattice.
536
00:46:52.560 --> 00:46:56.780
If you take a piece of ferromagnetic iron
537
00:46:56.880 --> 00:47:04.214
that can be represented by one
iron atom periodically repeated
538
00:47:04.215 --> 00:47:10.171
using the volume centred, the BCC cubic lattice.
539
00:47:11.000 --> 00:47:13.680
Sadly, you know, this definition is not unique.
540
00:47:13.780 --> 00:47:23.020
I can represent silicon not by taking, say,
two atoms and repeating them periodically
541
00:47:23.120 --> 00:47:32.000
with the FCC operation, but I could take
these four atoms that I am painting in green right
542
00:47:32.100 --> 00:47:36.680
now and I could repeat them
with the simple cubic operations.
543
00:47:36.780 --> 00:47:41.120
So in some ways I've created a super
cell that contains four atoms.
544
00:47:41.220 --> 00:47:43.460
I've changed my translation operation.
545
00:47:43.560 --> 00:47:51.680
The net result is an infinite system
that has been described equally well.
546
00:47:52.920 --> 00:48:01.080
OK, so with this, let me also introduce
the concept of the reciprocal lattice.
547
00:48:01.920 --> 00:48:08.000
Here we have a direct lattice that is a simple
cubic - apologies, the symbols are not coming out
548
00:48:08.100 --> 00:48:14.400
well-described, in general,
by the linear combination of a1, a2 and a3.
549
00:48:14.960 --> 00:48:20.240
These solid spheres could
be the Bravais lattice points
550
00:48:21.160 --> 00:48:26.000
or they could be, you know, real physical atoms
551
00:48:26.280 --> 00:48:31.620
but what I am going to define is
the reciprocal lattice
552
00:48:31.720 --> 00:48:38.360
and for that I'll introduce the concept
that may be are familiar with, of a plane-wave.
553
00:48:38.360 --> 00:48:39.460
So,
554
00:48:39.560 --> 00:48:49.100
in real space it's the exponential of the imaginary unit
i times a vector G, a wave vector G, times r.
555
00:48:49.200 --> 00:48:52.740
So what is the characteristic of a plane-wave?
556
00:48:52.840 --> 00:49:00.620
Well, it's modulation is going to be
along the direction of the G wave vector.
557
00:49:00.720 --> 00:49:06.580
This is because any point that lies in a plane
558
00:49:06.680 --> 00:49:11.000
perpendicular to G will have the same scalar
559
00:49:11.100 --> 00:49:16.940
product between any point - this point,
this point, this point and G, so the scalar
560
00:49:17.040 --> 00:49:23.740
product is the same so the value of the amplitude
of the plane-wave is the same at every point.
561
00:49:23.840 --> 00:49:28.420
And as we move along this direction,
562
00:49:28.520 --> 00:49:34.340
we are going to have a real and an imaginary
part that are both, you know,
563
00:49:34.440 --> 00:49:44.020
cosinusoidal and sinusoidal and
the wavelength is going to be very large
564
00:49:44.120 --> 00:49:49.200
if the modulus of G is very small,
because if the modulus of G is very small,
565
00:49:49.300 --> 00:49:58.060
I have to move a long way with r
before G dot r has changed by amount...
566
00:49:58.160 --> 00:50:02.500
by an amount of pi that would be
the phase inversion in the plane-wave.
567
00:50:02.600 --> 00:50:06.140
So G is small, wavelength is long.
568
00:50:06.240 --> 00:50:10.540
G is large, wavelength is short.
569
00:50:10.640 --> 00:50:13.080
And we are going to define the ...
570
00:50:13.180 --> 00:50:17.960
a set of G vectors, such that
571
00:50:18.800 --> 00:50:28.820
the plane-wave i G dot r is, we say, compatible
with my primitive direct lattice.
572
00:50:28.920 --> 00:50:32.840
So you see, if I have a plane-wave
that looks like this,
573
00:50:33.080 --> 00:50:39.720
this is compatible because it has the same
amplitude at every point of the Bravais lattice.
574
00:50:39.720 --> 00:50:41.600
That amplitude doesn't have to be zero.
575
00:50:41.700 --> 00:50:45.060
I mean, it could have been
this plane-wave shifted.
576
00:50:45.160 --> 00:50:50.060
It still has the same amplitude at all
the point of the Bravais lattice.
577
00:50:50.160 --> 00:50:58.240
But if I were to take for a moment
this orange plane-wave here,
578
00:51:00.640 --> 00:51:05.660
this is not compatible because it
doesn't have the same amplitude.
579
00:51:05.760 --> 00:51:10.260
So if my Bravais lattice was one-dimensional
580
00:51:10.360 --> 00:51:15.860
you understand, right away, that the
plane-waves - I'll do the real part
581
00:51:15.960 --> 00:51:19.040
or the imaginary part sines and cosines that are
582
00:51:19.140 --> 00:51:22.600
compatible are these, you know,
something like the fundamental
583
00:51:22.700 --> 00:51:32.080
harmonics, and then what I can have is, you know,
two wavelengths and then it becomes difficult.
584
00:51:32.180 --> 00:51:36.020
I could have three wavelengths
and any higher wavelength.
585
00:51:36.120 --> 00:51:39.700
So how do I define mathematically this condition
586
00:51:39.800 --> 00:51:45.700
if I want the amplitude to be
the same at every point?
587
00:51:45.800 --> 00:51:48.360
What I want is this equality here.
588
00:51:48.460 --> 00:51:52.940
That is, I want my plane-wave
to have the same amplitude
589
00:51:53.040 --> 00:52:00.700
at an arbitrary point r and I translate that point
by a Bravais lattice vector, capital R.
590
00:52:00.800 --> 00:52:05.060
So if I work out the algebra for a generic
591
00:52:05.160 --> 00:52:11.500
capital R, I want this to be true,
and you see it right away,
592
00:52:11.600 --> 00:52:18.360
that this is going to be true if G has the form
that is in itself a Bravais lattice form,
593
00:52:18.460 --> 00:52:24.560
as a linear combination of three
primitive lattice vectors: b1, b2 and b3,
594
00:52:27.480 --> 00:52:36.060
where b1, b2 and b3 have the property
that they are dual to the direct lattice vector.
595
00:52:36.160 --> 00:52:45.500
So bi scalar product aj - I have, you know,
three b vectors, I have three a vectors
596
00:52:45.501 --> 00:52:51.457
and the scalar product between each one
of the three b and each one of the three j
597
00:52:51.458 --> 00:52:57.080
needs to give me 0 or 1, depending
if he is equal or different from j.
598
00:52:57.180 --> 00:53:00.880
This is the property of a dual and you can work
599
00:53:00.980 --> 00:53:06.240
out the algebra that if I have
a1, a2 and a3, b1, b2 and b3
600
00:53:06.640 --> 00:53:09.260
are defined by this cross product,
601
00:53:09.360 --> 00:53:17.280
while at the denominator here this is nothing
else than the volume of the unit cell.
602
00:53:17.640 --> 00:53:21.620
OK, so this was
an introduction to Bravais lattice.
603
00:53:21.720 --> 00:53:25.880
Why is it important?
Because we need it to construct
604
00:53:26.800 --> 00:53:34.000
our input file for our simulations, because
we might want to study silicon and the simplest
605
00:53:34.100 --> 00:53:40.480
way to study silicon is saying that it
is given by an FCC Bravais lattice.
606
00:53:41.360 --> 00:53:45.480
So in the language of Quantum Espresso ibrav is
607
00:53:45.580 --> 00:53:53.380
equal to 2 and then given this translational
symmetry, you still need to specify that you have
608
00:53:53.480 --> 00:54:00.740
two atoms that would be at a certain
distance from each other.
609
00:54:00.840 --> 00:54:05.900
And in particular, if
this is the conventional unit cell, a,
610
00:54:06.000 --> 00:54:11.720
you could say put one atom in 000,
you could put the other in 1/4, 1/4, 1/4 a
611
00:54:14.040 --> 00:54:16.460
and then you could repeat all of this
612
00:54:16.560 --> 00:54:23.500
with the FCC lattice translations that are
given by these three Bravais lattice vectors.
613
00:54:23.600 --> 00:54:30.640
So that's the importance of Bravais lattices
for Quantum Espresso.
614
00:54:31.840 --> 00:54:34.160
In general,
615
00:54:35.160 --> 00:54:41.620
Bravais lattices have also relevance on how
we classify electronic states in a system.
616
00:54:41.720 --> 00:54:46.180
You know, remember for a moment
that the charged density, say,
617
00:54:46.280 --> 00:54:55.940
of density functional theory is going to be given
by the sum over all the occupied states
618
00:54:56.040 --> 00:55:02.380
of the square modulus of each one
of these Kohn-Sham electrons.
619
00:55:02.480 --> 00:55:09.840
So in general, what we have often
to do is to take sum over all
620
00:55:09.940 --> 00:55:16.940
the occupied states to get physical quantities
like the charge density or the total energy.
621
00:55:17.040 --> 00:55:19.000
In a solid the things are tricky because
622
00:55:19.100 --> 00:55:27.420
the system is infinite, so we can do
infinite sum and we'll do appropriate
623
00:55:27.520 --> 00:55:32.740
coarse graining of this infinite sum
using reciprocal space concepts.
624
00:55:32.840 --> 00:55:35.540
Let me start reminding you, you know,
625
00:55:35.640 --> 00:55:43.860
one of the fundamentals theorem of condensed
matter physics, the fact that if the Hamiltonian
626
00:55:43.960 --> 00:55:51.500
commutes with any translation given by a lattice...
627
00:55:51.600 --> 00:55:56.520
Bravais lattice vector R, OK,
so if applying the Hamiltonian first
628
00:55:56.620 --> 00:56:02.500
and the lattice translation afterwards or
the lattice translation first and the Hamiltonian
629
00:56:02.600 --> 00:56:08.060
later gives you the same result
if these two operators commute,
630
00:56:08.160 --> 00:56:15.740
what you can conclude, I will not demonstrate this
for you is that the eigenstates of the Hamiltonian
631
00:56:15.840 --> 00:56:21.780
can be expressed in a form
where they have also very well-defined
632
00:56:21.880 --> 00:56:25.340
symmetry properties with respect to the translation.
633
00:56:25.440 --> 00:56:29.020
Sadly, and interestingly these symmetry properties
634
00:56:29.120 --> 00:56:34.900
are not that the eigenstates
of the Hamiltononian are periodic.
635
00:56:35.240 --> 00:56:37.160
The charge density is periodic,
636
00:56:37.260 --> 00:56:43.280
it's a physical quantity but the eigenstates are
not periodic the eigenstates of a periodic
637
00:56:43.380 --> 00:56:51.540
Hamiltonian have the Bloch form, that is,
they are written in terms of a periodic function
638
00:56:51.640 --> 00:57:02.100
times, again, a plane-wave where the wave vector
of the plane-wave belong to the first Brillouin zone.
639
00:57:02.200 --> 00:57:08.700
And I remind you that if we have, say,
640
00:57:08.800 --> 00:57:18.780
a reciprocal lattice so suppose that we are in two
dimensions and we had a simple square lattice in
641
00:57:18.880 --> 00:57:23.940
direct space, we have a simple
square lattice in reciprocal space.
642
00:57:24.040 --> 00:57:25.600
These are my G vectors.
643
00:57:25.700 --> 00:57:27.129
This is the origin.
644
00:57:27.130 --> 00:57:33.759
Let's call this G1, G2, G3, G4, G5.
645
00:57:33.760 --> 00:57:40.057
These are the... you see the second
nearest neighbours, G6, G7, G8.
646
00:57:40.143 --> 00:57:44.640
So these are the reciprocal lattice vectors.
647
00:57:44.740 --> 00:57:49.000
So what is the Brillouin zone is the
Wigner-Seitz cell of the reciprocal lattice.
648
00:57:49.100 --> 00:58:01.614
So what I do is I bisect all the lines that connect
a lattice point to these nearest neighbours
649
00:58:01.814 --> 00:58:09.857
and I would have also these bisections
and then I take the internal convolution,
650
00:58:10.229 --> 00:58:13.180
and the internal convolution is this.
651
00:58:13.280 --> 00:58:21.660
So this is my Brillouin zone, and so my k
belongs to this first Brillouin zone.
652
00:58:21.760 --> 00:58:26.940
So Bloch theorem tells me that
at variance with a molecule,
653
00:58:27.040 --> 00:58:31.500
where I have molecular orbitals
that I classify with integrals,
654
00:58:31.600 --> 00:58:35.940
in an infinite system where I have
infinite number of electrons,
655
00:58:36.040 --> 00:58:43.360
I have also some continuum quantum number,
and this continuum quantum number is this k
656
00:58:43.460 --> 00:58:49.140
that is called the Bloch quasi-momentum
or sometimes the Bloch momentum
657
00:58:49.240 --> 00:58:55.860
and I classify every state in a solid
with an integer number that is the band index n
658
00:58:55.960 --> 00:59:00.680
and with that continuum number in the first
Brillouin zone, that is the Bloch quasi-momentum.
659
00:59:01.480 --> 00:59:11.480
And so it is very common both in experiments, you
can measure those with arcs, or in simulations to
660
00:59:11.960 --> 00:59:19.560
plot what are the possible energy states
for any possible integer number n,
661
00:59:19.660 --> 00:59:25.220
for any possible band index n,
any possible crystal momentum or
662
00:59:25.320 --> 00:59:28.700
crystal quasi-momentum, and that
is the band structure of a solid.
663
00:59:28.800 --> 00:59:38.060
OK, so we are looking at what are
all the possible energy levels for these solids,
664
00:59:38.160 --> 00:59:47.340
and those of you that don't need to take
this webinar will recognise here a classic
665
00:59:47.440 --> 00:59:49.960
ferroelectric perovskite with a band gap,
666
00:59:50.060 --> 00:59:57.880
with 2P states here and so in orange I have
667
00:59:57.980 --> 01:00:08.020
plotted the occupied bands, maybe in blue I'm
plotting the empty bands but these bands
668
01:00:08.120 --> 01:00:13.780
are the values of energy that are
allowed by my periodic Hamiltonian.
669
01:00:13.880 --> 01:00:16.520
Could be the true one,
it could be the Hartree-Fock, it could be
670
01:00:16.620 --> 01:00:24.420
the Kohn-Sham, for every possible integer
n and for every possible Bloch momentum k.
671
01:00:24.520 --> 01:00:28.360
And Bloch's theorem tells us here,
672
01:00:28.460 --> 01:00:35.700
maybe I'm just using x variable one-dimensional,
that my wavefunction psi is not periodic
673
01:00:35.800 --> 01:00:44.000
but is the product of a periodic wavefunction u
times a plane-wave and you see the periodic
674
01:00:44.001 --> 01:00:50.900
wavefunction u is this, you know,
atomic-like function that I'm plotting here.
675
01:00:51.000 --> 01:00:54.800
So this is periodic. So this is unk
676
01:00:56.200 --> 01:00:59.200
and it will vary as a function
677
01:00:59.600 --> 01:01:06.420
of n, as a function of k, but it is periodic
and the product of unk in, let's make it, blue
678
01:01:06.520 --> 01:01:12.980
times e to the i k x
that is in red is going to give me
679
01:01:13.080 --> 01:01:20.700
the wavefunction psi that you see is
this black wavefunction.
680
01:01:20.800 --> 01:01:29.600
It's not periodic but is the product of
something periodic times our red modulation.
681
01:01:29.720 --> 01:01:35.160
OK, I think to give everyone some breathing space,
682
01:01:36.560 --> 01:01:43.080
we'll stop for five minutes,
so we'll start again at 07.
683
01:01:43.400 --> 01:01:48.960
So I'll stop for five minutes the video
and you can all take a breath.
684
01:06:30.120 --> 01:06:33.800
OK, let's... let's go back.
685
01:06:34.040 --> 01:06:38.500
So we have looked at how the wavefunction
686
01:06:38.600 --> 01:06:46.940
looks like and as I move around,
say if this is the zone centre
687
01:06:47.040 --> 01:06:48.980
and this is the zone boundary,
688
01:06:49.080 --> 01:06:53.780
what I go is from a modulation
that has a very long wavelength
689
01:06:53.880 --> 01:07:04.720
to a modulation that starts to change phase
at ... for every time I cross a lattice vector.
690
01:07:05.280 --> 01:07:06.880
So,
691
01:07:08.320 --> 01:07:16.300
the important message here is that Bloch theorem
has told me what is the shape of psi.
692
01:07:16.400 --> 01:07:23.243
It's not periodic, but is the product
of a periodic part times a phase modulation,
693
01:07:23.429 --> 01:07:31.886
and so it's this periodic part, unk,
this blue part, that we need to know
694
01:07:31.957 --> 01:07:37.240
for every band that is for every index n and for every k.
695
01:07:37.340 --> 01:07:43.860
So for every continuum momentum
in the Brillouin zone is this
696
01:07:43.960 --> 01:07:49.940
blue periodic term that is going
to be expanded in plane-waves.
697
01:07:50.040 --> 01:07:54.740
Because this is periodic the plane-waves,
698
01:07:54.840 --> 01:07:58.500
and that was the important concept
of reciprocals space,
699
01:07:58.600 --> 01:08:08.740
the G are chosen such that E to the i G dot r has
the same periodicity or higher periodicities
700
01:08:08.840 --> 01:08:14.280
compatible with my direct lattice,
and so I can expand unk
701
01:08:14.600 --> 01:08:17.580
as a linear combination of plane-waves.
702
01:08:17.680 --> 01:08:19.540
This is what your code does.
703
01:08:19.640 --> 01:08:24.900
It approximates the periodic part
of the Bloch function
704
01:08:25.000 --> 01:08:28.580
as a linear coeffici... combination of plane-waves.
705
01:08:28.680 --> 01:08:34.420
The better you want it to approximate,
the more plane-waves you need.
706
01:08:34.520 --> 01:08:40.680
If your unk in your material happens to have
a lot of structure, you are not going to be able
707
01:08:40.780 --> 01:08:46.820
to simulate it with just a few plane-waves
because you need to go to plane-waves that have
708
01:08:46.920 --> 01:08:52.820
very high periodicity to be able
to describe all that fine structures,
709
01:08:52.920 --> 01:08:59.440
and for each of those plane-waves
you'll need to know what the coefficient c is.
710
01:08:59.540 --> 01:09:07.980
The coefficient c of the plane-wave G
for that band n and for that k-point, k.
711
01:09:08.080 --> 01:09:10.480
So in a typical plane-wave calculation,
712
01:09:10.580 --> 01:09:16.780
you might have a basis set that is of the order of
thousands, tens of thousands of plane-waves,
713
01:09:16.880 --> 01:09:22.420
you can have bands that are 10, 20, 100, 1000.
714
01:09:22.520 --> 01:09:24.000
The k can be a lot.
715
01:09:24.040 --> 01:09:25.400
We'll discuss the k in a moment.
716
01:09:25.500 --> 01:09:28.320
So a lot, a lot of coefficients.
717
01:09:28.420 --> 01:09:31.860
So that's why you need computers.
718
01:09:31.960 --> 01:09:34.600
What is good and what is bad about plane-wave
719
01:09:34.700 --> 01:09:42.080
basis set, exactly as in Hi-Fi it's
very easy to improve the resolution.
720
01:09:42.180 --> 01:09:48.201
OK, so when in 1981 and 1982,
Philips came out with the CD-ROM,
721
01:09:48.401 --> 01:09:56.671
what they were doing is they were sampling
a sound so that you could have basically
722
01:09:56.672 --> 01:10:00.640
all frequencies up to around 22 kHz.
723
01:10:00.740 --> 01:10:06.027
Relying on the fact that most rich people
that were buying CDs at the time
724
01:10:06.227 --> 01:10:13.029
wouldn't have the sensitivity of, you know,
a twelve year old that can go to 25 kHz.
725
01:10:14.480 --> 01:10:16.900
Also, it was good because then you can...
726
01:10:17.000 --> 01:10:21.540
could fit the Ninth Symphony
of Beethoven in just one CD.
727
01:10:21.640 --> 01:10:23.720
And that was a major point.
728
01:10:23.820 --> 01:10:27.700
I guess, you know, it's very important
over there in Northern Europe.
729
01:10:27.800 --> 01:10:33.840
But so by taking more and more plane-waves,
you are able to describe with finer and finer
730
01:10:33.940 --> 01:10:41.140
resolution your wavefunction so you can
make it systematically better and better.
731
01:10:41.240 --> 01:10:45.180
You need a lot of them and if you were to, you know,
732
01:10:45.280 --> 01:10:51.186
need to describe, say, in a piece of gold,
both the valence electron of golds
733
01:10:51.187 --> 01:10:55.840
that are very smooth, but maybe
the 1s electron of gold that are
734
01:10:55.940 --> 01:11:04.220
extremely localised around the nuclei, you have
an untenable number of basis elements
735
01:11:04.320 --> 01:11:08.780
and what you will discover, in a moment,
is that we are not going to
736
01:11:08.880 --> 01:11:18.140
study explicitly the core electrons in a solid
that are very sharply varying in space,
737
01:11:18.240 --> 01:11:23.440
but are also not very active in forming
and creating chemical bonds.
738
01:11:23.540 --> 01:11:27.220
We are only studying the valence or so-called
739
01:11:27.320 --> 01:11:29.920
semi-valence electrons and luckily,
those are very small.
740
01:11:29.960 --> 01:11:31.200
So
741
01:11:31.400 --> 01:11:34.220
I'll show you an example in a moment that
742
01:11:34.320 --> 01:11:42.520
it's very easy to calculate the gradient
of the Laplacian of our plane-wave expression.
743
01:11:43.400 --> 01:11:48.740
I'll comment later on that it's
very easy to calculate...
744
01:11:48.840 --> 01:11:50.980
because it's very easy to calculate the kinet ...
745
01:11:51.080 --> 01:11:54.580
the Laplacian is very easy to calculate
the kinetic energy operator.
746
01:11:54.680 --> 01:12:01.020
So remember, if we simplify for a moment and we
think our Hamiltonian as being a kinetic energies
747
01:12:01.120 --> 01:12:07.140
plus a potential and the kinetic energy
is just minus 1/2 Del square...
748
01:12:07.240 --> 01:12:11.120
.
the matrix element
749
01:12:11.280 --> 01:12:16.943
of the kinetic energy between two plane-waves is
super easy to calculate because basically
750
01:12:16.944 --> 01:12:21.600
if I have a plane-wave e to the i G dot r,
751
01:12:21.601 --> 01:12:28.143
its gradient is equal to i G times e to the i G dot r,
752
01:12:28.144 --> 01:12:34.720
so is here there are derivative
automatically and its Laplacian
753
01:12:36.240 --> 01:12:42.943
is going to be equal to minus G square e to the i G dot r.
754
01:12:43.114 --> 01:12:49.540
So twice, you know, I have taken the,
if you want, divergence of the gradient.
755
01:12:49.640 --> 01:12:52.860
So very, very easy to deal with this
756
01:12:52.960 --> 01:12:58.729
and it's also very, very easy to go
from real to reciprocal space.
757
01:12:58.730 --> 01:13:06.400
So I have a quantity that is written, say,
as in this case, the charge density that is
758
01:13:06.500 --> 01:13:16.580
periodic as a linear expansion
in plane-waves with coefficient rho of G.
759
01:13:16.680 --> 01:13:21.440
So either I know the coefficients rho of G,
760
01:13:21.720 --> 01:13:23.660
the coefficient of the plane-waves,
761
01:13:23.760 --> 01:13:31.020
and so I have the reciprocal space expression
of the rho G but using a fast Fourier transform,
762
01:13:31.120 --> 01:13:35.380
it's very, very fast to go from the rho G
763
01:13:35.480 --> 01:13:46.620
to the rho of r, where r belong to a direct space mesh,
764
01:13:46.720 --> 01:13:54.620
say this could be the Wigner-Seitz cell
of my crystal, belong to an equispaced
765
01:13:54.720 --> 01:14:00.220
Wigner-Seitz cell that is nothing else
than the dual of my set of G vectors.
766
01:14:00.320 --> 01:14:05.140
So fast Fourier transform, allow us if we have
767
01:14:05.240 --> 01:14:14.220
the rho sampled in real space on a grid
to give us the coefficients
768
01:14:14.320 --> 01:14:16.460
rho of G very efficiently.
769
01:14:16.560 --> 01:14:24.340
The cost is just a linear times a logarithm
of the linear cost, so it's just above linear.
770
01:14:24.440 --> 01:14:27.560
Or if I have the rho G to go back to the rho
771
01:14:27.660 --> 01:14:32.980
on a mesh, so I can go back from
real and reciprocal space very, very easily.
772
01:14:33.080 --> 01:14:37.180
And, you know, one of the fundamental
concepts that you'll use
773
01:14:37.280 --> 01:14:41.686
when you do a Quantum Espresso or any kind,
you know, it could be a VASP, it could be
774
01:14:41.687 --> 01:14:50.471
a CASTEP, it could be an abinit calculation,
"How do you select what is your plane-wave basis set?"
775
01:14:50.472 --> 01:14:55.140
Well, suppose this is my reciprocal space,
776
01:14:55.240 --> 01:14:57.080
so these are my G vectors.
777
01:14:57.180 --> 01:14:58.260
This is gamma.
778
01:14:58.360 --> 01:15:01.540
This is the set of G vectors.
779
01:15:01.640 --> 01:15:05.000
So to every one of these generic G vectors,
780
01:15:05.100 --> 01:15:10.860
this is G1, this is G2, I can give them names,
there corresponds a plane-wave.
781
01:15:10.960 --> 01:15:15.380
The larger the G is, so this one has
782
01:15:15.480 --> 01:15:21.780
a periodicity that is shorter than this one,
and this one has a periodicity that is longer...
783
01:15:21.880 --> 01:15:24.020
Sorry.
The opposite.
784
01:15:24.120 --> 01:15:25.640
So
785
01:15:25.840 --> 01:15:34.440
if I call this A, B, C and D, A has
a periodicity that is longer than B, that is...
786
01:15:34.720 --> 01:15:37.500
that has two wavelengths in the space
787
01:15:37.501 --> 01:15:43.320
that A has one wavelength and C has
three times, three wavelengths,
788
01:15:43.760 --> 01:15:47.140
and D has four wavelengths and so as I go
789
01:15:47.240 --> 01:15:53.700
to larger and larger G vectors I have plane-waves
that become more and more
790
01:15:53.800 --> 01:16:00.380
detailed and so they express better and better
791
01:16:00.480 --> 01:16:04.740
the fine features of my wavefunction.
792
01:16:04.840 --> 01:16:09.300
So how do I choose what is my basis set?
793
01:16:09.400 --> 01:16:19.280
I choose the radius of a sphere centred around
the origin and if I give you this radius,
794
01:16:19.560 --> 01:16:26.300
what I am saying is you are going to use all
the G vectors that follow inside of this sphere.
795
01:16:26.400 --> 01:16:30.380
And then I can choose a larger radius.
796
01:16:30.480 --> 01:16:33.940
and what I have is a basis set that contains
797
01:16:34.040 --> 01:16:40.420
all the previous basis set elements and more
so it has become systematically better.
798
01:16:40.520 --> 01:16:47.140
So choosing a cutoff radius is a very
economical and efficient way to express
799
01:16:47.240 --> 01:16:55.300
what is the basis set that I'm going to use and,
actually, rather than using the cutoff radius,
800
01:16:55.400 --> 01:17:00.400
we give the energy that corresponds to that radius.
801
01:17:00.500 --> 01:17:06.680
Because if let's say this corresponds to a vector GMAX,
802
01:17:07.200 --> 01:17:12.440
the kinetic energy of that GMAX vector...
803
01:17:17.840 --> 01:17:22.500
remember the Laplacian of e to the i GMAX...
804
01:17:22.600 --> 01:17:30.820
is going to be minus G square, and so this
kinetic energy is equal to 1/2 G square.
805
01:17:30.920 --> 01:17:39.220
So this quantity is called, in Quantum Espresso, ECUT.
806
01:17:39.320 --> 01:17:46.180
It could be the wavefunction cutoff
or it could be the charge density cutoff.
807
01:17:46.280 --> 01:17:48.500
Let me make a little detour.
808
01:17:48.600 --> 01:17:58.660
If I have that, say, a wavefunction
is written as u times e to the i k dot r,
809
01:17:58.760 --> 01:18:02.460
I dropped the index n just for simplicity,
810
01:18:02.560 --> 01:18:08.800
and if I look at the square modulus of psi k,
that is the charge density.
811
01:18:08.900 --> 01:18:15.900
So the charge density rho k for these orbitals
is equal to the square modulus of psi k
812
01:18:16.000 --> 01:18:22.020
and is equal to the square modulus of uk,
because when I take the square modulus,
813
01:18:22.120 --> 01:18:25.180
the i, ik, r modulation goes away.
814
01:18:25.280 --> 01:18:34.620
Well, you see, if my uk was written as
an expansion in plane-waves with,
815
01:18:34.720 --> 01:18:46.186
you know, certain G vectors of the coefficient,
c, I call them kG e to the i G dot r
816
01:18:46.360 --> 01:18:55.060
where G is a modulus that is
less than the G MAX cutoff.
817
01:18:55.160 --> 01:19:02.620
But now if I'm looking not at uk, but I'm
looking at the square modulus of uk.
818
01:19:02.720 --> 01:19:14.495
This is uk star times uk, so what I have here is
sum over G prime ck G prime star
819
01:19:14.695 --> 01:19:18.740
e to the minus i G prime dot r.
820
01:19:18.840 --> 01:19:29.380
This is uk star times uk, that is
sum over G ckG e to the i G dot r.
821
01:19:29.480 --> 01:19:31.871
And lo and behold you discover,
822
01:19:31.872 --> 01:19:38.040
I'll write it here that the charge density,
that is the square modulus of psi,
823
01:19:38.140 --> 01:19:46.043
that is the square modulus of u is written also
as a plane-wave expansion that is given by
824
01:19:46.243 --> 01:19:58.060
ck star G prime ckG e to the i G minus G prime dot r.
825
01:19:58.160 --> 01:20:00.020
So what do we discover?
826
01:20:00.120 --> 01:20:02.760
You see, we discover that if the cutoff
827
01:20:02.860 --> 01:20:07.980
for the wavefunction was this,
that is my G vectors were
828
01:20:08.080 --> 01:20:15.500
inside this cutoff sphere,
the charge density is described by
829
01:20:15.600 --> 01:20:21.640
plane-waves that can have
a wave vector that span all possible combination
830
01:20:21.740 --> 01:20:29.700
of G minus G prime and of course,
G minus G prime, if G and G prime live inside
831
01:20:29.800 --> 01:20:34.840
this blue sphere, so G minus G prime are going
832
01:20:34.940 --> 01:20:40.140
to live in a sphere that has twice the radius,
833
01:20:40.240 --> 01:20:47.860
because if I take G here and I take G prime here,
G minus G prime is going to be here.
834
01:20:47.960 --> 01:20:54.780
So the cutoff radius for the wavefunction
835
01:20:54.880 --> 01:21:02.420
is going to give me naturally a cutoff radius
for the charge density that is twice as large.
836
01:21:02.520 --> 01:21:05.400
And because we are not talking about radii,
837
01:21:05.500 --> 01:21:15.700
but we are talking with energies associated
with radii, it tells me that the energy cutoff
838
01:21:15.800 --> 01:21:21.460
for the wavefunction is related to the energy
839
01:21:21.560 --> 01:21:26.780
cutoff for the charge density by a factor of four.
840
01:21:26.880 --> 01:21:32.220
So if I decided that the energy cutoff for my wave
function is 30 Rydberg,
841
01:21:32.320 --> 01:21:38.660
the energy cutoff for my charge density
needs to be 120 Rydberg.
842
01:21:38.760 --> 01:21:41.260
So that was one statement.
843
01:21:41.360 --> 01:21:46.700
The second statement that is also important
to know - let's do this in black is that,
844
01:21:46.800 --> 01:21:51.440
as I said, Laplacian are easy,
and I'll give you an example of how am I going
845
01:21:51.540 --> 01:21:56.760
to solve the Poisson equation
using a plane-wave expansion.
846
01:21:56.860 --> 01:22:00.000
So Poisson equation is the Laplacian of V is equal
847
01:22:00.100 --> 01:22:06.740
to 4 pi of rho and suppose that I know that rho is
848
01:22:06.840 --> 01:22:11.120
given by a plane-wave expansion
where the rho G coefficient I know.
849
01:22:11.220 --> 01:22:12.960
So,
850
01:22:14.720 --> 01:22:17.040
I know what these are,
851
01:22:17.960 --> 01:22:19.740
my plane-wave expansion.
852
01:22:19.840 --> 01:22:26.340
And V, I don't know, I want to solve
this equation and I write V also
853
01:22:26.440 --> 01:22:35.029
as a plane-wave expansion, but the VG are not known.
854
01:22:36.520 --> 01:22:43.200
OK, so I want to solve this differential equation
and not in real space, but in reciprocal space.
855
01:22:43.300 --> 01:22:47.940
Given all the set of rho G s
what is all the set of VGs?
856
01:22:48.040 --> 01:22:50.940
And that is, you see, super easy because I
857
01:22:51.040 --> 01:23:03.628
insert these equations here and what I get is
sum over G VG times minus G square e to the i G dot r,
858
01:23:03.629 --> 01:23:11.171
is equal to 4 pi sum over G rho G e to the i G dot r.
859
01:23:11.400 --> 01:23:16.360
So plane-waves happen to be also orthogonal,
860
01:23:16.680 --> 01:23:20.386
and orthonormal, I won't go into the
nomalisation of periodic function here,
861
01:23:20.387 --> 01:23:29.814
but so if I take the integral of this
expression times the complex conjugate
862
01:23:30.029 --> 01:23:35.470
of any plane-wave I will have that
because of orthogonality
863
01:23:35.471 --> 01:23:38.940
that contribution is zero, and instead if
864
01:23:39.040 --> 01:23:45.540
my scalar product is with that plane-wave that has
wave vector G, that is not going to be zero.
865
01:23:45.640 --> 01:23:54.080
And so I get right away the solution that is
coefficient by coefficient.
866
01:23:54.520 --> 01:24:02.040
These plane-waves need to have
the same coefficient and so if I have the
867
01:24:03.560 --> 01:24:06.260
charge density coefficients,
868
01:24:06.360 --> 01:24:10.600
the coefficients of the potential... it could be, you see,
869
01:24:10.700 --> 01:24:17.740
the Hartree potential that solves the Poisson
equation... are just given plane-wave by plane-wave
870
01:24:17.840 --> 01:24:24.480
by the coefficient of the charge
density divided by G square.
871
01:24:24.640 --> 01:24:31.940
Now, what happens when you have a very, very
large system, a very, very large supercell.
872
01:24:32.040 --> 01:24:37.260
Remember that reciprocal space goes
as the inverse of direct space.
873
01:24:37.360 --> 01:24:41.900
So as I make my unit cell larger and larger,
874
01:24:42.000 --> 01:24:49.760
my G vector of my reciprocal lattice becomes, you
know, made by vectors that are closer and closer.
875
01:24:49.860 --> 01:24:52.400
If you remember, the expression of the dual said
876
01:24:52.500 --> 01:24:59.020
the volume was at the denominator,
so if the direct lattice vectors goes up linearly,
877
01:24:59.120 --> 01:25:04.740
the reciprocal lattice vector,
decreases inversely proportional to that.
878
01:25:04.840 --> 01:25:12.040
So in a big system, I will have a lot of G vectors
that are very small and if the G can be very
879
01:25:12.140 --> 01:25:17.580
small, it means that a little error
of the charge density is going to give me
880
01:25:17.680 --> 01:25:24.420
a large error, because it's divided by G square,
on the equivalent component for the potential.
881
01:25:24.520 --> 01:25:29.520
And that's why electrostatics
of large system is very tricky.
882
01:25:29.620 --> 01:25:32.660
Things that you're studying by first principle:
883
01:25:32.760 --> 01:25:40.140
a metal slab with vacuum or maybe with some
molecule looking here and smiling at you.
884
01:25:40.240 --> 01:25:43.040
Well, and you're solving this in reciprocal space
885
01:25:43.140 --> 01:25:46.640
so you're taking this charge density,
you're Fourier transforming,
886
01:25:46.740 --> 01:25:51.480
you're getting the rho G, and then out of that,
you're getting the VG are small.
887
01:25:51.580 --> 01:25:56.560
Error on rho G gives you a big error on the VG
and so it means that all of a sudden you could put
888
01:25:56.660 --> 01:26:01.520
the charge density here,
and this is a large system so in the next
889
01:26:01.620 --> 01:26:04.940
iteration, there will be an electrostatic
imbalance,
890
01:26:05.040 --> 01:26:11.880
that was to send back the charge density here
and this divergence in the infinite supercell
891
01:26:11.980 --> 01:26:18.960
limit gives rise to what are called sloshing
instabilities that make the convergence,
892
01:26:19.060 --> 01:26:26.160
the iterative convergence,
of first principle calculation slightly delicate.
893
01:26:26.960 --> 01:26:33.860
OK, so, as I said, we use plane-waves
but we didn't have to do that.
894
01:26:33.960 --> 01:26:36.940
A lot of quantum chemistry codes
895
01:26:37.040 --> 01:26:41.060
use Gaussian basis set because
the exchange integrals are easy to do.
896
01:26:41.160 --> 01:26:44.540
There are codes like CRYSTAL that do
897
01:26:44.640 --> 01:26:47.740
solids with Gaussian basis set.
898
01:26:47.840 --> 01:26:53.280
There are simplified codes that use
linear combination of atomic orbitals.
899
01:26:53.520 --> 01:26:56.360
There are codes that use a linear muffin-tin
900
01:26:56.360 --> 01:26:59.700
approach that has been introduced
by O. K. Andersen in the seventies.
901
01:26:59.800 --> 01:27:01.680
I don't know if you know what a muffin tin is.
902
01:27:01.960 --> 01:27:03.020
Muffin tins.
903
01:27:03.120 --> 01:27:08.040
Muffins have become very popular ever
since cupcakes have become very popular.
904
01:27:08.140 --> 01:27:11.860
So a muffin tin is something
to bake cupcakes in the oven.
905
01:27:11.960 --> 01:27:17.140
And so a muffin-tin approach
is an approach where seen from the above
906
01:27:17.240 --> 01:27:21.980
we partition space in atomic
spheres around the atoms
907
01:27:22.080 --> 01:27:30.100
and we use one set of basis sets inside
and in the interstitial valence region
908
01:27:30.200 --> 01:27:36.280
we use plane-waves and we have to match
the basis set at the boundaries very well.
909
01:27:36.960 --> 01:27:40.880
And that, in particular, LMTO
was very powerful when it was first
910
01:27:40.980 --> 01:27:45.500
introduced because it could allow to do
calculations that before were impossible.
911
01:27:45.600 --> 01:27:47.200
These days we even do better.
912
01:27:47.300 --> 01:27:49.700
We use linearised method
913
01:27:49.800 --> 01:27:56.900
plane-wave techniques that again use
atomic like basis sets inside the sphere
914
01:27:57.000 --> 01:28:01.640
and plane-waves in the middle
and allow ultimate accuracy.
915
01:28:01.740 --> 01:28:05.820
So, often there are a lot of reference calculation
916
01:28:05.920 --> 01:28:11.740
or calculation on very delicate properties
that can be done rather than using plane-waves and
917
01:28:11.840 --> 01:28:19.140
pseudopotential, using all-electron technique
that would be impossible to deal with,
918
01:28:19.240 --> 01:28:26.160
with plane-waves, but with LAPW
you can use atomic orbitals around
919
01:28:26.260 --> 01:28:31.700
the nucleus and plane-waves in the middle
so you can study the entire system.
920
01:28:31.800 --> 01:28:37.020
But again, Quantum Espresso like VASP, like
abinit, like CASTEP, like many other codes,
921
01:28:37.120 --> 01:28:41.740
we use plane-waves and that's how we are going to use.
922
01:28:41.840 --> 01:28:44.800
The other technical concept that you need to be
923
01:28:44.900 --> 01:28:50.940
familiar to do an ab initio calculation
is the one of Brilouin zone integration.
924
01:28:51.040 --> 01:28:53.820
As I said before, there are a lot of quantities
925
01:28:53.920 --> 01:28:58.980
that require summing over all the states in the system.
926
01:28:59.080 --> 01:29:01.840
Say, let's take silicon, an insulator.
927
01:29:01.940 --> 01:29:07.340
My charge density would be given by the sum over
928
01:29:07.440 --> 01:29:12.220
all the bands and by integral over all the k
929
01:29:12.320 --> 01:29:21.960
in the first Brillouin zone
of the square modulus of the unk.
930
01:29:22.080 --> 01:29:25.480
So we can't do continuous integrals in our
931
01:29:25.580 --> 01:29:30.120
calculations, so what we are going to do is we are
going to do sampling that means we take
932
01:29:30.220 --> 01:29:36.780
the Brillouin zone and we sample it
on a finite discrete mesh of k-points.
933
01:29:36.880 --> 01:29:41.320
So we don't consider all the blue states,
934
01:29:41.420 --> 01:29:47.700
but we consider all the black states
and by increasing the fineness
935
01:29:47.800 --> 01:29:53.340
of our k-point sampling mesh,
we do better and better.
936
01:29:53.440 --> 01:30:01.540
And the general technique is to use
equispaced meshes that are called
937
01:30:01.640 --> 01:30:07.140
Monkhorst-Pack meshes - I'll show them
in a moment - so we are doing something, actually,
938
01:30:07.240 --> 01:30:13.440
quite trivial. We are sampling there
the Brillouin zone on an equispaced mesh.
939
01:30:13.840 --> 01:30:15.380
There are a couple of subtleties.
940
01:30:15.480 --> 01:30:16.920
One that I'll describe later.
941
01:30:17.020 --> 01:30:25.880
If you have a metallic system you have a
Fermi energy that separates the empty states
942
01:30:26.720 --> 01:30:29.600
from the filled states.
943
01:30:30.200 --> 01:30:37.780
And if you want to calculate the charge density,
say, in copper, you have only to sum over
944
01:30:37.880 --> 01:30:44.920
the occupied states and so what you have
in a metal is actually a discontinuity.
945
01:30:44.920 --> 01:30:45.660
You might be, you know,
946
01:30:45.760 --> 01:30:52.420
following a band and integrate things but as soon
as that band crosses the Fermi energy
947
01:30:52.520 --> 01:30:59.620
it doesn't contribute any more to the...
to the charge density of the system.
948
01:30:59.720 --> 01:31:02.028
You know, the simplest case that could be a case of,
949
01:31:02.029 --> 01:31:07.200
you know, a very simple, model one-dimensional
system that there's a Fermi energy here.
950
01:31:07.300 --> 01:31:08.700
This is the band structure.
951
01:31:08.800 --> 01:31:12.100
So what I have when I do the integral,
952
01:31:12.200 --> 01:31:17.420
I need to count the states but as soon as
the states cross the Fermi energy,
953
01:31:17.520 --> 01:31:23.940
their occupation goes to zero, so you have
to integrate a discontinuous function.
954
01:31:24.040 --> 01:31:29.600
Discontinuous function are much nastier to
integrate carefully rather than smooth functions.
955
01:31:29.700 --> 01:31:38.300
So insulators are easy to sample in k-space than
metals and we'll get to metals in a moment.
956
01:31:38.400 --> 01:31:41.480
The other limit that is interesting is the limit
957
01:31:41.580 --> 01:31:49.000
when you study a very large system and maybe
the simplest way to think of this is, again,
958
01:31:49.100 --> 01:31:57.220
taking this band structure here in for a one-
dimensional system that was looking like this.
959
01:31:57.320 --> 01:31:58.660
You know, this could be, you know,
960
01:31:58.760 --> 01:32:06.120
the band structure of a simple monoatomic chain,
and as you know, I'm in real space.
961
01:32:06.220 --> 01:32:09.680
I could study this in the primitive unit cell,
962
01:32:09.780 --> 01:32:13.956
maybe one atom per unit cell,
but no one stops me from studying it
963
01:32:13.957 --> 01:32:18.700
in a double unit cell and a triple unit cell
and a quadruple unit cell.
964
01:32:18.800 --> 01:32:22.460
And in reciprocal space,
let's say if I double the unit cell,
965
01:32:22.560 --> 01:32:31.240
that would correspond to actually having the same
system but described in a different way.
966
01:32:31.340 --> 01:32:37.740
Rather than having this Brillouin zone,
now you have a Brillouin zone that is one half
967
01:32:37.840 --> 01:32:45.520
and rather than having one band, in black,
I have, let's say two bands, in orange.
968
01:32:45.620 --> 01:32:54.120
So the black band and the two orange bands
are describing exactly the same system.
969
01:32:57.520 --> 01:33:00.500
But in a different representation,
970
01:33:00.600 --> 01:33:06.520
the black one comes with having in direct space,
one atom per unit cell the orange comes
971
01:33:06.620 --> 01:33:14.820
from having two atoms per unit cell. And as I
increase my description in real space,
972
01:33:14.920 --> 01:33:21.580
what I get is that my bands
get folded and folded and folded.
973
01:33:21.680 --> 01:33:27.340
So, you know, I might end up
with having something that looks like this.
974
01:33:27.440 --> 01:33:30.840
So you see, if I study this system in a very,
975
01:33:30.940 --> 01:33:37.720
very big supercell, I have a very small Brillouin
zone and I don't have to do a lot of effort
976
01:33:37.820 --> 01:33:43.320
in sampling it because I just sample it at one
point, because when I sample it at one point,
977
01:33:43.420 --> 01:33:51.040
I get all these point here in the middle
that corresponds to taking all these point here.
978
01:33:51.200 --> 01:33:54.460
So if your unit cell is very large, you need...
979
01:33:54.560 --> 01:33:57.940
you can get away with very coarse sampling.
980
01:33:58.040 --> 01:33:59.860
Maybe just one point.
981
01:33:59.960 --> 01:34:01.540
It could be the gamma point.
982
01:34:01.640 --> 01:34:03.500
It could be the zone boundary.
983
01:34:03.600 --> 01:34:07.340
Or you could ask yourself the question, very smart
984
01:34:07.440 --> 01:34:14.500
of, you know, if you can only afford one point,
what is the best point that you can choose?
985
01:34:14.600 --> 01:34:16.240
And it might not be the origin.
986
01:34:16.340 --> 01:34:18.100
It might be something else.
987
01:34:18.200 --> 01:34:20.580
And that's called the Baldereschi point,
988
01:34:20.680 --> 01:34:26.760
and if you're actually dealing with an FCC
lattice, that's 1/4 1/4 2 pi over a.
989
01:34:29.760 --> 01:34:32.100
This was especially relevant when, you know,
990
01:34:32.200 --> 01:34:37.340
calculation was so expensive
that you could, you know, use...
991
01:34:37.440 --> 01:34:40.120
you could spare any CPU that you...
992
01:34:40.220 --> 01:34:43.080
that you could, so you wanted to sample those as
993
01:34:43.180 --> 01:34:46.200
little as possible, but it's
still very relevant these days.
994
01:34:46.300 --> 01:34:48.960
If you are studying your system in a big super
995
01:34:49.060 --> 01:34:55.320
cell, then maybe you can get away with only one
point, but then maybe the Baldereschi point is
996
01:34:55.420 --> 01:34:59.780
better than the gamma point because
it's actually more accurate.
997
01:34:59.880 --> 01:35:02.140
The gamma point has another advantage,
998
01:35:02.240 --> 01:35:10.660
that is at gamma point the wavefunction are real,
so those coefficients don't require
999
01:35:10.760 --> 01:35:18.500
a complex expansion, but you can just
get away with half of the informational
1000
01:35:18.600 --> 01:35:26.620
content and Quantum Espresso in the gamma only
formulation takes advantage, actually, of this.
1001
01:35:26.720 --> 01:35:28.920
But in general, this will be what we will do.
1002
01:35:29.020 --> 01:35:31.660
We will use equispaced meshes.
1003
01:35:31.760 --> 01:35:38.020
Those are called Monkhorst and Pack from the
paper in the mid 70s that describe those.
1004
01:35:38.120 --> 01:35:43.260
Quantum Espresso generates automatically for you,
so you don't have to worry very much.
1005
01:35:43.360 --> 01:35:45.460
But just so that you know what they are,
1006
01:35:45.560 --> 01:35:52.180
what I have on the right, the (4,4,4)
unshifted mesh, you see, is just
1007
01:35:52.280 --> 01:35:57.500
sampling the Brillouin zone with
1008
01:35:57.600 --> 01:36:00.840
four points in every direction and, of course,
1009
01:36:00.940 --> 01:36:05.940
this point here with the cross is equivalent
by symmetry with this point here.
1010
01:36:06.040 --> 01:36:09.760
And you can have the unshifted mesh or you could
1011
01:36:09.860 --> 01:36:15.660
have with the same fineness of sampling
what is called the shifted mesh
1012
01:36:15.760 --> 01:36:22.780
where this k-point sampling are
here in the middle of this division.
1013
01:36:22.880 --> 01:36:26.920
And Quantum Espresso would take these
1014
01:36:27.120 --> 01:36:28.480
as a 4 4 4 1 1 1, 4 4 4 0 0 0.
1015
01:36:30.880 --> 01:36:32.720
One thing that is very important.
1016
01:36:32.820 --> 01:36:37.260
Suppose for a moment,
we are studying, and this is a real space cell,
1017
01:36:37.360 --> 01:36:44.260
a carbon nanotube that, you know,
we might describe by something like this
1018
01:36:44.360 --> 01:36:48.686
In reciprocal space, you'll have also a Brillouin zone
1019
01:36:48.857 --> 01:36:56.900
but what happens in your system is that this
y direction is the direction of periodicity
1020
01:36:57.000 --> 01:37:03.140
and so is the direction in reciprocal
space where there will be band dispersion.
1021
01:37:03.240 --> 01:37:10.640
In reciprocal space, in the x and z direction,
there will be no dispersion because the system,
1022
01:37:10.740 --> 01:37:17.500
the nanotube, is effectively isolated and there's
not talking with these periodic images.
1023
01:37:17.600 --> 01:37:25.020
So the states don't have a k dependence
in the expectation value for the energy.
1024
01:37:25.120 --> 01:37:30.920
So the Monkhorst-Pack mesh that you would use
in this case is not maybe a sequence of better
1025
01:37:31.020 --> 01:37:35.480
and better Monkhort-Pack meshes,
be 2 2 2, 3 3 3, 4 4 4
1026
01:37:36.160 --> 01:37:38.340
until your quantities are converged.
1027
01:37:38.440 --> 01:37:43.360
You don't need to sample in two directions
so they will be 1 2 1, 1 3 1, 1 4 1
1028
01:37:44.920 --> 01:37:50.120
and, you know, 1 4 1 is going
to cost 16 times less than 4 4 4.
1029
01:37:50.220 --> 01:37:54.120
You have 64 point in one case
and only 4 points in the other.
1030
01:37:54.220 --> 01:37:58.140
So that is also very important to remember.
1031
01:37:58.240 --> 01:38:05.180
If your system is not truly three-dimensional,
but is separated by vacuum.
1032
01:38:05.280 --> 01:38:07.800
In this case, two direction, as is the case
1033
01:38:07.900 --> 01:38:12.180
for a carbon nanotube or in three
direction, if this was a molecule.
1034
01:38:12.280 --> 01:38:17.700
For a molecule, you always only
use gamma sampling, you know.
1035
01:38:17.800 --> 01:38:23.940
For a molecule, there is no point in sampling
1036
01:38:24.040 --> 01:38:30.160
anything because the molecular states
don't have any k-dispersion.
1037
01:38:30.160 --> 01:38:31.480
If they have a k-dispersion,
1038
01:38:31.580 --> 01:38:36.420
it means that you don't have a molecule,
but the molecule is feeling its periodic image,
1039
01:38:36.520 --> 01:38:40.900
so you are actually having something
that is an organic crystal.
1040
01:38:41.000 --> 01:38:43.380
And if you have, say, a layer of graphene,
1041
01:38:43.480 --> 01:38:49.720
you sample in the plane,
but you don't sample perpendicularly.
1042
01:38:50.280 --> 01:38:56.460
The last thing that makes, actually,
our calculations less expensive in solids
1043
01:38:56.560 --> 01:39:01.600
is exploiting symmetry operation,
the point group symmetry of the system.
1044
01:39:01.671 --> 01:39:07.700
You might have mirror planes,
you might have a rotation axis.
1045
01:39:07.800 --> 01:39:09.900
If you are studying, you know,
1046
01:39:10.000 --> 01:39:17.660
a silicon crystal you actually have 48 symmetry
operation on the Td point group.
1047
01:39:17.760 --> 01:39:24.260
And symmetry is used very
effectively in doing these sums,
1048
01:39:24.360 --> 01:39:29.260
because rather than, say, in this case,
if this was my Brillouin zone,
1049
01:39:29.360 --> 01:39:35.820
rather than summing over all the k
that belong to the Brillouin zone,
1050
01:39:35.920 --> 01:39:43.220
I can use, I will not demonstrate it,
but I can use this relation that tells me
1051
01:39:43.320 --> 01:39:48.300
that if I know the orbital at a certain k,
1052
01:39:48.400 --> 01:39:54.280
I can obtain the Bloch orbital at a certain...
1053
01:39:54.380 --> 01:39:55.920
sorry, this is unreadable...
1054
01:39:56.020 --> 01:39:59.600
at a certain
1055
01:40:01.520 --> 01:40:04.680
S of k that is at a k that is related
1056
01:40:04.780 --> 01:40:12.500
to the previous one by a symmetry operation,
simply counter rotating my points in real space.
1057
01:40:12.600 --> 01:40:18.860
That is to say that if this point here is k
1058
01:40:18.960 --> 01:40:28.820
and this point here is S of k, and they are
related by a minor operation, if I have psi k,
1059
01:40:28.920 --> 01:40:32.100
the psi Sk
1060
01:40:32.200 --> 01:40:37.300
is obtained by exactly folding back
1061
01:40:37.400 --> 01:40:40.760
the points according to this operation.
1062
01:40:40.860 --> 01:40:44.540
So my sum is broken apart.
1063
01:40:44.640 --> 01:40:50.180
The sum over all the k-points it becomes a sum
over only the k points
1064
01:40:50.280 --> 01:40:57.680
that live in what is called the irreducible wedge
of this segment of the Brillouin zone that once we
1065
01:40:57.780 --> 01:41:02.800
apply all possible symmetry operation
blankets the entire Brillouin zone.
1066
01:41:02.900 --> 01:41:10.520
Here I have only one mirror symmetry and so you
see the orange wedge plus the middle symmetry is
1067
01:41:10.620 --> 01:41:13.980
going to cover also the other
half of the Brillouin zone.
1068
01:41:14.080 --> 01:41:21.440
So here I do the sum over the irreducible wedge
and all other points that are outside
1069
01:41:21.540 --> 01:41:27.600
the irreducible wedge are not obtained
by doing additional calculation, by simply...
1070
01:41:27.700 --> 01:41:31.420
but by simply, using this symmetry relation here.
1071
01:41:31.520 --> 01:41:33.500
So you see, at the end of the day,
1072
01:41:33.600 --> 01:41:40.940
I have to solve the Kohn-Sham equation only
at the irreducible points of the Brillouin Zone.
1073
01:41:41.040 --> 01:41:46.900
And, as I said, if I'm studying
silicon is a saving of 48.
1074
01:41:46.943 --> 01:41:54.600
Of course, you know, a factor of 48 is not,
you know, life changing these days so
1075
01:41:54.720 --> 01:42:01.220
sometimes we don't bother with this,
but in very complex codes that use a lot of
1076
01:42:01.320 --> 01:42:04.360
calculation for, say,
calculating phonons and electrons
1077
01:42:04.460 --> 01:42:08.920
and electron-phonon interaction,
you know, these savings build up very quickly.
1078
01:42:09.960 --> 01:42:13.960
While symmetry, that is always there
if we don't have a magnetic field,
1079
01:42:14.060 --> 01:42:20.440
so if we don't have a vector potential in our
Hamiltonian, is time reversal symmetry.
1080
01:42:20.540 --> 01:42:23.880
So what I written here is the Hamiltonian,
1081
01:42:24.120 --> 01:42:27.500
not any more for the Bloch function psi.
1082
01:42:27.600 --> 01:42:32.280
but what happens if I apply H to psi,
1083
01:42:32.380 --> 01:42:38.860
that is if I apply H to unk, e to the i k r
1084
01:42:38.960 --> 01:42:47.380
so if I put explicitly the Bloch form it turns out
that the Hamiltonian for unk has this form here,
1085
01:42:47.480 --> 01:42:57.420
and you see right away that if I solve this
Hamiltonian for unk the solution for un minus k
1086
01:42:57.520 --> 01:43:03.320
is immediately given by taking
the complex conjugate of unk.
1087
01:43:03.420 --> 01:43:06.680
So always if this is my Brillouin zone...
1088
01:43:06.780 --> 01:43:08.660
OK, this is the gamma point...
1089
01:43:08.760 --> 01:43:15.540
I never need more than half of the Brillouin zone
because the points in the other half...
1090
01:43:15.640 --> 01:43:16.940
OK, if, you know,
1091
01:43:17.040 --> 01:43:22.660
the points on the left here have, say, negative k,
the points on the right have positive k,
1092
01:43:22.760 --> 01:43:30.780
I get from the one to the other, that is I go from
k to minus k, just taking the complex conjugate.
1093
01:43:30.880 --> 01:43:32.720
So you never have to do a full sampling
1094
01:43:32.820 --> 01:43:40.120
of the Brillouin zone but worst case scenario
you just have to sample half of it.
1095
01:43:41.400 --> 01:43:43.980
PWSCF, that is the self-consistent code
1096
01:43:44.080 --> 01:43:48.480
in Quantum Espresso, is very
sophisticated about k-points.
1097
01:43:48.580 --> 01:43:50.660
It gives you all these option.
1098
01:43:50.760 --> 01:43:54.260
It gives you the option of just choosing
1099
01:43:54.360 --> 01:44:01.340
the gamma point for a big super cell and then
your calculation becomes
1100
01:44:01.440 --> 01:44:04.700
twice as fast and twice as
1101
01:44:04.800 --> 01:44:08.060
less memory intensive.
1102
01:44:08.160 --> 01:44:12.980
Most of the time you lose
automatic k-point meshes,
1103
01:44:13.080 --> 01:44:18.140
maybe unshifted with k1, k2, k3
1104
01:44:18.240 --> 01:44:20.080
0 0 0 shifted 1 1 1
1105
01:44:20.360 --> 01:44:28.640
but you need to test when your calculations
are converged with respect to nk1, nk2, nk3.
1106
01:44:29.760 --> 01:44:33.500
And the calculation goes linear
in the product of these three numbers.
1107
01:44:33.600 --> 01:44:39.260
So doing a 2 2 2 mesh is 8 times
1108
01:44:39.360 --> 01:44:45.140
cheaper than doing a 4 4 mesh, that is
in itself 8 times cheaper than doing
1109
01:44:45.240 --> 01:44:51.000
an 8 8 8 mesh, that is in itself 64
times more expensive than a 2 2 2 mesh.
1110
01:44:51.120 --> 01:44:55.020
So you want to, you know, proceed
carefully but you need to check
1111
01:44:55.120 --> 01:45:02.740
that your physical quantities are converged
with respect to the k-point sampling.
1112
01:45:02.840 --> 01:45:06.900
Remember again, you know, if you
have a molecule you just use gamma.
1113
01:45:07.000 --> 01:45:12.460
If you have a nanotube, you sample only one
direction or you have a wire, a nanowire.
1114
01:45:12.560 --> 01:45:19.400
If you have a two-dimensional system or a metal
slab with a molecule absorbed, you don't need...
1115
01:45:19.500 --> 01:45:25.080
you don't want to sample in the
direction perpendicular to this.
1116
01:45:26.480 --> 01:45:30.520
OK, so this was k-point sampling,
so we are starting to look at all
1117
01:45:30.620 --> 01:45:34.560
the computational tools that we need
to run a proper calculation.
1118
01:45:34.660 --> 01:45:40.260
We have looked at the plane-wave cutoff,
the basis set, the k-point sampling.
1119
01:45:40.360 --> 01:45:41.740
Let me introduce
1120
01:45:41.840 --> 01:45:47.280
in very handwaving form the concept
of pseudopotentials, that is the idea that,
1121
01:45:47.380 --> 01:45:51.160
actually, when you look at a solid or you look
1122
01:45:51.260 --> 01:45:56.540
at a molecule, the electrons deep down in the core,
1123
01:45:56.640 --> 01:46:00.860
say, for aluminum, it would be
the 1s, the 2s and the 2p electrons,
1124
01:46:00.960 --> 01:46:09.329
they are so tightly bound to the nucleus that they
don't care if aluminum is in a molecule,
1125
01:46:09.330 --> 01:46:14.980
is in a bulk solid, is at the surface
or is in molten aluminum.
1126
01:46:15.080 --> 01:46:18.800
And so for this reason, what we want to do is only
1127
01:46:18.900 --> 01:46:25.120
study the valence electrons and
the core electrons that are very low
1128
01:46:25.220 --> 01:46:33.040
in energy, we consider frozen
in the atomic configuration.
1129
01:46:33.520 --> 01:46:38.480
Of course, if you want to study,
say, X-ray spectroscopy,
1130
01:46:38.580 --> 01:46:44.160
is where you excite electrons from the core,
well then you need to do something else.
1131
01:46:44.260 --> 01:46:46.500
We will not go into that but that are also
1132
01:46:46.600 --> 01:46:51.740
something that you can calculate
with standard density functional codes.
1133
01:46:51.840 --> 01:46:57.120
So how do we get rid of the need to study core electrons?
1134
01:46:57.220 --> 01:47:01.100
You see, core electrons have very large negative energies,
1135
01:47:01.200 --> 01:47:06.820
so a small error on a core electron messes up our
calculation because, you know,
1136
01:47:06.920 --> 01:47:10.780
binding energies are of the order
of a fraction of an electron volt.
1137
01:47:10.880 --> 01:47:18.100
So doing, you know, a .1% error in the 1s state
completely messes up our calculation.
1138
01:47:18.200 --> 01:47:23.200
So we don't want to deal with this
and the core electrons are very localised around
1139
01:47:23.300 --> 01:47:26.080
the nucleus, so they would need a lot of plane-waves.
1140
01:47:26.180 --> 01:47:31.460
So to do that, we introduce
the concept of the pseudopotential.
1141
01:47:31.560 --> 01:47:33.160
What is a pseudopotential?
1142
01:47:33.160 --> 01:47:34.786
Well, say for the case of aluminum,
1143
01:47:34.787 --> 01:47:41.540
a pseudopotential is a potential
that we create that reproduce the effects,
1144
01:47:41.640 --> 01:47:46.980
not only of the nucleus, minus 13 over r,
the coulombic attraction,
1145
01:47:47.080 --> 01:47:54.780
but it reproduces the effect of the nucleus
screened by the frozen core electrons.
1146
01:47:54.880 --> 01:47:58.700
And so we want to create that pseudopotential
1147
01:47:58.800 --> 01:48:06.860
so that when we throw the three valence electron
of aluminum in the field of those pseudopotential,
1148
01:48:06.960 --> 01:48:15.580
the energy of those three electrons
will be these energies here,
1149
01:48:15.680 --> 01:48:22.780
-7.8 and -2.7, that are exactly
the energy that the true...
1150
01:48:22.880 --> 01:48:25.760
that, sorry, that they would have had if they had
1151
01:48:25.860 --> 01:48:32.300
felt the full explicit electronic core configuration.
1152
01:48:32.400 --> 01:48:35.700
So we build pseudopotentials so that
1153
01:48:35.800 --> 01:48:42.940
valence electron thrown in that pseudopotential
field have the same energies that
1154
01:48:43.040 --> 01:48:50.580
valence electrons thrown in the through field of
the nucleus and all the core states would have.
1155
01:48:50.680 --> 01:48:55.700
But not only that, we want also that
the valence wavefunction is correct
1156
01:48:55.857 --> 01:49:03.880
and we are going to care for the valence
wavefunctions to be correct outside of the core.
1157
01:49:03.980 --> 01:49:08.460
So we want the valence wavefunction to be correct
1158
01:49:08.560 --> 01:49:14.540
outside of the core because that's where
chemical bonding takes place.
1159
01:49:14.640 --> 01:49:21.260
And so here you see what I have as a thin line is
1160
01:49:21.360 --> 01:49:29.140
a coulombic attractive potential if, say,
I have an atom sitting here at the origin.
1161
01:49:29.240 --> 01:49:38.100
In blue, sorry, in solid red, in thick red,
I have a pseudopotential that agrees exactly
1162
01:49:38.200 --> 01:49:45.100
with the true coulombic potential outside
the core and is different inside the core,
1163
01:49:45.200 --> 01:49:54.320
and this thick red pseudopotential has been
built such that the electron in the field
1164
01:49:54.321 --> 01:50:00.336
of the pseudopotential obtains
a ground state that is the thick blue curve,
1165
01:50:00.536 --> 01:50:04.280
and that thick blue curve has two characteristics.
1166
01:50:04.380 --> 01:50:10.243
One is identical to the thin blue curve
that was, you know, this could be the
1167
01:50:10.244 --> 01:50:13.780
6s, 7s state in an atom.
1168
01:50:13.880 --> 01:50:20.843
So outside the core, the pseudo wavefunction
and the true wavefunction are identical
1169
01:50:20.844 --> 01:50:30.086
and also the energy of the blue state, thick,
is equal to the energy of the thin blue state.
1170
01:50:30.480 --> 01:50:31.620
How to do that?
1171
01:50:31.720 --> 01:50:35.570
How to build this pseudopotential is very complex,
1172
01:50:35.571 --> 01:50:41.486
but you gain a lot because all of a sudden
you don't have to study all those core electrons
1173
01:50:41.600 --> 01:50:48.260
that are very high energy, so very little
error are going to be very harmful
1174
01:50:48.360 --> 01:50:53.820
and they are very sharply peaked,
so they require an impossibly large cutoff.
1175
01:50:53.920 --> 01:50:56.900
How do I build these pseudopotentials?
1176
01:50:57.000 --> 01:51:03.700
There were a number of techniques developed,
but really at the end of the 70s Don Hamann,
1177
01:51:03.800 --> 01:51:09.960
together with Michael Schluter and Chiang,
developed these ideas of norm-conserving
1178
01:51:10.060 --> 01:51:16.333
pseudopotentials and they developed a
construction protocol for pseudopotential,
1179
01:51:16.533 --> 01:51:23.595
where as I said, the eigenvalues,
and they agree for the real states
1180
01:51:23.795 --> 01:51:29.100
and for the pseudo states,
the wavefunctions agree
1181
01:51:29.200 --> 01:51:33.300
beyond core radius.
1182
01:51:33.400 --> 01:51:42.300
The concept of norm-conservation told us
that the integral of the pseudo wavefunction
1183
01:51:42.400 --> 01:51:48.520
and the integral of the pseudo and the real
wavefunction both squared, that is the real
1184
01:51:48.620 --> 01:51:55.300
and the pseudo charge are going to be identical
within the core, the integral,
1185
01:51:55.400 --> 01:51:58.014
and then because the pseudo and the
1186
01:51:58.015 --> 01:52:03.100
real wavefunction agree outside the core,
they're also going to agree outside the core.
1187
01:52:03.200 --> 01:52:05.786
And then there is another subtlety
1188
01:52:05.787 --> 01:52:13.240
that I'll just mention in passing,
there is going to be an agreement for,
1189
01:52:13.340 --> 01:52:18.420
and I'll explain what this is in a moment,
what the logarithmic derivatives are,
1190
01:52:18.520 --> 01:52:21.280
and there are first energy derivatives.
1191
01:52:21.380 --> 01:52:23.671
I'll get to that in a moment.
1192
01:52:24.000 --> 01:52:28.229
But in order to make this
norm-conserving pseudopotential,
1193
01:52:28.230 --> 01:52:33.380
it turns out that what you have to do is,
1194
01:52:33.480 --> 01:52:39.880
you know, you can't really do it with just a local
pseudopotential, that is a pseudopotential
1195
01:52:39.980 --> 01:52:41.860
that depends only on position.
1196
01:52:41.960 --> 01:52:49.744
You need to act differently on the wavefunction,
that is, you need to act differently
1197
01:52:49.944 --> 01:52:55.074
on the wavefunction that looks like an s orbital,
or on the wavefunction that looks like a p orbital
1198
01:52:55.274 --> 01:52:57.560
or a wavefunction that looks like a d orbital.
1199
01:52:57.660 --> 01:53:01.700
That is to say, in order to act, your pseudopotential
1200
01:53:01.701 --> 01:53:10.643
takes a scalar product on the spherical
harmonics and you are going to act differently...
1201
01:53:11.240 --> 01:53:15.140
so this would be the matrix element...
1202
01:53:15.240 --> 01:53:19.640
you are going to act differently on the component
1203
01:53:19.740 --> 01:53:27.020
of the wavefunction that has a non-zero projection
onto an s spherical harmonic, onto a p spherical harmonic
1204
01:53:27.120 --> 01:53:30.629
onto a d spherical harmonic
and so on and so forth.
1205
01:53:30.630 --> 01:53:36.560
So, proper norm-conserving pseudopotential will act
1206
01:53:36.660 --> 01:53:39.060
and it's, you know,
1207
01:53:39.160 --> 01:53:44.840
represented here differently on the s component
of the wavefunction, on the p component,
1208
01:53:44.841 --> 01:53:51.120
on the d component, on the f component,
and at a certain point for higher order
1209
01:53:51.220 --> 01:53:56.980
angular momentum components, it will stop
doing that because it's not necessary anymore.
1210
01:53:57.080 --> 01:53:59.380
So maybe for the case of gallium,
1211
01:53:59.480 --> 01:54:02.780
we have an s component,
a p component, a d component,
1212
01:54:02.880 --> 01:54:10.140
and we say for anything that has a higher order,
angular symmetry will just act in the same way
1213
01:54:10.240 --> 01:54:15.660
and that will be the local component
of the pseudopotential that might look like this.
1214
01:54:15.760 --> 01:54:21.460
And again, you see outside a certain core
radius, they all look identical.
1215
01:54:21.560 --> 01:54:24.580
And if I study in this case, the gallium atom,
1216
01:54:24.680 --> 01:54:32.460
you see, I'm going to reproduce the s,
the p and the d wavefunction of the atom
1217
01:54:32.560 --> 01:54:35.540
exactly outside the core.
1218
01:54:35.640 --> 01:54:37.520
So this is inside the core.
1219
01:54:37.620 --> 01:54:39.786
This is outside of the core.
1220
01:54:39.986 --> 01:54:42.880
And inside the core, they are different.
1221
01:54:42.980 --> 01:54:44.160
They are much smoother,
1222
01:54:44.160 --> 01:54:45.600
so you see they are much easier.
1223
01:54:45.700 --> 01:54:50.220
You see how sharp the s wavefunction was here.
1224
01:54:50.320 --> 01:54:51.400
So the 3s,
1225
01:54:51.500 --> 01:54:57.186
so they are much smoother
but they have still the same energy.
1226
01:54:57.520 --> 01:55:02.240
And the last concept was this
concept of a logarithmic derivatives.
1227
01:55:02.340 --> 01:55:04.460
I won't go into the detail of this
1228
01:55:04.560 --> 01:55:10.300
but it comes from scattering theory. The way at
the end of the pseudopotential are built
1229
01:55:10.400 --> 01:55:18.020
comes from the idea that a real core,
that is the nucleus plus the real core electrons
1230
01:55:18.120 --> 01:55:25.760
and a pseudopotential, that is a potential that
represents all of this in a, you know, a unique potential.
1231
01:55:25.860 --> 01:55:30.220
So there are no core electron
explicit in the pseudopotential,
1232
01:55:30.320 --> 01:55:33.080
so these two object, the real core
1233
01:55:33.180 --> 01:55:42.100
with the nucleus and the pseudopotential should
scatter an incoming plane-wave in similar ways.
1234
01:55:42.200 --> 01:55:44.040
And what does it mean in similar ways?
1235
01:55:44.140 --> 01:55:51.080
It means that you can take the plane-wave
and after the scattering you can do a Born
1236
01:55:51.180 --> 01:55:56.500
expansion in spherical waves and you can
represent the effect of the scattering
1237
01:55:56.600 --> 01:56:01.100
in terms of the phase shifts
of all these coefficients.
1238
01:56:01.200 --> 01:56:09.720
And the norm-conserving pseudopotential are built
to have, say, the correct phase shift,
1239
01:56:09.820 --> 01:56:14.340
that's the logarithmic derivative at one given energy,
1240
01:56:14.440 --> 01:56:16.360
but they are also transferable.
1241
01:56:16.460 --> 01:56:23.020
That means that they have the correct derivative
and having not only the absolute value,
1242
01:56:23.120 --> 01:56:30.220
but also the correct derivative means that those
phase shifts are going to be actually good
1243
01:56:30.320 --> 01:56:33.000
also in the neighborhood and, actually, you see,
1244
01:56:33.001 --> 01:56:39.280
you can make this phase shift to be good
in a very large neighbourhood of energy.
1245
01:56:39.380 --> 01:56:43.200
So that's why the norm-conservation concepts
1246
01:56:43.201 --> 01:56:47.671
of Hamann, Schluter and Chiang was so powerful.
1247
01:56:48.086 --> 01:56:56.284
Later on in 1985, David Vanderbilt extended
this concept and it developed a formalism
1248
01:56:56.484 --> 01:57:03.390
by which you could impose this condition
and not at only one value of energies,
1249
01:57:03.590 --> 01:57:05.400
but at multiple values.
1250
01:57:05.500 --> 01:57:09.180
So you could reproduce better and better
1251
01:57:09.280 --> 01:57:20.620
the scattering from the real core and the
pseudopotential in a wider region in energy.
1252
01:57:20.720 --> 01:57:23.960
And these are the extended norm-conserving
1253
01:57:24.060 --> 01:57:28.740
pseudopotential that are now becoming
more and more used because
1254
01:57:28.840 --> 01:57:33.600
you can make them more accurate than standard
norm-conserving because they have better
1255
01:57:33.700 --> 01:57:38.500
scattering properties
but are not so complicated as
1256
01:57:38.600 --> 01:57:42.680
what had been developed in between,
also by David Vanderbilt,
1257
01:57:42.681 --> 01:57:48.580
of the ultrasoft pseudopotentials or by Peter
Blochl, the projected augmented wave method.
1258
01:57:48.680 --> 01:57:52.260
These two approaches are somehow similar.
1259
01:57:52.360 --> 01:57:55.080
They tend to make the calculation much less
1260
01:57:55.180 --> 01:58:03.620
expensive, but they also tend to make
the entire algorithm more expensive.
1261
01:58:03.720 --> 01:58:07.780
I'm looking at my notes,
I think, today I'm going long.
1262
01:58:07.880 --> 01:58:13.740
So I'll mention I think we'll have 20 more minutes to go.
1263
01:58:13.840 --> 01:58:16.940
If I don't faint. If I faint, that's it,
1264
01:58:17.040 --> 01:58:24.500
if I don't faint, you know, and you also know
that all of this is recorded and will be posted.
1265
01:58:24.600 --> 01:58:30.580
So if you have to go, you can also
watch it, watch it later.
1266
01:58:30.680 --> 01:58:37.386
One thing that is interesting about the ultrasoft
pseudopotential, I just draw a graph here.
1267
01:58:37.387 --> 01:58:46.220
Suppose that this was the charge density,
the radial charge density of a 3d orbital,
1268
01:58:46.320 --> 01:58:50.040
well, the ultrasoft and similarly the projected
1269
01:58:50.140 --> 01:58:57.140
augmented wave, does very smart,
introduces some very smart formalism
1270
01:58:57.240 --> 01:59:03.986
to make the wavefunction, the pseudo
wavefunctions that we have in our calculation,
1271
01:59:04.114 --> 01:59:08.286
even smoother than what norm-
conserving pseudopotential can do.
1272
01:59:08.370 --> 01:59:10.720
And when they were introduced in the early 90s,
1273
01:59:10.820 --> 01:59:18.780
it means that materials that were very challenging
for plane-wave calculation, materials with 2p electrons
1274
01:59:18.880 --> 01:59:22.280
or 3d electrons... you see
2p electrons or 3d electrons
1275
01:59:22.380 --> 01:59:26.040
don't have to be orthogonal to 1p
1276
01:59:26.140 --> 01:59:29.900
because there is no 1p, or to 2d
because there is no 2d.
1277
01:59:30.000 --> 01:59:35.560
So 2p and 3d electrons are very
compressed towards the nucleus
1278
01:59:35.560 --> 01:59:37.120
so they need a lot of plane-waves.
1279
01:59:37.220 --> 01:59:41.220
They are very tough to describe and ultrasoft
1280
01:59:41.320 --> 01:59:48.940
pseudopotential manage to describe those
with very, very smooth pseudo wavefunction
1281
01:59:49.040 --> 01:59:56.843
but in that, what you lost is that the charge
density was not any more the square of the
1282
01:59:56.844 --> 02:00:01.580
wavefunction, but was actually
the square of the wavefunction
1283
02:00:01.680 --> 02:00:10.060
plus some so-called augmentation terms
because the charge density is a physical observable,
1284
02:00:10.160 --> 02:00:14.598
and so if I make my wavefunction
smoother than it should be,
1285
02:00:14.798 --> 02:00:18.000
I'm not going to have any more that the square
1286
02:00:18.100 --> 02:00:21.380
of the wavefunction gives me the true charge density.
1287
02:00:21.480 --> 02:00:25.580
I need to have these augmentation terms.
1288
02:00:25.680 --> 02:00:28.060
And so because of this,
1289
02:00:28.160 --> 02:00:37.137
it means that the relation between the cutoff
for the charge density, being 4 times the cutoff
1290
02:00:37.337 --> 02:00:45.438
for the wavefunction, has become broken.
The cutoff for the charge density cannot change
1291
02:00:45.638 --> 02:00:49.580
because of any computational trick that I do,
1292
02:00:49.680 --> 02:00:54.540
because that's a physical quantity and
that depends on the system I'm studying.
1293
02:00:54.640 --> 02:00:59.940
With the ultrasoft pseudopotential,
I can make the cutoff on the wavefunction
1294
02:01:00.040 --> 02:01:04.440
smaller and so it means that
the ratio between these two
1295
02:01:04.540 --> 02:01:09.500
is not any more 4, but it could
be 6, 7, 8, 10, 12.
1296
02:01:09.600 --> 02:01:14.260
It's something that, in principle, you also have to test.
1297
02:01:14.360 --> 02:01:17.420
Luckily, if you go to the Materials Cloud,
1298
02:01:17.520 --> 02:01:21.680
we have tested for you a standard set of pseudopotential.
1299
02:01:21.780 --> 02:01:30.843
We call them standard solid-state pseudopotential, SSSP,
and we have prepared them for the case that you want to
1300
02:01:30.844 --> 02:01:38.340
use PBE as your exchange-correlation potential
or PBEsol as your exchange-correlation potential.
1301
02:01:38.440 --> 02:01:43.673
No one has generated, as far as I know,
a set of pseudopotential for the SCAN functional.
1302
02:01:43.674 --> 02:01:45.085
That would be very good.
1303
02:01:45.086 --> 02:01:47.920
There are some, but not an entire periodic table.
1304
02:01:48.020 --> 02:01:53.983
So there are going to be inconsistency when you
do SCAN calculation between your pseudopotential
1305
02:01:54.183 --> 02:02:00.000
that has been generated using density functional
theory for the atom, but with an exchange-correlation
1306
02:02:00.100 --> 02:02:06.020
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.
1307
02:02:06.120 --> 02:02:11.214
But this is to say that, you know, given the
pseudopotential that you are using,
1308
02:02:11.215 --> 02:02:14.313
those can be hard, like if you use
1309
02:02:14.314 --> 02:02:16.400
norm-conserving pseudopotential
1310
02:02:16.500 --> 02:02:20.757
for oxygen, they might require a cutoff
of 80 Rydberg or 100 Rydberg.
1311
02:02:20.914 --> 02:02:24.420
If you use ultrasoft pseudopotential or PAW,
1312
02:02:24.520 --> 02:02:29.280
you might get away with 30 Rydberg,
and there is a huge difference in the cost
1313
02:02:29.380 --> 02:02:32.686
but you need to make sure that also your cutoff
1314
02:02:32.687 --> 02:02:40.520
of the charge density is correct. There is not going
to be a factor of four between the [inaudible].
1315
02:02:42.200 --> 02:02:44.960
Let me go to something else as well.
1316
02:02:45.060 --> 02:02:51.640
The case I wanted to delve a little bit
deeper in this concept of methods.
1317
02:02:51.740 --> 02:02:54.560
That is what happens when you have, actually,
1318
02:02:54.660 --> 02:03:00.820
a Fermi energy separating the occupied
states by the unoccupied states.
1319
02:03:00.920 --> 02:03:06.843
And, as I said, that makes our integrals,
that here have been discretised as
1320
02:03:06.844 --> 02:03:11.180
Monkhorst-Pack meshes, much more tricky
1321
02:03:11.181 --> 02:03:16.637
because now we have to say if
that state is occupied or not.
1322
02:03:16.638 --> 02:03:24.800
OK, so f is equal to one, occupied,
or is equal to zero if the state is unoccupied.
1323
02:03:25.280 --> 02:03:27.800
So this is really a sharp discontinuity.
1324
02:03:27.900 --> 02:03:29.980
I was doing my happy integral here,
1325
02:03:30.080 --> 02:03:35.420
but things dropped to zero and this
one cuts and this is very difficult.
1326
02:03:35.520 --> 02:03:45.780
So during the years there has been the concept
of trying to smooth these discontinuities
1327
02:03:45.880 --> 02:03:52.740
by introducing a fictitious temperature in the system,
you know. Suppose that we had a real temperature,
1328
02:03:52.840 --> 02:03:56.980
now, the occupation of my states
would be Fermi-Dirac occupations,
1329
02:03:57.080 --> 02:04:02.200
and while at zero temperatures,
the Fermi-Dirac is a step function that goes
1330
02:04:02.300 --> 02:04:08.500
from one or two, if you are, say, spin
degenerate to zero at finite temperature.
1331
02:04:08.600 --> 02:04:10.860
It's not. It's smoother,
1332
02:04:10.960 --> 02:04:16.186
and so it means that this, you know,
discontinuity has now become
1333
02:04:16.187 --> 02:04:17.643
much smoother.
1334
02:04:18.240 --> 02:04:21.220
So if you have a temperature in your system,
1335
02:04:21.320 --> 02:04:29.943
your functional is not any more the total
energy, but is the Helmholtz free-energy,
1336
02:04:33.680 --> 02:04:39.360
that is equal to the total energy
minus the temperature
1337
02:04:39.460 --> 02:04:41.020
times the entropy.
1338
02:04:41.120 --> 02:04:43.203
Well, this would be the... say, if I have
1339
02:04:43.204 --> 02:04:49.980
a Fermi-Dirac occupation would be
f log f, one minus f, log of one minus f.
1340
02:04:50.080 --> 02:04:57.530
So this is the proper extension to a finite
temperature formalism and often
1341
02:04:57.730 --> 02:05:02.860
the temperature, to remind ourselves
that it's not a real physical temperature,
1342
02:05:02.960 --> 02:05:05.680
we call it a smearing, sigma.
1343
02:05:05.780 --> 02:05:11.057
That is just there to smoothen our discontinuity.
1344
02:05:12.320 --> 02:05:19.580
Quantum Espresso, to confuse you, calls this
quantity, the total energy also for a metal,
1345
02:05:19.680 --> 02:05:22.557
but remember, it's not a total energy it's a
1346
02:05:22.586 --> 02:05:23.757
total....
1347
02:05:24.320 --> 02:05:26.843
total free energy.
1348
02:05:27.680 --> 02:05:35.160
And so, all of a sudden, what happens is that
now you have an additional thing to figure out.
1349
02:05:35.260 --> 02:05:42.540
That is how much temperature or how much smearing
you have to give to your system.
1350
02:05:42.640 --> 02:05:43.943
And things are delicate.
1351
02:05:43.944 --> 02:05:46.486
Here I am looking at a physical quantity.
1352
02:05:46.487 --> 02:05:53.843
I think this is an iron bulk system
with a few atoms in the primitive cell.
1353
02:05:53.844 --> 02:06:00.820
I am displacing one. I am measuring...
calculating the force on an atom.
1354
02:06:00.920 --> 02:06:04.860
So you see, what I am plotting here is
1355
02:06:04.960 --> 02:06:13.240
the force on an atom as a function
of the smearing. You know, .02 Rydberg
1356
02:06:13.360 --> 02:06:18.540
is equal to .27 electron volts,
1357
02:06:18.640 --> 02:06:25.340
that is roughly equal to... well, let's say 3000 K.
1358
02:06:25.440 --> 02:06:27.900
So if we think in physical terms,
1359
02:06:28.000 --> 02:06:34.186
you know, already .02 Rydberg of
smearing is an enormous temperature.
1360
02:06:35.440 --> 02:06:42.400
So what is the right result here?
You see, what I'm seeing is that the force that I
1361
02:06:42.500 --> 02:06:51.871
get changes if I change the smearing, and
changes if I change the k-point sampling.
1362
02:06:52.840 --> 02:06:57.000
What I haven't shown you is that if we go to very,
1363
02:06:57.100 --> 02:07:05.560
very high smearing, we would always get
one curve, because at very high smearing,
1364
02:07:05.920 --> 02:07:13.020
my discontinuities become so smooth
that I can integrate it very well.
1365
02:07:13.120 --> 02:07:15.420
And you see a remnant of this here.
1366
02:07:15.520 --> 02:07:18.560
At .02, you have three or four curves
1367
02:07:18.660 --> 02:07:26.140
that are almost identical, the ones with
10 or 20, 24, 36 point meshes.
1368
02:07:26.240 --> 02:07:34.940
But as I decrease the smearing, you know,
all this curve starts to become different
1369
02:07:35.040 --> 02:07:43.770
Even the best curve that I have, the best curve
is the green diamonds here and the black dots here.
1370
02:07:43.771 --> 02:07:48.540
You see these two best curves are identical
1371
02:07:48.640 --> 02:07:55.380
in the high temperature, high smearing
region, and then they start to separate.
1372
02:07:55.480 --> 02:07:59.460
And as you go to very small smearing,
things go actually berserk.
1373
02:07:59.560 --> 02:08:06.800
So it's not that, you know, getting
to zero smearing, gives you a good result.
1374
02:08:06.900 --> 02:08:09.460
You are here. You are all over the place.
1375
02:08:09.560 --> 02:08:13.380
What you would want ideally is the limit
1376
02:08:13.480 --> 02:08:20.460
of infinite k-point sampling and then
zero smearing, but you cannot reach that.
1377
02:08:20.560 --> 02:08:28.228
So, you know, what you need to figure out roughly...
let me remove these very confusing pen marks...
1378
02:08:28.229 --> 02:08:35.259
is finding, if you want the plateau around here
and, you know, you need to look at what is,
1379
02:08:35.459 --> 02:08:38.040
you know, the error that you can tolerate here,
1380
02:08:38.140 --> 02:08:43.040
and this is a matter of, you know,
around .001, .002 electron volt
1381
02:08:43.140 --> 02:08:46.500
so, it's not that big error.
1382
02:08:46.600 --> 02:08:50.080
But you need to look at the plateau where,
1383
02:08:50.180 --> 02:08:55.620
you know, with the highest sampling,
you roughly have the same results.
1384
02:08:55.720 --> 02:08:57.360
And that's a function of smearing
1385
02:08:57.460 --> 02:09:02.560
things are not changing too much.
If you go to larger smearings
1386
02:09:02.660 --> 02:09:05.540
it's easier to converge
1387
02:09:05.640 --> 02:09:10.500
but you converge to the wrong result
because you have a lot of temperature.
1388
02:09:10.600 --> 02:09:12.740
If you go to smaller smearing,
1389
02:09:12.840 --> 02:09:17.220
you have the right temperature,
but your sampling error is killing you.
1390
02:09:17.320 --> 02:09:21.457
So here, sampling is killing you.
1391
02:09:24.480 --> 02:09:26.280
And here,
1392
02:09:26.880 --> 02:09:29.786
temperature is killing you.
1393
02:09:31.600 --> 02:09:38.980
It's like being in a place where
lava falls into ice water,
1394
02:09:39.080 --> 02:09:44.940
You need to find that comfort zone and it
is very thin in between these two.
1395
02:09:45.040 --> 02:09:49.520
But, you know, to make things a little bit
1396
02:09:49.680 --> 02:09:51.180
more reassuring,
1397
02:09:51.280 --> 02:09:56.740
I would say that, you know, this interval of
smearing tends to be conservatively good.
1398
02:09:56.840 --> 02:10:03.260
So never go below, never go above
and make sure you have enough k-points.
1399
02:10:03.360 --> 02:10:13.580
And the reason why this is actually a good range
is because we have developed 20, 25 years ago
1400
02:10:13.680 --> 02:10:19.040
what are called generalised
entropic formulation that I'll just
1401
02:10:19.140 --> 02:10:25.280
briefly go through. That is we
have been able to construct
1402
02:10:25.440 --> 02:10:32.940
formulation that are not leading to a Fermi-Dirac
entropy, are leading to a different entropy.
1403
02:10:33.040 --> 02:10:36.260
And they are such that A,
1404
02:10:36.360 --> 02:10:44.460
the free energy as a function of the fictitious
temperature sigma, is actually quite flat.
1405
02:10:44.560 --> 02:10:51.571
That means that if A is flat as a function
of temperature, I can do my calculation
1406
02:10:51.572 --> 02:10:56.820
at high temperature where my
k-point sampling is converged,
1407
02:10:56.920 --> 02:11:02.700
but this number that I get is actually very,
very close to the number I can get at zero temperature.
1408
02:11:02.800 --> 02:11:03.980
How do we do that?
1409
02:11:04.080 --> 02:11:08.500
Well, we write the generalised Helmholtz free energy,
1410
02:11:08.600 --> 02:11:14.420
and we look at the minimum
conditions of these functionals.
1411
02:11:14.520 --> 02:11:16.800
So I minimise it.
1412
02:11:16.900 --> 02:11:24.660
I have two Lagrange multipliers on the
total number of electrons on orthogonality
1413
02:11:24.760 --> 02:11:30.080
and I have the derivatives, functional
derivatives, with respect to the wavefunction
1414
02:11:30.180 --> 02:11:34.100
and I have the derivative with respect to the occupations.
1415
02:11:34.200 --> 02:11:36.700
And in particular, these derivatives
1416
02:11:36.800 --> 02:11:42.400
gives me a constructive relation between
the expectation value of the Hamiltonian,
1417
02:11:42.500 --> 02:11:48.220
that's not the eigenvalues, the expectation
value of the Hamiltonian on the psi i
1418
02:11:48.320 --> 02:11:55.120
minus the Fermi energy, the chemical
potential, mu, and this is equal to T times
1419
02:11:55.220 --> 02:11:58.640
the derivative of the entropy
with respect to the occupation.
1420
02:11:58.740 --> 02:12:01.680
So this can be written also
1421
02:12:02.129 --> 02:12:08.614
epsilon i minus mu divided by T
needs to be equal to the derivative of the
1422
02:12:09.720 --> 02:12:13.114
entropy with respect to the occupation.
1423
02:12:13.400 --> 02:12:15.557
So from this
1424
02:12:16.760 --> 02:12:24.720
we can do something very simple, that is,
you know, these ideas of smearing,
1425
02:12:24.820 --> 02:12:32.860
this idea of an occupation being
a smooth function of the temperature,
1426
02:12:32.960 --> 02:12:39.180
And this variable is epsilon minus mu versus T
was introduced already in the mid eighties by Fu and Ho
1427
02:12:39.280 --> 02:12:45.340
that said, well, you know, let's take a broadening
of a Dirac delta, it could be a Gaussian,
1428
02:12:45.440 --> 02:12:49.720
and we calculated the occupation as this integral.
1429
02:12:49.820 --> 02:12:53.260
You see, for epsilon is very large,
1430
02:12:53.360 --> 02:12:58.600
this integral is going to give me one...
Sorry, this is...
1431
02:12:58.920 --> 02:13:01.800
Yeah, so.
1432
02:13:03.600 --> 02:13:07.940
If epsilon is very large,
this integral is going to give me zero,
1433
02:13:08.040 --> 02:13:14.060
and if epsilon is very negative,
this integral is going to give me one.
1434
02:13:14.160 --> 02:13:16.120
If epsilon is equal to mu this integral
1435
02:13:16.220 --> 02:13:22.200
is going to give me one half, and,
you know, they proposed different forms.
1436
02:13:22.300 --> 02:13:25.529
Let me not get into this now.
1437
02:13:27.120 --> 02:13:28.640
And...
1438
02:13:32.840 --> 02:13:37.020
So, remember what we had before,
1439
02:13:37.120 --> 02:13:44.600
so we had this constitutive relation between the
expectation value of the Hamiltonian and the entropy.
1440
02:13:45.120 --> 02:13:48.860
We have defined the occupation as an integral
1441
02:13:48.960 --> 02:13:54.760
of a broadening function g that we decide.
This is a broadening function.
1442
02:13:54.860 --> 02:13:59.729
And there is a very beautiful
relation that Alessandro De Vita
1443
02:13:59.730 --> 02:14:05.920
and Mike Gillan developed in the early 90s
that then tells me, well, if that is the case,
1444
02:14:06.020 --> 02:14:12.060
so if f is written as an integral of a broaden
1445
02:14:12.160 --> 02:14:19.920
occupation, the entropy is the integral of minus
1446
02:14:20.120 --> 02:14:24.420
whatever variables t times
the broaden occupation g(t).
1447
02:14:24.520 --> 02:14:28.020
So there is if you want a very simple way to go
1448
02:14:28.120 --> 02:14:34.480
from broadening g to entropy S so
you can choose your broadening.
1449
02:14:34.580 --> 02:14:40.380
You can take your Dirac delta and you can
broaden it in very many different ways.
1450
02:14:40.480 --> 02:14:43.920
And this construction will tell you what is
1451
02:14:44.020 --> 02:14:52.020
the form of the entropy that needs to go
into your Helmholtz free energy here.
1452
02:14:52.120 --> 02:14:56.840
And so we use this concept to build
1453
02:14:57.960 --> 02:15:01.140
broadening function
1454
02:15:01.240 --> 02:15:05.520
delta tilde that had a number of quality,
1455
02:15:05.620 --> 02:15:12.420
that is, of course, is normalised,
that for which this integral was giving zero.
1456
02:15:12.520 --> 02:15:16.420
That means that the entropy
at zero temperature is zero
1457
02:15:16.520 --> 02:15:19.860
but in particular, we built
1458
02:15:19.960 --> 02:15:28.620
distribution functions, broader delta
that add a second moment that is equal
1459
02:15:28.720 --> 02:15:34.440
to zero and adding this second moment equal
to zero was demonstrated by Stefano de Gironcoli
1460
02:15:34.540 --> 02:15:43.200
that is equivalent to having the free energy,
not be, say linear or quadratic as a function
1461
02:15:43.300 --> 02:15:51.300
of temperature, but actually to have only cubic
or quartic terms as a function of temperature.
1462
02:15:51.400 --> 02:15:54.440
And this is an example of the broadening function.
1463
02:15:54.540 --> 02:15:57.620
This is called cold smearing in Quantum Espresso.
1464
02:15:57.720 --> 02:16:00.720
There are other variation on the theme.
1465
02:16:00.820 --> 02:16:07.500
There is the Methfessel and Paxton who's actually
introduced first and doesn't go through this.
1466
02:16:07.600 --> 02:16:12.540
If you want conceptual derivation
of a free energy functional,
1467
02:16:12.640 --> 02:16:17.680
the Methfessel-Paxton has slight defects
that can give rise to negative occupancies.
1468
02:16:17.780 --> 02:16:20.800
That is something that we do not like but both
1469
02:16:20.900 --> 02:16:25.820
of these broadening distribution
have a zero second moment.
1470
02:16:25.920 --> 02:16:33.840
That means that the free energy does not
depend on temperature up to the cubic order.
1471
02:16:33.940 --> 02:16:42.871
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
1472
02:16:42.914 --> 02:16:48.180
and you see how the free energy
is quadratic in temperature.
1473
02:16:48.280 --> 02:16:53.460
The total energy is actually quadratic
with the opposite coefficient.
1474
02:16:53.560 --> 02:16:59.320
It was realised very early on that if we were to take
the semi-sum, this is what Mike Gillan called
1475
02:16:59.420 --> 02:17:06.360
the corrected energy, E minus 1/2 T S,
we would have, you see, with very large smearing...
1476
02:17:06.460 --> 02:17:12.700
now these are electron volts... a very good
estimate of the zero value term.
1477
02:17:12.800 --> 02:17:19.940
Remember that this is the no-go zone
where we can't do good sampling in there.
1478
02:17:20.040 --> 02:17:24.460
But this is just good for the total energy.
1479
02:17:24.560 --> 02:17:28.020
The functional is still the free energy,
1480
02:17:28.120 --> 02:17:33.180
so the stress and the forces that would
come out of would not be corrected,
1481
02:17:33.280 --> 02:17:38.143
but the generalised smearing formulation,
like the cold smearing, the Methfessel-Paxton,
1482
02:17:38.144 --> 02:17:38.900
do both.
1483
02:17:38.901 --> 02:17:46.960
They both give you a free energy that in the case
of the cold smearing, say, you see, it goes like this.
1484
02:17:47.060 --> 02:17:52.900
You see how much flatter it is than the free energy.
1485
02:17:53.329 --> 02:17:55.700
The free energy, you see it's exploding,
1486
02:17:55.800 --> 02:18:00.380
and remember, this is the region
where we are doing our calculations.
1487
02:18:00.480 --> 02:18:02.980
But because it's a variational functional,
1488
02:18:03.080 --> 02:18:10.243
we can calculate also forces and stresses
and all derived quantities.
1489
02:18:11.200 --> 02:18:13.660
So let me summarise, what are, you know,
1490
02:18:13.760 --> 02:18:19.540
all the ingredients that you'll find in your
homework to do a self-consistent calculation.
1491
02:18:19.640 --> 02:18:23.800
You need to figure out what is the
translational symmetry of your system?
1492
02:18:23.900 --> 02:18:25.580
What is your Bravais lattice?
1493
02:18:25.680 --> 02:18:32.780
What is the group of atoms that we call
the basis for the crystal structure?
1494
02:18:32.880 --> 02:18:36.560
That is, the atoms that are periodically repeated
1495
02:18:36.660 --> 02:18:41.640
by the Bravais translation, and they
will give rise to the entire real crystal.
1496
02:18:41.641 --> 02:18:45.340
So for silicon, you would have a basis of two atoms.
1497
02:18:45.440 --> 02:18:51.771
For gallium arsenide, also a basis
of two atoms and an FCC Bravais lattice.
1498
02:18:51.857 --> 02:18:56.020
For aluminum, we would have one
atom in the FCC Bravais lattice.
1499
02:18:56.120 --> 02:19:01.260
For BCC iron, you would have one atom
and then BCC Bravais lattice,
1500
02:19:01.360 --> 02:19:04.280
but, of course, if you want to study a defect
1501
02:19:04.380 --> 02:19:10.843
in silicon, then you might need to construct
a much bigger primitive unit cell with maybe
1502
02:19:10.844 --> 02:19:15.300
sixty four atoms where then you
remove one of these atoms.
1503
02:19:15.400 --> 02:19:18.660
So you need to construct your crystals.
1504
02:19:18.760 --> 02:19:21.471
You need to get your pseudopotentials from somewhere.
1505
02:19:21.472 --> 02:19:27.458
From the SSSP, from the Quantum Espresso
website, from the PseudoDojo website.
1506
02:19:27.658 --> 02:19:33.060
You need to find out what is the correct
cutoff energy for your wavefunction.
1507
02:19:33.160 --> 02:19:36.460
The SSSP website gives you suggestions.
1508
02:19:36.560 --> 02:19:40.020
What is the correct cutoff for the charge density?
1509
02:19:40.120 --> 02:19:43.529
Four times larger for now, concerning pseudopotential
1510
02:19:43.530 --> 02:19:48.100
but may be more and more for ultrasoft pseudopotentials.
1511
02:19:48.200 --> 02:19:51.386
What is the appropriate Monkhorst-Pack mesh?
1512
02:19:52.033 --> 02:19:55.580
With this, the fictitious temperature that you can use.
1513
02:19:55.680 --> 02:19:57.760
And then you need to make sure that you have
1514
02:19:57.860 --> 02:20:03.580
a good technique to get to self-consistency.
1515
02:20:03.680 --> 02:20:06.700
And I'll discuss that in a moment.
1516
02:20:06.800 --> 02:20:09.966
If you are desperate, you can go to the Materials Cloud.
1517
02:20:09.967 --> 02:20:15.271
In the the work tool, there is something called
the Quantum Espresso input generator
1518
02:20:15.286 --> 02:20:18.460
that allows you to very simple...
1519
02:20:18.560 --> 02:20:26.180
very simply upload a structure, choose a table
of pseudopotential for PBE or PBEsol,
1520
02:20:26.280 --> 02:20:33.540
either efficiency, that is less accurate
but less expensive, or accuracy, more accurate.
1521
02:20:33.640 --> 02:20:36.820
And then you just have to decide if your system is
1522
02:20:36.920 --> 02:20:42.620
insulating or metallic and if
it's magnetic or non-magnetic.
1523
02:20:42.720 --> 02:20:48.800
And you decide one of the standard three
meshes of k-points and Quantum Espresso
1524
02:20:48.900 --> 02:20:55.520
generates... the Quantum Espresso input
generator generates for you the...
1525
02:20:55.760 --> 02:20:57.180
the input file.
1526
02:20:57.280 --> 02:21:03.100
And if you want to plot band structure,
there is another tool called SeeK-path,
1527
02:21:03.200 --> 02:21:11.314
that gives you also the appropriate crystallographic
approved path in the Brillouin zone.
1528
02:21:11.920 --> 02:21:15.720
So let me conclude for today...
I think I'm going a little bit longer.
1529
02:21:15.820 --> 02:21:18.420
let's try to do this in eight minutes...
1530
02:21:18.520 --> 02:21:26.900
with how to get to self-consistence and exactly
as in Hartree and Hartree-Fock, in Kohn-Sham.
1531
02:21:27.000 --> 02:21:31.420
We need to get to the self-consistent
solution of these equations.
1532
02:21:31.520 --> 02:21:38.380
So, you know, once we set up our problem
with all these computational parameters,
1533
02:21:38.480 --> 02:21:43.740
we need to construct, you know, the Hartree
operator, the exchange-correlation operator.
1534
02:21:43.840 --> 02:21:50.500
So we need to, you know, have actually a starting
electronic density to then iterate,
1535
02:21:50.600 --> 02:21:56.387
and so this is typically given by sums of atomic densities.
1536
02:21:56.587 --> 02:22:03.180
At these points we can construct the Hartree
operator, the exchange-correlation operator.
1537
02:22:03.280 --> 02:22:06.820
So we have the Hamiltonian. For that Hamiltonian,
1538
02:22:06.920 --> 02:22:09.300
we solve the Kohn-Sham equation.
1539
02:22:09.400 --> 02:22:14.380
This Kohn-Sham equation
gives us a set of, you know, phi 1,
1540
02:22:14.480 --> 02:22:20.000
phi N for a molecule or phi 1k,
1541
02:22:20.120 --> 02:22:26.820
phi Nk for a solid with k belonging
to the Monkhorst-Pack mesh.
1542
02:22:26.920 --> 02:22:28.880
Now with these
1543
02:22:29.240 --> 02:22:33.400
phis, we calculate the new charge density.
1544
02:22:33.500 --> 02:22:35.740
This was maybe rho 0.
1545
02:22:35.840 --> 02:22:38.860
Now we have rho 1 and with rho 1
1546
02:22:38.960 --> 02:22:42.460
we construct the Hartree
and exchange-correlation operator.
1547
02:22:42.560 --> 02:22:48.460
We solve, we get a rho 2. We iterate, iterate, iterate
1548
02:22:48.560 --> 02:22:52.343
until things stop changing.
1549
02:22:52.771 --> 02:22:59.143
And, you know, this process
can get to self consistency.
1550
02:22:59.271 --> 02:23:04.220
More often than not, this process
would actually not converge
1551
02:23:04.320 --> 02:23:08.440
but you can get a little bit smarter
and introduce a damping term.
1552
02:23:08.540 --> 02:23:10.340
You can say that
1553
02:23:10.440 --> 02:23:17.880
the new charge density that enters
into this loop is not equal to the old
1554
02:23:18.720 --> 02:23:24.657
but is actually just equal to a fraction of the old
1555
02:23:24.771 --> 02:23:30.400
plus the remaining of the previous one.
1556
02:23:30.500 --> 02:23:37.740
So you sort of mix the charge density
in an iterative way by making beta very small.
1557
02:23:37.840 --> 02:23:41.820
You slow down the process of
1558
02:23:41.920 --> 02:23:49.860
updating what was your initial starting condition,
so you much more smoothly move around
1559
02:23:49.960 --> 02:23:53.340
and if you're lucky,
it becomes easier to converge,
1560
02:23:53.440 --> 02:23:56.020
although you might also be slower.
1561
02:23:56.120 --> 02:24:02.460
And this is actually what is
taken from the review of modern physics
1562
02:24:02.560 --> 02:24:09.120
of Mike Payne on the total energy pseudopotential
plane-wave method from 1992.
1563
02:24:09.220 --> 02:24:14.529
This is exactly what I said, just in graphical form.
1564
02:24:15.760 --> 02:24:19.780
So but, you know, if you don't converge
1565
02:24:19.880 --> 02:24:27.620
and you decrease your beta and nothing happens,
you just might be unlucky in your system.
1566
02:24:27.720 --> 02:24:29.440
You know, it's very tough to converge.
1567
02:24:29.540 --> 02:24:32.940
Maybe it's very large, maybe
it has sloshing instability.
1568
02:24:33.040 --> 02:24:36.940
Maybe it has two levels, one empty
and one full that, you know,
1569
02:24:37.040 --> 02:24:42.020
keep oscillating back and forth and
you know which one to fill.
1570
02:24:42.120 --> 02:24:47.960
And so I wanted to introduce the complementary
approach that maybe might be less efficient,
1571
02:24:48.060 --> 02:24:53.320
but is much more robust to converge
total energy plane-wave calculation
1572
02:24:53.420 --> 02:24:58.314
that is based on variational recipes.
1573
02:24:58.443 --> 02:25:01.940
That is the idea that I've given before that rather than
1574
02:25:02.040 --> 02:25:06.620
iterating to self-consistency, the
Hamiltonian and the Kohn-Sham equation,
1575
02:25:06.720 --> 02:25:11.420
you could do a direct nonlinear
minimisation of the energy functional.
1576
02:25:11.520 --> 02:25:16.440
And to give you an example, I've written
here the kinetic energy for a simple case.
1577
02:25:16.540 --> 02:25:17.880
I haven't put k-points.
1578
02:25:17.980 --> 02:25:26.540
I've just put, say, orbitals with just a discrete
index n expanded in plane-waves.
1579
02:25:26.640 --> 02:25:32.020
And, you know, in the energy,
as I already told you, there is,
1580
02:25:32.120 --> 02:25:39.260
you know, for the every plane-wave,
we have a one half G square component.
1581
02:25:39.360 --> 02:25:45.800
And so the expression of the kinetic energy for a
1582
02:25:46.440 --> 02:25:49.440
given orbital and then summed over all
1583
02:25:49.540 --> 02:25:53.920
the orbital is, you see, sum over all
the orbitals, one half sum over all
1584
02:25:54.020 --> 02:25:58.520
the plane-waves, the square modulus
of the coefficient of the plane-waves.
1585
02:25:58.620 --> 02:26:03.260
So the kinetic energy is quadratic
in the coefficient of the plane-waves.
1586
02:26:03.360 --> 02:26:08.280
But there is not only the kinetic energy, there is
also the potential energies in the system.
1587
02:26:08.380 --> 02:26:16.340
And now the potential energy is not diagonal like
the kinetic energy in the basis of plane-waves.
1588
02:26:16.440 --> 02:26:22.260
So the matrix element between G and G prime
is not zero as in the kinetic energy.
1589
02:26:22.360 --> 02:26:30.520
So the Hamiltonian is... has kinetic energy
on the diagonal and potential energy everywhere.
1590
02:26:30.620 --> 02:26:35.814
This is the Hamiltonian in a basis of plane-waves
1591
02:26:36.640 --> 02:26:42.560
and the matrix element of the
potential are nothing else.
1592
02:26:42.660 --> 02:26:51.260
You see, this is just the coefficient of the
Fourier transform of the potential in real space
1593
02:26:51.360 --> 02:26:56.860
with respect to a plane-wave
exponential of i G prime minus G.
1594
02:26:56.960 --> 02:27:00.800
So if we put the two together and if for a moment
1595
02:27:00.900 --> 02:27:07.120
the Hamiltonian and the total energy of the system
was just given by a kinetic energy term
1596
02:27:07.220 --> 02:27:11.660
and a potential term without considering
that in Kohn-Sham we have
1597
02:27:11.760 --> 02:27:16.260
the external potential, the Hartree potential,
the exchange-correlation potential.
1598
02:27:16.360 --> 02:27:19.820
But you see, there is a little star here that's difficult to see.
1599
02:27:19.920 --> 02:27:23.160
This would be the expression of the energy and we
1600
02:27:23.260 --> 02:27:27.220
can minimise it with respect
to all these coefficients.
1601
02:27:27.320 --> 02:27:29.460
So finding the ground state
1602
02:27:29.560 --> 02:27:37.471
can be seen not as an iterative process,
but as a minimisation of this form that is ...
1603
02:27:37.472 --> 02:27:43.786
here it looks like a quadratic form on the
coefficients, cG, but whilst we have
1604
02:27:43.787 --> 02:27:50.140
not only the, say, external potential,
but we have the external, the Hartree
1605
02:27:50.240 --> 02:27:55.520
and the exchange-correlation in this
potential, the charge density enters
1606
02:27:55.620 --> 02:27:59.360
so this coefficient enters inside the Fourier
1607
02:27:59.460 --> 02:28:05.460
transform of the potential, making this
a non-linear minimisation problem.
1608
02:28:05.560 --> 02:28:14.043
But, you know, it's still a problem that, you know, people
have studied and so we can try to minimise this stuff
1609
02:28:14.086 --> 02:28:18.200
in the space of all the possible coefficients.
1610
02:28:18.300 --> 02:28:23.840
So I've written it here again,
my approximate, here, energy function because I
1611
02:28:23.940 --> 02:28:29.880
don't have Hartree and exchange
correlation so I can minimise it in the space
1612
02:28:29.980 --> 02:28:35.040
of all these coefficients, so it's a
a multidimensional minimisation.
1613
02:28:35.140 --> 02:28:36.940
Remember that you might have
1614
02:28:37.040 --> 02:28:43.100
thousands or tens of thousands or
hundreds of thousands of coefficients.
1615
02:28:43.200 --> 02:28:47.540
How do I minimise in a space of coefficients?
1616
02:28:47.640 --> 02:28:50.760
Well, maybe I'm somewhere here
in the space of coefficient,
1617
02:28:50.860 --> 02:28:53.320
that is, I have certain values and I need
1618
02:28:53.420 --> 02:29:01.060
to evolve to get to the minimum of the energy
in this multidimensional space of, you know,
1619
02:29:01.160 --> 02:29:07.020
10,000 dimensions and each
coefficient has its own dimension.
1620
02:29:07.120 --> 02:29:09.740
And I need to get to the global minimum.
1621
02:29:09.840 --> 02:29:14.280
But the process is similar to what
you would see in, say, this case,
1622
02:29:14.380 --> 02:29:19.560
one-dimensional, two-dimensional system,
and maybe we have a two-dimensional space
1623
02:29:19.660 --> 02:29:29.600
of coefficients and we want to get to the minimum
of this energy in the space of c1 and c2.
1624
02:29:30.200 --> 02:29:31.360
There are a lot of techniques
1625
02:29:31.460 --> 02:29:37.857
but one thing that you want to know is what is
the derivative of the energy with respect to c1 or c2
1626
02:29:37.858 --> 02:29:43.340
because that gives you the right
direction, the right gradient, where to move.
1627
02:29:43.440 --> 02:29:46.860
So in general, this is written in a compact form.
1628
02:29:46.960 --> 02:29:48.880
We need to know what is the derivative
1629
02:29:48.980 --> 02:29:53.240
of the energy. In this case,
a functional derivative with respect...
1630
02:29:53.340 --> 02:30:00.440
maybe these are just real coefficients
and, actually, it turns out these derivatives is just
1631
02:30:00.540 --> 02:30:03.780
given by the application
of the Hamiltonian to psi i.
1632
02:30:03.880 --> 02:30:05.580
So if I have a wavefunction
1633
02:30:05.680 --> 02:30:12.210
the application of the Hamiltonian gives me
the direction where to move minus.
1634
02:30:12.410 --> 02:30:18.740
And so I can start to find ways
with these derivatives to go to the bottom.
1635
02:30:18.840 --> 02:30:22.700
And there are sort of two main
ways to go to the bottom.
1636
02:30:22.800 --> 02:30:26.280
One is to say, well, if I'm here and I feel
1637
02:30:26.380 --> 02:30:32.580
a force like this, I can use
that force as an acceleration.
1638
02:30:32.680 --> 02:30:40.240
So, you know, I'll start moving in the space
of coefficients towards the bottom,
1639
02:30:40.320 --> 02:30:41.957
but because I'm feeling the force of
1640
02:30:41.958 --> 02:30:49.060
some acceleration, by the time I'm here, I have
zero acceleration, but I have still a velocity.
1641
02:30:49.160 --> 02:30:53.700
And so, you know, without any damping,
I would oscillate forever.
1642
02:30:53.800 --> 02:31:02.220
And so I need to put some damping here to make
sure that at the end I converge to this ground state.
1643
02:31:02.320 --> 02:31:06.060
The other possibility would be to say that
1644
02:31:06.160 --> 02:31:08.640
I instead move with a velocity,
1645
02:31:08.740 --> 02:31:13.300
this is, you see, acceleration versus velocity.
One dot or two dots.
1646
02:31:13.371 --> 02:31:20.180
So if I have a force that is proportional
to the velocity, again I'll start moving,
1647
02:31:20.280 --> 02:31:23.560
but a force that is proportional to the velocity
1648
02:31:23.660 --> 02:31:29.260
means that I have critical slowing down.
As I get closer and closer to the minimum.
1649
02:31:29.360 --> 02:31:33.660
I never reach it because
my driving force goes to zero.
1650
02:31:33.760 --> 02:31:37.020
So the first recipe that we call Car-Parrinello
1651
02:31:37.120 --> 02:31:41.800
and the second recipe that we call conjugate
gradients have the defect, one,
1652
02:31:41.900 --> 02:31:45.900
that Car-Parrinello overshoots
the minimum. That has to be damped down.
1653
02:31:46.000 --> 02:31:51.560
Conjugate gradient has a critical slowing down
and it has to be pushed to get to the minimum.
1654
02:31:51.660 --> 02:31:55.820
But these are the two main techniques:
1655
02:31:55.920 --> 02:32:06.420
conjugate gradient comes again from some
mathematical analysis in which if you have, say,
1656
02:32:06.520 --> 02:32:08.700
a potential energy surface,
1657
02:32:08.800 --> 02:32:15.320
that is very asymmetric and you use
the gradients to get to the bottom.
1658
02:32:15.420 --> 02:32:21.920
So you see here, I'm moving along the gradient
until the derivative along that direction is zero.
1659
02:32:22.020 --> 02:32:26.700
Then I move along the gradient and move
along the gradient. I zigzag to the bottom
1660
02:32:26.800 --> 02:32:33.020
and it's much better once I've taken the first
step not to move in the
1661
02:32:33.120 --> 02:32:37.920
steepest descent direction,
but to move in what is called a conjugate
1662
02:32:38.020 --> 02:32:43.360
direction that takes the steepest direction,
but tries to make sure that as you move
1663
02:32:43.460 --> 02:32:47.700
in the direction, you don't lose the
minimisation that you have done before.
1664
02:32:47.701 --> 02:32:49.100
You see fine, we're here.
1665
02:32:49.200 --> 02:32:51.420
And moving down the steepest direction,
1666
02:32:51.520 --> 02:32:53.700
it's true that I would improve around that,
1667
02:32:53.800 --> 02:32:57.640
but I would do more so with respect to what
I done on the previous iteration.
1668
02:32:57.740 --> 02:32:59.880
So conjugate gradients...
1669
02:33:00.120 --> 02:33:03.420
Thus the conjugate gradients techniques
1670
02:33:03.520 --> 02:33:12.520
have what have been used in a lot of codes,
CASTEP, especially, to get in a very robust way
1671
02:33:12.620 --> 02:33:17.100
to the ground state and thus give
that when you have to the ground state,
1672
02:33:17.200 --> 02:33:20.580
we'll see that tomorrow, you can calculate forces,
1673
02:33:20.680 --> 02:33:27.300
and then maybe you let atoms move
using molecular dynamics techniques.
1674
02:33:29.029 --> 02:33:29.871
But
1675
02:33:30.000 --> 02:33:31.020
exactly
1676
02:33:31.120 --> 02:33:33.260
if I'm
1677
02:33:33.360 --> 02:33:38.020
looking at a system during a time evolution.
1678
02:33:38.120 --> 02:33:45.380
So if the atoms move, every time an atom move
I have a new potential energy surface,
1679
02:33:45.480 --> 02:33:47.560
and so you have to minimize.
I have to minimize.
1680
02:33:47.660 --> 02:33:48.940
I have to minimize.
1681
02:33:49.040 --> 02:33:52.720
And so getting to the minimum
1682
02:33:53.360 --> 02:33:57.100
at every step is an expensive operation.
1683
02:33:57.760 --> 02:33:59.600
And that's the idea that came
1684
02:33:59.700 --> 02:34:02.859
to Car and Parrinello in the mid 1980s,
1685
02:34:02.860 --> 02:34:09.386
that is to say, well, let's try to find a way
by which we don't have to reach the minimum
1686
02:34:09.443 --> 02:34:16.140
of this, you know, total energy as the function
of the coefficient space at every step,
1687
02:34:16.240 --> 02:34:22.440
but with each of the coefficient,
how to evolve themselves so that they follow,
1688
02:34:22.540 --> 02:34:26.140
maybe not exactly, but they follow very closely
1689
02:34:26.141 --> 02:34:32.600
what would be their minimum at every
instant, at every time step.
1690
02:34:32.601 --> 02:34:35.820
So this could be time zero,
this could be time zero plus delta T,
1691
02:34:35.920 --> 02:34:40.700
this could be time zero plus two
delta T and so on and so forth.
1692
02:34:40.800 --> 02:34:46.460
And Car and Parrinello did that through the idea
of extended Lagrangian, that is,
1693
02:34:46.560 --> 02:34:52.580
rather than having a standard Lagrangian
with the kinetic energy minus
1694
02:34:52.680 --> 02:34:59.820
the potential energy of our system, that is
the electronic energy, the Kohn-Sham energy.
1695
02:34:59.920 --> 02:35:03.780
What they did is that they introduced
1696
02:35:03.880 --> 02:35:07.500
dynamics, also for the coefficients
1697
02:35:07.600 --> 02:35:13.000
of, sorry, of the orbitals. You see here
there are the velocity of the orbital.
1698
02:35:13.100 --> 02:35:18.460
So this is like a physical
kinetic energy for the atoms
1699
02:35:18.560 --> 02:35:24.580
and this is a fictitious kinetic energies
for the plane-wave coefficients.
1700
02:35:24.680 --> 02:35:29.720
But by making the mass of these
fictitious plane-wave coefficients
1701
02:35:29.820 --> 02:35:36.180
that are now thrown in into my dynamical
algorithm by making this mass very, very small.
1702
02:35:36.280 --> 02:35:42.860
It's like making them the flies in the cow.
The cow that are the atoms can move,
1703
02:35:42.960 --> 02:35:47.900
and these coefficients don't affect the dynamics,
1704
02:35:48.000 --> 02:35:55.360
but just follows the atoms and very quickly
rearrange themselves in that blue
1705
02:35:55.460 --> 02:35:59.340
line at the bottom of the
evolving potential surface.
1706
02:35:59.440 --> 02:36:05.500
So this would be the extended Car-Parrinello
Lagrangian. It gives rise to equation of motion,
1707
02:36:05.600 --> 02:36:11.100
not only for the atoms, but also for the
coefficients that learn how to follow.
1708
02:36:11.200 --> 02:36:13.980
And once you look at the,
1709
02:36:14.080 --> 02:36:18.160
actually, constant of motions,
you discover that because you are
1710
02:36:18.161 --> 02:36:23.229
integrating the extended Lagrangian,
you have a constant of motion where the
1711
02:36:23.230 --> 02:36:28.140
kinetic energy of the atom sees
the potential energy of the atom is.
1712
02:36:28.240 --> 02:36:30.220
But there is also these
1713
02:36:30.320 --> 02:36:33.940
kinetic energies of the coefficient of the plane-waves,
1714
02:36:34.040 --> 02:36:37.440
and that's, you know,
it's just a computational trick,
1715
02:36:37.540 --> 02:36:42.820
but by making this effective mass
very small, it can be negligible.
1716
02:36:42.920 --> 02:36:46.680
And also the blue curve that would be this
1717
02:36:46.780 --> 02:36:50.760
expression becomes closer
and closer to the red curve
1718
02:36:50.860 --> 02:36:53.300
that would be the right expression.
1719
02:36:53.400 --> 02:36:57.180
The only cost that you have to pay
if you make a mu too small
1720
02:36:57.280 --> 02:37:03.780
is that the coefficient of the plane-waves
change so fast that you're time step
1721
02:37:03.880 --> 02:37:07.940
for the Car-Parrinello integration
needs to become very small.
1722
02:37:08.040 --> 02:37:09.960
And that's a bit of a technicality.
1723
02:37:10.060 --> 02:37:17.460
But if you're able to do this calculation well,
things can become very efficient.
1724
02:37:17.560 --> 02:37:23.460
You have a very flat physical constant of motion,
1725
02:37:23.560 --> 02:37:29.860
because the total constant
of motion is very well preserved.
1726
02:37:29.960 --> 02:37:32.485
And this, you know,
1727
02:37:32.486 --> 02:37:38.671
unphysical quantity is very small
and very flat. These kinetic energies,
1728
02:37:38.672 --> 02:37:41.680
you know, we have hundreds
of thousands of coefficients
1729
02:37:41.780 --> 02:37:46.600
of plane-waves, but their total kinetic
energy can be a tiny fraction,
1730
02:37:46.601 --> 02:37:52.020
one tenth, one twentieth of the true
physical kinetic energy of the ions.
1731
02:37:52.120 --> 02:37:58.600
So the coefficient of the plane-waves follow,
but they don't give back energies to the ions.
1732
02:37:59.760 --> 02:38:04.000
So in a nutshell, you can do
1733
02:38:04.160 --> 02:38:06.820
Born-Oppenheimer
1734
02:38:06.920 --> 02:38:16.600
molecular dynamic simulation, where you move
the ions. You converge the electrons back to...
1735
02:38:19.200 --> 02:38:26.700
to the ground state with, say, iterative procedure
or conjugate gradient procedures,
1736
02:38:26.800 --> 02:38:31.540
but you're always, sort of, you know,
you have always a little bit of a systematic error,
1737
02:38:31.640 --> 02:38:36.060
so your constant of motion will not be very constant.
1738
02:38:36.160 --> 02:38:38.260
Instead, with Car-Parrinello techniques,
1739
02:38:38.360 --> 02:38:47.220
your constant of motion is exact but you need
to make sure that you have integrated them.
1740
02:38:47.320 --> 02:38:51.900
Things can go well in a semi-conductor or insulator.
1741
02:38:52.000 --> 02:38:55.680
The issue, and then I'll conclude, is what happens
1742
02:38:55.780 --> 02:39:01.360
if you have a metal. And if you happen,
what if ... what happens if you have a metal
1743
02:39:01.460 --> 02:39:07.800
in a variational algorithm is that you
start to have three degrees of freedom,
1744
02:39:07.900 --> 02:39:15.860
that is, you need to evolve your wavefunction
but you need also to evolve your occupations.
1745
02:39:15.960 --> 02:39:25.686
And you need that to make sure that your
occupations are calculated on the eigenvalues of
1746
02:39:28.960 --> 02:39:36.243
the Hamiltonian, so, you know, what is
the representation that diagonalises
1747
02:39:36.244 --> 02:39:41.020
the Hamiltonian, so that you can get
the actual eigenvalues of the system.
1748
02:39:41.120 --> 02:39:44.300
So in a metal, as opposed to a semiconductor
1749
02:39:44.301 --> 02:39:49.860
or an insulator, you don't have
to evolve only the wavefunctions.
1750
02:39:49.960 --> 02:39:55.640
In a semiconductor or insulator, the total energy
is invariant for unitary transformation.
1751
02:39:55.740 --> 02:39:57.820
So you don't care about that.
1752
02:39:57.920 --> 02:40:02.080
But instead in a metal,
you need to make sure that you evolve
1753
02:40:02.081 --> 02:40:08.020
also the occupations that can change
from zero to one, or they can be fractionally.
1754
02:40:08.120 --> 02:40:11.300
And you always need to know what is
1755
02:40:11.400 --> 02:40:17.360
the eigenvalue of the Hamiltonian that is unique
to being the representation where you can
1756
02:40:17.460 --> 02:40:22.340
calculate eigenvalues because
the problem is not invariant anymore.
1757
02:40:22.440 --> 02:40:25.829
So you have, if you want,
two additional degrees of freedom,
1758
02:40:25.857 --> 02:40:29.780
occupation numbers and unitary rotations.
1759
02:40:29.880 --> 02:40:35.480
And these occupation numbers can be
very effectively combined with the unitary
1760
02:40:35.580 --> 02:40:43.340
rotation to, say, in a metal, the additional
degree of freedom is actually a matrix
1761
02:40:43.440 --> 02:40:45.500
that is the representation.
1762
02:40:45.600 --> 02:40:49.200
Say, if this was a standard finite temperature problem,
1763
02:40:49.201 --> 02:40:57.813
it would be the representation of the Fermi-Dirac,
operator on the basis of wavefunction psi i.
1764
02:40:57.814 --> 02:41:02.260
And, in general, this matrix is non-diagonal.
1765
02:41:02.360 --> 02:41:09.680
And so to deal with this non-diagonal matrix,
you need to write the charge density in a covariant
1766
02:41:09.780 --> 02:41:14.400
form like this,
and you need to write the total energy
1767
02:41:14.500 --> 02:41:19.460
and the external potential
in a covariant form, as a double form.
1768
02:41:19.560 --> 02:41:21.340
And of course, you could always
1769
02:41:21.440 --> 02:41:30.360
choose a special representation in which
this f is diagonal, but then you would need also
1770
02:41:30.460 --> 02:41:37.394
to evolve the unitary matrices
that connect the psi H psi
1771
02:41:37.594 --> 02:41:41.220
with its eigenvalues, epsilon i.
1772
02:41:41.320 --> 02:41:43.100
So there is no winning.
1773
02:41:43.200 --> 02:41:46.620
I mean, you have all this complexity.
1774
02:41:46.720 --> 02:41:50.040
And the nice thing about this covariant
1775
02:41:50.140 --> 02:41:55.740
representation is that it allows you to define
this functional G
1776
02:41:55.840 --> 02:42:01.160
that is defined as the minimum
of the free-energy functional with respect
1777
02:42:01.260 --> 02:42:05.820
to this matrix representation of the occupations.
1778
02:42:05.920 --> 02:42:13.140
And you see this functional G is a functional
where occupations have disappeared.
1779
02:42:13.240 --> 02:42:16.240
Not only that but is a functional that is
1780
02:42:16.340 --> 02:42:20.520
invariant for a unitary
transformation of the orbital.
1781
02:42:20.620 --> 02:42:27.320
So if I take the psi and I send them
into psi U, U is a unitary matrix,
1782
02:42:27.420 --> 02:42:29.940
I have done a unitary transformation.
1783
02:42:30.040 --> 02:42:32.280
What happens is that the values of this functional
1784
02:42:32.380 --> 02:42:40.100
doesn't change because if I do a transformation
of the psi into psi U,
1785
02:42:40.200 --> 02:42:48.620
what I can do is a transformation of the U into
U dagger f U and the value has not changed.
1786
02:42:48.720 --> 02:42:56.300
So if I rotated psi, the f counter rotates
and the minimum doesn't change.
1787
02:42:56.400 --> 02:43:03.420
So this trick gives me a functional G that looks
like the functional of a semiconductor or insulator.
1788
02:43:03.520 --> 02:43:09.020
It doesn't have explicit occupation numbers
and it's invariant for unitary transformation,
1789
02:43:09.120 --> 02:43:12.400
so I can do conjugate gradient iterative
1790
02:43:12.500 --> 02:43:20.614
techniques and so on in this and I can get to
very stable molecular dynamics again.
1791
02:43:20.757 --> 02:43:24.814
This is called... we called it ensemble DFT.
1792
02:43:24.900 --> 02:43:27.460
It's coded in VASP.
1793
02:43:27.560 --> 02:43:34.380
This was actually done by Freysoldt
and collaborators. It's coded in CASTEP.
1794
02:43:34.480 --> 02:43:36.560
It's not coded in Quantum Espresso.
1795
02:43:36.660 --> 02:43:38.300
That's somehow very sad
1796
02:43:38.400 --> 02:43:46.360
but I think in a few months we'll be able to make
a release of this together with the GPU library
1797
02:43:46.460 --> 02:43:53.800
serials that Anton Kozhevnikov and
collaborator are developing at CSCS.
1798
02:43:54.200 --> 02:44:02.540
So, I think this, and that's the second lesson,
I've been only 45 minutes longer.
1799
02:44:02.640 --> 02:44:04.640
I think there were a lot of technicalities
1800
02:44:04.740 --> 02:44:08.680
so I figure out that only people that
want to do actual calculations
1801
02:44:08.780 --> 02:44:12.620
were going to follow this
and be interested in this.
1802
02:44:12.720 --> 02:44:15.620
Tomorrow we'll go back to 2 hrs
1803
02:44:15.720 --> 02:44:21.940
and a much more conceptual discussion. This will
be recorded, will be put on the Materials Cloud.
1804
02:44:22.040 --> 02:44:25.500
If you want to do the calculation,
go to the Materials Cloud.
1805
02:44:25.600 --> 02:44:29.420
In the fireside section, there is a link to GitHub.
1806
02:44:29.520 --> 02:44:34.629
There is a complete handout that tells you how
to download and install Quantum Espresso
1807
02:44:34.771 --> 02:44:42.020
or how to download and install the Quantum Mobile,
the virtual machine with Quantum Espresso
1808
02:44:42.120 --> 02:44:50.685
preloaded into it and one by one there are
homeworks that tells you and, you know,
1809
02:44:50.686 --> 02:44:54.216
teach you how to set up
the cutoff for your calculation,
1810
02:44:54.416 --> 02:44:56.820
how to check on the k-point sampling,
1811
02:44:56.920 --> 02:45:02.820
how to do a simple calculation,
an equation of state for sodium chloride.
1812
02:45:02.920 --> 02:45:09.020
If you are lost, the Quantum Espresso input
generator always gives you a good start
1813
02:45:09.120 --> 02:45:14.760
to, you know, have a good chance
for an input file for Quantum Espresso.
1814
02:45:14.860 --> 02:45:20.260
The SSSP library gives you tested
pseudopotential and tested cutoff.
1815
02:45:20.360 --> 02:45:27.300
And there is always a mailing list where you can
ask questions for the Quantum Espresso users.
1816
02:45:27.400 --> 02:45:32.760
And if you're interested also in more on this,
we have put on the lower section, also,
1817
02:45:32.860 --> 02:45:35.300
you can enrol in a mailing list for
1818
02:45:35.400 --> 02:45:42.200
more fireside chats. With this I'll conclude.
The panellists have collected all your questions.
1819
02:45:42.300 --> 02:45:45.720
We'll put them in a PDF in the next few days.
1820
02:45:45.820 --> 02:45:50.220
We'll find them and we'll post them next week.
1821
02:45:50.320 --> 02:45:54.020
Other than this, I look forward
to seeing you tomorrow.
1822
02:45:54.120 --> 02:45:57.660
And this we'll also post on the lower section.
1823
02:45:57.760 --> 02:46:01.380
Thanks a lot and thanks especially exactly
1824
02:46:01.480 --> 02:46:06.060
to Davide, Francisco, Luca and Norma
that have prepared the handout
1825
02:46:06.160 --> 02:46:12.629
and then been sitting here all of these
2 hours and 45 minutes listening to this,
1826
02:46:12.929 --> 02:46:16.640
and, you know, see you indirectly tomorrow.