 The lowest energy here, eigenstate and energy eigenvalue of a many-body fermionic Hamiltonian. And actually, we've subsequently since generalized our algorithm for the calculation also of excited states. So later on, we'll deal with that problem as well. So I'll spend quite a long time on what I call introductory remarks, because I have to set the scene as far as actually describing the problem and how it's actually written down, at least in quantum chemistry. So the notion of orbitals, second quantization, and so on. And then I'll describe this algorithm that we've developed, SCI QMC. Oh, that would be great. Thank you. So if you're wrong with the microphone, we can use this thing to change it in just a moment. So the bulk of the lectures will be concerned with this algorithm, SCI QMC, but which actually comes with a slight tweak to it that we call the initiator for SCI quantum Monte Carlo method. And the initiator is a systematically improvable approximation to it. So what we will show is that SCI QMC can generate in a Monte Carlo form the exact ground state, eigenvector for a fermionic system, but is expensive. And then this tweak, the initiator method can provide systematically improvable approximations to that. And so as the number of walkers increases, we can demonstrate that we can systematically converge onto the exact solution that SCI QMC provides. And then I'll provide some actual online demonstrations today so we're connected to my PC sitting under my desk in Stuttgart. And then we'll actually run some calculations. It's very important to get real-time experience with this method because it's a lot of things that I talk about seem as though these are going to be incredibly expensive calculations. And in fact, they are very inexpensive. And that's just important to see how it goes and also to see the nature of the fluctuations, how quickly things converge, et cetera, et cetera. So I will do two sets of online demonstrations. One with the basic algorithm and one with what's called a semi-stochastic adaptation of the algorithm, which greatly reduces the stochastic errors. And then we have a series of more advanced applications that will hopefully come onto probably tomorrow on calculation of properties, density matrices, excited states, self-consistent field optimizations, what the quantum chemists call multi-configurational self-consistent fields, and so on. So we'll see. I don't know how much exactly of this we will go through. And of course, the other thing is to mention is if you have any questions, then just ask them as you see fit. OK, so let me just point out that the grand challenge problem that we would like to be able to solve is the many electron Schrodinger equation. This is here in non-relativistic form. And n being the number of electrons that we want to treat. And this, of course, is the realistic Coulomb interaction which we generally have to face with in ab initio calculations. Of course, if you're doing lattice models, Hubbard models, or such like, you have analogous interactions. And so this is just indicated here as 1 upon R1, 1 upon Rij. But you could replace that with anything else you like. And this is the external potential. And again, in ab initio calculations, we would assume that we're working in the Born-Oppenheimer approximation. So you have clamped nuclei instantaneously. So the nuclei are assumed to be fixed for the purposes of these calculations. And they generate a Coulomb field which somehow keeps the electrons bounded. And we're then interested in solving the correlation problem that arises therein. And of course, there's a kind of a frustration involved because the nuclei suck the electrons in. And the repulsion, the R1, the Coulomb repulsion tries to keep the electrons apart. The kinetic energy leads you to itineracy, delocalization, and so on. So there's a huge amount of what you might call frustration implied by the Schrodinger equation. And that's what makes it so interesting, actually. So this is the eigenvalue problem that we would like to address. And in general, at least from the point of view of, let's say, where we want to start, we're going to be interested in the lowest energy eigenvalue of the Schrodinger operator. But it's not the absolute lowest energy eigenvalue that we're after. We're actually after the lowest energy eigenvalue corresponding to the eigenstate that is fully anti-symmetric. So it's relatively easy to show that there will always be an eigenvalue corresponding to a symmetric wave function which is a lower energy than the fermionic one. That's the bosonic ground state. And in Coulomb systems, the bosonic ground state is way, way, way lower in energy than the fermionic one. And that causes what's called the fermion sign problem that if you generally perform an unconstrained Monte Carlo simulation of the Schrodinger operator, you actually immediately collapse the projector techniques, immediately take you towards this very low-lying bosonic ground state. You can think of the bosonic ground state as where all your particles are sitting nicely in the 1s core orbital, which is usually a mile down. And it's actually preventing that instability that is the major problem. So the fermionic ground state is usually way above. So the thing to point out is that the many particle wave function that we're dealing with is obviously a very high-dimensional object. So that, if you like, is the primary source of difficulty. And that would be difficult even if you were dealing with a bosonic system, actually. So you have high-dimensional functions to somehow deal with. And for fermions, for electrons, say, each of these coordinates refers to three spatial coordinates and one spin coordinate. So our electrons are spin-half particles. So they will have one spin degree of freedom, which has got two values, let's say alpha and beta, that it can manage. And yes, so this is obviously a very high-dimensional function. And that is the primary source of difficulty, how to very accurately estimate these wave functions sitting in very high-dimensional spaces. The second difficulty comes from the fact that we have anti-symmetry to deal with so that if you hear exchange, you see the labels i and j, you end up changing the sign of the wave function. And that is the anti-symmetry. And that holds for every pair of exchanges. And then the other thing to point out in this preliminary discussion is that you can see there are no fundamental constants in the Schrodinger equation as we've written there. So h bar has been set to 1. The mass of the particle, the electron mass, has been set to 1. And the electron charge in absolute value has been set to 1. And this actually sets us our unit of energy, atomic unit of energy, which is one heart tree, 27 electron volts approximately. It's a huge energy, in other words. In chemistry, 27 electron volts is you would dissociate the strongest bonds out there, in other words. You can think of strong bonds like N2 as being 10 EV. And of course, often in chemistry, we're actually interested in describing very subtle energy changes. And so we're often interested in actually calculating these energy eigenvalues, or at least differences between energy eigenvalues, which are out there in the second or even third decimal place. And that's an extraordinarily difficult thing to do given the complexity of this problem. So one can set up very relatively sophisticated calculations starting up with Hartree-Fock theory that will capture even 98% to 99% of the energy, and yet is chemically totally useless by and large. The remaining 1%, you have to still be able to calculate to 80%, 90%, almost 100% accuracy before you're able to say anything useful about chemistry. So that is the range of difficulty that we're facing here. OK. So of course, an enormous amount has been researched about this problem. And there are different schools of thought, which are out there. And so when we talk about ab initio, what we really mean is starting from, let's say, the Schrodinger equation. In other words, starting only from these fundamental constants, and then we want to obtain some estimate of our ground state energy. So that is the target. Of course, you have to make approximations. So ab initio does not mean it's missing an approximation, but what it really implies is that we're not going to introduce any semi-empiricism into the problem. Now, if whether this is a feasible task in general, that, I think, is the jury is out there for sufficiently complex systems, it probably isn't. But anyway, that is our aim. So one school of thought, which I will call quantum chemical, the physicists would probably call this many body theory. But essentially, they're the same. And that really starts off with making two types of approximations to begin with. So in the first place, we introduce finite basis sets. So in other words, even at the level of solving a single particle, one electron problem, we already approximate that problem. And so in the simplest case, you could just imagine a basis set that is a discretization in the form of a lattice. And you can make the lattice finer and finer and finer. And that will eventually convert you to the solution of the single electron, a Schrodinger equation. That's actually not a very efficient lattice basis set. And how to construct accurate solutions whilst minimizing the number of basis functions is a problem that the quantum chemists, among others, atomic physicists, et cetera, have studied, have long studied. And there are different classes of basis functions that are out there. We will, in the calculations that I'll be showing you later on, be using Gaussian basis sets, which are somewhat non-ideal from one respect. The Gaussians don't have the correct cuss behavior at the nuclei. They don't decay correctly into the vacuum. They decay too quickly. Slater functions would be better. But Gaussians have got very nice mathematical properties. And in particular, it means that the relevant integrals that we need to be able to do the many body calculations in a Gaussian basis set, they're analytically solvable. Whereas in other basis sets, they're not, apart from plane waves. So Gaussian basis sets will form our choice of basis sets. And so when we talk about Gaussian basis sets, we're often dealing with functions whose radial behavior decays as a Gaussian. And then the angular behavior is given by the usual spherical harmonics. So you can kind of build complete basis sets out of those type of functions. And the point about Gaussians is that they're atom-centered by and large as well. They don't have to be. But generally speaking, you end up placing your Gaussians at the centers of your atoms. So that is the type of very commonly used basis set in quantum chemistry. If you're doing condensed matter physics in periodic systems, you would typically use plane waves, which is another way of filling your space. Those are usually not tied to the atoms. And plane waves are useful in some respects. But unfortunately, it's a rather inefficient basis. It's a nice basis in the sense that you can take it to completeness relatively easily. But it's inefficient because the sides of the basis sets tend to get gigantic very quickly. So quantum chemistry in plane waves is something that we've also looked at, but it's not really very practical. But I'm going to assume that the basis set problem is kind of solved. You may or may not agree with that, but that is the point of view that at least I'm going to take. And there are problems with basis sets, but that's a separate issue. But the real question comes in the second set of approximations, which are the many body approximations that we have to be able to treat. That is to say, once you've assumed that you've got some representation, degrees of freedom, which express how the single particles can move, then you have to deal with the problem of the fact that the electrons do not move in an independent fashion. If they moved in an independent fashion, then already Hartree-Fox theory or mean field theory would give you the exact solution, and you wouldn't have to go any further than that. But in fact, they don't move in an independent fashion, and that is where the many body approximation comes in. And so within quantum chemistry or many body theory, there are then entire hierarchies of theories which try to deal with this problem of correlation. And they almost always start off from mean field theory, which is Hartree-Fox. And then they build up solutions based on mean field theory. So you've got, for example, many body perturbation theory. You've got diagrammatic methods of which couple cluster theory is probably the most systematic way of doing it. And so I won't be talking about these. They're probably going to be dealt with in other lectures, actually. But if you take these methods to essentially infinite order or as far as it can go, and if they don't diverge, then you end up with what's called full CI. And that's just the exact solution of the Schrodinger equation in the basis set that you hypothesize to begin with. So it's not the exact solution of the original Schrodinger equation, of course. But of course, as the original basis set itself is then taken to completeness, then the full CI solution does indeed provide variational upper bound. So you gradually get closer and closer to the final exact solution. So that is, in a sense, what most people would call the ultimate limit. Unfortunately, I should say the rate of convergence to the true exact solution as you expand your single particle basis set turns out to be exceedingly slow. And that's because, as we'll see later on, the nature of the two particle wave function, as two electrons get closer and closer together, and generally speaking, it's two opposite spin electrons as they get closer and closer and closer together, the wave function develops a cusp. In fact, it becomes exactly a cusp at coalescence. And those cusps are extremely difficult. They're non-analytic objects, essentially. They're extremely difficult to represent in this kind of analytic expansion. And that causes all sorts of problems. And again, we'll touch upon how to improve the rate of convergence in those cases. OK, so this is a really beautiful hierarchy. I don't know how many Nobel Prizes have been awarded for this, but it's certainly quite a number, let's put it this way. But the problems are yet far from resolved. And actually, I would say that the most interesting problems out there all hinge in how to go from Arctic Fog through to full CI. And that's a really good area for research, if you're interested. So the nice thing about this hierarchy is that it's systematically improvable. So in other words, if you have an exceedingly accurate experiment that you're trying to somehow hit, which really describes very, very subtle features of the problem, then in principle, you can get there by ramping up your basis set, by ramping up the level of the many-body approximation that you're dealing with. And as long as you're dealing with the relativistic, so non-relativistic problem, you will always get there in the end. And even, actually, I should say, extension to relativity itself doesn't really cause any serious problems either. So even that limit is approachable. So systematically improvable, but the problem is that it's very expensive to do. So in practice, you cannot actually do it. So it's like one of these promised land that's out there that is just never really achievable. It's expensive. And how to make these expensive calculations practical that is the interesting thing. And a lot of it actually boils down to smart approximations, smart algorithms. In other words, given a computer, how should you actually best use that resource that you have? That is far from obvious, actually, as a proposition. And OK, so that's one thing. So density functional theory, I won't be saying much about, but this is obviously another huge school of thought. It's, again, probably going to be covered by other lectures in this workshop, sort of. OK, well, I would simply say that although it's formally an exact theory, for the ground state at least, it's predicated on the exchange correlation functional. So we have to hypothesize an exchange correlation functional, which somehow captures all the subtle many-body effects that these theories here try very hard to actually describe. So these exchange correlation functionals, the problem with them is that it's essentially introduces an uncontrolled approximation into your theory. So if it works for your system of interest, then that's fine. You have a nice scaling theory. It's generally quite cheap. You can scale it up to many particles, et cetera, et cetera. But where it fails, then I'm sorry, you're stuck. And it usually fails, actually, in precisely the systems that, as I happen to know myself, I started off in density functional theory. And so it usually fails most frequently where you have electronic configurations that are not dominated by a single determinant, but are inherently multi-determinant. And those are like strongly correlated systems. And you get these strongly correlated systems quite ubiquitously, actually. It's not just in very exotic materials, but every time you break a chemical bond, you typically go into a state in which the Hartree-Fock determinant is not dominant. There will be at least another one out there. And these methods then tend to struggle. But I should say it's very widely used as a method. And then finally, we have quantum Monte Carlo, which is a relatively niche technique. So out there, you have many thousands of people doing quantum chemistry, many tens of thousands of people doing density functional theory. But probably the number of people doing quantum Monte Carlo can be counted in the hundreds, something like that. It's a very niche technique, so to speak. It's never really caught on in the same way quantum chemistry and DFT have sort of caught on within electronic structure theory. And so we'll see. I have my own views as to why that is. But intrinsically, the ideas behind QMC are very, very powerful. And in a way, the way they've been implemented historically has been to not to invoke basis sets. So usual quantum Monte Carlo calculations work in continuum space, real space. So your particles are not discretized in any practical sense. They can move around in continuous space. So they have one advantage over these quantum chemical methods, that they don't invoke any basis sets. And so the famous CUS problem, which arises when you're using finite basis sets, doesn't arise so much in these continuum space methods. But unfortunately, because we have the problem that we're dealing with fermionic wave functions, so that our functions are necessarily negative in some places and positive in others, you have to impose some constraints in the sampling. And that introduces fixed node approximation. And it turns out that these constraints introduce another uncontrolled source of error. They're extremely difficult to remove the constraints or to push them, if you like, towards the exact constraints where the nodes of the ground state wave function are. And no one has any practical means of doing that. And so this is another source of uncontrolled error in conventional QMC, which actually I don't think is that different when you look at it in practice from DFT calculations. So people have done very extensive benchmarks of these. And despite the fact that these calculations tend to be a hell of a lot more expensive than DFT, you don't actually see the benefit at the end of the day in terms of dramatically reduced errors, for example. That's the sort of thing. If you're spending 1,000 times more compute power, then you would certainly expect to see some benefit in that regard. But by and large, that's not true. But now, one of the interesting things is that I've actually highlighted these three as three different subjects almost. But actually, the interfaces between these are very, very interesting. And a lot of the research that currently goes on is looking at the interface between quantum chemistry and many body theory and density functional theory. So if you like how to improve these exchange correlation functionals by borrowing ideas from many body theory, that's a big area of research. And that's fine. But another interface that's very interesting, essentially rather unexplored, is the interface between quantum Monte Carlo and quantum chemistry. So in other words, using stochastic ideas to actually solve these types of problems. And from an algorithmic perspective, this turns out to be really interesting, actually. Plenty of space, if you like, virgin territory out there to investigate. And so what we're going to be doing is actually looking at QMC in the context of solving the full CI problem, which is probably the most interesting one, actually. But one can also think of there's also quite nice work from actually originating from my own students who've gone off and done their own things, so to speak. But how to solve couple cluster equations using QMC as well. And that's another very interesting. It becomes technically much harder, but potentially very rewarding. OK, so let's now get to some more detail. So the first thing to introduce is the notion of slated determinant space, which is the appropriate Hilbert space if you're dealing with fermions. So fermions have this anti-symmetry property that is absolutely fundamental. So if you're doing calculations on fermions, the one thing that you cannot allow to be violated at all costs is their anti-symmetry. That's the name of the game. And in a way, the way the quantum chemists have dealt with this is to introduce functions which are manifestly every single one of them anti-symmetric. And then you're going to use that as a basis to expand the wave function that you're going to be interested in. So this is entirely analogous, by the way, to setting up second quantized algebras where you've got creation and annihilation operators which obey the fermionic commutation or anti-commutation relations, they're entirely analogous things. So what we suppose is that we have in our hand a set of spin orbitals or 2M spin orbitals or M spatial orbitals, OK? So you can think of it if you're used to dealing with lattices, lattice models. You can think of a spatial orbital simply as one lattice site if you're dealing with a Hubbard model. Each lattice site would correspond to an orbital. You can think of it if you're a quantum chemist, that would correspond to, let's say, having an S function on each lattice site. But of course, this is a much more general notion. You can deal with symmetry adapted functions which obey symmetry properties of the underlying Hamiltonian or whatever. So we're just going to keep it abstract for the moment. And so we're going to say we have 2M functions u1, u2 up to u2M which somehow somebody has given us. Now that can either be the specification of the problem as it would be in the case of the Hubbard model or it could be a quantum chemistry calculation where somebody has done a Hartree-Fock calculation and those Hartree-Fock orbitals are then our set of spatial orbitals. And for each orbital, we then associate either an alpha spin or a beta spin, and that then becomes our spin orbital. So that is how we generate 2M functions. And so when you're doing your Hartree-Fock calculations, you have to specify your basis set to begin with. That specifies the number of spatial orbitals. And of course, if you do a more refined Hartree-Fock calculations, then the number of spatial orbitals increases and you get a better resolution of your single particle description. OK, now we want to do a calculation with N electrons. So we have M spatial orbitals, so M can be let's say 100. And now we want to do a calculation with N electrons. Let's say we're doing the nitrogen molecule. We'll deal with that later on. So nitrogen has actually got 14 electrons. And it's got 14 electrons, but four of the electrons are core electrons. So they sit either in the 1s orbital of this nitrogen atom or the 1s orbital of that nitrogen atom. They're a long way down in energy. And we can assume that they're frozen. Of course, they're actually not frozen. And in practice, you could correlate those if you like. But let's, for the moment, deal only with the valence electrons of nitrogen. And that would mean that we have a 10 electron problem. And furthermore, we're going to say, well, actually, it's a spin unpolarized system, the nitrogen that's sort of spinning around us right now. So each one is spin unpolarized. So we're going to say that of those 10 electrons, five of them are spin up and five of them are spin down. And so that is the problem that we face. Now we're going to ask, well, how many of these slated determinants could we actually construct given our M spatial orbitals? So we have M spatial orbitals, 2M spin orbitals. And what I show here is a representation. Each one of these lines is a representation of our spatial orbitals. And each one of our spatial orbitals, if you like, can take, can accept either an alpha electron or a beta electron. So here goes 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. If you're doing Hartree-Fock theory, this is how you would actually fill up your orbitals to construct the lowest energy eigenstate. So again, I'm assuming that these lines, as I've drawn them, these are the Fock energies. They've actually been given to us by our Hartree-Fock calculation. And that's, by the way, no real bottleneck at all. OK, so that's a specific representation. And this thing here represents one slater determinant. So essentially, you would be creating products of these orbitals and then anti-symmetrizing them. That produces you a given slater determinant. And this normalization, 1 over the square root of n factorial, is the normalization that you need for these determinants. And now you see you can start exciting electrons out of the Hartree-Fock determinant. So for example, you can create a double excitation where you excite these two spin down electrons, leaving two holes in the original Hartree-Fock determinant and creating two particles in the virtual manifold. So we'll often talk about the occupied orbitals, meaning occupied in Hartree-Fock. These are the occupied orbitals. And these orbitals, which have a higher energy, are the virtual orbitals. And the correlation problem that we're going to face is how to excite electrons from a Hartree-Fock determinant out into this virtual space, which create higher energy determinants. And how these higher energy determinants remix with the original Hartree-Fock determinant to produce you a lower energy overall solution. So this one you see is a slater determinant with two holes in the original Hartree-Fock and two particles. And this one is four holes and four particles and so on. And so you can imagine you have a kind of, I think of it as like an onion shell at the center of your Hilbert space is your Hartree-Fock determinant. And then you have a set of determinants, which are single excitations. So you've just excited one electron. And then you have a set of determinants, which are double excitations, d i j a b. So this notation means you've emptied orbitals i j and you filled orbitals a b. And then you have a set of determinants that are triple excitations all the way out to n-fold excitations, where you've taken every single one of these electrons and you've put them into the virtual space. And at least the first approximation is simple problems. This hierarchy that I've pointed out is a reasonable hierarchy describing the importance of the determinants. But I should say that that's a very, very crude first order approximation to it. So there will be determinants, which are, let's say, sitting out here, which are quadruple excitations of the original Hartree-Fock, that are actually much more important than a whole set of determinants, which are sitting closer in energy on the double excitations. And one of the problems that we face is that we don't actually know, at the end of the day, which are going to be the important determinants from the point of view of correlation, which are the unimportant determinants from the point of view of correlation. The only thing that we know is that we have a hell of a lot of determinants to deal with. Now, how many of them do we have? Well, from a strictly mathematical point of view, we have n alpha, let's say, spin up electrons. And there are m orbitals that we can distribute our n alpha electrons in them. And so the number of ways that we can do that is this binomial coefficient, m choose n alpha. And the same we can do for our beta electrons, m choose n beta. So the total Hilbert space that we're dealing with is this product. And so this grows extraordinarily quickly, gets out of hand extremely quickly. So for this problem, here, with 100 spatial orbitals and 10 electrons spin unpolarized, you end up with 10 to the 16. So if you wanted to store the one vector, if you like, of this problem, you would require a vector with 10 to the 16 words associated with that's about 80,000 terabytes. And so that's before you've even done anything actually to manipulate the problem. And so that, if you like, is the nub of the many part, many body problem, the underlying Hilbert space just gets out of hand extremely quickly. And that's how we deal with that problem, ultimately, that is what we try to face. So in full CI, we don't try to confront that problem, but we just state it. So we're looking for psi 0, which is the variationally the very best lowest energy wave function, which is simply going to be a linear superposition over all of these, let's say, 10 to the 16 objects. So in full CI, the dynamical variables of our theory are these CI coefficients. So these are the objects that we will try to vary if we can in such a way as to get the lowest energy solution to this problem. And because these DI's are already anti-symmetric, it will give us the lowest energy anti-symmetric solution. And that is a statement of the problem. So in a way, it's easily stated. It's a ground state eigenvalue problem in an exponentially large space. That is, in essence, the problem that we face. Now, these objects here, this is the matrix element of the Hamiltonian between a pair of slated determinants. Now, I'll mention a little bit about that. To somebody who's just coming into this, might say, well, crikey, that's going to be a very difficult problem to actually evaluate. Fortunately, those matrix elements come from a computational perspective, turn out to be very easy objects to evaluate. In fact, you can calculate them in order one operation. So that's a good thing. But the point is that there's no rhyme or reason to these matrix elements. They can be positive. They can be negative. They can be large. They can be small. They can be all over the place. And the very fact that they can be positive or negative, in turn, implies that these CI coefficients can themselves be positive or negative. In fact, you can show that, for example, if all these Hamiltonian matrix, if all the off diagonal matrix elements were all negative, then the corresponding CI vector would be entirely positive. And actually, there'd be no sign problem in that. And you could do straightforward Monte Carlo calculations on that. But unfortunately, that's not the case. We have to treat the general problem in which these can be positive or negative. And actually, somebody might ask, well, what about complex Hamiltonian matrix? And the answer is, yes, they can even be complex. And that's also a problem that we have to be able to deal with as well. But I should point out that the sign problem that we face here is not the fermion sign problem. So we're not having to deal with the fact that we've got a system that's unstable with respect to a bosonic solution, which, as I say, bosonic ground state is much, much lower in energy. So it's a sign problem, but it's a weaker sign problem in general. How large can we actually go if you're trying to do this? The answer is about 10 to the 10 determinants. You can probably do about 10 to the 11 determinants nowadays if you really push the boat out. So that means this problem, 10 to the 16, is completely out of range for the foreseeable future. And this is also a very modest problem, by the way. Of course, we want to deal with more than 10 electrons, et cetera, et cetera. OK. Now let me just point out some properties of full CI wave functions, some of which are obvious and some of which are less obvious. The first one is that this is a variation of the minimum energy achievable within the given one particle basis. So it produces a useful, actually, upper bound on the exact solution, number one. Number two, less obviously, is that this total energy is invariant orthonormal transformations of the orbitals. So that means that if we just go back to this statement, so I emphasize that we were dealing with Hartree-Fock orbitals, which is a specific linear combination of our single particle basis sets. But actually, the full CI energy is invariant to that. So I could have chosen a different basis set, a different linear combination. And that wouldn't change my energy. It would change my solution, of course. So the CI vector would change as I change my single particle representation. But the energy doesn't change. And a very interesting question is, from the point of view of full CI, what is the best single particle representation? In other words, which choice of orbitals produces the most compact form of wave function? And I don't think the answer to that is actually known. There are many conjectures onto it. And we can test these various conjectures. Sometimes they work, sometimes they don't. But anyway, that's a very interesting question. Maybe I'll come on to that later on. So that's one thing. Now, the other thing is that our wave functions are pure spin eigenfunctions. In other words, if you apply the s squared operator to it, you'll always get a function of the type s, s plus 1, times the original function. I should point out this is true in the case of non-degenerate wave functions. If you have a degenerate wave function, then you don't have to be a pure spin eigenfunction, even if you're the exact solution. But if you're non-degenerate, then the corresponding solutions are pure spin eigenfunction. So in principle, if I give you a full CI wave function and a non-degenerate, it has a well-defined spin associated with it. And you could discover what that spin is through this by applying this s squared operator. Another very important property is that our energies are size consistent, which is a kind of a special case of size extensivity. So what that means is that if I've got a system which actually consists of two widely separated fragments, so a and b, but these are actually widely separated, and I do a full CI calculation on that system, I get some energy. And then I go away and I do full CI calculation on each individual fragment in the absence of the other one, and I get the corresponding energies of the fragments, ea and eb. Then if I add ea and eb together, I get what the original full CI told me I should get. And that's actually an extremely important property, obviously in chemistry. If you're trying to dissociate a bond, you absolutely need the size consistency condition fulfilled. But it's a very, very tough constraint on the approximations. So it's very easy to construct size inconsistent theories, actually. And one type of size inconsistent theory is that if you take your original Hilbert space and truncate it at some excitation level. So if you take your original Hilbert space and take only, let's say, the set of determinants up to double excitations of the Hartree-Fock and you calculate the energies in that restricted truncated Hilbert space, you find your energies are not size consistent. So if you do this calculation, let's say for a nitrogen molecule, you will not find that the energy of two nitrogen atoms at a long separation will not equal the corresponding energy of one nitrogen atom multiplied by two. So this size consistency is obviously an obvious property that should always be satisfied, but there are many approximations which violate it and they're usually not useful approximations. So that's good. So these are four desirable properties, but the undesirable property, which are full CI wave functions, but this is actually more general of any kind of orbital expansion of this type, is the fact that the wave functions do not satisfy the electron-electron cusp condition. So the wave functions do not linearly, so as two electrons coalesce, the wave functions do not linearly decrease in magnitude, as that occurs. They always have a rounded... So there's two electrons, so if you imagine you fix one electron and now you scan, you ask for the amplitude of the wave function as you're moving the second electron through it, the wave function, so this is psi of r12, the wave function itself for alpha-beta correlation goes through a cusp, so that when the two electrons are actually on top of each other, you have a non-analytic behavior at that point. And that's actually how the energy doesn't diverge. So I should point out this value here is non-zero, there is a non-zero amplitude for two opposite spin electrons to be actually exactly on top of each other. So how is it that the energy of that isn't infinite? The answer is, well, the Coulomb energy is infinite, but the kinetic energy is minus infinite, and the way this has been set up is that those two infinities exactly cancel each other and you end up with a finite number. So the point is that whenever you're dealing with finite basis sets, you don't get this behavior, you cannot resolve that. And in fact, the wave functions typically display this kind of behavior, so they round it off, and as you increase your basis set, you describe this correlation hold better and better, but it's a very, very slow convergence, so that is the problem. Slow convergence with respect to basis set, and this means one needs to use large basis sets and extrapolate to the complete basis set limit. Okay, now what about these Hamiltonian matrix elements? So the first thing to point out is that because our Hamiltonian contains at most two-body interactions, at least for electronic problems, you can easily show that if you got two slated determinants that differ by more than two spin orbitals, then the orthogonality properties of the single particle orbitals means that those Hamiltonian matrix elements are all zero, which is great. So there's some sense of locality in our Hamiltonian in this Hilbert space. So here is one slated determinant. In red is, for example, the occupied orbitals, and here's another one, so you can see I've basically emptied orbitals i, j, and I've filled orbitals a, b, spin orbitals. And so these two, I say that these two slated determinants differ in two spin orbitals. Semantically you might say they differ in four spin orbitals, but that's what I mean when I say they differ in two spin orbitals. Now what you can show is that the corresponding Hamiltonian matrix element only contains the Coulomb operator, because only the Coulomb operator is two-body, the other stuff is one-body, so you don't have to deal with it. And the corresponding matrix element turns out to be the difference between two numbers, which is the four index integrals i, j, multiplied by the Coulomb kernel r, 1, 2, multiplied by the corresponding orbitals a, b. So it's i, j, a, b minus i, j, b, a. Now these are integrals which the quantum chemistry code, the Hartree-Fock solver or whatever, has already generated for you as a byproduct of those calculations. So whenever you do your Hartree-Fock calculation, unbeknown to you, actually the code has generated a long list of these integrals. Actually you can see that there will be m to the 4 of them in general. Now I may not print them out, in which case you have to, you know, get them to print it out, and so that's most of our work in interfacing with various quantum chemistry codes has been exactly that. But once you read these integrals in, and you can more or less store them in memory, m to the 4, and then there are various symmetries, which you can use to reduce that m to the 4 by a factor of 8 roughly. Once you've read those in, you can then on the fly calculate the Hamiltonian matrix elements. Essentially by looking up, for a given pair of slated determinants, it's relatively easy to perform a line-up operation. So you line up the orbitals that match, the holes that match, and the orbitals and holes that don't match, then give you the i, j, a, b indices, and then you can just look up the relevant 4 index integrals and calculate the matrix elements. This is a bit of technology, but the main point is that these Hamiltonian matrix elements are basically order one operations. They're extremely fast. So just to give you an idea, in our code, we can generate these on the fly on the order of 10 to the 7 per second per core. So these moves, the Monte Carlo moves that we'll be doing later on, which fundamentally involve this kind of electron hole where excitation are exceedingly efficient to actually implement on a computer. Where the slated determinants differ only by one spin orbital, the corresponding matrix elements is a bit more complicated. So essentially here you take, so if they differ by one spin orbital, let's say i going to a, then you have to trace out the orbitals, the remaining n minus one orbitals where they actually match. So you have an order n operation here and where they're actually the same operator, same determinant on both sides, then you have a double loop of order n squared in the calculation of the Hamiltonian matrix elements. So the most expensive matrix elements are the diagonal matrix elements, the energy of a given determinant. So those are actually ones that we can afford to store as I'll mention later on. So that's how we try to minimize the impact of the more expensive matrix elements later on. Okay, now one of the interesting things about our slated determinants space is that it turns out to be a giant, network is giant because it's got let's say 10 to the 16 nodes on it. If you're dealing with the problem that we were dealing with earlier on, it's a giant network but it's a network which is on the one hand local in connectivity but on the other hand it has small world character and what that means is the following, is that obviously the number of nodes around grows exponentially in N and M, that's a feature of this binomial, of this byproduct of binomial coefficient. It's growing, fundamentally it's growing exponentially in the number of orbitals and electrons and the connectivity around each node of our network, each slated determinant, is roughly uniform because essentially given any determinant there are roughly N squared, M squared, double excitations away from it. You can pick two electrons to excite from and two holes to excite to and that scales as N squared, N squared. Roughly speaking to first order we have a uniform graph, so to speak, a uniform connectivity. The details of course in practice enter because there are symmetry restrictions which restrict the, which means some matrix elements are zero, etc. This is an upper bound to the connectivity of it. But N squared, M squared is still actually quite a large number. So out of each slated determinant you have a large number of connections coming out. So you can just plug in if N is 10 and M is 100 you can see what N squared, M squared is getting to probably about a million something like that. So you have a very high dimensional network in other words. In fact it's essentially close to infinite dimensional most practical purposes. Now these infinite dimensional networks are actually very interesting because it means that you can set up if you like particle flows relatively easily in Ombom. So if you have networks with very restricted dimensionality then you know things don't necessarily flow through this network very very easily. But where the connectivity is very high you actually end up with very few bottlenecks. And so that's a great thing because as we'll be doing later on where we're going to be setting up walkers on this network these walkers will be able to move around our network with great efficiency. Relatively few bottlenecks. That's not to say that you can actually set up some pathological systems with bottlenecks in them again that's a type of system where the algorithm won't work so well it gets trapped. But otherwise there are very few ergodic bottlenecks. And of course the other thing to point out is that you need only at most N over two steps to go from any one slater determinant to any other slater determinant in this network. So it's a relatively well connected network as a result. Another two is basically is the number of double excitations which you need to connect any pair of slater determinants N being the number of electrons. So in this network where we had ten to the sixteen ten to the sixteen slater determinants N is ten, N over two would be five. So in only five moves you can go from one point in your Hilbert space to any other point ten to the sixteen elements in it. So that's also a good thing. And then finally I should point out that the FCI problem is a linear problem with only one minimum in the problem. So the other stationary points are all saddle points. And that again from an optimization perspective is nice. The other many body methods such as couple cluster theory end up with non-linear problems to solve and so you end up with non-linear problems and they're obviously always going to be harder. So the idea is to combine QMC with quantum chemistry and actually this combination, so we're going to solve our full CI problem by concepts that come from QMC. And actually this mixture takes the form of a very very simple algorithm that I call the game of life. It's a bit like the game of life actually. So you have walkers that give birth that die that annihilate each other and if you play this with a simple set of rules but definite and you run it for long enough you actually converge onto the full CI solution. And what's nice is that these solutions are truly emergent solutions. So you don't start off by knowing anything about the ground state and you actually learn something about the ground state as you run the calculation which is nice. So in particular the nodal surface of the wave function you don't need to input, you don't need to know anything about the nodal surface. The nodal surface if you like emerges in the course of the simulation. So that's the full CI quantum Monte Carlo and this I which we call the initiator introduces a kind of a Darwinian twist to this to this algorithm which greatly improves the rate of convergence. It does so at the cost of an approximation but the approximation is itself systematically improved. So that's essentially the main content of the method but once you've got these emergent many electron wave functions you can then do things with them like calculate density matrices. Yeah, later on, later on I'm just talking and then we'll see specific examples later on. And then once you've got these things you can calculate your reduced density matrices which are specific operators and then from these operators you can calculate other properties and so on. So in principle you have a fairly complete method which can generate you more or less exact results for fairly sizable systems. Not as big as we would like unfortunately but we're working on it. Okay, now the key, really key step and it's a surprisingly difficult one to accept is to actually forget about the notion of these amplitudes. This CI vector which if you think about it always ties you to this gigantic Hilbert space and that more or less kills you before you've even started. So if you think about amplitudes then you have to think about approximations to the amplitudes in essence that's what we do in many body theory. But instead what we do is we think about walkers. Now walkers is a term that's used in quantum Monte Carlo. It's essentially means when you're talking about walkers it means you're talking about delta functions. So when we talk about a walker we mean delta functions. And in the context of our method because our Hilbert space is finite, it's very large it might have 10 to the 16 components to it but it's still finite, our walkers are then Kronecker deltas. Kronecker deltas are a lot easier to think about than Dirac deltas. And Dirac deltas which somehow sit in continuum spaces is what you normally use in conventional quantum Monte Carlo. Now there are a number of advantages I would say crucial advantages to working in finite Hilbert spaces despite the fact that they're very large compared to Kronecker deltas. So I'll talk about this a little bit later. In some ways it is but I'm going to come to that in a second. So yes, in some ways it is. The crucial point about walkers is that we're dealing with distributions and then those distributions describe our wave function. So in general we're going to deal with a large population of walkers. NW is the number of walkers and these populations will be on the order of millions or even billions of walkers, large populations. And actually large enough for interesting statistical mechanics in other words to arise among our walkers. So these are not going to be so our walkers will have a certain interaction between them and that's So each walker will be associated with a slated determinant. So each walker will have a slated determinant and I should point out that you can think of slated determinants as binary strings. So here is the length of the binary string. It's equal to 2m and a given slated determinant will just be a binary string. So one for an occupied orbital and zero for an unoccupied orbital. So actually from the point of view of a classical computer there's an ideal architecture for it actually. So there's an ironic sense in which fermions are ideal for classical computers. Don't tell your quantum computing friends about that by the way. They won't like it. Now the other thing is that apart from this binary string which each walker will have each walker will also have a sign minus one. So those are the two defining characteristics of our walkers. Binary string and sign. And so henceforth we're only going to be concerned with essentially a population dynamics of our walkers. But whenever we want to come back to quantum mechanics and our wave function we're then going to say that the CI amplitude on any given determinant will be proportional to the number of walkers instantaneously on that determinant. And of course this will be a strongly fluctuating entity. So this is essentially you're summing over all your walkers and if a particular walker is on that determinant you add in its sign. But the way the algorithm is going to be constructed with every iteration you'll only ever have walkers of one sign on any given determinant. And that's actually crucial part of this story. So here is a pictorial example so along this line is some description of our configuration space. So you can imagine that there are 10 to the 16 components along that line. So that is our description of our system. And I have now distributed 11 walkers on this space. 6 positive, 5 negative. So NW is 11 and the first thing to point out is that you can of course have more than one walker on a given determinant. In fact the important determinants will actually accumulate walkers on. That's actually a good thing in our algorithm. Whereas the unimportant determinants will have no walkers on them and you won't be wasting any CPU memory or effort in describing those. But they'll still be accessible. So there is no sense in which we're going to be truncating our space. The whole space is available for exploration. But at any instant in time we will have basically a pretty coarse snapshot of our system. So here you see each walker is one of these delta functions. That's basically it. Now the other thing to point out is that you can think of a vector Ni which is analogous to the Ci vector. So this vector Ni would have the same dimension as the original Ci vector. But that's not something that we actually store. And you can see that the L1 norm of that vector Ni is the number of walkers. So when we're dealing with simulations as we'll be later on, where we're dealing with a fixed number of walkers essentially what you have to think about is we're doing calculations on vectors that have a fixed L1 norm. Which is kind of an unusual thing because that's not what we do when we normally do quantum mechanics. What we normally do in quantum mechanics is work with vectors that have a fixed L2 norm which is the fixed normalization. Instead the normal normalization which is given by the sum over the squares of the number of walkers and then square rooted, this object determines the magnitude of each of our delta functions. So if you're dealing with a normalized distribution such that this object has been set to 1 then you've got a different your L1 norm will then be some number, it won't be equal to 1 of course. But that is the so you can either fix your L1 norm and then you end up with different L2 normalization or you can fix your L2 norm and you end up with different L1 normalization. But in essence what our method is going to be doing is essentially fixing the L1 norm. And so what we want to do is to interpret the imaginary time, the imaginary time Schrodinger equation which so the imaginary time Schrodinger equation is dh by do I use tau there or do I use beta which is some kind of imaginary time parameter is h psi and then we have a minus sign here. So that's the imaginary time Schrodinger equation. If you were to propagate this for large beta that would convert you onto the ground state solution. So psi of beta is e to the minus beta h acting on some initial solution and in the infinite time, some initial trial vector and as beta goes to infinity this projects out the ground state solution of the problem. So if you take this imaginary time Schrodinger equation and you substitute the full ci expansion and you end up with is this expression here. So you have dci by d beta this would if you like be the instantaneous force on a given ci coefficient if you're just dealing with a conventional ci calculation and that instantaneous force would be nothing other than minus the Hamiltonian operator acting on this should be j rather than i Hamiltonian acting on your ci vector and that would be the the imaginary time Schrodinger equation for our ci coefficient. So now what we're going to do is we're going to introduce a parameter that I call the shift into this problem and the shift acts only on the diagonal so it modifies this operator so it only acts on the diagonal basically if we're in a stationary point so the dci by d beta is 0 so we're at the stationary point there then this object is a corresponding eigenvector and the corresponding shift would be the eigenvalue and if you were to take the lowest energy one of those that would give you the ground state eigenvalue so this shift parameter that we're going to use as a population control parameter at convergence will directly give you the exact energy for the problem so in our population dynamics so that's quite an interesting point of view which of course in Monte Carlo is well known and so we're going to reinterpret this this equation in terms of a population dynamics of a set of positive and negative walkers in such a way that in the infinite beta limit the infinite time limit the expectation value of our number of walkers on any given determinant will become exactly proportional to the full ci solution of the problem and so that means that if you have a ci coefficient that has a tiny tiny value so it's an unimportant one it will only occasionally be visited by one of our walkers and as in its time average will actually acquire its exact value that's the idea so this we can sort of think of in a differential way the first thing that we do is we actually define a matrix K the original Hamiltonian matrix and we for convenience it's kind of it's nice to actually subtract out the Hartree-Fock energy from the diagonal values so that means that the all the all the diagonal matrix elements of K are either zero if they're the Hartree-Fock determiner or the positive and I've more or less said what I've just said on this slide that if you're dealing with a stationary distribution so dci by dt equals zero that implies you have a an eigenstate of K of course if you have an eigenstate of K you also have an eigenstate of the Hamiltonian and the corresponding shift would then be the exact correlation energy E minus E zero of this Hartree-Fock energy and the important point is that any initial distribution so if you start off with any initial distribution random vector set of coefficients and you simply iterate if you were to carry out this matrix vector product many many many times that would actually lead you to the exact solution ground state solution that's just a totally robust algorithm if you could implement it so the reason why you cannot implement this is of course if you're just dealing with your amplitudes then you have to store this entire vector and actually you need two copies of it one the vector one for the force and that pretty much kills you so it's really the memory bottleneck ultimately that is the real killer in many body theory ok so instead what we want to do is to think of the processes implied here by this force update in terms of a population dynamics so to do that we're going to separate out the diagonal term so k i i minus s times c i and we're going to separate this step as a diagonal death step ok so you imagine you have a set of walkers maybe one maybe more on your slater determinant on determinant i you're then going to kill walkers each walker with the probability that's going to be proportional to i minus the shift so let's for the moment set the shift parameter to zero for the sake of argument then you will kill the walkers with the probability given by k i i so that is the first step in the population dynamics rather than applying k i i to c i we're going to discretize c i in terms of walkers run through each walker and then say ok the first walker will you die with the probability k i i or not if you die you remove it if you don't die you leave it then you go to the next one you do the same thing so we want to essentially develop this notion of population dynamics instead of vector or matrix acting on vector in other words this is the crucial point so the diagonal death step in a way is the easiest one to imagine now what about this one this is saying k i j times c j so what this would be doing would be sending flux amplitude from determinant j on to determinant i at a rate given by k i j times c j so again we're going to now think of that as a kind of birth step so we sit on determinant j so if we're a walker on determinant j we now try to create a walker on determinant i with the probability given by k i j that's the idea so that's how we construct our population dynamics but you will say well hang on k i j can be positive or it can be negative how are you going to handle that well the answer is that if k i j is negative then i will construct i will construct then i will give birth to a walker of the same sign as the parent and if the Hamiltonian matrix is positive i will give birth to a walker of the opposite sign and this inversion in sign comes because we're dealing with the ground state problem so we have this minus sign so these off diagonal terms give rise to these spawning processes and then after all birth so these spawning processes, birth process and death processes have taken place there's an all important annihilation step that i'll talk about in a second so the overview of our algorithm is so i call it as a random game of life, death and annihilation we typically start off with one or maybe ten walkers, something like that on our it doesn't really matter on our Hartree-Fock determinant an initial value of the shift that we set to zero and a time step tau that's fixed again our code chooses the time step on the fly so to speak in the build up phase of the simulation so you don't have to worry about the time step either and at each step each walker tries to give birth each walker then tries to die and then after all birth and death processes have taken place we have this all important annihilation step and what that annihilation step is a pair wise removal of walkers of the opposite sign that have ended up on the same determinant so if you had let's say a positive walker here and a negative walker was spawned onto it then at the end of that iteration assuming the positive walker itself hadn't died in the meantime those two annihilate each other and you remove the problem of the simulation and this is if you like this is the crucial thermionic step in this algorithm and then finally if you're working in constant number of walker NW mode that should have a W here you then adjust the shift and go back and iterate and the shift is adjusted according to some logarithmic criterion so as a pictorial example so here is a spawning event so here is one walker attempting to spawn another walker here now there's a choice to be made because this walker could have actually spawned anywhere so we'll discuss how that choice is made but in essence the probability of spawning depends on the Hamiltonian matrix element connecting the two and the generation probability of this determinant from that one the death event is when this walker then attempts to die and then finally the annihilation event is two walkers of the opposite sign on the same determinant at the end of the iteration they annihilate each other and you've ended up with another distribution so somebody asked about the nodal structure of the wave function well essentially the nodal structure is the sign structure of this vector some determinants have positive amplitude some determinants have negative amplitude and at the end of this iteration you have a different distribution which is potentially a different nodal structure and the whole aim of FCI QMC is to just generate many many many distributions in such a way that in the long time limit you actually converge onto the exact solution okay so that sounds a bit like magic how can you do that and in fact the rules of this game are more or less unique because you can't mess with them so if you mess with them you'll end up with the wrong solution in other words and so these are the rules so and you can actually derive them from the underlying imaginary Schrodinger equation so let's first of all deal with the spawning rule so the probability to spawn so let's suppose here you're a walker on this determinant it's connected via the Hamiltonian to these ones shown dotted line but of course the connectivity in practice is very large, it's more than four but I've only indicated four for the sake of just convenience now the first thing that we do is we select another determinant let's say this one with some probability and so P gen is the generation probability of determinant J is this one from the original determinant I this one now in our method this generation probability can essentially so the can essentially be arbitrary subject to two constraints the first is that if two determinants are connected via the Hamiltonian then the generation probability has to be non-zero so you cannot if you like ignore certain links in your Hamiltonian if the Hamiltonian matrix element connecting two determinants is non-zero you must have a non-zero probability of generating that move that's the first thing the second is that whatever is the algorithm you choose to make your move you have to be able to calculate the probability of making a particular move so whatever is the algorithm that does the move must also return the number which is the probability that it actually made that move and not some other move so if you of course make the moves with uniform probability in this case that generation probability would just be a quarter and that would be relatively easy to do and actually in our early work we actually use uniform generation for simplicity but it's an inefficient way of doing it so we've and I'll talk about improvements of that later on but the crucial point is that these generation probabilities are normalized so that if I sum over for any given determinant i sum over all connected j's the generation probabilities sum to exactly one and so if you whenever you I'm sure immediately after that you're all going to go and start programming this that is one constraint that you have to ensure that your generator, your probability, your algorithm that generates your excitation is normalized so if you sum over all connections to any given determinant you end up with one okay so you make this move and you end up and your subroutine that makes this move or attempts to make this move generates you this number which is the probability of making the move then you calculate the Hamiltonian matrix element between these two so that's again easy to do you take the absolute value of that and you multiply it by a time step tau and that then gives you your spawning probability so that'll give you a number like 0.69 okay so then you generate uniform random number between 0 1 and 69% of the time you would accept that move 31% of the time in that case you would reject it if you reject it it means that spawning operation is aborted and you move on to the next walker if you accept it then you're going to actually create a walker on J and the sign of the walker that you're going to create then depends on the sign of the Hamiltonian matrix element connecting the two and the sign of the parent so if the Hamiltonian matrix element is negative the child will have the same sign as the parent and if it's positive it'll have the opposite sign so that more or less defines the spawning process now you'll note that the generation probabilities if you're doing things uniformly doesn't unfortunately scale very favorably because around each determinant you've got roughly n squared m squared connectivity so that your probability to selecting any one of them if you're doing it uniformly will go as 1 over n squared over m squared so those are probabilities that are coming crashing down so if n squared m squared is 10 to the 6 then 1 over it is 10 to the minus 6 and that is one thing the time step, tau is to somehow bound this spawning probability so if your generation probability is let's say 10 to the minus 6 and these Hamiltonian matrix elements let's say are typically on the order of 10 to the minus 2 something like 10 to the minus 3 then this product will be on the order of let's say 10 to the 4 and you don't want to be spawning 10 to the 4 walkers so that's where the time step comes in and that time step, tau has to be sufficiently small so that the spawning probability remains kind of bounded you don't want to be spawning more than one walker you can spawn more than one walker so if for example p s turned out to be 11.2 then what that would mean is you spawn 11 walkers with probability 1 and then 20% of the time you would spawn another walker as well but generally speaking we want to avoid having what we call walker blooms where one walker creates a whole mass of other walkers there because they then take time to effectively re-equilibrate so we typically choose our tau so that the spawning probability rarely exceeds 1 you can tolerate it occasionally exceeding 1 but you don't want it to grow out of hand there is actually a fundamental limit on the time step that will come to in a second and that's to do with the probability of death so the probability of death so this walker here now attempts to die and it dies with the probability that's given by the energy of the determinant minus the Hartree-Fock energy so by convention subtract out the Hartree-Fock energy there will always be a number that's either 0 or bigger than 0 and if we are in variable shift mode you then subtract out the value of the shift from it as well but let's suppose we are in sort of fixed shift mode and s is 0 so then the probability of death is basically given by the energy of the Hamiltonian so the higher so if you happen to be born on a high energy determinant then I'm sorry you will die with a high probability if on the other hand you've happened to have been born on a low energy determinant then you are likely to survive for a number of time steps so what the death step does is to if you concentrate the walkers on the low energy determinants of your Hilbert space whereas the spawning process is more or less quite agnostic as to where you end up spawning so you typically end up creating walkers on high energy determinants and they die relatively quickly those that are spawned onto low energy determinants live long enough to be able to themselves have children and propagate the process so that is roughly speaking this game and now you can see that when we're going to variable shift mode the shift which will equal the correlation energy the correlation energy is always negative so as you go into variable shift mode you enhance the probability of death in other words that is in essence what's going on and that's how the system does its population control so typically you have an exponential process, birth processes that are giving rise to increasing numbers of walkers then you're getting how you control this explosion is typically through the death step and that reduces makes s negative which enhances the probability of death and this balance between the two fixes more or less what the probability of death is is the potential process which is the annihilation process which also limits the particle growth as well okay so that's the description in plain vanilla form of the algorithm and now if you want to calculate the energy which is like the most basic quantity we want to do we typically do it called a projection formula now you're probably used to in quantum mechanics dealing with expectation values of the type psi operator psi that's great but they're very difficult operators to handle, very difficult expectation value to handle so instead we deal with this projected thing where you see psi the wave function that is being sampled is let's say only in the set not in the bra and in the bra instead is a trial wave function in this case the Hartree-Fock determinant d0 and that d0 h psi divided by d0 psi so that's the overlap between the exact wave function and your trial wave function on the denominator now obviously if psi were the exact solution to the problem so h psi would give you the energy times psi then this would in fact give you the exact energy because in that case the numerator and the denominator would exactly cancel leaving you just the energy so that's why this is a useful object if psi were exact the corresponding energy would be exact but it's much easier to handle now you can actually work it out that what you end up with is the Hartree-Fock energy so that's when you sum over j j is actually equal to the first term where j is equal to the Hartree-Fock determinant d0 plus a set of corrections that runs over the double excitations of d0 in principle also the single excitations but in Hartree-Fock theory the singly excited determinants have a zero Hamiltonian matrix element with the Hartree-Fock determinant so their contribution drops out but you could easily include the singles in here as well so you have d0 hdj so that's the Hamiltonian matrix element connecting the Hartree-Fock with the double excitations times cj over c0 and that gives you your the contribution of the corresponding double now this might to the uninitiated you might think well if that's true then all you need is to do a calculation truncated at your doubles and that would give you your exact energy but that's actually not the case because these cj's the amplitude at the double is itself connected to the rest of the space to the triples and the quadruples and the triples and the quadruples are connected to the fifth and sixth order slated determinants and so on so you have a problem in which all the determinants get involved in determining the amplitudes of the double determinant but only the ones of the doubles enter this expression for the energy so you have this ratio cj over c0 and now we replace that ratio with the number of walkers on the respective determinants so again that's in the plain vanilla form of this algorithm now you can see the problem immediately is that you got this n0 here and you got the walkers on the Hartree-Fock determinant explicitly entering this expression and what does that mean is that first of all that number better not be 0 because you won't get anything out of this and the other is that you want that number to be quite large so that fluctuations in that number are not too severe so the fluctuations in the denominator are particularly damaging so in other words whenever you have been able to accumulate enough walkers on the reference determinant and then everything will work fine ok so let's actually just see a simulation in practice how much time by the way do we have alright ok then we'll have coffee now and then we'll... yeah no problem