1
00:00:00,548 --> 00:00:01,340
NICOLA MARZARI: OK.

2
00:00:01,340 --> 00:00:05,330
Welcome, everyone,
lecture 12, second part

3
00:00:05,330 --> 00:00:09,350
of molecular dynamics, also
after the spring break.

4
00:00:09,350 --> 00:00:13,110
So I decided to spend two or
three slides in reminding you

5
00:00:13,110 --> 00:00:16,430
what was going on
from last time.

6
00:00:16,430 --> 00:00:18,000
And I've actually
been intrigued.

7
00:00:18,000 --> 00:00:20,030
I'm starting with
Heraclitus because this

8
00:00:20,030 --> 00:00:22,220
is a beautiful citation.

9
00:00:22,220 --> 00:00:24,470
So this is, you know,
fifth century BC.

10
00:00:24,470 --> 00:00:28,940
But there is the intuition
that a macroscopic property,

11
00:00:28,940 --> 00:00:33,320
like the river, actually comes
from a microscopic property.

12
00:00:33,320 --> 00:00:35,570
That is the water molecules.

13
00:00:35,570 --> 00:00:38,540
Though he wouldn't know
what a water molecule

14
00:00:38,540 --> 00:00:41,570
is, but comes from what the
water molecules are doing.

15
00:00:41,570 --> 00:00:44,550
So Heraclitus was saying, you
cannot step twice in the same

16
00:00:44,550 --> 00:00:45,050
river.

17
00:00:45,050 --> 00:00:47,540
Because the molecular
conformation

18
00:00:47,540 --> 00:00:51,260
is always different, but the
river is always the same.

19
00:00:51,260 --> 00:00:55,260
Anyhow, reminder of what
we have seen last time--

20
00:00:55,260 --> 00:00:57,890
first of all, of course,
the standard thing

21
00:00:57,890 --> 00:01:01,340
is that we have looked at
Newton's equation of motion.

22
00:01:01,340 --> 00:01:04,220
Really, molecular
dynamic simulations

23
00:01:04,220 --> 00:01:07,730
involve the integration
of these couple

24
00:01:07,730 --> 00:01:10,620
set of second-order
differential equation.

25
00:01:10,620 --> 00:01:12,650
So if we have n
particles, there are

26
00:01:12,650 --> 00:01:15,930
going to be n
trajectory defined by,

27
00:01:15,930 --> 00:01:20,210
for each particle, how the
vector position evolves

28
00:01:20,210 --> 00:01:21,230
with time.

29
00:01:21,230 --> 00:01:25,400
And the force acting
on each particle

30
00:01:25,400 --> 00:01:27,440
can be obtained
from the derivative

31
00:01:27,440 --> 00:01:28,910
of the potential field.

32
00:01:28,910 --> 00:01:31,860
And all the particles
interact with each other.

33
00:01:31,860 --> 00:01:34,010
So this is still a
many-bodied problem,

34
00:01:34,010 --> 00:01:36,170
although a much
easier one to solve

35
00:01:36,170 --> 00:01:39,230
than the electronic problem.

36
00:01:39,230 --> 00:01:42,320
And we sort of remind
you how we integrated

37
00:01:42,320 --> 00:01:44,250
this equation of motion.

38
00:01:44,250 --> 00:01:46,970
But sort of one of the
fundamental concepts that we

39
00:01:46,970 --> 00:01:49,430
have seen is how
we can use actually

40
00:01:49,430 --> 00:01:52,580
the integration of the molecular
dynamics equation of motion

41
00:01:52,580 --> 00:01:55,310
to obtain thermodynamical
quantities.

42
00:01:55,310 --> 00:01:59,180
And you see what I'm
writing here on the left

43
00:01:59,180 --> 00:02:02,630
is really the Boltzmann
idea, the sort

44
00:02:02,630 --> 00:02:05,930
of basic idea of
statistical mechanics.

45
00:02:05,930 --> 00:02:10,310
That is, if you have a certain
variable that I'm generically

46
00:02:10,310 --> 00:02:15,410
calling here A, what we
know is that the expectation

47
00:02:15,410 --> 00:02:18,500
value for that
thermodynamical variable

48
00:02:18,500 --> 00:02:21,690
is given by the integration
over all the phase

49
00:02:21,690 --> 00:02:23,960
space, integration
over all the position,

50
00:02:23,960 --> 00:02:29,420
integration over all the momenta
of the value of that variable A

51
00:02:29,420 --> 00:02:31,670
for every position
and every momentum

52
00:02:31,670 --> 00:02:35,580
weighed with the appropriate
probability distribution.

53
00:02:35,580 --> 00:02:37,130
So this would be,
say, for what we

54
00:02:37,130 --> 00:02:41,990
call a canonical distribution,
a distribution in which

55
00:02:41,990 --> 00:02:46,190
each microstate, its single
position and momentum,

56
00:02:46,190 --> 00:02:50,120
is weighed with the
exponential of the energy

57
00:02:50,120 --> 00:02:52,850
with a coefficient
beta that is just 1

58
00:02:52,850 --> 00:02:56,990
over the Boltzmann constant
times the temperature.

59
00:02:56,990 --> 00:03:02,030
And this fundamental sort of
statistical mechanics approach

60
00:03:02,030 --> 00:03:05,360
can actually be translated into
a molecular dynamics approach.

61
00:03:05,360 --> 00:03:09,320
That is we can actually
obtain the same expectation

62
00:03:09,320 --> 00:03:14,930
value for the operator A without
really doing these integrals.

63
00:03:14,930 --> 00:03:18,200
These integrals are obviously
impossible to perform

64
00:03:18,200 --> 00:03:20,780
explicitly as soon as
you have many particles

65
00:03:20,780 --> 00:03:25,160
because the dimension of
your phase space explodes.

66
00:03:25,160 --> 00:03:26,810
And we have seen it
over and over again

67
00:03:26,810 --> 00:03:30,270
that, as the
dimension increases,

68
00:03:30,270 --> 00:03:33,470
it becomes increasingly
difficult or computationally

69
00:03:33,470 --> 00:03:35,720
impossible to perform
these integrals.

70
00:03:35,720 --> 00:03:39,840
But what we can actually
do is evolve our system.

71
00:03:39,840 --> 00:03:41,690
But, remember,
our system is just

72
00:03:41,690 --> 00:03:44,660
a point in phase space, the
collection of all position

73
00:03:44,660 --> 00:03:45,710
and momenta.

74
00:03:45,710 --> 00:03:49,550
We evolve it with the molecular
dynamics equation of motion.

75
00:03:49,550 --> 00:03:51,870
And in particular,
say for this example,

76
00:03:51,870 --> 00:03:55,750
we need to evolve it in what
we call the canonical ensemble.

77
00:03:55,750 --> 00:03:58,250
That is we need to evolve
under the condition

78
00:03:58,250 --> 00:04:01,130
that the temperature of
the system is constant.

79
00:04:01,130 --> 00:04:04,040
But in this trajectory,
in this evolution,

80
00:04:04,040 --> 00:04:08,840
we can average for each
instant of this trajectory what

81
00:04:08,840 --> 00:04:12,380
is the instantaneous value
of your observable A.

82
00:04:12,380 --> 00:04:15,050
And the average
over the trajectory

83
00:04:15,050 --> 00:04:17,279
gives us an expectation value.

84
00:04:17,279 --> 00:04:20,149
And these two quantities
are equivalent.

85
00:04:20,149 --> 00:04:22,230
They are obtained in
two different ways.

86
00:04:22,230 --> 00:04:24,470
Here, as an average
over the trajectory,

87
00:04:24,470 --> 00:04:26,990
that is much, much simpler
to do computationally.

88
00:04:26,990 --> 00:04:30,260
And here is an average
of all phase space.

89
00:04:30,260 --> 00:04:33,950
And there is a subtlety about
this equivalence that is called

90
00:04:33,950 --> 00:04:36,110
the hypothesis of ergodicity.

91
00:04:36,110 --> 00:04:39,470
That is, in order for this
equivalence to be valid,

92
00:04:39,470 --> 00:04:43,670
we really need to think that
our system spans the phase

93
00:04:43,670 --> 00:04:46,820
space with the appropriate
probability distribution

94
00:04:46,820 --> 00:04:50,150
without being confined
in corner of phase space.

95
00:04:50,150 --> 00:04:53,000
And I'll give you an
example later on of what

96
00:04:53,000 --> 00:04:54,890
means being confined.

97
00:04:54,890 --> 00:04:58,010
So in a way, we sort of
work under the hypothesis

98
00:04:58,010 --> 00:05:02,300
that really this trajectory is
spanning an appropriate sample

99
00:05:02,300 --> 00:05:03,693
of trajectories.

100
00:05:03,693 --> 00:05:06,110
And so when you do a molecular
dynamic simulation-- again,

101
00:05:06,110 --> 00:05:08,840
reminder of last class--

102
00:05:08,840 --> 00:05:13,400
you really perform an ideal
computational experiment

103
00:05:13,400 --> 00:05:15,920
in which you start
deciding what are

104
00:05:15,920 --> 00:05:19,700
the position and the velocity
of all your particles.

105
00:05:19,700 --> 00:05:23,330
And then this is sort of the
computationally important step.

106
00:05:23,330 --> 00:05:26,400
You integrate the
equation of motion.

107
00:05:26,400 --> 00:05:28,970
And this is the sort
of step in which

108
00:05:28,970 --> 00:05:33,590
you need to pay particular care
and attention on the algorithm

109
00:05:33,590 --> 00:05:36,110
you use to integrate
your equation of motion

110
00:05:36,110 --> 00:05:39,980
and, also, on the accuracy
that is needed to integrate

111
00:05:39,980 --> 00:05:41,870
this equation of motion.

112
00:05:41,870 --> 00:05:50,060
And once you have put all
your particles in motion

113
00:05:50,060 --> 00:05:53,240
and you are sort of integrating
them and evolving them

114
00:05:53,240 --> 00:05:56,160
according to the
appropriate equations,

115
00:05:56,160 --> 00:06:00,020
well, you just need to make
sure that you wait long enough.

116
00:06:00,020 --> 00:06:02,270
That is that you
equilibrate your system,

117
00:06:02,270 --> 00:06:04,160
so that you
basically lose memory

118
00:06:04,160 --> 00:06:05,630
of your initial condition.

119
00:06:05,630 --> 00:06:09,620
You see, your system
starts from a choice

120
00:06:09,620 --> 00:06:13,010
that has been made
by you and by hand.

121
00:06:13,010 --> 00:06:16,220
And in data, it's
never representative

122
00:06:16,220 --> 00:06:19,100
of what could really be a
microscopic configuration

123
00:06:19,100 --> 00:06:19,760
of your system.

124
00:06:19,760 --> 00:06:21,980
You're putting molecules
certain position,

125
00:06:21,980 --> 00:06:23,750
and you're giving
them velocities.

126
00:06:23,750 --> 00:06:27,510
But that's just your assumption.

127
00:06:27,510 --> 00:06:31,290
And so you need to lose memory
of this initial condition.

128
00:06:31,290 --> 00:06:35,930
And that's accomplished
by evolving the system

129
00:06:35,930 --> 00:06:37,190
long enough.

130
00:06:37,190 --> 00:06:41,480
And how long is the time
needed to reach equilibrium

131
00:06:41,480 --> 00:06:43,640
is something that
depends on your system,

132
00:06:43,640 --> 00:06:45,350
depends on the
temperature at which you

133
00:06:45,350 --> 00:06:46,850
are performing a simulation.

134
00:06:46,850 --> 00:06:49,280
But it's something
that you can check

135
00:06:49,280 --> 00:06:51,500
in the course of
your simulation.

136
00:06:51,500 --> 00:06:53,420
And again, I'll
show you an example

137
00:06:53,420 --> 00:06:57,620
on how you check for
appropriate equilibration.

138
00:06:57,620 --> 00:06:59,870
And once you have
reached equilibrium,

139
00:06:59,870 --> 00:07:03,770
well, then it becomes
straightforward to just

140
00:07:03,770 --> 00:07:06,350
calculate the results
of your experiments.

141
00:07:06,350 --> 00:07:10,700
That is calculate, say, what are
the thermodynamical averages.

142
00:07:10,700 --> 00:07:12,920
And so you start
accumulating what are

143
00:07:12,920 --> 00:07:14,310
your quantities of interest.

144
00:07:16,850 --> 00:07:23,440
And again, sort of reminder of
last class, the initialization

145
00:07:23,440 --> 00:07:27,040
step involves choosing
the initial position

146
00:07:27,040 --> 00:07:28,720
and the initial velocities.

147
00:07:28,720 --> 00:07:30,670
And we need position
and velocity

148
00:07:30,670 --> 00:07:34,480
because really we are dealing
with second-order differential

149
00:07:34,480 --> 00:07:36,250
equations.

150
00:07:36,250 --> 00:07:38,590
It's more critical to
choose appropriately

151
00:07:38,590 --> 00:07:41,500
the initial position
and more difficult.

152
00:07:41,500 --> 00:07:43,850
If you are studying a
system that is crystalline,

153
00:07:43,850 --> 00:07:45,310
obviously it's very simple.

154
00:07:45,310 --> 00:07:48,970
You start just with sort of a
periodic arrangement of atoms.

155
00:07:48,970 --> 00:07:52,360
But say if you are starting
something more complex,

156
00:07:52,360 --> 00:07:56,620
like, say, liquid water,
it becomes more critical.

157
00:07:56,620 --> 00:07:59,740
Because your liquid,
as a structure,

158
00:07:59,740 --> 00:08:02,215
that is not immediately
evident to you

159
00:08:02,215 --> 00:08:04,570
when you try to put
molecules around.

160
00:08:04,570 --> 00:08:07,510
Anyhow, there are sort of
obvious things to avoid.

161
00:08:07,510 --> 00:08:12,250
That is you want to avoid that
atoms overlap too closely.

162
00:08:12,250 --> 00:08:14,860
Because if you put two
atoms too close together,

163
00:08:14,860 --> 00:08:16,600
they will have an
initial repulsive

164
00:08:16,600 --> 00:08:18,130
force that is very strong.

165
00:08:18,130 --> 00:08:22,940
And they'll just zoom away
very quickly from each other.

166
00:08:22,940 --> 00:08:26,020
It's much more critical to
equilibrate with respect

167
00:08:26,020 --> 00:08:29,950
to position than with
respect to velocities.

168
00:08:29,950 --> 00:08:34,194
Velocities, or if you want
the phrase space part,

169
00:08:34,194 --> 00:08:38,070
the momentum part
of your phase space,

170
00:08:38,070 --> 00:08:39,929
is much easier to integrate.

171
00:08:39,929 --> 00:08:43,500
Basically, because in any
collision between atoms

172
00:08:43,500 --> 00:08:46,480
or between molecules, the
velocity, if you want,

173
00:08:46,480 --> 00:08:48,190
get reversed.

174
00:08:48,190 --> 00:08:52,800
And so in the momentum, in the
velocity part of phase space,

175
00:08:52,800 --> 00:08:57,060
your representative points
can really jump around.

176
00:08:57,060 --> 00:09:00,360
Say, if you had an ideal
gas of hard spheres

177
00:09:00,360 --> 00:09:03,870
that only collide when
they touch each other--

178
00:09:03,870 --> 00:09:05,760
and in an elastic
collision, they

179
00:09:05,760 --> 00:09:08,100
just reverse their velocity--

180
00:09:08,100 --> 00:09:11,160
well, in that case, your
representative point

181
00:09:11,160 --> 00:09:14,520
is really jumping
around discontinuously.

182
00:09:14,520 --> 00:09:19,020
So in velocity space, you
really move around very easily.

183
00:09:19,020 --> 00:09:21,840
And you span all
the phase space.

184
00:09:21,840 --> 00:09:25,410
It's really in the position
part of the phase space

185
00:09:25,410 --> 00:09:28,500
that you need to be careful
and where your system

186
00:09:28,500 --> 00:09:30,600
can become entangled.

187
00:09:36,170 --> 00:09:39,500
Once you have set up
your initial position

188
00:09:39,500 --> 00:09:41,150
and your initial
velocity, you have

189
00:09:41,150 --> 00:09:44,720
everything you need to integrate
your equation of motion.

190
00:09:44,720 --> 00:09:47,240
And to integrate
equation of motion,

191
00:09:47,240 --> 00:09:49,310
there are a variety
of algorithms.

192
00:09:49,310 --> 00:09:51,770
And I have mentioned,
and you'll see again,

193
00:09:51,770 --> 00:09:55,790
the standard or simple
Verlet algorithm.

194
00:09:55,790 --> 00:09:59,330
That is actually a very robust
and very accurate algorithm.

195
00:09:59,330 --> 00:10:02,220
So it's what you'll
see in the next slide,

196
00:10:02,220 --> 00:10:05,630
but there are sort of different
varieties of Verlet algorithm.

197
00:10:05,630 --> 00:10:09,800
Sometimes you see them called
leapfrog Verlet or velocity

198
00:10:09,800 --> 00:10:13,130
Verlet depending on the
specific integration scheme.

199
00:10:13,130 --> 00:10:16,310
And then there is another class
that I'll also mention briefly

