Lecture 13: Molecular Dynamics I

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

Topics covered: Molecular Dynamics I

Instructor: Prof. Nicola Marzari

PROFESSOR: So what we'll do today is really switching to the computational description and characterization of thermodynamic properties. As we have seen in a lot of the previous lectures that we have been focusing on finding out good energy models, that is being able to calculate, what is the energy of a system given the coordinates of its ions?

And now what we want to understand, really, is, what is the evolution of the system in the presence of a finite temperature? Temperature really means that there is an average kinetic energy available for all the atoms and the molecules in your system. That means that there is really a microscopic dynamics. And we want to follow that.

One of the most important things-- and this goes back to, really, the first lecture that you have seen in this class for the quantum part-- is that there is the fundamental de Broglie relation that gives us an estimate of what is the quantum wavelength of any kind of object. Remember, we have said that the wavelength times the momentum should be equal to the Planck constant.

And that really means that for an object like an electron, the de Broglie wavelength that comes out from that relation is comparable to the distance between the atoms. So if the wavelength of an electron is of the order of the angstrom, 10 to the minus 10 meters, it means that the wavelike properties of the electrons needs to be studied. And the electron will behave as a wave. And that's why we use the Schrodinger equation to describe electrons.

The wavelength of the nuclei is much, much smaller. That means that nuclei really didn't-- we don't really need to take into account to a good approximation the quantum nature of the nuclei. The wavelength is so small that they can be treated practically, correctly, as classical particles.

And that's why, actually, nuclei, in a way, still follow the rules of classical dynamics. That is, they evolve in time following Newton's equation of motion. Again, this is an approximation that involves treating the nuclei as classical particles and not as quantum particles. But especially at high temperature, this is a very good approximation.

And so a lot of approaches in computational and material science have to do with solving Newton's equation of motion for a system of interacting nuclei and evolving those nuclei in time to study the thermodynamic properties of the system. And so that's why I actually started with a reminder of Newton's equation of motion. That is really an ordinary second-order differential equation.

That is, when you have a particle of mass m in a force field that is this F of r, well, we can calculate its trajectory-- that is, its position-- as a function of time and its velocity as a function of time by integrating this differential equation. And in the classical case that is always presented, say, is the evolution of a particle in a constant force field, like a particle under the gravitational law, that, to a large extent, is constant if you don't move very much.

And because, in that case, the force is really a constant, the integration of this becomes trivial. It's really a second derivative that needs to be equal a constant. And integrating that gives us a parabola.

But the reason why I'm focusing on this is to remind you of the initial conditions in our differential equations. That is, if we have a second-order differential equation, the general solution-- let's say in this case for this force field would be a parabola-- is uniquely identified once we specify the position at a certain instant and the velocity at a certain instant.

That is, among all the possible parabolas in the world, you selected the one that your projectile is going to follow once you have the position and the velocity at a certain instant. And that's why actually dynamics and molecular dynamics and other physical properties basically depend essentially on the position of the particle in your system and the velocity of the particle in your system.

That's why there is no physical observable that depends, say, on the fifth derivative with respect to time. It's really because the dynamics of that system follows Newton's equation of motion.

OK, with this, what we really need to deal in general is with a system of interacting particles. And for a moment you could actually think that this i, interacting particles, could be the planets orbiting around the sun. So each of the planet i will have its own trajectory. Our r of t will have its own mass.

And they all live in a force field. That's really where the many-body complexity, again, comes out. So many-body complexity, that is much easier to deal with than the many-body complexity of the electron problem because, again, now we describe each object just with a set of six variables, so its 3-space coordinate and its 3-momentum coordinate.

But still, each and every object interacts with each other. When you look at the planets orbiting around the sun, you, at first approximation, could only consider the gravitational attraction to the sun. But in reality, also the planets interact with each other.

And so, again, even if this is much simpler than the Schrodinger equation, as soon as you start to have more than two bodies interacting, it really becomes too complex to solve analytically. And so all the molecular dynamic techniques have, as their only object, basically, the integration in time of this equation of motion in an efficient form.

There are a couple of comments that one wants to make. In general, we deal with conservative fields. That is, the force depends on position only. That also means that the work that you make in going from one point to another point, the integral of the force times the displacement is going to be independent from the trajectory. So the energy, if you want, is just a function of the position.

The other thing that descends immediately from this equation of motion is that the total energy of the system is conserved. And that's very easy to prove. And I'll do it in a moment.

And that's what's often called a microcanonical evolution. When you read about microcanonical ensembles, it means that you are actually under the condition in which the total energy of your system is conserved. And again, what is the total energy of your system? Well, if you have one particle, it's just, trivially, the kinetic energy 1/2 mv squared plus the potential energy V of r of that particle, and just doing everything in one dimension.

And you ask, why is this going to be conserved? Well, you can just prove it by taking the time derivative of this quantity, OK? So we have taken kinetic energy plus potential energy. We look at the time derivative. And this is just simple algebra.

What you get is 1/2 m times 2 times v times the derivative with respect to time, that is the acceleration, plus the derivative of the potential with respect to time. And this can just be written as mass times velocity times acceleration plus-- and we can rewrite the derivative of the potential with time as the derivative of the potential with respect to position times the derivative of the position with respect to time.

And now you see where we are getting to, because this term here is equal to minus the force, while these terms here is equal to the velocity. And so it goes away with here. And so what we have is that m, the derivative of the total energy with respect to time, is equal to the mass times acceleration minus the force. That is 0 because the particle follows Newton's equation of motion.

So because of Newton's equation of motion and because the particle follows exactly that evolution, we have that the sum of the kinetic energy plus the potential energy needs to be conserved during the evolution. And that's actually the fundamental sanity check that you will always have and you can always do during a molecular dynamics simulation. That is, you can check that your integration is correct by making sure that the total energy of your system is concerned.

