 Massimo, se mi avvisi quando siamo online. Siamo online. Siamo aperti e adesso lo sta partendo lo streaming. Un attimo di... Ok. Ci siamo. Quando vuole? Buongiorno a tutti. Eh... Quindi... vediamo più di 1-2 minuti, perché vedo che ci sono persone still joining the meeting. In the meanwhile I just recapped our ideas that don't... just let's keep the zoom chat clean while the speaker give the talk. So don't write anything in the chat. Keep all your questions in your mind. At the end of the talk we will... I think we might have time for one or two questions on voice. So raise your hand and we allow you to the microphone so you will speak. Or... and even if you have other questions or you don't want to speak or we don't have time to speak write in the zoom chat. For people who are on the Slack channel please use Slack because it is much easier for us to answer and even if we don't... in case we don't have time during the session to answer we can still answer on Slack. So on Slack is sure that for sure you will get an answer. And now in the session now we will give priority to questions coming from the streaming. So at the end of the talk there will be room for one or two questions on voice from people here and other questions from streaming and all participants can write on Slack every time. So maybe it's 8.33 I see we can start now. So I introduce Stefano the Gironcoli who today will speak about advanced functionals. So higher accuracy functionals, hybrid functionals, strongly correlated materials and weekly bound systems. So please Stefano let's start. Ok good morning good afternoon everybody I share my screen yes you should be seeing my screen and let me ok I hope you can see my the first slide of my talk ok so this lecture is about exchanging relational functionals that is a vital ingredient of any DFT code but everybody knows there are many many functionals that are available and the most there are several generations several let's say groups of functionals of increasing complexity starting from LDA then generalize grid and corrected functionals and then meta-GGA and then more and more complex functionals. So I will very quickly just remind to everybody the main ingredients of the simpler functionals and how they perform and then I will move actually and tell something more about the more advanced functionals that will be showcased in the hands on session afterwards. So let as I said there are several generations starting from the simplest all of them have let's say the simplest only contains information about the density this is the name of this local density approximation then if you want to include some more in homogeneity coming from the system then in the 90s the generalized gradient approximation was developed by the U and co-workers and complicated functional in which include also kinetic energy densities as ingredients have been proposed giving rise to the so-called meta-GGA methods but also hybrid functionals and then we will talk about also the Vander Waals interaction later. For the most basic normal let's say functionals like a DA GGA there are a number of resources that one can find there are millions of books about the FT this is just one review early I mean some early review of the properties and features of LDA and I think it was a good review back then and I think it's also still give some basic understanding of the method. So, Koncham and Hoenberg-Kontheorem and Koncham in a nutshell Hoenberg-Kontheorem say that given a density you can write in principle it exists a general function that will contain all the information that is of the system that is specific of the electronic system not the interaction with the external potential so contains the kinetic energy and the interaction energy of the system Koncham introduced the non-interacting fictitious system that contain the information about the kinetic energy of the non-interacting system and using this let's say tool as a this ingredient as a tool one can rewrite the interacting Hoenberg-Konfunctionals in terms of three terms two terms that are known the kinetic energy of the non-interacting system via the Koncham construction and the classical electrostatic energy of the density that is called the heart of the energy and what is left which is therefore defined by this difference by this equation is the so-called exchanging correlation function and given this definition then the functional becomes the explicit functional for the total energy of the system can be written adding the coupling with the external potential and from that one get the Koncham equations which basically are simple to solve is a self-consistent set of equations as we saw the previous days and it's a simple mean field like approach with the advantage that it is formally exact the problem is that the exchanging correlation functional is not known exactly and so one has to actually make approximation the still simple approximation to exchanging correlation functional are possible the reason why are possible the term that they tried to mimic the term that they want to describe is relatively small if you take a system this is just an example for a manganese atom isolated manganese atom and if you break down the energy contribution you have the kinetic energy the direct interaction of the system with electrons with the external potential which is the highest the biggest contribution then the kinetic energy second the Hartree valance valance electrostatic interaction is next in size and the exchange and correlation would be even smaller are just a small fraction of the total energy so the hope is that the reasonable approximation for this relatively small part of the energy of the energy functional will still be useful and so the simplest approximation is the one that goes under the name LDA or local spin density approximation actually if you include the possibility of dealing with magnetic systems and the idea is simply that you have a system and every little tiny volume of the system could be let's say one can hope that the exchange and correlation energy associated to that small volume is the same or similar to the exchange and correlation energy that you would have in a system which is homogeneous everywhere and has the density that you find in that point so basically you can extract from one single system the homogeneous electron gas as a function of its density the energy exchange and correlation energy per particle and then replace this model functional everywhere in the space everywhere in the volume of the system that you want to describe and this in the end is not too bad I mean the LDA and local spin density approximation can describe a number of properties the DFT is a ground state theory so it is if it is if you use a good approximation should describe a number of properties that are related to ground state to the ground state of the system and so the energetics the phase stability different thermodynamics the equilibrium geometry we saw yesterday how to relax the system and many other properties that you can obtain by applying perturbation to the system so you can calculate apply electric fields calculate static electric constants calculate stresses compute piezoelectric constants or elastic constants etc etc so there are a number of quantities that are that are accessible that are within the the reach of DFT and if you apply a diabatic approximation then also lattice dynamics lattice dynamics and many properties other properties like in general termodynamic properties where you address the temperature issue for the ionic degrees of freedom but you think that the electrons are basically in their ground state can be addressed LDA is good I mean is not perfect this is why we don't better function but provide good trends chemical trends so one can compare these are historical historical slides let's say look at the first ionization energy and you see that LDA and local spin density in particular the spin the possibility of addressing different spin component to get nice trend that are parallel the experimental the experimental results this is for SP and transition metal elements then you can also calculate excitation energy from one angular momentum to another and for instance in transition metals you can look for how much the relative stability of the diffused electrons the S electron in the S shell and the localize the orbital so many things are qualitatively correct not always quantitatively correct but the picture is pretty was pretty successful and actually the fact that we are all using DFT is related to the fact that overall even the simplest approximation is not too bad and then people studied different phases and again one get results that are in reasonable agreement with experiment and then the problem is maybe describing exactly the correct transition pressures and in general the results are reasonable so these are just as an example the first row of dimers and comparing three quantities that define the binding curves the minimum position the binding energy so the depth of the of the binding well and the curvature around the minimum that is associated to the frequencies and typically what one obtain and this is then true across also for more complicated system the geometry is not too bad is pretty good a few percent typically and typically the LDA bone length is maybe a little bit underestimated the binding energy is overestimated instead so for instance boron 2 is bound by almost four electron volts in LDA in local spin density while it should be a little bit more than three in experiment oxygen is strongly overestimated binding energy of oxygen is strongly overestimated so there is a general tendency of LDA to make bones a bit shorter a bit more bound than they should and correspondingly a bit more stiff than they they would be in the experiment so this is a general let's say trend to which I summarize here so lattice constants are typically good a little bit too small typically binding energy or cohesive energies are too high and correspondingly then the stiffness the response function is a bit overestimated band gaps so expectation energy estimated from difference of eigenvalues are typically too small but this is not a defect of LDA this is actually a feature of DFT which are completely understood that this is not something that if you stay inside LDA you should expect to get right so there are other methods like GW that can address the problem of the excitation energy and time dependent DFT can address some of these some of these problems but let's say normal ground state DFT is not meant to describe properly the band structure the band gap of insulators so one should not be too harsh let's say on the function as if you don't get the the band gap right but I mean LDA was a very nice very successful method approximation for its complexity I mean for how complex it is but people got an IP I mean people wanted to have something better and so in the years the generalized gradient approximation were developed in which in addition to the local density you include an information about the local gradient local reduced gradient so this is just again a simple picture of the density of a beryllium atom super simple system in the radial coordinate so for an isolated beryllium atom and what you see is that the density decay when you go away from the look and the gradient of the density divided by n to the fourth third that is a dimensional quantity actually diverge because the gradient decay exponentially but also the density decay exponentially and since the power is fourth third actually this this ratio diverge so basically you see that the region where most of the electron live has a relatively small gradient so the green curve is more or less bound but as you move toward the surface the density decays the gradient may become bigger and bigger so the these functionals that include a dependence on the local gradient may treat differently what happens at the surface of the external region of a molecule of a crystal with respect to what happens in the more densely densely and homogeneously described regions inside the molecule and so being able to describe this different different situation in different in different way with different portion of the function allows to get a better description so this is the hope and this is why all the gga were developed the historically the one of the first is pw91 per the one functionals introduced in the early 90 more recently more nowadays basically people use rather pve and pve actually is nothing but another functional that has been written with the same ingredient of pw91 but trying to simplify the functional form so pw91 and pve are expected to be actually very similar they are basically the same functional but with a different with a different functional representation but the ingredient in them are the same and the typical ingredient is written as a term that is the ldaa contribution coming from the uniform electron gas times an announcement factor depends on the density through the rs bigger size radius which is related to the density the magnetization if you take spin polarized version of the theory and the radius and this functional this announcement factor as a function of the radius so it is goes to the limit of the homogeneous electron gas when s goes to zero but has the ability to to have different behavior when s is bigger and while the ldaa limit is well defined because the ldaa limit is what you get in the homogeneous electron gas limit so for large value of s is somehow not bound is not constrained by or is much less constrained and so this is why there are millions of gradient corrected functions there is only one ldaa because there is only one energy external correlation energy of the homogeneous electron gas if you compute it correctly there is only one function but there are actually many ways in which you can explore the largest limit but the one that is historically in practice per the 191 PBE and PBE sol which is some revised PBE so there are a few functionals but they are there are millions of functions that are being proposed but let's say the one that people use in practice are a few what these functionals the reason why they were they were introduced was to solve let's say the over binding problem of ldaa so here is just a table of small molecules and the last of the atomization energy of small molecule how much energy in kilocalpher mole takes to take a molecule like H2 and break it and for instance H2 hydrogen dimer cost 109 kilocalpher mole to break and local density was actually over binding so it was giving a number which was significantly higher the new functionals improve in the sense that they give lower binding it is 105 in this case but if you look for instance the oxygen yeah the oxygen case or the nitrogen case for instance nitrogen oh sorry was 267 in ldaa which is much more than the 229 in of the experimental results and moving to pw91 or pbe improve it and it gives 240 so there is a general improvement in the energy so that mean absolute error in ldaa was for this set molecules like 30 this error was reduced was cut by a factor of 4 on average by moving to the gradient corrected functionals and similar results can be obtained for other properties for solids let's say magnesium so the lattice parameters are typically good the the lattice parameter sorry no this is the sorry the energy of adsorption I think and no this is the lattice parameter sorry the lattice parameter as I said for ldaa is typically slightly underestimated and this is improved in a little bit larger lattice parameter in gga but overall the energy is much better so if you put the two met together ldaa is good the lattice parameter in both the lattice parameter are very good in ldaa tends to over bind gga gives better cohesive energies and the bike modulus in ldaa is typically slightly larger than an experiment gga has a bit on the other side tends to be a little bit less bound bike modulus elastic constant slightly lower gga is very important for magnetic systems in particular one of the first examples that made let's say the success of gga is that if you do local spin density calculation iron normal iron is come out come out to be non magnetic and instead gga correctly give the magnetic ground state so that was one of the first success of the gga function right so ldaa tends to predict bonds a little bit too small gga tends to give always small error sometimes a little bit larger sometimes a little bit small so there is no more no not too systematic and but typically errors are pretty small and ok if you calculate compare different phases like FCC and then you can calculate for different elements which one is more stable and the prediction using gga is typically accurate so gga really improve the the the energetics both for let's say structures that are very different like bcc and FCC but also for more like the the relative stability of close packed structure like FCC and exactly close packed so yeah they predict typically if you work with the gga you get pretty reasonable picture of your system and as I said iron in LDA was wrongly predicted to be non magnetic but in gga this is sold still and this we come to the to the topic of today there are there are errors there are problems with these functions overall they are fantastic approach but this doesn't mean that there are no issues as we said as we saw one of the main motivation for introducing gga was to improve the energetics to get more more accurate energies and if you are doing thermochemistry or anyway if you compare different different phases you would like to have accuracy a high accuracy on the on the computer energy and we saw that LDA was had a big error like 30 kilocal per mole gga improve and reduce it to say 8 still chemical accuracy that is what we principle would like to obtain which is energy accuracy within the kilocal per mole is still far and trends bond lengths type of bonding are accurate but typically I mean if you have covalent, ionic bonds so the strong bonds are typically well described but weak bonds are not very well described typically there is another problem that if you consider the system of electrons and you look at the how the electron electron interaction is described we all know that Hartree-Fock is defined by requesting that the proper symmetry of the wave function is contained is described and so this actually define the energy the electrostatic energy the Hartree energy and the exchange energy in Hartree-Fock is built from the Hamiltonian and it is by construction cancel the self interaction that is included in the in the Fock operator cancel exactly the self interaction that is spuriosly included in the Hartree-Fock this is exactly in Hartree-Fock but Hartree-Fock does not have correlation and so for other reasons Hartree-Fock is not desirable but for this feature Hartree-Fock does the right thing so the fact that one electron should not interact with itself LTA and GGA replace the Hartree-Fock exchange energy with a functional an approximated functional that therefore does not cancel exactly the self interaction in the in the system and this in some system where the self interaction is important may give problems therefore may have problems when you look at molecular dissociation limit because when you dissociate the molecule then the fact that the electrons I mean you deal with smaller number of electrons and then the details of how these electrons interact becomes more important or when you have localized orbitals of the self interaction is more stronger when the involved orbitals are more localized and so in transition metal oxide in materials involving rare hertz or in other in any other system in which the atomic the atomic feature of the of the atom survive in the solid in all these systems then the discreteness of the of the atomic system is more important and the fact that the self interaction is only dealt with in a kind of average way is that more evident when you have many many electrons many so the fact that maybe you are not treating the self interaction super correctly is not so a big deal because anyway you have many many contributions but when when you deal with the system which a few electrons so are bound to a transition metal then what happens to those few electrons is more important and so the self interaction cancellation maybe a problem we will deal I will talk about that in a moment van der Waals interaction is another issue that is not well described in the in the normal usual LDA no by the way I mean van der Waals is a correlation effect is due to the fact that you may have local fluctuation in the density and and this local this local fluctuations binds can give a contribution to the binding energy of the system and however this is a non local correlation effect I mean fluctuation in one part of the system internet fluctuation in another part of the system and their interaction gives you a contribution of the energy LDA and GGA focus on what happens around one point in the system and therefore do not account for this non local term and so this for a while this van der Waals interaction where since they are relatively weak this van der Waals interaction where neglected and for many system for many properties this is not a bad approximation but van der Waals is everywhere everywhere and especially in organic molecules it may be an important part of the binding energy so in more recent years there have been a lot of effort in trying to account for it ok so let try to then see what options one has to deal with this problem so these are the two main problems I want to talk about in the rest of the presentation and it's what you do what you could do to address self interaction correction self interaction cancellation and what you could do to address van der Waals ok so the first one actually self interaction error uh has been dealt with in various ways with a method that was originally proposed by Zunger and Perdyu Alex Zunger and John Perdyu called self interaction correction then hybrids, functionals and the DFT plus method so we will go quickly through the self interaction correction because it's not implemented in quite an espresso but historically I think is important so the idea is as I said that in standard DFT electrons move in an effective potential in which electrons interact with an average the average density of the energy, of the potential plus an average potential where instead the actual system each electron would see the other electrons where they are and not their average distribution if you have a very broad band, if you have a very a very diffused system this fact, the fact that the you interact with an average of the other electrons including actually some contribution from yourself is a relatively small small term the broader the more diffuses the orbital the smaller is the self interaction and so the smaller is going to be the self interaction instead instead if you have if you are dealing with orbitals that are strongly strongly localized then the fact that you are means within the self interaction is more relevant and so this is why self interaction error is more relevant in transition metal oxide or rare hurts or whenever the orbitals are very very localized for instance and in particular you see problems when you change the locality the localization of your orbit because that the fact that you are comparing two situations in which the localization is more or less the same even if you make an error you may hope some compensation between the two systems but when you let's say take magnesium oxide manganese oxide and you interact with lithium lithium is a metal and so the electrons of lithium are delocalized strongly delocalized when you combine it and you make a lithium manganese salt or lithium manganese oxide then the lithium gives away its electron which is strongly bound to manganese in a localized orbital and therefore the localization properties of these electrons is completely different and this may make a big error so the gga give you a formation energy or reaction energy for this reaction of 3.6 electron volt while experimental is 4.1 and for the iron phosphate it's 2.8 versus 3.5 for vanadium phosphate is 3.3 versus 4.6 so you see that you may have big errors error of the order of electron volts or more because the localization of the orbitals is very different on the two sides of the reaction and therefore the fact that the interaction is not well represented is is bad what is the problem the problem is that when you have localized orbitals if you make a model, if you make a tie binding model you have localized orbitals and then you have op-in terms that give you the scribe the the op-in, the the localization of the electrons across this levels this op-in term actually comes from the overlap between the orbitals and so this is actually something that the LDA or GGA described correctly because the wave function the shape of the wave function is quite a robust properties what is wrong is the electrostatics because if you have two electrons in two separate orbitals at some distance this there is one interaction if you bring two electrons on the same site of localized orbitals of course the electrostatics is very deep but at the LDA level you describe actually a band structure of one electron moving in an average potential and this average potential should describe the two situations when the two electrons are far apart and where the two electrons are together and these two situations are very different and the LDA is not able to describe them correctly so the problem of this self interaction correction this self-energy error is not that the op-in term is wrong but that the specific interaction the fact that the electrostatic interaction strongly depends on how many electrons you have at that particular moment on your atom is the problem for instance if you do calculation for iron oxide and you do simple GGI calculation this will be one of the example in the in the unsung session the simulation fill the several levels and here is a plot of the density of state project density of state something that we saw yesterday and basically you have the majority spin the blue curve are the the orbitals of the majority spin for iron which are completely filled but iron has six electrons and so it has one electron in the minority spin and this electron in the minority spin is goes in the red states but you see that these red states are broad and only partially occupied by the single electron that is there but what you see is that from this picture the iron oxide is a metal has a non zero density of state at the fermi energy which is Mach 20 here this is wrong I mean this is what GGA leave but this is at variance with what experimental is found that iron oxide is actually an insulator and this is due to the fact that the self interaction energy is means represented for the localizer with us over over so the methods that have been proposed to overcome this problem by then originally proposed by and and then also proposed by so and so on and or the TFT plus methods so we will very briefly discuss the self interaction correction method and then we will go to the hybrid functions and the FFT plus you that are implemented in the quantum space ok this is the original the original paper and so it's I think good to have a read if you can it's all paper 81 but the idea is that the functional that describe sorry kinetic energy plus R3 plus the functional extension correlation functional except that the approximate functional that we typically use does not does not properly describe the self interaction and so they basically say let us add some term and let's make this term such that when we have one orbital only when orbital is when we have only one orbital that is occupied let's make that this correction give us the result that there is no self interaction so it's basically a dock solution for restoring what you would have naturally in Artifox it's a bit complicated because orbital dependent it's conceptually is helps and one can show that it solves some of the problem but it's a bit complicated so this is not implemented in quantum espresso but I think is useful to know and read maybe this paper to have an idea the in order to address this problem let's say in the in the with the aim of addressing the typical problems in DFT people have introduced various methods and one one important formal concept that has been developed and has been very fruitful in the development of functional is this the adiabatic connection formalism we saw that we have Holmberg and Kohn theorem that say that we can define a functional f that contain kinetic energy and interaction we are seeing that we can have a no feticious non interactive system the conscious system in which you have a kinetic energy functional and plus potential the conscious potential external potential so actually there is a whole family of and so these two things the conscious sorry the original Holmberg functional and the conscious functional are the same except that the conscious functional is specializing the case in which there is no interaction so basically whatever is the interaction that you choose you can define a Holmberg Kohn functional f in the special case in which the interaction is zero you get the Kohn functional and you can think that you can introduce a parameter lambda that switch from zero where you have the Kohn and you know everything to one where you have the actual system that you are interested in which is the interactive system and you may have to do some approximation there but in principle there is this parameter that switch continuously adiabaticamente from one to the other ok there are I mean I will not enter in the details of this procedure but basically there are variational principles for every value of lambda there are variational principles and so you can calculate the derivative the derivative since it is variational it allows a kind of element theorem to be used and so eventually you can from this formally write the term associated today to the interaction as the sub of two terms which in DFT we call Hartree and exchanging relation functional and this is nothing but average between going from non interacting to interacting of the expectation value of the interaction on the ground state way function correspond to that interaction style this looks all complicated but this is just to say that in principle one can if you have an approximation for different if you have good approximation for different values of lambda you could actually make a refined approximation in which you make an average of these approximations this was the idea that basically put forward in the early 90s in which you said ok we have this object that has as a function of interaction strength from go from non interacting to interacting in the and let's say that it has a certain dependence let's be hope that it is smooth and so that we can approximate it with a straight line and then the value of the integral is the rule is one half the value in zero one half the value in one and the nice thing is that the value in zero is actually Hartree-Fock I mean is the expression for lambda equal to zero is something that you can explicitly formally write correctly is the Hartree-Fock energy that is computed with the what is called the exact exchange in DFT and for lambda equal one you may need to do some other approximation like for instance LDA but then you combine one more accurate description that you can compute that lambda equals zero with some approximate description write at lambda equal one and this was the half half functional proposed originally by Beck you can at that point once you open a gate then the herd goes out and so in the next years many other functionals were proposed in which you have a fraction of exact exchange so alpha zero like 20% of Hartree-Fock exchange plus a bunch of linear combination of LDA, GGA functional for what concern the reminder is very very popular and basically is 20% of Hartree-Fock plus some gradient corrected function Purdue and co-worker also played with this idea and eventually they found out that actually they write they had some argument to say that the right fraction of Hartree-Fock exchange that should be included actually one quarter is not 20% but it's 25% and for the reminder they propose to use PBE basically so this is PBE zero and a variant of PBE zero in which you also decompose a short range in a long range part mix differently with some empirical parameter that however is more or less fixed is this head scuseira hand of functional, HSC which has let's say a parameter localization parameter omega make one fifth of an angstrom if I remember correctly or actually five angstrom I mean the inverse of of five angstrom and the fraction of exotic change is again one quarter so you see that the ingredient that enters Fock DFT in this hybrid function in which you mix a fraction of Arctic Fock and some other GGA like functional is the fact that you have to do the Arctic Fock calculation I mean you have to be able to compute the the Arctic Fock exchange unfortunately and this is something that can do all the time and the many codes in the chemical community that is based on localized orbitals do implement do implement hybrids hybrids or let's say Arctic Fock energy in a plane wave setup is difficult it's complex but it's in particular very expensive so yeah so this is the expression for the Arctic Fock exchange energy okay some over double sum over occupied bands for each pair of bands you compute a auxiliary density V of K prime V prime star times V of KV and then you calculate the interaction energy between these auxiliary densities and you sum over all of them this is formally well defined I mean it's written there but computationally is expensive now but this is the ingredient that you need if you want to do an Arctic Fock calculation or an exact exchange calculation in DFT that goes also under the name of optimized effective potential or you want to use hybrid functionals either half half B3 lip, Bb0 or brain separated HSC or other that may have been proposed so the problem is that in a plane wave setup the way you the most efficient way if you cannot exploit localization so there are other methods that try to use a Vanier representation but let's say in a plane wave cold you need as I mentioned you need to build these auxiliary functions these auxiliary densities and calculate the corresponding the corresponding interaction so basically for any wave function you need to bring it which we have called written as Fourier component in the basis set so you need to bring this wave function to real space then for each pair of wave functions and therefore which means for each k additional k points as additional band you need to compute a charge density and then you need to solve the Poisson equation and to solve the Poisson equation you need Fourier transform back to real to reciprocal space solve the Poisson equation and then accumulate these results sorry solve the Poisson equation for this auxiliary orbital bring back to reciprocal space multiply by the wave function and accumulate the result so it's a long procedure but the important thing is that you have to do it for every orbital for every k points and each of them each for each of them you have to do a lot of FFT so basically you need to do all this procedure and we will see in a moment that this is actually scales how it scales and how much it costs but there is one technical things that is important is the fact that this integral sorry let me show you this integral that is written there is a double sum over k and k prime but in reciprocal space can be written as sum over g reciprocal space and so that if you write done the this auxiliary densities and then you basically can compute this auxiliary quantity the integral of the square of the of the density the external correlation energy in the Artifox style is an integral in reciprocal space over the brillant zone of this function where this function is a diverging function because you have the one over g square contribution that we have always having in half the energy in this case q is an integral over the brillant zone so q plus g has may have the value of zero so when the when q plus g is equal to zero this quantity or is going approaching zero this quantity is diverging this quantity is diverging it's not really a problem because it is an integrable diverging so you can show that this one can show Gigi Baldereschi back in the 86 shows that this is actually something that you can deal with but not in the trivial way so not as a simple sum over k points so you cannot just sum this quantity over a set of q and hope that the divergence is integrated correct so what they have done is actually remove some term that display the same the same divergence and has the property that can be integrated analytically so basically you remove a term zero times some decay factor and then you integrate this term that you removed analytically and so in this case you can do the integral over q exactly and here in this analytic form you can integrate the divergence and what you obtain by removing and adding back the term is that now the quantity that remain and you want to make a summation of is non divergent anymore and so this term can be sampled accurately if you do a grid of k points so basically the in practice the integral of the of the difference become a sum of g vector and q point and instead the divergence has been integrated and gives you some analytic term that can be computed so this is the idea of the divergence and just to motivate the reason why one do all this is because the energy is much better the energy is much better Hartree-Fock would be terrible so pure Hartree-Fock would be terrible but if you look at a PBA functional we saw that for nitrogen and oxygen the error was reduced from the value that you got for for LDA for instance but was still significant if you go to PBA zero so an hybrid functional the error compared to the experiment is further improved so we are approaching the chemical accuracy limit that is what what we want however the Hartree-Fock scaling is terrible because just to comparison when you do a SCF calculation and we calculate kinetic energy and interaction with an external, local external potential we have to do kinetic energy is diagonal in reciprocal space so just take products of the wave functions times the g square and pw very fast then you go to real space you can transform to real space get the wave function in real space multiply by the potential which is the number of points in real space and then come back in reciprocal space so overall using the dual so called dual space formalism in miltonian with a local potential would be super fast would be basically cost two FFT and a few operations the potentials that we use are not local and so we need to add a number of projectors and the number of projectors is proportional to the number of bands that you have and so this cost more this cost of playing waves times two because you have to you have to do a few operations for each projector but I mean it's a small number times a quantity that now scales quadratically because the number of playing waves and the number of bands both are proportional to the number of electrons that you have so in general is quadratic but with a relatively small coefficient the foc operator as I said for each pair of bands for each k-point you need to multiply make an FFT to go to real space multiply by the build density Fourier transform back in reciprocal space solve the the Poisson equations Fourier transform in real space and combine the results so this require a number of Fourier transforms which is scales with the number of bands and the number of k-points in your grid and the point is that both depends on the number of bands but there is an additional dependence on two points and instead of multiplying by the number of playing waves you have to multiply by how much it costs to perform an FFT and summation over the point in real space therefore typically from one order two order and possibly even more higher order of times slower than a typical normal normal calculation so yes that's it we saw here that hybrids do have better energetics and address at least partially the problem of self-interaction correction but it is hopefully expensive some exercises in the hands on will address this ok ok so an alternative to full hybrid functionals has been proposed independently again in the 90's and it goes under the name of LDA plus U in similarly to the self-interaction correction approach the idea is that you have a functional at that time was LDA but nowadays could be any DGA and you add a correction and this correction is inspired by the Hubbard model community and so this correction is basically a term in the energy that is proportional to some parameter U and times the product of N and one minus N one can look at the original paper to understand how why it was defined like this it has to do to make some kind of simplify Artifoq and then remove some average mean field approximation of the same Artifoq simplify Artifoq description but what is N and what is N of i sigma N of i sigma is the projection of I mean is the occupation of some atomic orbitals so at the A plus U make a model which is applied to some specific orbitals and these specific orbitals are the orbitals that are most localized in for instance in a transition metal would be the D orbitals of the transition metal element and so one has to compute how much is the occupation of these orbitals and then from the occupation of these orbitals one can write this term the energy and with U U is a parameter and this is basically a cost function so basically what these how about you method do apply a penalty the higher is U the bigger is the U parameter the bigger is the penalty 2 partial occupation of the orbitals so if an orbital is partially occupied let's say 1 half then there will be a cost an additional energy of U divided by 8 and so the bigger is U the more it will cost to partially occupy orbitals and so the system will be pushed to try to either fully occupy or fully vacate fully empty the given orbitals and this is the practical effect of LDA plus U what is the justification of these LDA plus U ok different different autor will give different justification but I think what is relevant and is also useful to understand what U parameter is is to focus on what happens when you have a finite system an isolated atom in an isolated atom what you have and you look as a function of the occupation of the levels if you do a calculation with the atomic code and as a function of the occupation of the levels what you obtain is the black curve you get a continuous curve that has a certain curvature and basically has a parabolic a parabolic behavior the true exact functional actually if you look in the literature what the exact function as a function of fractional occupation of levels should be is not a continuous quadratic function but is a piecewise linear function because when you have a fractional fraction of electrons so this doesn't mean that you have a fraction of electron in your system this means that you have a full electrons a fraction of time a certain fraction of the time and so it is the fraction of the value that that fraction of the value that it would have when it has of the energy that it has when it has n and when it has n plus 1 plus 1 electrons and so the exact result should be a piecewise linear function and the fact that the the ionization energy and electron affinity of the systems are different is related to the fact that whenever you cross an integer number of electrons then you are starting occupying different orbitals and also the the slope of this of this curve is suddenly changed which is related to the gap problem that I mentioned before but let's not enter into that but if you look now at the difference between what LDA or GGA would give this parabolic object and what the exact description should be this straight line piecewise linear function then the correction that is needed is exactly something that looks like a penalty function that tends to penalize partial occupations and the you the curvature the spurious curvature that you describe with LDA and GGA that one can show is related to the self-interaction error that we have is associated to the U parameter and so the U parameter is therefore associated to the spurious curvature that you have in the LDA and GGA and that you shouldn't have because of the exact function would be piecewise linear if you do this calculation so if you add this this function and this will be also done in the hands on then basically this penalty break the allow the system to break the system and to the allow the system to break the symmetry break the symmetry and to push some of the orbit I mean to select some of the orbiters one orbiter in particular and to fully occupy it while the other four orbiters of the minority spins are completely empty and the system that you get is insulating as is the case in the experiment and similar things happens to rare hearts like but how do you compute this U parameter this is the in principle this U parameter has been used as an empirical parameter by many people but in principle we saw that there should be there is a motivation this U parameter describe in the atomic case is described the wrong LEA-GGI curvature and this is something that you could do easily in atoms because you have control of how many electrons you put in the system so you can see that the energy is quadratic should not be quadratic you calculate the derivatives and you extract it in solid we don't have a control of how many electrons we put on an orbital because the consham equation will decide how to distribute the electrons but still you would like to see the derivative of the energy the quadratic term the second derivative of the energy associated to the occupation of the levels but you have to take care of the fact that due to hybridization since the atoms are not isolated when you change the position of a level in effect the occupation of the interacting levels but this is not due to interaction this is just the fact to the hybridization of the levels so you have to obtain the effect of change the occupation of a level but without taking into account the fact that you don't have any way when the orbitals hybridize so this is something that we Matteo Coccioni and myself we work on some years ago and we devise a methodology in which we wrote a code to implement this idea in a very efficient way much more efficient than we did back then and it will present this scheme in the afternoon so in the exercise so there is a way but the point is there is a principal way to compute the new parameters and looking at how the occupation of the levels changes when you perturb the the local orbitals the orbitals that are localized and there is a code that now is doing it efficiently in the distribution ok so I think I spent a lot of time talking about the the the this self-interaction energy and the way to address it I have a couple of minutes to talk about van der Waals ok I'm sorry maybe I will try to be fast but I need to tell a few things so van der Waals function also is as I said another big issue especially if you are interested in organic systems or whenever you have physical solution or whenever you have two dimensional materials and interacting with a substrate so I mean this is a subject that recently became more and more important because people are more interested in systems in which van der Waals is relevant and the real point is that it is a correlation effect it is a correlation effect and it is not included in LDA it is not included in GGA it is not even included in MetaGGA that I didn't cover today but also Archifok since it is not a correlation I mean the hybrids that mix a fraction of Archifok may do better for the self-interaction energy but they do not do better for the van der Waals because it's a correlation effect and it is a truly non-local effect in the sense that if you have two fragments the local di Santa Tanya associated to one fluctuation of density in one fragment will induce a response in the in the other fragments with decay law which is inverse power third power of the distance and then when this interaction when this induced dipole and then go back to the original fluctuating dipole gives you the famous R2-6 dependence that we have for instance in Reynolds also or the model potential that describe the van der Waals in principle this is a density functional I mean this is a correlation effect in the density function the problem is that the functional and so should be taken care of of the term Vx in the external correlation function but LDA and GGA contains information over what happens around one point in space and instead as we saw it is actually depends on of at least two points so semi local density functional do not describe van der Waals interaction neither do hybrids or DFT or self interaction correction and if you do for instance graphite and you calculate the geometry of graphite what you obtain is that for the lattice parameter in plane lattice parameter this is pretty much determined by covalence see and it's pretty sharply and pretty accurately this time by any GGA but if you look at the C over A ratio the GGA is very very bad and it's a very large C over A ratio it is very very shallow so it's completely wrong how do you address it well for many years we simply neglected we were interested in more tightly bound systems and so we were not caring too much or let's say the system that we were working on were not needing it and so people live without it but more recently it became an issue and so there are basically two ways practical ways that have been used one is to add some empirical dispersion correction before is also out is has been used and are pretty successful model also has been proposed and both of these two methods are implemented I mean I've been implemented in the quantum espresso and so it will be shown in some session another options I mean this is ok I mean people have done this and there is some accuracy I mean these are accurate methods but somehow it's a bit ad hoc in the sense that this should be a correlation effect should not be an extra term that you have to in the Hamiltonian but should go inside the exchange correlation function and so there are being people in this in this direction it was the I mean the VandervaasDF Langret Lundquist and their collaborators have been working on these things and so the VandervaasDF functionals of several generation of a few generation of this functionals the approach to the to the problem so there are a few functionals that try to address the Vandervaas issue inside the DFT and this is also has been implemented in the also this function has been implemented in quantum espresso and you will be comparing the results of these two ways so just super quickly the green idea is that you write the DFT plus a dispersion correction and this dispersion correction is modeled in a linear jobs like way so with coefficients of interaction between the fragments so which each atom is considered as a different fragment and so the energy is something that is proportional to some empirically computed coefficient divided by r6 there is a dumping function because when two atoms are too close you don't want the energy to explode so f dump goes to zero when the the distance gets too small and there is some empirical parameters there and it has been these parameters have been fitted on a number of systems and they reproduce pretty well a number of situations so that they are can compare reasonably well the energies and geometry of weakly interacting pairs of fragments Tachenko Schaefer do similar things they write the energy again as a dumping term times a coefficient times r6 except that the coefficient are re-computed in terms of some local redistribution of the charge entities so basically this is again empirical, this is again add extra term to the energy so it's not addressing directly the exchange regression functional but the coefficient that is c coefficient that enters in the expression of the energy depends on the local redistribution of the charge density according to some rule of the partitioning of the charge density this is also implemented and they have again very good this coefficients can be computed and reproduced pretty well what one expect from quantum chemical calculation then I finish with this brief description of this truly non local exchange correlation functional that have been proposed by and back in 2006 was I think the first full functional which was named BDVDF and basically is a functional in which the energy is written as the product of two density in two different places times a kernel a kernel coupling but this kernel actually depends in addition is not I mean if this kernel was just a function of R or minus prime that would be like the Coulomb interaction would be just just a different interaction the point is that this kernel this interaction the point R and the point R prime does actually depends on what happens in R and R prime and also on the distance between the two so it's a more complicated object there are different schemes to define what this kernel is they all of them have basically a kind of local density like idea behind them only that local density try to address the local exchange and correlation energy what is doing what they are doing here is trying to address the local polarizability because they exchange and they have Vandervaas interaction and they can be described by the polarizability of a fragment and so actually generally the general idea is to have a local polarizability approximation and then try which is done needs some some assumption and then try to compute it exactly and then from that in q the technical problem is a six dimension functional the comes expensive in order to deal with that okay first of all we see that it is working I mean it described correctly the T shape dimer in internazione o il sanzo che abbiamo visto prima. Quindi è comparabile con i risultati buoni, come potete dire. Questa issue di internazionale sistimazione è stata soltata da cose tecniche che sono implementate nel colpo. Non è necessario che entri nel dettaglio, perché non è qualcosa che puoi vedere. Non puoi fare nulla, perché è già implementato nel colpo. Ma ok, non è una cosa tecnica, ma è soltata. Ok, andiamo con questo. E' diversi funzionali entranti nel modello che usate per questa sensibilità del colpo local. E come potete coprire con questo colpo non local e con questo colpo semi local. Ma tutto questo è definito in funzionali. E ok, questo funziona nel senso che, per esempio, se guardate a Glycin e compare il polniumorfo di Glycin, Glycin è un moleccule più organico con 10 atomi e ha 5 o 6 arrangementi differenti. Si può avere. E l'ordere di questo arrangement è che si tratta solo se inclusiva il funzione di Vandervals e il VdVdf in questo caso o il VdVdf combinano con alcuni specific exchange funzionali. Ok, Alan, ok, questo era solo per mostrare che, se inclusiva queste quantità, il VdVf è un variante di funzionali che è molto bello perché è una forma molto semplice che si tratta per questo senso ma ha un complico perché, al contrario, il trick che è usato per fare il lavoro l'integra efficatamente non funziona immediatamente quindi bisogna fare un speciale tratmo ma ok, ho scoperto questo è stato fatto, è implementato in quattro espresso non il VdVf originale ma il VdVf, questa versione modificata può esplorare il trick di Pérez di Roman Pérez e questo è perché lo chiamo il VdVf e è basicamente equivalente al VdVf quindi il blu e il verde sono errori su un set e vedete che sono piccoli e praticamente the same e meglio che il verde che è il risultato del VdVf2 e queste sono altre molestine e ancora il VdVf il blu e il verde è lontano il C over A e il VdVf4 per esempio può essere computato e è in graffino uno ottenere lontano quindi sto qui, mi sono scoperto ho mai che ho aspettato ho spentato un tempo ricordando cosa è il VdVf2 che è il successo ma ci sono due problemi principali l'errore di interaczione che può essere attraverso da hiberi ma questo anche improvva la termochemistria le energetiche che sono molto spensive o da DFT plus U e ho spentato un po' di tempo attendendoci di questo questo è stato molto veloce ma il problema più grande che non è attraverso da alcuni di l'ultimo è l'interaczione del VdVf che è un effetto non local che ha stato attraverso in molte mani iniziamo da Grim e non locali tutti questi sono implementati in il metodo quantum espresso e alcuni degli esercizi nel il metodo e l'interaczione quindi mi sono piaciuto grazie il tema è molto molto grande ci sono davvero molto interessanti e argomenti anche ricordando in quei giorni ci sono molte domande di qualsiasi funzione di DFT grazie per questa overview perché credo che era davvero utile perfettamente credo che ci sono un po' di domande perché stiamo facendo un paese perché è troppo troppo lungo sicuramente ci vedo o se ci sono tutte le domande in stream nel canale youtube Yuri Timrov io vedo che hanno già risposto di qualsiasi domande e per tutti i participatori che sono nel canale ricordate di scrivere nel canale come stiamo facendo e avremo un risposto grazie stiamo facendo qualsiasi domande e avremo un risposto grazie a tutti alla 10.30 da Yuri Timrov vedo un po' hi Yuri ne vedete la mia screen? no perché ancora abbiamo un screenshot da i CTS assistenti ma magari possiamo la massimo e la removing dai oggetti grazie, potete puoi partire Grazie. Grazie. Siamo andati un po' di minuti per riconvogliere qui. Sì. Siamo andati a testare. Sì. Quindi magari possiamo iniziare ora. Quindi ora abbiamo l'answer in sessione. E' da Yuri Timro. Sì, come sempre, posto tutte le vostre questioni su Slack. Sì, abbiamo molte persone che possono rispondere e potete rispondere per sicuro. E per i tutor, se ti mette necessariamente, si può invitare i participatori a aggiungere le roome per dare più dettagli per aiutare a raccogliere più problemi tecniche. Quindi, grazie a Yuri. Procediamo con D&Sol. Grazie per l'introduzione. Ciao a tutti. Sono da FEDERAL DOLOSAN in Svisiola. È un grande piacere di essere con te oggi o domani. Quindi, il giorno 4 abbiamo un'altra funzione. Abbiamo anche 4 tutor disponibili. Quindi, vi ricordo che potete ascoltare questioni in room break per quelli che sono in Zoom e per quelli che sono on YouTube o potete postare questioni come sempre. Allora, iniziamo. Abbiamo 3 esercizioni. Il primo è per iron oxide. Ho già visto questo nel lecture. Hanno usato DFT plus U metodo. Il secondo esercizio è DFT con funzioni di hybrid. Ci vediamo come questo funziona per silicon. L'ultimo esercizio è usare funzioni di van der Waals e vediamo come questo funziona per graphite. Poi puoi vedere l'esercizio sotto che questi esercizioni possono essere abbastanza piccoli. Questo è il motivo per cui possiamo usare classi allocati per questa scuola. E per quelli che sono on YouTube, potete girare sui laptop o sui workstationi. Depende di quello che hai. Allora, iniziamo con DFT plus U. Un ricordato per il lecture come professore di Jerome Cooley. In DFT plus U, potete corrette l'energia total, il funzionamento di energia, così che abbiamo l'espressione di DFT, l'energia total di DFT e il termo corretto di hybrid U. Questa terma corretta è già la contribuzione di hybrid minus il termo double counting. Questo è l'espressione finale. E questa corretta di hybrid è un'energia così che abbiamo il parametro di hybrid U perché questo è il termo di hybrid U minus il termo exchange J, quindi lo chiamo 1 effettivo di hybrid U. In questa formulazione di DFT plus U abbiamo il simbolo cronico delta m m' e m' sono quanti quanti magneti. Oggi vediamo un obiettivo che è N, che dipende in indice I, sigma, m, m, prime. L'indice I risponde ad atomi, quindi le label ad atomi. Sigma è l'indice spin, che può essere spin up o spin down. Questo obiettivo è chiamato l'occupazione matrix. È definito come l'indice con la funzione psi, che dipende in spin. V, che è l'indice electronega e K è il punto importante. Il prodotto scolar è tra psi e phi che sono le funzioni localizzate. Questo può essere, per esempio, le funzioni atomiche che sono centrate su atomi. E lo stesso prodotto scolar complice conjugato, non complice conjugato, in modo diverso. E anche l'indice m, invece di m prime. E f, sigma, v, K sono le occupazioni electronegiche, per esempio la ferma di Dirac o altro. E se computiamo la traccia di questa occupazione di matrix, noi sumiamo l'indice m e il spin. Abbiamo un obiettivo che dipende solo niente, ma il numero di electronegiche siano sulla traccia che stiamo cercando di correggere. Per esempio, se mettiamo la traccia di un obiettivo alla traccia di D, quindi la traccia di D, o la traccia di F, la traccia di F. Quindi questa traccia di N ci darà come molte tracce siano siano sulla traccia di D. Finalmente, per sapere questa espressione per l'energia, abbiamo modificato l'enegro con la traccia di D. È scritto sotto per ottenere la derivazione funcionale dell'energia, rispetto alla funzione. Quindi abbiamo ottenuto l'enegro con la traccia di D dove il termine blu è il standard di contribuzioni di D. L'energia canetica e il potenziale con la traccia di D. Però ora abbiamo un potenziale extra in red. Questo è l'energia derivativa di E di U. E non ho scritto l'espressione per il potenziale, ma è molto semplice. E in pratica dobbiamo soltare questa modificata con l'enegro con la traccia di D che interessa con la traccia di D o le elezioni di F. Ok, questo era un ricordato. Ora, in regardi i stati localizzati, perché nel slide precedente mi raccomando ci sono alcune funzioni misteriose fai, che sono funzioni localizzate. Quindi, che sono queste? Queste funzioni sono i stati che vogliamo correzzare. Per esempio, in iron oxide che abbiamo visto in la lectura il professore di Geronclay ha detto che le elezioni non sono corrette in LDA o GGA. Quindi, vogliamo correzzare le elezioni. Quindi, queste funzioni fai corrispondono alle orbitose atomiche delle elezioni. E come questo è implementato in Quantum Espresso o come ha detto a la coda che vogliamo correzzare le elezioni di iron oxide. Per ora, in versione 6.7 di Quantum Espresso, questo è molto codito. Quindi, non può essere utilizzato in l'input file. Quindi, questo è molto codito. In la routine chiamato SET HUBBARD N e SET HUBBARD L. Quindi, dove N serve per il numero principale di Quantum N come in Quantum Mechanics e L è il numero orbitoso. Quindi, in caso di iron vogliamo correzzare la cellula di 3D. Quindi, il numero principale di Quantum è 3. E il numero orbitoso è 2. Perché questo corrisponde alle elezioni di D. E questo è codito in queste due routine. In versioni future di Quantum Espresso, chiamiamo di cambiare questo e permettere al user di specificare le molte fold che vogliamo correzzare dall'input. Questo sarà molto meglio e molto più user-friendly. Questo è molto codito e se vuoi modificare questo bisogna modificare queste due routine fortran e compilare il codito. Ma non bisogna fare questo, solo in alcune case rare. Ok, quindi esempio 1 l'iron oxide iniziamo con l'input file e sono già tutti gli esperti per la calculazione Pw per la calculazione Pw quindi non discutirò molte cose solo cose per la calculazione Pw In Yellow questo è quello che dovresti sapere il tipo di calculazione è Cf Abbiamo usato l'experimentale lattice parameter questo è sale dm1 Abbiamo usato smear perché l'iron oxide può essere metallico è un metallo piccolo sullo livello di LDA o GGA Quindi sappiamo che è metallico quindi bisogna usare smear e lo specificiamo in questo caso il smear marzari o il smear cold e anche bisogna specificare il valore del parameter per questo tipo di smear in questo caso è 0.02 W typically this is the default value and we should use this or slightly lower it depends on the system and one has to carefully check this parameter next ok in blue I actually did some mis explanation before so in blue we show the block corresponding to the magnetization this will be described in more detail next week when you will hear about introduction to the magnetism so here we specify and spin equals to 2 which means this is spin polarized collinear calculation and we need to specify two parameters starting magnetization depending on for type 1 and for type 2 we need to specify value in the range from 0 to 1 o from 0 to minus 1 in this case we just take something like 0.5 and minus 0.5 so this value can be changed but in some cases systems are not really sensitive to this but in some cases systems can be sensitive to the starting magnetization in particular from my knowledge in 2D materials the solution to which you converge can depend quite a bit on this starting magnetization so one has to check carefully because there are many in magnetic system there are many local minima ok next we have in yellow also highlighted atomic species we have two types of iron iron 1 and iron 2 and we use PW projector augmented waves to do potentials why do we have two types of iron iron 1 and iron 2 this basically means that we have two sub lattices in this cartoon they are shown in red and blue balls so iron 1 is for example red and iron 2 is blue so we specify these two types of iron because red atoms have spin down and blue ones have spin up so this is why we need to specify two types now coming back to Harvard specific input files so they are actually not highlighted but they are written here first of all we need to add the keyword called LDA plus U equals to true so this is to tell to the code that we want to do DFT plus U calculation next we need to specify parameter LDA plus U kind equals to zero this parameter can take values zero one or two if it's equals to zero it means we use DFT plus U simplified formulation in the do-derive framework if you specify one it's more full formulation according to Liechtenstein which can also be non collinear and you can have also J parameter in that case and more recent addition in quantum espresso is 2 which is DFT plus U plus V where V is the inter-side correction this was not mentioned in the lecture I have some slide about this I will discuss in more detail later but for the time being just keep in mind that we do simplified DFT plus U calculation so we set zero then U projection type this is very important point because in DFT plus U as I said before we have localized functions phi so this functions can be different I mean and depending on which functions you choose the value of U also can be adapted I will explain it can be computed but basically these functions phi are controlled by U projection type this is the type of the Hubbard manifold or the Hubbard projections in quantum espresso there are a few options one is atomic corresponding to atomic like functions so this corresponds to atomic orbital thread from the pseudo potentials in pseudo potentials you have information about atomic orbitals of an atom of iron so basically we take the orbitals of iron in this case another option is ortho atomic so as you can guess ortho means orthogonalized so what we do in that case we take atomic orbitals of iron and we orthogonalize them between different sites also including oxygens and it's even possible to use vanier functions this is less developed for the time being and this is working in progress so typically we have atomic and ortho atomic options and even though here it's written atomic but for production purposes it is recommended to use ortho atomic because it's from our experience it gives more accurate results and finally we need to specify the value of the HABARDU parameter effective HABARDU parameter in this case this is a trick we essentially put it equals to zero because we want to do PB sol calculation in this case so without HABARDU correction but we put it some very small number just to activate the DFT plus U machinery because when you use DFT plus U the code prints extra information in the output file so we want to see this extra information when we don't have U correction so that's why we just activate all flags but essentially we are not applying any U correction so we are just trying to cheat the code and we specify two values of HABARDU parameters one corresponding to iron one and second corresponding to iron two in this case they are equal but if you say for example have a system which contains different elements from the periodic table for example iron, nickel cobalt and other elements of course each element needs different HABARDU correction so we need to specify them here but in this specific case of iron oxide we have two types of iron but both are iron so if these atoms are crystallographically equivalent the HABARDU should be the same but if for some reason there is some complex distortion low symmetry and these two types of iron become crystallographically inequivalent then we need to apply most likely different HABARDU parameter and the question is how to find those parameters it's a big problem it was problem until sometime ago but now we have first principles method which allows us to compute those HABARD parameters for any element any material at hand so this is no longer an issue I will talk about this in more detail in a moment ok another interesting point is that I would like to let you know about the so called quantum express input generator because I've just discussed the input for iron oxide and during the first three days you learned a lot about how to set up the input how to choose key points cut-offs that you need to perform convergence tests that was the routine maybe like 10-20 years ago and earlier but nowadays we want to improve and make it more automatic so for this reason it was created quantum express input generator the link to this tool is below so the goal of this tool is that it provides to you a first initial guess of the input you can do you can upload a crystal structure this can be the input file that you just prepared but let's say you don't have a lot of experience in setting up all those parameters so you do something that you learned but then you want to verify if that is a good initial setup or not you just upload it here and this tool will provide you on the input the suggested parameters for key points, for cut-offs, for potentials or if you don't have input file you can upload CIF file or any other formats you can choose then yeah this is exactly the parser the file format that you upload then the pseudo-potentials I will mention the next slide you already heard a lot about pseudo-potentials that there are norm conserving ultra-soft project argumented wave but again which one to choose if you are a beginner you can be lost and don't know what to do but nowadays we have standardized library it's called SSP so basically no longer problem you just go there and pick up pseudo-potential for the element you need next you need to know preferably the type of your system is it metallic or insulating is it magnetic or magnetic so of course it's not easy to tell but maybe you know it from the literature or some other source but you need to provide this information as an initial guess and then you specify the density of k-point grid so it's not the grid directly like 2x2x2 or 4x4x4 but it's rather the the distance in k-points in the reciprocal space in the brilliant zone so in this case for example we want to tell that please tell me which k-point grid I need to use in order to have the distance between nearest k-points of 0.21 over angstrom and this will give you the accuracy of 0.2EV on energies and this corresponds to the fine grid but then there are other two options which give more denser k-meshes so this can be controlled here and of course depending on your structure you have different size of the unit cell so for the small cell you need more k-points in the brilliant zone because the brilliant zone is larger instead if your unit cell is large in real space so the brilliant zone will be smaller you need less k-points so with this QE input generator it can provide you the good initial guess for the input parameters but of course you should verify that you just need to still do some checks convergence checks on k-points decrease it a little bit increase and see that it's really the good suggestion by the tool and speaking of SSP library the link is below but again this is very well tested studio potentials so what was done people from from EPFL they took different studio potential libraries existing and they compared all of them so many properties were computed like phonons total energies and many many other properties for each studio potential all these properties were compared with all electron calculations that do not use studio potentials at all so that is a reference and by doing this comparison we pick up the best performing and the most accurate studio potentials so if you need some studio potential for example iron we see magenta color and oxygen awesome magenta so actually I don't have it on my slide but this corresponds to some specific studio potential library which is PS library generated by Andrea Del Corse using the atomic code that is included in quantum express but if you take some other element for example nickel and again oxygen so oxygen is from one library and nickel from another library all right let's move on coming back to DFT plus U before we discuss the input for the SCF calculation cell consistent field calculation and this is the input for the NSCF non cell consistent calculation we consider this input because we want to compute projected density of states so as you learned already first we need to do SCF calculation and then NSCF calculation and here many things are the same exactly the same as an SCF just in yellow highlighted we need to add a number of bands in this case 35 so we didn't specify this parameter in the SCF input because in the SCF the code will compute just occupied states for insulators and for metals it computes occupied states plus 20% of empty states so since iron oxide is a metal we add smearing so we have not only occupied states but a little bit of empty states and partial occupied states but here we since I already checked before that I want to include a bit more that was already computed so we add a bit more empty states so this is why we have 35 but for each system this should be checked must be checked and also here we use K point grid because in the SCF we use 3x3x3 and here we use 6x6x6 because NSCF is much cheaper so we can use denser K point grids but notice that these parameters are not converged this is just for the demonstration purposes so of course you need to carefully check the convergence of the calculation respect to K points and even cutoffs we are also very low ok so SCF and SCF and the last point is calculation of the PDOS this is the input file we specify type of smearing gaussian in this case the broadening for this gaussian function which is in Riedberg units not electron volts and the energy range in which we want to plot the project density of states and delta is the step and once this is done we will use the GNU plot script it's called plotPEDOS dot GP to plot the spectrum ok so now let's do perform calculations ok so now we are back to the virtual machine we are in day 4 example 1 dot DFT plus U so any questions so far if anything you can either unmute yourself or you can raise your hand also yes ok if not let's move on so as I mentioned before we will use HPC facilities for this exercise because it's quite heavy in this repository in day 1 in example 1 you click on readme.md and these are the instructions so you can proceed on your own and I will try to proceed as well so we have 3 parts one is standard DFT calculation second DFT plus U and third calculation of Harvard U so let's do the first step ok ok let's use HPC now let me remind in case you missed it or you forgot if you click on Mozilla browser here there is a button called how slash 2 ok you click on it then you choose remote access usage and you scroll down and these are the commands remote mpi run to submit the calculation with the default number of course which is 20 and we will use remote sq to monitor the execution of the calculation so please do not open the output file when the job is running submit the calculation and monitor the execution with remote sq and then we copy file to our virtual machine so let's do that so we type remote mpi run then pw.x minus input and the name of the input file this is pw. pw.x ok so the job was successfully submitted now let's see what's happening type remote underscore sq press enter ok so the job is running already 10 seconds it should be very fast calculation this is standard dft calculation still running ok while waiting there is a question on youtube how is the double counting correction implemented so there are different types of double counting and in quantum express so one of them is implemented it's called fully localized limit next missing ssp from different libraries a mixing sorry mixing from different libraries is it ok and different norm conserving of pw yes this is ok you can mix potential from different libraries and if even different types of potentials you can use norm conserving and with ultra soft what you cannot do is to mix functionals for example for iron you use LDA and for oxygen you use PBE so that is not as possible but as long as you use one functional you can mix you those ok this calculation should be finished ok now let's copy the file from the cluster here on the local machine so this is done by rsync underscore from underscore hpc and you do quotes star dot out enter ok copied files and we have this we can open the file output file pw dot if iron oxide dot scf dot out ok you are already familiar with the output let me just highlight what is new with respect to in dft plus you with respect to dft so I just remind that this is dft calculation but we activated the dft plus you machinery just to print some extra output but the you correction is actually zero so what is in particularly interesting is that at the end of the consistency the dft plus you machinery prints some information about the occupation matrix so let me recall what is this occupation matrix is this object n i sigma mn prime so this is a matrix ok respect to mn prime in the case of iron we are considering 3d electrons so the orbital quantum number is 2 so the size of the matrix with respect to mn prime is 5 5 ok so this is 5 by 5 matrix and we can take compute the trace of this matrix so take the sum of the diagonal elements and this will give us number of electrons in this 3d shell what we can do we can also diagonalize this 5 by 5 matrix and obtain its eigen values which can also tell us what's happening in the system so let's see so in the output we see the trace of the occupation matrix this matrix is called ns n a is the atomic index so we take the trace of this 5 by 5 matrix and we have for the spin up for the spin down contribution in the total which is spin ups plus spin down for the spin up we see 4.98 so this is the number of electrons in the spin up channel of the 3d electrons of iron and this is number of electrons which is 1.8 in the spin down channel for the d electrons of iron and if you take the sum you get 6.78 electrons overall in the 3d shell so we know that roughly almost 7 electrons on the 3d shell not exactly 7 but around 7 but this is not very accurate estimate because we use atomic projectors so it's better to use ortho atomic projectors which will give more physically meaningful value for the number of electrons and you can compare determine the oxidation state of iron oxide I mean if since oxygen is 2 minus so iron is 2 plus so you can count electrons and see what you should expect so for atom 1 this is information about number of electrons, spin ups, spin down channel the total and these are the eigen values this is what they said before when you diagonalize 5 by 5 matrix you obtain 5 eigen values so on the diagonal you have 5 eigen values so we see that all of them are almost 1 so we are fully occupied and also you have eigen vectors and so this for the spin 1 which is spin up and for spin 2 which is spin down you see that in the spin down channel instead of 1 we have some fractional occupations so we see that those are partially occupied so when it's 1 it's fully occupied when it's 0 it's fully empty so here we have partial occupancy in the spin down channel and and then for the second atom because we have second iron we have the same information number of electrons it's actually reversed so in the spin up channel we have 1.8 in the spin down 8 but the sum is the same 6.7 as before and now again here we have some reversed picture so in this spin 1 we have fractional occupancies in the spin down fully occupied ok now let's continue now we do the nscf calculation same way so we do wx read the input now add nscf and we execute ok let's check if it runs properly yes it works ok still running ok one question by colleague Romain Nikita about we can see LDA plus you parameters what to the potentials are PbSol is it ok yes we are using PbSol functional because this is three dimensional crystal so for bulk systems we should use PbSol gga functional PbE so because PbSol I mean Sol stands for solids so PbSol was optimized for solids while PbE that was optimized for molecules so PbE should be well historically PbE was created before PbSol that's why people applied also PbE in solids but PbSol is now developed so we should use PbSol for solids Yuri we have also one raised and you want to answer now or? ok so let me pancake here i now and and we i may at with poor Good morning and sorry for the interruption. My question is why did you write remotes in the terminal in the beginning? I didn't get the question Ivan, did you understand? Why did you something in the terminal at the beginning? I didn't get just the middle word. The command, the command remote. Okay, thank you now it's clear. Okay, remote underscore MPI run. So this is the special command created for this school because we are using the cluster. So when we execute this command, what happens underneath is that the calculation is submitted in the cluster. So not on your computer, so via SSH. The message is passed to the cluster one or three clusters that we use. And there, there are local scripts on the cluster which are submitted. And the calculation is run there. So, I mean, if you're connected directly to the cluster, you're using MPI run or other analogs on the HPCs. But here we have a special command for this school, which does more than just that. It simplifies our life a lot. So just with this simple command you submit calculation. That's it. So this was explained yesterday. All right. Thank you, Samira. And thank you, Yuri, of course. We have questions from streaming. Do you want me to read them or you read yourself from the. No, please, please, let's answer one more and then maybe we answer others later. Okay. So the first one from streaming is mixing SSSP pseudo potential from different libraries. Is it okay? Or if different non conserving PAW on Ultrasoft? Yeah, yeah. This actually is already answered. It's possible to mix it. Yeah. Okay. The next one is how is the double counting correction implemented? Also answered. Fully localized limit. Maybe the last ones. How we can read the eigenvectors for each eigenvalue by row or column. And the last one, how should they give FE1, FE2 for inequivalent atoms? And the last one, what difficulties can NBE calculation have in magnetic systems applied the FT plus you. Okay. Regarding the eigenvectors. I mean, I never really use those data information about eigenvectors. And that's why I cannot answer. Is it by row or by column? But this can be easily verified by checking the code. So I will check this and I will answer on Slack. Okay. So this I will answer later. Yeah. Maybe Yuri, because the guy who posted this question is on YouTube streaming. Oh, yeah. So I will. So he cannot get the answer by Slack. So maybe we can answer later on during the end son, or maybe even in tomorrow or in voice because. Yeah. Yeah, I know the name of the routine. I mean, if you're familiar with Fortuna, you can check. Right on YouTube. P W slash S R C slash write underscore NS dot F 90. To check. So P W S S R C. Right underscore NS. Right underscore NS. Yeah. There you can see how it's printed dot F 90, right? Yeah. Okay. So I, I, I just written in the YouTube stream chat. And so you can check that. Yeah. Okay. Next question. Should I give iron one iron two for inequivalent atoms? Yes. So in this case, we use iron one iron two because of spin. Because we want to have one spin up spin down, but these are still equivalent. It happens that we have just to iron in the, in the hour simulation cell. But if you take some other polymorph of iron oxide, which has bigger cell, so you have more than two iron atoms. It can happen that not only you need to define two types of iron because of spin but it didn't happen that even there are some crystal graphic inequivalents tra i rinforzi. Quindi potete definire magari più tipi di rinforzi. Per giocare a questa strada di diversi tipi di rinforzi, potete controllare il spin e anche l'equivalenzia tra i atomi. Quindi potete avere i rinforzi 1, i rinforzi 2, i rinforzi 3, i rinforzi 4. Il primo è di corrispondere al spin, il spin down, l'equivalenzia tra i rinforzi e i rinforzi. Potete anche spin, il spin down, ma c'è un altro formato di atomi in equivalenzia. L'ultima domanda, quale difficilità può... Qual è l'NBE? Che può essere? Alcalzatione? May be they mean the NAB. Ah, NAB. Ho detto, non ho capito. Soprattutto non ho capito la domanda. So, per esempio, William Sennin, potete repostare la questione su Youtube. Ok, c'è anche un'altra questione. Il pbe sol è meglio per i materiali 2D di pbe. Quindi nel nostro gruppo usiamo PBE per i materiali 2D, perché PBE è per le molte, PBE per le soli, 2D è in-between. Quindi c'è una domanda, è PBE o PBE soli? Well, from my knowledge people use typically PBE, but you should compare PBE and PBE sol, just to be safe. Okay, and Stefano answered eigenvectors are by column from lowest to highest, the same ordered eigenvalues. Okay, thank you. Okay, just to move on because I'm super late. So we did SCF and NECF. Now the last step is a calculation of the PDOS. So we type remote MPI run. It's called projwfc.x minus input. And we type the name of the input file, which is projwfc, ion oxide input, emit, remote, SQ, already finished. Yes, takes a few seconds. Now we need to copy files. Let's do this way rsync from HPC dot. We should copy everything because I don't remember the names of files and there are many actually. So we should copy all files from the cluster to my local machine. So we see that the PDOS calculation produced many files, each of them corresponding to specific type of orbitals. We have the electrons of iron, we have P electrons of iron, D electrons of iron. So what we want, we want to plot a figure which we saw in the lecture for iron oxide. So we want to plot D states for iron one, D states for iron two and P states of oxygen and see how this looks. So for this we have a script. Where is it? It's called plot underscore PDOS GP. So we just type glue plot and the name of the script. So let's do that. The name of the script. Okay, so this produced file called iron oxide underscore PDOS. Let's visualize this. So we use this common atrial. So you should obtain this. I guess many of you already obtained it half an hour ago because you're faster. But this is in fact what we saw in the lecture. So the vertical axis is the project density of states horizontal axis is the energy minus Fermi energy. We see in red iron 3D majority spin. So we have in fact they are fully occupied because we saw that Eigen values are one, essentially one. So we see they are all occupied. And for the other iron minority spin in green, we see that they are around the Fermi level. Because we saw that Eigen values were not one, but they were some fractional numbers between 0 and 1. And for oxygen we have blue PDOS, which is also mainly occupied. So this is a problem of standard DFT with LDA or PBE or PB. So now we want to apply the plus you correction. So to do that. So you follow not this one. So in the extractions of this exercise, you do step two, you add the value of the Abardu parameter, which is 4.6. And all the rest is the same, but you just repeat steps that you saw before. SCF and SCF, project density of states. Okay, so and I will not do that in the interest of time. So let's see what is the final result. So on the left is our DFT result. On the right is the DFT plus result. So when we apply the Abardu correction. So what happens is that we apply the, we remove self interactions for the D electrons. So they were all delocalized in PB. So but by putting plus you correction, we are localizing them more as they should be. And when we do the self consistent cycle, what happens is that the calculation converges to the insulating ground state. So we have a bang gap now. And in fact, this agrees with the experimental observation that iron oxide is a insulator. So this is good news. So plus your correction works. Now the question, very interesting question, and it's the main problem of the DFT plus your method. Why and how to choose the value of the Abardu parameter? Why did we use 4.6 electron volts? So 4.6 EV, we use just an example for demonstration purposes. But what we need to do in practice, we need to find the value of the new parameter somehow. So what people, many people typically do, they use empirical value by trying to reproduce, for example, the value of the bang gap. For example, you know that in iron oxide, there should be bang gap. I don't remember the experimental value. Let's say, I don't want to say some random number, but some value and we want to reproduce that value as close as possible, but we can try different values of Abardu parameter. We can pick up that one, which gives us the best bang gap. So this is what many people do. However, this is not very satisfactory approach. Because, first of all, DFT plus U is not really the method for computing very accurately the bang gap. DFT plus U is the method to correct the energetics of the system, to remove self-interaction and to have correct energetics of the system. Bang gap is a spectral property because you need to go to, you need to have conduction states, so on and so forth. So trying to reproduce experimental bang gap is already not really the best way to go. But the better choice is to compute Abardu using the first principle methods. There are different methods. Professor Dejaron clearly explained that there is linear response method, but there are also others like constrained random phase approximation and some other methods. In this slide, I would like to tell a few words more about linear response. And also, as Professor Dejaronko mentioned recently, we recast this theory or formulation in a more computationally efficient way, using density function perturbation theory that was initially developed for phonons, in which it will be discussed tomorrow by Professor Baroni. So the FPT was created in around in the late 80s. The first paper was in 1987 for phonons. Then in 1991, very extended physical review, B paper by Professor Giannazzi, Professor Baroni, Stefano Dejaronkoly and Andrea Del Corso, or maybe I'm wrong. And then there is a review paper in 2001, if I'm not wrong. But you will learn tomorrow about the FPT. This is a very powerful theory. And we use this intelligent idea also for this problem of computing Abardu parameter. It's quite theoretically involved, so I will not discuss here. I must say that you is defined in this way where chi and cannot are so-called response matrices, which can be computed by computing, by obtaining the response of our occupation matrix to some perturbation with the strengths of the perturbation lambda. In this FPT framework, this object is defined as the sum over q points over this complex object delta n of q. Unfortunately, I don't have time to explain in detail. But the message to pass is that there is some q point sampling of the brilliant zone and each q points q vector, if you want, corresponds to the wave vector of the perturbation, so-called monochromatic perturbation. So we have several monochromatic perturbations and for each of those perturbation we compute the linear response of the system by summing up all responses. We obtain the total response from which we obtain the value of Abardu parameter. So physically what we do in our system, we apply some perturbation to the electronic occupations. For example, in 3D shell of Pyron. So we perturb slightly the electronic occupations in the 3D shell and see how the system responds to this. And from this response, we can obtain the value of the Abardu parameter and how to compute it with quantum express. So there is a recent code called HP standing for Abard parameters that implements this idea using density function perturbation theory. The input is very simple. You just provide prefix and output directory if needed and the q point grid which I mentioned. In this case, it's simply one by one by one to make it fast. But in practice, you need to converge the value of Abardu parameter with respect to this Q mesh. So Q mesh and also K point mesh. So you value must be converged. So this is why we use 4.6. Because if we compute Abardu for iron oxide using Q mesh 111, K mesh 333, we obtain 4.6. But of course this is not converged. So this cannot be used for production purposes, but this is just for demonstration purposes. Sorry, there is a question which I think is pertinent regarding the calculation of Abardu parameter. That is, if there is an easy way to determine a range for it, there is a range for the U value. So yeah, if there is an upper bound for... Is there any upper limit of U? No, no, there is any constraint. In fact, this method is not suitable for closed shell systems. So when the D electrons, this shell is fully occupied. In that case, this method gives an physically large Abardu parameters. It can be 20V, 50V, an physically large. So this method is not suited for closed shell systems, only for open shell systems where the electrons are partially occupied, like in iron oxide. And also since we're talking about questions, there was a question on a Slack or in Zoom, I don't remember. So because some participants says that the U parameter is applying for to 6S states in some system. And that was done using Empirical U to reproduce the experimental bandgap. So question was, can we apply U to S states or even to P states? So typically, U correction was introduced, it was explained going from the Hubbard model when there is a problem of strong correlation. So historically LDA plus U or DFT plus U was used for correcting D electrons and F electrons. But then later people realized that actually when you introduce the Hubbard correction in the framework of DFT, it actually cures self interactions and not correlations. Correlation problem is very complicated problem, strong correlations. You need multi determinant theory and there are symmetry arguments, it's very complicated problem. So DFT plus U is not the theory to cure strong correlations, it's really curing self interactions. And self interactions are present not only for D and F electrons, they're also present for P electrons and S electrons for all types of electrons. They are larger for D and F electrons and they are smaller for S and P. So in principle, why not applying U to P states or even S states? And in fact in the literature, some groups apply U correction to the P states of oxygen for example. There is no general consensus on this point, but this is what people do. But putting U on P states is, well this has to be discussed. For the S electrons it's even more questionable. But then people do what they want. You apply to S6S, you try to reproduce experimental bang up. Yeah, I mean, I don't know what to say. Yeah, so that's basically the comment. To move on, we can compute a how about your parameter. Let's do that quickly. So this is explained in point three, you can skip the step of cleaning temporary directory. All we need to do is to run the SCF calculation as before and running the HP code, HP.x, read this input. So since we are running on the cluster, we don't need to add this name of the output file. And instead of this symbol, we need to put minus INP or minus IN. So we do remote MPI run, W.x minus INPUT. And the name of the input file PW.ionic site SCF input, remote SQ, to check. The HP calculation will be much longer than the PW calculation. The PW is the ground state calculation while HP is the linear response calculation. It can be 20 times longer. It depends on the system, depends on the density of the Q point grid and other things. Okay, so SCF is finished. Now let's submit the HP calculation. Sorry, remote MPI run, HP.x. And the name of the input, it's HP.ionic site.in. Sorry Yuri, may I say, just interrupt you one minute. Okay, for the people that is running on the CISA cluster, please, before launching starting jobs, remote jobs, please run the update CISA reservation name command as it is in the instructions in the auto. Because if you don't do that, you run the jobs in the queue of all the users, not in your reservation for the school. And so first of all, you remain in queue without the job executing. Second thing, you occupy nodes that are not reserved for the school. So please, just read the instructions in the auto remote and follow the instruction for CISA. This must be done every morning before starting the end zone. Sorry for the interruption and the inconvenience. Thank you. Okay, the calculation is finished. I repeat here. So we do our sync from HPC. Okay, and we can check the output file. Well, HP calculation will produce file called ironoxide.aborparameters.dots. So we can open it. And in this file, you see information about theabordew parameter, which was computed. And we see that, in fact, for iron one, the U value is 4.655. And for iron two, it's exactly the same value. So we obtain exactly the same value for two irons and because they are crystallographically equivalent, they just differ by spin, but this doesn't mean that you should be different because of spin. There can be spin resolved, theabordew, but this is also work in progress. For the time being, we have global U for both spins. Okay, now a couple of questions before we go to the second example about hybrid functionals. We have an issue with the raise done and we answer it. Thank you. So on this issue of calculating U for orbitals other than 3D and 4F, I was thinking about these projectors that you have or these inner products between the bands and the atomic wave functions. And my question is, is there a problem if the wave functions that we use are actually pseudo wave functions which don't have nodes? Okay, no, it's not a problem. So we use, in fact, pseudo wave functions. So for the conchia wave functions, we have pseudo wave functions and we project them on the 3D atomic orbitals of iron. So no, there is no any problem there. So we don't need the information about core electrons for the conchia wave functions because in conchia wave functions, we have just valence electrons which make the chemistry. So we project those on 3D electrons. Excuse me, if you want to project the conchia orbitals over atomic orbitals which have nodes, but you are actually approximating those with pseudo atomic orbitals that don't have nodes. Okay, yeah. That was the question, if that could become a problem or not. Ah, okay, sorry, yes. We were referring to atomic orbitals. No, no, it's not a problem. No issue. Okay, thank you. You're welcome. Then other comments? Yes, we have another. Sir, my question is, did we need to modify the input before doing the second calculation, the calculation of u parameter? Because when I check now the input that we have in the day for, there isn't the value of u and the input. Okay. Yeah, so in this example, there is a note that before running the SCF calculation, you need to initialize back how about u to 10 to the power minus 8. Yeah. To start from the PbSol ground state. I think what you did, you forgot to do this. So I guess in your SCF input, you still have 4.6, correct? I think so. I will check. And so what happened is that in your test, you use 4.6. So this means actually your ground state calculation is dft plus u. In this ground state, you compute the hybrid u and it's different. And this is expected. So this is another important point that actually I was not planning to discuss, but this is important. So we need to do self-consistency when we compute the hybrid u. So what does this mean? When we compute the hybrid u from the dft ground state, so with u equals to zero. What happens is that our ground state is completely wrong. It's metallic. And so from this completely wrong ground state, we obtain u 4.6. If we repeat this step again, we put in our SCF the value we just computed 4.6. We redo the ground state calculation with dft plus u. And we open the gap. So now our ground state is way better than before. We still have insulating ground state which is much closer to the truth. And starting from this improved ground state, we compute again the hybrid u and you obtain something maybe 5 or 6 EV. I don't know what did you obtain. So what we do next? We use this updated value of u and plug it in again in the SCF calculation. So we redo again the dft plus u ground state calculation and we compute u again. We do this so-called self-consistency loop until the hybrid u does not change. So this way we obtain the self-consistent value of the hybrid u parameter. If we do just one shot, it's called one shot calculation. You just do dft and on top you compute u. That is not the value you should use for production purposes. You need to do it self-consistently. Moreover, there is an external loop of self-consistency for the geometry. So to go even further, what we do, we not only compute u self-consistently, we also relax, optimize our structure in the external loop. So we do this double loop for the structure and calculation of u until the value of the hybrid u does not change and until geometry does not change. So this is very important because your hybrid u is consistent with the geometry. And that value of u that you obtain and that geometry that you obtain should be used for the production purposes. There is a paper we published in this year in January in Physical Review B that discusses this point. So you see dft plus u is quite tricky when you start working on it. It's not just picking up some random u or empirical u that in calculation publish something. There is a lot of details that you should be aware. Thank you, sir, and if you can let us the title of the paper. Yeah, I will write it on Slack and on YouTube. It's Physical Review B 2021 where I'm the first author so you can Google it. Ok, next. Thank you, Samira. We have four, which might be reduced to condensate to three questions from streaming and another raised hand. Would you like to proceed to answer? Yeah, the point is that we have a little time. We have 15 minutes, correct? And then we have another half an hour. Can we use that half an hour? Well, we have, yes, no, we have until 12.30. Ok, so yeah, we have 55 minutes. Ok, ok, because dft plus u was initially for half an hour, but now it's one hour and 15 minutes. Ok, so let's answer a few questions. Just do as you want, if you want to. This is my favorite topic so I can continue forever. Ok, let me explain the last thing and then again questions and we move on. Ok, because I mentioned about dft plus u plus v. This is extended formulation where we have not only the Hubbard U correction, which is on site. So we are correcting electronic interactions sitting on iron atom. But also we are adding another contribution, actually subtracting another contribution, which depends on the v parameter, which depends on two indices, i and j. And this is the so-called inter-site interaction. So in the case of iron oxide, this extra contribution will correct interaction between 3D electrons on iron and 2P electrons on oxygen. So this is why it's called inter-site or interatomic interactions. And in this cartoon we see that in dft 3D electrons are overly localized. In dft plus u we are putting plus u correction, so it's a penalty where localizing electrons so they become less spread and sharper. In dft plus u plus v in fact we apply plus u to make it sharper, so to localize 3D electrons. But we also take into account the effect of ligands. We know that there are oxygens, oh this is not the right structure, but we know that there are oxygens around iron. So and 2P electrons of oxygen interact with 3D electrons of iron, so we take into account that interactions and we correct that. So this extended functional gives more flexibility and this functional allows us to correct not only the localization problem of 3D electrons, but also correct the interactions with ligand states. So this general formulation is very powerful and from our experience dft plus u compared to hybrids because some people asked if dft plus u is better than hybrid functionals or GW. So dft plus u is worse than hybrids, sometimes they are quite compare well, but sometimes hybrids are superior to dft plus u. And if you use this extended formulation dft plus u plus v, this extended formulation is gives much closer results to hybrid functionals. But it is important to remember that we should not take hybrid functions as a reference. Because hybrid functions are also approximate, you have a fraction of exact exchange which is 25% in PB0 and HSE06, but hybrid functions also fail in solids because this fraction of exact exchange is not proper, not correct. So there are more advanced hybrid functionals, they are called the electric dependent developed in the group of Giulia Galli in Chicago and other words by chronic and other groups. So they choose the fraction of exact exchange ab initio, so they compute it. Actually I will discuss about this soon, but to summarize standard hybrids like PB0, HSE06, they are also not perfect for solids and taking them as a reference is not always a good idea, but in some test cases that we considered from our experience dft plus u plus v is in good agreement with hybrids or sometimes even superior. Ok. And to finish, not only u but also v can be computed using the HP code in quantum expression using the latest version. Ok, before we move on, let's answer a few questions on this first part. Yeah, we have two questions from streaming, which are related, we can plot the density of states of subbands, maybe subbands of d orbitals, I mean dxy, yz, zy, and another related, how can we plot dEg and dT2g states. Ok, so for the first one. Yes, we can plot different contributions, you just need to open this file. Let me try to do that. I'm not sure this works. Ok, so one, then five. Ok, if you open this file produced in the p dot calculation, it doesn't want to fit. So we do view zoom out. Ok, in the p dot file, you see lots of numbers. There are many columns. The first column is the energy. Second and third are local density of states for the up, spin up, and for the spin down. So this is what we plot actually in this example, plot for the spin up and spin down. But if you are interested in really contributions to this local density of states of this dxy, dyz, dzy, etc. So these are, I highlighted them. For the spin up channel, we have five, of course, because this is the D shell. And for the spin down, we have five as well. Also at the end of the p dot output file, where is it? Should be proj output. At the end of this file, there is useful information, so called loading charges. For atom one, you see lots of interesting information and what you're asking is essentially there. So for the spin up, the electrons, you see dz squared, dz, dyz, etc. So you see those numbers. And you can identify and plot them using the file I showed before. And regarding t to g e g, this is not explicitly written anywhere. And actually I have also troubles to identify what is what. So you need to plot these contributions and analyze a little bit the results to identify where is t to g, where is e g. But you don't have something explicitly written t to g or e g in the output. But I think this information is enough to obtain what you need. And then what else? Another one from streaming hp.pix code. Calculate abbar parameters for the atoms that have abbar parameter assigned in the SCF calculation, right? What if I wrongly assign abbar parameter for all in the SCF file? So the last question. What if I wrongly assign abbar parameter for oxygen? What does it mean? I mean, if you wrongly assign, well, you will obtain wrong results. I mean, maybe I didn't get the question right. I mean, yeah, maybe if you want to write, if you want to write more, maybe write reformulate it on YouTube so Yuri can answer more precisely. For example, for iron we put 4.6 and 4.6. And for oxygen we put abbar d3. So we put 3 because oxygen is type number 3 in the atomic species list. And for oxygen we put the value we compute where we want. Ok, so next. Is it better to use dft plus u plus v to evaluate the band center? Not sure. I understand the question and don't understand what does it mean evaluate the band center. So again here maybe. Yeah, then the question can they have a new parameter, how about you value be different according to the applied potential. Yes. The hubbub you parameter is depends on the potentials depends on the hubbub manifold, if it's atomic or orth atomic depends on the oxidation state. So yeah, when you report value of the hubbub you in your paper. Remember to report which pseudo potentials did you use, which functional, which hubbub manifold, etc, etc, so that people can reproduce your results. Otherwise, if you just report u plus 5 e v, and that's it. That's not enough. Not satisfactory. Ok. How dft plus u is better from tb mbj. Fortunately, I don't know. I don't have experience with that functional. Is there a way to determine the g factor of transition metals of radicals using quee that is not related to this hands on. Okay, so there is a question about hybrids but let's discuss it in the next part. So to raise them if you want. So we have a half an hour and we have still hybrid functionals and wonder was so let's continue. Actually they are shorter so they should be 15 minutes each. Okay, so let's go to hybrid functionals. For example two and apologies for those questions that I couldn't answer now. So presentation hybrid functions. Okay, this is the slide to remind what you saw in the lecture. Because in hybrid functional there is this problem of divergence for the q plus G when it is zero in the summation. So there is a Q mesh in hybrid functionals. So what the meaning of this Q mesh in the expression for the exact exchange energy. There is a double sum over bands and key points. And that is very expensive having this double sum over key points. To reduce the computational cost. One of the key points sums is replaced by some over Q points which is has less Q points than in the original key points summation. For example, if you have a double sum of our key points and each of these sums runs from one to 100. 100 times 100 terms. That is a lot. So what we do, we say the following. The first key point sum we keep from one to 100. And the second key point summation we replace by Q points from one to 10, for example. So basically we are reducing the second sum over key points using this trick of Q points. Okay, so this is the meaning of Q points and there is the singularity problem. There are different ways how to treat this singularity problem. And the one was mentioned in the lecture is the so-called GG Baldereschi. I think it was, if I'm not wrong, it was created in Lausanne in Switzerland and published in this work. So we will use this correction method to treat divergence in this example. The input file for silicon with hybrid functionals is shown on this slide. In yellow I highlight input parameters related to hybrid functional. What we need to specify is the keyword called input underscore DFT equals to the name of the hybrid functional that we want to use. In this example we use PB0 functional, hybrid functional, so we write PB0. If you want to use B3 lip, you write here B3 lip and HSE respectively. It is important to notice that when you select studio potentials, you need to use those which are generated with the functional, which is as close as possible to your hybrid functional. So this means that if you use PB0, you need to use studio potential generated using the PB functional. You should not use LDA functional in this case because you're using PB0. Unfortunately, nowadays we don't have studio potentials generated with hybrid functional. So this is why people use studio potential generated with PB sol functional or LDA. Next we specify the Q point grid, these parameters NQX1, NQX2, NQX3. So we use just one Q point, so we write one by one by one. This is really the simplest approximation we can do, which just say there is one Q point, just gamma. Then we need to set up X gamma extrapolation to true. This means that we are taking care of the divergence singularity problem when Q goes to zero. If you set it to false, you basically don't take care of that and it will be much slower convergence. And then EXX divergence treatment, we specify the GG by the rescue method. This was explained in the lecture, there are other methods. For example, a spherical cutoff of the Coulomb potential. And there is also for anisotropic systems, so-called beginner sites truncation. So it's not spherical in all directions. It depends on the, if it's anisotropic in one direction it's longer, another it's shorter. So it takes into account this. Okay, but spherical cutoff and beginner sites cutoff, they are not supporting non cubic systems if I'm not wrong. Okay, this is just to remind that there are several types of hybrid functions and different treatments of the divergence in the hybrid function. Digibaldreschi spherical cutoff and cutoff for anisotropic supercells. And if you specify none, then set Coulomb potential simply to zero. But this is not really the best way to go unless you are using this specific type of hybrid, which is GAUPB. Which is not that popular actually. Okay, and now we need to do the exercise. We need to run a hybrid functional calculation with PB0 for bulk silicon. We need to compare the total energy using different Q point meshes and see how the total energy converges as a size of the Q point mesh. So I remind that we're using here a simplification when we replace the one summation over Q points by summation over Q points. But this is an approximation and we need to be sure that we are not making a big error there. So we need to make some convergence tests using or not using the extrapolation at the gamma point. So let's do this exercise during 10 minutes and the remaining 15 minutes we use for wonder wells. So hybrid functional calculations are extremely expensive. Currently in quantum espresso there is efficient implementation for the case of norm conserving pseudo potentials and even is an expert in this implemented very efficient. And recent method called ASE, correct? He's an expert so if you have questions. Yes, the method was proposed by Linlin and then we implemented it in quantum espresso is based on a projection of the exit exchange operator on the manifold. Yes, and this speed ups calculation by which factor more or less. Yeah, it depends on the systems we registered up to 7. It depends on cases if you compare with the very old previous implementation is even an order of magnitude. Whereas in other cases it can be 3, 4, 5, but still it is a significant let's say in fact now it is the default method. So thanks to even you can now run every functional calculations very effectively and efficiently. Instead for ultra soft to the potentials and project augmented wave. Is this supported your implementation even? Yes, it is supported. Yes, it is. Okay, okay, but just notice that. In the case of ultra soft and PAW to the potentials in hybrid, there is less, let's say, more limited support. I mean, some things you cannot do like it's first of all, it's a bit quite slow. And also certain things cannot be done. I think it's like stress forces, things like that I don't remember by heart, but for non conserving there are more things you can do. Then you can do with ultra soft and PAW. So I would suggest probably to try first with non conserving and. And then you can also try to explore ultra soft and PAW. Okay. Yeah, so you use HPC for this case. As usual remote MPI run. X minus input in the name of the file. And you submit a while waiting. One important comment for for hybrid functional equations is there are some tricks. At this point, I should not do that will go to pw dot X input, there is parameter called ecot folk. Where is it? I mean, I have some issues if I try to use some, you see some comments. Certo, c'è l'ultimo, l'ho visto, è la seconda linea. Ok, grazie. Quindi questo parametro è chiamato if cut foc. Il value default è the same as if cut raw. Quindi questo è per... Se ricordate, non considerate il potenziale, per esempio, il value foc è 4x più grande che il genetico e il giocatore della funzione. L'ultra soft e Pw è 8x più grande. Quindi il value foc è anche... è uguale ad il value foc per default. Ma puoi significativamente speedup la calculazione usando il value foc. Per esempio, puoi iniziare a setterla ad il value foc per WFC. Ma io credo che se non ero, potete usare magari un factor di 1.5 o 2 o qualcosa di così per avere un'accursione riuscita. Quindi il value foc è... può dare valori tra il value foc per WFC, il value minimo, e il value default è il value foc. Ma guardate attenzione, se è una coppia molto dispensa, puoi ridurre questo parametro di value foc. Il value foc è 2. E questo era il punto gamma, quindi è finito. Quindi dobbiamo fare la calculazione per altre mesche di Q-point. Non abbiamo tempo per fare questo. E questo è il motivo di aver visto questo risultato. E infatti, vediamo che usando il value foc per Q-point è vero. Perciò la convergenza più forte per mesche di Q-point. Quindi già per gamma questo punto è più vicino al value converso di Q-point. Se usiamo la tecnica di extrapolazione gamma se non, sei più fuori con il value converso di Q-point. Per Q-point 2x2x2, vedete l'improvenzione quando usate la extrapolazione, ma, diciamo, è più piccola. E poi per Q-point 4x4x4 è quasi negligibile e per Q-point 8x8 sei già in convergenza. Q-point 8x8x8. Quindi questo è why it's recommended to use extrapolation technique to allow us to use scarcer Q-mashes. Ok, quindi questo è tutto per i functioni di Q-point. Let's try to answer some questions if we can. We have Pallon in call. Please. Yes. Thank you very much for your talk. I have a question concerning the e-cat fork. So you said that it will speed up the calculation, right? I want to know if it doesn't affect the total energy and something like that. Ok. Yeah, so it affects. So this why you need to also check carefully the value of the e-cat fork. So when you reduce it to the value of e-cat WFC, this cutoff for the wave functions. So this is the minimum value of the e-cat fork. But be careful because perhaps in your system this simplification can lead to some quite big errors in the total energy and other properties. So this is why always before production calculations, before publishing the results, as everything else you always have to perform convergence tests and checks. So this parameter has to be used with care. Ok, ok, thank you. Welcome. Ok, then from Slack there was a question about the hybrids. How can we modify the exact exchange mixing parameter in the hybrid functional? Also can we modify the screening parameter present in the HSE functional? So the mixing parameter we can modify from the input. Recall the name of it. Ivan, do you remember the name of this mixing? So there is some parameter. Screening parameter. Ok, let me open it here. Oh, thank you, Pietro. So that's it. I'm not wrong in something, yeah. Here it is, yeah. E6 fraction. Ok, here, this one, no? Screening parameter. Ah, this is screening parameter for HSE like hybrid functions. Ok, there are two questions. Exate exchange mixing parameter, I guess it is E6 fraction and the other one is screening. Exactly, so apologies, yeah. So the fraction of exact exchange is, I mean, it's called exact exchange fraction. Here it is. So we can change it. And screening parameter is specific to HSE. So since we're talking about this, let's open some papers for the electric dependent hybrid functions from the Galli group. So what they said that this parameter alpha is actually related to the dielectric constant of your material. It's one over epsilon. So, and of course, in different solids the dielectric constant is different. But using global value of 0.25 is not always working well for all solids. So for those solids where one over epsilon is close to 0.25 it should give quite good results. But for other solids where the electric function one over epsilon infinity is quite far away from 0.25 then standard hybrid most likely will not give satisfactory results so you can tune this fraction of exact exchange by computing it like one over the electric constant. You can go further and do even more advanced things. And this is done in another paper by the same group. So instead of using equation five where you have alpha one over epsilon infinity you can use more sophisticated expression in equation seven which is then used in equation six which is this screened Coulomb potential. So in this expression for W it's not simply one over epsilon infinity it's more complicated function where you use a complementary error function it's sort of a modeled expression for the screen Coulomb. It has interesting properties that if you take the limit when R goes to R0 or R prime so basically two points are close to each other in this limit it gives you simply one over R minus R prime but if you are going in the limit one R minus R prime it goes to infinity it gives you one over epsilon infinity basically the first term survives in the limit when R minus R prime goes to infinity and in a limit when R minus R prime goes to zero so this term survives. I don't remember the details I worked on this five years ago but basically if you take the two limits you should reproduce the expected behavior of the screening of the Coulomb potential but these features I think they are not implemented in quantum express currently but these electric dependent hybrid functions are very popular nowadays becoming more and more popular for solids so you can check the literature and let's say last 15 minutes for Wonderwalls but before we go maybe one more question if there is any. Yes, Dorje Esperas. Hello, thank you for your explanations for answering our questions I was thinking how to extract the Havar potential matrix you saw us at the beginning could we can we do this using the STF matrices present in the output or is something more I mean the occupation matrix and the project source Apologies, I didn't understand fully your question so you're asking about occupation matrix in the output Yes can we construct the Havar potential a Hamiltonian term for example to use it in a tie binding using this occupation matrix Oh, okay I don't have experience with this but I think something can be done of course but I don't have experience with doing what you suggest constructing a binding Hamiltonian on top of DFT plus the results but I think yeah, yeah I would say yes Okay, may I ask this in a slack Thank you Thank you too Okay, so let's move on for Wonderwalls functions in graphite so this example last one of this hands on we consider graphite because it has a layered structure we have layers and between these layers we have Wonderwalls weak interactions the experimental value of the interlayer distance is 3.333 x angstrom if you try to use a standard LDA or GGA exchange collection functional you will not obtain satisfactory results because LDA underestimates and GGA overestimates interatomic distances and we will compute those distances in this example so that's why here we need to use Wonderwalls corrections so how to do that with quantum espresso you need to specify there are two ways how to do that depending on which type of about correction you use it's empirical or or other so you need to specify parameter input underscore DFT Wonderwalls slash DF this is just one kind of Wonderwalls correction there was a question earlier about Wonderwalls kernel in earlier versions of quantum espresso you need to generate the Wonderwalls kernel not in this case in others so like RVV10 and then do the calculation but in the latest versions of quantum espresso this kernel is generated on the fly and so you don't need to worry about that and similarly to hybrid functions when you choose pseudo potential it must be the closest GGA because there are currently no pseudo potentials for non-local functions so this is why also you need to choose pseudo potentials wisely so in the previous slide I showed that we can activate Wonderwalls correction with input DFT keyword but also you can do that using a different keyword Wonderwalls underscore correction and there are different types like Grime D2 there are several ways you can specify it it's known as DFT slash D or Grime slash D2 then there is Grime D3 or DFT D3 Tkachenko Scheffler where this coefficients 6 coefficients are computed from first principles while in the first two they are semi-impedical and this XDM and in this example unfortunately we don't have time to do it now but you can do it after the ending of this hands-on session what you can try to do is to perform these 7 calculations by in the first 5 you activate the Wonderwalls correction in different ways the first three using the keyword input DFT this is our non local flavor and the fourth and the fifth one are using the keyword Wonderwalls correction DFTD and DFTD3 these are semi-impedical so in the first three cases it's really the functional that you are changing because input DFT keyword controls the functional while the fourth and the fifth one are just corrections that are done at the end of the CF calculation so that's why we have two different keywords which might be confusing for some of you and point 6 and 7 are standard PBE and LDA calculations and so what you need to do using what you learned yesterday for the structural optimization so-called BC Relax you need to do the BC Relax calculation in these seven ways and then for the converged optimized structure you need to determine the interlayer distances using X-Cris then for example by doing so you should obtain something like that so the vertical axis is the interlayer distance and the horizontal axis is the method we are using actually the functional or correction, Wonderwalls correction in the last two are standard the PBE LDA functionals and the height of each bar corresponds to the value of the interlayer distance in graphite and the horizontal dashed line is the experimental value so we can see that there are big variations depending on the correcting scheme for example we see the first two overestimated RVV10 which is actually Stefano De Gironcholi ha lavorato on this about ten years ago and this is also supported in the phone encode which we will discuss tomorrow so the third is RVV10 which is seems to be the closest to the experimental value and DFTD is largely underestimating it DFTD3 is second best result and we see as we said before PBE largely overestimates interlayer distance while LDA underestimates but we see interestingly that DFTD is even worse than LDA this is quite interesting and surprising at least to me but notice that these results are not carefully conversed because of course we use K-point mesh low cutoff etc etc so these results should be considered as not fully conversed so you should really converse them very well to be sure on the performance of each correcting scheme but at least it gives you an idea that you see big variations depending on which wonder wells correction you choose so this is why as far as I know there is no one single wonder well asking which is the best and everybody should use only that one so what people typically do they test several ones and see which ones perform the best for the system of interest and they continue working with that but please correct me if you have more experience and this was not totally correct so the remaining five minutes we can use to discuss maybe answer some questions there is one question from streaming what is the difference between input DFT and wonder well score in the relaxation and if we mention them in the same time in the input what will happen ok so I think I already explained the difference between input DFT in wonder wells correction because in the first one input DFT you are really changing the functional because input DFT controls the functional we use this input DFT parameter also to specify hybrid functional you remember we said input DFT equals PB0 so it's really the functional that you are changing while here 4 and 5 wonder wells correction these are post processing corrections at the end of the pb calculation you just apply some post processing correction to your total energy and forces and stresses ok yeah so Ivan Giurotto said there are many questions on YouTube apologies that we cannot answer all of them yes some of them we are trying to answer the questions we can maybe with as a short video we are doing it ok if someone else want to raise and for some boys question please one sorry there is also one question on Slack or there is someone wants to ask now on zoom ok so on Slack there is so we use wonder wells to get the interlayer distance right would the wonder wells correction also affect the in-plane lattice parameters in my experience semi-impedical wonder wells also shrinks the in-plane lattice parameter given a worse result than the in-plane is that a common effect would it also be the case with wonder wells functions actually I don't have experience with that maybe step from the journal we can answer if he's still around if not fine we will reply later ok I know it was explained but in case of HSE functional which pseudo potential is best to choose so I would recommend to use pseudo potentials from the SSP library but all of them are ultra soft and project augmented wave however I explained before that for hybrids wish well I would go for non conserving first so it's a yeah it's interesting question maybe you can try to start you can start with SSP library and try to use ultra soft and Pw if this is not working well you can try norm conserving pseudo potentials there are two libraries for non conserving pseudo potentials that I know are very good one is called pseudo dojo developed by belgian group pseudo dojo so this is a norm conserving pseudo potential so you can try them and then there is another library norm conserving pseudo potential library it's called SG15 ok and you can also try them out for hybrids next does the e-cut fork parameter pseudo potential no sorry is it pseudo potential dependent e-cut fork parameter the answer is no e-cut fork is agnostic to pseudo potentials so it just related to Fourier transform it doesn't matter well of course for changing pseudo potentials you can you should check e-cut fork maybe in some cases e-cut fork is smaller in others it's larger so e-cut fork indeed is sensitive to pseudo potentials but it's not explicitly dependent somehow in order to choose in order to choose the best correction that represents best the interaction between a substrate and a metallic slab shall I proceed like in example 3 following the best correction that brings the lattice parameter of the pure metallic material closer to the experimental value my answer is yes but again you can also ask professor dogeronically maybe he has a answer I would try several flavors and see which one is the best same for hybrids you have PB0 you have HSE well a priori you cannot tell which is the best commonly people use HSE but also you can try PB0 and others so the same for GGA functions there is PBE, PBSOL and tons of other functions some are the most popular but of course you can test different ones and choose the best ones for your system ok so I think we should stop here we can continue the discussion on slack if you still have some questions we will try to answer all of them or the majority of them and yeah thank you very much thank you very much Yuri thank you so let me just say one more thing of course we meet again tomorrow at 8.30 for Stefano Maroni we explain the functional distribution theory and another recommendation if you don't have the last version of zoom you might find useful to update it because it might be that the breakout rooms might work better on the last version so yes please you can continue posting questions on the slack channel and thank you again very much Yuri we meet again tomorrow morning bye bye have a nice thank you massimo maybe we can close the connection or put stream server now no e come ivan scusa mi e magari possiamo