200
00:10:16,310 --> 00:10:18,320
just for reference
that are called

201
00:10:18,320 --> 00:10:22,200
the Gear predictor-corrector
integration algorithm.

202
00:10:22,200 --> 00:10:25,970
So this, more or less,
relies on all the zoology

203
00:10:25,970 --> 00:10:28,730
of integration algorithm.

204
00:10:28,730 --> 00:10:33,380
And obviously, one of
the important qualities

205
00:10:33,380 --> 00:10:37,580
of your integration algorithm
is that it must be robust.

206
00:10:37,580 --> 00:10:41,960
That is it shouldn't be
too sensitive on the choice

207
00:10:41,960 --> 00:10:43,550
of discretization.

208
00:10:43,550 --> 00:10:45,950
Remember, we have sort
of differential equation

209
00:10:45,950 --> 00:10:48,860
that are basically perfect
analytical equation.

210
00:10:48,860 --> 00:10:51,020
But in order to integrate
them on a computer,

211
00:10:51,020 --> 00:10:52,880
we need to discretize them.

212
00:10:52,880 --> 00:10:56,930
And again, sort of one of the
key features of your simulation

213
00:10:56,930 --> 00:11:01,280
is figuring out what is
the discretization interval

214
00:11:01,280 --> 00:11:03,080
or what we call the time step.

215
00:11:03,080 --> 00:11:08,300
That is how much
time you can spend

216
00:11:08,300 --> 00:11:13,040
before you calculate another
point in your trajectory.

217
00:11:13,040 --> 00:11:15,050
And there are sort of
a lot of subtleties

218
00:11:15,050 --> 00:11:22,130
that we won't go into on what
really makes a good integrator

219
00:11:22,130 --> 00:11:24,570
for the equation of motion.

220
00:11:24,570 --> 00:11:27,120
In particular, one of the
properties that we would want

221
00:11:27,120 --> 00:11:29,090
is that it's time reversible.

222
00:11:29,090 --> 00:11:31,010
And that's what the
simple Verlet is.

223
00:11:31,010 --> 00:11:34,880
That is, if you invert
all your velocities,

224
00:11:34,880 --> 00:11:39,200
your trajectory is recovered
in the opposite direction.

225
00:11:39,200 --> 00:11:42,620
And not all algorithms
actually are time reversible.

226
00:11:42,620 --> 00:11:45,140
And then there is something
a little bit more subtle.

227
00:11:45,140 --> 00:11:47,810
That is a sort of
criterion of robustness

228
00:11:47,810 --> 00:11:50,210
for an integrator
is that it preserves

229
00:11:50,210 --> 00:11:53,690
the volume of your
representative system in phase

230
00:11:53,690 --> 00:11:58,670
space, but we really want
to go into these details.

231
00:11:58,670 --> 00:12:02,000
There are sort of, you
know, several classes

232
00:12:02,000 --> 00:12:04,610
of molecular dynamic
simulations depending

233
00:12:04,610 --> 00:12:07,220
on what is the thermodynamical
ensemble that you

234
00:12:07,220 --> 00:12:08,540
want to represent.

235
00:12:08,540 --> 00:12:10,610
That is, if you
straightforwardly

236
00:12:10,610 --> 00:12:13,400
integrate Newton
equation of motion

237
00:12:13,400 --> 00:12:16,100
we have seen in the
last class, you actually

238
00:12:16,100 --> 00:12:18,980
conserve the total
energy of the system.

239
00:12:18,980 --> 00:12:21,890
That is you conserve the
sum of the kinetic energy

240
00:12:21,890 --> 00:12:23,990
plus the sum of the
potential energy.

241
00:12:23,990 --> 00:12:26,180
And as you know from
statistical mechanics,

242
00:12:26,180 --> 00:12:30,860
that is what is called the
microcanonical ensemble.

243
00:12:30,860 --> 00:12:35,180
Sometimes it might be
more worthwhile to perform

244
00:12:35,180 --> 00:12:39,060
your simulation at a
constant temperature.

245
00:12:39,060 --> 00:12:43,460
And so that is what you would
call the canonical ensemble

246
00:12:43,460 --> 00:12:45,260
in statistical mechanics.

247
00:12:45,260 --> 00:12:49,730
And so you need to find a way
to impose a constant temperature

248
00:12:49,730 --> 00:12:50,750
in your system.

249
00:12:50,750 --> 00:12:54,000
Your system, in principle,
would conserve the energy.

250
00:12:54,000 --> 00:12:57,450
And so it would have a
fluctuating temperature.

251
00:12:57,450 --> 00:12:59,510
And again, from
statistical mechanics,

252
00:12:59,510 --> 00:13:03,200
the larger your system is,
that is the closer you get

253
00:13:03,200 --> 00:13:05,600
to the thermodynamically
limit, to the limit

254
00:13:05,600 --> 00:13:07,760
of an infinitely large system.

255
00:13:07,760 --> 00:13:12,650
It really doesn't matter because
an infinitely large system that

256
00:13:12,650 --> 00:13:16,460
is simulated with, say, a
microcanonical algorithm

257
00:13:16,460 --> 00:13:21,320
will also have 0 fluctuations
on the average temperature.

258
00:13:21,320 --> 00:13:26,990
And so you would obtain the
same results in a microcanonical

259
00:13:26,990 --> 00:13:28,490
or in a canonical simulation.

260
00:13:28,490 --> 00:13:33,290
But as always, we deal with
finite systems maybe replicated

261
00:13:33,290 --> 00:13:35,340
with periodic
boundary condition.

262
00:13:35,340 --> 00:13:38,690
And so it might be more
convenient to choose

263
00:13:38,690 --> 00:13:41,270
a canonical simulation,
so that we can fix

264
00:13:41,270 --> 00:13:42,650
the temperature of our system.

265
00:13:42,650 --> 00:13:44,870
Suppose that, again,
you want to study water.

266
00:13:44,870 --> 00:13:46,670
You want to make sure
that you are studying

267
00:13:46,670 --> 00:13:49,670
your system at a
temperature that is

268
00:13:49,670 --> 00:13:52,190
consistent with liquid water.

269
00:13:52,190 --> 00:13:54,410
Suppose that you
were studying water

270
00:13:54,410 --> 00:13:56,810
with a microcanonical
simulation.

271
00:13:56,810 --> 00:14:01,550
Again, what would you do is
choose your initial positions,

272
00:14:01,550 --> 00:14:05,990
choose your initial velocities,
and let your system evolve.

273
00:14:05,990 --> 00:14:09,290
But it's only after you
have let your system evolve

274
00:14:09,290 --> 00:14:13,640
for a long time, that is after
your system has equilibrated,

275
00:14:13,640 --> 00:14:17,570
that you can actually
measure, start averaging, some

276
00:14:17,570 --> 00:14:19,190
of the macroscopic properties.

277
00:14:19,190 --> 00:14:21,450
And so in the
microcanonical ensemble,

278
00:14:21,450 --> 00:14:23,960
you could start, after
a while, to average

279
00:14:23,960 --> 00:14:26,390
as a macroscopic
property the temperature.

280
00:14:26,390 --> 00:14:29,090
The temperature is nothing else
than the average kinetic energy

281
00:14:29,090 --> 00:14:29,980
of the ions.

282
00:14:29,980 --> 00:14:31,730
Well, suppose that you
are studying water.

283
00:14:31,730 --> 00:14:33,320
You have set up your
initial position,

284
00:14:33,320 --> 00:14:34,580
your initial velocity.

285
00:14:34,580 --> 00:14:35,990
You let this evolve.

286
00:14:35,990 --> 00:14:38,780
And after a while, you discover
that the average temperature

287
00:14:38,780 --> 00:14:40,550
is 150 kelvin.

288
00:14:40,550 --> 00:14:42,140
And your system is frozen.

289
00:14:42,140 --> 00:14:42,740
OK.

290
00:14:42,740 --> 00:14:46,520
Well, what you have to do is
then restart your simulation.

291
00:14:46,520 --> 00:14:49,700
It's not going anywhere if you
want to study liquid water.

292
00:14:49,700 --> 00:14:51,230
And how could you change this?

293
00:14:51,230 --> 00:14:55,550
Well, you could try using the
same initial configuration

294
00:14:55,550 --> 00:14:59,690
that you had used before,
but much larger velocities,

295
00:14:59,690 --> 00:15:04,100
hoping that once the system
evolves deterministically

296
00:15:04,100 --> 00:15:07,100
from that initial point,
your average temperature

297
00:15:07,100 --> 00:15:10,850
is going to be a temperature
in which your water is liquid.

298
00:15:10,850 --> 00:15:13,520
But maybe you do it again, and
your velocity is too large.

299
00:15:13,520 --> 00:15:16,100
And now, the temperature of
your system is 500 kelvin.

300
00:15:16,100 --> 00:15:16,700
OK.

301
00:15:16,700 --> 00:15:18,980
So you see the
problem, in this case,

302
00:15:18,980 --> 00:15:21,900
in starting with a
microcanonical approach.

303
00:15:21,900 --> 00:15:24,830
So it would be much
more convenient to have

304
00:15:24,830 --> 00:15:27,590
an integration scheme
that actually allows

305
00:15:27,590 --> 00:15:30,530
you to fix the temperature
instead of allowing

306
00:15:30,530 --> 00:15:32,510
you to conserve the energy.

307
00:15:32,510 --> 00:15:35,660
And there are several ways
in which something like this

308
00:15:35,660 --> 00:15:36,720
can be done.

309
00:15:36,720 --> 00:15:39,380
And again, you'll see examples.

310
00:15:39,380 --> 00:15:41,810
But generally speaking,
we say that we

311
00:15:41,810 --> 00:15:45,320
couple our dynamical
system with a thermostat.

312
00:15:45,320 --> 00:15:49,340
That is nothing less than, if
you want, a computing device

313
00:15:49,340 --> 00:15:53,690
that allows us to control,
in this case, the temperature

314
00:15:53,690 --> 00:15:55,590
of our simulation.

315
00:15:55,590 --> 00:15:57,080
And if instead of
the temperature

316
00:15:57,080 --> 00:15:59,720
we wanted to control
the pressure,

317
00:15:59,720 --> 00:16:02,990
we would couple our
system with a [INAUDIBLE]

318
00:16:02,990 --> 00:16:04,590
and so on and so forth.

319
00:16:04,590 --> 00:16:07,400
But if we focus for a moment
on the thermostat that's

320
00:16:07,400 --> 00:16:10,580
the simplest, well, there
are sort of several ways

321
00:16:10,580 --> 00:16:14,280
in which we can control the
temperature of our system.

322
00:16:14,280 --> 00:16:21,200
And I've mentioned here sort
of three classes of devices

323
00:16:21,200 --> 00:16:25,610
that we use to couple our
system with a thermostat.

324
00:16:25,610 --> 00:16:30,290
And one approach could
be the Langevin dynamics

325
00:16:30,290 --> 00:16:32,360
or the stochastic approach.

326
00:16:32,360 --> 00:16:36,740
That is you try to sort of
simulate a constant temperature

327
00:16:36,740 --> 00:16:39,800
by letting your system evolve
with the microcanonical

328
00:16:39,800 --> 00:16:41,220
equation of motion.

329
00:16:41,220 --> 00:16:43,340
But every now and
then, you give little

330
00:16:43,340 --> 00:16:48,350
kicks to random atoms in which
those atoms either accelerate--

331
00:16:48,350 --> 00:16:50,580
so they start picking
up more velocity.

332
00:16:50,580 --> 00:16:52,760
And so the temperature
of your system rise.

333
00:16:52,760 --> 00:16:57,080
And so you drive your system
to your required temperature.

334
00:16:57,080 --> 00:16:59,870
Or instead of sort of kicking
and accelerating them,

335
00:16:59,870 --> 00:17:03,470
you slow them down sort of
with a kick in the direction

336
00:17:03,470 --> 00:17:05,910
opposite to its velocity.

337
00:17:05,910 --> 00:17:08,250
And so you cool the system down.

338
00:17:08,250 --> 00:17:12,410
And this is really sort of
analogous to the thermal bath

339
00:17:12,410 --> 00:17:13,460
if you want.

340
00:17:13,460 --> 00:17:16,609
You have all your
molecules moving around,

341
00:17:16,609 --> 00:17:20,130
and collisions thermalize them.

342
00:17:20,130 --> 00:17:24,000
So this would be a stochastic
or Langevin approach.

343
00:17:24,000 --> 00:17:28,010
Another approach would be
a constrained approach.

344
00:17:28,010 --> 00:17:31,310
Often, the simplest one
is the velocity rescaling.

345
00:17:31,310 --> 00:17:34,820
Suppose that, again, you want
to keep a constant temperature.

346
00:17:34,820 --> 00:17:37,640
Well, constant
temperature means that you

347
00:17:37,640 --> 00:17:44,480
have fixed overall average
kinetic energy for your ions.

348
00:17:44,480 --> 00:17:48,020
And so when integrating
the microcanonical equation

349
00:17:48,020 --> 00:17:50,390
of motion, at
every time step you

350
00:17:50,390 --> 00:17:54,110
have the velocity on every atom.

351
00:17:54,110 --> 00:17:57,020
And if the average
kinetic energy

352
00:17:57,020 --> 00:18:00,320
is below your target
temperature or above your target

353
00:18:00,320 --> 00:18:02,330
temperature, the only
thing you need to do

354
00:18:02,330 --> 00:18:07,760
is change your velocity in
order to have the constraint

355
00:18:07,760 --> 00:18:11,240
that, at every time step, the
total kinetic energy of the ion

356
00:18:11,240 --> 00:18:12,210
is fixed.

357
00:18:12,210 --> 00:18:15,500
So if you want, a
constrained thermostat

358
00:18:15,500 --> 00:18:20,000
involves keeping the total
kinetic energy of your system

359
00:18:20,000 --> 00:18:21,270
fixed.

360
00:18:21,270 --> 00:18:22,880
And in order to
do this, you need

361
00:18:22,880 --> 00:18:25,790
to have actually an
integrating algorithm

362
00:18:25,790 --> 00:18:28,700
that uses the
velocities to integrate

363
00:18:28,700 --> 00:18:30,380
the equation of motions.

364
00:18:30,380 --> 00:18:32,705
So something like the leap
frog Verlet that we want,

365
00:18:32,705 --> 00:18:36,830
see, or the velocity Verlet that
we want, also see, would work.

366
00:18:36,830 --> 00:18:40,410
But the simple Verlet that
I'll show you again in a moment

367
00:18:40,410 --> 00:18:41,900
wouldn't actually be useful.

368
00:18:41,900 --> 00:18:44,660
Because in the simple
Verlet, you never

369
00:18:44,660 --> 00:18:47,720
control or never
use the velocity

370
00:18:47,720 --> 00:18:50,670
of the system to integrate
the equation of motion.

371
00:18:50,670 --> 00:18:52,970
So you are not able
actually to control

372
00:18:52,970 --> 00:18:54,660
the velocity of your system.

373
00:18:54,660 --> 00:18:59,000
And so you can't thermostat
your dynamical ensemble

374
00:18:59,000 --> 00:19:02,630
using the simple Verlet.

375
00:19:02,630 --> 00:19:05,300
And sort of the last
approach that is probably

376
00:19:05,300 --> 00:19:08,840
the more generic and
the more powerful

377
00:19:08,840 --> 00:19:13,370
is actually adding
to your sort of set

378
00:19:13,370 --> 00:19:17,990
of equation of motion an
additional dynamical variable.

379
00:19:17,990 --> 00:19:20,270
So up to now, the
dynamical variables

380
00:19:20,270 --> 00:19:23,570
are all the position and all
the velocity of your particles

381
00:19:23,570 --> 00:19:24,810
in the system.

382
00:19:24,810 --> 00:19:29,090
Well, you could add an
additional dynamical variable

383
00:19:29,090 --> 00:19:34,310
that I tend to sort of compare
it to some generic wind going

384
00:19:34,310 --> 00:19:35,380
through the system.

385
00:19:35,380 --> 00:19:38,590
And that sort of additional
dynamical variable

386
00:19:38,590 --> 00:19:42,880
is something that either
slows down your particle

387
00:19:42,880 --> 00:19:46,570
or accelerates your
particles in order

388
00:19:46,570 --> 00:19:51,320
to fluctuate appropriately
around your target temperature.

389
00:19:51,320 --> 00:19:54,070
And this is the
most elegant sort

390
00:19:54,070 --> 00:19:56,890
of solution to the problem
of the thermalization,

391
00:19:56,890 --> 00:19:58,430
and it's also the most powerful.

392
00:19:58,430 --> 00:20:01,000
And we'll see examples
both in this class

393
00:20:01,000 --> 00:20:03,740
and especially the
next class of applying

394
00:20:03,740 --> 00:20:06,430
this extended approach.

395
00:20:06,430 --> 00:20:10,240
And often, it goes under the
name of Nosé that was the first

396
00:20:10,240 --> 00:20:14,680
one to develop the equation
of motion for a canonical

397
00:20:14,680 --> 00:20:16,750
simulation at
constant temperature,