OK, so what do we do in a practical molecular dynamic simulation? Well, the idea is that, as in the motion of the planets, we want to integrate the equation of motions with this caveat that, again, each particle interacts with each other particle. And so it's really a many-body problem, even if it is a classical many-body problem. And so it's really easier to integrate than a many-body Schrodinger equation.

And in general, the way particles interact with each other can be described by empirical force fields or, in the best-case scenario, can be described using quantum mechanics. That is the [? functional ?] theory [INAUDIBLE].

And so what you will have to do in that case is that for every configuration of your ionic system, you need to solve an entire total energy self-consistent problem. And that's why actually the first principle in molecular dynamics becomes a very, very expensive challenge. And we'll see in one of the later classes how this is solved.

There are a lot of systems that can be described accurately, even with very simple potential. Even more than in the case of the total energy, certain thermodynamic transitions, like the transition from a solid to a liquid, can be thought of as broadly independent of the details of the interacting potentials.

So even studying a system with a Lennard-Jones potential that, after all, was developed especially to do molecular dynamic simulations of rare gas atoms, or even just studying a system that is made by hard spheres that don't really interact when they travel in space, apart when they really hit each other-- so what is called a contact interaction that is 0 all the time apart from the instant in which the two spheres collide-- can sometimes give you a very accurate qualitative description of the thermodynamic transition in a problem.

And molecular dynamics actually was first introduced really to study systems, like liquids, for which it's very difficult to develop analytic techniques. You see, if you want to study a solid, you can actually choose a potential. And if the potential is simple enough, you can analytically calculate the potential energy of a periodic repetition of all your atoms.

But if you study liquid, disorder is really an essential part of the state of the system. And so the entropic contribution to the free energy of the system is fundamental. And really there are no reasonable analytic ways to address a lot of the questions that you have. And that's why people really moved to solving on the computer the Newton equation of motion.

And this is a short history, really, of molecular dynamics, reminding you what is the history of computers, really. In 1952, one of the first computers, MANIAC, was operational in Los Alamos. And that's when Metropolis, Nicholas Metropolis, in collaboration with actually a husband and wife couple, Rosenbluth and Rosenbluth and Teller and Teller, developed what is called the Monte Carlo method.

That was a very powerful method to address exactly this entropic free energy problem. And that you'll see in one of the later classes. But really, 1956 sees the first molecular dynamic simulation by Alder and Wainwright.

And as I said, in looking at the dynamics of hard spheres, these are the milestone papers. Vineyard is studying radiation damage in copper, what happens when you have a copper atom that is really put in motion by a very energetic radiation kicking it and starting running around.

And really this is the classic paper of Aneesur Rahman studying, with a Lennard-Jones potential, the dynamics of liquid argon. Again, a Lennard-Jones potential is actually very accurate to describe rare gases that ultimately repel each other at very close distance and, at large distances, attract each other with this weak dipole-dipole interaction that goes as 1 over r to the 6th.

And then really the first ab-initio molecular dynamics and the first theory for ab-initio molecular dynamics was developed in the mid '80s by Car and Parrinello. OK, now, there is a universal way to think of these dynamics of your classical system composed of N particles. And as we have said, since they follow Newton's equation of motion, what we really need to track, the only quantities that we need to track, is their position as a function of time and their velocity as a function of time.

And generally speaking, the kinetic energy is just a function of the velocities. And the potential energy is just a function of the positions. So if you have N particles, you really have to follow in time six n variables.

And if we define for a moment a six-dimensional space, what we really have is that an instantaneous state of your dynamical system, that is all its position and all its velocity, can be represented by a single point in this 6N dimension of space. So if you have got just a hydrogen molecule with two nuclei, what you really need is a point in six dimensions.

If you have a water molecule, you need a point in nine dimensions. And the evolution-- sorry. If you have a water molecule, you need a point in 18 dimensions. And if you have a hydrogen molecule, you need a point in 12 dimensions.

And the dynamical evolution of your system is described exactly by the evolution of this point. So we often think at this phase space that we'll have, for an N particle, 6N dimension, of which-- let's try to draw a bit dimension-- 3N of them, those represented in green, will really represent the positions. And the other 3N will represent velocities.

OK. And obviously this is all orthogonal dimension. And your point will be somewhere. The certain instant in time T will have a set of velocities, will have a set of positions. And the Newton equation of motion will make it evolve in this phase space.

So really your molecular dynamic integration consists in following in time this evolution. And that evolution is, in principle, analytically and uniquely defined once you know the position and the velocities at every instant. And-- well, we'll see that later.

OK, so what do we do once we have a computational algorithm that allows us to evolve this point in phase space? And I sort of summarized here the three main goals that could follow from a molecular dynamic simulation in the order that we'll see them.

Really, one of the most common ones is to use molecular dynamic simulations to calculate thermodynamic properties, that is to calculate what we call ensemble averages. Suppose that you have a system like a solid. And you have calculated, say, its lattice parameter, like you're doing in lab 2 or lab 3.

And suppose that you start heating up this solid. What happens? Well, the atoms, in the classical sense, start oscillating more and more. And you know from experience that most solids start expanding when you warm them up.

So a molecular dynamic simulation could tell you how much a solid expands when you heat it up because basically you start evolving your system from 0 temperature. And then you increase its temperature. And really, temperature, remember, is just the average kinetic energy of the ions.

And as you increase the temperature, the ions oscillate more and more. And if you let your unit cell expand or contract in time following what is the stress overall of all the atoms in your unit cell, you see that your system increases the dimension of the unit cell and starts oscillating around a value that's really an average, depending on temperature.

And all these sort of statements can be actually made formally in a way that we'll see in a moment that is the language of thermodynamics. That is, if we have a finite temperature, we'll have a certain probability associated with every microscopic configuration.

OK, so if you are at 0 temperature, really what you have is that you have a probability 1 to be in the lowest energy state possible and probability 0 of being in any other energy state. If you start increasing your temperature, states that are energetically similar to the ground state but not really the ground state start to have a certain probability of being occupied.

