And we can start to think. Sandro, if you remove the chair or whatever it is, that would be even better. Chair? What is it? Yeah, yeah, now it's OK. Now it's out of the window. Oh, excellent. OK. I actually don't see myself, but will somebody pin my video now? I see all now. Oh, OK, good. Now you can see me. Great. This is an experiment, by the way, huh? Yes. Massive partito. Yeah, this is actually better. OK, I think we can gradually start. So welcome back everyone to the school. Today we have a, it's my pleasure to introduce Professor Sandro Scandal, a pioneer contributor of condensed matter physics with studies in material science and atomistic simulation. Sandro is a long term member of the ICDP. Initially the position of the scientist then became the head of the condensed matter session and now acting as a research director of the center. Sandro today will speak about Abinicio molecular dynamics with the, and actually we have the opportunity to test this innovative lecturing methodology where Sandro is hidden somewhere at the ICDP. Nobody, apparently there is a legend on going here in the center that nobody knows where he effectively is and actually will be, the game of today will be to find him out where he's hiding in the center. But meanwhile, we'll enjoy his talk and so Sandro stages yours. Welcome, also welcome to the school from the organizer. Thank you. Thank you, Ivan. Thank you everybody for connecting so early in the morning.  So this is actually a new setup that was created by our scientific fab lab, the sci-fab lab here. It's still rudimentary, but I'm actually testing with this lecture. So you can consider it also as a test of the system and I hope everything works today. So the lecture today will be about Abinicio molecular dynamics. You'll have some tutorials later on, hands on. But since I guess most of you come from the background of electronic structure calculations, the idea is that I will start with some introductory concepts in the field of molecular dynamics. So my lecture will be structured in more or less two parts. First, I will give you some brief introduction to the basic concepts of molecular dynamics. And then I will move on to the concept of potential. And this is where the Abinicio part comes in. And I will also try to show you some simple methods to do Abinicio molecular dynamics, among which the famous Carr and Parinello method. I mean Roberto Carr and Michele Parinello they developed this method actually here in Trieste almost 40 years ago, 35 years ago. But it's still a method that is used in practice for Abinicio molecular dynamics. Okay, so let me start. It's some very basic concept. So what is molecular dynamics, first of all? Molecular dynamics essentially consist in simulating the dynamics of molecules or in fact the name molecular dynamics is an historical in some sense mistake because we are not just dealing with molecules, we're dealing with atoms, we're dealing with any particles as a matter of fact, as long as it moves with Newton's equations. So it should be called actually atomic dynamics more than molecular dynamics.  You have some forces acting on the atoms. Let me introduce some notations here. I will consistently use capital letters for ions, atoms and small letters for electrons because electrons will come into play later on, of course because we're doing Abinicio molecular dynamics. For the time being, it's just classical, okay? I'm just dealing with the classical system of particles. Each particle will have a force attached to it. And the idea is that we want to essentially solve Newton's dynamics for these particles, right? So mass times acceleration, this is the second derivative and this is equal to the force acting on that particular particle, okay? So we want to evolve, we want to solve this problem. And in solving this problem, let me write away mention that there are two aspects which are quite complimentary. One is the aspect of how we're solving this as a second order differential equation because this is what it is in practice, right? It's a second order differential equation in time. And there is another aspect which is in fact probably the most interesting one in some sense for us, for our community, for the Abinicio community which is how I calculate the forces acting on the particle. And this is where of course the Abinicio theory comes into play. Now, when you set up a simulations, you of course you have to define the atoms that you have in your system. For example, this could be a four atom molecule and you want to study the dynamics of this for a certain amount of time. But in most cases you do Abinicio molecular dynamics to study behavior of bulky systems, right? Liquids or collections of molecules, water for example, or amorphous systems. So your goal is to study systems that are actually macroscopic in size. And I don't need to explain you this concept, right? You're already familiar with the fact with the concept of a simulation cell where you of course you restrict your system into a given box. And you're also familiar with the concept of periodic boundary conditions. But in the concept of periodic boundary conditions in molecular dynamics is a bit more, I mean it's a bit less obvious because you're dealing with systems as move and you're dealing with systems at finite temperature. You're dealing with liquids, right? So there are no periodic boundary, there's no periodicity in a liquid. So what you actually do in practice, if you have a say a macroscopic system of particles now, let's imagine that this is 10 to the 23rd particles which is roughly the size of a macroscopic system. You don't want to study 10 to the 23rd particles in the computer. So you need to, and they're all different in some sense, right? There's no periodicity in principle if this is a liquid. So what do you do in practice?  It can be a hundred atoms, it can be a thousand, it can be a million atoms depending on the size of your computer, the one that you have available. So the first question that you ask yourself if you do molecular dynamics is, how do I truncate my system, right? Which kind of boundary conditions I'm going to apply I'm going to apply to my system. Now, you're familiar, as I said, with periodic boundary conditions, but this is not obvious at all because periodic boundary conditions in principle introduces porous periodicity in your calculations, right? If you're really, if you're studying a liquid with periodic boundary conditions, you're not studying a liquid. You're studying something which is periodically, in fact, you're studying a crystal if you wish because even if your simulation box contains, a hundred, a thousand atoms, these thousand atoms repeat periodically, right? So there are some artifacts that are introduced by the introduction of periodic boundary conditions. In fact, in the early days of molecular dynamics people were not actually using periodic boundary conditions. They were using open boundary conditions. They were putting the atoms in a selecting, say a chunk of the system and then just put vacuum outside and just calculating the forces between the particles inside this cluster, let's call it cluster. Now, there was a debate, especially historically, whether open boundary conditions were better or worse than periodic boundary conditions. periodic boundary conditions, from the computational point of view is likely more complicated because you have to compute forces coming from infinitely many particles. While with open boundary conditions, you're only looking at forces coming from finite number of atoms. So from a computational point of view, open boundary conditions are slightly easier. Of course, open boundary conditions means that you're exposing a number of atoms to a surface, to the vacuum, right? Which is clearly unphysical if your goal is to study a bulk system. I mean, and this is actually quite serious, more serious than perhaps you may imagine. Just to give you an example, I mean, suppose you take a thousand particles, right? And you put them in a cube with open boundary conditions, right? Suppose that the system is simple cubic, right? So you have a thousand particles, you have a cube, and the cube is upside, say 10 times 10 times 10, okay? And you're doing this with open boundary conditions. Now, can you guess the number of particles that are at the surface of this big cube? Well, I mean, I'm sure you're now trying to calculate it, but you have six phases and each phase is about a hundred particles. So you immediately realize that more than half of the particles that are in a cube of side 10, so a thousand particles are actually at the surface of this cube, right? You thought that a thousand was a big number, but you realize that you try to take a cube made of a thousand particles, more than half are actually at the surface of this cube.  So there is actually a half of the forces acting on those particles. So you see, you immediately see that open boundary conditions really need to go to very big simulation boxes if you want to achieve, you know, something that looks like a bulk state. So it's something where the vast majority of the particles in your system is actually in a bulk state. I mean, surrounded, completely surrounded by other particles. So the concept of periodic boundary conditions been, in fact, taking over as the main way to do boundary conditions in molecular dynamics. Now there are even more exotic choices. People have experimented something like twisted boundary conditions. These are actually quite common in the field of when you simulate wave functions because you're also introducing some K point sampling, but that's a different, that's a different story. So we discussed periodic boundary conditions. We social concluded that what we want to do is to study essentially systems in inside some simulation box, which is representative of a large amount of particles. We put our particles inside our simulation box and we now want to evolve them according to Newton's equations. And we use periodic boundary conditions. That is, if a particle gets to the surface of this box, it will automatically reenter from the other side of the box. We use periodic boundary conditions, okay? So as a matter of fact, what you're studying is not just, you know, completely dissolve the system, you just study a liquid, but you're studying a system which is periodically repeated in space. And what you have to make sure is that the size of your simulation box is not introducing some spurious physical, I mean, unphysical effects in the results of your simulation. So you normally try to do bigger sizes, bigger and bigger sizes to try and see whether there is anything in your calculations that depends on the size, which is of course something you have to fix, of course, in your simulations. Okay, so let me now go back to the solution of the second order equation, which is one of the first steps you have to consider, you know, how am I going to solve this second order differential equation? So in order to do that, oh, well, let me actually share the screen because it's probably faster if I do that. Let me go straight to this one here. Okay, I hope you can see this. Yes, I'm sharing the screen. So as you can see, meanwhile, I will erase this blackboard. As you can see, I mean, the goal is to essentially evolve the position of a given particle, which is our capital I with time. And the way you do that is by discretizing space, is discretizing time, sorry. Discretize your time in small steps. And I'll come back to the choice of the time step in a moment. And essentially what you want to do is to guess, well, try to come up with the best guess of what will be the position at time t plus delta t, once you know the properties of your system at time t, for example, position and velocity, and of course, the force because you're calculating the force, right? So based on the force, you want to guess what will be the coordinate of the position atom I, for example, in this particular example, at time t plus delta t, once you know the position and velocity at time t, right? And this is going to be an iterated procedure.  Now, there are a number of subtleties in doing this. The most important concept that you have to consider is that you want to minimize that particular error that you're making when you truncate this sum, right? Because you know the position, you know the velocity, you know the force, but you don't know what comes next, what comes at the third order in delta t, right? This is the third derivative of the position with respect to time. So if you truncate your estimate just after the force, you're going to make an error, which is of order delta t cube. So you want to try and minimize as much as possible that error, okay? So there are a number of tricks and there are several algorithms that try to improve this, but let me just give you a flavor of how these algorithms work with a given example. This is called the Verlet algorithm. The Verlet algorithm starts from the consideration that if you actually do the same exercise with minus delta t, and I'll show you in a moment why it's interesting to look at minus delta t, what you'll find, of course, is something where the terms that change sign are the velocity term, the ones, of course, which are odd in delta t, and of course the term in delta t cube, whatever it is, right? You don't know what it is, but you certainly know that this term is going to change sign when you reverse the delta t. Of course, the term in delta t fourth will have the same sign and will have the same sign. Now, if you sum these two expressions, you just sum the left-hand side and the right-hand side, and what you see is that two things. First, the velocity disappears, right? You don't need any longer to know what the velocity is. You just need the positions. Of course, you need to know also the position at time t minus delta t. In other words, your estimate of the position at time t plus delta t depends on the position at time t, position at times minus delta t, the force at time t. And then the interesting thing is that whatever remains, the residual of this calculation is now one order higher in delta t. It's now over delta t to the fourth. So essentially, you've been able to reduce the error that you are making. So this is a very popular, it's probably the most popular algorithm in molecular dynamics because it only requires the knowledge of the positions at the previous step and the two steps before the one you're calculating. And of course, the force because the force is what guides you in the definition of the trajectory. And the error you're making is to order delta t to the fourth. Now, you can easily imagine that you can make this algorithm following similar procedures, similar derivation. You can actually, by including memory of the system at time t minus two delta t or t minus three delta t, you can make that error delta t fourth even smaller by going to delta t to the fifth delta t to the sixth and so on and so forth. So it can be systematically improved following essentially the same derivation as the one that I've just shown you here for the first order, essentially. So it is, I mean, you have in principle the possibility to make it systematically better if you want. But the Verlet algorithm is already an extremely stable algorithm. Another reason why Verlet algorithm is quite popular, especially in molecular dynamics is that molecular dynamics simulations are never too long in time because abynition molecular dynamics is much heavier computationally speaking with respect to standard molecular dynamics. So with the standard molecular dynamics, you can run systems for, even for microseconds and there have been simulations that have been running for milliseconds. With abynition molecular dynamics instead, you rarely go above, say, one nanosecond. The typical abynition simulations is between one picosecond and one nanosecond. So it's actually quite short in time and errors that come from the integration of the equations of motions are less serious than errors that you do in the case of classical molecular dynamics where you evolve your trajectory for a much longer time. Now, let me actually say something about delta t now. How do we actually choose delta t? So let me actually go back to my screen here. And so I have my time here, my position here as a function of time. And I'm discretizing my trajectory, my time, right? The distance is delta t and I'm, of course, I'm evolving my trajectory in this way. And so on and so forth, right? So I have a discrete set of points which describe more or less the trajectory of my particle, right? And these are discrete points, of course. So let me, sorry, let me highlight the fact that what I have is not, of course, a continuous set of points, but it's set of points that define only where the delta t, every delta t, right? So this is what I call my trajectory. Now, the question is now, how do I choose delta t? Of course, I would like to make, there are two contrasting arguments that I have to take into consideration when I choose delta t, right? Based on what we just discussed, I want to make delta t as small as possible because the smaller is delta t, the better will be my integration. The more precise will be my estimate of the position of time t plus delta t, right? I'm making this error which I've now discussed, which is of order delta t to the fourth in the case of the Verlet algorithm, I'm making it as small as possible. So on one hand, I would like to make delta t as small as possible, tiny. On the other hand, there are some obvious considerations coming from the fact that I want to simulate a system for a given amount of time. And the given amount of time is determined by the physics of my problem, right? I want to simulate, I don't know, I mean the vibrations of a molecule, I need to be able to simulate my system for at least the several periods of vibration of my molecule, right? So the physical problem I'm interested in will set the actual time scale I need to consider. So if I make delta t too small, the total number of steps that I have to, you know, iterate my trajectory will be, of course, inversely proportional to delta t. So if delta t is too small, right? Obviously the total trajectory, let me call it t, the total, you know, length of a trajectory will be, of course, the number of steps times delta t. So if delta t is too small, the number of steps will be too large. And so the calculation will be too long, it will take too long, right? So after these two conflicting criteria that I have to somehow try to accommodate. Well, it turns out that at least with the Verlet algorithm as a rule of thumb, if you consider the typical fluctuation of your system, okay? Which can be here, for example, represented by this roughly, right? Let me call it tau. This is a sort of characteristic time of my system. For example, if I have an oscillator, it will be the period of the system. This is more or less the time that describes the oscillations, the average oscillations of my system. Well, as a rule of thumb, in order to be able to integrate this properly, certainly my delta t must be much smaller than tau, obviously, right? Otherwise I won't be able to describe a fluctuation. But as a rule of thumb, it turns out that delta t can be, if I want, of the order of, I can make it something of the order of 30 times smaller. It can be 25, it can be 40, depending on how accurate I want this oscillation to be described. But the Verlet algorithm, in particular, is quite stable if you choose a delta t, which is only 120th or 130th of the actual period of oscillation of your system, okay? So in practice, your system will be a complicated system and may have different time scales. Take water, for example, right? Liquid water, there are several time scales. There is a time scale of diffusion, which is normally a fraction of picoseconds. And then there is a time scale, the fastest one is the time scale of the vibrations of the OH bond, for example. This is actually extremely fast. So there are different time scales. So how do you choose actually your delta t in a molecular dynamic simulation? Well, you need to integrate all your degrees of freedom. So if you want to do that, you need to be able with your delta t to describe the fastest ones. Once you describe the fastest ones, you, of course, be able to describe also the slowest one with the Verlet algorithm. So as a matter of fact, your delta t must be a fraction of the order of 20, 30 times smaller than the fastest oscillation of your characteristic time of your system, okay? You have to pick the fastest oscillation time in the case of a water molecule, for example, that would be the stretching of the OH bond. And then that defines your minimal time scale that you have to be able to reproduce in terms of oscillations. And then your delta t will be correspondingly set, right? And in some sense, those of you might be familiar with the simulation of water. You might now understand why a number of water molecules, water models assume that the molecule is a rigid object. 90% of the molecular models for water using molecular dynamics assume that the water molecules is a rigid object, right? Of course, by doing that, you're introducing a lot of approximations because you're not looking at the dynamics of the water molecules. But in practice, you immediately see that this gives you a huge advantage in terms of delta t because you can actually now, you can forget about the faster oscillations of the stretching of the bonds and you can only integrate, you only need to integrate the dynamics of the relative motion of the water molecules and rotations which is much slower. And therefore, this allows you to use a delta t which is much longer than the one you would have to use if you wanted to integrate the short one. Okay, good. So let me now move on and let me again share the screen and show you, let me see. Okay, right. So this is a standard protocol of a molecular dynamics simulation. You first initialize your atomic positions, you need to give some initial positions and velocities because you remember, I mean, in the algorithm, you need to know, well, either the velocities or you have to give the positions at time t and time zero and time delta t which is sort of equivalent to given velocities. That initializes your simulation and then you start running your simulation. Now, of course, if you want to study your system at the given temperature, you have to choose atomic positions and velocities that are as much as possible close to what you think is going to be thermal equilibrium in your system, right? For example, it doesn't make sense to start, if you want to study a solid crystal at, I don't know, room temperature, it doesn't make sense to start from the crystalline positions because this is not a good representation of your system at finite temperature. At finite temperature, the atoms will be already displaced from their position. Similarly, the initial velocities, you can actually choose them. This is simpler. This is simpler. I mean, choosing velocities, you can actually do it quite straightforward way because you can just sample the Boltzmann distribution. So you can choose the initial velocities that are somehow chosen to represent, to describe the Boltzmann distribution, okay?  And how much do you need to weigh? Well, if you are interested in reaching thermal equilibrium, there are, I mean, different goals that you can give yourself. I mean, if you're interested in looking at valuation of frequencies of a molecule, for example, you're not interested in reaching thermal equilibrium. You just want to know the period of the oscillation of the molecule. But in most cases, you're interested normally in study systems of thermal equilibrium. So what you do is you let your system evolve until the system reaches some sort of thermal equilibrium. And certainly you want the system to lose memory of what was the initial configuration because that initial configuration was chosen, you know, out of a guess. So you don't want that guess to affect the final result of your calculation, okay? And then once you reach the equilibrium, you start accumulating averages of the observables you're interested in. Here, this is just the generic definition of an observable. Let me give you one example. I mean, it could be the temperature, for example, the instantaneous temperature of your system. Temperature in the specific case of this equation, A for temperature would be just P squared, the sum of all P squared, the velocity squared. But you can also calculate other averages if you want. And you can average it over the length of the simulation. Now, it's interesting at this point to make somehow a connection with statistical mechanics. I think this is an interesting consideration that perhaps there's no time to go into any detail now, but you're also, you're of course familiar that in the realm of statistical physics, averages are calculated in another different way, right? They're calculated as, I mean, as averages over the Boltzmann factor, right? Using standard techniques in statistical physics. So the assumption here is that the molecular, in fact, it's rather an assumption on the side of statistical mechanics. We know that the real system behaves like a Newton, follows Newton dynamics. But if you average your system over long times, the assumption is that this averaging will eventually be closer and closer to the averaging, the perfect averaging that you would achieve if you were doing perfect statistical sampling of your observable. Okay, so there is an hypothesis here. And there is a, those of your interest, there is also the assumption, I mean, the ergodic hypothesis hidden behind this equation between the two averages. Because in writing the right hand side, you're assuming that the system is free to sample all possible configurations that are available in phase space. In molecular dynamics, you never know because you might be stuck in some, you know, local minima, for example. And you don't know whether another minimum exists as long as you're exploring with your dynamics only a certain minimum, okay? I just, I mean, this is not, I don't want to enter into this, there's of course a huge, you know, literature on this. And I don't want, I just wanted to give you some hints that this is sometimes a problematic aspect that you need to consider. So let me again show you in practice what happens. This is a simulation in which the goal is to simulate your system at room temperature, right? 293 Kelvin. So you start from different initial conditions.  But as you can see here, your initial guess is never perfect, right? It's a guess. So it is a fact that if you've evolved molecular dynamics simply using Newton's equations, you're actually, this is now the instantaneous temperature that you calculate in your system. So let me just show it here one second. Sorry, sorry, time here to erase the, okay, so in this particular case, what I'm calculating was a instantaneous temperature which is essentially, whoof, now I have to remember the factors now. Well, there's certainly a sum over the particles of the, they are kinetic energies, the I squared, right? So this is, I know this. I mean, the equipartition theorem, this is three over two N number of particles KBT, right? This is from equipartition theory. So the instantaneous temperature is now one over two over three number of particles KB, okay? So if I calculate my instantaneous velocities here and I compute this object, I divide by the number of particles and by KB, what I get is, I can define it this way, right? And notice this value depends on the velocities at any given time. So it's not a constant, it will actually vary in time. So let me go back now to this curves. These are actually curves that show you the instantaneous temperature along different runs where essentially the difference was only the choice of the initial conditions. By the way, the choice of the velocities, as I said, is typically straightforward. You just sample the Boltzmann distribution at the temperature you want. So here it was 293. But in terms of positions, there's no way you can guess with exactly what would be the positions that sample the E to the minus beta V, the potential part. So you can only guess it, right? You can try to displace the particles a little bit and try to guess what would be the equilibrium values. I mean, the value of the displacements which gives you the final temperature that you want. So you can see different choices of the initial guess lead actually to different values of the average temperature that the system reaches after quite some time, after the equilibration. What you see here is also the fact that it takes time, normally of the order of 10,000 steps in this particular case to reach what we can call equilibrium. At the beginning, everything will fluctuate because you started from a configuration which is not exactly representative with your statistical system at the given temperature. But after a while, the system will reach a state which is stable in some sense where the average is calculated and the temperatures will be different. So this of course raises an important question and this is how can I actually simulate my system? Oops, sorry, stop sharing here. How can I actually simulate my system at the given temperature? You don't want to do a tie-in and errors. You want to simulate your system at room temperature. You don't want to start with different guesses and then choose the one that gets closer to room temperature. Now comes the concept of simulating doing molecular dynamics in different ensembles. Okay, so far I've essentially described molecular dynamics as a solution of Newton's equations. Now, Newton's equations as you know, conserve the energy, the total energy of the system that is kinetic energy plus, sorry, I'm using here total energy in a different meaning than in a standard of the initial. Here for me, total energy is the sum of the kinetic energy and the potential energy, right? In classical mechanics, if you evolve the system with Newton's equations, that energy is conserved, right? So we call these simulations, we call them micro-canonical simulations, okay? So a simulation where energy is conserved, we call them micro-canonical, micro-canonical, okay? But sometimes what you want is simulation in which instead temperature is, I wouldn't say conserved, but I would say where the simulation reaches samples, the macro, the canonical ensemble, okay? I'm using here conserved, but it's not actually conserved because temperature will always, what I mean by that is a simulation in which I'm actually sampling the canonical ensemble at a given T that I can give at the beginning of the simulation, okay? At given T. This is not what we did a minute ago. We did this and we didn't know what the final temperature was. Now we want to do something else. We want to do a simulation, a canonical simulation at the target temperature that I impose to my system. Now the trick to do that, and there are several tricks to do that. I'm just showing one because there are plenty of tricks to do molecular dynamics in the canonical ensemble. Let me just describe very qualitative terms how one goes about doing that, right? You have your equation, which is the Newton's equation, but what you do is in addition to your force, you add the friction term, right? The friction term is proportional to the velocity of your particle, right? And you modulate your friction in such a way that you reach exactly the same temperature the temperature you want, okay? So if your instantaneous temperature is above the target temperature, then you switch on the friction. Vice versa, if the temperature is below the target temperature, you actually change the sign to this friction and you accelerate your particles, right?  You give an acceleration which is proportional to the velocity. So the trick here is how can I modulate gamma in such a way that I reach a target temperature? Now, it turns out, and there's of course a plenty of literature on this and I can point you to the correct literature, is that the best way to do it is to actually do it in a slightly different ways than I described it here. And by saying that this actual gamma is the derivative of an object which evolves with this equation essentially, where this is the difference between the instantaneous temperature and the target temperature, okay? So now this is a slight variant of what I described before, but it still goes in the right direction because if the instantaneous temperature is above the target temperature, I want gamma to be positive, to increase, right? Because I want the friction to start acting on my velocities, on my acceleration here. Vice versa, if the instantaneous temperature is below the target temperature, this becomes negative and gamma gets reduced until it becomes negative and accelerates the particles, right? Now, it turns out that this mathematical framework is an extremely nice mathematical framework that you can show it samples exactly the canonical distribution at the given temperature. And there is actually a nice mathematical proof that was derived by Nose in the, I believe in the 70s, probably early 80s or late 70s, okay? So this sort of algorithm actually goes under the name of the Nose thermostat. So these are called thermostats in general because these are algorithms that allow you to sample the canonical distribution while you're doing a molecular dynamics. Now, we have to be careful.  It's nothing to do with the physics of your system. It's something that you're introducing to make your calculation going, sampling the canonical ensemble, okay? So if you do this, you'll be able to reach the canonical ensemble, but you forget about extracting dynamical information from your trajectories, right? So the frequency will be completely, well, completely, there will be slightly different with respect to the real frequencies of your system. The diffusion coefficients will not be the real ones, okay? So in practice, what you do is you first start with this framework with the canonical simulation.  Hopefully the temperature in the micro canonical ensemble will be very close to the temperature. If the system was equilibrated, was at equilibrium, will be very, very close to the temperature you fixed in the simulation done with the canonical ensemble, okay? This is thermostats. And of course, there are many different applications of this concept of thermostats. You can actually extend this, not only to temperature, but you can extend it also to pressure, for example. You can this, I mean micro canonical simulations are, the cell is fixed, everything is fixed, energy is conserved. There's no chance the system can say, change the volume or change the pressure. Well, actually the pressure will fluctuate, but the volume will essentially remain fixed. So there are techniques which essentially are all inspired by this concept of thermostats that you can apply to sample other thermodynamical ensembles. For example, the constant pressure ensemble. I'll show you an example right now of what the extreme again. Okay, let me see. Yes, this one, for example. This is something that was introduced by Kelly Parinello and Ysraman in the early 80s. And you see the similarity between what I showed before at the blackboard, at the temperature and here now the pressure. So instead of considering now a friction part that you consider the cell as varying in time. So H here is the vectors, actually the matrix that defines the three vectors that define the simulation box. And you're giving now the freedom to the simulation cell to vary in shape, you see. And with a dynamics, and what you see at the bottom is the dynamics of the cell now of the H matrix, which is driven by the unbalance between the instantaneously calculated stress tensor, okay, the equivalent of the instantaneous temperature in constant temperature simulations and the pressure, the target pressure that you want to achieve. In other words, if the system will reach an equilibrium only when and actually the shape of the cell will reach equilibrium, only when the internal stress is equal to the target pressure, okay? So this is a mechanism to drive your simulation towards the target values of the, the thermodynamic variables you want to fix, which is temperature in the previous case and it is pressure in this particular case, okay? Let me stop sharing here. Okay, now I just wanted to, you know, give you one example of quantities that you may wish to calculate out of a molecular dynamic simulation. A temperature is one, of course, and I just showed it before. You may wish to compute the instantaneous temperature and to monitor the instantaneous temperature while you're evolving your trajectory with time. And this is actually a very good way to determine whether you've reached the equilibrium because the fluctuations of the temperature tend to remain, say, constant around the given value which is the equilibrium value of the temperature. By the way, you will never get rid of fluctuations in your simulations, right? Even a instantaneous temperature will always vary in time. What you need to make sure is that its average remains constant. This is equilibrium, but the instantaneous value of the temperature will, of course, always fluctuate. And the fluctuations will actually be, can be shown, will be proportional to the inverse square root of the number of particles of your system, okay? So the oscillations of the instantaneous temperature. This is actually quite clear if you think about the single particle, right? Think about just a single harmonic oscillator. In a single harmonic oscillator, there's essentially only one type of dynamics at this one, right? And the kinetic energy, the average kinetic energy of your system will go from a maximum value when the particle is at the bottom of the potential to a zero when the particle is at the top of the potential, right? If it oscillates like that. So you see that with one particle, temperature, temperature, whatever is the meaning of temperature for one particle, which of course, you have to be very careful, right? Because temperature is only defined in the thermodynamic limit, of course. But I mean, let's imagine that you do a simulation with a single particle, an harmonic oscillator, you immediately see that the value, at least of the kinetic energy, let's not even call it temperature, oscillates between zero and the maximum value forever, right? So in some sense, the fluctuation of this temperature is of the same order as the average value, which is half of it, okay? So you see the fluctuation is essentially equal to the instantaneous value of the, so to the average value of the temperature for one particle. If your system is composed of many particles, these fluctuations will reduce, reduce, reduce, and in the thermodynamic limit, you will essentially have a constant temperature for the average and also for the instantaneous temperature. But it's case like one over square root of the number of particles in your system. Okay, let me give you another example of observable that you may wish to calculate. I'm not sure whether this will be considered in the tutorial, but it's actually a very common exercise in especially when you study liquids. And this is the so-called pair correlation function. So you look at, if you have a system of particles, I mentioned this because this is extremely common type of observable. You have a liquid, let me say you have a liquid, and well, you sit on a given particle now. You sit on this particle, for example, and you look at all the distances that you have between this particle and all the other particles, right? Like this. And you construct a probability distribution that at the given distance R from a given particle, okay, from this, you have another particle, okay?  But if you're average over several steps, what you will obtain is something typically of the type of this type, okay? And this is quite universal. I'm showing this because this is another universal behavior for essentially all systems. What this is showing is that, for example, there is essentially zero probability that two particles come in contact, right? Well, if in contact, you mean really coming in contact at zero distance, right? You have always repulsion at some point between two particles. You see more probability for particles to be at the given distance, right? This distance is normally the distance that particles take in the solid, nearest neighbors in the solid, right? In the liquid, there is always a reminiscence of the fact that you have a shell, that you have some type of bonding regardless of the system you're studying. And so you always have a maximum of the probability at the distance, which is, you know, of the order of a few Armstrong's, which is the typical distance of two particles in when they approach each other. And then you have fluctuations, which eventually tend to a constant. The constant of course comes from the fact that you're feeling the disorder system. If you look sufficiently far away from yourself, you're essentially going to see an average distribution which corresponds to the density of your system. It won't be able to recognize fluctuations, systematic fluctuations due to the fact that there is a memory of yourself in that particular position, okay? So this distribution normally is then to a constant, to one, to the density of your system, okay? So this is called pair correlation function, and it's one of the most popular pair correlation function. It takes also other names, okay? It's called geobus. Now, I don't know how we're doing with time, so it's 9.15, maybe, I think it's probably, I think it's probably, I mean, I don't like to talk for one hour and a half, so is it okay if we take, say, a five-minute break so you can go to the toilet if you need, you can get a glass of water or whatever? So can I suggest that we stop here briefly for five minutes and we reconvene in five minutes, which for me is 9.30, let's do 9.35, okay? Yeah, fine, okay. Sorry, sorry, not 9.35. Yeah, it's fine, why not, yeah. 9.35, okay? Yes. Good, okay, so what do I do? I'll write it here. Reconvene at 9.35 a.m., of course, Central European time. So in five minutes, approximately six minutes. Central European summertime. Okay, whatever. It's GMT plus two, okay? See you in a moment. Okay, I think we should start now. Can I start? Yeah, sure. Okay, good. So first, I guess there was a question on the chat about the difference between Nose and Nose Hoover. I've been using the word Nose just because of historical reasons. Nose was the one who actually came out with the proof that the thermostat samples, the canonical ensemble. But in fact, the scheme as I presented is actually called now as Nose Hoover. Hoover gave some important contributions in the definition of the algorithm. So it's in fact, there are several different flavors of this thermostats. I've only shown this qualitative, very sort of poor main description of the thermostat. But you're right, it's called actually in practice, it's called Nose Hoover. Okay, so this is the second part now of my lecture. And this is now going to focus on the potential, right? We are in a Venetian community or a place where a community that wants to learn how to use a Venetian method. So this is the place where a Venetian comes into play now. Let me just give you a brief historical, excursions of how the field of potentials, force fields, inter-atomic potentials, was developed over the years. The beginning of molecular dynamics was essentially just done what we call classical molecular dynamics. That is, it was the first calculations were done with inter-atomic potentials, which were written as some of pair interactions in the 60s at least. And then in the 70s, people started to introduce some environmental variables into the definition of the potential. It was only in the 90s, I would say, with the breakthrough made by Karen Parinello, that the people have started to use the abinitial potential, concretely in molecular dynamics. And it's the last 20 years have been something that in fact, slightly goes back to the beginning in the sense that the people are now constructing inter-atomic potentials based on the abinitial data to what we call abinitial-based molecular dynamics. But let me just, I mean, briefly go back to the beginnings because this is very instructive also from a number of points of view. So what is the concept? The concept is that suppose you don't know anything about the potential. Well, the first thing that comes to your mind is to say that the energy of interaction between two particles, well, to first approximation, perhaps, it can be written as a sum of pair interactions that is atom I interacts with atom J with a given potential. In fact, if I can now stop sharing and I'll go back to my, if you take any pair of atom, okay, any pair of atom in the periodic table and you look at the energy of interaction between two atoms, well, it's quite, one can come out with a quite universal description of this interaction, which certainly has a repulsive wall at short distances, atoms, any pair of atom beyond the given distance will start to repel each other for a number of reasons, power repulsion, nuclear repulsion and so on and so forth. So there is always a steep repulsion wall at very short distances, much shorter than chemical bonds, for example. Then at long distances, we also know that at least neutral systems, they all, that's actually an exact statement and two neutral systems, they always attract with something that goes like one over D to the sixth, right? This is the famous Vander Waals coefficient. In fact, it's the long distance interaction between any pair of neutral objects. It goes like this, like minus one over D to the sixth. So these are quite universal statements, right? Any pair of atoms would do that. Now, if you just join these lines analytically, you will have to realize that there will be a point where there is an equilibrium, right, at some distance. And of course, this corresponds to a given, so energy of equilibrium between any pair of particles. So again, this is any pair of particles, any pair of atoms in the periodic table, with probably one or two exceptions. Okay, where the difference is coming to play is when, of course, if you choose different atoms, is what is the distance and what is the depth of this well? I'm sure you're familiar with the fact that covalent bonds, atoms that interact with the covalent bond, like carbon, for example, the distance is very short, okay? The distance in carbon is about one point, something Armstrong's. And the depth of this energy is also quite large. It's of the order of electron volts in the case of chemical bonds, right? So this universal curve has, it's short and it's very deep in the case of covalent interactions. But take helium, for example, right?  So the depth of this potential is less than millivolts, and in fact, even less depending on the rare gases. And the average equilibrium distance can be several Armstrong's, 10, 20, 50 Armstrong's, the equilibrium distance between two pair of completely neutral and, I mean, like rare gases. So while this universal curve in qualitatively is the same, the actual quantitative aspects of this can be quite different across different atoms in the field. So let me go back now to my slide here. So if I assume that the two particles in my system interact with the pair potential, I write the sum of the energy. I mean, the total energy has a sum of pair interactions, making sure I don't count interaction twice. And then, and then of course the next step if I want to improve this potential is to say, well, how is the interaction between I and J affected by the presence of atom K, for example, right? Clearly, if I take atom K in this figure here, and I move K, say from the position you see it now to the position that is exactly in the center between the connecting line between I and J, you will immediately realize that the interaction between I and J is going to be affected. Suppose it was a covalent bond, for example, between I and J. If I place an atom in the middle of a covalent bond, of course the interaction between I and J is going to be affected. So you clearly see that the pair interaction approximation is really in approximation. And in principle, you need to correct it for the presence of third atom. So this is what people scientists did in the 70s. They started to construct the potentials that were essentially correcting the pair interactions by including three body terms that were supposed to correct the pair interactions for the effect of a third atom K in this particular case, depending on the position, of course, of atom K with respect to atom I and J. Well, you immediately realize that that's never ending story, right? Because even if you now correct it with the presence of atom over third atom, you can say, well, but the interaction between I, K and J is going to be affected by atom, whatever, L, right? So there's going to be a fourth atom which is going to affect, depending on where the fourth atom is, the interaction between I, K and J. So you must, in principle, include also a term that corrects the third order contribution with the fourth order contribution, four bodies, five bodies, six bodies. In fact, it has been shown that this is not even in some cases, it's not even converging as an expansion, okay? So people have come up with quite different approaches in the later on, especially in the 80s. And especially, people have introduced the terms that somehow consider the fact that the two atoms, when two atoms interact, you have to consider in an average way the position of all the other atoms. This is the underlying concept, underlying the embedded atom models. You're embedding your pair interactions into something, into an average field which is the contribution, which essentially gives the contribution of all the other particles to the two-body interaction between these two specific atoms. I don't want to go into any detail about this. I mean, there's plenty of literature about how people have constructed empirical potentials for that. And of course, you have several parameters that you need to fix in deriving your empirical potential. Historically, this was done, especially in the 70s and 80s, by fixing those parameters to experimental values, such as, for example, melting temperatures, equilibrium, lattice spacings, dynamical properties, and so on and so forth. But let me now come to the ab initio potential and let me actually stop sharing this because I'm used to the black part now. So the breakthrough came when, I wouldn't say people realized because that was obvious to everybody, but when people started to become, when it started to become possible to say that the potential could be evaluated entirely ab initio instead of, right? So what we want to do in our new, in our molecular dynamics, let me remind you is that we want to solve Newton's equations and we want to calculate the force, but we know that the force is essentially the derivative with respect to the coordinate of a particle of some sort of energy, which is not only a function of our eye, of course, but it's a function of all the particles in my system. Oops, sorry, I need to do this. R1, R2, blah, blah, blah, blah. RI, the particle I'm interested in, and RN, the final, I mean, if I have N particles in my system, okay? This is a very important concept, okay? So the energy of my system is, of course, a function of the position of all the particles in my system. If I move one, the total energy will change, right? And when I calculate the force, I need to calculate the derivative with respect to one of these variables. And this will give you the force acting on the eye. As I've seen before, as you've seen before historically, people have described this as a sum of pair interactions or some embedded scheme or so on and so forth. But the truth is that this one was able to evaluate this total energy of the system for a given configuration of the atoms ab initio, that will be the solution to the problem. And in fact, this is what we do when we do ab initio molecular dynamics. We're actually calculating the total energy of the system. We're evaluating the force as you did in nothing last week. And then we use the force that was calculated with ab initio approximation with ab initio accuracy to evolve the molecular dynamics of my system, okay? Now let me somehow elaborate a little bit on the complexity of this task, right? So I have my box here and I have my particles here. Now I'm not going to call them atoms any longer, okay? I'm going to call them ions or nuclei because in an ab initio framework, the electrons are everywhere, right? I have the particles here. So this is, for example, a nucleus with charge z i and mass m i, okay? And position r i, of course. So the points will be defined by, will be our, what we used to call atomic positions, okay? But then inside this box, I have electrons everywhere, right? Well, let me just to be a little bit more, let me assume that electrons are somehow localized close to the atoms, for example, okay? So these are electronic distributions. This could be, for example, a typical electronic distribution in a system of particles where the particles are placed, the nuclear are placed in these positions, okay? So in order to calculate this energy here, I need to solve my quantum mechanical problem for a given configuration of the atoms. So this of course means that every time, so that every time I solve my quantum mechanical problem, I have to find what is essentially the wave function of the electrons. So this is now the electronic wave function and it's the ground state wave function. And this is going to be a function of all the electronic coordinates or only one if we are in density functional theory. If you think in the many body framework, this is actually the collection of all the position of the electrons. And notice that I'm using now small letters to describe positions of electrons, okay? I'll be consistently using, as I said, capital letters for ions, atoms or nuclei and small letters for electrons, okay? So for example, if I'm describing this particular wave function, the orange wave function inside this simulation box, I will describe it with, well, let me use a collective index for all the electronic particles. But then I have to remember that this is parametrically dependent also on all the RI positions, right? Because this wave function will be different depending on where the positions of the nuclei are, okay? So it's a wave function expressed in terms of, as a function in terms of the electronic positions, one or N, depending on whether many body are being filled, but also parametrically on the position of the nuclei. What do I mean by that? What I mean by that is if I change the configuration of the nuclei, I displaced them, this will be different. And of course the wave function will be different because for every given configuration of the nuclei, I will obtain a different wave function. And what I know of course is that this wave function is the solution of the Hamiltonian problem. And the Hamiltonian is of course, a parametric function of the coordinates, okay? So this wave function will be the ground state, well, let me put the ground state here just to be clear, of the Hamiltonian problem. And the Hamiltonian of course is defined for a given configuration of the nuclei because it depends on where the nuclei are. I mean, the Coulomb potential or the pseudo potential will be placed at the position of the nuclei. So this Hamiltonian depends parametrically on the position of the nuclei. Its ground state for the electrons will be given by the solution of this equation. So as a consequence, the ground state will depend parametrically on the position of the nuclei. I hope this is clear, right? So I'm solving the quantum mechanics of the electrons. You have the initial problem. And I'm solving the quantum mechanics of the electron parametrically. That is by changing every time the parameters in this Hamiltonian, which is where the atoms are. Let me actually show it to you a bit more with real equations so that you can, here you go, okay? So this is essentially what I just explained. The box at the bottom right of this slide gives you the parametric Hamiltonian, which depends parametrically on the position of the nuclei. And the energy, the total energy of the system is the ground state energy of electronic problem for a given choice of the nuclear position. Now, when I now do a binational molecular dynamics, what I have to remember is that every time I solve this problem, I need to solve it once, I will have to calculate forces. And I'm normally using Helman-Feinman theorem to calculate forces. I don't need to go through Helman-Feinman theorem. I guess I think you already saw that last week. So there are ways, efficient ways to calculate from the electronic ground state, its derivatives with respect to the coordinates. And then I need to evolve. I remember what we discussed at the beginning of the lecture. I need to involve my particles, right? So the next steps, if I now do a one molecular dynamic step, the nuclei will change a little bit, right? It will be displaced. This one will go here. This one will go here. This one will go here. So there will be small displacements where small depends of course on delta T, delta T I chose. But it will be slightly displaced. Now, the next step of the molecular dynamics, I need to solve again the electronic problem, right? I need to go back here, change the Hamiltonian, find the ground state again, determine the ground state, determine the ground state energy and recalculate the force, okay? So I hope you get the flavor of why doing that initial molecular dynamics is so complex and so time consuming because every time in principle you move the particles, you have to solve the electronic problem again. You have to find the ground state again of your system. Ground state energy, determine the forces and then evolve delta T, delta T, delta T, delta T, okay? So you easily imagine that the complexity of this task is much, much higher than, for example, determining the forces based on a pair potential, right? The pair potential in nanoseconds any, you know, with any computer you can calculate it for an arbitrary large number of particles. Here, of course, if your system has, I don't know, you know, a thousand atoms, you're simulating a system of a thousand atoms that means calculating the ground state energy and forces of in a simulation box of 10,000, of a thousand atoms, right? So this may take actually quite a long time. So this is why abinition molecular dynamics, of course, is limited, much more limited than classical molecular dynamics in terms of sizes of the systems and length of the simulations. Instead of doing, as I said at the beginning, you know, micro, even milliseconds in the case of classical molecular dynamics with abinition molecular dynamics, you rarely go above one nanosecond of total dynamics. In fact, most simulations are of the order of 10 to 100 the picoseconds in terms of total duration of the molecular dynamics. And even in terms of sizes, you rarely go above 1,000 particles in your abinition simulation because that has to become quite consuming, time consuming and resources consuming in practice, okay? Now, having said this, however, I'm sure you're already seeing that we can take advantage of some facts here.  Yeah, that's right. But if your particles move only very little, that is, if this distance is small, and it is small because we are evolving the particles with a small delta T in a continuous dynamics, then of course, the solution of this problem is going to be enormously simplified by the fact that I already know the solution of the problem at the previous step, right? In fact, as a guess, the wave function of the system, the entire system at the previous step is an excellent starting point for the solution of the ground state problem at the next step. It's much better than starting from random orbitals of some, even some educated guess of what could be the ground state, okay? So binational dynamics takes a lot of advantage from the fact that you already know the solution at the previous step. And so when you evolve it, you can take advantage of that solution, okay? To speed up, and there are various algorithms that allow you to do that in an efficient way. Let me rephrase this problem in a slightly different sort of mathematical framework because this is very useful to now understand what is a carparinello molecular dynamics. So let me, oops, okay, what time I have? Yes, a few minutes, almost done. Okay, let me now show you this. This is now the energy of my system, the ground state energy of my system as a function of the, sorry, the energy of my system as expressed as psi H of R, psi as a function of psi. So what am I doing here? I'm fixing R i, I'm fixing the configuration and I'm simply imagining what's happening to the energy that is to the expectation value of the Hamiltonian. If I vary the wave function, okay? So I'm here in the realm of quantum mechanics, if you wish. There will be a ground state, okay? And this expectation value we know already that it has some sort of minimum here, right? I mean, that's the variational principle of quantum mechanics. If I evaluate the energy as an expectation value of the Hamiltonian on an arbitrary wave function, this arbitrary, this energy will always have to be variationally higher than the value I obtained for the energy at the ground state, okay? But this is my ground state. And I mean, this is just very qualitative description of what, mind you, I'm using only one coordinate, but this is a Hilbert space, right? Because the wave functions are a Hilbert space. So you have to imagine this as an infinite dimensions in not just as one, I'm just collapsing the wave function as if it was a one dimensional object, but it's an Hilbert space. But there is a minimum, right? At the corresponding to the ground state in quantum mechanics, if you fix the Hamiltonian and now I'm fixing the Hamiltonian, this corresponds to a minimum. And of course, I mean, I'm just rephrasing concept that perhaps are very obvious to you, right? That quantum mechanics finding the ground state is equivalent in a variational picture to identifying the minimum of this object in the space in the Hilbert space of wave functions. And this is the ground state. Of course, a ground state which will be determined by whatever is the Hamiltonian I've been using here. I'm using the Hamiltonian for a given value of the nuclear positions. And of course, there are several ways to find the minimum of a quantum mechanic. I mean, variationally, I can do diagonalization in a different framework, but I can think of sort of seeing this as an optimization problem, as a minimization problem in which, for example, if I start from a guess of the wave function here, let me assume that the goal now is to find the ground state. I can start from a guess, a wave function which is a guess here, and I can try to identify where the minimum is. What do I do? Well, I need to know the derivative of this, right? If I want to determine where I need to move. But the derivative of this with respect to psi, that is the derivative of the... This is just quantum mechanics. It's of course, it's a functional derivative, right? So this is, forget about the factors, right? I'm really using, you know, hands waving arguments. Is H psi, right? The functional derivative of this object with respect to psi is H psi. Now it depends on... I mean, forget about the details, but essentially qualitatively, you can say that the force, right? This is the force in principle. I mean, if you see it from a classical point of view, the derivative of your energy with respect to the variable you want to optimize, which is the wave function, is H psi. Which is this force acting, let me use force, but it's not really a force because we're acting, we're in the quantum space now. We are seeing the quantum problem as a classical problem here, right? Optimization problem. We have a force which pushes the electrons towards the ground state in this direction. And I can come up with many algorithms. Once I know the force, I know the gradient. I can use conjugate gradients. I can use deepest descent, for example.  Well, then I need to evolve my wave function from the guess following the gradient. Sorry, this is minus the gradient, right? The gradient goes up. This is minus the gradient, okay? The one that pushes me towards the minimum, of course, I have to go again. Sorry, sorry, sorry. Yes, yeah, so the force is minus the gradient. So if I imagine an algorithm which evolves with time, now it's a fictitious time because I'm just finding the ground state of a quantum problem. And I'm introducing some sort of fictitious, I mean an algorithm time, which says my velocity of this particle, this object that I'm optimizing, goes with the force acting on it, right? This is called steepest descent.  And then in several steps, if you do discretization, you will eventually find the minimum following this gradient. Now, this is, of course, there are much, much more clever ways to optimize, to find the minimum of a quantum problem. You can diagonalize, you can do lunches, you can do conjugate gradients. I mean, quantum espresso has several ways to actually find the minimum.  Because of the following.  Suppose I now do, instead of doing this, I do this. Again, let me remind you that we are not in the classical space here. We are working with Hilbert spaces, okay? So I'm just imagining some infinitely dimensional field that evolves with this second order equation following what I call the force, which is the gradient of a potential along this least trajectory. If I do this, what will I obtain? Well, something a little bit different, right? Because this is now second order dynamics in an algorithm time. And so my particle or wave function, depending on how I want to call it, will start, well, it will start oscillating back and forth from here, boom, boom, boom, boom, boom, right? So imagine this in a quantum space, right? This is the ground state. I'm just letting my wave function evolve around the minimum in this infinitely dimension Hilbert space, okay? This is what I would get. Sounds silly, right? It sounds quite silly. Why should I do that? Why would I like to do that? Well, there is one reason I may want to do that. It's that on average, if you think about it, if I evolve the wave function in this, on average, the wave function will be very close to the ground state on average. It will be oscillating, but on average it will be very close to the ground state. So if I do this and I evaluate the average value of the wave function, that's not going to be very far from the actual ground state. Okay, we are ready now to move to Carparinello. And Carparinello can easily be understood if I now add the third dimension, which is now the position of the particles arrive. Once again, I'm collapsing everything into a single dimension, but it's a three-end dimensional space, right? That of the coordinates. In which sense? In the sense that once I found the ground state, and now remember our problem, our problem is to evolve the wave function finding the ground state every time we move the atoms, okay? If I now introduce the variable, the nuclear coordinates, the atomic coordinates, once I find the ground state, I need to now move the particles, right? So I will now be sitting at the position, which is slightly displaced in our eye, right? I will have to move in this out of the blackboard in a different space, everything will change. This will change, the Newtonian will change. So I will have to find another minimum, right? So well, let me now move it a little bit here. So this will be a slightly different parabola. It will be slightly displaced with respect to the original position. Its minimum will also be slightly different with respect to the original minimum, right? We'll be somewhere here, for example, now it's a little bit different, but in other words, this reflects a bit more at what I said before. If I now want to minimize the wave function for a slightly different choice of the nuclear coordinates, I need to, I mean, the ground state will be slightly different, but the original ground state, which was this one here, now projected this one will not be very far from the new ground state, which is this one here. So as you can see, I mean, I can see the problem of finding the minimum of the quantum problem when the particles evolve in their positions as a problem of identifying the minimum of something which looks like a parabola, which actually moves with time according to the position of the particles. So it's a bit like, you know, having a particle of which you find the minimum here or perhaps in a bowl, for example, and then you move the bowl and the bowl changes a little bit, the minimum changes, the position of the minimum changes and you want to continue to keep your bowl, your little particle at the bottom of this cup. So Carparinello comes into play when you say that you are actually evolving this and this is now the magic of Carparinello. Carparinello said, well, why don't we do this simultaneously? Why don't we evolve the Newton's equations for the particles and these fictitious, very strange equations for the electrons simultaneously? In other words, why do we need to find every time we move this parabola, always the ground state and let's stay actually out of it and let's start to oscillate around the minimum, okay? While this parabola moves in space, according to the evolution of the nuclear coordinates. So this is essentially the basic principle of Carparinello. Instead of looking every time for the minimum, which is very time consuming, because if you want to look for the minimum every time, every time you fix the value of the nuclear coordinates, you must find the minimum. So we have to solve a quantum mechanical minimization. Here instead, every step, you don't find the minimum. You just do H times psi. It's a simple matrix times vector operation. You don't have to diagonalize H. You don't have to find the minimum eigenvalue. You just do H times psi, which is numerically extremely faster than finding the ground state every time. Now if you do this, you know, if you let these oscillations here be fast enough, okay, the particle will adiabatically follow the minimum of your parabola. But this has to be fast enough to follow adiabatically the minimum of your parabola, okay? Let me now just, you know, summarize it in a year. Time is almost over. So let me see, no, sorry. This is non-file, and I skipped that, sorry. Here you go, here you go. And now I'm essentially summarizing what I just said here at the blackboard, okay? Perhaps a little bit better from a point of view, right? I'm evolving electrons with this fictitious second order dynamics, which keeps them, and you see now the blue line. This is the Carparinella evolution close, but never actually in the ground, in the instantaneous ground state of the particles. The red line instead is the so-called Borno-Penheimer molecular dynamics, okay? The molecular dynamics where every time you move the particles, you end up in the, you find the ground state of your system. Now, very briefly, can you tell the differences? Yes, there are some differences. And here, for example, on the right side, you see the, no, sorry, on the left side, you see the forces on the atoms calculated with the electrons in the ground state, and the forces calculated with the Carparinella approximation. They are not exactly the same, of course, because your wave function is not exactly the minimum. But what you can see on the right side is that the difference is actually an oscillating function. So on average, the force on the particles is quite similar to the force. And you can make this as small as you want if you say make the dynamics of the electrons faster. So, I think I will just stop here, oh, sorry, yes. I'll leave you with some basic texts about molecular dynamics, a little bit about the history, the basic papers, and also a very nice and comprehensive book that was written by Dominique Marx and Jörg Hütter about a definition of molecular dynamics where you'll find essentially everything you want to know. It's, I mean, it's a bit old, but it contains really everything you need to know to run a definition of molecular dynamics. And with this, I will stop here, and I'm not sure how it works here. Thank you very much. Thanks, Andrew. It works that normally we leave people from Zoom to raise their hand and make a question directly if they want. Meanwhile, I can post you on the chat some questions from the streaming, like one that is here now. Approach of molecular dynamics to set the system in phase transition water. Okay. All right, this is a very complicated question. I mean, now, starting phase transitions, unfortunately requires a huge effort in terms of statistical sampling, because you have to determine essentially the free energy of different phases. And determining free energies is not straightforward. It requires a lot of statistical sampling. So this is can be normally, you know, hardly done with ab initio molecular dynamics. So if your problem is to, you know, study specifically a given phase transition, then I would probably do ab initio molecular dynamics. But if your goal is to study, say, an entire phase diagram, including liquids, my suggestion would be to instead use, you know, classical potentials. And I didn't discuss it, of course. I just gave you a flavor at the very beginning, but you can now optimize classical potential based on the ab initio data. And you can do it also now with, you know, machine learning techniques. There are now a plethora of potentials being developed based on ab initio data and say machine learning, neural networks, and so on and so forth. So my suggestion would be, if you need a lot of statistical sampling to do it with potentials that are optimized based on the ab initio trajectories and ab initio data. Okay, there is another question from Slak that I can forward you here on Zoom. There you go. Okay, I go to the molecular dynamic calculation where you put on the table for consideration. Oops, it's crawling too fast. Oh, the question is the CPMD versus Borno-Penheimer. Yes, that's a very good question. So it has to do with your, essentially, just as a summary, it's a long discussion, but I can summarize it as follows. For insulators, it turns out that normally carpanella molecular dynamics works better than Borno-Penheimer molecular dynamics. For metals of resistance with small gaps, small electronic gaps, it's the opposite. You better do Borno-Penheimer molecular dynamics. It essentially has to do with the curvature of this parabola that I was showing at the end of my lecture. For an insulator that the curvature of that parabola is very high. So the parabola is very steep. And so it's easy to follow. Electrons have an easy job in following adiabatically the evolution of the nuclei. For metals, the parabola is very flat because actually the curvature of that parabola is the square root of the energy gap of your system. It's been shown. And so for metals, you better use Borno-Penheimer. It's a bit of a trade-off. It depends also on the efficiency of your code. Of course, if you have an extremely efficient code for Borno-Penheimer molecular dynamics, you better use Borno-Penheimer molecular dynamics all the way. One more. Thank you. There are two more that I already posted. Yes. OK. Then if I can just one last slide, which addresses that last question, if I can... Here, here, here, here. OK. These are all the things I would have loved to tell you, but I didn't have the time to do. So the first question is actually the one that was in the chat that addresses the first one. I'm assuming here that you can do classical molecular dynamics on the nuclei. This is true for 99% of the atoms, but there are exceptions. And helium is a very notable exception. Helium is a quantum fluid or extremely important quantum effects. So you have to be careful. Hydrogen, of course. The lighter is the particle, the less you... The more you have to be careful about doing quantum dynamics for the particles. And let me just take this opportunity to mention three more topics that for you to think about. One is how... I mean, I've been discussing essentially ground state dynamics, but it's extremely interesting to also look at the dynamics of particles when the electrons are in excited states. Things of photochemistry, for example. I didn't discuss how to accelerate sampling, in case you want to jump from different... I just alluded to it when I mentioned the Ergodic theorem. I mentioned in these things because they're all extremely fast expanding areas of research. So you might be interested in somehow looking at what's happening in these areas. And then only at the end I mentioned this concept of now building machine... Using machine learning, building classical potentials that is interatomic potentials, the way we used to do it in the 60s and 70s, but with parameters that are now determined based on a binational data. Okay, there are two, I think, one from the streaming and one from the slack, and then we can close. So one is how to do melt quenching. And the other one is... I posted it on Zoom if you want to answer. Where is it? Sorry. The first one is what? Sorry, I didn't understand the question. Is this one here? How to do melt quenching? I'm not sure what... I guess this refers to the... to how we produce amorphous systems. We normally produce them in simulations at least, but also in the real life by quenching melts. And I think the big question here is how to somehow approach the time scales that are required to produce correctly an amorphous system. I mean, ordinary glass is a typical example of a system that is obtained by slow, even on real-time scales dynamics. It takes minutes, if not seconds. I mean, I was born in Venice. I know very well how to produce arts, masterpieces with glass. And it takes minutes, in fact, even hours to obtain the final product. So this is clearly outside of the bridge for molecular dynamics, as a matter of fact, not to mention abinition molecular dynamics. So there are techniques that one can use to speed up the quenching of melts to obtain amorphous systems. Yes. Okay. I think we can take the last question from the Slack channel. Is it possible to perform molecular dynamics by using implicit methods which the solution can converge at any values of time step? Sorry, sorry, say it again. Is it possible? I can post it here, so you can read it. Yes, please. Go ahead. By using implicit methods, which solution can converge at any values of time step. I think this is what we call Born-Oppenheimer molecular dynamics. So every step, we find the ground state of the electronic wave function. I guess, I mean, the issue here is to which precision you want to obtain the ground state of your function. Because of course, numerically, you will never get exactly the ground state. And that difference can be important. In fact, one of the main problems of Born-Oppenheimer molecular dynamics is the fact that your solution, the error you're making in determining the ground state is normally biased in a given direction in wave function space. And that introduces a systematic error in the determination of the ground state and the determination of the ground state energy. Something that the Kalparanello method doesn't have because the Kalparanello method has the wave function oscillating back and forth. So you're sort of averaging out all possible errors you're making because of the fact that you're not exactly in the ground state. Okay, great. I think we can close it here. I would like to quote a comment from the stream saying a great lecture, Professor Skandolo, very clear and captivating. We relied the interactive aspect of the lecture. So thanks a lot. I think we reconvened 1035 with the Elson session. Thank you very much. I'm glad the system worked properly. Yes, yes, yes, it did work. I will say it did work. Thank you. Okay, thank you. See you later. Bye-bye. Bye-bye. Ricardo, did you say? Yes. Great. Yes, yes, yes. You talked to Stefano, didn't you? With Stefano, yes, yes, yes. What do you say? Shall we start? As you wish. We'll start in a minute. Okay, okay, fine. I'll share the screen. Yes. Okay, so we are gradually ready to start. So here we have today, Ricardo, who is the student of Stefano Baroni, I think in the last year of his program, I think, and driving this Elson session. So Ricardo has been also tutoring in the previous days. You might have met him in some breakout room. Anyway, now I leave to him the stage and having the Elson on a Benishi monoclonal. Thanks Ricardo. Okay. Thank you very much, Ivan. Actually I will be the tutor of the tomorrow Elson. So you will meet me tomorrow at the GPU Elson. So now we are going to see what the Carpernello code can do in Quantum Express. So just a small introduction. So the Carpernello code is able to do, to compute many different types of stuff. For example, the direct minimization of the total DFT energy that was the conjugative gradient that previously Sandra told you, then is able to do Born-Oppenheimer molecular dynamics that will be the first step that we are going to do later. Then of course the Carpernello molecular dynamics where the electronics degrees of freedom becomes some sort of a classical classical coordinates and they are integrated with the Verlet algorithm. Then it is possible, but we are not going to cover this in this tutorial to do variable cell calculation. That is, you let the cell to become a dynamical coordinates. So you can fix the pressure and simulate the system at a given pressure. Then the Noseuver thermostat that we are going to use to equilibrate the system at a fixed temperature. Then an important difference that this code can do only gamma point calculation. That means you do not have to converge the K point grid and you have to check of course for the size effect that there can be present. So you can start with a small system and then you can try the same simulation with a larger system and see if the quantities that you compute from the equilibrium simulation are the same or not. And then also other stuff that I don't say. So let's see what we are going to do today. So we will start with Born-Oppenheimer molecular dynamics. That means if what I am going to say and you calculate this with standard minimization of the DFT energy. Then we will use the ground state calculated within this Born-Oppenheimer simulation to start the Carparinello. We will see that to start the Carparinello we need the equivalent of the velocities of the atoms also for the wave function because wave function are classical quantities in this Carparinello formatism. So we need initial state that is the ground state and also the velocities of the wave function and we will see how to compute them when we initialize the simulation. Then we will see how to equilibrate the system at a given temperature. We need to set a temperature if we start with the Newton equation as you saw before the temperature is not given we have a fixed energy and the temperature is what we have. So we have to equilibrate with the Nozare-Hubert-Ternostat. After that we will see that the system is equilibrated at a temperature and we are able to run a longer microcanonical simulation that means fixed number of particles volume and energy whereby energy means kinetic energy of the ionics degrees of freedom plus the potential energy that in our case will be the potential energy calculated with the density function of theory. At the end of the day we will calculate the pair correlation function that was explained by professor Scandro before and also another dynamical quantities that is the mean square displacement that is able to tell us for example if the system is liquid or solid some atomic types are diffusing and other stuff like that. So for today we are running a very small system we will have a cubic cell with eight water molecules it is very small and very far from the thermodynamic limits but it is enough for our tutorial. So this is more a cap when I say trajectory I mean the time series of position and everything that is calculated at every time step. We are using of course like every day the bonapenibra approximation and we are saying that electrons are always in the ground state and we are not treating nucleus quantum particles and that's it. So let's begin with the bonapenibra as was explained before we have the Newton equations so we compute this with the Helman Feynman theorem and to do that we have to select all the parameters that you selected in the previous day that is the plane wave cutoff which sort of potential you have to use but we don't have to select the k points because we use only gamma points and then there would be additional parameters that are the time step in this case. So let's start by looking at the input file for our first step of our tutorial that is the bonapenibra molecular dynamics so you see this input file looks familiar, it is almost the same as the PW code you set the title of the simulation the type of calculation you have many different types of the calculation and now I'm reading the input description that you can find in the quantum express website and also in your virtual machine there is a shortcut in the Firefox browser so as you see many different calculations that are possible today we will cover only this one the other out of this tutorial then we start we have to tell the code that we are starting from scratch because we will see that later we are using also the restart files that the code will write in the file system so it is common with the company to start from previous previous runs with the same code so this after will be not from scratch but will be another thing then this number means a number directory write and basically it tells the code the name of the directory where it is going to write the restart files that we will read later with the next step then you have to tell of course how many steps I want in my molecular dynamic simulation here for example we are doing 100 steps then this variable is important since usually the time step is very small to compute the quantities that you want for your research for example let's say the fairy correlation function the g of r you don't need to save every every time step of the dynamics and if you do that you may end up with very large files that are useless because consecutive steps are very correlated there is not so much information in two consecutive steps so you can say I want to save one step every 10 so I have a lighter trajectory and enough information for compute all my quantities ok this is telling the code that I want to write a restart file every 1000 steps here we don't reach 1000 steps but suppose you are running a much longer simulation let's say that it runs for a day on a super computer and suppose at a certain point of the day the super computer crashes but as it can happen it happens to me a lot of time and what you do you don't want to lose the computation of all the days so you tell the code to write the restart file every 1000 steps that can be let's say an hour one hour whatever so you later you can restart from the last restart file then this is not necessary this test to write the forces in the trajectory files later we are going to see the trajectory files write the forces then the time step ok the time step is in units of atomic units and here there is a very delicate point this is taken from the documentation this is very important that the atomic units of the carparenello calls are different from the pw calls so there is a factor too so here in carparenello the time atomic unit is this value here in pw is twice that so 4.8 ok this is very important so here we set this to 5 it is a really small value if you think about it and it is useful at the beginning of the simulation to set a small value especially if you choose an initial state or initial atomic position that you guessed and you are not sure if they are good enough or not so you could for example end up with two molecules that are too near each other and if you choose a time step that is too big maybe this atom displays too much and you will not sample with enough accuracy the trajectory at the beginning is useful to use a very small time step for the bonopenheimer standard time step ok prefix is just the initial part of the name this is standard then we are using a cubic system with a size of 13 radius 24 alpha of two times hydrogen and oxygen ok the cutoff you have to converge as in the previous exercises here I choose a pretty low cutoff just to run faster the exercise but you have to check that the forces and the energy are converged within this cutoff so this is something that you have to check then ok now we tell what to do with the electronic degrees of freedom we tell the code that we want to do at every time step of the dynamics are conjugate gradient minimization of the electronic degrees of freedom here the codes for the first step we start from scratch so we set I think a random wave function or some educated west then at the next consecutive time steps it will use the previous solution for the ground state to start with and then it will minimize this initial guess that is a lot better than a random electronic wave function for example so you can see if you look at the output that the first step take very little time very long time for the first step because you start from scratch then from for all the other steps it will be much faster then here we say hey I want the Newton equation and I want to integrate them with the valet algorithm for the ionic degrees of freedom that means the r i all these quantities then how do I select the initial velocities of the atoms here I call the sample from the Maxwell-Boltzmann distribution that is something like a Gaussian where you have the square of the velocity of the atom times 2 k dt so you use this distribution that is not normalized of course you have the normalization factor in front you sample the velocities of the ions from this distribution with a temperature that is 600 kelvin then you have the rest of the input this is not necessary I don't want to move the cell the cell is fixed during all the simulation then the pseudo-potentials and all the atomic positions as in the pw code then you can see there is an additional ok let's go to the director maybe day 8 this is the first input files at the end of the files you have all the atomic positions all the atomic positions then at the end of the files you can have a special card that is specific to the caponello code that is able to tell the code to change some parameters of the simulation during the run this can be useful especially with the first step for the simulation because you don't have to write a lot of input files so for example here after 10 steps I use a bigger time step why? because suppose that I am in this situation where two atoms are too near each other after 10 steps I suppose I'm guessing it is not a rule or anything I'm supposing that the situation is fixed so I can use a bigger time step with the Born-Oppenheimer molecular dynamics maybe you can use also time step bigger than that ok this is typical with the Born-Oppenheimer then at the end I set again a small time step because this is the time step that I'm going to use later for the caponello simulation I will explain later why we are using a smaller time step with caponello simulation if you don't do this and you change the time step between restart different restart files that means you start you do the first input with let's say the time step equal to 20 atomic units then you stop the simulation you write the start file if you want to restart the simulation but with a time step for example of 5 atomic units in place of 20 you have to tell the code with an input that is this one you have to tell the code that you changed the time step why? because in the restart files there is no the velocity of the electrons of the ions is not saved you have only the position of two consecutive time steps so if the code does not know what was the old time step it cannot rescale the position of the time step before the last one to have let's say correct velocities in the algorithm so you have to tell this if you change the time step between different runs of the code ok so I think I saved all we can try to run this to run this is very simple ok this is here we have the input file so just check that you are connected to your cluster with maybe you have to run the VPN stuff I already run it I just I have to type run and here I run ok this one just a trick if you press ok it's going if you press Ctrl R you can look at the history of your bash for example here press Ctrl B plus R then you start typing and you can search in your history this is very useful when you want to keep up with this command so this should be very fast and let's see if it works remote remote school ok it's running so now we are running Bono Penheimer ok I see how do you select the initial position of the water molecules ok this is a very good question actually I guessed it in this particular cases I selected I put one molecule oriented randomly in the space and I used some symmetry group operations to to move them to move it around the simulation cell I do stuff like that so it is some sort of art to select the initial position usually I suggest to pick the crystal one if the crystal is simple or to maybe also start with classical molecular dynamics pick some equilibrated position of your system and we use this starting initial position there are also a lot of programs that help you to build the initial state for the positions so let's see if it's finished oh yes so now I copy back copy back from the cluster I think from HPC there is there are a lot of files you see in the list you see ok the trajectory files these those are all the files that the code generates at every time step then you have the restart files that are the files that we are going to read later to start the next step so let's have a look to this just have a look for the moment later we will plot it so you see there is here let's start with the cell file you see this is the prefix H2O the cell file is very simple is the cell vectors since the cell is constant they are always the same you see this number is the number of steps this is the time in picoseconds not atomic only these are picoseconds and then you have all the cell vectors you see the cell is simple we can see the position that are the most interesting one the position file ok you see the first line let's say the header for each time step is the same as the other file number of steps and time in picoseconds this is different time in picoseconds maybe I will write it time in picoseconds in H2O dot whatever files then you have all the atomic position in units of ball radius so the positions are in ball radius and you see now we don't see because we are not looking at graphically but the order of the position is the same as the input file so if you put the atoms in order you will find the same order note that I have written this question this sign here because this changed in version 6.6 before they were sorted by type and then by input file so this changed in the last version it is the most included way of writing them ok the velocities is the same but with velocities here time is in the same atomic unit as before so between velocities in carparinello and pwcode there is a factor 2 then there is the evp file that has some thermodynamic information on all the fun pieces for example number of step times this we will see later what it is it is the kinetic energy of the electrons but now since we are not running carparinello it is 0 temperature of the cell if you are moving the cell for example then the temperature of the ions that is something that we are really interested in then you find a thought that is the energy of the system the dft energy ok entropy a cons it is they write it later ok yes a cons is the potential energy plus the kinetic energy of the nuclei that is the energy that is conserved in a classical molecular dynamic simulation then a cons it is the constant of motion of the carparinello Lagrangian that now is not conserved because we are not running carparinello ok volume is constant pressure we are not computing the pressure so it is 0 you have to tell the code to compute the pressure if you want it so maybe we can see also what happened I prepared a little script to run to run some conversion tools that are not specific to quantum espresso and just are you able to see the result of the simulation and then I can open the trajectory that the script generated with a program just to visualize the atoms so what happens ok this is my system you have a different view of the same system top bottom right left everything you see our 8 water molecules and there is not much trajectory in here we have 9 frames so you see the atoms move really a little so that's it we are ready to go to the next step ok this is a slide just to remember you what what does the code write in the rest of the files so remember the rest of the files is specified by this n d w number directory writer that is the number appended here so this is 15 and the w and inside you find the all the data that the is the rest of the simulations ok now we move to the next point that is a little bit delicate so ok let's move to this before so to do molecular dynamics it is useful to express the system in terms of the ground that is an expression that together with the equations can define the equation of motion for the system so in particular this is the Lagrangian of the couple method ok my important points are these are the equations that we need to integrate with the algorithm so you see that there are some parameters but before the parameters let's focus on the equation itself it is a second order differential equations so to integrate a second order differential equations you need the initial conditions that means you have to provide the zero derivative of the function at the initial time step at t equals that means the position of the atoms are i and the initial wave function that is the ground state then you have to give also the first derivative of these two quantities so the velocities and the derivative of the wave function now the problem is what is that ok we could in principle set this to zero don't worry about the parameters i will explain them later so this how to guess this one can think ok i put a to zero it is ok also to put it to zero at the beginning but it is not really physical so what we want to do is to find a way to estimate the wave function velocities for the initial condition of the second order differential equation so as i said the nuclear for the nuclei this is simple but for the electrons what to do it is not possible to do this that means this is not possible because this can have different gauges because the minimization of the ground state is independent nothing nobody told the conjugate gradient algorithm to choose the same gauge for the two minimization they are independent completely so if you do this difference let's say for example that this here there is a different phase and you are subtracting something that is not correlated and you get on the random numbers so what we do is to do a trick to use the parallel to fix the gauges practically and we use the parallel and transport gauges that means we want to compute this ok let's write something trivial the projector over the occupied main fault and we say ok the wave function is equal to the projector to the occupied main fault time to wave function i've written nothing here now let's do the derivative i do the derivative of the projector then over the wave function and the parallel transport gauge that will not i will not explain here in the details but it is equivalent to say put this term i define this term as equal to zero in this way we are able to compute this because this object is well defined also with if we use the finite difference method to compute the derivative that means it does not depend on the gauge that the minimizer take quantity gradient choose to compute the ground state exactly we compute this and we differentiate it numerically and we use it to compute to compute a new state at time t using the previous step at times t minus dp and this object ok in this way you see that if a particular gauge is selected for this wave function for this is the same one and this will have the same gauge same gauge as as t minus dp so we are ok we are ready to start the ballet algorithm a little note that right now this is implemented only for non conserving pseudo potentials so if you don't use non conserving we will have to start with the velocity of the wave function equal to 0 ok now since we are moving to carparinello I will start to explain you the parameters of the Hamiltonian and the Lagrangian so you see there is this mu that is a mass for the electrons you can imagine the electrons I like to imagine the electrons are some sort of classical fluid like you have the the potential you have above let's say you have some water inside it and this water has a mass ok then there are other stuff that maybe it's not very important but you have to know that it exists that are auto normality constraints that are needed because you want the wave function to be orthogonal to each other so there is a standard method in Lagrangian mechanics that distributes some Lagrangian multipliers that are here only to ensure that the wave functions are auto normal these are quantities that are computed by an iterative algorithm and produces some sort of additional forces on the wave function that make sure that the wave function are also auto normal each other ok then ok spin at this this is the electronic fictitious kinetic energy that is called in C in the output files that is the equivalent of the kinetic energy of the nucleus but for the electronic fictitious freedom but here is by far more complicated because you see you don't have here you have for example we have in our case we have 24 atoms but here we have a lot more degrees of freedom that is the number of plane waves that is called NPW for example times the number of of bands so you have a lot of degrees of freedom inside here under the kinetic energy of all these degrees of freedom is called the Tron fictitious kinetic energy in the code is written like in C in the output files also so that's it then the mass parameter how to choose the mass parameter ok as professor scandal told you before what we want to do is to have this potential we have that's for simplicity that the electrons are like like water inside the bowl the potential moves and we want that the water moves with the bowl and that's not so suppose for example that you use a very large mass for this electron for this water inside the bowl it will be very difficult to move the bowls and if this is an ion you will have that the mass of the ions is like the mass of the ion plus an additional mass that depends on the mass of the electrons that are attached to the ions I repeat again electrons are like some classical fluid that is attached to the ions by the potential so if you choose an electron the mass that is too large you have the issues but this is not all you will have another issue that is the issue of of the coupling between the electronic degrees of freedom and the ionic degrees of freedom maybe I don't know so what is the issue it is that we have two system, different system we have the ionic system and the electronic system we want that the electronic system stays always cold that means the electrons are near the ground state and keeps oscillating near the ground state but what happens if the two system exchange energy since the electrons are very cold and the ions are hot because they have a temperature like a few hundred degrees or Kelvin's or more if there is an overlap between the frequency of the ionic system and the frequencies of the electronic system that are very fast because the electrons are light if you choose a low enough fictitious mass if there is an overlap let's suppose that those are the frequencies this is the ionic frequency and this is the electronic system frequencies if there is an overlap the energy we start to go from the ionic system to the electronic system so what happens is something like that the electronic system goes up and eventually the electrons will go out of the ground state and you will end up with very expensive random numbers so what you want is that the highest frequency of the ionic system is smaller by big enough amount of the smaller frequency of the electronic system and it turns out that that the smaller system of the electronic system is proportional like I call it the smaller frequency of the electronic system is something that goes like the gap the electronic gap of the system over the fictitious parameter over the square root so you see if you put here a too large number eventually this small frequency will overlap with the frequencies of the ionic system energy will exchange and the electrons will go away from the ground state so as a rule of thumb I start with an emas that is the name for the moo parameters in the input file of about 15 units of atomic units that is a reasonable enough value that usually allows to to have a good allows the electrons to follow the minimum very well then there is another thing on the other side that since you are integrating the electronic degrees of freedom with the verlet algorithm you want a long enough simulation because you want to simulate for example one picosecond let's say if you choose the lower the mass that you choose the faster will be the oscillation of the electronic degrees of freedom because they are lighter if you have a lighter object it will run faster so with the verlet algorithm there is a part as you saw before that is something like that you have a mass then you have delta t squared so if you divide this by two let's say to have the same magnitude of the displacement we have to divide also this by the square root of this so the lower the electronic mass the lower the time step and since the electronic mass is the lighter element in the system you have to choose a smaller time for the wall simulation so you have to balance between a small time step and an accurate enough simulation this usually is a good guess to start with so let's go through the input file ok calculation restart let's see sorry for being long but there are a lot of concepts here but I find it very interesting by the way so CP molecular dynamics so you see we have to change a few few things we want to restart that means we want to read the ground state and the electronic velocities from the restart files we had we calculated before with the previous run and with the velocities calculate with the projector treat let's say we have to tell the code number directory read that is the number that identifies the folder where the code is going to read the restart files then you have to tell number directory write and I increased this by one I do not overwrite the previous directory so if I want I can restart again from the previous calculation then I want more steps and since Carparinello is usually cheaper for a single step I can put here a much bigger number then this is very important this switch on the Carparinello equation of motions so you should tell the code electro dynamics equal to Verlet that is the electrons name list we are using the equation of motions of Carparinello then the last very important thing is to set an emas to a reasonable value you can start from this and then see if it is possible to increase it then we will see how to check that the emas is good enough or decrease it if the simulation does not go well as I said before if you choose a too large emas it can happens that this becomes too low and the electronic system and the ionic system exchange energy electrons heats up and everything goes bananas as the code will tell you ion velocity is equal to default test the code to read the ion velocities from the restart file remember to not forget here random otherwise velocities will be overwritten with new one velocities and also the autopilot card is not necessary anymore ok everything is the same so now we can submit also this I used the reverse search trick that I like pretty much so remote and here I am close under if you press control R more time you see the command so this is true that true Carparinello multiple bananas so I run it I hope that I'm not saying too much okay some somebody asked what is the K nuclear maybe I did not explain it better so this is okay the DFT is this one the DFT of the previous slide okay this is K nuclei I called this part this part is the kinetic energy of the nuclei and this is the kinetic energy of the electrons okay now we are going to see in more details the output because we will have more so let's see if something happened all the output and why some parts some energy change some other energy does not change so okay it worked sync from APC you see there is a new folder of the new restart directory and the trajectory is updated an important thing if you keep the same prefix for the trajectory every time you run the code the code is going to append at the end of the file the new quantities so it is not going to overwrite the file but it is going to append new atomic position new atomic velocities, new cell vectors, new thermodynamic quantities at the end of the file so don't worry the time step the number of times is always increased between each restart file so when you run the new the new input the first step will be the number with the same number of the last step of the previous run so I already showed you how the output appears in the terminal but now let's use a new plot I prepared a small script that it's nothing special, you can see it it is it plus some columns of of the output files so I will run it without the plot termo I don't know if you know but if you press tab in the terminal the terminal will try to complete what you are writing so you can type faster ok I hope that this is clear enough these plots so what we are seeing here we are seeing the output of the grid the thermodynamic output of the code so you see here there are the few steps that we ran before with the conjugate gradient algorithm then here at this point we switched off conjugated we did the restart the restart file the legal instruction we did the restart file and from this point on we are now running the Karparino algorithm so here I am plotting this is the constant of motion this line between maybe I can zoom you see this line here this line here is the constant of motion of the Lagrangian each Lagrangian has many constant of motion and the integration time step is small enough this constant of motion will be constant within numerical accuracy that means that the parameters of your simulation are good with respect to the verlet integration so it does not mean that your simulation is good it means that the time step the delta t is small enough then I printed the electronic kinetic energy that is this quantity here and I shift it from this point so the electronic kinetic energy is only from here to here this one and you see that is approximately a constant so if the simulation is good this electronic kinetic energy will be more or less almost the same during the simulation if the temperature does not change that means that there is no exchange of energy between the ionic system and the electronic system that in turn means that probably your the forces that you are calculating are good enough so you have to check around the Carparinello equation of motion that this kinetic energy of the electrons is always the same if it is always the same means that there is no energy transfer between the nuclei and the electronic use of freedom so you can trust the forces that the code is computing for you then ok the constant motion is nothing physical it is not like the energy of a classical simulation because you have this term inside it so the constant of motion maybe I wrote it before yes a counter this is called a constant in the output file it is the potential energy it is the potential energy plus the kinetic energy of the ions is something physical plus this part this part is not there is no physics inside this number because it is only an energy that arises from the way that you are using to minimize say minimize but it's not a minimization you saw before it is a kinetic energy that arises from the way that you are using to compute the forces at the end of the day so it is not related to a physical quantities also because if you change a move if you change the fictitious mass of the electron this number here the kinetic energy changes so it is not something physical so this constant of motion is something useful to know but it is not something physical the physical energy is the sum of the nuclei kinetic energy plus the dft potential energy the dft energy that is this line here under the constant of motion by the way you see that it is evident from this plot that the constant of motion is this plus the electron you can see that the plot are mirrored around the constant of motion so if your simulation is good or if your simulation is good the physical energy the physical total energy of the system will be more or less a constant ok then let's see if I have something more ok the potential energy also in this case when the system is equilibrated the potential energy will oscillate around a constant value that is ok in this case the system is the simulation is very short so it is difficult to say if later this energy will stay more or less constant or it will go up or down or whatever but since now everything seems ok so then here I wrote something about the forces you can ok the forces as I told you before and also as Sandro Scandro told you before forces are not exactly the ones that are given by the PW code for example because the ground state keeps oscillating in the minimal of the potential and and it is not always exactly the same but on average the forces are about the same of the one that the code pw.px calculates but there is another fact that maybe you can imagine that if let me put another if let's say that is the potential here you have the electronics that has a mass so you have a mass that is added to the ionic mass how is this reflected in the forces is reflected in the fact that the forces calculated by the cp are a fraction of the forces calculated with the true Born-Oppenheimer method and this is something that is usually like 99 something so it is nothing that you have to worry about but be aware that if you put this too large these numbers can go to something like 0.8 so you have if you are computing the statistical properties of the system those does not depend on the mass of the ions plus the electrons let's say but if you are computing dynamic properties like diffusion coefficient you have to be careful because if you choose a move such that the forces calculated by a factor of this you are computing wrong diffusion coefficient so be aware of this fact so here I was suggesting to compare the forces a part of the factor 2 with an electronic mass of 50 atomic units you should not see nothing strange so you should see that the forces are more or less the TW1 ok so this is done now let's move to the next step how we can choose the mass of the electrons I thought before you can start with the mass of 50 or something like that as an educated guess then you can check first of all if there is no transfer between the ionic degrees of freedom that means that in the plot of the kinetic energy of the electrons the kinetic energy of the electrons remains more or less constant then as a last step you can check the ratio between the forces calculated with the Carpanello method and the forces calculated with the Born-Oppenheimer method if the mass is low enough the forces will be similar within oscillation typical of the Carpanello method but that are not harmful to the simulation ok then let's start with this thermostat so as Sandro thought you before if you start with a niche condition choose and randomly say you will never get the temperature that you want so let's switch on the thermostat ok thermostat let's go what did I do switch text beam water 8 3, no say over no say over so you will have a direct rewrite read and number as always we don't want to have a write restart because we may want later to restart from a previous step of this so we increase this this has the point of course to a valid restart the directory then I want more steps to me enough here there is no a precise recipe it depends really on the system the time needed to equilibrate at a given temperature you will have always to look at the output files and see if the system is equilibrated or not so to switch on the no say over thermostat we have to put in the name list ions and the temperature in Kelvin very nice so what is this ok you see I don't know that means you copy the velocity from the previous restarts the temperature of the thermostat then there are two additional parameters that are very technical I will not explain them if the restarts you have to study the way the no say over thermostat works but just you have to know that there are additional parameters that are explained here for example you see ions you have no zip that is explained here it is just a parameter that you can tune to have a better equilibration at your temperature and also this is the number of chains in the no say over thermostat it is a very technical thing just you know that there is so if it is difficult to equilibrate your system at a given temperature you can try to change this parameter and see what happens this usually you can choose so it overlaps more or less with some peaks in the vibrational spectrum of your ionic system ok so let's submit this to the cluster so remote mpi run and it is step number 3 no say over.in it will take a little bit thank you sorry you are waiting there are questions on youtube if I particularly want to ask if I can use any kind of self-dopotential for the capparinello as far as I know yes but be aware for example there with function velocities feature that we used at the end of the Bernofenheimer step that we did first before it is implemented only for non-conservative cells so there are some features of the code usually the more advanced they are the less the number of pseudo-potential types that are implemented so for example everything implemented and but yes there is no you can use whatever pseudo-potential you think is good enough for your system for ultra-soft pseudo-potential there are additional parameters that maybe I explained here in the input file that is this number here it uses a different algorithm with respect to the PW code and you have to specify this number for the ultra-soft that are some sort of grid that is around each atom to contain the additional charges if I'm not wrong so be aware if you use ultra-soft there is this additional input that you have to provide and then yes you can use whatever you want there is another question if you can use the optimized gamma-trick everything is gamma-trick because we have only gamma-point we don't have k-points so you have to check for convergence with respect to the system size in place of the convergers with k-points and the codes use only gamma-tricks ok then you can leave the default basically the recommendation is if your system does not equilibrate fast enough you can try to plot the vibrational density of your ionic system and try to choose a value that overlaps with it for example in my case I put 10 wrong yes I put 10 ok this is a frequency in terabytes that has to overlap in a good way with your system to enable a faster exchange of energy between the thermostat degrees of freedom and the ionic degrees of freedom but it's really technical this one it's better to do some trial and error and see what happens after reading a little bit of theory of course so you know more or less what you are doing ok I think that maybe the code finished they are all very fast run the mode way way yes yes sync from apc the new trajectory no say over ok now we do the same step as before to see what happened let's plot the same interesting so you see the data is always attended at the end of the file here at the beginning you have the plot that we had before this is exactly the plot that we had before but then we switched on the no say over thermostat from this plot what you can see this is the physical total energy that is the kinetic energy of the ions plus the potential energy for the ions that is in that in this case is the density functional theory total energy you see that it is no more an approximate constant of motion because I remember you that the exact constant of motion is the one that includes also the electronics degrees of freedom plus of course the thermostat degrees of freedom that are the additional variables that you introduce in the equation of motion to fix the temperature you see that it changes because we switched on the thermostat and more or less after a certain time you see that it starts to oscillate around some values so it looks like it looks like the system wants to equilibrate around this physical total energy then you can see that the electronic kinetic energy went a little bit down why did it went down because the temperature of the system is lower by the way later I'll show to you also the temperature also and you see also that the constant of motion of the Lagrangian is a constant of motion so we are happy that our time step is small enough it looks like everything is going well this is the potential energy it is equilibrating around another different values of course because we switched on the thermostat and ok let's see the temperature the temperature is both column 5 5 yes you remember that at the beginning we chose a temperature of 600 Kelvin for the initial initialization of the temperature with same thing from the Maxwell-Boltz distribution I wrote it in the first into 5 yes it was here 600 Kelvin but looks the system choose a temperature that looks like the double of twice the 600 Kelvin why? because as I said before the initial state was guessed by randomly so apparently there was it was a configuration where there was too much potential energy so when I started the verlet algorithm the system took all this potential energy and this went into kinetic energy so temperature then I switched on the thermostat at this point you remember we did we did 100 steps of conjugate gradients of this the very very beginning 100 is here you see on the bottom left corner there is the x and y values for the plot I switched on the verlet or then at 1100 that is here I switched on the nozzle over the thermostat to see the system goes down with the temperature very fast and then it equilibrates around the value that I selected 600 Kelvin it is the right value now I'm going to switch off the thermostat and let's see what happens so why wasn't it seen before so let's go ahead so now our simulation is working very well but now I'm going to explain you a little trick that can be useful if the simulation is more problematic what do I mean with more problematic I mean that maybe there is some energy exchange between electronics degrees of freedom and the ionic degrees of freedom so the kinetic energy of the electrons start to go up and maybe I don't know maybe it is too high or the forces are starting to become random numbers so it may happen that you need to minimize again the ground step so what we do we do the same thing that we did at the beginning but in this case we need the only we need only one step because one step is enough for minimizing not that we stole the one step in the input file but in reality the code does two gg minimizations if the norm conserving pseudo potential is used it does one minimization for t-minus dT and one for t then it computes the projector and then it computes again the c-minus dT with functions so you select you say one step inside the the code does two steps of conjugated gradient to compute the projector over the occupied manifold so this is easy to do and also very fast because we have only one step so let's see the input file waterA2.4 conjugated gradient 15 as always I increase by one number directory read and number directory write I do again conjugated gradient for the electrons so no Verlet algorithm no Carparinello we are doing Born-Oppenheimer for only one step and ion temperature not controlled so I switched off also the so we do another minimization the time step I left the same as the Verlet algorithm with Carparinello so I don't have to use the ugly change time step input of the code everything is the same the electronic mass is ignored when you conjugated but I left here nothing happens ok here everything is the same ah I forgot to tell you this part of the input file the atomic position file is the same that I used at the beginning but if you use the restart feature of the code this part will be simply be ignored if you have to put it it is sufficient to copy and paste the atomic position that you have from the first input file but the code will not use this position to do the computation it is going to use the position from the folder number with the restart file ok so this is a very fast calculation ok there are some methods that run are slow ok I will go slower but now it is almost finished so don't worry what I was about to ok submit the minimization this is a very fast step in now we are using I am doing this to use this trick we compute a new gram state and we compute a new gram state velocity let's say so we are happy and why this is better than a velocity of 0 because this is somewhat consistent with the dynamics of the ions it is like with the water and the bowl so suppose that everything is moving or accelerating it is accelerating in this direction the water is like this and here I choose a velocity that is consistent with the acceleration if I don't choose this maybe I can go the other side like this and then like this and it is not ok you will have at the end of the day more kinetic energy for the electrons also you can start with 0 a velocity of 0 imagine that you are on a running car with a bowl of water on the top of the car the water is consistent with the velocity of the car if you put some water at the rest on a car that is running the water will start going out of the bowl ok it will gain kinetic energy ok s