398
00:20:16,750 --> 00:20:18,610
or actually Nosé-Hoover.

399
00:20:18,610 --> 00:20:22,420
That was the other person that
made this equation a little bit

400
00:20:22,420 --> 00:20:24,430
more tractable.

401
00:20:24,430 --> 00:20:29,090
And again, a reminder
of the simple Verlet,

402
00:20:29,090 --> 00:20:32,440
the only algorithm that we
have seen, that is sort of

403
00:20:32,440 --> 00:20:34,790
summarized in here.

404
00:20:34,790 --> 00:20:36,550
And what we are doing
is we are looking

405
00:20:36,550 --> 00:20:40,690
at this trajectory,
this vector that

406
00:20:40,690 --> 00:20:43,370
evolves as a function of time.

407
00:20:43,370 --> 00:20:46,810
And we are just doing
a Taylor expansion.

408
00:20:46,810 --> 00:20:50,380
That is, if we know
the position at time t,

409
00:20:50,380 --> 00:20:54,550
the position at times t plus
delta t, the Taylor expansion,

410
00:20:54,550 --> 00:20:57,970
involves the derivative of the
position with respect to time.

411
00:20:57,970 --> 00:20:59,878
The second derivative
of the position

412
00:20:59,878 --> 00:21:01,420
with respect to
time, that is nothing

413
00:21:01,420 --> 00:21:03,880
than the acceleration,
the third derivative

414
00:21:03,880 --> 00:21:05,540
and so on and so forth.

415
00:21:05,540 --> 00:21:08,950
And I've only expanded
it up to the third order.

416
00:21:08,950 --> 00:21:11,800
But the very nice thing
is that we can expand it

417
00:21:11,800 --> 00:21:14,860
both forward in
time and backward

418
00:21:14,860 --> 00:21:18,910
in time around the
instantaneous position.

419
00:21:18,910 --> 00:21:22,600
And so using, you see, the
instantaneous velocity,

420
00:21:22,600 --> 00:21:26,050
the instantaneous acceleration,
and the instantaneous third

421
00:21:26,050 --> 00:21:31,090
derivative, we are able to
correlate the present position

422
00:21:31,090 --> 00:21:34,570
with the next position in
time and the previous position

423
00:21:34,570 --> 00:21:40,420
in time where delta t is this
fundamental discretization unit

424
00:21:40,420 --> 00:21:42,490
that we call the time step.

425
00:21:42,490 --> 00:21:45,520
And that is the
unit that we choose

426
00:21:45,520 --> 00:21:49,660
to integrate our equation
of motion accurately enough.

427
00:21:49,660 --> 00:21:54,700
And just by taking these two
expression and summing them,

428
00:21:54,700 --> 00:22:00,760
we see that all the
terms cancel out.

429
00:22:00,760 --> 00:22:05,230
And just by sort of moving left
and right the appropriate term,

430
00:22:05,230 --> 00:22:08,410
we have just, by summing
these two Taylor expansion,

431
00:22:08,410 --> 00:22:10,420
there's a fundamental
expression that

432
00:22:10,420 --> 00:22:14,650
is really the core of the
simple Verlet algorithm.

433
00:22:14,650 --> 00:22:17,050
And you see, this is the
recipe that we use, actually,

434
00:22:17,050 --> 00:22:20,530
to integrate our equation
of motion in time.

435
00:22:20,530 --> 00:22:24,760
That is what we are saying
is that the position at times

436
00:22:24,760 --> 00:22:29,680
t plus delta t is related
to the position of time t

437
00:22:29,680 --> 00:22:34,150
and the position at
times t minus delta t.

438
00:22:34,150 --> 00:22:38,800
And so by knowing this position
and the previous position

439
00:22:38,800 --> 00:22:41,860
and by knowing what
is the acceleration

440
00:22:41,860 --> 00:22:43,450
at the present time--

441
00:22:43,450 --> 00:22:45,680
and the acceleration
is nothing else

442
00:22:45,680 --> 00:22:48,010
than the force
divided by the mass.

443
00:22:48,010 --> 00:22:51,830
Here is where Newton's equation
of motion comes into play.

444
00:22:51,830 --> 00:22:54,340
We are able to
predict what is going

445
00:22:54,340 --> 00:22:58,180
to be the next position
with an accuracy that

446
00:22:58,180 --> 00:23:03,310
is of the order of the fourth
power of your time step.

447
00:23:03,310 --> 00:23:07,930
So if you want, if you are
going to have your time step,

448
00:23:07,930 --> 00:23:13,310
the accuracy that you obtain
is going to be 16 times higher.

449
00:23:13,310 --> 00:23:15,550
So one of the fundamental
sort of checks

450
00:23:15,550 --> 00:23:18,490
that you want to do in your
molecular dynamic simulation

451
00:23:18,490 --> 00:23:22,210
is that your time
step is small enough.

452
00:23:22,210 --> 00:23:25,220
But obviously, it doesn't
need to be too small.

453
00:23:25,220 --> 00:23:29,050
And again, we'll see an example
on how you can check for this.

454
00:23:29,050 --> 00:23:32,710
And you see the simple
Verlet algorithm never uses

455
00:23:32,710 --> 00:23:34,510
the velocity of your system.

456
00:23:34,510 --> 00:23:37,690
It's able to calculate the
new position just using

457
00:23:37,690 --> 00:23:41,410
the present acceleration
of the present force.

458
00:23:41,410 --> 00:23:43,930
And so it's not
an algorithm that

459
00:23:43,930 --> 00:23:48,680
allows you to control the
velocity of your particles.

460
00:23:48,680 --> 00:23:50,560
So the velocity of
your particle is

461
00:23:50,560 --> 00:23:53,930
something that is given to
you and you can, of course,

462
00:23:53,930 --> 00:23:54,880
calculate.

463
00:23:54,880 --> 00:23:56,920
Suppose that you
want to calculate

464
00:23:56,920 --> 00:23:59,650
what is the kinetic
energy of your system.

465
00:23:59,650 --> 00:24:01,870
Suppose that you
want to find out what

466
00:24:01,870 --> 00:24:03,430
is the average temperature.

467
00:24:03,430 --> 00:24:06,670
Well, you need to
calculate the velocity.

468
00:24:06,670 --> 00:24:12,500
And again, you can use a simple
finite difference formula.

469
00:24:12,500 --> 00:24:15,580
And you see, what we use is
actually the position of t

470
00:24:15,580 --> 00:24:18,940
plus delta that's written
here and the position at t

471
00:24:18,940 --> 00:24:20,950
minus delta that
is written here.

472
00:24:20,950 --> 00:24:23,500
And these are
symmetric derivative.

473
00:24:23,500 --> 00:24:25,390
And you could actually
do the homework

474
00:24:25,390 --> 00:24:29,590
of inserting these
quantities and finding out

475
00:24:29,590 --> 00:24:34,330
what is the order of error
in a power of delta t

476
00:24:34,330 --> 00:24:36,710
for your velocity
and why we use,

477
00:24:36,710 --> 00:24:39,620
actually, this
symmetric formula,

478
00:24:39,620 --> 00:24:42,800
say, compared to a left
or a right derivative.

479
00:24:42,800 --> 00:24:47,150
You could also calculate the
velocity as rt plus delta

480
00:24:47,150 --> 00:24:51,500
t minus rt divided by delta t.

481
00:24:51,500 --> 00:24:55,310
And it's easy to do it here.

482
00:24:55,310 --> 00:25:06,990
If you want to calculate what
is the velocity at times t,

483
00:25:06,990 --> 00:25:11,130
you obviously need to know
what is the derivative here.

484
00:25:11,130 --> 00:25:13,870
This is the definition
of velocity.

485
00:25:13,870 --> 00:25:17,520
And you can see that the
most accurate definition

486
00:25:17,520 --> 00:25:22,620
of a derivative is obtained
by taking the value at rt

487
00:25:22,620 --> 00:25:28,620
plus delta t and t minus delta t
and measuring this slope rather

488
00:25:28,620 --> 00:25:32,770
than, say, measuring this
slope or measuring this slope.

489
00:25:32,770 --> 00:25:37,020
So this slope here is given
by the symmetric derivative

490
00:25:37,020 --> 00:25:40,290
that I have written
down here and is

491
00:25:40,290 --> 00:25:42,960
the sort of most
accurate representation

492
00:25:42,960 --> 00:25:46,380
of the true slope that is
the green position here.

493
00:25:46,380 --> 00:25:49,290
And you can also do this not
only sort of graphically,

494
00:25:49,290 --> 00:25:55,440
but analytically, just writing
out the explicit expressions.

495
00:26:00,290 --> 00:26:04,220
I won't go into the details of
other integration algorithm.

496
00:26:04,220 --> 00:26:06,890
But just to mention
another very common one

497
00:26:06,890 --> 00:26:11,270
is what is called the Gear
predictor-corrector algorithm.

498
00:26:11,270 --> 00:26:14,270
This is actually
not time invariant.

499
00:26:14,270 --> 00:26:19,640
And it tends to be probably
more accurate than the Verlet

500
00:26:19,640 --> 00:26:23,660
algorithm if you have
an analytical expression

501
00:26:23,660 --> 00:26:26,090
for your forces, so
if you can actually

502
00:26:26,090 --> 00:26:29,720
derive from a potential
field your force

503
00:26:29,720 --> 00:26:32,240
and then obtain
your acceleration.

504
00:26:32,240 --> 00:26:37,100
It tends to be much worse
than the Verlet algorithm

505
00:26:37,100 --> 00:26:40,310
if you have noise in your data.

506
00:26:40,310 --> 00:26:41,930
And that's actually
the case when

507
00:26:41,930 --> 00:26:45,710
you perform an ab initio
molecular dynamic simulation

508
00:26:45,710 --> 00:26:48,590
because the forces that
you obtain in an ab initio

509
00:26:48,590 --> 00:26:53,750
simulation are exact only,
if at each time step,

510
00:26:53,750 --> 00:26:57,510
you have reached the
perfect self-consistency.

511
00:26:57,510 --> 00:26:59,690
And so ab initio
molecular dynamics

512
00:26:59,690 --> 00:27:05,270
has some inherent
amount of noise

513
00:27:05,270 --> 00:27:09,620
that calls for the use of a more
robust integration algorithm,

514
00:27:09,620 --> 00:27:11,450
like this simple Verlet.

515
00:27:11,450 --> 00:27:13,880
But if you are doing
classical molecular dynamics,

516
00:27:13,880 --> 00:27:15,830
this is often superior.

517
00:27:15,830 --> 00:27:18,290
And really it's sort of
an extension of the Verlet

518
00:27:18,290 --> 00:27:19,130
algorithm.

519
00:27:19,130 --> 00:27:22,880
That is, instead of
using just the first

520
00:27:22,880 --> 00:27:25,310
and the second derivatives
of your trajectories,

521
00:27:25,310 --> 00:27:28,520
it uses the third, the fourth,
and the fifth derivative,

522
00:27:28,520 --> 00:27:30,240
and so on and so forth.

523
00:27:30,240 --> 00:27:33,260
And then sort of it
estimates the error

524
00:27:33,260 --> 00:27:36,740
and corrects, once you
have taken your time step,

525
00:27:36,740 --> 00:27:40,460
for the difference
between your predicted,

526
00:27:40,460 --> 00:27:42,560
say, acceleration
or third derivative

527
00:27:42,560 --> 00:27:46,370
or fourth derivative that comes
from your Taylor expansion

528
00:27:46,370 --> 00:27:49,010
and actually, say, the
acceleration that you

529
00:27:49,010 --> 00:27:50,000
calculate.

530
00:27:50,000 --> 00:27:51,920
That is the force
that you calculate

531
00:27:51,920 --> 00:27:53,750
once you move in that position.

532
00:27:53,750 --> 00:27:57,830
And again, using higher
order expressions--

533
00:27:57,830 --> 00:28:00,200
that is using third and
fourth derivatives--

534
00:28:00,200 --> 00:28:04,610
can give you more accuracy and
allow you to use longer time

535
00:28:04,610 --> 00:28:08,540
step to integrate appropriately
your trajectory if you

536
00:28:08,540 --> 00:28:10,550
have no noise in your forces.

537
00:28:10,550 --> 00:28:12,560
But if you have
noise in your forces,

538
00:28:12,560 --> 00:28:15,830
it can actually be more
catastrophic in the errors

539
00:28:15,830 --> 00:28:18,040
that it obtains.

540
00:28:18,040 --> 00:28:18,950
OK.

541
00:28:18,950 --> 00:28:22,610
This was a reflection of some
of the things that we have seen.

542
00:28:22,610 --> 00:28:25,460
One of the sort of
very simple concept

543
00:28:25,460 --> 00:28:29,000
I wanted to sort
of remind you that

544
00:28:29,000 --> 00:28:31,670
takes place whenever
we try to do

545
00:28:31,670 --> 00:28:35,060
numerical integration of
a differential equation

546
00:28:35,060 --> 00:28:37,910
is the concept of
Lyapunov instability.

547
00:28:37,910 --> 00:28:41,330
That is sort of what you have
heard over and over again

548
00:28:41,330 --> 00:28:45,360
called the butterfly
effect in weather forecast.

549
00:28:45,360 --> 00:28:51,060
That is, if a butterfly beats
his wings in northern Africa,

550
00:28:51,060 --> 00:28:53,600
we have a tornado
in the Caribbean.

551
00:28:53,600 --> 00:28:56,390
And the reason
for that is saying

552
00:28:56,390 --> 00:29:00,260
that, in any of this
differential equation, when you

553
00:29:00,260 --> 00:29:04,790
try to integrate them,
an infinitesimal error

554
00:29:04,790 --> 00:29:06,890
at the beginning
of your integration

555
00:29:06,890 --> 00:29:10,190
will exponentially
increase to the point

556
00:29:10,190 --> 00:29:14,810
that what comes at the
end of your simulation

557
00:29:14,810 --> 00:29:18,440
can be radically different
for infinitesimal differences

558
00:29:18,440 --> 00:29:20,120
in your initial condition.

559
00:29:20,120 --> 00:29:22,520
That is to say
suppose that you have

560
00:29:22,520 --> 00:29:24,530
an exact analytical trajectory.

561
00:29:24,530 --> 00:29:26,660
That is you are
able to integrate

562
00:29:26,660 --> 00:29:31,850
your equation of motion starting
from a certain initial position

563
00:29:31,850 --> 00:29:34,320
and initial velocity here.

564
00:29:34,320 --> 00:29:39,470
So in the limit of an
infinitesimal time step,

565
00:29:39,470 --> 00:29:43,310
in principle you are going
to do a perfect integration

566
00:29:43,310 --> 00:29:45,560
and recover this trajectory.

567
00:29:45,560 --> 00:29:51,830
In practice, it means that,
for larger and larger steps,

568
00:29:51,830 --> 00:29:53,390
what you will have
is that you'll

569
00:29:53,390 --> 00:29:58,370
be able to sort of follow this
trajectory for a certain time.

570
00:29:58,370 --> 00:30:03,080
But after a while,
you'll start diverging.

571
00:30:03,080 --> 00:30:06,950
And the smaller
your time step is,

572
00:30:06,950 --> 00:30:08,780
the more accurately
you are going

573
00:30:08,780 --> 00:30:13,160
to reproduce this
trajectory for a while.

574
00:30:13,160 --> 00:30:17,390
But at the end, the
exponential explosion

575
00:30:17,390 --> 00:30:18,620
will always take place.

576
00:30:18,620 --> 00:30:23,720
That is you start having
an exponentially increasing

577
00:30:23,720 --> 00:30:28,400
difference between the exact
trajectory and the trajectory

578
00:30:28,400 --> 00:30:31,070
that you are calculating.

579
00:30:31,070 --> 00:30:34,290
Luckily, this doesn't
really matter.

580
00:30:34,290 --> 00:30:36,320
This is very important.

581
00:30:36,320 --> 00:30:40,370
And this is sort of true in
molecular dynamic simulation

582
00:30:40,370 --> 00:30:45,020
of our small sample periodically
repeated in, which we want

583
00:30:45,020 --> 00:30:47,210
to calculate thermal averages.

584
00:30:47,210 --> 00:30:51,110
And it also applies to
the case of the butterfly

585
00:30:51,110 --> 00:30:53,660
beating its wings and
creating a tornado.

586
00:30:53,660 --> 00:30:58,070
What really happens is that all
these systems on the long time

587
00:30:58,070 --> 00:31:02,860
scales are basically
inherently chaotic.

588
00:31:02,860 --> 00:31:05,970
So all these
trajectories diverge

589
00:31:05,970 --> 00:31:10,710
from the ideal trajectory,
but it doesn't really matter.

590
00:31:10,710 --> 00:31:13,380
When we do a molecular
dynamic simulation,

591
00:31:13,380 --> 00:31:18,330
we are not trying to reproduce
an exact trajectory starting

592
00:31:18,330 --> 00:31:22,170
from the atomic timescales and
propagating for picoseconds,

593
00:31:22,170 --> 00:31:24,480
nanoseconds, or microseconds.

594
00:31:24,480 --> 00:31:30,180
Because what we want to find out
are either average properties