And the more temperature you have, the more very costly configurations can be accessed by your system. And really, the average properties at a finite temperature need to be an average over all these possible states weighed with their appropriate thermodynamic weight, weighed with the probability of being in that state.

And so this is what an ensemble average does really. It tries to span all the possible states in a system with a probability that is proportional to their thermodynamical weight. And then for each state, you can calculate a certain property, like what is its average potential energy, what is its lattice parameter, and so on. And out of that, you can calculate the average thermodynamic property.

And we'll see, actually, an example in a moment. So we'll see a practical application of this. The other thing that you can do with molecular dynamic simulator is actually study the evolution in real time of your system. So you could actually study a chemical reaction.

What happens if you have, say, an explosive molecule that is decomposing, OK? And in particular, whenever you have a chemical reaction that tends to be a bond-breaking and bond-forming reaction, you'll probably need to use a quantum mechanical approach.

But in principle, you can set up your system in an initial configuration and see what happens as you evolve in time. You could look at a catalytic reaction. Say, maybe you are interested in studying fuel cells. And you want to see how a hydrogen molecule decomposes when it arrives on a platinum surface.

Well, you can take your hydrogen molecule and project it delicately towards the platinum surface. And follow, with your molecular dynamic techniques, the dynamics of all these nuclei as the molecule chemisorbs and dissociates in time.

And then really the last application of molecular dynamics is more of an optimization algorithm. There are problems in which it's very complex to find what is the lowest possible solution, what is the optimal solution.

That is, you could try to find out what is the schedule of all your flights in an airplane company. And obviously that's a complex optimization problem because you can't move one plane at a time and figure out what is the best possible solution because there might be a completely different choice of planes and itineraries that actually gives you a best performance overall on the net.

And so you can always take one optimization problem and express a cost function. That is depending on what is your interest. Minimize the amount of oil, let's say, used in running all these routes.

So when you have a cost function of a very complex problem, you have really an object that depends on many variables and that has a lot of possible minima. OK, so this would be, really, the problem of optimizing the air routes of an airline. What you have is a lot of possible variables. And there are a lot of reasonable solutions.

But you really want to find the one that is the best for you, something like this. And molecular dynamics can actually help to solve this problem because it's a problem that is absolutely analogous of having a system of interacting particles that has a certain potential energy surface and trying to explore that potential energy surface, moving around, until you find the lowest minimum.

You could think of this as a mountain range. And you have your skiers that are really your particles moving around. And you want to sort of let the skiers ski around as much as possible until they find the minimum energy state.

And so you can use a thermodynamic analogy. You heat up your system. That really means that every particle has a lot of kinetic energy. So having a lot of kinetic energy can easily overcome any kind of potential barrier.

And you have a lot of these particles moving around. And then you cool them down slowly. And as you cool them down, they'll start piling up in the lowest possible energy state.

And once you have reached, very slowly, 0 temperature, you look at where all your particles have ended up. And so you have at least a good sampling of the relevant minima in your potential energy surface. So actually molecular dynamics-- and this is called simulated annealing, in this context-- can actually be used as an optimization technique.

And we don't have to go as far as the case of an airline or what is called the traveling salesman problem. But often what we have is that we have a system that is fairly complicated to characterize.

Say we want to describe maybe an interface between a perfect silicon crystal and the SiO2 substrate because we are in the electronic industry. Well, I mean, we have a problem because how do we construct the interface between a crystalline solid and an amorphous system?

I mean, we can try putting atoms in a reasonable position. But it's never going to work very well. And so what we can do is actually use molecular dynamic techniques to let the atoms evolve according to their interaction and to the Newton equation of motion and somehow find, by themselves, a slightly more favorable green minima instead of the initial configuration in which we have put them.

OK, so, in principle, if you had infinite computing power or infinite mathematical prowess, you would have really solved the problem of characterizing on a computer any material that you wanted. The reason why you can't really do that in general and indiscriminately has to do with four problems that arise in the approach that I've mentioned.

And I've mentioned them here, I guess, in order of importance. And the first and most dramatic problem is the one of time scales. That is, atoms move around and vibrate at time scales that are characteristic of atoms.

That is, say, in a molecule, it will take tens of femtoseconds or hundreds of femtoseconds to have a periodic oscillation. So these are femtoseconds. 10 to the minus 15 seconds is the order of magnitude of the dynamics of the atoms at room temperature.

Now, a lot of interesting problems take place in times of seconds, minutes, or hours. A classical problem is how a protein folds itself in its native state. And you know, this is a process that can take seconds.

And so we have, really, the issue that all our calculation takes place on the timescales of 10 to the minus 15. And often we need to reach physical properties that are of the order of seconds. So we need to span a 15 order of magnitude.

And that means that if it takes a time x to perform one single molecular dynamic operation, we have a 10 to the 15 cost to get to the reasonable time scale for a lot of macroscopic processes. And that's really a scale that we can't bridge.

If we do quantum mechanical simulations for a small system, we can easily reach the times of tens of picoseconds or hundreds of picoseconds. If we do classical molecular dynamic simulation, we can go maybe three, four, five orders of magnitude better. And so we can start accessing nanoseconds and maybe start getting closer towards the microsecond.

But there is still an enormous gap to reach. And this is really a difficult problem to overcome because we have dynamics taking place at different levels.

Say, think of a protein. We have the atomic dynamics of an atom vibrating around. And then we have dynamics of subunits maybe interacting with the different amino acids along the chain.

And then we have structural motifs that want to move around and then maybe the first coil-up in a certain configuration. But that's not really very good. And then after a while they sort of open up and recoil in a different configuration.

And so the more you look at your system on a larger and larger length scale, the more you discover that there are dynamics of more and more complex groups that are slower and slower. And integrating from the atomic motion up to the mesoscale and the macroscopic scale of these dynamics is really an impossible problem to do by brute force.

