 Well, if we were in Trieste, we could go to a nicer restaurant, have a fish, but we can only imagine. Well, we'll be compensated by a nice talk, okay? Well, that's the thing, you're in bad shape. Okay. Yeah, so maybe it's a good time to start, you agree? Okay, so good morning, good afternoon, and good night to you all over the world. And my name is Lucia Reining, and I'm a CNRS researcher at the Ecole Polytechnique close to Paris, and I'm chairing the session. And we have the great pleasure to have as a first talk, Michele Parinello, who will speak about metadynamics of paths. Okay, so good morning. Thank you for the invitation. I was telling the organizer, I really hate giving Zoom talks, but total energy is special for me. I was there since the very beginning, some 40 years ago, and so reminds me of my youth when I was in Trieste, and also the Trieste kind of philosophy and science style, I would say. Okay. So let's now start, and I will start with one of the main preoccupations that I have had in the last decade is that we're trying to address the problem in molecular dynamics, as you know, molecular dynamics is a very powerful tool, but it has like all things, all human things, the limitation. And one important limitation which always we are trying to overcome is that on the time scale, that simulation, molecular dynamics simulation can explore. And in fact, there are many physical phenomena here, I illustrate two such phenomena, one is protein folding, the other crystallization, they take place on a time scale which cannot be reached by brute force, molecular dynamics simulation. So that's what we are going to talk about today. And before going into the details, I just want to establish some nomenclature, so one concept that we're going to use over and over during this talk is that collective variable, which for physicists, that's a common, it's not such an alien concept, if you can see the Landau's theory, what they call about all the parameters, or even chemistry, similarly reaction coordinate or whatever. I mean, this is not a new concept. But so imagine that we have by hook or by crook obtained some, what we call a collective variable, so a variable which encodes the flow degrees of field of the system. Then we do the marginal of the probability distribution, probability distribution is the Boltzmann distribution of S and we get possibly a distribution like the one that I tell you, this would be, would correspond to a two-state system, so the system can be with probability here, with some probability there, but with very low probability in the middle. And we can express the probability in terms of free energy, and in the space of free energy, which is the log, you find that there is a two-state system. This is just two ways of representing the same phenomenon. And here you talk about transition state, which is a low probability situation. So if you, so what happens if you start from here, and since the probability of going here is so low, you remain stuck. And similarly, if you start safe from there. So what we want to do, our strategy is to make some transformation of the distribution function such that the probability in the middle, which allows you to go from here to there, is increased. And we can do that by designing an algorithm that takes, in a controlled way, the original probability distribution into a target probability. For instance, we start with these two systems, so we can make the final probability look flat. And this would be the case for those who know these things, what happens in umbrella sampling. Or in metadynamics. We do something a bit more sophisticated. We don't want totally flat, because then we sample with the equal probability states, which are really not relevant. So these states are not relevant, and these states are not relevant. And so instead, we do what we call well-tempered. So we start, this is the red, this is the original probability. We want to make a transformation such that the final probability is distorted, and so we have a higher probability of passing this transition state region. And I mean, the mechanism evolved a lot over the years. The first version was done with the, okay, sorry, the first version was the simplest, and then we evolved in well-tempered metadynamics, and then we go from, we lie, I'm sorry, the name was done, we lie, and then the most recent version is OBS, which has the present evolution and what we use nowadays. And this is obtained going from one initial distribution to a final distribution is to add or construct a function of the collective variable, which is added to the Hamiltonian that changes the probability distribution. And of course, you have to be able to go from the original probability to the bias probability and from the bias probability back to the original. And this operation is called re-weighing. The OBS, which is the latest, would require a lecture to be explained because there are many subletters and many tricks that make it very powerful. But the gist of it is we write the bias in terms of the target and the probability distribution. And on the fly, like a metadynamics, rather than changing V, we change POS by adding Gaussians wherever we go. In this way, we can do this nice transformation and observe all sorts of phenomena. And I thought that since many people here have electronic structure minded and application minded, I would like before moving on to the next chapter, some example of applications to real problem of the methods that we have developed. The first one is the study of the phase diagram for instance of silicon. Silicon is one of the most important elements. Solid to liquid transition is accompanied by a non-metal to metal transition. So liquid to silicon is a metal, crystalline silicon is a semiconductor. And there are two problems. One is the time scale because crystallization takes place on microscopic time scale. You need to have a good potential that describes in the realistic manner the interaction. In fact, this was one of the motivation why a long time ago with Roberto Carr, we developed ab initio molecule analysis. Of course, ab initio is of course very powerful, but it's also very expensive. So the issue is, can we do an ab initio quality calculation at much more cheaply? And that's what we started doing some 10 years ago, we had better and nearly to collect a large number of configurations. On this configuration, we train a neural network, which has a special structure. So you start the input at the positions of the ions. They are symmetrized manipulator, so you can respect the symmetry of the problem. Then you feed it into a neural network. This neural network is equal to this, is equal to this. And each one represents the contribution to the total energy of one single atom. And I believe you have talks about this in this work. Okay, so when you train a new work, you need data and you have to choose the data on which to train the network. And the issue is with neural network, that neural network are very good at interpolating. So if I have data, I have a good representation of this function, but very poor at extrapolating. And so if the yellow is your real curve, as soon as you move from the region where you have data, you're in trouble. In the context of what we want to study, which is a phase transition between one phase and the other, which is where the system goes through a transition state, the problem is that even if I do a standard MD, I only find a configuration in A and B. So, I mean, to study the transition, it's very important that I have data on the transition state, otherwise the barrier will be too high or too low. I mean, it will be totally wrong. So the way we do it is we use meta dynamics. So we start, we use a potential, in this case, a cylinder variable potential, some simplified potential. We generate configuration here and there, but also we do meta dynamics. That means that we are not stuck here or there, but we can move from here to there and pass through the transition state. You see here, there's a nucleus of solid phase in a sea of liquid phase. So meta dynamics enters directly in the construction of the potential, not only as a sampling tool. And we have studied the silicon. That's a given sampler of things that you can do in this way. So these are all ab initial quality things. But we studied the phase diagram of gallium. Gallium is a particularly complicated system with a special bond. We don't have time to go into the details, but we can really understand all these four phase transitions. We can use, we can, we can study chemical reaction in solution, in the full complexity. This is the case of you read the compositions in water, which we, and it's a very recent application is to this nice system. Pt, I don't, maybe I should have talked about this, but I mean, now it's too late. And the phosphorous as many, many solid phase, as you know, but also in delicate state is the odd behavior. It has two phases. So it has a molecular phase here and the polymeric, which is insulating and the polymeric metallic phase. And through these tools that we have developed, we can study these transitions naturally. We find that there is a critical point here, a secondary phase transition, which is associated to a metal, metal transition in the liquid state. Okay. So after this introduction, which is also in the kind of establish a bit language and motivation, I move on to the main subject of my talk, which is this meta dynamics of past. Okay. For doing that, I will recall how you do, you sample Boltzmann distribution, which is written here, using molecular dynamics. Okay. So you solve Hamilton's equation and you add the thermostat to keep the temperature constant and you express your trajectory as a sequence of positions. So position from position N, you move to position N1 and so on and so forth and similarly the velocities. From this trajectory, you can calculate average properties, simply summing or also correlation function, which in the practice are written because you have an observable, a time i, then you want to see how it correlates with the same operator, different operators at time tau, and that's how you write things. Okay. So that's all standard. Now, what we want to do is rather than sampling configurations as you do, we want to sample entire segment of trajectories and in a certain limit, which is called more small cosky dynamics, we can express in finite form the probability of finding a certain trajectory from zero to tau. And this action that appears here is written, it's an integral over the, or here's a bracket missing, over the trajectory r of t from zero to tau. Okay. So from here you can sample trajectories. Okay. So computer doesn't like continuous object and so we have to discretize as just as you do in standard the molecular dynamics. And so the trajectory is nothing else, but a succession of points where the index here is really a label for the time step where you are. So when you discretize your trajectory, the velocity becomes a finite difference. Estimator of the velocity, this is the force and so on. And these are the ionic position. Okay. So to, if you want to visually see what happens, we have the system repeated in times. As you move along here, you move along the time. Each snapshot, each instantaneous configuration is linked to the other by this expression and by the forces. Okay. So that's, that's what. Okay. So we want to look for probability to start a certain point, say A in the, in the basin away. So I put this here and the probability now becomes like this, where this is something complicated basis, nothing but the classical Boltzmann. It only the system is repeated at times as you do, for instance, in party integrals. But you, you can treat it as a standard classical problem. So by doing this, you turn what is a dynamical problem into a static one. And so you can apply all the techniques, including, including meta dynamics. Okay. You might say, okay. So that's you, you pay a big price because at the origin, you had one system. Now you have to multiply it by n if you do n steps, your system, but it is trivially parallel. And so you have two advantages by doing this approach. One, you turn something dynamical and something static or look that looks static. And then you turn a serial problem up because when you do molecular dynamic, you start with the configuration, you do all your calculations, then you update, then you update, then you update. And this is naturally serial. And so it's not cannot profit that much from parallel machine. In, in, in this path approach, instead, you, it's a really, it's a trivially, it's a trivially parallel formulation. Okay. Let's see how, how it goes. So I have this model system. It has one minimum here, one minimum there. There are two barriers. So in this case, we took symmetric. Sorry. Okay. Sorry. So this is one of these paths, which we represent as a, as a, it really is a bit like a polymer. So the initial position say it's red and the final position is black. If I don't do anything as I sample parts, all the parts will be confined in this minimum, just like the configuration, because these are the most likely parts. If I start with a, the system remains in place for a long time until the rare fluctuations takes you over the barrier, but nothing happens. So now we know we can treat this dynamical problem as a classical thing. And so we can apply the non-sampling method I briefly discussed earlier on. And so what we use as a CV in this case, what do we want to see happening? We want to see the initial and final position to be very distant. And so I use a CV, the distance between the first bit to the last bit. And if I do that, then you see what happens. The system can go from one side, makes a lot of fluctuations, then you make a successful transition. Then it goes also from the other side. So there are methods that call part sampling where you start with a guess for the part and you remain trapped here. But this method allows you to go forth and back as you please. Okay, so that's the basic thing. Okay, now let's do something a bit more physical. The case of ammonia, as you well know, ammonia has this umbrella motion where you go from one side to the other. And this is basically, it's a minimalistic part in which I put the initial configuration, the transition state, and the final configurations. In reality, if I want to be realistic, I have to divide the part in more steps. So we divide it into 200 replicas. In this case, it was sufficient. We use this new method, opens. This is a sequence of configurations that belong to the known, when the trajectory where no transition is observed. This is a trajectory in which a transition has been observed. Okay, I can look at the distribution of trajectories. And there is an interesting aspect that happens. If you look at the action here, when I'm close to the transition state, the force is close to zero. So this spring becomes operative. It's not shifted to a new position. And so here I have an accumulation of data, of transition state data, which can be very important in many cases. Okay, now we've done the calculation. How about calculating dynamical properties? We know how to do this. As I move along here, I move along the dynamics. In this case, discretized time axis, and I can calculate the correlation function, just as we did for molecular dynamics. The same thing is like here. Since I do metadynamics, I have to do reweighting, but this is a technical detail. Basically, once I have a sequence, I can calculate the other correlation function. In particular, I can calculate the probabilities. If at time zero, I'm in A, at time T, I'm in B. And the delivery of this formula, so we shall see, gives the rate with which the system goes from A to B. Let's look at this thing. For ammonia, we can calculate this correlation function. A different temperature from low to high temperature. This correlation function has an asymptotic linear behavior. The slope of this curve, this linear part of the curve, gives you the rate. Then you can do an arranged plot, a different temperature. That below is the arranged plot in linear. I can extract the transition state energy, and we find this estimate from the feet, and then if you do the exact calculation, you find a very, very close estimate. Then let's make our life a bit more difficult, and we look at the case of ammonia in water, the same mode, the same umbrella, but now this time in water. We prepare everything. We do open, the CV is the distance between the initial and the final configurations. Then we do the simulation. We can do the vacuum, and the water, of course, the drag of the water brings down the rate. As you would expect, we get two different rates, one in the vacuum and one in the in solution, and the ratio of these two rates is similar to what you would have predicted by transition state theory knowing the difference in free energy in the transition state. So the things seem to be fitting together very well. Okay, so let's look a little deeper into what happens, what is the role of water. Okay, if we look at put together the reactive, the non-reactive part, we see that there is a dominant structure here, and it is the hydrogen bond that this water molecule makes with the hydrogen. And if I look at the precondition function, you see that there is a strong peak around the sodium which represented this distance. Now let's look at the reactive part. Now in the reactive part, well that's very difficult to understand from this picture what happens, but for sure this hydrogen bond is almost disappeared, plays no role. So breaking this bond is crucial for the event to take place. A more detailed version of this effect can be seen if we can look at the water, oxygen and hydrogen distribution around the in this plane. That's what happens, that's what you expect here. You have nitrogen, your hydrogen and the oxygen part of the part. Okay, non-reactive trajectory represents it's just a check and then the standard equilibrium md gives the same thing. So you can sample A both in the ordinary way or through the part dynamically get the same result. Now what is more interesting is the case of the reactive trajectory. So for the reaction to take place, the solvation shell around the nitrogen has to change nature, has to become almost symmetric. If you look at this, if you make your nitrogen flat and you force it to be flat and you look at the solvation around the flat and the ammonia, you see that this is the same symmetric configuration that allows. So you have to have a fluctuations of the water solvent such that the energy of the transition state has the optimal solvation energy. So that's how much time do I have? Actually, the speaking time would just be over and then we have another... Okay, so that's fantastic. So I can skip this thing, which is the marginal relative to the thing. Yes, I'll show you this little worm. That's a trajectory and I have a system with two states. That's a mini protein that goes for the back. But I think I got through you. I want to thank the group that I had over the last few years. I let you guess which picture was taken in 2020. And I thank you for your attention. I'm ready to take questions. Thank you very much. So we actually do have time for questions. And I can see two questions from participants and two raised hands from an analyst. So let me start with a participant question. The first question that came in was a quite general question from a participant called Juhid Gupta. Juhid Gupta, do you want to ask your question yourself? I don't know if the person is online. Yes, ma'am. Yes, ma'am. I want to ask how to do the phase transition in the density function theory. Katia. The question was how to calculate a phase transition using density functional theory. So a very general question. Okay, so the way we do is to first determine the potential. As we said, training a neural network to have a potential which has a DSP like accuracy. So that's the potential part. Then we have to define an order parameter that can distinguish between density and solid. We have, there are many possibilities, but one of the most physically transparent CD that we can think of is the blood peak. It's a lack of experimentalist look to see whether there is a phase transition. So we do a metadynamics of a simulation using the CD so the system can go from liquid to liquid, from liquid to solid. And then we see for the focus of the point in which the difference between solid is. Okay, thank you. Okay, so I suggest to move to a question of a panelist. So the next one would be Francesco Mauri. Hi, Miguel. A very nice talk. A question. So the whole example you show, there is a correspondence between the pati integral MD and the constrained classical MD. Probably because of a renewed type of reaction. You have example in which you get something more than from standard DFT. Well, I mean, I would say the role of the salvation was not clear before our calculation. And I can imagine a more complicated system where it is difficult to know exactly the real condition of state. So this is the very beginning. So we do it, we know that there is the traditional state. In many cases, this is not the case. In the application that we have made so far, we have taken segment trajectories which are not very long because that would require massive parallelism. So in the beginning, we have focused our attention in situations where the transition between A and B is very fast. And these are the situations in which the method has been tested. In a situation more complicated, there were, you say here, a mesa and then the system goes wandering in the mesa before falling into the other. That's more complicated. But I think also this could be tackled. But if you ask me if I have example now, no, I don't because we just published the paper. Okay, so now we start to have more questions from my participants. So next one would be from Seyyed Vaid Hosseini. If you want to ask your question yourself. Hello, can you hear me? Yes. Yeah, this is Bahid from University of Tehran. I'm wondering to know how to include time reversal, rotational symmetry and translational invariance of the solids, especially for molecules in the correlation functions for calculating the properties by molecule dynamics. There are two issues here. Okay. One is time reversibility. Translational invariance, yes. Translational invariance, so that's easy. In what case? Because there are two parts to my talk. We're talking about standard. We are talking about the last path. A standard molecule dynamic calculation, a standard. Okay, the linear momentum, it's easy to enforce because you calculate the force, the way it is done, you calculate the force on particle i due to particle j. And then you put the minus sign and then you get the force on day from i. And this automatically concepts the linear momentum. If you start with the zero linear momentum, you continue to have linear momentum. The question of the linear momentum, we have to distinguish two cases. Whether you have periodic boundary conditions or non-periodic boundary conditions. So if you have periodic boundary conditions, you break momentum invariance. And so that's not the concept quantity. The other, I think, I never really worried about it because I don't do planetary motion and things like that, where conservation of angular momentum is important. So it's something folks, I don't think there is a major issue there unless you want because the question that you use, if you use Newton's question, that's angular momentum is going to be, that's the accuracy of your... It is conserved. Yes. Thank you. Thank you very much. Thank you very much. Okay, I think we can take one last short question maybe from a panelist, Guillermo Mazzola. Thank you very much. So a short technical question. So if you study reaction involving tunneling, do you observe the tendency of these parts to perform some kind of instantonic trajectory, like with all the replica on the left or the other on the right, because no one wants to sit on the top? And if you do, how do you avoid this? Well, instanton they come in a quantum context, I believe, and we are not treating the quantum part. We could do it in a bit. So suppose you can use path integral molecule, path integral description of the quantum system and look for the path that go from A to B. And I believe, and I get wrong here, but I believe that one approximation that is used to discuss this thing is about instanton, as you say, in this context. But that's nothing to do with what the classical dynamics, because there ain't no tunneling. Okay, no, because it's an example of ammonia where I saw that. Let me see. No, no, actually it's something different. Let's go back, if I may. So I think if the answer will be long, we should maybe... And here, yeah. Here, I see the chairman. So this is the order parameter, which would be the height of the the tetrahedron there. So you see, you find the lots, lots of data in the transition region, which is also false, so there are zero. Suppose you were at just enough the time to go on top of the barrier. Ideally, we would say they are forever, but then your noise is in there. Okay. Okay, so I think that unfortunately, in order to save at least five minutes, maybe, of coffee break, we have to stop here, there are more questions in the chat, maybe Michele or someone else can answer that. And so thank you for the good discussion, and thank you very much again, Michele, for having presented these exciting ideas. And I would say we recombine in five minutes for the next talk. Okay. It's closed now. Lucia, what do I have to do now? You simply have to pull away to a slide. It means that the next speaker will technically have the possibility to put his slides. Okay, okay. I think it's enough that there is an unshare button or stop share button. It's already done, it's already done. It's already done, excuse me. Okay. So we are in place. Five minutes is okay, right? I mean, it's all right. What are the things, what are the things? Since we are not in a huge heart for coffee break. Maybe it's because to have room for eventual delays. We put all the 10 minutes pause to accommodate for eventual delays. No, I mean, the coffee break has to be shortened. Yeah, yeah, okay. It's like a dictatorial way to five minutes, and we will do what we can, right? We could try the screen too. If you have now, since there's no coffee break, if you have to bake your coffee, then it's a problem. I mean, there's no queue, so there's no problem. Also to test the screen-sharing of the new speakers. So Zvastika, David, Lorenzo, Massimiliano. Okay, I will be waiting for one minute. You can do this, right? You do this. Okay. Okay, so who shot? Sorry, Cola. Should I try to share my screen? Yes, please. Okay, so I would share this here. Okay, can you see it? Yes, fine. Okay, perfect. So if I do presentation. Okay. I can. Okay, so you can see my cursor, my mouse pointer. Yes, great. Okay, thank you. Then I will answer it. Okay, okay. Could I please share and check? Yeah, Zvastika. Yeah, maybe you can try last because you will be... So you leave it on directly, okay? Okay, okay. Lorenzo, are you there? Yes, so I try sharing the screen. Can you hear it? Perfect. Can you see? Can you see also the video now that I put the screen? Yes, yes, yes. Okay, okay. So it works in this very technological way of including video in this lens. Like running separately, right? Yeah. Okay. Okay, so Massimiliano. Yeah, okay. You see it? Yes, can you put it in? Yeah, sure. Okay, perfect. Does it work? Yes, perfect. Yes. Okay, stop share. Yeah. Okay, so Zvastika, it's your turn now. And so then you can leave it on directly. Okay, perfect. Yeah, perfect. So your name is Zvastika Chatterjee. Yeah, so it's like my passport has been here. So when I was registering for this, it was the passport name that was asked for. So I need to change my passport name actually. Okay, so yeah. So I'm academically Zvastika Chatterjee, yeah. Okay. Then I would say we start in one minute, right? Okay. Nicola, one thing. Shall we do you have a prepared conclusion after? No, no. Eventually you should thank all the ICTP organization. I can say two words regarding the meeting after at the end, and you can eventually thanks the secretariat or technical staff. Very important. And we should say something about the meaning next year in Shanghai, right? Yeah, yeah. We should do it. Do we have a organizer here? Can you say something? No, I see Gong is not there, and I don't know who else is involved. Maybe you can prepare this a little bit in background, and we can maybe start. Okay. You want to finish this discussion? It's over. Thank you. All right. So then I would suggest that we start the session again with a talk about the composition of the Earth on Aqua given by Zvastika Chatterjee. Yeah. Thank you so much. So am I audible? Thank you. Okay. So the topic of my presentation has just said is composition of the Earth's inner core, and this work has been done in collaboration with Professor Tanushree Shahadash Gupta, Professor Sujay Ghosh, and Dr. Tilabdash. Okay. So the Earth as we know it today, it has an outer crust followed by the mantle. Then we have the liquid outer core, and then we have a solid inner core. Now, since the discovery of this solid inner core something about 80 years ago, its composition has been highly controversial. Okay. So however, there is a general consensus among the scientific community that the major constituent of the Earth's inner core is iron with some amount of nickel. So there have been several studies that have been done in the past, where people have tried to figure the structure of iron under extreme conditions of temperature and pressure. And most of these studies, they indicate that it is the HCP, the hexagonal closed pack structure of iron, which is stable at Earth's inner core pressure and temperature condition. Now, I would like to point here is that the BCC phase of the iron, on the contrary, has been shown to be dynamically unstable at Earth's inner core PT condition. So we are going to get back to this at the end of the talk. All right. So if I now calculate the density of HCP iron at inner core pressure and temperature and compare this density with the density obtained using geophysical models, then I will find that the density of HCP iron is higher than the true density of the inner core. So what it basically signifies is that there must be some light elements which are alloying with iron inside the core. And there have been actually several studies that have been reported in the past regarding this. And out of this, the several light elements that have been proposed, carbon is one of the most dominant candidates. So there are several reasons why carbon is one of the important candidates, as far as light elements are concerned. Firstly, carbon has very high cosmic abundance. So it is basically the fourth most abundant element in the solar system. It has a high affinity for iron. Then there are theories and formation of the Earth's core, which suggests that the Earth's core was formed out of a magma ocean which was rich in carbon. So there is a finite probability of some amount of carbon getting incorporated into the inner core in the process of this formation of this core. Now further evidence comes from the study of C1 chondrites. So what are C1 chondrites? They are basically thought to have similar elemental composition as the original solar nebula. So the original solar nebula is basically from where our planets have been formed. So if I study the composition of C1 chondrites and compare it with the composition of our mantle, we find that the carbon content of the mantle is actually lower than the carbon content of C1 chondrites. So what it indicates is that there is some missing carbon, and this missing carbon might have gone into the core when it was being formed. Further evidence comes from, once again, iron meteorites. So what are iron meteorites? They're basically thought to be parts of cores of planet decimals, which unfortunately broke apart instead of growing further into a proper planet. And the study of these iron meteorites suggests that there are abundant iron carbides there. In fact, the iron carbide in the form of Fe3C, that is cementite, has been abundantly found in these iron meteorites. So what happened is that previously in the earlier times, people thought that even our Earth, the core of our Earth is made out of Fe3C. However, this very belief was broken by First Principles calculations. So it was shown using First Principles theory that as you increase pressure, so if you cross 600 gigapascal, then it is basically, so basically iron carbide converts from a magnetic to a non-magnetic state. And this non-magnetic Fe3C, it has a much higher or stiffer bulk modulus as compared to the inner core. So what it suggests is that you can no longer consider Fe3C, that is cementite as an inner core constituent. So then after people were in search of better and better iron carbides and recent high pressure experiments, they suggest that Fe7C3 could be a possible candidate. Now once again, the structure of Fe7C3 is controversial. So there have been two different structures of Fe7C3 that have been proposed. One is orthorhombic, the other is hexagonal. So in the following, we will discuss the thermodynamic properties, the thermodynamic stability of, relative thermodynamic stability of these two phases, that is orthorhombic and hexagonal, and also look at its physical elastic properties. Okay, so when we are talking about thermodynamics, the very first parameter we generally talk about is the enthalpy. So enthalpy is simply U plus PV, so it's the internal energy of the system, plus the pressure volume work done, you have to create room for the system. So you basically have to push out the atmosphere and create the system. So it is basically the amount of energy you spend to create the system in the absence of temperature. So calculation of enthalpy shows that the enthalpy of the hexagonal phase is lower than the orthorhombic phase. As a result, we can see that in the absence of temperature, it's the hexagonal phase which is stable. Now we look into the dynamical stability. So here we have plotted the phonon density of states and the absence of imaginary frequency suggests that both the orthorhombic as well as the hexagonal phases, they both are stable, dynamically stable at inner core pressures. Finally, when we are talking about the inner core, it is very important to talk about temperature because you are basically considering very high temperatures of the order of 5000 to 6000 degrees. So the appropriate parameter to study here is basically the Gibbs energy. So Gibbs energy is nothing but the enthalpy minus the t times s. So you can basically borrow some energy from the environment and this assists you in creating the system. So greater is the temperature of your surrounding, greater is the probability of the system being created or easier is the system being created. So when we calculate the relative Gibbs free energy between the hexagonal and the orthorhombic phases, we find that once again, the hexagonal phase of FV7C3 is the stable phase. Right. So next we just discussed that apart from iron, we also have some amount of nickel, about 5 to 15 percent of nickel inside the core. So we next try to see what is the effect of nickel on the relative stabilities of these two phases and we find that nickel has no effect on the relative stabilities of the hexagonal and the orthorhombic phases. Okay. So having understood that, we also know that our overlying mantle, that is the mantle that lies just above the core is very rich in silicate minerals. So there is a very high chance that some amount of the silicon can diffuse into the core. Now apart from this, there are also other arguments which suggest that silicon could be a constituent, a major light element constituent of the inner core. So what we do next is that we try to dope silicon into this FV7C3 system. Now when we're talking about doping silicon, we can dope silicon in two different ways. We can dope silicon at the iron site as given by the first equation or I can dope silicon at the carbon site as given by the second equation. So when we do the thermodynamic calculations, we find that it is the second equation that is when silicon sits at the carbon site, that is the most preferred configuration. Now this is also expected because silicon and carbon, they are very similar charges. So it is easier for silicon to come and sit at a carbon site rather than to come and sit at an iron site. Okay. So once we know where to dope silicon, so we know that silicon goes into the carbon site, the next question that we ask is, will this system be stable or break up into HCP iron, silicon doped HCP iron plus iron carbide, this carbon? The question that we are asking here is important because there are several studies which have been reported in the past which says that silicon doped HCP iron could be the stable phase inside the earth's inner core. So once again, we do thermodynamic calculations, we calculate the Gibbs free energy of the left inside of the equation and the right inside of the equation, find out the difference and we see that no, fortunately our silicon doped Fe7C3 is stable. So it does not prefer to break up into HCP iron Fe7C3 plus carbon. Okay. So now that we understand that where silicon gets doped and whether our system is stable or not, we next go over to calculate the thermodynamic properties of this 3.2 weight percent silicon doped Fe7C3. So what do we see? We find that when silicon goes into the carbon side, once again the hexagonal phase of Fe7C3 is stable as it was in the absence of silicon doping. So this is where we have calculated the enthalpy and the pressure range is from 200 gigapascal which is the pressure at the core mantle boundary up to 364 gigapascal which is the pressure at the center of the earth. We next once again check the dynamical stability of all the four phases that is orthorhombic and hexagonal Fe7C3 and the silicon doped counterparts and we find that all four phases are dynamically stable at 364 gigapascal pressure that is at the center of the earth. Finally we calculate the Gibbs energy and this is where things turn interesting. So initially we found that the Gibbs energy of the HCP phase was lower than the orthorhombic phase or rather the HCP phase was more stable but upon doping of around 3% of silicon we find that this stability relationship just reverts. Now we see that instead of the hexagonal phase it's the orthorhombic phase which becomes stable at 364 gigapascal and in the temperature range of 2000 to 5500 Kelvin. Okay so apart from being thermodynamically stable what is most important to be one of the major constituents of the earth's inner core is that it should satisfy the physical properties of the inner core and when you talk about physical properties we are mostly talking about the elastic parameters why because all we know about the inner core is through our seismic studies seismology. So here we will calculate the longitudinal wave the poisons ratio and the shear wave velocity. Now to determine that we have basically determined the elasticity tensor to through stress strain relationship and then from that we could determine the bulk and shear modulus from which further gives us the relations for vs vp and sigma which is the poisons ratio. So when we plot the vp and vs that is the phase wave and the shear wave velocity of all the four phases that we just proposed and compare it with respect to the preliminary reference earth model which corresponds to the vp and vs of our earth which has been given by these square red squares we find it is the silicon doped fe7 orthorhombic fe7 c3 which matches the best with the earth's model. So in the case of vp the match is quite perfect whereas for the case of the shear wave velocity it's still a way off but still as compared to the rest of the phases it fares the best. Now if you look into the density of the phases of the hexagonal and orthorhombic fe7 c3 and the silicon doped counterparts we find that as compared to the preliminary reference earth model which is this blocks these squares red squares it is our silicon doped orthorhombic fe7 c3 which is closest. Now the difference we can see that it's there's still some difference the differences of the order of around 2.6 percent. Now the other two lines which is shown here which is the green and the green are experimental curves and they basically have been determined at much lower temperatures so the pressures are corresponding to the earth's inner core but the temperatures are as low as 300 degrees as a result they tend to show a better match however for our case it is the silicon doped orthorhombic fe7 c3 which gives the best match. Finally when we talk about the earth's inner core we all know that it has a very high poisons ratio so the poisons ratio is 0.4 for it's very close to the poisons ratio of the of a rubber so you can consider the earth's inner core to be like a rubber ball so it's important that the phase you calculate the phase you propose should also have a high poisons ratio so once again our silicon doped orthorhombic fe7 c3 has a poisons ratio which is close to the preliminary reference model as given here so it's around 0.421 whereas that for the earth it is 0.44 okay so to conclude we find that hexagonal fe7 c3 is thermodynamically more stable than orthorhombic at earth's core temperature and pressure nickel does not affect the stability field however interestingly introduction of silicon at carbon sites it makes the orthorhombic fe7 c3 thermodynamically more stable than the hexagonal now the calculated seismic parameters of silicon doped fe7 c3 are in good agreement with the reference earth model now next let us ask a very fundamental question is it possible in any way to achieve stability of bcc iron at earth's core pressure temperature condition so we all know it has been reported widely that the bcc phase of iron is unstable at inner core pressures this is what we obtain here so the black curve here shows the uh phonon density of states of bcc iron and the presence of these imaginary modes they suggest that the bcc iron is basically unstable at 364 gigapascal which is the pressure at the center of the earth now let us look at the effect of nickel doping so we all know there is some amount of nickel also there inside the earth's core so if we now dope the nickel in a very systematic manner that is we place nickel atoms at the centers of the bcc motors as in model one and calculate the corresponding phonon density of states of this model at 364 gigapascal once again we find that this model is unstable because we tend to get imaginary frequencies now let us try another way of doping nickel so now we have dope nickel in this model 2 such that the nickel-nickel distance interatomic nickel-nickel distance is not constant or rather nickel is somewhat randomly distributed in the lattice so in that case we find that the phonon density of states suggests that this system surprisingly becomes stable so what we conclude from here is that nickel doping in an appropriate fashion it can stabilize the bcc phase of iron-nickel alloy now we have repeated this experiment with various different configurations of nickel and also various concentrations of nickel so we have carefully distributed these nickel atoms in a random fashion in the bcc matrix and then we have calculated the phonon density of states at 364 gigapascal and we find that for all possible all the different configurations considered our configuration our bcc iron-nickel alloy is found to be stable at earth's inner core pressures as shown in this dispersion curve at this phonon density of states okay so now next if you really want bcc iron to be the stable phase inside the earth's inner core apart from dynamical stability you should also look at the thermo dynamical stability so we have to check whether our bcc iron-nickel alloy is stabler than the hcp and the fcc phases of iron-nickel alloy so this is what we have done here here we have plotted the difference between the bcc phase and the corresponding xcp and fcc phase as a function of pressure this is the enthalpy that we have calculated here and we see that for both the hcp as given by the black curve and the fcc as given by the red curve the bcc bcc iron-nickel alloy is more stable as as compared to the other two phases now as I just saved that for the inner core conditions it's important to introduce temperature so here we are talking about the Gibbs free energy so difference in Gibbs free energy also suggests that as compared to the hcp and the fcc it is the bcc phase of iron typical iron-nickel alloy that we are talking about which is thermo dynamically stable is thermodynamically also stable finally we look into the sound the phase and the shear wave velocities of the hcp fcc and hcp and bcc phases of iron now when we compare our vp and vs of bcc iron-nickel alloy with respect to that of hcp iron-nickel alloy we find that yes it's it fares better than hcp iron-nickel alloy in a sense that the vp of our proposed bcc iron-nickel alloy is off by around 1% whereas for the hcp it's off by around 10% whereas for the vs both bcc and hcp they fare quite badly when we calculate the density then what we find is that since the bcc phase of iron it since it's the packing fraction is quite lower then the density deficit that we were seeing in the case of hcp iron is lowered in the case of bcc iron so when we are considering bcc iron-nickel alloy the density deficit is around 3% whereas for the hcp iron it goes as high as 7%. So in case we are able to stabilize the bcc phase of iron then the budget of light elements in the core would be hugely affected. So to summarize nickel doping in an appropriate manner can bring in dynamical stability of bcc iron under earth's co-condition. Now the dynamically stable configuration of nickel doped bcc iron is also found to be thermodynamically stable with respect to fcc and hcp at inner core conditions. Since the packing fraction of bcc iron is lower it is expected to have a lower density than hcp or fc phases at inner core condition and the density of bcc iron is found to be greater than the earth's core density by about 3% indicating that once again there is a necessity of some light element to lower the density but however the amount of light element required or the budget of light element is usually reduced when we are considering bcc iron-nickel alloy. Thank you. Thank you very much. There are lots of questions unfortunately we have gone a little bit over time so let me summarize one question that several participants have asked and which is maybe to give just a few words about the technical ingredients of your calculations including the code you have used. Maybe you could stress whether there is any particular problem linked to calculations at extreme conditions like high pressure for example. Okay so the code we have used here is VAS and simulating high pressure is definitely not very challenging what is more challenging is to simulate high temperature because you essentially need to do a molecular dynamic simulation if you are really trying to figure out calculate things in a very perfect manner and doing molecular dynamic simulation for such kinds of systems which are huge and for example when we are talking about the bcc iron phase with nickel doping we have to give away all the symmetry of the system so it was a huge system with very low symmetry and then you have to go up to you have to simulate very very high temperatures so computationally it's very expensive. Okay thank you very much so I think that unfortunately we have to move on but as I said you have a lot of questions there's also a question of Shobana Narasimhan which is not in the chat so maybe you want to discuss this separately and I would say that we thank you again very much and we move on to the next speaker who is David Jakob. I cannot share the screen still. If you don't share your screen please thank you. Okay so David Jakob got more metal insulated transitions very good. Thank you Lucia. First I would like to to thank the organizers for giving me the opportunity to to represent our work here. This work has been done in collaboration with Gianluca Stefanocchi and the following I'm going to show you how to describe the mod metal insulator transition in the framework of steady state density functional theory. First let me briefly remind you of the foundations of density functional theory. DFT is based on the Humbacons theory which guarantees the existence of a universal density functional F which when you add it to the energy of the caused by the by an external potential by the external potential of a system and you minimize it with respect to the density you obtain the ground state energy and density of the system and I would like to emphasize that in principle this is an exact theory for the ground state energy and density. Now unfortunately the exact form of this universal functional is unknown and there we need approximations. Now a practical scheme to implement density functional theory is the Koncham scheme where you introduce a fictitious non-interacting system that reproduces the electronic density exactly. So this allows you on the one hand to decompose the universal energy functional into known contributions namely the kinetic energy and the Hartree energy and an unknown part which we call the exchange correlation energy and which we need to approximate and this is the starting point for a vast number of approximations. Now the other advantage of the Koncham system is that it allows to efficiently solve the DFT problem in terms of one-body Schrödinger equations which where you the Koncham wave functions are subject to an effective external mean field potential right and so this is all still exact in the sense that if we have an exact functional for the exchange correlation part then the Koncham scheme will give you the exact density and energy of the system right. Now it is very tempting and actually often done to interpret the Koncham wave functions and energies as a quasi particle energies and wave functions of the system and to use this to for example calculate band structures and spectra of a solid and yeah this is very often done of course we know that but one should keep in mind that this is not really strictly justified to do this right. There's no theory that states that that that the Koncham energies are actually the quasi particle energies and in fact there is there is a we know how to calculate exactly the for example the band gap and within the density function of the theory and this theorem states that the exact gap is given by the Koncham gap plus the so-called exchange correlation part and sorry sorry the derivative discontinuity is called so it's the exact gap in DFT is given by the Koncham band gap and the derivative discontinuity and so when you then use a very good functional so a very good approximation for the functional for example here exact exchange plus RPA for the correlation part and to use this for example for silicon you actually reproduce using this formula the exact or the experimentally measured band gap of silicon but you also see that the Koncham gap only makes up 50 percent of the total gap so this of course means that the band structure the Koncham band structure cannot be the actual quasi particle band structure right and this is even so this is even in for relatively weakly correlated systems like silicon like and other band insulators and semi-conductors the case right now the worst case scenario is when we actually have strong correlations so we consider for example the mod insulator and what happens there is that the Koncham gap is zero so the Koncham systems actually metal why and the the gap is completely given by the by the derivative discontinuity so so meaning that the Koncham band structure and spectra they can be sometimes good approximations maybe to the real band structures and spectra but one has to be really careful and cannot really take this always as the case now related problem is to use dft to calculate the density and current of a nano-scale system driven out of a equilibrium by combining dft with the nano-albitica formalism or also known as the non-equilibrium greens function formalism and the problem there is that dft is a ground state theory and it's not there's no basis to to use this basically in the in the case of a system driven out of equilibrium and my collaborator Stefan and Gianluca solved this problem by introducing a new flavor of density functional theory called steady-state density functional theory or brief idft and what you do there is you you introduce a new density namely the current in addition to the electronic density and and and then you have then you have a one back home like theory which maps which establish a one-to-one map between the density and the current on the one hand and the external potential and the applied bias on the other hand you can then also introduce a Koncham this system that is an interact a non-interacting system which reproduces the exact density and the current of the our nano-scale region and then well there you have for the Koncham system a one-to-one map with an effective external potential and effective bias so these are the Koncham potential and the Koncham bias of these systems and then you can within these you can then in this Koncham scheme calculate the density and current with the formulas known from the Landau-Bürtigarh DFT and the only thing that you have to count is that now the Koncham potential is a functional of both the density and the current and in addition you have a Koncham bias which also is a functional of the density and the current which and you can split this Koncham bias into the exact into the actual applied bias in the real system and then excited change correlation part now this was so Gianluca and Stefan have shown that this so with this approach you can actually describe Coulomb blockade and condo effect through nano-scale junctions which is something which is out of region for the Landau-Bürtigarh DFT now I will show you now how to apply this approach in order to calculate spectral functions with it so in order to do that we introduce a special setup and as basically it's an STM experiment an idealized STM experiment where we have our system S where that we want to have to measure the spectral function of and we have the tip electrodes which is weakly coupled to our system to be measured and now in the system in the limit of vanishing coupling to the tip electrode this our system is in quasi equilibrium and so that the derivative of the Lesser-Greens function with the voltage vanishes and therefore becomes the Lesser-Greens function basically becomes can be described in the equilibrium from the equilibrium density of states or a spectral function sorry and so so basically now also the density matrix of our system and therefore also the the electronic density it becomes the equilibrium one and if we apply now the Mia-Wingreen the asymmetric Mia-Wingreen equation in this ideal STM limit the DIDV then is we see that the DIDV is basically given by the spectral function of our system so so if we can now choose the coupling as we like we can couple it to some orbital some state in our in our system S only and then this approach this setup allows us to measure the spectral function of that orbital m from the DIDV and the idea is now to use IDFT to calculate the DIDV and therefore the spectral function so now in the ideal STM limit the equations for the density and the current decoupled right so so basically the the electronic density is given by the equilibrium spectral by the equilibrium conchamp spectral function so we so the first step is to to do a equilibrium DFT calculation simply and then once this is done we solve the IDFT equations for the current which is simply given now also in terms of the equilibrium conchamp spectral function but the only difference is now we have here the the conchamp bias and we have to solve this equation simply by root finding right because this conchamp bias also depends on the current right so we have to solve this equation additional equation and therefore we we will find the current as a function of V and we just take the numerical derivative and therefore with this equation calculate the spectral function of some orbital in our system yes and this actually works so we have here applied this to calculate the spectral function of the Andersen model and you see how when we crack up the the interaction in the Andersen model here we see this emergence of the Habat side peaks and the the narrowing of the condors and we can also derive an exact relationship between the many body spectral function and the conchamp one in terms of the exchange correlation bias and yeah and now the idea is that we could also apply this to measure the spectra to calculate spectral function of bike systems so to that we example we just call connect our so to say it Duncan STM tip to to a site in our in our bike material to a single side and then we can apply the asymmetric Mia Wingering we can still apply the asymmetric Mia Wingering formula right because we have now we we only have to assume that the tip is non-interacting electrons and therefore all this what we derived before still applies so the spectral function of that side is still given by the DIDV and we can make this relationship between the spectral function and the conchamp spectral function of the of of a single side in our lattice so so we want to apply this approach to to calculate the mod transition in in the Habat model and just to remind you what the Habat model is you have a lattice of a tight binding lattice where you have a strong onsite interaction you right and in this we know that in this lattice as you crank up the interaction the the quasi particle weight of the initially metallic system decreases until you reach a critical you where disappears and your the system becomes an insulator and our reference theory is a dynamical mean field theory which neglects non local correlations but it becomes exact in the limit of infinite dimensions and and gives a very good picture of the mod insulator transition in terms of the spectral function here you see initially the uncorrelated metallic state which there then the the Habat side bands start to appear as you crank up the interaction and finally the the the quasi particle peak disappears and we have the insulating state in so in order now to what we have to do now is to find approximations to our xc bias to our exchange correlation bias in that may describe this transition and first we have we derive from thermo liquid theory and dynamical mean field theory number of conditions that we can apply to to to our exchange correlation bias that our exchange correlation bias should fulfill these conditions so for example a very simple condition is that the xc bias at zero current should vanish and otherwise we would have a finite current at zero bias and this is unphysical and then from field some rules and for us that the the actual spectral density at the Fermi level and the Korn-Chem one should should be the same and this this implies that the derivative the first derivative to the with respect to the current should vanish at then additionally at half filling from the derivative of this expression we can derive also that the second derivative of the xc bias should vanish at zero current and finally we can relate the third derivative of the xc bias to the quasi particle weight of our to the quasi particle weight and yeah so we use this condition to to cook up an xc bias functional and so we know from our work in the Anderson model that in order to to produce a gap system what the what the xc bias needs to have is a is a discontinuous jump at zero current at zero current right which is actually related to the derivative discontinuity right and determines the the magnitude of this jump determines the the gap and on the other hand in the metallic phase this this discontinuous jump at the at zero current must disappear right so this is a Fermi liquid condition right the first derivative should vanish and so here we we make then this ansatz here for our xc bias where we have a part which is describes the mod insulating part of our functional on the Anderson model would be the column blockade part equivalently right which has this a discontinuous jump basically and times a pre-factor we call it the condo pre-factor because in the Anderson model it gives rise to the condo peak um which basically erases this discontinuous jumps and makes it flat right and we have this ansatz here which is acousticians and and this basically i to the three in here and we can derive then from the third derivative condition we can derive an expression for this lambda k parameter from and related to the quasi particle right so now we can basically we have a part for the Coulomb blockade part of our system which is relatively straightforward and this condo of pre-factor which is linked to the quasi particle weight which for example can be calculated in energy right um now we apply this and of course it works because we we have cooked up the functional that has all the conditions to describe the mod-metal insulator transition right so here this is for the beta lattice and you see a reasonably good correspondence between the exact energy spectra with the blue dashed lines here and our id of t spectra right um we can also apply this for the for example in the cubic lattice and again it works quite well apart from a few details which are not captured by our function right so to summarize so we can calculate many body spectral functions from id of t and the really nice thing is that this is basically a dft like computational course because you do just the dft calculation and equilibrium and then for each spectral function that you want to that you want to solve you just have to to solve this to find the root of the current equations in id of t and this is a very light computational overhead over conventional dft so in the Anderson model we can describe with this condo effect of Coulomb blockade and correspondingly in the bike systems where we can calculate the mod-metal insulator transition with dft friendly but that's all thank you very much thank you very much um we have uh one question from a participant and i think uh when we don't have questions from panelists so let us just take this one question because we are running over time again so please sushil Kumar can you ask your question hello hello hello hi we can hear you we can hear you am i audible yes please ask your question yeah can we calculate the electronic transmission function using id of t yes but the electronic transmission function in id of t is just an intermediate quantity what you really want to calculate is the current right actually the transmission function does not really mean mean that much in id of t because it's just the transmission function of the conchamps system right which per se is just a fictitious system and the the the the meaning of the conchamp system is just that to reproduce the density and the current through the device okay thank you okay thank you very much so um i think we have to move on and um um please can you unshare your slides yes and the next talk will be by uh Lorenzo monacelli about time dependence of consistent harmonic approximation hi so i hope you can hear me and today i'm going to first of all i want to thank the organizer for giving me the opportunity to be here in this online meeting and to tell you something about what i have done in my last years that is developing a new theory to compute nuclear quantum dynamics in strongly harmonic systems and this theory is called that time dependence are consistent harmonic approximation so in the beginning of this small presentation i'm going to tell you why it is interesting to actually compute nuclear quantum dynamics and then i'm going to introduce a bit what is time dependence of consistent harmonic approximation and i will going to spend some time telling you the linear response regime which is very interesting because many experiments can actually be simulated in this way so why do we want to simulate nuclear quantum dynamics so let's think about that standard experiment so uh usually we when we do an experiment we have our system at the beginning is an equilibrium at some temperature and then we send a probe to our system that perturbs our equilibrium state into some dynamical state and then we measure a property of this system that is changing in on time due to the effect of our probe we can formalize all these in a practical way so for example at time t equals zero we are at equilibrium which means we are in the ground state here i'm talking about nuclear uh dynamics so this psi here means refers to the nuclear wave function while all electronic properties are inside are encoded inside a Hamiltonian that is the bond of an amera approximation okay so at the beginning our system is in the ground state so it solves the static Schrodinger equation then we perturb the system through an external perturbation and our system is no more in the ground state of the total Hamiltonian which is composed by the polynomial Hamiltonian plus our perturbation and it evolves in time according to the Schrodinger equation and if we want to measure a property of the system that depends on nuclear opposition like a an observable a in this case we just need to bracket the observable with the wave function the time dependent wave function so uh this works very well when our system is in a pure state which holds at t equals zero but we can extend very easily even a finite temperature by replacing the wave function with a density matrix and everything remains the same apart that the now the equation we need to solve is no more the Schrodinger equation but the polynomial equation and this is how we compute the average of our observable so so far we have made no no approximation on the extent of perturbation so these include experiment even beyond linear regime so experiment uh very complex experiments uh if the perturbation is small and by small i means that if we turn off our perturbation the system doesn't change with respect to the equilibrium so at long time we will find again the same equilibrium as before so it doesn't change temperature or doesn't do any space transition or something strange then we can treat our system in linear regime and the good thing about linear regime is that how the full time dependence of our observable is fully determined once we know the response function that is a property of the system in equilibrium okay so our aim is to measure the so if we are able to simulate the response function we are able to simulate the outcome of most of experiments so this is the aim okay i can give you some example for example if you want to simulate the vibrations the nuclear vibration in solid if you want to measure them with infrared spectroscopy it's equivalent you can you can get the result of the experiment by simulating the dipole dipole response function if you want to to get the result of a Raman spectroscopy you need to simulate the polar visibility polarizability response function or if you want to get the result of x-ray or Newton scattering experiments you need the dynamical structure factor and the dynamic structure factor can be seen as a simple response function of this observable here so in the end what we really need to compute is response functions and how do we compute response functions so the way in which people nowadays compute response functions are these two ways here so we have an exact expression for the response function given by the cubo formula however solving numerically the cubo equation is really hard therefore we need some kind of approximation so what people usually do is the perturbative expansion so we for example use a auxiliary Hamiltonian that around which we do a perturbation theory this is usually for nuclei is the harmonic approximation and then we treat an harmonicity by expanding this equation perturbatively at the lowest order through usually this is done through Feynman diagrams the good point of this method here is that we fully account for quantum effects of nuclei if it is a quantum theory but unfortunately this theory doesn't work very well if the harmonic phonons are not well defined for example these outputs if our system is close to a transition and we have a dynamically stability at the harmonic phonons for the isometric phase and therefore in this case we need something else and one way the other way that is heavily used by people is with molecular dynamics with molecular dynamics we can compute the correlation function as also the talk of Professor Parinello talk a bit about this and if you have the correlation function you can link correlation function with the response function of the system through the two equation dissipation theory the good thing of this is that this equation here is exact in principle so you are including an harmonicity to all orders here you do not rely on an expansion or a perturbative expansion the problem is that this is classical molecular dynamics therefore you are neglecting quantum effects of nuclei we would like to have a tool that with the same computational cost of molecular dynamics for example is able to tackle is able to tackle systems with highly harmonicity so close proof transition and be able to account for quantum effects so this is the aim of the theory we developed so this is the aim of time dependence of consistent harmonic approximation this theory is an approximation as the name itself says and the approximate the only approximation we do is to constrain the wave packet to be a Gaussian so the most general Gaussian as possible then if the whole packet is a Gaussian we can derive a dynamical equation for the Gaussian very similar to the Schrodinger equation but the Hamiltonian now is depends is non-linear because depends have consistently on the wave function and this Hamiltonian is harmonic because harmonic Hamiltonians are the only one in which if we have a Gaussian it remains Gaussian during the wall the wall evolution and the parameters of these Hamiltonian are constrained so that we minimize the action so this is our first principle theory by meaning that it's not empiric we derived by first principles the equation and we can indeed also propagated not only pure states but also mixture of states so the theory works also for finite temperatures emulation by just replacing the wave function with a density matrix and the Schrodinger equation with the von Neumann equation for the density matrix so I'm going to show you a little video on how it performs on a one-dimensional system so here we have the external potential that depends on time the face of all your time here we have the dynamics of the exact way function that as you can see it is not a Gaussian while here we have the result of the time dependence that was is the harmonic approximation that is evolving and on the top right corner we have the average position of the nucleus in this during the dynamic as you can see the dynamics is very good actually and even if this is a point model it's not a realistic system I will show you this system in the next slides but we already can have the idea that the dynamics is really very good and this simulates a non-linear experiment that because we have dynamics even beyond linear but indeed as I told you at the beginning of this presentation linear response is a very important and very useful so we can derive the linear response theory and computer response function within time dependence and consistent harmonic approximation this is done quite easily simply we replace the density matrix that is the full dynamical density matrix to the equilibrium one plus a small perturbation and we linearize the von Neumann equation this is very similar to what we do usually also in time dependent density functional theory time dependent on the fork and all these time dependent self consistent theory and here we have a term that is the free evolution which is composed by two terms this term here is the free evolution according to the self consistent harmonic Hamiltonian with the equilibrium density plus another term the dates into account the fact that Hamiltonian itself depends self consistently on the density okay so we have harmonic evolution plus something that is non-animal okay and indeed the solution of the equation is very similar to time dependent Ft we have the green function here plus the perturbation applied to the ground state density matrix so for those of you who are familiar with time dependent Ft this is very similar okay and we can to get an idea of which processes which harmonic processes we are accounting with this theory we can derive the one phonon green function which is the green function of phonons is that if we apply a perturbation a dynamic perturbation on an atom is placing an atom in position b for example and we look at another atom inside the system in position a on how it moves this is correlation this response function here is the one phonon green function and in the framework of time dependent show we can prove that this is equivalent to a particular diagrammatic expansion where these straight fillet lines are the is the interactive green function of time dependence across its harmonic approximation while the dashed lines here are the green function obtained the harmonic like green function obtained with the self-consistent harmonic Hamiltonian with the equilibrium probability distribution so you can see that we correct this green function with a bubble like diagrams in which we include the three phonon scattering vertices and this bubble here is the two phonon the interacting two phonon propagator which is written as a Dyson equation itself where it depends on four phonon scattering sequences so this theory actually includes an harmonicity up to any order because inside the self-energy we have another Dyson equation which includes four phonon processes up to infinity and we have both three and four phonon processes inside so this is actually usually more much more than what done with imperfection and however if you look at the equations it seems that they are very difficult to calculate because in practice to compute the infunction we need to invert the propagator but if you look at this propagator here you can write it as a dense matrix n square times a square where n is the number of atoms in your simulation cell that can be hundreds or thousands also and you will immediately see that inverting a dense matrix of this big requires an order of n to the sixth operation which is a computational cluster or closer to couple clusters theory so I mean it seems useless and you do you need to leave this inversion for each value of the frequency that this is the cost of computing that diagrammatic expansion with distribution theory however likely since this is a set consistent theory the application of this linear operator dense big operator can be written in a much more compact way because it actually is a set consistent non-linear operator which is much smaller and so the actual application of l to a vector cost n square so we this means that we can use sparse linear algebra to compute the response function and in particular we can use an algorithm that in which we are inspired from time dependent DFT that is the largest algorithm to compute the response function this is inspired as I told you to what don't in time dependent DFT for example in the turbo launches approach maybe some of you are familiar with and so what what what this algorithm does is to find a basis in which this propagator the three propagator is can be written in a three diagonal form and if we can find this basis then computing the inversion computing the response function is very easy because we can write as a continued fraction function of this parameter so as you see it's new iteration of the algorithm which means adding a line on this matrix correspond to a new pole of this continued fraction so each time we do an iteration we apply a new pole until we converge I will show you a simple movie of this on high pressure hydrogen this is the green function as a function of this is the Raman spectrum the full line spectrum response function of high pressure hydrogen at 150 gpa and zero kelvin as a function of the number of step of the launch or as you can see at the beginning it's very ugly but as the number of steps increase the system starts getting something really nice and very similar to what we can get to an experiment and this is the full unharmonic spectrum of high pressure hydrogen and as you can see we have a strongly non-lorenzian shape of this vibrant peak here and you can reproduce system with a lot of unharmonicity include temperature effect and so on and just to do a comparison this is what you would get with the harmonic form so as you as you can see already from here it's really the effect of an harmonic in the system is extreme this is comparison of the experiment and data which is quite good and it's slightly different with the theory is due to the exchange correlation because this has been all computed with using for a full ab initio treatment of electrons so this is dft plus time dependent shuffle and the whole thing is that you can do simulation of systems with hundreds of atoms with full ab initio treatment of an electron and a complete also quantum treatment of nuclei so to conclude my presentation the time dependence of consistent harmonic process approximation is a theory to study nuclear quantum dynamics which can be applied in system with in which dynamicity is very strong in which perturbation theory fails and in system where quantum situations are very important are usually system with light atoms for example hydrogen or hydrates or a lot of lots of water and and ice and the computation of the linear response is really very efficient it's so efficient that actually is cheaper than ab initio molecular dynamics so you get quantum effects in your spectra with the same computation cost or even cheaper than ab initio molecular dynamics and to conclude I would say that the code is available is open source to this kind of calculation you can download at www.sstha.eu after this the website is currently under construction so we are uploading everything and probably the code will be available before the end of the week or early next week but so stay tuned and check out our website so that you can do this apply this in your simulation thanks for the attention thank you very much we have time for a few questions so I can see one question from a participant this is Seyed Vahid Hosseini okay okay I can see it is about spin orbit coupling you I can answer directly to the question so yeah you can include spin orbit coupling by directly on the density function of theory level so if you let me so the electronic degrees of freedom are all included inside the Hamiltonian so inside the external potential so everything that comes from electronic degrees of freedom can be accounted with the level of theory on which you treat electrons what you actually need to do this kind of calculation is a program let's say that you give an atomic position an atomic configuration and it gives you total energy and forces so if the level of theory you use to compute the energies and forces includes spin orbit coupling then you will include spin aspect of spin orbit coupling in your calculation within a bonon kind of approximation okay we have two more questions popping up one from Yehidu hello the computational cost of time dependence and consistency harmonic approximation is is uh is comparable with ammunition molecular dynamics it's a bit different because time depends this theory is highly parallel because what we do usually is a stochastic sampling of the wave function so we we and since it's a Gaussian our wave function we can extract the nuclei configuration very efficiently in parallel and then we need a way to compute energies and forces for this configuration this is the computational cost and this energies and for the number of calculations you need to do is more or less the same of ammunition molecular dynamics the difference is that since you have all together you can trivially parallelize them and indeed it is I answer also to the next question because it's very simple as I told you the the we can replace the wave function with the density matrix so then it is a theory at finite temperature so yeah maybe we have time to for the next question which is about this response function object maybe I could add something to the question in the sense that you introduce this response function that you have to invert maybe you can remind us the physics of of this object because here the question is I think whether you use the konsham potential as input so the konsham potential is somehow so once you arrive so when you do linear response theory you assume you so you already know the equilibrium so I've not spoke about this but actually there is another theory that is the self-consistent harmonic approximation without time dependency that is how to get that is very simple that is the equilibrium if you want the solution of our equation which is how to get the equilibrium density matrix once you get the equilibrium density matrix you have the equilibrium self-consistent Hamiltonian corresponding to that density and so you already have done all the calculations to compute the response function so you don't need to any other ab initio cost for your response function so everything you because then you already know know everything about who is this Hamiltonian and how is change when we change the response function so the konsham is within within within the the electronic part but the theory doesn't see it really thank you I think we have to stop here so maybe you can have a look at the other questions by the side and unshare your okay thank you very much again and we will move to the last talk now this is massimiliano stengel yes screening yep okay is it all right fine okay so good afternoon everybody first of all I'd like to thank the organizer for organizers for having me here it's always a pleasure for me to be back in Trieste even if this time I'm not really there but still I'm particularly excited this time because this talk will be about phonons and phonons are a long time tradition of the Trieste community so motivating phonons is not really hard it's actually pretty easy because phonons are involved in many physical properties of materials phonon frequencies are important for structural characterizations then they give information about lattice instabilities through imaginary frequencies that appear in the calculated phonon spectrum as you can see here in the example and then they can be integrated over the Brillouin zone to get thermal properties of crystals another important aspect of phonons is electron phonon couplings which give provide a whole new range of physical properties like electron mobility optical absorption and so on and so forth so in this talk I will mostly focus on phonon frequencies although the same arguments I'll be discussing will be 100 percent relevant also for electron phonon interactions I don't have time to go through this specific topic but if you want to know more about it you can find it in this reference so phonon frequencies are given together with eigenvectors via the usual Lagrange equation of motion which is written in terms of the dynamical matrix and of course the main ingredient here is the force constants matrix which can be calculated accurately and pretty efficiently at any cue point from first principles via the established methods of density function and perturbation theory but for most applications actually we want to know what is the physics behind the numbers that we calculate and it is particularly crucial in insulators to separate short range and long range contributions to the force constants so this talk will be mostly about this topic about the separation between short range and long range in a variety of contexts so let me start from the very beginning so from the theory from the fundamental theory that lies behind the long range character of the force constants and for this I'll take you through a journey back to the early 70s and in particular to this classic paper by Pico and Martin which has essentially pioneered a modern ab initio lattice dynamics in this paper the the force constants matrix as you probably know is written in terms of three as the product of three quantities two of them are the pretty simple are the bare point dipoles that are induced by the displacement of a nucleus and the third one is where the complexity of the electron-electron interactions lies is and is the screen Coulomb interaction so the remarkable result of Richard Martin's paper is that the long range force constants stem from the non-analytic behavior of the screen Coulomb interaction for small q vectors so to know to know better about where this problem lies let me write down the screen potential produced by an external charge which is this operator w which is written in terms of two ingredients which are which is the per Coulomb kernel which is the usual one of our interaction and we have the so-called irreducible polarizability that describes the charge response of the electrons to a screened potential and it contains the electron-electron interactions only at the level of the exchange correlation kernel so within the local or semi-local approximations to density functional theory we have that the irreducible polarizability is analytic and so all the non-analyticities can only arise within the bare Coulomb kernel and is what where we need to act if you want to operate this separation so the basic idea is to do a range separation in the Coulomb kernel so separate the a short range analytic part which is related to local fields and is usually a big operator that describes how electrons react to an external field on a short length scale and then we have the non-analytic part which is hopefully lives in a much smaller space and describes the non-analytic portion of the Coulomb kernel and this leads automatically to an exact separation between short range and long range and the long range part is written pretty much like the full the full dynamical matrix except that here hopefully we work in a in a very small space in a in a highly reduced dimensionality so it is easy to treat the material dependent properties that enter this formula can easily be calculated in this case our analytic function of q and our second derivatives of the total energy with respect to scalar potential or phonons and can be readily calculated within density functional perturbation theory so how does this work in practice let me give you the example of the three-dimensional case which is relatively simple in that only the head of the Bercoulomb kernel is non-analytic and this means that the long range potentials are macroscopic that is uniform on the scale of the unit cell of the crystal and the long range electrostatics is governed by analytic scalar functions of q which are the charge written in terms of the born dynamical charges at the lowest order and then we have the susceptibility chi which is depends on the macroscopic dielectric susceptibility of the crystal and so by combining all these ingredients in the formula I showed you in the previous slide we get to the popular to the very well known dipole dipole formula which is first established by Cochran and Coley and then this was the demonstration essentially by by Richard Martin and this governs the interactions between dipoles the dipole-dipole interaction which scales as the inverse third power of the distance between atoms so we try to see actually what happens this is only the leading order in the long range interactions and so it is reasonable to wonder what happens if we extend this order by adding a higher order multiples and this is what we did in a recent work where we considered the next lowest order which is given by dynamical quadruples and gives an interaction that decays as the inverse fourth power of the distance why do we consider this well the first reason is that we know how to calculate quadruples so why not check how they do in in calculations and the second perhaps more fundamental reason is that quadruples are related to piezoelectricity and actually this is a another classic paper by Richard Martin where he wrote the clamped ion piezoelectric tensor as the sublattice sum of the dynamical quadruples that now how many years later maybe 50 years later we know how to calculate well so given this relationship one would expect that the dynamical quadruples are important maybe important in the interpolation in piezoelectricity so to check this idea which we looked at rhomboidra barium titanate barium titanate is a perovskite in this case we consider the ground state for electric configuration where the polarization is along one one one and here on the left you see the calculated form by structure the red dash lines refer to the calculation performed just with the dipole dipole formula and the black lines refer to our higher order formula including dynamical quadruples and at first sight the result is pretty disappointing because you can see that the two curves do not well they match pretty closely across the the whole blue one zone except for a small region that I want to point your attention to here which are the acoustic phonon bands near the gamma point and with the standard dipole-dipole interpolation we get a persistent imaginary frequencies in the transverse acoustic band close to gamma while these frequencies completely disappear this is an artifact actually that stems from the fact that the electrostatics associated with the acoustic phonons is not described correctly unless you take quadruples into account so after this long introduction about the three-dimensional case let me move to the main topic of this talk which is can we do the same in two dimensions so two-dimensional the two-dimensional case this sketch by the way is a sort of very rough way of plotting the the the charge density perturbation associated with a polar optical phonon traveling through the crystal so in the three-dimensional case the fields are only modulated along the longitudinal direction because I have full crystal periodicity in the planes that are perpendicular to q while in 2d the the problem looks much more complicated because we have an extremanisotropy in the sense that the crystal is extended in plane but is microscopic out of plane and any polar phonon generates strongly non-uniform fields that propagate in vacuum for a distance that goes like the inverse of q there are pretty convincing first principles so this study by Francesco Mauri and the workers is a solved this problem of capturing these fields and in understanding how they influence long-range interactions between atoms but this understanding was based on two-dimensional dielectric models and we want to do actually something that is more comparable to what Richard Martin did for the 3d solids so derive these interactions for a fundamental theory so let's let's start from the beginning how what are these problems how these problems manifest themselves if i calculate phonons in a 3d code so let's take an isolated hexagonal abolon nitride layer in a vacuum supercell with a length l which is treated a converges parameter as usual and let me simulate an optical phonon of this type propagating along a some q vector in plane so this is the dispersion the red curve is the dispersion given by the 3d code which as you see departs quite strongly in a region close to gamma from the real curves which are associated to the isolated layers suspended in an infinite vacuum and this is precisely due to the stray fields i mentioned in the previous slides which induce a cross talk between images which become important when these stray fields propagate far away in vacuum and so forth small q vectors so this problem has been solved by the use of the Coulomb cutoff method which cuts the days of the Coulomb interaction along the z direction and this eliminates the cross talk between images and has been successfully applied to the phonon problem again by by Francesco and co-workers and this of course is the method that i use to calculate the black dispersion curves in the previous slide but the question is this works nicely for the actual dfpt calculations but how about the understanding of the long range forces that is what we are aiming at in this talk so for that sake let me write down the reciprocal space representation of this cutoff Coulomb kernel which is written in terms of the standard Coulomb kernel times some factor that depends on the exponential of q times the length l and what i want to point your attention to is that this factor makes the Coulomb kernel non-analytic at any gn which is the out of plane component of the reciprocal space vector this is a problem because i cannot just extract one part of the reciprocal space representation of the Coulomb kernel and just treat it separately because in this case i would have a huge matrix and i would have to calculate very complex response functions to describe the longer range forces and it is not what i want so we use a trick to work around this problem and the trick is to write the the short range Coulomb kernel as by using the image charge method so we take the original charge and we construct a vertical array of image charges with alternating sign and this is short range because if you look at a column of neutral charges seen from a distance this has no charge at all and so this is a representation of the short range part of the kernel and we call the longer range part the remainder what is this short range part sounds a bit exotic but you can see that it corresponds exactly to a phonon calculated in a three-dimensional code but just at the zone boundary along the vacuum direction and in fact the fact that you calculated the zone boundary introduces the 180 degrees phase shift between images which is exactly consistent with the alternating sign we don't get quite the correct curves yet but this discrepancy is now pretty small and this is how we define the long range part of the Coulomb kernel as the difference between these two curves so we can derive analytically this long range Coulomb kernel and it is written as in this formula where we have various ingredients we have the phi is 2d the 2d version of the scalar potential perturbation which is now non-uniform is mediated by hyperbolic cosine functions we have some range separation function which is a sort of e-walled filter and then we have the usual small momentum behavior of the 2d Coulomb kernel so what the the important thing about this slide is that the long range electrostatics becomes a 2d problem that involves only two by two dielectric matrices which is this matrix here so it is only marginally more complicated than the problem in 3d where we have scalar functions why hyperbolic functions well you can take the expansion of hyperbolic functions in a vicinity of q equals zero and you see that one the cosine the hyperbolic cosine is inversion even so it mediates interactions between multiples that are even with respect to a reflection with respect to the z equals zero plane and the opposite is true for the hyperbolic sine and the interesting thing is that unless the material breaks in version symmetry we don't have any cross talk between these two subsystems and so we have two components of the multiples at any order so this is the exact formula with in the mirror symmetry case as I said we have two components inversion even and inversion odd the inversion even part is consistent at lowest order with earlier models and the inversion odd terms are new and describe the interaction between out of plane dipoles which was missing in earlier theories we have some other features like an introduction of dynamical quadruples but let me point out that this is an exact formula which is we can push its expansion up to an arbitrary order in the multipolar expansion in the present work we stopped at the dynamical quadruple level so regarding the results this is what we get for boron nitride for the same mode as before we don't get qualitative huge qualitative improvements over the state of the art which is the work by so yeah at all we get a slight improvement due to screening and the quadruples are surprisingly ineffective at correcting the phonon frequencies in practice they are only relevant for the off diagonal elements so what we get actually is maybe the the original result which is a linear dispersion for the zio phonon branch which is not captured if you don't take these out of plane electrostatics into account so I'm already over time I leave you the summary for you to read and I thank you for your attention thank you very much you are actually almost perfectly on time and we have a question concerning the origin of the sharp phonon kink at the gamma point maybe you can just say a word about the physics of this sharp kink uh do you mean in the longitudinal branches like this one or in the zio phonon which my question is the question of said by rossini yeah and I guess it refers to the earlier slide but yeah okay so this this kink is due to the fact that the longitudinal phonons carry some macroscopic fields so in 3d you have an analytic behavior non-analytic behavior at q equals zero that gives you an lotio splitting but in 2d you don't have an electrostatics is weaker so you don't have an lotio splitting but instead you have a linear dispersion of the longitudinal phonon approaching gamma and this is uh this was a result that was known earlier okay we have a new question also uh from mod or by duraman who maybe wants to explain the question his or herself are you online so it's about the the fact that you're only talking about optical phonons if I get the question correctly so you want to know about the the question is about the acoustic branches so I actually talked about the acoustic branches in the context of the bokeh barium titanator rhombohedral barium titanator 3d solid uh in the case of 2d materials uh acoustic branches are also important uh but we found that uh the electrostatics is not did not play a substantial role in the in the systems that we simulated possibly in some systems it is important and if you find some issues like some imaginary frequencies like near gamma that we would we would be happy to know uh because it might be a candidate for our uh higher order interpolation method I think um to have ample time for the poster discussions we should stop here there's one one more question in the chat that maybe you can ask answer directly and thank you very much again I thank you all the speakers of the session I wonder whether the organizers want to say something now yeah yeah Lucia thanks so I think that uh I should add so I'm Francesco Mauri and so I would like just to thanks everybody uh and we'll go in detail of the tanks later uh so when we so we were uh thinking to postpone by two years this this uh this um this meeting because of a covid uh crisis uh at the end uh also thanks to I think also Nicola I I think was very enthusiastic in trying this kind of online conference and our worry was to involve the people because I think the total energy beside being a conference is a is a moment where the community meets showing the advances and in a community manner and so in doing the trying to reproduce with this online conference such a community participation we decided to open a a call for a contributed co-talk and poster and I think the response surprised us we have many very good offers and we managed to do a very nice uh program out of your proposition so so I think the experience could be improved but I think I was pretty positive and uh and I think we should take this this as also a message for the uh I hopefully physical conference in a couple of years that perhaps adding also some contributed talks to the invited talks so I would like to finish thanking all of you so I start thanking the participants so we have over 1000 applications uh more or less 700 people connected so much more than usually with total energy we have 300 250 people in physical presence so the participation was very large for many also countries we have many from India for example and so we thank the people that offer for contributed talks and for poster the key speakers I would like to thank also the ITCP organization I think Nicola will tell something more on that part and I thank the two organizer my co-organizer so uh Tanusri and Nicola and especially Nicola because I think 90% of the job was done by Nicola so I thanks a lot so we just help him in uh choosing the the the the the the speakers and organizing the outline so giving the idea but I think uh really Nicola put a lot of effort and the the success is also uh thanks to him so now with that I conclude we'll uh next year we'll have in in China a mini and hopefully in presence and also in two years the next uh Maxi okay thanks also obviously also all the people that help sharing the the session okay I just want to I just want to echo uh with Francesco that Nicola really deserves a very big thank because as Francesco has said bulk of this actually goes the credit goes to Nicola okay thank you thank you very much but uh this is not true because actually 90% of the work was behind the scenes so we really have to thank the the ACTP people so the sigma team Victoria Adriana and Monica and the IT people of course Massimo Sabrina and uh Valter that really really did 90% of the work I think so thank you very much okay and don't forget maybe at the poster session now so we go yes okay yes we are doing the the the greetings now because this is the last common session but yes we will be spread around the the poster rooms in 10 minutes no so and I hope to see you next year in the mini total energy in Shanghai and hopefully in two years here again in Trieste thank you okay so thanks everybody bye bye thanks bye bye this was very nice yes yes I I still find it a bit it was little time for discussion but yeah the discussion time is short but uh well I mean one has to learn right yes nice and to see how I mean nice physics going on and I mean and honestly even in when when we have