595
00:31:30,180 --> 00:31:31,500
of the system--

596
00:31:31,500 --> 00:31:36,210
well, it doesn't really matter
that the trajectory diverges.

597
00:31:36,210 --> 00:31:39,450
Because even this
trajectory that diverge

598
00:31:39,450 --> 00:31:44,130
still conserve the fundamental
macroscopic variables.

599
00:31:44,130 --> 00:31:47,700
If you are doing a simulation
at constant temperature, what

600
00:31:47,700 --> 00:31:50,550
it really matters is that
all these trajectories

601
00:31:50,550 --> 00:31:54,630
correspond to evolutions
at constant temperature.

602
00:31:54,630 --> 00:31:58,960
And then you are sampling
appropriately your phase space.

603
00:31:58,960 --> 00:32:02,550
And even if you want to do
molecular dynamics of an event

604
00:32:02,550 --> 00:32:05,850
instead of using it to
obtain an assemble average,

605
00:32:05,850 --> 00:32:07,950
even if you want
to, say, find out

606
00:32:07,950 --> 00:32:12,870
how an atom diffuses across
a barrier, what is really

607
00:32:12,870 --> 00:32:16,260
your macroscopic
observable is the average

608
00:32:16,260 --> 00:32:20,880
over tens, thousands,
millions of these events,

609
00:32:20,880 --> 00:32:24,000
all sort of slightly
different from each other,

610
00:32:24,000 --> 00:32:27,420
but all corresponding to
the overall constraint

611
00:32:27,420 --> 00:32:30,910
of having the same
macroscopic observable.

612
00:32:30,910 --> 00:32:35,010
So it doesn't really matter
that we replicate exactly

613
00:32:35,010 --> 00:32:38,760
the trajectory,
neither for an ensemble

614
00:32:38,760 --> 00:32:42,420
average of thermodynamical
quantities, nor for the case

615
00:32:42,420 --> 00:32:45,240
in which we simulate
a specific event.

616
00:32:45,240 --> 00:32:47,430
Because in any case,
what we need to do

617
00:32:47,430 --> 00:32:50,850
is repeat the simulation
over and over again

618
00:32:50,850 --> 00:32:54,510
to have appropriate statistics
and to have an average.

619
00:32:54,510 --> 00:32:57,630
Of course, the only case where
having the exact trajectory

620
00:32:57,630 --> 00:32:59,850
matters is that
if you are sending

621
00:32:59,850 --> 00:33:02,790
a satellite in outer
space because you are not

622
00:33:02,790 --> 00:33:07,260
sending 1 million satellites and
seeing what they do on average.

623
00:33:07,260 --> 00:33:10,980
You want that satellite to
do something very precise.

624
00:33:10,980 --> 00:33:13,680
And so you need to integrate
things very accurately.

625
00:33:13,680 --> 00:33:16,980
And you need to keep
correcting your integration

626
00:33:16,980 --> 00:33:19,260
as time goes by.

627
00:33:19,260 --> 00:33:22,350
And we can actually
sort of formalize

628
00:33:22,350 --> 00:33:26,040
this sort of concept of
Lyapunov instabilities

629
00:33:26,040 --> 00:33:30,030
in terms, basically, of
an exponential divergence

630
00:33:30,030 --> 00:33:33,210
in some metric of
our actual trajectory

631
00:33:33,210 --> 00:33:38,350
from the ideal analytic
trajectory of our system.

632
00:33:38,350 --> 00:33:38,850
OK.

633
00:33:38,850 --> 00:33:41,460
So now, everything is in place.

634
00:33:41,460 --> 00:33:43,950
We are doing our molecular
dynamic simulation.

635
00:33:43,950 --> 00:33:47,190
First of all, how do
we check that something

636
00:33:47,190 --> 00:33:49,890
is working, that what we
are doing makes sense?

637
00:33:49,890 --> 00:33:53,820
Well, the first thing
is just making sure

638
00:33:53,820 --> 00:33:57,180
that what you are doing
makes sense computationally.

639
00:33:57,180 --> 00:33:59,940
It might still not
make sense physically,

640
00:33:59,940 --> 00:34:02,580
but your sort of
order [INAUDIBLE]

641
00:34:02,580 --> 00:34:06,420
is making sure that you are
integrating your equation

642
00:34:06,420 --> 00:34:08,219
of motion appropriately.

643
00:34:08,219 --> 00:34:10,409
Because, remember, you
need to choose a time step.

644
00:34:10,409 --> 00:34:14,190
You need to choose a
discretization interval.

645
00:34:14,190 --> 00:34:18,239
And the size of that
interval depends really

646
00:34:18,239 --> 00:34:20,730
on what is the
typical atomic motion.

647
00:34:20,730 --> 00:34:23,610
If you are studying a
system that is very cool

648
00:34:23,610 --> 00:34:26,219
and is made of very
heavy atoms, you

649
00:34:26,219 --> 00:34:28,170
can use very large times steps.

650
00:34:28,170 --> 00:34:31,380
If your system is very hot and
made of very light atoms that

651
00:34:31,380 --> 00:34:35,790
move much faster, you need
to use very small time steps.

652
00:34:35,790 --> 00:34:39,630
But on average, there is really
a small range of time steps

653
00:34:39,630 --> 00:34:42,090
that is really acceptable.

654
00:34:42,090 --> 00:34:45,300
And they are of the order
of the femtoseconds, 10

655
00:34:45,300 --> 00:34:47,880
to the minus 15 seconds.

656
00:34:47,880 --> 00:34:52,170
So for something as light
as a hydrogen atom bound

657
00:34:52,170 --> 00:34:54,480
to an oxygen, say,
in water molecule,

658
00:34:54,480 --> 00:34:56,670
you might need to
use a time step as

659
00:34:56,670 --> 00:34:59,970
small as 1/2 femtosecond.

660
00:34:59,970 --> 00:35:02,710
For something in which
maybe you are studying,

661
00:35:02,710 --> 00:35:06,630
say, liquid caesium,
a free electron metal

662
00:35:06,630 --> 00:35:08,190
that is very heavy,
you might have

663
00:35:08,190 --> 00:35:10,530
a time step that
can be appropriate

664
00:35:10,530 --> 00:35:13,380
and as large as 20 femtoseconds.

665
00:35:13,380 --> 00:35:15,420
But this is the
range, basically,

666
00:35:15,420 --> 00:35:18,210
of the order of femtoseconds.

667
00:35:18,210 --> 00:35:21,330
And the very important
thing is that you can always

668
00:35:21,330 --> 00:35:26,160
check that your time
step is accurate enough.

669
00:35:26,160 --> 00:35:29,490
Because if you are doing a
microcanonical simulation,

670
00:35:29,490 --> 00:35:33,300
what you need to have is
that the total energy--

671
00:35:33,300 --> 00:35:37,170
that is the sum of the potential
energy and the kinetic energy--

672
00:35:37,170 --> 00:35:40,920
is conserved or, if you
want, that the kinetic energy

673
00:35:40,920 --> 00:35:44,340
of the ions and the
potential energy of the ions

674
00:35:44,340 --> 00:35:47,280
mirror each other exactly.

675
00:35:47,280 --> 00:35:50,610
And again, in the limit of
infinitesimal time step,

676
00:35:50,610 --> 00:35:55,830
this total energy line
should be perfectly flat.

677
00:35:55,830 --> 00:36:01,140
As you increase your time
step, well, this line

678
00:36:01,140 --> 00:36:05,870
will start to have some
microscopic oscillation.

679
00:36:05,870 --> 00:36:10,780
So it will start fluctuate and,
depending on your algorithm,

680
00:36:10,780 --> 00:36:13,960
could also have a
systematic drift.

681
00:36:13,960 --> 00:36:16,360
A systematic drift,
generally speaking,

682
00:36:16,360 --> 00:36:21,730
is very bad because it means
that your system is really

683
00:36:21,730 --> 00:36:22,930
cooling down.

684
00:36:22,930 --> 00:36:27,190
That is a systematic drift
means that your total energy is

685
00:36:27,190 --> 00:36:30,740
going to get closer and closer
to the potential energy.

686
00:36:30,740 --> 00:36:33,550
So the kinetic
energy of your ions

687
00:36:33,550 --> 00:36:35,750
is decreasing and decreasing.

688
00:36:35,750 --> 00:36:39,370
And so systematic drift
means that, after a while,

689
00:36:39,370 --> 00:36:42,890
your system will
slow down to a halt.

690
00:36:42,890 --> 00:36:45,010
And so you need
to check for that.

691
00:36:45,010 --> 00:36:48,070
And then even if there
is no systematic drift,

692
00:36:48,070 --> 00:36:50,740
a simple algorithm
like the Verlet

693
00:36:50,740 --> 00:36:53,470
tend to be very good
at that, at not having

694
00:36:53,470 --> 00:36:55,000
any systematic drift.

695
00:36:55,000 --> 00:36:58,690
You need to make sure that your
fluctuation that will always

696
00:36:58,690 --> 00:37:02,200
exist as long as you are
using a finite time step

697
00:37:02,200 --> 00:37:04,560
are actually small enough.

698
00:37:04,560 --> 00:37:07,590
So you want the fluctuation,
your total energy,

699
00:37:07,590 --> 00:37:10,810
to be much, much smaller
than the fluctuation

700
00:37:10,810 --> 00:37:14,170
in your potential energy,
say two orders of magnitude,

701
00:37:14,170 --> 00:37:16,690
three orders of
magnitude smaller.

702
00:37:16,690 --> 00:37:23,960
So that locally your system is
varying its potential energy,

703
00:37:23,960 --> 00:37:28,030
but its total energy
is actually unchanged.

704
00:37:28,030 --> 00:37:30,730
And again, thanks to the
power of the Verlet algorithm

705
00:37:30,730 --> 00:37:34,810
that sort of scales as the
fourth power of your time step,

706
00:37:34,810 --> 00:37:37,000
you can do very simple tests.

707
00:37:37,000 --> 00:37:41,860
You half your time step, and
you should see your total energy

708
00:37:41,860 --> 00:37:44,560
curve becomes much flatter.

709
00:37:44,560 --> 00:37:47,380
And so you can basically
choose sort of what

710
00:37:47,380 --> 00:37:49,300
is the optimal time step.

711
00:37:49,300 --> 00:37:55,450
And you can even use a sort of
real space picture in which you

712
00:37:55,450 --> 00:37:56,800
can look at the trajectory.

713
00:37:56,800 --> 00:37:58,870
After all, what
you have are atoms

714
00:37:58,870 --> 00:38:01,360
moving around and oscillating.

715
00:38:01,360 --> 00:38:05,530
And so what you want to do in
a real trajectory of an atom

716
00:38:05,530 --> 00:38:08,440
moving around is
you want to have

717
00:38:08,440 --> 00:38:13,960
enough points to
describe correctly, say,

718
00:38:13,960 --> 00:38:16,180
a period of an atom.

719
00:38:16,180 --> 00:38:23,530
And the rule a thumb is that you
want between 20 to 100 points

720
00:38:23,530 --> 00:38:26,470
in describing a full
period of an atom.

721
00:38:26,470 --> 00:38:29,050
And that, again,
sort of brings us

722
00:38:29,050 --> 00:38:33,970
back to the order of magnitude
of femtoseconds for the time

723
00:38:33,970 --> 00:38:37,780
steps of your systems involved.

724
00:38:37,780 --> 00:38:44,170
So once you have set up your
computational part correctly,

725
00:38:44,170 --> 00:38:49,070
you still need to check out that
your simulation makes sense.

726
00:38:49,070 --> 00:38:52,150
And in particular, the
sort of trickiest part

727
00:38:52,150 --> 00:38:56,450
is making sure that you reach
thermodynamic equilibrium,

728
00:38:56,450 --> 00:39:00,550
so that you lose memory of
your initial conditions.

729
00:39:00,550 --> 00:39:02,380
And again, that's
actually something

730
00:39:02,380 --> 00:39:05,470
in which the Lyapunov
instabilities, if you want,

731
00:39:05,470 --> 00:39:06,430
help.

732
00:39:06,430 --> 00:39:12,880
No matter how bad our choice
is, if we let the system evolve

733
00:39:12,880 --> 00:39:17,230
long enough, our
trajectory will sort of

734
00:39:17,230 --> 00:39:22,790
diverge from our
initial exact trajectory

735
00:39:22,790 --> 00:39:25,280
that was somehow poorly chosen.

736
00:39:25,280 --> 00:39:28,250
And it will become more
and more representative

737
00:39:28,250 --> 00:39:32,150
of the thermodynamical ensemble
that we are simulating.

738
00:39:32,150 --> 00:39:35,300
But it's still very
tricky to understand

739
00:39:35,300 --> 00:39:37,130
if we are reaching equilibrium.

740
00:39:37,130 --> 00:39:40,910
I've shown you here an example
in the case in which, say, we

741
00:39:40,910 --> 00:39:46,040
are looking at a solid,
something like a slab of metal

742
00:39:46,040 --> 00:39:49,790
in which we have a, say,
bulk atom in the metal.

743
00:39:49,790 --> 00:39:52,490
And then as we move
towards the surface,

744
00:39:52,490 --> 00:39:56,120
we have sort of second layer
atom and first layer atom.

745
00:39:56,120 --> 00:39:58,160
And what we are
averaging here, what

746
00:39:58,160 --> 00:40:02,120
is the sort of thermodynamical
observable we are interested,

747
00:40:02,120 --> 00:40:04,160
is the mean square displacement.

748
00:40:04,160 --> 00:40:07,430
That is how much
the atom fluctuates

749
00:40:07,430 --> 00:40:09,420
around the initial condition.

750
00:40:09,420 --> 00:40:11,420
And again, you'll see
more and more example

751
00:40:11,420 --> 00:40:15,470
in what comes later of a
mean square displacement.

752
00:40:15,470 --> 00:40:17,600
But you see, somehow
it's written here

753
00:40:17,600 --> 00:40:18,710
in a reverse scale.

754
00:40:18,710 --> 00:40:22,760
But if we start at time
t equal to 0 and we

755
00:40:22,760 --> 00:40:25,700
measure what is the mean
square displacement,

756
00:40:25,700 --> 00:40:30,030
well, at the beginning,
the atoms will keep moving.

757
00:40:30,030 --> 00:40:35,750
If you think, this peak
here, this sort of 1/3 to 1/4

758
00:40:35,750 --> 00:40:38,840
of a picosecond, is
a measure of what

759
00:40:38,840 --> 00:40:42,470
is the period of
oscillation of an atom.

760
00:40:42,470 --> 00:40:45,590
Because at the beginning,
the atom starts moving.

761
00:40:45,590 --> 00:40:48,860
And from the point of
view of calculating

762
00:40:48,860 --> 00:40:53,420
what is its average
displacement, it keeps moving.

763
00:40:53,420 --> 00:40:56,870
At the very beginning of
the simulation, atoms move.

764
00:40:56,870 --> 00:41:00,140
And then they start feeling the
potential from the other atom.

765
00:41:00,140 --> 00:41:03,800
And they sort of settle in
an oscillatory movement.

766
00:41:03,800 --> 00:41:06,650
And then you are
starting to average over

767
00:41:06,650 --> 00:41:09,350
more and more oscillations.

768
00:41:09,350 --> 00:41:12,170
And you see sort
of, in the bulk, you

769
00:41:12,170 --> 00:41:15,440
fairly quickly converge.

770
00:41:15,440 --> 00:41:17,360
These thin lines
are for the bulk.

771
00:41:17,360 --> 00:41:20,570
You very quickly
converge to your sort

772
00:41:20,570 --> 00:41:24,140
of final exact results.

773
00:41:24,140 --> 00:41:27,320
But somewhere else
for the surface atoms,

774
00:41:27,320 --> 00:41:28,640
it takes more time.

775
00:41:28,640 --> 00:41:32,120
It takes a few seconds
before you finally

776
00:41:32,120 --> 00:41:37,670
converge to your
sort of final result.

777
00:41:37,670 --> 00:41:39,850
But again,
thermalization sort of

778
00:41:39,850 --> 00:41:46,700
can be checked by looking at the
time evolution of the quantity

779
00:41:46,700 --> 00:41:48,470
that you want to average.

780
00:41:48,470 --> 00:41:51,320
Say, you could look at the
temperature of your system,

781
00:41:51,320 --> 00:41:54,980
and you want to see
the quantity that you

782
00:41:54,980 --> 00:42:01,520
are sort of monitoring
converge to a final result that

783
00:42:01,520 --> 00:42:06,218
doesn't change as you increase
the length of your simulation.

784
00:42:06,218 --> 00:42:07,760
And there are
quantities, again, that

785
00:42:07,760 --> 00:42:11,360
are sort of easy to
converge and quantities that

786
00:42:11,360 --> 00:42:13,790
are very difficult to converge.

787
00:42:13,790 --> 00:42:17,990
Suppose for a moment that you
want to study actually, again,

788
00:42:17,990 --> 00:42:20,810
the structure of liquid water.

789
00:42:20,810 --> 00:42:24,560
Well, you could start from
a sample that is frozen--

790
00:42:24,560 --> 00:42:28,340
that is ice-- and
all of a sudden

