 Okay, so let me start thanking the organizers for the invitation and you see the title and the outline of my talk is given here. I will talk mostly about our applications, development of post-DFT methods in the last couple of years. I will give a short introduction, then we'll talk about MP2 and RPA using the Gaussian and plane-based method, so that's the basic methods that we are using and then we'll go to recent developments and if I have time also to some outlook and challenges. So of course this was done by many gifted people. All of this is in collaboration, it was in collaboration with Yostfandi Wanderle and most of the work was done by three brilliant PhD students, Manuel Guidon, Mauro Delbein and Jan Wilhelm. What do we want to do? So we want to do liquids and solutions and we know that after many years of doing water and ions in water and similar systems, there are still some problems with GGADFTs and the question is can we go beyond GGADFT, can we use non-local correlation methods and improve for something like, for example, the density of water where we need a balanced description of hydrogen bonding and van der Waals forces. Another shortcoming that we recently found or it was given to our attention by some experiments at PNNL that, for example, the structure of the solvation shell of ions in water is not very accurately described by GGAs. Or for example, of course, well-known level alignment ions in solution at solid-liquid interfaces are problematic for GGAs. So, for example, if you look at the density of water, we all know if you do this for standard GGA functionals, you get a too low density. If you add now van der Waals corrections in some ad hoc manner, you get better densities but you get a widespread of densities. Depending on the functional you're using, typically are between 1 and 1.1 grams per cubic centimetre, but you can even have outliers like these mini-sorta functionals that are 20% of too dense functionals. So, is there a possibility to have here much better or more predictive values? Another point that as I already said is ions in water. For example, here this is potassium in water. If you do PBE D3 or also a hybrid functional PBE 0 D3, the peak of the radial distribution function is off by almost half an angstrom. So, could this be something that can be done with post-Hartifoke methods or post-density functional GGA methods? Another thing, energy levels in liquid water. We all know the band gap of water and that then is only described in these functionals also in hybrid functionals and this then affects also the levels of ions in water. So, to do this type of simulations, what we want to do is electronic structure calculations and we want to use non-local correlation because that includes the van der Waals correction and we hope to have van der Waals correction and for example hydrogen bonding on the same level so this should add to better accuracy. Of course, if you go to quantum chemistry, there are many methods doing this but we only want to do the most simple ones. For example, Mp2 or direct RPA may be in connection with a double hybrid system. Still, it's a challenge because we want to do system sizes that we can do sampling for liquids so at least 200 atoms in the order of 500 correlated atoms. We are using a Gaussian basis set or a localized basis set so that means for us the order of 4,000 basis functions because these methods are much more demanding on the basis than GGA functionals are so we are typically in the order of 4,000 more basis functions. On the other hand, we want to do liquids so of course we have periodic boundary conditions but the gamma point is hopefully enough so we concentrate on gamma point only calculations. Simulations are not enough. We have to do sampling. Sampling can be done either by molecular dynamics or Monte Carlo. These methods are expensive so you better do the best you can meaning multiple time step schemes or very accurate bias potentials so that you can get a long way without doing a calculation on the highest level. In all calculations showing this has been done. Then, still we need smooth energy surfaces. Monte Carlo will not work if you have a racked energy surface or accurate analytic forces if you want to do MD. So if we have these methods, they have to be analytic. We can cut corners but not all corners can be cut. We need in the order of 20,000 for water. This is now water at ambient temperature. We need in the order of 20,000 energy calculation. So this means 20,000 RPA calculations gives us one density of water. What does that mean? Nowadays, if you go to praise, you probably get easily from praise. CPU budget of one million no dollars. If you want to use this for this type of calculations and you want to spend it in three months, you don't want to wait too long, and you run on 512 nodes that means roughly one RPA or MP2 calculations in six minutes. So that's our goal. Now how do we do this? All of these calculations need at some point electrostatic integrals, two electron integral types of occupied times virtual orbitals interacting with occupied times virtual orbitals. Now we use a Gaussian auxiliary basis and an RI approach, resolution of identity or density fitting approach that these type of integrals are always calculated from this formula here where we have here three center type integrals and then here the inverse of the Coulomb matrix. This integral is calculated using the Gaussian and plane wave approach, meaning by Fourier transforms, and then these integrals can be calculated as a tensor contraction over these three index tensors where the three index tensors themselves are calculated using these GPW integrals that meaning we are calculating these integrals by putting these functions as these are our auxiliary basis on the grid doing an FFT, calculating the real space potential of this orbital and then integrating over the product of two Gaussian functions. This is the GPW approach. We are very good at doing this and then calculating this quantity, and this quantity then has to be transferred with your molecular orbitals to the BIAS. Once we have this, we can plug this scheme into the MP2 energy, you see you have these type of integrals, you calculate them using this tensor contraction, that is in fact the most expensive part. If you look at it, it's a contraction over the auxiliary basis that makes it an n to the fifth step, so we don't reduce the scaling, we just make it faster. So it's still an n to the fifth step for MP2. The RPA formula given in this form includes the Q matrix that is a matrix of the auxiliary basis and it has tensor contraction over the particle whole states and this makes it an n to the fourth step, so RPA calculations in this form can be done on n to the fourth steps. Now, applying these methods to our system, what we are doing is isobaric isothermal Monte Carlo simulations are liquid water, 64 molecules, 192 atoms. The basis set for the experts is a correlation-consisted triple-seat surveillance basis set, roughly three and a half thousand basis functions. The RI basis is 8700 basis functions. A little bit on scaling, these methods are very expensive but they are also very easily to, or rather easily to, parallelize, so parallelization up to several ten thousands of cores for this specific system is very efficient. Now, you do the scale, the sampling and in the end, you can look at how much CPU time I have burned to get my density of water. So you start with a GGA, you need roughly, we need roughly 100,000 CPU hours to get the density. Hybrid functions, depending on how you do it, is a factor of 3 to 20 in our code, more expensive to get the density and then RPA and MP2 are even more expensive, typically 50 to 100 times more than a GGA calculation. So that's what you have to invest to do this type of calculations. Now what do you get? Here is the Monte Carlo profile, the Monte Carlo Cycles profile of the density, MP2, some oscillations of course and fluctuations but it converges to a density of 1.02, that's 1.02 much, it's 1.02 grams per millilitres or not bad. RPA gives you a slightly different than lower density 0.99 grams per millilitre. So very good results, seems to be that really this non-local correlation and being doing everything on the same level is a way to get very good results for these systems. What about the sodium in water again? We get a very or rather different radial distribution function and the peak is now much closer to experiment again. Energy levels in liquid water, they are better, still not perfect, this is certainly something where we don't get away easily with just using MP2 and direct RPA. This was the standard implementation that has been done a couple of years ago, mostly by Maugo and in recent years we have made some progress on these codes. First, we have forces and stress tensors for MP2, both restricted and unrestricted. We have also implemented GW, once you have RPA, the step to GW is not that big. This has been done by Jan Wilhelm in the last couple of years and very recently implementations where we reduced the scaling of the RPA and GW from N to the fourth to N to the third. Now, some information how expensive are MP2 forces or RPA forces, these are now non-variational methods, so you really pay for forces. Typically, a factor of four. So the force is four times more expensive than the energy in these type of calculations. It's the same for close-cell and open-cell calculations. The additional factor going from close-cell to open-cell is roughly a factor of three. You have to do three independent energy terms instead of one in MP2 going from close to open-cell. With this, it was possible to calculate an MD trajectory of water and from the MD trajectory to calculate the infrared spectra of water. So what you hear see is the experimental and the MP2 infrared spectroscopy calculated from a velocity auto-correlation function on an MP2 MD trajectory. With some tricks, one of them is that it's scaled by 0.95 that all the quantum chemists will say, of course MP2 scaling factor for infrared frequency is 0.95 and we also use this. This seems to be a universal law because MP2 is not the final answer. There are still many types of correlation to be added. The GW code, just to show you, here the scaling with system size and here the parallel scaling with system size, if you look here in this double logarithmic plot, it's also n to the fourth. So GW in this implementation is n to the fourth like the underlying RPA calculation. Roughly the same timing. So these calculations with 20 quasi-particle energies took about, for a system of 64 waters, take about a couple of minutes on a machine of this size here, few thousand processes. Now last and newest development is a cubic scaling RPA. So here the idea is that we are using another resolution of identity that is no longer based on the Coulomb metric but now on the overlap metric. Integrals can be calculated through this contraction, tensor contraction. The advantage is now that these tensors here, these three index tensors are overlaps between three different atomic functions and so they are extremely sparse. Only when all three functions overlap, they are non-zero. So these are extremely sparse quantities here and we can make use of this sparsity. What you need is the three center overlaps, the two center overlaps, and you need here the integrals and I forgot to say, of course, whenever I write here an integral, in principle I mean a periodic local function. So it's an atomic function but it's periodic. Fully periodic means all these integrals, in fact these Coulomb integrals are in fact evolved integrals. And for this type of evolved integrals over Gaussian functions we now have a semi-analytic integration so we don't do this anymore with Fourier transforms, we directly calculate these integrals analytically. And this gives you a big speed up for this type of calculations, reducing the pre-factor again, not the scaling of the method. So we are using this sparsity of the three center overlap integrals and we need an additional trick. That is we are using also now a time, or we are going to the time domain and we are calculating the Q matrix in the direct IPA now in the time domain or with other words we are doing a Fourier transform. If you do a Fourier transform you can separate the energy part here and you can derive here a quasi-density matrix or in fact it's a Green's function that only depends on the occupied orbitals and one that only depends on the unoccupied orbitals. Now this is the final formula. If you go through it you make use of all the sparsity. This is cubic scaling. Cubic scaling because in the end you have to calculate the logarithm of this matrix. This matrix we don't assume it's sparse so you do a logarithm of a full matrix that's n cubed. That's where n cubed is coming in. Here is the comparison. The cubic scaling code compared to the linear scaling roughly 128 waters is the break-even point in the current implementation. The interesting thing is also that in fact the cubic part is not yet setting in and for all of these systems up to 1000 water molecules in fact the code is quadratically scaling. Even better the quadratically scaling part becomes linear scaling once this density matrix is also become sparse. The only cubic scaling part left working with this matrix is all the rest is in fact becoming linear scaling. The same can be done for GW so this is now a cubic scaling GW calculation. You already heard before of these nano ribbons so what Jan together with the people at EMPA he has calculated these nano ribbons the GW calculations and you can see here for these quasi 2D systems the scaling is even better and the break-even point is at roughly 100 carbons and it could go up to 1,700 nano ribbons with 1,700 carbon atoms. Almost finished. I also almost finished in time. A little bit of outlook and challenges that we have basis set convergence how to calculate properties for these type of functionals and can we do sustainable code development once we go in these type of complex systems. So basis set convergence I just wanted to say if you ask a quantum chemistry we'll tell you it's solved just use an F12 algorithm unfortunately nobody ever thought or I don't know of anybody of thinking of F12 algorithms for RPA and GW especially not once the tires are low scaling as the implementations we have. Another possibility would be maybe we should go for double hybrids with only long range not short range and then the basis set convergence is not a problem. The other thing is we need an additional basis set Gaussian basis sets are already something not that easy to handle especially if you come from a plane wave a community where you have one number and suddenly you have 1,000 different basis sets which one to pick and now comes additional RI basis sets so what do we do? We can automatically generate RI basis sets but they are typically two times bigger than the ones that we do by hand for the same accuracy but maybe it's nice to go this way anyway CPU time is always there and human time is very valuable. Properties as I said derivatives and everything we have non-variational wave functions it's hard to do and question is how do we do properties any property a response property is already a second derivative and so on so we have to do second order response it's something we have to think about also periodic boundary conditions so for example how do you calculate the dipole moment of mp2 in periodic boundary conditions Sustainable code development because I'm at the end I can just see everything means increased code complexity if you decrease the complexity of your algorithms you usually increase the complexity of your code that's something that is almost universal the same is true also for hardware the more sophisticated you want to go in your hardware usually your codes become more and more complex something that is difficult to be sustainable in the future I don't know I don't have a solution for this so with this at the end I'm happy to take questions thank you