And I would say that it's still the conceptual problem for which less progress has been made. It's really the most difficult one. And you'll hear a lot of talks nowadays in science about multiscale modeling and about bridging timescales. And it's obviously a very important challenge. But I still think we are at the stage in which there are very little solutions.

Somehow related to the time-scale problem, there is the length scale problem. Suppose you want to study a system close to a phase transition. Suppose that you want to see a gas evaporate. You want to study a liquid-to-vapor transition.

What happens exactly at the temperature where the liquid and the vapor are in equilibrium? Suppose that you are studying water boiling. At normal condition, it boils at 373 Kelvin. And at that temperature, why there is a phase transition? Basically because the free energy of the liquid and the free energy of the gas are the same.

You go a little bit above the temperature, the free energy of the gas is lower. So the system wants to be all gas. You go a little bit below the boiling temperature, the free energy of the liquid is lower. So the system wants to be a liquid.

But as you get closer and closer to the temperature, these two free energies becomes comparable. That means that the cost of transforming from liquid to the gas and from gas to liquid becomes smaller and smaller. And it is 0 at the transition temperature.

That means, in practice, that since it doesn't cost anything to transform from one state to the other, you have lost all possible length scale because you can nucleate a bubble of gas of any size with 0 cost. So if you think of the transformation from a liquid to a vapor as the nucleation of bubble of gas exactly at the temperature of the transition, that cost has become 0. And we can nucleate bubbles of any size.

So the phase transition is a moment where we lose length scales. The fluctuation can have any size. And so any kind of finite simulation will break down because it's really only in the infinite system that we represent appropriately all these sizes popping up in our system.

And so that's why it becomes very difficult to study accurately, say, a boundary at a phase transition. We say that the length scale starts diverging. And that's why one needs to be careful.

This problem is not as serious as the problem of time scales. But it's still one that you want to consider very carefully and one you need to spend some time thinking, actually.

That is, often we use a periodic boundary condition to describe an extended system. And so the question that you always have to ask yourself is, how large should I make my unit system that is periodically repeated before it's really describing accurately the infinite system that is the goal of my simulation?

Suppose that you want to study water, OK? Or suppose that you want to study any liquid. In principle you want an infinite system. In principle you can deal only with a finite number of degrees of freedom.

So the two possibilities that you have, you could study a droplet of water, OK? But then obviously you encounter the problem that, when you study a droplet, you have a very large influence coming from the surface effect.

A droplet is a droplet because of the surface tension. And before reaching the thermodynamic limit of that droplet behaving as a bulk liquid, you need really to go to enormous sizes.

If instead you study, say, the same number of molecules in periodic boundary conditions, you eliminate the presence of the surface. And so you can reach the thermodynamic limit of the system behaving as a bulk much faster.

But still you need to always ask your question, is your simulation still large enough? And there is a size that is large enough. And that's the other very important concept. That is, whenever your simulation cell is so large that the dynamics of one water molecule doesn't affect the dynamics of a water molecule that is far away from yours as much as it is its periodic image, your system has become infinite.

That is, whenever you have a molecule in a liquid, you really dynamically interact not with every other molecule in the infinite liquid but only with a sphere of neighbors. And you need to understand how large is that sphere of neighbors because those and only those are the molecules that affect yourself in the center of that sphere.

And if you decide that for a water molecule at room temperature the sphere of neighbors with which there is a dynamic interaction has only a radius of 10 angstroms, whenever your unit cell is larger, is a cube, say, larger than 20 angstroms inside, you are studying, in practice, an infinite system because every molecule in that system sees all the independent neighbors that it needs to see, OK?

So even a simulation with a finite number of molecules can reproduce exactly the behavior of an infinite system provided a molecule interacts with other independent molecules up to a maximum range of dynamical correlation. And we'll actually define this quantity in the lecture that follows.

So again, length scales are another problem and a problem in which periodic boundary conditions tend to help a lot. But again, suppose that you want to study something like the strength of a metal.

Well, you have now a problem in which a lot of length scales becomes intertwined because in the strength of a metal, you might want to consider what is the dynamics of impurities that strengthen that metal. Like, if you have steel, you have carbon atoms in between the ions. And so you need to understand what is the atomic dynamics.

But then there will be dynamics of units on a larger length scale because in a metal there's going to be these location and slip and glide planes. And so there is coordinated motion of thousands of atoms that you want to consider.

And then, in reality, a real metal is never a crystal. But it's always a polycrystalline material. So you have grains. And inside that grain you have a consistent ordering of the atoms. But then grains are mismatched. And there are grain boundaries.

And so you have a lot of length scales to consider if you really want to predict what would be the mechanical response of a real system. OK, so this is the challenge, if you want, number 2.

The challenge number 3 is the one that we have described more in detail up to now that is determining how accurate can be your prediction of the interaction between nuclei. And you have seen everything that we could tell you about empirical potential quantum mechanical simulations.

And in principle, even if it's obviously very expensive to make very accurate predictions of the interaction between atoms, in principle it's the less conceptually challenging of the problems. I mean, we can just try all our electronic structure methods at will. And we can try to see how accurate we become.

And the last challenge that is also not really a conceptual challenge and it can be solved but is just a very expensive one to overcome is this one in which we treat the nuclei as classical particles. And that tends to be, in general, very good. But it's sort of an approximation that breaks down the more you go towards lighter and lighter nuclei because, remember, the de Broglie relation.

Your wavelength times your momentum is equal to the Planck constant. So the lighter you are, the smaller your momentum, the longer your wavelength, the more quantum you become.

So something like hydrogen that has a nucleus that is just one proton, so it's fairly light-- or it's the lightest of all possible nuclei-- tend to still have significant quantum properties even at relatively high temperatures. You see, the more you increase the temperature, the more the average momentum of your particles become and the more classical they become. But the lower the temperature, the more quantum they become.

And something like hydrogen is still significantly quantum at room temperature. And what does that mean? Well, it means that that particle, it's really a delocalized system that can actually tunnel through barriers. That's if you want one of the important differences between a particle that behaves according to classical mechanics and one that behaves according to quantum mechanics.