791
00:42:28,340 --> 00:42:31,430
you increase the temperature
with your thermostat.

792
00:42:31,430 --> 00:42:33,950
You bring it, say,
to 300 kelvin.

793
00:42:33,950 --> 00:42:36,530
And then you let it evolve.

794
00:42:36,530 --> 00:42:40,130
And you could sort
of look in time

795
00:42:40,130 --> 00:42:44,210
as the sort of relative
position or the diffusion

796
00:42:44,210 --> 00:42:48,170
coefficient of the water
molecule evolve in time.

797
00:42:48,170 --> 00:42:50,270
And you would discover
that, depending

798
00:42:50,270 --> 00:42:55,220
on the temperature, maybe you
need 10, 20, 30 picoseconds

799
00:42:55,220 --> 00:43:00,920
to reach a situation after which
your expectation value do not

800
00:43:00,920 --> 00:43:04,490
change if you increase your
simulation longer and longer.

801
00:43:04,490 --> 00:43:08,030
This is sort of the
concept of equilibration.

802
00:43:08,030 --> 00:43:10,250
If you were to do
the opposite, if you

803
00:43:10,250 --> 00:43:13,400
wanted to find out what
is the structure of ice

804
00:43:13,400 --> 00:43:18,470
by taking a sample of liquid
water and cooling it down,

805
00:43:18,470 --> 00:43:22,070
your equilibration time
would be almost infinite.

806
00:43:22,070 --> 00:43:26,280
Because if you take liquid
water and you cool it down,

807
00:43:26,280 --> 00:43:30,560
it will actually not freeze
in the perfect crystalline

808
00:43:30,560 --> 00:43:32,720
structure of ice.

809
00:43:32,720 --> 00:43:36,710
But it will much more
likely sort of slow down

810
00:43:36,710 --> 00:43:40,520
in a sort of glassy transition
in which it freezes.

811
00:43:40,520 --> 00:43:45,530
It stop diffusing, but it's
only sort of locally coordinated

812
00:43:45,530 --> 00:43:46,920
as ice should be.

813
00:43:46,920 --> 00:43:49,130
But it's not in a
crystalline state.

814
00:43:49,130 --> 00:43:50,100
What does that mean?

815
00:43:50,100 --> 00:43:53,840
It means that what we
actually ended up with

816
00:43:53,840 --> 00:43:58,100
is a situation that is not
at thermally representative

817
00:43:58,100 --> 00:43:58,730
situation.

818
00:43:58,730 --> 00:44:02,510
We are not in an equilibrium
state for a system,

819
00:44:02,510 --> 00:44:08,120
but we have sort of and ended
up in a corner of phase space

820
00:44:08,120 --> 00:44:11,210
out of which we are
never going to get out

821
00:44:11,210 --> 00:44:14,570
unless we wait forever.

822
00:44:14,570 --> 00:44:16,760
That is, even after
their water is

823
00:44:16,760 --> 00:44:19,790
sort of frozen in
this glassy state,

824
00:44:19,790 --> 00:44:23,450
if you keep simulating it
longer and longer and longer,

825
00:44:23,450 --> 00:44:27,020
every now and then you'll have
a fluctuation that changes

826
00:44:27,020 --> 00:44:29,630
the position of two water
molecules and two more water

827
00:44:29,630 --> 00:44:31,310
molecules and two
water molecules.

828
00:44:31,310 --> 00:44:34,040
And after a very long
time, you might actually

829
00:44:34,040 --> 00:44:36,620
find that your system
has become crystalline.

830
00:44:36,620 --> 00:44:39,800
But it's never going to
happen in the timescale

831
00:44:39,800 --> 00:44:41,280
of your simulation.

832
00:44:41,280 --> 00:44:43,430
So one needs to be
careful when one

833
00:44:43,430 --> 00:44:47,180
does one of these simulation to
make sure that, actually, you

834
00:44:47,180 --> 00:44:51,060
are sort of finding yourself
in a representative state.

835
00:44:51,060 --> 00:44:53,480
And that's why some
of these issues,

836
00:44:53,480 --> 00:44:56,450
like simulating
a solid to liquid

837
00:44:56,450 --> 00:44:59,690
or a liquid to solid transition,
require particular care

838
00:44:59,690 --> 00:45:04,160
and somehow require actually
knowing something of the system

839
00:45:04,160 --> 00:45:05,750
that you want to discover.

840
00:45:05,750 --> 00:45:09,800
And again, there is
an enormous debate

841
00:45:09,800 --> 00:45:16,340
on what are some of the
phases of supercooled water.

842
00:45:16,340 --> 00:45:20,030
That is, if you take
water, and it's very pure,

843
00:45:20,030 --> 00:45:23,870
and you bring it down below the
regular freezing temperature,

844
00:45:23,870 --> 00:45:27,260
well, you can keep
water liquid up to,

845
00:45:27,260 --> 00:45:33,350
I think, 40 degrees
below the freezing point.

846
00:45:33,350 --> 00:45:35,270
And you can do that
experimentally.

847
00:45:35,270 --> 00:45:37,100
And you can do it
with simulation.

848
00:45:37,100 --> 00:45:41,570
And then, all of a sudden, your
system will ultimately freeze.

849
00:45:41,570 --> 00:45:46,490
And it will end up in phases
that neither the people doing

850
00:45:46,490 --> 00:45:49,190
simulation, nor the
people doing experiment

851
00:45:49,190 --> 00:45:50,490
can still agree with.

852
00:45:50,490 --> 00:45:54,980
So there is sort of dark areas
in the phase diagram of water

853
00:45:54,980 --> 00:45:56,210
at those low temperatures.

854
00:45:56,210 --> 00:45:59,600
And people this discuss about
high density liquid and low

855
00:45:59,600 --> 00:46:00,710
density liquid.

856
00:46:00,710 --> 00:46:04,290
And a lot has been achieved
in the last 10 years.

857
00:46:04,290 --> 00:46:07,190
But it's still not
formally known.

858
00:46:07,190 --> 00:46:11,330
And this is why one really needs
to pay particular attention

859
00:46:11,330 --> 00:46:15,060
to a number of these problems.

860
00:46:15,060 --> 00:46:15,560
OK.

861
00:46:15,560 --> 00:46:18,800
So suppose that
everything is going well.

862
00:46:18,800 --> 00:46:21,200
You have initialized
properly your system.

863
00:46:21,200 --> 00:46:24,320
You have sort of
weighted enough,

864
00:46:24,320 --> 00:46:26,960
so that you have
actually lost memory

865
00:46:26,960 --> 00:46:29,040
of your initial condition.

866
00:46:29,040 --> 00:46:31,280
And if you want,
you should always

867
00:46:31,280 --> 00:46:35,570
throw away the part of your
molecular dynamic trajectory

868
00:46:35,570 --> 00:46:38,810
that corresponds to this
initial thermalization time

869
00:46:38,810 --> 00:46:42,380
because it doesn't really
represent high quality data.

870
00:46:42,380 --> 00:46:45,920
Only after you have
waited long enough,

871
00:46:45,920 --> 00:46:48,890
say, in this
previous simulation,

872
00:46:48,890 --> 00:46:52,130
only after you
have waited, say, 4

873
00:46:52,130 --> 00:46:57,110
picosecond, you can start
taking averages of what is, say,

874
00:46:57,110 --> 00:47:01,340
the mean square displacement
for a surface atom.

875
00:47:01,340 --> 00:47:03,830
So you throw away
from your simulation

876
00:47:03,830 --> 00:47:05,990
your thermalization
time, and you

877
00:47:05,990 --> 00:47:08,630
start accumulating averages.

878
00:47:08,630 --> 00:47:11,750
And again, the averages
that could be relevant,

879
00:47:11,750 --> 00:47:14,330
if you are seeing a
microcanonical simulation,

880
00:47:14,330 --> 00:47:16,700
you could look at the
average kinetic energy.

881
00:47:16,700 --> 00:47:19,220
Because that average
kinetic energy

882
00:47:19,220 --> 00:47:22,730
gives you the average
temperature of your system.

883
00:47:22,730 --> 00:47:28,130
Or you could look at the
average pressure in your system.

884
00:47:28,130 --> 00:47:30,980
Or I show you an
example in a moment.

885
00:47:30,980 --> 00:47:33,260
You could look at
the potential energy

886
00:47:33,260 --> 00:47:37,080
or what is sometimes
called the caloric curve.

887
00:47:37,080 --> 00:47:40,040
So if you are actually
simulating a system

888
00:47:40,040 --> 00:47:42,110
at different temperatures--

889
00:47:42,110 --> 00:47:44,420
say, again, suppose
that you have water,

890
00:47:44,420 --> 00:47:47,300
and you are still
in the frozen phase.

891
00:47:47,300 --> 00:47:52,340
Well, the potential energy
will oscillate around

892
00:47:52,340 --> 00:47:54,030
a certain value.

893
00:47:54,030 --> 00:47:58,310
But remember, in order
to sort of melt water--

894
00:47:58,310 --> 00:48:03,410
that is in order to go
at equilibrium from being

895
00:48:03,410 --> 00:48:05,570
an infinitesimal
amount of temperature

896
00:48:05,570 --> 00:48:09,020
below the melting temperature
and an infinitesimal amount

897
00:48:09,020 --> 00:48:11,000
above the melting
temperature-- you

898
00:48:11,000 --> 00:48:13,700
need really to give
to your system what

899
00:48:13,700 --> 00:48:15,740
is called the latent heat.

900
00:48:15,740 --> 00:48:17,900
In this case, the
latent heat of fusion

901
00:48:17,900 --> 00:48:20,990
would be the heat
released when it cools

902
00:48:20,990 --> 00:48:25,040
or the heat that you
give in order to melt it.

903
00:48:25,040 --> 00:48:28,850
And so your average
potential energy

904
00:48:28,850 --> 00:48:32,540
will have a discontinuity as
a function of temperature.

905
00:48:32,540 --> 00:48:34,220
And we'll see it in an example.

906
00:48:34,220 --> 00:48:36,710
And that's really a measure
from your molecular dynamic

907
00:48:36,710 --> 00:48:40,610
simulation of this fundamental
thermodynamical quantity.

908
00:48:40,610 --> 00:48:44,480
Or you could look, say, in
liquid water at the diffusion

909
00:48:44,480 --> 00:48:47,780
coefficient, so the means
square displacements.

910
00:48:47,780 --> 00:48:49,130
And we'll see this.

911
00:48:49,130 --> 00:48:51,740
Or you could look
at the structure,

912
00:48:51,740 --> 00:48:54,870
the geometric structure,
of your liquid.

913
00:48:54,870 --> 00:48:58,010
And again, sort of molecular
dynamics techniques

914
00:48:58,010 --> 00:49:00,410
were actually initially
developed for system

915
00:49:00,410 --> 00:49:03,260
like liquids in which
it was much more

916
00:49:03,260 --> 00:49:07,760
difficult to have analytical
expression for something

917
00:49:07,760 --> 00:49:09,170
like the structure.

918
00:49:09,170 --> 00:49:10,970
I mean, if you are
studying a crystal,

919
00:49:10,970 --> 00:49:13,100
the structure is
somehow trivial.

920
00:49:13,100 --> 00:49:15,920
And you can obtain it
from X-ray experiments,

921
00:49:15,920 --> 00:49:20,670
but you can also obtain it
from total energy calculation,

922
00:49:20,670 --> 00:49:23,640
minimizing the energy with
respect to your volume.

923
00:49:23,640 --> 00:49:26,190
If you are studying
a liquid, again, you

924
00:49:26,190 --> 00:49:30,180
can have information on your
sort of local structure.

925
00:49:30,180 --> 00:49:33,660
That is how far the
molecules keep on average

926
00:49:33,660 --> 00:49:37,020
from each other in the
simulation from diffraction

927
00:49:37,020 --> 00:49:38,160
experiments.

928
00:49:38,160 --> 00:49:43,380
But you don't have an
immediate computational way

929
00:49:43,380 --> 00:49:45,840
to calculate data
like you have done it

930
00:49:45,840 --> 00:49:48,690
for a solid in which you
were, again, minimizing energy

931
00:49:48,690 --> 00:49:49,950
with respect to volume.

932
00:49:49,950 --> 00:49:53,790
What we really need is to do a
molecular dynamic simulation.

933
00:49:53,790 --> 00:49:56,670
And this is sort of the
first of the examples

934
00:49:56,670 --> 00:49:59,430
I want to show you
today of application

935
00:49:59,430 --> 00:50:00,750
of molecular dynamics.

936
00:50:00,750 --> 00:50:05,220
That is one in which we
calculate ensemble average.

937
00:50:05,220 --> 00:50:08,100
And in particular, we,
say, calculator, say,

938
00:50:08,100 --> 00:50:10,410
the structure of liquid water.

939
00:50:10,410 --> 00:50:13,560
I'll also give you an
example, a couple of examples,

940
00:50:13,560 --> 00:50:16,350
of a chemical
reaction that is using

941
00:50:16,350 --> 00:50:18,330
molecular dynamics to
look at the evolution

942
00:50:18,330 --> 00:50:20,220
in real time of a system.

943
00:50:20,220 --> 00:50:24,930
And I'll sort of hint at
this more general idea

944
00:50:24,930 --> 00:50:27,660
that we can use these
thermodynamical analogies

945
00:50:27,660 --> 00:50:30,360
and sort of molecular
dynamic analogies

946
00:50:30,360 --> 00:50:34,560
to actually solve complex
optimization problems.

947
00:50:34,560 --> 00:50:35,070
OK.

948
00:50:35,070 --> 00:50:39,750
And so this is an example from
a molecular dynamic simulation,

949
00:50:39,750 --> 00:50:42,690
actually an ab initio
molecular dynamic simulation,

950
00:50:42,690 --> 00:50:45,420
and one in which one
tries to figure out

951
00:50:45,420 --> 00:50:49,450
what is the local
structure of the liquid.

952
00:50:49,450 --> 00:50:53,460
And so how do you actually
define structure in a liquid?

953
00:50:53,460 --> 00:50:56,580
Well, sort of one of the
fundamental quantities

954
00:50:56,580 --> 00:50:59,670
is what is called the
pair correlation function.

955
00:50:59,670 --> 00:51:02,160
Suppose that we are
looking, again, at water.

956
00:51:02,160 --> 00:51:04,600
We have the water molecules.

957
00:51:04,600 --> 00:51:08,260
So say we have an oxygen
atom on each molecule.

958
00:51:08,260 --> 00:51:11,340
And so what we can
sort of ask ourselves

959
00:51:11,340 --> 00:51:16,230
is, what is the probability
of finding another oxygen

960
00:51:16,230 --> 00:51:21,120
atom at a certain
distance from one oxygen?

961
00:51:21,120 --> 00:51:23,640
So you could sort
of imagine yourself

962
00:51:23,640 --> 00:51:27,730
sitting on a water molecule
sitting on that oxygen

963
00:51:27,730 --> 00:51:31,110
and sort of
monitoring on average

964
00:51:31,110 --> 00:51:35,670
what is the distribution of
all the other oxygen atoms.

965
00:51:35,670 --> 00:51:38,970
And this would be one of these
pair correlation function.

966
00:51:38,970 --> 00:51:41,580
This is specifically
for the case of water.

967
00:51:41,580 --> 00:51:43,590
I'll tell you in a moment
what the difference

968
00:51:43,590 --> 00:51:46,740
between these curves is,
but what we are plotting

969
00:51:46,740 --> 00:51:49,320
is the probability
of finding an oxygen

970
00:51:49,320 --> 00:51:52,260
atom at a certain
distance if I'm

971
00:51:52,260 --> 00:51:55,200
sort of sitting on one
of these water molecules

972
00:51:55,200 --> 00:51:56,770
and running around.

973
00:51:56,770 --> 00:52:00,450
And as you can see, the
probability below roughly 2.3,

974
00:52:00,450 --> 00:52:03,150
2.4 angstrom is 0.

975
00:52:03,150 --> 00:52:08,982
So you never find an oxygen atom
closer than 2.3, 2.4 angstrom.

976
00:52:08,982 --> 00:52:10,440
And this is because,
obviously, you

977
00:52:10,440 --> 00:52:13,260
know, molecules
by Pauli repulsion

978
00:52:13,260 --> 00:52:16,950
and electrostatic repulsion,
sort of repel each other very

979
00:52:16,950 --> 00:52:18,930
strongly at close distance.

980
00:52:18,930 --> 00:52:21,120
This, in the
Lennard-Jones potential,

981
00:52:21,120 --> 00:52:23,400
was the 1 over r
to the 12th term.

982
00:52:23,400 --> 00:52:25,320
When you get very
close, the repulsion,

983
00:52:25,320 --> 00:52:29,970
it's so strong that you never
get very, very, very close.

984
00:52:29,970 --> 00:52:34,920
Actually, the only time you
sort overcome this short range

985
00:52:34,920 --> 00:52:38,070
repulsion is in something
like nuclear fusion

986
00:52:38,070 --> 00:52:42,330
in which you can manage to
get, say, hydrogen atom to get

987
00:52:42,330 --> 00:52:45,160
so close together that
actually their nuclei melt