If you are behaving as a classical particle, you need to have enough kinetic energy to overcome a barrier. If you are a skier and you are in a valley and you want to go in the other valley, you need to have enough velocity to be able to overcome the mountain.

If you are a quantum skier, you can actually have some finite probability of tunneling through. And so there are tunneling problems that are important not only for hydrogen but sometimes even for much more massive particles or even for groups of particles, provided that the energy barrier that they need to overcome is small enough.

And I've shown you in the first quantum class an example in a perovskite crystal in which we have this cubic structure with an octahedral cage of oxygen that can configure itself in a ground state that is ferroelectric, in which the oxygen cage displaces with respect to the cubic symmetry and induces a permanent electrical dipole in your system. Well, there are certain materials that have this structure.

But since the octahedral cage can have roughly six possible choices on where to move offsite, if you decrease the temperature enough, the oxygen cage will start actually tunneling between all the six possibilities. Instead of finding itself at 0 temperature in a broken symmetry configuration in which the oxygen cage has chosen one of those six possibilities, in reality, because of quantum tunneling, it keeps tunneling between all the six possibilities.

And on average your system looks cubic. And on average you don't have really any ferroelectricity. But the polarizability of the system is going to be different from that one in which the octahedral cage was sitting squarely in the cube. So there is a dielectric response that is different.

So there are cases that are fairly exotic, I would say, for most applications, in which tunneling effects of the nuclei could become important. But there is one other important consequence that comes from quantum mechanics that tends to be important in a lot of cases. And that consequence has to do with the quantization of vibrational excitations.

That is, if you look at the motion of any system-- a molecule, a solid-- in classical terms, well, that molecule or that solid can have any amount of kinetic energy. In a hydrogen molecule, the atoms can have 0 kinetic energy if we treat it classically. Or they can move just a little bit. So they have a very minute amount of kinetic energy. And we can just increase the kinetic energy that those atoms can have at will.

In reality, because the system is a quantum system, the vibrational modes of a molecule or the vibrational modes of a solid are actually quantized. And you can't have arbitrary amounts of vibrational energy. Actually, in its quantum ground state, a molecule will not be at rest. But it will be in what is called the zero-point motion state. That is what is the quantum mechanical ground state.

You can think of the analogy with the stability of the atoms. In its ground state, the electron around the proton is not collapsing onto the proton because of the Coulombic repulsion. If the electron was a classical particle, it would collapse on the nucleus. And the hydrogen atom would be just a proton and an electron sitting in the same place.

But because the electron is quantum, it can't really collapse. And its lowest energy state is a state that is the 1s orbit for the hydrogen in which the electron is around the atom but doesn't collapse onto it. And the excitation of that electron to the next energy state is a quantized excitation.

You can't go to a level of 1s that has an energy for the hydrogen of minus 1 Rydberg to another state that has energy minus 0.99 or minus by 0.98 or any kind of amount. The next state in which the electron can live is a state that has an energy of minus 1/4 and then minus 1/9 and so on and so forth.

The same thing happens for the vibrational level of molecules, OK? A quantum molecule in its ground state is not in a state where the atoms do not move. But it's in a state where the atoms have a certain amount of quantum kinetic energy. And it's in a state called zero-point motion.

And if you want to increase the kinetic energy of this molecule, you can't do it continuously. There is a quantized cap. So you need to provide enough energy to heat up the system to the next state and then to the next state and then to the next state.

And that, if you think, affects deeply how you are going to heat up a system. That is, if you look at the specific heat of a solid, that is basically the quantity that tells you how much heat the system is going to be able to take, and you look at what happens at very low temperature, you discover that actually the system is almost unable to take any heat because all the vibrational excitation are much larger than the average amount of kinetic energy that you have if you are looking at what it takes to bring something from 1 Kelvin to 2 Kelvin.

OK, but this is still relevant if you have a system like water, even at 300 Kelvin. In water, say, at 300 Kelvin, what you have is you have these water molecules. Remember what you have-- your oxygen, hydrogen, and hydrogen.

And the modes of vibration of these molecules are stretching modes of the hydrogens. And then there are scissor-like models in which the two bonds do this. And then these are intermolecular modes, just sort of relevant to one molecule. And then there are all the molecules interacting.

Well, at room temperature, all the internal modes of the molecule, the stretching and the scissor modes, have so much vibrational energy in their ground state that it's almost impossible to excite them to the next vibrational state with the average temperature that you have at 300 Kelvin. So even a system as fundamental as water has really a hybrid behavior in which the vibrational interaction between molecules can actually be populated arbitrarily by the temperature that is surrounding you.

But it's almost impossible to excite a molecule from a zero-point stretching motion to the next one. And so, say, the specific heat of water will be slightly exotic because really part of its vibration are quantized.

And so there are a number of cases in which this classical nuclear approximation is not exact. And luckily they can be studied with specific techniques that we'll see in some of the later-application classes during this course.

But remember this-- if there is one thing that you want to remember about the behavior of nuclei-- that this vibrational excitation are quantized. And we actually see this quantization in physical properties like the specific heat.

OK, so let's discuss the first of the possible applications of molecular dynamics that I have mentioned. And that's using molecular dynamics to really calculate average properties of a system at finite temperatures.

And what I have written here is really the summary of statistical mechanics. So that is what, if you want, Boltzmann would tell you. You have a dynamical system that is described by this point moving around in phase space, OK? So with a certain possible set of positions and with a certain possible set of momenta.

And suppose that we want to calculate, say, what is an average property at a certain temperature. And we call that property A. It could be the lattice parameter of the system. Or it could be the diameter, say-- it's even simpler-- of your water bubble.

OK, you have this water bubble. And you want to see what is the diameter of this water bubble at a certain temperature. And this diameter is going to increase with temperature. And what you do, you just let this system evolve. And during its evolution, you keep measuring what is its diameter. And you let this evolve a lot. And you measure what is its width.

Well, the way Boltzmann would set this problem is say, we really need to integrate over all possible configurations of the system. That is, we need to go through all possible configurations of the system. That is, we need to integrate over the full phase space.

And for every possible point in phase space, we need to calculate the diameter of my water bubble. And then we need to take into account what is the probability of being in that configuration. And if you are in a system at a constant temperature, the probability of finding yourself at a given temperature is given by exponential of minus beta is just 1 over the Boltzmann constant times the temperature times the energy.

That means that configurations that have a high internal energy are very unlikely. And configurations that have lower energy are more likely, to the point that if the temperature is 0, that is this beta is diverging to infinity, only the lowest energy configuration possible is represented.

So this is the statistical mechanic points. To calculate an average quantity at a finite temperature, we need to sum over all possible configurations in space and momentum. For each possible r and p configuration, we'll have a certain value of the diameter of that bubble.

We weigh that with the statistical weight of that configuration. And we actually normalize this expectation value. And this is what we obtain.

So if, for a moment, you think at a phase space that is two dimensional and in this square we represent all the possible positions and velocities of our system, what Boltzmann tells us is that for every possible point, we will have a certain value of A. And the sum of all the possible values of A weighted with the probability of finding them gives us the expectation value.

And so what you really need to do is an integral over all space of this. And the cost of these integrals diverges as the number of dimensions diverge because, again, every time we increase the size of the system by one particle, we have six more dimensions over which to integrate.

And again, if you think that we do this just by numerical integrations by taking points on a grid and maybe we just discretize our system with 10 points along one dimension, every time you add a particle, you have 10 to the 6 points that multiply your system. So this quantity explodes.

The point of view of molecular dynamics is different. What we say is that, well, let's just let this system evolve. Let this system evolve in time. And let it evolve at a constant temperature.

And so now our point in phase space is going to go around, following a trajectory. And if we let it go around long enough, it will spend more time in the parts of phase space that are more favorite, the ones that are more likely, according to the Boltzmann factor.

And so if our evolution, really, is distributed in the trajectory according to this thermodynamic weight, thermodynamic factor, then our expectation value, the thermodynamic average for our observable A, is just trivial, obtained by integrating over the trajectories all the values of A. That is what I was saying before.

If you want to calculate the diameter of your bubble, what you do in the molecular dynamic sense, you just let the bubble evolve according to the equation of motion at a constant temperature. And you keep monitoring its diameter. And you will see that the diameter will have an average and that this is a much more efficient way then considering all possible configurations of the bubble, weighing them with how expensive that configuration is, and doing it for all the possible configurations.

The molecular dynamic algorithm automatically, if you want, keeps spanning the most likely configurations. And so it does, for you, the job of going in the place where it's most likely to be.

So in order to calculate a sample average, we can actually transform a problem into integration over a trajectory. And that is almost always doable under what is called the criterion of ergodicity. So the equivalence between the two formulations that I'd written before-- that is calculating a thermodynamic observable as an integral in phase space or as an integral over trajectories-- are equivalent, provided your trajectory is really able to go all over the place in phase space.

So if, again, this were phase space, you could have, in principle, the problem of a known ergodic system if there is something that prevents your trajectory to reach certain so-called forbidden segments in your phase space. And I can give you an example.

Suppose that you are studying something like silicon and you are in the crystalline state. And now you heat it up. And once you go above the melting temperature, your silicon will become liquid, OK?

Now, suppose that you cool it down again below, now, the melting temperature. OK, in principle, the liquid silicon should become again crystalline silicon because, at that temperature in which we have brought it back to, the thermodynamic stable state is really the crystalline solid. If you do that in practice with our molecular dynamic simulation, you see that you start from the crystal. You heat it up. It becomes a liquid with no problems.

And then you cool it back down. And instead of going back to a crystal, it goes back to a glassy, amorphous state. Silicon will find itself in a configuration in which locally all atoms are happy to be fourfold coordinated. But if we look at the long-range order, it's completely lost. The system is not crystalline anymore. It's glassy.

What's happening there? Well, it happens that it's just, for a liquid cooling down, so much easier to trap itself into a set of configurations in which atoms are really fairly happy as they are, because what is important for silicon atoms is to be fourfold coordinated.

But what is important to be in a crystal is that this fourfold coordination extends and becomes a long-range order. But for silicon, it's very inexpensive to mess up that long-range order, like for many other systems.

And in principle, same thing happens for polymers. You lower it down. And they sort of get trapped in a glassy configuration. A lot of polymers might want to stay actually in a crystalline state. But really that is impossible to reach.

And so if you want to study, say, what are the properties of your silicon and you have obtained or you are trying to obtain your silicon structure from a liquid that you cooled down, you'll actually-- in your molecular dynamics, you'll keep orbiting around in this glassy mass state and will never reach that region of phase space that is actually very, very small because it needs to have all atoms long range in very specific positions.

And so you are really spanning a system in a nonergodic way. That is, you spend a lot of time in a region of phase space that are not really representative but out of which is very difficult to escape.

OK, so there are problems in, again, brute force molecular dynamic approaches that you have to be careful about. OK, so let's see how we actually go and try to do this on our computer experiment. And I'll start with the simplest case that is what is called the microcanonical case.

That is, we perform a simulation where the number of particles is constant, the volume is constant, and the total energy of the system is constant. This is what is called a thermodynamic ensemble. Again, these are all extensive variables. And we are fixing them.

Sometimes it's appropriate, instead, to do a Legendre transformation and go into a different thermodynamic ensemble. Say, you could want not to be at a constant volume. But you might want to study a system at a constant pressure where the volume changes.

And so we need to describe how we'll actually be able to do that. That is, how will we be able to move to a system in which different intensive thermodynamic variables are constant?