988
00:52:45,160 --> 00:52:46,950
into a helium.

989
00:52:46,950 --> 00:52:50,640
But you see, at this standard
condition of pressure

990
00:52:50,640 --> 00:52:52,950
and temperature,
really the probability

991
00:52:52,950 --> 00:52:57,780
of finding some oxygen atom
up to 2.3 angstrom is 0.

992
00:52:57,780 --> 00:53:00,330
And then all of a
sudden, the probability

993
00:53:00,330 --> 00:53:03,390
of finding an oxygen
atom skyrockets.

994
00:53:03,390 --> 00:53:04,510
What does this mean?

995
00:53:04,510 --> 00:53:08,820
Well, it means that,
actually, to have,

996
00:53:08,820 --> 00:53:13,350
if you add an oxygen
atom in a water molecule,

997
00:53:13,350 --> 00:53:19,080
it's very likely to have other
water molecules surrounding you

998
00:53:19,080 --> 00:53:22,800
at a distance that is a
tad shy of 3 angstrom.

999
00:53:22,800 --> 00:53:23,370
OK.

1000
00:53:23,370 --> 00:53:26,830
This is what the chemists would
call the first solvation shell.

1001
00:53:26,830 --> 00:53:27,510
OK.

1002
00:53:27,510 --> 00:53:31,440
The water molecule sort of
runs around in the liquid.

1003
00:53:31,440 --> 00:53:35,040
And there's a certain
number of water molecules

1004
00:53:35,040 --> 00:53:38,070
that are at a certain
average distance that

1005
00:53:38,070 --> 00:53:41,460
is around 2.6, 2.8 angstrom.

1006
00:53:41,460 --> 00:53:43,740
And so this is a
solvation shell,

1007
00:53:43,740 --> 00:53:47,100
a shell of water molecules
that, on average, surround

1008
00:53:47,100 --> 00:53:48,540
this molecule.

1009
00:53:48,540 --> 00:53:52,560
And because there is a sort
of average shell of water

1010
00:53:52,560 --> 00:53:57,000
molecules around, once you
move outwards from this oxygen,

1011
00:53:57,000 --> 00:54:00,720
then there is going to be a
depletion of water molecules.

1012
00:54:00,720 --> 00:54:03,960
Because there is one shell
at a certain distance.

1013
00:54:03,960 --> 00:54:06,240
You move further,
it becomes much more

1014
00:54:06,240 --> 00:54:09,120
difficult to have
water molecules

1015
00:54:09,120 --> 00:54:10,770
at an intermediate distance.

1016
00:54:10,770 --> 00:54:13,770
And this is really sort
of a general description

1017
00:54:13,770 --> 00:54:15,480
for all kinds of liquids.

1018
00:54:15,480 --> 00:54:18,640
And then you sort of
move farther away.

1019
00:54:18,640 --> 00:54:20,990
And then there is a
second solvation shell.

1020
00:54:20,990 --> 00:54:23,020
OK, there is structure
in the liquid.

1021
00:54:23,020 --> 00:54:24,400
You have a water molecule.

1022
00:54:24,400 --> 00:54:27,250
On average, it's surrounded
by a first shell.

1023
00:54:27,250 --> 00:54:30,730
And then, on average, that's
surrounded by a second shell.

1024
00:54:30,730 --> 00:54:33,010
And as you move
farther and farther,

1025
00:54:33,010 --> 00:54:35,390
the boundaries become
farther and farther.

1026
00:54:35,390 --> 00:54:38,590
So you can sort of spot
a third solvation shell.

1027
00:54:38,590 --> 00:54:42,590
But once you move
farther and farther away,

1028
00:54:42,590 --> 00:54:45,550
if you look at what is the
probability distribution 10

1029
00:54:45,550 --> 00:54:51,220
angstrom, 100 angstrom, 1
meter away, it becomes flat.

1030
00:54:51,220 --> 00:54:54,790
That is, if you look what is the
probability of finding a water

1031
00:54:54,790 --> 00:55:00,280
molecule 1 meter away or 1
meter plus 1/10 of an angstrom,

1032
00:55:00,280 --> 00:55:01,870
there is really no difference.

1033
00:55:01,870 --> 00:55:06,010
Because you have lost
long range order.

1034
00:55:06,010 --> 00:55:07,750
You have lost structure.

1035
00:55:07,750 --> 00:55:11,590
This would be very different,
say, in a perfect crystal.

1036
00:55:11,590 --> 00:55:16,220
So if you study something like
silicon at 0 temperature--

1037
00:55:16,220 --> 00:55:20,060
and for a moment, forget about
zero-point motion and quantum

1038
00:55:20,060 --> 00:55:20,560
effects.

1039
00:55:20,560 --> 00:55:23,110
Think just at a silicon
crystal as being

1040
00:55:23,110 --> 00:55:25,780
made of classical
particles at 0 temperature.

1041
00:55:25,780 --> 00:55:28,280
They are perfectly periodic.

1042
00:55:28,280 --> 00:55:30,610
And so what you would
find, you could still

1043
00:55:30,610 --> 00:55:33,580
define pair
distribution function.

1044
00:55:33,580 --> 00:55:36,750
And what you would
find is Dirac's delta.

1045
00:55:36,750 --> 00:55:40,030
That is the probability
of finding a silicon

1046
00:55:40,030 --> 00:55:43,830
atom at a certain distance
from new silicon atoms,

1047
00:55:43,830 --> 00:55:45,880
say, sitting on
the origin is going

1048
00:55:45,880 --> 00:55:51,400
to be 0 until you hit the
distance of your first shell

1049
00:55:51,400 --> 00:55:52,420
of neighbors.

1050
00:55:52,420 --> 00:55:55,360
And then it's sort of infinite
or properly normalized

1051
00:55:55,360 --> 00:55:56,530
at Dirac's delta.

1052
00:55:56,530 --> 00:55:59,140
And then it's 0 again, and
then another Dirac's delta,

1053
00:55:59,140 --> 00:56:01,600
and 0 again, and then
another Dirac's delta.

1054
00:56:01,600 --> 00:56:05,140
And as you move farther
and farther away,

1055
00:56:05,140 --> 00:56:07,630
the density of
your Dirac's delta

1056
00:56:07,630 --> 00:56:09,460
becomes higher and higher.

1057
00:56:09,460 --> 00:56:13,000
Because, again, if you are
at 1 meter of distance,

1058
00:56:13,000 --> 00:56:15,940
there are an enormous
amount of possibilities

1059
00:56:15,940 --> 00:56:19,660
of finding sort of a
combination of lattice vectors

1060
00:56:19,660 --> 00:56:23,450
that give you a silicon
atom in that position.

1061
00:56:23,450 --> 00:56:28,720
And so if we were to smooth
out those Dirac's delta, that

1062
00:56:28,720 --> 00:56:30,820
is if we were to have
some temperature,

1063
00:56:30,820 --> 00:56:33,700
we would start seeing something
a little bit like this.

1064
00:56:33,700 --> 00:56:37,360
The first Dirac's delta
would become slightly larger.

1065
00:56:37,360 --> 00:56:39,790
The second Dirac's delta,
the third Dirac's delta,

1066
00:56:39,790 --> 00:56:41,830
we move farther
and farther away.

1067
00:56:41,830 --> 00:56:44,170
There are so many
close to each other

1068
00:56:44,170 --> 00:56:46,300
that they all start overlapping.

1069
00:56:46,300 --> 00:56:50,020
And also, for a crystal
at final temperature,

1070
00:56:50,020 --> 00:56:52,210
you would start to have
this pair correlation

1071
00:56:52,210 --> 00:56:54,790
function that goes towards 1.

1072
00:56:54,790 --> 00:56:58,390
In a liquid, it's a
much more dramatic.

1073
00:56:58,390 --> 00:57:02,470
And this quantity tells you
also something fundamental

1074
00:57:02,470 --> 00:57:06,640
about the appropriate size
of your simulation cell.

1075
00:57:06,640 --> 00:57:09,640
Remember that you are studying
something in periodic boundary

1076
00:57:09,640 --> 00:57:10,580
condition.

1077
00:57:10,580 --> 00:57:13,720
So the question is, if we
want to study liquid water,

1078
00:57:13,720 --> 00:57:16,240
how large should
our unit cell be?

1079
00:57:16,240 --> 00:57:18,070
I mean, does it
need to contain 10

1080
00:57:18,070 --> 00:57:20,020
to the 23, the
Avogadro's number,

1081
00:57:20,020 --> 00:57:22,610
of water molecules in
a macroscopic sample?

1082
00:57:22,610 --> 00:57:23,110
No.

1083
00:57:23,110 --> 00:57:26,350
And this is a
fundamental concept.

1084
00:57:26,350 --> 00:57:30,880
It needs to be large enough
so that the correlation

1085
00:57:30,880 --> 00:57:33,400
between your molecule
in the origin

1086
00:57:33,400 --> 00:57:37,870
and that molecule
periodically repeated is 0.

1087
00:57:37,870 --> 00:57:44,410
That is, by the time this pair
correlation function reaches 1,

1088
00:57:44,410 --> 00:57:47,450
it means that you have
lost any structure.

1089
00:57:47,450 --> 00:57:50,860
So somehow you
don't have any way

1090
00:57:50,860 --> 00:57:55,180
of recognizing that there
is a molecule in the origin

1091
00:57:55,180 --> 00:57:59,180
or that molecule is 0.1
angstrom displaced away.

1092
00:57:59,180 --> 00:58:01,660
So a molecule that is on
the origin and a molecule--

1093
00:58:01,660 --> 00:58:03,490
suppose that you
continue this curve

1094
00:58:03,490 --> 00:58:06,190
and maybe becomes
15 angstroms away.

1095
00:58:06,190 --> 00:58:09,010
Well, those two
molecules are not

1096
00:58:09,010 --> 00:58:12,320
aware of each other's existence.

1097
00:58:12,320 --> 00:58:14,260
And so if you are
studying your system,

1098
00:58:14,260 --> 00:58:17,020
such data, the distance
between a molecule

1099
00:58:17,020 --> 00:58:20,650
and your periodic
image is 15 angstrom.

1100
00:58:20,650 --> 00:58:23,290
Actually, from
all point of view,

1101
00:58:23,290 --> 00:58:26,080
you are studying a
system that behave

1102
00:58:26,080 --> 00:58:28,840
as it was an infinite system.

1103
00:58:28,840 --> 00:58:32,290
Your molecule sitting on
the origin moves around.

1104
00:58:32,290 --> 00:58:37,510
And it sees independently
all the molecules with which

1105
00:58:37,510 --> 00:58:39,160
on average it interacts.

1106
00:58:39,160 --> 00:58:41,920
And the molecule that
is 15 angstroms away

1107
00:58:41,920 --> 00:58:45,250
is so far away that
there is no correlation.

1108
00:58:45,250 --> 00:58:47,600
And so for all
practical purposes,

1109
00:58:47,600 --> 00:58:52,270
this molecule on the origin
sees an infinity liquid system.

1110
00:58:52,270 --> 00:58:55,030
It doesn't recognize
that there is a boundary

1111
00:58:55,030 --> 00:58:57,470
and things are
periodically repeated.

1112
00:58:57,470 --> 00:58:59,740
And this is absolutely
very important.

1113
00:58:59,740 --> 00:59:04,570
And that's why we can actually
simulate a realistic liquid

1114
00:59:04,570 --> 00:59:06,640
with actually very
few molecules.

1115
00:59:06,640 --> 00:59:09,880
It turns out that, for
something like water,

1116
00:59:09,880 --> 00:59:15,250
you need only 32 water molecules
at sort of regular conditions

1117
00:59:15,250 --> 00:59:19,120
in the liquid, 300
kelvin, 350 kelvin,

1118
00:59:19,120 --> 00:59:23,050
to actually simulate
an infinite system.

1119
00:59:23,050 --> 00:59:25,780
You don't need to put
more water molecules.

1120
00:59:25,780 --> 00:59:28,330
It won't change anything.

1121
00:59:28,330 --> 00:59:32,500
Of course, as you sort
of get closer and closer

1122
00:59:32,500 --> 00:59:37,000
to a phase transition, remember
that, at the exact temperature

1123
00:59:37,000 --> 00:59:39,730
of the phase transition,
the solid and the liquid

1124
00:59:39,730 --> 00:59:40,880
are in equilibrium.

1125
00:59:40,880 --> 00:59:43,180
So the free energy of the
solid and the free energy

1126
00:59:43,180 --> 00:59:45,220
of the liquid are identical.

1127
00:59:45,220 --> 00:59:48,580
So it doesn't cost anything
in terms of free energy

1128
00:59:48,580 --> 00:59:50,950
to transform the solid in
the liquid and the liquid

1129
00:59:50,950 --> 00:59:52,060
in the solid.

1130
00:59:52,060 --> 00:59:56,920
And that means that you can have
fluctuation, bubbles of solid

1131
00:59:56,920 --> 01:00:02,610
and bubble of liquid,
that have arbitrary sizes.

1132
01:00:02,610 --> 01:00:06,000
So at the exact temperature
of the phase transition,

1133
01:00:06,000 --> 01:00:10,130
you have a divergence in
this characteristic length.

1134
01:00:10,130 --> 01:00:13,590
That is you can have
fluctuations of all sizes.

1135
01:00:13,590 --> 01:00:16,490
So if you were to study
the system, exactly

1136
01:00:16,490 --> 01:00:19,700
the phase transition, you
need an infinite sized system.

1137
01:00:19,700 --> 01:00:22,610
But as you move away from
the phase transition,

1138
01:00:22,610 --> 01:00:26,000
very quickly the correlation
function dropped down.

1139
01:00:26,000 --> 01:00:30,350
And it's something like 25
kelvin above the melting point.

1140
01:00:30,350 --> 01:00:33,800
It's so short that you
need only 32 molecules.

1141
01:00:33,800 --> 01:00:36,650
And water is actually
an excellent example

1142
01:00:36,650 --> 01:00:40,750
to sort of show you also what
goes wrong in the simulation.

1143
01:00:40,750 --> 01:00:41,360
OK.

1144
01:00:41,360 --> 01:00:44,020
Because what I've
plotted here are actually

1145
01:00:44,020 --> 01:00:46,810
a result of accurate
or very accurate

1146
01:00:46,810 --> 01:00:49,190
ab initio molecular
dynamics simulation.

1147
01:00:49,190 --> 01:00:52,190
And we know that sort of
density functional theory

1148
01:00:52,190 --> 01:00:55,250
actually describes
very well the hydrogen

1149
01:00:55,250 --> 01:00:58,250
bonds, so describes the
interaction between water

1150
01:00:58,250 --> 01:01:00,030
molecules very well.

1151
01:01:00,030 --> 01:01:02,300
But you can see it here.

1152
01:01:02,300 --> 01:01:04,850
You see, these are the
pair correlation function.

1153
01:01:04,850 --> 01:01:07,130
You see this sort of
set of curves here,

1154
01:01:07,130 --> 01:01:14,270
the more structural set of
curves, taken at 325, 350, 375

1155
01:01:14,270 --> 01:01:15,260
kelvin.

1156
01:01:15,260 --> 01:01:19,310
There is actually very, very
little difference between them.

1157
01:01:19,310 --> 01:01:20,630
And that's a bad sign.

1158
01:01:20,630 --> 01:01:23,420
So what we are seeing is
that this sort of structure

1159
01:01:23,420 --> 01:01:27,010
is not changing as we
increase the temperature.

1160
01:01:27,010 --> 01:01:28,760
And then we increase
the temperature, say,

1161
01:01:28,760 --> 01:01:31,820
from 375 to 400 kelvin.

1162
01:01:31,820 --> 01:01:33,990
And you see this
different set of curves.

1163
01:01:33,990 --> 01:01:35,480
You see the red and the green.

1164
01:01:35,480 --> 01:01:38,420
So there is really a
discontinuous change

1165
01:01:38,420 --> 01:01:41,720
in going from 375 to 400 kelvin.

1166
01:01:41,720 --> 01:01:43,610
What has happened
in our simulation?

1167
01:01:43,610 --> 01:01:48,500
Well, it turns out that, at 400
kelvin, our water is liquid.

1168
01:01:48,500 --> 01:01:52,520
At 375, it's a glassy system.

1169
01:01:52,520 --> 01:01:54,620
So it's an amorphous glass.

1170
01:01:54,620 --> 01:01:58,010
And you see it doesn't
change with temperature.

1171
01:01:58,010 --> 01:02:02,660
And you need to get to 400
kelvin for it to become liquid.

1172
01:02:02,660 --> 01:02:04,250
What is the bad thing here?

1173
01:02:04,250 --> 01:02:08,090
It's not that the water is
in an amorphous glass state.

1174
01:02:08,090 --> 01:02:11,270
It's reasonable, as I said, even
if it's not the result that we

1175
01:02:11,270 --> 01:02:16,850
want, to find an amorphous
states for a lot of liquids

1176
01:02:16,850 --> 01:02:20,330
that have been cooled down
below the melting temperature

1177
01:02:20,330 --> 01:02:23,630
especially sort of liquids
that tend to have sort