But if you take a system of particles in periodic boundary conditions and you have 10 water molecules in a certain unit cell and you evolve them with the Newton equation of motion, you can obviously know that the number of molecules and the volume of the unit cell are not going to change. And if you integrate your equation of motion appropriately, the total energy of the system, that is the kinetic energy plus the potential energy, is not going to change.

And so as we evolve this point in phase space, its trajectory is going to span what we call the microcanonical ensemble. And you let this point evolve long enough. And all of a sudden, you will have a meaningful trajectory, that is a meaningful assembly collection of microstates of representative configurations in space and in velocities of your system, out of which you can calculate averages of any of the observables that you want to calculate and that you want to characterize.

So what do you do in these molecular dynamic simulations? Well, there are really four important steps that you need to take into account. And we keep thinking of water, liquid water, as an example. That is, you need to start from something. That is, you need to put your water molecule in a certain position. And you need to give them initial velocities.

And you are almost never going to be good enough in choosing a state of the system that is somehow meaningful. You'll always, when you start building something complex, put molecules in places that they don't want to stay. And one needs to be very careful, then, in all of this.

But suppose that we have, in any case, an initial configuration of position and velocity. Then what we need to do is, over and over again, integrate the equation of motion, that is, from the current position and the current velocities, predict the new position after a certain time.

And once you have the new position at a certain time, you can recalculate what are the forces acting on the particle in the new position. And with these forces, you can, again, move the particles by a certain amount and so on and so forth.

And that is where the molecular dynamics algorithm comes into play. We need to integrate our equation of motion. And with this, our system will start to evolve.

We have put that point in phase space in motion. And it'll start moving around. And if we are integrating correctly, the total energy of the system is going to be conserved.

Now what we need to do is we need to make sure that our characterization, our thermodynamic averages, do not depend on our initial condition. Again, we have started from a configuration that we have built by hand. And because of that, it's most likely a very unfavorable configuration of the system.

OK, so we need to let the system evolve and somehow start figuring out really in which condition it would want to stay. Think of the previous example of silicon that is glassifying. Well, we had liquid silicon. We have to cool it down.

In a way, what's happening is that we are never equilibrating it because the system goes into a glass phase. And so it never loses memory of its initial condition. It's always sort of trapped in this pseudo liquid, now glass, state. While in reality, at a temperature below the freezing or the melting point, equilibrium would be given by a crystalline state with the ions vibrating.

So there are cases in which equilibration is so difficult that we actually never reach it. Silicon never goes to a crystal. That's why we need to be careful. But in most cases awfully we are lucky. And we lose memory of initial condition.

Say, the opposite is actually very easy to achieve. If we start from the solid and we warm it up, it becomes a liquid. And if we wait long enough, that liquid will have lost all memory of having come from a solid.

OK, so it's actually very easy in one direction and very difficult in the other direction. But once this equilibration time has passed, well, at that point we have truly everything in place. Our point in phase space is moving around. And it's really going into the regions of high probability and distributed, according to Boltzmann distribution.

And so at that point, we can start an average among any of the physical observables that we want to calculate. And so there are these four phases. And I'll briefly mention how we actually choose them.

Initialization. Well, again, we need position and velocities. You need to be careful, in particular, to avoid overlap between your molecules.

If you put two atoms very close, they are going to repel each other with an enormous strength. And if you put them close enough, they are going to shoot in opposite directions very quickly and completely kick off and mess up everything on the side.

So what is difficult is really, always in configuration space, finding the right structures. This is the difficult part of choosing initial condition. Choosing velocities is much easier. One possibility is just giving zero velocity to the system or a very small velocity and then slowly adding things.

And what we usually see is that configuration in velocity space is much faster. And this is the reason. Think of that trajectory for a moment of one particle. And think of what happens when you have a collision.

Suppose that you have particles, for a moment, that are hard spheres. In a collision, the velocity of the particle all of a sudden changes sign. So in the velocity part of the phase space, your point makes a sudden jump, OK, because the velocity that was plus 10 has become all of a sudden minus 10.

So in velocity space, for the limit case of a system of hard spheres, your representative point hops everywhere continuously. In reality, I mean, you have interaction at all times. And you don't really hop. But you can sort of move very quickly from one region to the other region.

There is not anything similar to the glass entanglement in configuration space. So it's very easy to thermalize velocity. It's very difficult to thermalize and equilibrate position. And that's where, all, you need to be careful.

And if you really want to be accurate, what you can also do instead of giving zero velocity to your system, you can use your concept of statistical mechanics. That is, a system of classical particles will have a distribution of velocities that depends on the temperature at which it is. And this is 0 in Celsius. So it's 273 Kelvin.

And the distribution of velocity in a classical system is actually given by the Maxwell-Boltzmann distribution written above that has basically this shape. So it means that if we are, say, at 0 Celsius, 273 Kelvin, we'll have a lot of molecules with a certain velocity. And then, in principle, we can have molecules with higher and higher velocity. But the probability of finding them decreases.

And as we increase the temperature, the probability distribution of the velocities for the molecules changes. And so as usual, you can define both a median or an average in your system. That is, you can look at what is the most likely velocity, that is the top of this curve, at 0 Celsius. Or you could also define what is the average velocity. And since that this curve is asymmetric, it's going to be really shifted a little bit.

And these are actually the numerical expressions for these two velocities, the most likely and the average velocity. And so I've given you the example of an oxygen molecule at room temperature.

We are talking about 1,000 meters per second. So these velocities, if you want, are similar to the speed of sound in your material because, at the end, sound propagates at a certain speed because that's the velocity of the collision from one molecule to the next that are really carrying sound away.

So if you want to be sophisticated, you give to your system a set of velocities that could be, say, all identical to the most likely velocity. But again, equilibration and thermalization in velocity space is very easy. So once you have this, that is, once you have a set of positions and the set of velocities, you need to use Newton's equation of motion of saying, of predicting, where your system is going to go.

And so that's when the molecular dynamic algorithm comes into play. So what you need to do is, as we say, you need to use an integrator, the algorithm that tells you what your next position and what your next velocities are going to be.

And there are a number of integrators. There are subtleties. But often you find the name of Verlet algorithm. This is one class of algorithms. And the other common class is the Gear predictor-corrector, just to give you the terminology. And we look, in a moment, at the Verlet algorithm.

They're all fairly simple. But what is really important is that this algorithm should be robust. That is, they should give you a trajectory that is very close to the trajectory that you would get from the perfect analytical integration of the equation of motion.

And in particular, what is the most important thing is that they should be accurate in, say, conserving the constant of motion. As we have seen from the Newton's equation of motion, the sum of the kinetic energy and the potential energy should be a constant.

If your integration, if your numerical evaluation, of that equation is poor, you will see that the constant of motion is not any more constant. So again, this is the sanity check on the accuracy of your integration.

And there are more subtle elements of the integration algorithm that can become important. But we won't go into that. And then as I said, straightforward integration of the Newton's equation of motion gives you the microcanonical ensemble in which the total energy's conserved.

But sometimes you want to study systems in which the temperature is fixed. Or you could want to study systems in which the pressure is fixed. And so you need to augment your molecular dynamic integrator with some feature that allows you to control those intensive variables.

So if you want, say, to control the temperature of the system, what we say usually is that we couple the system to a thermostat. That is, we make sure that the temperature in your evolution doesn't change. And you'll see examples in the next class of this.

And in the case, in particular, of something, say, like the temperature, you have a whole set of possible thermostats. That is, you could use stochastic approaches. That is, you have your evolution of your system. And your system doesn't conserve, if it's microcanonical, the temperature. Sometimes it goes towards higher temperature. Sometimes it go towards lower temperature.

And what you could do is either every now and then absorb some energy or give back some energy in your system-- and that would be a Langevin dynamics-- in order to control the average kinetic energy of your system. Or you could somehow forbid, in the microcanonical evolution, your system to choose a state in which the temperature changes.

So if your new positions are such that your velocity gives an increase in temperature, you just shorten how much your atoms move in order to make the velocity smaller and in order to make the kinetic energy conserved. Or there are more complex dynamical ways of imposing temperature that goes under the name of extended Hamiltonian, or extended system, of which we'll see an example that is [INAUDIBLE].

But let me give you really what is the simplest and still one of the most used numerical integrators. And it tends to be very robust and very accurate. So it's still used in a lot of molecular dynamic simulations. And there is often no reason to do anything better than this.

So in a way, we are in a much luckier situation than in quantum mechanics. This is good enough in almost all cases. And this is how we look at the problem. That is, what we need to do always is we need to calculate the position of each particle at a time t plus delta t once we know the position, the velocity, and the forces on that particle at time t.

So suppose I'm at a certain instant in time that I call t. I know everything about my system-- position, velocity, and forces. What I really need to do is predict what are the new positions. And when I'm in the new positions, I want to calculate the forces again and evolve the system.

And now, what do we do? Well, let's do a simple Taylor expansion. That is, this is r of t is a function of the variable t. So we can use a Taylor expansion. And r of t plus delta t is going to be r of t plus the derivative of r with respect to t times delta t.

And that's the velocity. Plus 1/2 of the second derivative of the position with respect to time. And that's the acceleration. Times delta t squared, and so on for the third derivatives and so on for the fourth, fifth, ad infinitum.

This is just the Taylor expansion as a function, OK? And so it gives us the position of t plus delta t, knowing all the derivatives-- first, second, third, fourth, nth-- at times t. And once we have the derivatives, we can not only predict what the position would be at t plus delta t, but also what the position would be at t minus delta t. And it's written here.

And if you think, the terms with an odd power change sign in delta t. And the terms with an even power in delta t do not change sign. So we have just done a Taylor expansion. And now what we do is we actually sum these two quantities together, OK?

And this is what we obtain, OK? We sum the two quantities together. All the odd terms disappear. This sums to 0. This sums to 0. The fifth order, seventh order, and so on sum to 0.

So what we have is that we have rt plus delta t plus rt minus delta t is equal to 2 times rt plus alpha-- sorry-- acceleration times delta t squared. It's just written here. And then rearrange the terms.

But you see, this is the Verlet algorithm. Now what we have when we look at this expression, we have an expression that tells us what is the position at my new time step, knowing the position at my current instant, knowing the position at the previous instant.

We are keeping configurations from our trajectories. So we know where we were before. We know where we are now. And we know where we are going to be just by adding to this the term acceleration at the present time times delta t squared.

We make an error because we are not introducing terms that are higher order in delta t squared. But you see, the term in delta t cubed has canceled out. So the error that we make is just of the order delta t fourth. That is, that means that the smaller we make the delta t, the smaller we make the time step, the more accurate our integration is.

And so this is, if you want the most critical parameter in your molecular dynamic simulation, you need to make sure that the steps that you take are small enough, considering what is really the variation of your force as a function of time because remember that now that we use Newton's equation of motion, the acceleration at a time t is just given by the force divided by the mass. And because the force is just the gradient of the potential, it's just given by this.

So at a certain instant in time t, you have the acceleration acting in your particles just by taking the force and dividing by the mass. So at every instant, you calculate force. And the force lets you evolve the system. And it keeps giving you the new positions.

And so you know what is the trajectory and configuration space. And you can very simply find out what is your velocity at every instant t, if you need to calculate what is the velocity just by doing a final difference calculation of the tangent basically. The velocity is the tangent of the trajectory with respect to t.

OK, this concludes class. What we'll see in the next lecture is how we actually take this Verlet integration and how we use it to do the characterization of the thermodynamic properties of a real system.

So today is Thursday. Next week is spring break. So we'll see each other again on Tuesday, 29th. And I will be at a conference. And all the TAs will be at a conference next week.

So it will be slightly difficult communicating via email. So try to ask a lot of questions between now and Sunday for your problem set that is due the Tuesday after spring break. And that will be problem set 2. Otherwise, enjoy spring break.