1178
01:02:23,630 --> 01:02:25,700
of pseudo-covalent bonding.

1179
01:02:25,700 --> 01:02:29,840
You study carbon
below its melting pot.

1180
01:02:29,840 --> 01:02:30,750
You study silicon.

1181
01:02:30,750 --> 01:02:32,460
You study something like water.

1182
01:02:32,460 --> 01:02:35,720
They tend to want to be sort
of tetrahedrally coordinated.

1183
01:02:35,720 --> 01:02:40,730
And they tend to sort of form
a lot of glassy structure.

1184
01:02:40,730 --> 01:02:44,930
What is really bad is
that the melting of water

1185
01:02:44,930 --> 01:02:48,620
doesn't seem to be taking
place at 273 kelvin,

1186
01:02:48,620 --> 01:02:51,020
but it takes place
at 400 kelvin.

1187
01:02:51,020 --> 01:02:51,710
OK.

1188
01:02:51,710 --> 01:02:55,250
So if you want, if you were to
use densely functional theory

1189
01:02:55,250 --> 01:02:58,490
and sort of molecular dynamics
to predict what happens,

1190
01:02:58,490 --> 01:03:00,600
we would all be frozen to death.

1191
01:03:00,600 --> 01:03:04,940
You need to be at 400 kelvin
to see something liquid.

1192
01:03:04,940 --> 01:03:07,340
And you see, then
you can also compare.

1193
01:03:07,340 --> 01:03:08,090
This is very good.

1194
01:03:08,090 --> 01:03:10,910
There are a lot of
experiments that

1195
01:03:10,910 --> 01:03:13,160
look at what is the
structure of water.

1196
01:03:13,160 --> 01:03:15,350
You can do this with
neutron scattering.

1197
01:03:15,350 --> 01:03:17,580
You can do this with
X-ray scattering.

1198
01:03:17,580 --> 01:03:19,760
So you can probe
the local structure.

1199
01:03:19,760 --> 01:03:22,430
And you see this is
what the experiment is.

1200
01:03:22,430 --> 01:03:26,510
So it's sort of good, but
it's really not perfect.

1201
01:03:26,510 --> 01:03:29,060
And so what's going on here?

1202
01:03:29,060 --> 01:03:31,400
What's wrong in this?

1203
01:03:31,400 --> 01:03:34,820
Well, in many ways,
we actually think

1204
01:03:34,820 --> 01:03:38,090
that density functional
theory is good enough

1205
01:03:38,090 --> 01:03:39,170
to simulate water.

1206
01:03:39,170 --> 01:03:39,770
OK.

1207
01:03:39,770 --> 01:03:43,310
So the energy for any
given configuration

1208
01:03:43,310 --> 01:03:46,760
is given accurately
by quantum mechanics.

1209
01:03:46,760 --> 01:03:49,450
But there is something
very important

1210
01:03:49,450 --> 01:03:52,490
that I sort of hinted in
some of the discussion.

1211
01:03:52,490 --> 01:03:56,330
300 kelvin, 400 kelvin
is a temperature

1212
01:03:56,330 --> 01:04:02,270
at which molecules can still
display quantum effects

1213
01:04:02,270 --> 01:04:05,120
in their vibrational
degrees of freedom.

1214
01:04:05,120 --> 01:04:08,300
That is when we discussed
at certain point

1215
01:04:08,300 --> 01:04:12,710
the fact that a simple
molecule like hydrogen

1216
01:04:12,710 --> 01:04:17,780
is actually only discrete
vibrational energy

1217
01:04:17,780 --> 01:04:19,880
and discrete vibrational states.

1218
01:04:19,880 --> 01:04:22,280
And at 0 temperature,
it's in the state

1219
01:04:22,280 --> 01:04:24,110
called zero-point motion.

1220
01:04:24,110 --> 01:04:27,380
That is a purely quantum
mechanical effect

1221
01:04:27,380 --> 01:04:32,060
and doesn't really come
from describing correctly

1222
01:04:32,060 --> 01:04:36,020
the energy of the electrons
in the potential of the ions,

1223
01:04:36,020 --> 01:04:38,990
but comes from the fact
that the ions themselves

1224
01:04:38,990 --> 01:04:40,640
are quantum particles.

1225
01:04:40,640 --> 01:04:43,520
And the excitation
of a quantum particle

1226
01:04:43,520 --> 01:04:46,760
from one state to
another vibrational state

1227
01:04:46,760 --> 01:04:49,020
is in itself quantized.

1228
01:04:49,020 --> 01:04:51,650
And so what turns out
to happen in water

1229
01:04:51,650 --> 01:04:57,000
is that the hydrogen oxygen
stretching frequency.

1230
01:04:57,000 --> 01:04:57,770
OK.

1231
01:04:57,770 --> 01:05:01,040
So the frequencies by
which one hydrogen beats

1232
01:05:01,040 --> 01:05:04,280
against another oxygen
are actually very high,

1233
01:05:04,280 --> 01:05:08,210
are of the order of 3,000
centimeters to the minus 1.

1234
01:05:08,210 --> 01:05:10,280
So that's a typical
sort of measure

1235
01:05:10,280 --> 01:05:11,720
of vibrational frequency.

1236
01:05:11,720 --> 01:05:15,050
They are, in particular,
sort of 10 times larger

1237
01:05:15,050 --> 01:05:17,210
than the average
kinetic temperature

1238
01:05:17,210 --> 01:05:18,930
at room temperature.

1239
01:05:18,930 --> 01:05:24,740
So that means that, in reality,
those vibrational modes cannot

1240
01:05:24,740 --> 01:05:25,760
be excited.

1241
01:05:25,760 --> 01:05:28,830
They are frozen in
their ground state.

1242
01:05:28,830 --> 01:05:32,780
So in all this liquid, you have
a variety of vibrational modes,

1243
01:05:32,780 --> 01:05:35,330
some that are very
soft and that involve

1244
01:05:35,330 --> 01:05:39,470
sort of how one water molecule
correlates with other water

1245
01:05:39,470 --> 01:05:40,280
molecules.

1246
01:05:40,280 --> 01:05:43,760
But you have some very
stiff intermolecular

1247
01:05:43,760 --> 01:05:48,050
vibrational degrees of
freedom that you can't excite.

1248
01:05:48,050 --> 01:05:52,890
And so the classical
description of the ions is good,

1249
01:05:52,890 --> 01:05:54,020
but it's not perfect.

1250
01:05:54,020 --> 01:05:57,140
And that's what brings
us this error here.

1251
01:05:57,140 --> 01:06:00,230
If we were to study-- sadly,
there are no experiments, say,

1252
01:06:00,230 --> 01:06:01,820
at 3,000 kelvin.

1253
01:06:01,820 --> 01:06:06,230
But if you were to do an
experiment of, obviously,

1254
01:06:06,230 --> 01:06:09,320
water at a sort of
fixed volume, so

1255
01:06:09,320 --> 01:06:13,070
at very high pressure
at 3,000 kelvin, then

1256
01:06:13,070 --> 01:06:17,420
the ions would all sort of
behave as classical particles

1257
01:06:17,420 --> 01:06:20,370
and probably our agreements
would be much better.

1258
01:06:20,370 --> 01:06:22,520
But at this low
temperature, some

1259
01:06:22,520 --> 01:06:26,130
of these vibrational degrees
of freedom are actually frozen.

1260
01:06:26,130 --> 01:06:29,390
And they can't
exchange temperature.

1261
01:06:29,390 --> 01:06:31,580
They can't exchange
vibrational energy

1262
01:06:31,580 --> 01:06:34,400
with the rest of the
degrees of freedom.

1263
01:06:34,400 --> 01:06:39,800
It's the same concept by which
the specific heat of solids

1264
01:06:39,800 --> 01:06:42,140
follows the Dulong
and Petit law.

1265
01:06:42,140 --> 01:06:45,740
And it's much smaller than
its ideal sort of value

1266
01:06:45,740 --> 01:06:46,820
at low temperatures.

1267
01:06:46,820 --> 01:06:50,570
Because even in a solid,
where vibrational frequencies

1268
01:06:50,570 --> 01:06:54,140
tend to be much lower than 3,000
centimeters to the minus 1,

1269
01:06:54,140 --> 01:06:58,190
it tends to be of the order 300,
400 centimeters to the minus 1.

1270
01:06:58,190 --> 01:07:01,400
So they tend to be all
populated at room temperature.

1271
01:07:01,400 --> 01:07:03,900
But in the regime
of low temperatures,

1272
01:07:03,900 --> 01:07:07,790
much lower 300 kelvin, some
of those degrees of freedom

1273
01:07:07,790 --> 01:07:09,020
are also frozen out.

1274
01:07:09,020 --> 01:07:11,450
So they can't
exchange temperature

1275
01:07:11,450 --> 01:07:14,150
when you increase the heat.

1276
01:07:14,150 --> 01:07:16,610
And so the specific
heat of the system

1277
01:07:16,610 --> 01:07:20,090
is much smaller because a lot
of those degrees of freedom

1278
01:07:20,090 --> 01:07:24,450
can't actually talk
with your heat bath.

1279
01:07:24,450 --> 01:07:24,950
OK.

1280
01:07:24,950 --> 01:07:26,480
So this was an example.

1281
01:07:26,480 --> 01:07:28,160
Let's see another
example in which

1282
01:07:28,160 --> 01:07:32,750
we look at ensemble averages.

1283
01:07:32,750 --> 01:07:35,490
And this is the example
of the caloric curve.

1284
01:07:35,490 --> 01:07:39,590
So that is plotting
the potential energy

1285
01:07:39,590 --> 01:07:41,430
as a function of time.

1286
01:07:41,430 --> 01:07:45,210
This is done for a sort of
classical molecular dynamic

1287
01:07:45,210 --> 01:07:46,350
simulation.

1288
01:07:46,350 --> 01:07:53,060
So you see it's using
6,400 particles, so

1289
01:07:53,060 --> 01:07:54,710
a very large amount.

1290
01:07:54,710 --> 01:08:00,680
This is actually a slab of
lead in the 1, 0, 0 direction.

1291
01:08:00,680 --> 01:08:04,040
And you see, if we simulate
this at 600 kelvin,

1292
01:08:04,040 --> 01:08:06,120
well, this is as a
function of time.

1293
01:08:06,120 --> 01:08:08,000
This is a very long
simulation, you see.

1294
01:08:08,000 --> 01:08:10,820
We are talking about
the nanoseconds.

1295
01:08:10,820 --> 01:08:11,600
OK.

1296
01:08:11,600 --> 01:08:17,370
So it's of the order
of 1 million time steps

1297
01:08:17,370 --> 01:08:19,850
if we were using times
step for a femtosecond.

1298
01:08:19,850 --> 01:08:21,979
And you see, this
is how the sort

1299
01:08:21,979 --> 01:08:29,930
of average potential energy
fluctuates at 607 kelvin.

1300
01:08:29,930 --> 01:08:32,390
And then we increase
the temperature,

1301
01:08:32,390 --> 01:08:35,930
and you see it fluctuates
around another value.

1302
01:08:35,930 --> 01:08:38,840
Increase the temperature,
it fluctuates another value.

1303
01:08:38,840 --> 01:08:44,720
We increase the temperature a
little bit more from 625 kelvin

1304
01:08:44,720 --> 01:08:49,439
to 630 kelvin, and what
happens in this solid?

1305
01:08:49,439 --> 01:08:52,550
You see, at the beginning,
during, if you want,

1306
01:08:52,550 --> 01:08:55,850
your thermalization
time, you oscillate

1307
01:08:55,850 --> 01:08:57,680
around a certain value.

1308
01:08:57,680 --> 01:09:01,370
But then something
happens, and your system

1309
01:09:01,370 --> 01:09:05,899
goes into a new state in which
the average potential energy

1310
01:09:05,899 --> 01:09:08,479
oscillates around a new value.

1311
01:09:08,479 --> 01:09:10,580
What your system
has done, it's gone

1312
01:09:10,580 --> 01:09:12,740
through a melting transition.

1313
01:09:12,740 --> 01:09:15,410
And so the difference
between, if you want,

1314
01:09:15,410 --> 01:09:19,580
the potential energy just
below the melting temperature

1315
01:09:19,580 --> 01:09:24,080
and the potential energy just
above the melting temperature

1316
01:09:24,080 --> 01:09:29,930
is this latent heat of
fusion, the energy that

1317
01:09:29,930 --> 01:09:34,250
is needed to melt the crystal
or the energy that is released

1318
01:09:34,250 --> 01:09:36,560
when the crystal solidifies.

1319
01:09:36,560 --> 01:09:39,620
And this is actually a sort
of very specific system

1320
01:09:39,620 --> 01:09:43,910
in which there is a sort of
more exotic phase transition.

1321
01:09:43,910 --> 01:09:46,310
And that's why molecular
dynamics simulations

1322
01:09:46,310 --> 01:09:47,090
were needed.

1323
01:09:47,090 --> 01:09:49,220
It was a system
in which there was

1324
01:09:49,220 --> 01:09:52,340
what is called incomplete
surface melting.

1325
01:09:52,340 --> 01:09:55,310
That is, if you look at
the surface of lead 1,

1326
01:09:55,310 --> 01:09:59,570
0, 0, it starts melting
at a temperature that

1327
01:09:59,570 --> 01:10:03,830
is lower than the temperature
at which the bulk melts.

1328
01:10:03,830 --> 01:10:06,380
And this is because
surface atoms are really

1329
01:10:06,380 --> 01:10:08,360
different from bulk atoms.

1330
01:10:08,360 --> 01:10:10,920
So they have their
own phase transition.

1331
01:10:10,920 --> 01:10:13,400
And it's temperature for
that phase transition

1332
01:10:13,400 --> 01:10:16,140
that is different from the
bulk melting temperature.

1333
01:10:16,140 --> 01:10:18,320
So you heat up lead.

1334
01:10:18,320 --> 01:10:22,370
All of a sudden, you see that
the surface becomes liquid.

1335
01:10:22,370 --> 01:10:28,430
And then you need to reach
another critical temperature

1336
01:10:28,430 --> 01:10:31,140
before also the bulk melts.

1337
01:10:31,140 --> 01:10:31,640
OK.

1338
01:10:31,640 --> 01:10:35,420
So a couple of examples we have
seen from ensemble averages.

1339
01:10:35,420 --> 01:10:38,810
As I said, we can also
use molecular dynamics

1340
01:10:38,810 --> 01:10:42,210
to actually follow
a chemical reaction.

1341
01:10:42,210 --> 01:10:45,440
And so this is another
example from ab

1342
01:10:45,440 --> 01:10:47,600
initio molecular
dynamic simulation.

1343
01:10:47,600 --> 01:10:50,150
Because what we really
want to investigate

1344
01:10:50,150 --> 01:10:55,460
is how we can take a molecule
of methane-- you see C and H4.

1345
01:10:55,460 --> 01:10:57,440
There are four hydrogen atom.

1346
01:10:57,440 --> 01:11:01,670
One is sort of the pink and
below, beneath the yellow.

1347
01:11:01,670 --> 01:11:04,550
We want to figure out
how to break this apart.

1348
01:11:04,550 --> 01:11:07,310
This is a fundamental
technological problem

1349
01:11:07,310 --> 01:11:11,060
because we have enormous
reserves of methane is

1350
01:11:11,060 --> 01:11:13,500
what people call natural gas.

1351
01:11:13,500 --> 01:11:17,390
But it's extremely expensive
to transport because you really

1352
01:11:17,390 --> 01:11:22,250
need to liquefy it and carry
it in sort of cryogenically

1353
01:11:22,250 --> 01:11:24,170
frozen containers.

1354
01:11:24,170 --> 01:11:26,760
So it's much more expensive
to transport methane

1355
01:11:26,760 --> 01:11:27,980
than to transport oil.

1356
01:11:27,980 --> 01:11:30,500
And that's why we are
actually not using

1357
01:11:30,500 --> 01:11:35,100
sort of all the methane natural
gas reserves that exist.

1358
01:11:35,100 --> 01:11:37,670
So one of the important
things that we could do

1359
01:11:37,670 --> 01:11:44,150
is actually transform methane
in some higher and heavier

1360
01:11:44,150 --> 01:11:46,430
hydrocarbons, like
methanol, that we

1361
01:11:46,430 --> 01:11:48,800
can transport much more easily.

1362
01:11:48,800 --> 01:11:50,780
And that, in principle,
is very doable.

1363
01:11:50,780 --> 01:11:53,600
Methane has an enormous
amount of energies

1364
01:11:53,600 --> 01:11:54,950
ready to be released.

1365
01:11:54,950 --> 01:11:58,160
But sadly, in order to
release the energy, in order

1366
01:11:58,160 --> 01:12:01,910
if you want to burn methane
or in order to break it apart,

1367
01:12:01,910 --> 01:12:06,770
you need first to break apart
the first carbon-hydrogen bond.

1368
01:12:06,770 --> 01:12:08,100
And that's tough.

1369
01:12:08,100 --> 01:12:10,520
Once you break
that hydrogen bond,

1370
01:12:10,520 --> 01:12:13,760
all the other hydrogen go
and release a lot of energy.

1371
01:12:13,760 --> 01:12:16,160
But breaking this
is very difficult.

1372
01:12:16,160 --> 01:12:19,460
And we want to do this in
a very controlled manner

1373
01:12:19,460 --> 01:12:24,290
in order not to burn methane,
but construct some higher

1374
01:12:24,290 --> 01:12:26,760
hydrocarbon, like methanol.

1375
01:12:26,760 --> 01:12:29,960
And so one sort of very
suggestive hypotheses

1376
01:12:29,960 --> 01:12:33,690
is actually using
nanostructured catalyst, that

1377
01:12:33,690 --> 01:12:36,900
is using well-known
catalyst transition metals,

1378
01:12:36,900 --> 01:12:40,350
like platinum in this
case, but catalysts

1379
01:12:40,350 --> 01:12:43,590
that have a very exotic
surface, so that they

1380
01:12:43,590 --> 01:12:44,820
become more reactive.

1381
01:12:44,820 --> 01:12:49,090
They can help break part of
this carbon-hydrogen bond.

1382
01:12:49,090 --> 01:12:52,770
And so this is where ab
initio molecular dynamics come

1383
01:12:52,770 --> 01:12:56,040
into play because we
can study the system

1384
01:12:56,040 --> 01:13:06,430
and try to figure out what
happens when we take a methane

1385
01:13:06,430 --> 01:13:11,000
molecule, and we throw
it to a platinum surface.

1386
01:13:11,000 --> 01:13:14,230
And as you can see,
nothing happens.

1387
01:13:14,230 --> 01:13:16,750
And we are very, very
happy, not so much

1388
01:13:16,750 --> 01:13:18,970
because we have sort of
simulated an interesting

1389
01:13:18,970 --> 01:13:23,230
system, but because actually we
can measure the kinetic energy

1390
01:13:23,230 --> 01:13:25,060
of the molecules coming down.

1391
01:13:25,060 --> 01:13:27,400
And we can measure
when it comes out.

1392
01:13:27,400 --> 01:13:30,970
And we find out that it
has the same velocity,

1393
01:13:30,970 --> 01:13:32,440
but opposite direction.

1394
01:13:32,440 --> 01:13:36,550
So our system is conserving
energy very accurately.

1395
01:13:36,550 --> 01:13:38,500
But you know, what
you can do now

1396
01:13:38,500 --> 01:13:40,780
with your molecular
dynamic simulation

1397
01:13:40,780 --> 01:13:43,930
is try to figure out if
this is the whole story.

1398
01:13:43,930 --> 01:13:46,300
That is, well, as you
see, we have sort of

1399
01:13:46,300 --> 01:13:49,480
thrown this molecule
into the trough here.

1400
01:13:49,480 --> 01:13:52,990
Because we thought
that maybe this area

1401
01:13:52,990 --> 01:13:54,850
is the most
catalytically active.

1402
01:13:54,850 --> 01:13:58,420
It's where there are more
platinum atoms able to sort

1403
01:13:58,420 --> 01:13:59,950
of rip this thing apart.

1404
01:13:59,950 --> 01:14:01,690
Well, it turns out
that, actually,

1405
01:14:01,690 --> 01:14:05,920
the most catalytically
active part of the surface

1406
01:14:05,920 --> 01:14:09,550
is the under-coordinated
platinum atoms sitting here

1407
01:14:09,550 --> 01:14:10,510
on the ridge.

1408
01:14:10,510 --> 01:14:15,400
And you see, when the
molecule falls down on this,

1409
01:14:15,400 --> 01:14:18,760
it very easily
breaks apart this.

1410
01:14:18,760 --> 01:14:22,390
And so doing our
ab initio MD, we

1411
01:14:22,390 --> 01:14:26,380
discover that these
under-coordinated atoms

1412
01:14:26,380 --> 01:14:29,020
are actually much more
active than the sort

1413
01:14:29,020 --> 01:14:32,020
of surface-like atoms down here.

1414
01:14:32,020 --> 01:14:34,900
So under-coordinating
a transition metal

1415
01:14:34,900 --> 01:14:36,730
is actually very helpful.

1416
01:14:36,730 --> 01:14:40,180
And this is actually
what bacteria that

1417
01:14:40,180 --> 01:14:42,790
eat methane for a living do.

1418
01:14:42,790 --> 01:14:46,690
So what those bacterias
do, they have an enzyme

1419
01:14:46,690 --> 01:14:49,720
that is called a methane
monooxygense that

1420
01:14:49,720 --> 01:14:54,130
actually is able to break apart
this carbon-hydrogen bond.

1421
01:14:54,130 --> 01:14:57,370
And what it has, it
has sort of two iron

1422
01:14:57,370 --> 01:15:01,810
atoms in the middle of a sort
of porphyrin like structure.

1423
01:15:01,810 --> 01:15:04,930
And those under-coordinated
transition metals

1424
01:15:04,930 --> 01:15:06,940
do all the trick.

1425
01:15:06,940 --> 01:15:11,650
Other examples in which a
molecular dynamic simulation

1426
01:15:11,650 --> 01:15:14,770
helps us understanding
the structure

1427
01:15:14,770 --> 01:15:18,320
of microscopic system--

1428
01:15:18,320 --> 01:15:19,370
as shown here.

1429
01:15:19,370 --> 01:15:23,080
In this case, what we
are doing is taking

1430
01:15:23,080 --> 01:15:26,330
a liquid droplet of aluminum.

1431
01:15:26,330 --> 01:15:29,980
So this is aluminum above
its melting temperature.

1432
01:15:29,980 --> 01:15:33,340
That is sort of 933 kelvin.

1433
01:15:33,340 --> 01:15:36,910
It's been thermalized,
and so it's a liquid.

1434
01:15:36,910 --> 01:15:39,160
And then we have a
surface that could

1435
01:15:39,160 --> 01:15:42,340
be, say, two different
crystallographic orientation,

1436
01:15:42,340 --> 01:15:46,860
the 1, 1, 0 surface and
the 1, 1, 1 surface.

1437
01:15:46,860 --> 01:15:51,610
And what we do is we drop this
liquid droplet on the surface.

1438
01:15:51,610 --> 01:15:55,390
And we sort of see what is
the evolution of the system.

1439
01:15:55,390 --> 01:16:00,340
And on the 1, 1, 0 surface, that
is one of these metal surfaces

1440
01:16:00,340 --> 01:16:03,700
that undergo a pre-melting
transition where the surface

1441
01:16:03,700 --> 01:16:04,810
itself--

1442
01:16:04,810 --> 01:16:06,190
you can see it here--

1443
01:16:06,190 --> 01:16:07,810
becomes liquid-like.

1444
01:16:07,810 --> 01:16:10,210
What happens is that you
have a liquid droplet falling

1445
01:16:10,210 --> 01:16:12,350
on a sort of liquid surface.

1446
01:16:12,350 --> 01:16:17,080
And so it wets completely the
surface, and this moves around.

1447
01:16:17,080 --> 01:16:21,760
You throw this same liquid
droplet on aluminum 1, 1, 1,

1448
01:16:21,760 --> 01:16:24,130
that is a surface
that doesn't matter

1449
01:16:24,130 --> 01:16:26,260
and that actually overheat.

1450
01:16:26,260 --> 01:16:30,250
You can sort of bring a
crystalite of aluminum,

1451
01:16:30,250 --> 01:16:32,440
and you can look at
the 1, 1, 1 surface.

1452
01:16:32,440 --> 01:16:35,140
And you can bring it above
the melting temperature,

1453
01:16:35,140 --> 01:16:37,600
and it will not
melt until you reach

1454
01:16:37,600 --> 01:16:40,510
sort of a critical
temperature of 1,000 kelvin.

1455
01:16:40,510 --> 01:16:44,560
Well, this liquid droplet
actually doesn't wet away,

1456
01:16:44,560 --> 01:16:47,170
but it forms, you
see, a little disk.

1457
01:16:47,170 --> 01:16:49,150
And you can actually
sort of measure

1458
01:16:49,150 --> 01:16:51,070
what is this contact angle.

1459
01:16:51,070 --> 01:16:53,950
And so you can figure out
what is the surface tension.

1460
01:16:53,950 --> 01:16:56,290
And you can understand
a lot basically

1461
01:16:56,290 --> 01:17:01,360
about the free energy of the
liquid and the solid just

1462
01:17:01,360 --> 01:17:04,240
by those contact angles.

1463
01:17:04,240 --> 01:17:08,090
Another example for which
I have another movie

1464
01:17:08,090 --> 01:17:10,360
is sort of related
to what you have seen

1465
01:17:10,360 --> 01:17:12,280
in the sort of previous slide.

1466
01:17:12,280 --> 01:17:14,350
And it's sort of
looking at actually

1467
01:17:14,350 --> 01:17:20,470
what do adatom on aluminum
1, 0, 0 surface do.

1468
01:17:20,470 --> 01:17:24,100
And so you see what we do is,
in ab initio molecular dynamics,

1469
01:17:24,100 --> 01:17:26,170
we put an atom here.

1470
01:17:26,170 --> 01:17:29,980
And we let the system evolve
at a certain temperature.

1471
01:17:29,980 --> 01:17:33,310
In this case, it's
evolving at 500 kelvin.

1472
01:17:33,310 --> 01:17:36,890
That is roughly 1/2 of
the melting temperature.

1473
01:17:36,890 --> 01:17:39,460
And again, you can
see some features

1474
01:17:39,460 --> 01:17:42,610
of an ab initio molecular
dynamic calculation in which we

1475
01:17:42,610 --> 01:17:44,740
can simulate very few atoms.

1476
01:17:44,740 --> 01:17:46,240
And so actually,
in this simulation,

1477
01:17:46,240 --> 01:17:51,610
you see our surface area is
given just by 16 atoms, this 4

1478
01:17:51,610 --> 01:17:52,690
by 4 here.

1479
01:17:52,690 --> 01:17:55,340
And then we have periodic
boundary condition.

1480
01:17:55,340 --> 01:17:59,830
So this atom and this atom are
doing exactly the same thing,

1481
01:17:59,830 --> 01:18:02,440
but we are actually
ourselves in the condition

1482
01:18:02,440 --> 01:18:08,390
in which there is no correlation
between what this atom does

1483
01:18:08,390 --> 01:18:10,010
and what this atom does.

1484
01:18:10,010 --> 01:18:13,580
So these atoms sort of
behave independently

1485
01:18:13,580 --> 01:18:16,040
from what this atom does.

1486
01:18:16,040 --> 01:18:17,870
Just because, at
this temperature,

1487
01:18:17,870 --> 01:18:22,070
the correlation length between
these two is very, very short.

1488
01:18:22,070 --> 01:18:24,710
And what I'm not showing
here in this simulation

1489
01:18:24,710 --> 01:18:29,120
is actually the lower
layers of out aluminum slab

1490
01:18:29,120 --> 01:18:32,000
just because they are not
doing anything interesting.

1491
01:18:32,000 --> 01:18:36,830
And so we want to understand
how the atoms move around.

1492
01:18:36,830 --> 01:18:38,580
And everyone would
have sort of guessed

1493
01:18:38,580 --> 01:18:41,930
that adatoms on the
surface hop around,

1494
01:18:41,930 --> 01:18:44,390
sort of the classical
model of diffusion,

1495
01:18:44,390 --> 01:18:47,690
from one side here in the
middle to another side

1496
01:18:47,690 --> 01:18:48,710
here in the middle.

1497
01:18:48,710 --> 01:18:50,240
But you see, you
do your ab initio

1498
01:18:50,240 --> 01:18:51,950
molecular dynamic
simulation, and you

1499
01:18:51,950 --> 01:18:54,470
discover that something
very intriguing happens.

1500
01:18:54,470 --> 01:18:58,070
Actually, adatoms do
not diffuse by hopping,

1501
01:18:58,070 --> 01:19:00,800
but they diffuse
by sort of diving

1502
01:19:00,800 --> 01:19:04,865
into the surface layer and sort
of pushing another atom out.

1503
01:19:04,865 --> 01:19:06,740
And you'll see it again
in a cleaner fashion.

1504
01:19:06,740 --> 01:19:10,820
You see the diffusion takes
place by an atom diving in.

1505
01:19:10,820 --> 01:19:11,450
Look at this.

1506
01:19:11,450 --> 01:19:14,750
This sort of, in a
moment, will dive in

1507
01:19:14,750 --> 01:19:19,160
and pushes this sort
of second atom out.

1508
01:19:19,160 --> 01:19:21,800
And by the time it
reaches here, now you

1509
01:19:21,800 --> 01:19:25,560
see it appearing in
periodic boundary condition.

1510
01:19:25,560 --> 01:19:28,820
And again, the molecular
dynamic simulation

1511
01:19:28,820 --> 01:19:32,630
has given you an insight on
the microscopic process that

1512
01:19:32,630 --> 01:19:35,150
provided density functional
theory is accurate.

1513
01:19:35,150 --> 01:19:37,460
And it, for this
particular system,

1514
01:19:37,460 --> 01:19:42,950
is actually very different from
any kind of physical intuition

1515
01:19:42,950 --> 01:19:44,890
would have told you.

1516
01:19:44,890 --> 01:19:46,790
And it's 1 more minute.

1517
01:19:46,790 --> 01:19:49,700
Let me finish this movie because
something more interesting

1518
01:19:49,700 --> 01:19:54,630
takes place once you increase
the temperature to 700 kelvin.

1519
01:19:54,630 --> 01:19:56,840
This is what we have done
in the simulation here.

1520
01:19:56,840 --> 01:20:00,350
And what you actually see
is that the first surface

1521
01:20:00,350 --> 01:20:04,080
layer starts really
behaving as a liquid.

1522
01:20:04,080 --> 01:20:04,760
OK.

1523
01:20:04,760 --> 01:20:07,280
So we actually have
a phase transition

1524
01:20:07,280 --> 01:20:10,340
that is induced by the
presence of that atom in which

1525
01:20:10,340 --> 01:20:14,240
this first layer, if you
monitor long enough what

1526
01:20:14,240 --> 01:20:15,770
the mean square
displacements are,

1527
01:20:15,770 --> 01:20:18,890
how the atoms are, well,
the atoms in this layer

1528
01:20:18,890 --> 01:20:21,200
move around like a liquid.

1529
01:20:21,200 --> 01:20:25,100
But they don't move
around below the surface.

1530
01:20:25,100 --> 01:20:28,190
And the surface
remain crystalline.

1531
01:20:28,190 --> 01:20:30,410
And I'll finish the movie.

1532
01:20:30,410 --> 01:20:33,290
We actually then, in our
computational experiment,

1533
01:20:33,290 --> 01:20:37,520
all of a sudden we can remove
atoms from our liquid sample.

1534
01:20:37,520 --> 01:20:40,100
That is we can sort
of take them away.

1535
01:20:40,100 --> 01:20:44,330
And we actually see very quickly
that the system crystallizes

1536
01:20:44,330 --> 01:20:45,050
again.

1537
01:20:45,050 --> 01:20:46,920
And you'll see it in a moment.

1538
01:20:46,920 --> 01:20:48,050
See, we have removed atoms.

1539
01:20:48,050 --> 01:20:49,600
We have actually
removed two atoms.

1540
01:20:49,600 --> 01:20:52,040
So we have created a vacancy.

1541
01:20:52,040 --> 01:20:54,020
And the system has
become crystalline.

1542
01:20:54,020 --> 01:20:57,590
Now, the atoms vibrate again
around the crystalline surface.

1543
01:20:57,590 --> 01:20:59,390
And there is a vacancy here.

1544
01:20:59,390 --> 01:21:02,600
And from this, we also
discover how the vacancy moves.

1545
01:21:02,600 --> 01:21:07,310
And you'll see that there is an
atom from below that hops out.

1546
01:21:07,310 --> 01:21:10,460
You see this hops out,
leaving a hole here.

1547
01:21:10,460 --> 01:21:13,520
And an atom from the
first layer pops down.

1548
01:21:13,520 --> 01:21:17,400
So the vacancy that was
here has moved here.

1549
01:21:17,400 --> 01:21:20,330
And again, the way
it has moved is not

1550
01:21:20,330 --> 01:21:23,000
sort of the usual way in
which one could imagine

1551
01:21:23,000 --> 01:21:27,140
a vacancy moving by a
surface atom moving around,

1552
01:21:27,140 --> 01:21:31,460
but it actually involves a
second layer atom jumping out.

1553
01:21:31,460 --> 01:21:31,970
OK.

1554
01:21:31,970 --> 01:21:33,890
With this movie,
I'll finish today's

1555
01:21:33,890 --> 01:21:35,660
molecular dynamic simulation.

1556
01:21:35,660 --> 01:21:39,780
We'll sort conclude the
MD part on Thursday.

1557
01:21:39,780 --> 01:21:41,780
And as a reminder,
today was also

1558
01:21:41,780 --> 01:21:45,770
the deadline for your
second assignment, lab two.

1559
01:21:45,770 --> 01:21:48,950
And you can give it to
me or to one of the TAs.

1560
01:21:48,950 --> 01:21:51,790
And otherwise, we'll
meet on Thursday.