 Okay, so what we're going to have a look at is just this algorithm in its simplest form implemented for this molecule, BE2, in this basis set. You don't have to worry about it. The only thing that's relevant is that the Hilbert space is about 346,000 determinants. 346,485 determinants, to be exact. And so that is a non-trivial problem. You can exactly diagonalize it. And in blue is the exact ground state eigenvector for this problem. So along this axis here is some discretization. It's a discretization. So there are, in fact, 346,485 impulses on this plot. The vast majority is so small that they appear as zero, but they're actually not zero. If you look at the CI vector, there's plenty of 10 to the minus 8s, 10 to the minus 10s, 10 to the minus 12s, and so on. But the ones that are not substantially norm zero you can see in blue. This one out here at an amplitude of minus 0.8 in a bit is the Hartree-Fock determinant. And then you have these ones out here and these ones out here. And the main point is that I really don't care what these guys are, where they are, or any such thing. And that's one of the nice things about, so I want to start off with minimal information regarding the ground state solution. So what we're going to do is to start one walker that's in red on the Hartree-Fock determinant and then just execute the algorithm. So let's see what happens. As we, let's just concentrate on this top plot as this simulation proceeds. So the red now is our instantaneous distribution of walkers. And as time goes on, and this is the number of iterations on here. So what you begin to see is the distribution in red actually settles down on the original exact eigenvector. And in particular you can see the signal that's growing on the Hartree-Fock determinant. And that indicates more and more and more walkers are actually condensing onto it. And that's actually what you want. So if we, I'll just rerun the simulation from the start again. Now on this plot we see the energy. And you can see the projected energy shown in green is fluctuating all over the place. And at some point it begins to settle down. And actually at that point if I just go back, if we just watch it, you begin to look for a red signal which you begin to see, you begin to see it about now. You see that red signal is growing here. That is when the projected energy begins to dampen down these wild oscillations. And you can see it settles down on the red line here, which is the exact solution to the problem. So that's what I was mentioning before is that until you end up with a large number of walkers on your reference determinant, your projected energy is showing very large fluctuations. And in fact, if the sign on the reference determinant is oscillating sometimes as positive, sometimes as negative, you even go through essentially, you get these singularities there. They're completely meaningless. But once you get enough walkers on your reference determinant, then you get this nice signal emerging. And the energy essentially settles down on the exact solution. Now what's shown in blue is the shift, which initially is held zero. So in this phase of this simulation, we're not trying to control the population of walkers. In fact, the population shows a very interesting dynamics. What happens is that's shown down here is first of all, the population grows. And then despite the fact that we're not controlling the population, the population stabilizes. Here it goes. It's completely flat. And then it begins to grow again. Now this feature of it going flat, as we call this the annihilation plateau, that's actually the point at which you have enough walkers that they start bumping into each other and start annihilating each other. And at that point, the wave function actually properly emerges. And as it comes out of this plateau phase, you get another exponential rise. But although it's not obvious in this plot, it's a lot slower than this. Now when I hit the target number of walkers, which in this case was set at 2 million, which was just the input, we'll vary it in a second. I now vary the shift in order to control the population at 2 million. And the shift now relaxes and oscillates about the exact energy. So I have two measures of the energy in other words. I have one coming from the population control and I have one coming from the projected energy. And if those two don't agree with each other, then we know that we don't have the exact solution for the problem. Now what's shown here in this plot is the amount of annihilation that's actually going on in this eigenvector. And the interesting thing to see is that these regions where the exact solution is sparse, is that sparsity arises to a large extent because of annihilation. In other words, you get positive and negative walkers being spawned onto those determinants, but they're being spawned at roughly the same rate. So they annihilate each other and nothing is left. And that's what makes these fermionic eigenvectors so difficult to intuit. So many of these determinants out here are very strongly coupled to the important ones. But they're just unimportant in the final solution because of annihilation. So that was the beryllium dimer. Here's another example with a water molecule. This is now a much bigger Hilbert space, 252 million determinants. This is again exactly diagonalizable but now with considerable effort. And here you see we start our simulation off with I think with one walker and it grows very rapidly, in this case to 26 million walkers. When it hits 26 million walkers, once again the simulation goes into this annihilation plateau. The growth just automatically attenuates for quite a long period of time. You can see maybe 15,000-20,000 iterations. And then once again it starts to come out of this annihilation plateau. But this rate of exponential growth is now what I call coherent exponential growth. It's a lot slower than the initial growth phase that was essentially incoherent exponential growth where you were just spawning walkers but they weren't annihilating each other. You can see in the dotted line is the amount of annihilation that's going on. You see initially until you hit the plateau there is not much annihilation and then you get to a critical amount of annihilation and once you get to that critical amount of annihilation when you combine it with the normal death rate. So you see this is the death rate, this is the annihilation rate and that's the birth rate and if you add these two together you get exactly that and you end up with no growth whatsoever. So in this case the annihilation plateau is at 26 million walkers and then so here I think we'd set our target number of walkers to be about 50 million or so or 40 million which is about one-tenth of the total Hilbert space. So what we have as an interesting result is that the amount of storage I need to get the exact solution is about one-tenth of the original Hilbert space which is actually quite an interesting, I mean if I say that to a numerical analyst they will say well how can you possibly do that and that is what the Monte Carlo method is, this sampling technique is buying us. So I was actually very excited when we first thought of this algorithm because we had this nice population dynamics that could give us the ground state solution of these huge Hilbert spaces and so the question was well what's the snag? And the snag is that the annihilation plateau, the number of walkers that you need to hit before you can extract anything meaningful unfortunately scales with the original Hilbert space. The pre-factor will be different from one system to another. So for example this is actually a selection, this is more or less all known exact results in molecular quantum chemistry and so we looked at all of them. We'll actually just do this one in a second as an online demonstration is the neon atom. Here you have a plateau, Hilbert space of 6.6, 6.7 million determinants and the number of walkers that we need to get to the plateau is 210,000. So if you like the fraction of the space that you need to instantaneously occupy with walkers is about 3% of the space. If you now deal with the carbon dimer in this basis set that has a Hilbert space of about 27-28 million determinants and that needed 15 million walkers before you could hit the plateau. If I terminated my simulation or if I went into variable shift mode before I hit that 15 million then I get nonsense out of the simulation. Why do I get nonsense is because I haven't condensed enough walkers on the Hartree-Fock determinant to be able to get a signal and so that means that I needed to occupy about 54% of the space before I could converge and so on. Nitrogen molecule, typically for these molecular systems you need about 40, 50, 60% of the space instantaneously occupied by walkers before you can get to this plateau. But nevertheless we were still able to get some system that would be more or less beyond the exact diagonalization, for example the oxygen molecule or the CO molecule. These were Hilbert spaces that were now really very big. This in this case is about 5.4 billion determinants and for the oxygen molecule we needed about 2.6 billion walkers before we could get to the plateau. So the big question was how, maybe I should do a demonstration at this station. It's kind of interesting. I will just log on to my, just to give you an idea of timings and things, let's look at the neon atom and just, so I have two files here, FCI dump, this is just a symbolic link to it, which is the file that my Hartree-Fock Solver has provided me that is this list of integrals in it. So if you look at it, you have some header information. So this is a problem with 23 spatial orbitals, so 46 spin orbitals and 10 electrons. There's some information about the symmetry and there is a long list here of these four index integrals. So 1, 1, 1, 1 refers to the four index integral, you know, spatial orbital 1 for electron 1, spatial orbital 1 for electron 2 and so on and these are the, these are the integral. So this is a file which I'm going to assume is the input to our calculation. So if you've got a molecule and you want to say, okay, I want to do an FCI QMC calculation on it, then the first thing you need to do is to take mole pro for example, run it, get the Hartree-Fock orbitals and get this file out of mole pro, which is possible to do. And now we want to run our calculation. So there is a second input file which specifies the input parameters and the only thing that's going to, that we're going to be concerned about here is the total number of walkers that I will set to 500,000. Now I'm setting it to 500,000 because I know in this problem the plateau is actually at about 200,000 or so, but this is just as a demonstration of what's going to go on. And so now let's just start the calculation and I'm going to start it with one walker on the Hartree-Fock. Okay, so we will run, so this is a quad core, so we will run it on the four processes that we have. So now we've launched a four processor job. Our code is called NISI. It stands for N-Electron-CI. It was the original CI code that I started developing. Okay, so let's now have a look. The first thing that we're going to look at is the number of walkers, which started off at one. And just in the few seconds that we've been talking, it's got to actually, let me just have a look. I think this is an initiating mode. Yes, I think I am. I'm going to stop this job and rerun it. This wasn't quite the calculation I wanted to do. I want to display the, let's rerun it. Okay, nice. So just in these few seconds, you see we started off with one walker and we've actually got to about 150,000 or so, so this problem is the plateau. Now it's a bit lower than the plateau that I was showing before, which was at 210,000. And the reason is that we're using a certain spin symmetry in this calculation that we weren't in the previous one, and that has the effect of lowering the plateau. But basically, if you plot this on a log scale, you see you get super rapid growth and now the population is stabilized at the plateau. Let's just update it. You can see now it's slow growth. Now what about the number of walkers on the reference determinant? That's now shown in green. Now you can see that until I hit the plateau, I didn't have many walkers on the reference determinant. One, two, but then when we go into the plateau, so this is now the population in the whole space is mysteriously not growing very quickly, but the population on my reference determinant is now racing up. That's kind of showing that you're getting an accumulation of walkers onto the reference determinant. Now if we plot the energy, you can see wild oscillations to begin with, and then it settles down on some very well-defined value. You can now plot, if we now put the exact energy on, which is actually nothing, that might be the nitrogen molecule. Let me just find out what the exact energy is. By the exact energy, I mean if you did an exact diagonalization, so you have this exact energy to more than six figure accuracy, with any luck they will be rather close to each other, so you can see this is our oscillations converging onto the exact result. Okay, now suppose we were to do the calculation, but we didn't allow the system to grow to the plateau. So I'm now going to reduce the number of walkers to 50,000, or let's say, I've never done this calculation by the way, so I'm quite curious to see what will happen, but let's now repeat the calculation. So now we will go into variable shift mode at 50,000 walkers, and that happens to be below the plateau. So here you can see the number of walkers shot up, we went into variable shift mode, and actually we've now dropped below our 50,000 that we had set, but this is the number of walkers that's on the Hartree-Fock determinant, let's put it on a log scale, and it's essentially fluctuating between one, it probably gets to zero, and then it doesn't appear on this plot because there's a minus infinity and so on. So you can see that if I go into variable shift mode before we get to the plateau, we don't get a signal appearing on the Hartree-Fock determinant. And if I then try to measure the energy, I'll probably get total nonsense. Let's just have a look. I mean, look at this scale. Ironically, it's actually, the mean may not be so bad, but I wouldn't rely much on that. I mean, the scale is absolutely huge, and so this is the problem that basically I don't have any reasonable description of my wave function. So the message is that if you're running with FCI QMC, let's go back to our talk now, you've got to hit the plateau, but of course the plateau may be at such a high number that you may not be able to get there. And then the problem is, well, how can you... It's not, actually. So it's usually, if there's any determinant that has a significant weight, it's usually the Hartree-Fock, but of course there are many problems in which there are other determinants that have comparable weight to the Hartree-Fock. And that's one of the things that we'll deal with later on. But anyway, the question is, how can we reduce the number of walkers while maintaining full CI accuracy? In other words, we want to be able to run walkers with, you know, reasonable number of walkers and still end up with the exact solution. And this actually turns out to involve one minor tweak to the algorithm that we call survival of the fittest. And this is the initiator step. And it occurs immediately after we've tried to spawn. And so this was actually... So the way we came up with this idea, I had a graduate student, Deirdre Cleveland, who really traced, if you like, trajectories of walkers over time. And then we realized that the damaging problem was whenever determinants... So occasionally it's, let's say, occupied by a positive walker and then it dies and at some time it becomes occupied by negative walkers. So the signal emanating from that determinant would change sign. Because of course, whenever a determinant changes sign, all the subsequent progeny also change sign. And that becomes, if you like, a source of incoherent noise. And we wanted to cut that out. And the question was, how do you cut it out? And so what we say is... So we introduce two if statements into the algorithm. And so it introduces the concept of an initiator. So let's suppose here we're a parent determinant and now we want to spawn on this determinant. Suppose we've selected that to spawn on. And we spawn a child here. So we've already decided we want to spawn a child here. Now we ask two if statements. So the first if statement asks, is D currently empty or not? In other words, is there currently a walker on D or not? If there is a walker on D, then the algorithm proceeds as normal. We don't change anything. But if D is empty, so in other words, we are just about to bring to life a new determinant, we then ask how many walkers were on the parent. And we only allow that new child to survive if the parent determinant had a number of walkers that exceeded a critical number, typically three. Otherwise you kill that new child. So in other words, determinants that instantaneously have three or more walkers on them are allowed to bring to life new determinants. Otherwise those children are killed. Whereas if you're not an initiator, so your amplitude is less than let's say three, then you can continue to spawn on the already occupied determinants. But you can't spawn on new ones. And actually I'll do an online demonstration, so I'll skip that thing. So the first thing that arises is the question, well okay, that's a nice rule you're introducing, but surely you're going to mess up your solution. And the interesting point here is that our initiator method remains exact. And it remains exact in the limit of large walker limits. So I have a guarantee that as I crank up my number of walkers, my algorithm will revert to being the original algorithm. The original algorithm we know is exact if you're above the plateau. So therefore we have this guarantee that in the large walker limit the initiator algorithm remains exact. And so what's the proof of that? The proof is that in the large walker limit all determinants acquire an occupation. And so therefore even if you're a non-initiator, the children that you're spawning will survive on those determinants. And in other words the first if statement, the survivor of the fittest if statement just doesn't arise. You pass the first if statement. And so therefore the large walker limit of our initiator algorithm becomes the original algorithm. And because we know the large walker limit of the original algorithm is full CI, we have that the large walker limit of the initiator algorithm is the full CI limit. Now what's interesting is that whereas the original algorithm before you get to the plateau you cannot extract any noise, any meaningful information out of it. In this problem it would the initiator method even with a very small number of walkers you actually get a very, very precise solution. So this is an example, we can run this in a second of the nitrogen molecule with 540 million determinants in it. So this plateau is at about 240 million walkers. So with the original algorithm if I truncated the calculation before I got to the 240 million walker limit I would get nonsense out of it. But instead here only with 10,000 walkers first of all I get a very reasonable energy, it's only a couple of millihartries off the exact one. And secondly the arrow bar is nice and bounded. And then as I increase the number of walkers here at 100,000 I'm really pretty close and then by the time I get to a million walkers I've essentially, I'm indistinguishable from the exact solution. So let's now run this neon calculation but now let's do the initiator method. So now I'm going to switch on the initiator keyword. It's frozen. And now we're going to keep the same 500,000, sorry 50,000 simulation that with the full method gave us nonsense. So if we now plot the number of walkers so now you see we've gone, we've hit 50,000 we've now gone into variable shift mode so we're stabilizing. And now you see here that's already the population on the Hartree-Fock determinant. So that is, if we go on to a log scale probably about 5,000 walkers. So you see with the initiator method we actually grow the population on this reference determinant. It essentially a similar rate to which we are growing the full population. Let's see. Network seems to have slowed down. It's frozen. Let's see if I could fire it up again. OK, well let's continue our talk. We'll come back to this in a second. OK, so one thing that one might ask is is this algorithm sensitive to this critical parameter? So we set it to pick three as defining the initiator threshold. But in fact the algorithm is remarkably insensitive to that parameter as well. So you can set it to be one or two or three or four or ten and in all cases you can see convergence in the large walker limit to the exact value. So that's a good thing. You don't want algorithms that show excessive sensitivity to some of these internal parameters. Let me see if the connection has come back. Well, let's hope it comes back. OK, so how does this initiator method actually work? So this is going back to the beryllium dimer again. We do now a simulation with 2,000 walkers and in this problem which the Hilbert space was about 350,000 walkers the plateau was if I recall about 200,000 walkers in that region. So now we're doing a simulation with many, many fewer walkers than the value of the plateau. So here we start off the simulation I think with ten walkers and we grow the population to 2,000 and then we stabilize the population. And in green you have the number of determinants that are instantaneously occupied. So obviously it's below 2,000 it's I don't know 1,900 in that order 1,800 in that order. And of those instantaneously occupied determinants in red it's shown the number of initiators so the number of determinants that actually have a population of more than three walkers. And that's shown in here. So there are about 80 or so initiators instantaneously. And the first thing to point out is that you can see this is a fluctuating number. So you can be an initiator at this time step because you've got four walkers and if on the next time step you drop below four walkers to three then you cease to be an initiator. So the initiator, the distribution of initiators is itself a strongly fluctuating object. Okay, now the first thing that might appear what this algorithm is doing is just selecting the important determinants in the system. This might be a reasonable hypothesis as to how this method is working. But that's actually not what is happening. So if you take those initiated determinants and then you just diagonalize the Hamiltonian in that space of let's say 80 determinants. So you set up your 80 by 80 Hamiltonian and diagonalize it and get the lowest energy eigenstate. You find that the energy of that lowest energy eigenstate is way off the exact one. So the exact correlation energy in this problem is about 150 milli heart trees. Whereas the correlation energy that is embodied only in these 80 or so determinants is only one third of that. So that's one thing. Now you can say, well what happens if I take all 1,800 or so determinants and I diagonalize in the space of 1,800 determinants that I'm instantaneously occupying? Well, of course it's a variational calculation so therefore you'll gain a bit more correlation energy that's shown in purple. But it's still only half the correlation energy of the problem. And yet if you look at the projected energy that's shown in green, you can see that it's nicely oscillating about the exact value. So in other words, the projective dynamics that we are doing is actually sampling the correct wave function although the best instantaneous wave function is completely off. So that can be kind of shown here. So in blue is the exact wave function and now we've actually reordered the determinants so that the first set is the Hartree Fock and then the double excitation and so on. But it doesn't matter. So in blue I show the exact amplitudes and in green is an instantaneous snapshot that we get from our Monte Carlo. So this, I don't know if you can see it, this thing here which has got an amplitude of in these units two is a determiner with two walkers on it. This one here has got minus one walkers on it and so on. And so instantaneously you have your walkers and they're kind of coarse grain description of our wave function but nobody would say that our green is in any sense an accurate approximation of our blue. So instantaneously our sampled wave function isn't very good. We've only got 2,000 walkers. But now if we do something that we don't normally do in our simulation but it's to histogram the space. So we introduce a vector that is now the length of the Hilbert space and every time a walker enters it for the time that it's living there we keep a running average and so we keep a running average the time average our distribution. So that as I say is something we don't do usually because it would essentially obviate the whole method but then if you actually time average our populations this is what you get. It's a staggeringly good representation of the exact wave function. So you can see the fine details of the original wave function have now been resolved. And so I come back to this that the way this method works is by generating many distributions of walkers which on average give you the exact solution not instantaneously. Okay let's see if we can get okay ah we got it. So let's plot. So this is the calculation that we did with 50,000 walkers and there you go. So you see this is the exact solution and now this is in our initiator method and you could continue this simulation for as long as you like and you will get a completely stable signal oscillating about the exact number. So that's with 50,000 walkers. Let's now do a calculation with 500 walkers. So I'm now going to do a calculation with 500 walkers. Let's see what happens now. Again I haven't done this calculation before. So now this is a really gross undersampling of this Hilbert space. By the way the Hilbert space 6.6 million determinants 6.7 million determinants. So now we've only got 500 walkers in the whole space about 100 walkers on the on the Hartree-Fock determinant. Again it's stable. This is the crucial thing. It's showing fluctuations but it's stable. And if I now plot the energy let's see where the energy is, there you go. Again we're getting oscillations now quite large not surprising. So this by the way is about 50 milli-Hartree. So we've got oscillations about 50 milli-Hartree but stably about about the solution. We're not getting so the usual problem with fermion Monte Carlo is the longer you oscillate the wilder become the the long time fluctuations until you cannot extract any signal from this. So that is the and actually you can take this down to really tiny numbers of walkers and as long as the Hartree-Fock determinant has a well-defined sign you find that you get a well-defined solution. So that is the initiator method and now you can ask well okay well that's very well enough but what ultimately is the scaling of the algorithm? So this is actually quite an old calculation that we did on a set of atoms and anions and anions are really very strongly correlated beasts they're difficult to describe and so we did these calculations on I think basically from lithium first row and then to sodium and along this axis is the size of the Hilbert space that we actually span for those set from roughly a thousand to ten to the sixteen in this case and on this axis we show the number of walkers that we needed to converge the calculation and our rule of thumb when we do our calculations is we put in enough walkers so that we get 50,000 walkers on the Hartree-Fock determinant that's our rule of thumb and so why 50,000 and the answer is that you can show that the fluctuation in the number of walkers basically follows a square root law so that if you've got 100 walkers on any given determinant then the fluctuations on that determinant will be on the order of 10 50,000 walkers then the fluctuations on that determinant will be on the order of 100 which is 1% relative and 50,000 you're below the 1% you're below the 1% limit and so this is the number of walkers you need in order to end up with 50,000 walkers on the Hartree-Fock determinant and the point is that you can see that if you go up from about 10 to the 3 to 10 to the 16 determinants actually the number of walkers that we needed goes up from about 50,000 well that was the minimum we needed because that was the criterion that we said through to about 10 million walkers so with only 10 million walkers in this case we were able to very very accurately resolve these systems that had 10 to the 16 slated determinants in them and so you can then actually fit a straight line through this log-log plot and you get this exponent so the slope of this line turns out to be about 0.16 so in other words the method remains exponential scaling but because it goes as NSCI is self-exponential but to a small exponent these points here by the way refer to the original what we call the full method without initiators and actually the slopes here are as it has a slope 1 so the original algorithm has the same scaling as exact diagonalization but the initiator method reduces the exponent of the exponential scaling and we don't really have a theory for this 0.16 by the way turns out to be fairly uniform although there are difficult systems in which that number actually goes up so if you're dealing with a Hubbard model for example in a difficult u-regime the exponent is unfortunately we think close to 1 but for many molecular systems at least it turns out to be quite small so the thing is that the set of initiators is dynamic and self-selecting so from the perspective of the user you don't have to make any choice of what's going to be an initiator and what's not and that's a very nice feature of the algorithm it basically does the selection itself and so the method really is a black box method no knowledge of the wave function is assumed and the key parameter is the number of walkers so you have to go on increasing the number of walkers until your energy basically ceases to vary that's the thing and so we usually aim for enough walkers so that you've accumulated about 50,000 on the Hartree-Folk determinant there's a fairly conservative criterion and you usually get to convergence well before that although there will again be systems where this is actually not sufficient at the same time we're also aware of those systems but for most systems if you can run enough walkers so that about 50,000 end up on the Hartree-Folk determinant then you're okay now the CPU time is proportional to the walker number and that's actually a non-trivial statement because if you implement the algorithm in a naive way you would think the annihilation step actually would scale quadratically with the number of walkers so for every given walker you have to search through all the other walkers to find out if there's anything that's co-incident on it and so that's the naive implementation of annihilation and that would obviously scale quadratically it turns out that there are very nice either using sorting algorithms or better still hash algorithms that we could actually use in practice now in which you can in order n steps instead of order n squared steps determine how many walkers are co-incident on any given determinant and so that makes the CPU time proportional to walker number on the downside so what's the downside of the algorithm and the main downside is that I cannot at priori tell you how long it would take what CPU resources you would need to solve a specific problem so if you do normal many body theory apart from the number of iterations you might need to converge a problem you can work out more or less beforehand how much time you're going to need to solve a problem here you cannot do it because I cannot tell you how many walkers you would need in practice to solve a problem so this means that there is some amount of experimentation that's involved but in practice I mean we do calculations at 200 million walks actually that's an old slide we've gone to well beyond that now even a billion walkers and that's perfectly feasible you need more than a quad core to do it but you can do it on a machine actually ironically the biggest bottleneck that we have with our calculation right now with our method right now is the input output of actually dumping the information of our distribution of walkers on to disk and then reading it back in for a subsequent calculation that can take hours of time actually when you're generating gigabytes worth of I have on a parallel machine okay it's nice and parallel algorithm although with some of the tweaks that I'll be talking about the more sophisticated versions don't parallelize quite as well but nevertheless we can typically parallelize quite successfully to several thousand cores before you begin to move off to move off the ideal scaling so that's a nice feature again with most sophisticated quantum chemical algorithms couple cluster and so on is very hard to get parallelism out at a thousand or several thousand processor level okay so let me see how far we've got eleven thirty how much time do we have half an hour or so okay so I'm going to now talk a little bit about some more technicalities of the algorithm and the first thing is that what we actually store when we're talking about a determinant is we store binary strings so so given determinant you can think of as a list of ones and zeros determining what is occupied and what's unoccupied and the nice thing about this binary representation is that of course there's sort of a unique and invertible map that takes you from your binary string onto a given integer so you have a list of these NIs which is ones and there is an integer associated with that so that is if you like a label or given this capital I label or index I can work out by inverting the map which orbitals are occupied and unoccupied and that's that's something that can be done very efficiently so this at the heart of our algorithm is this is going backwards and forwards between this so when I have a list of determinants I have a list of these I's but any given capital I I can easily decode and get the occupied orbitals in that determinant so the first thing is that this gives us a compact storage so if I have a system with 2M spin orbitals then the number of integers I need let's say 64 bit integers I need to encode the information is this is the ceiling of 2M over 64 so if I have a problem with 200 spin orbitals then essentially 4 integers 464 bit integers is sufficient for me to encode the information about any determinant I like irrespective of the number of electrons in the system so that's compact storage which is one thing now the other thing is that a random list a random selection of determinants can always be uniquely ordered and in fact in the version before the one that we actually currently are on we would always maintain an ordered list of determinants so every time I spawned a new walker and I wanted to add it to the list of currently occupied walkers I would scan the list and I would insert that determinant in the appropriate point so there is a bit of shuffling of memory but the crucial point about keeping these ordered lists is that it makes a subsequent searching exceedingly efficient because you can search ordered lists in log time and this is crucial for the annihilation steps so every time I spawn a walker and I don't know if that determinant is occupied or not I can easily search within our list of gigantic list of I don't know 100 million walkers I can easily find if there is a walker there as I said nowadays we actually use the hash algorithms to do that so it means it basically obviates a binary search of a long list but you know those are their fairly analogous algorithms so to further uses of our bit string representation is that our excitation generation is extremely efficient so given a bit string I know which bits are ones and which bits are zeros so when I come to do an excitation I can select a pair of electrons and a pair of holes with extreme rapidity in fact we can do this essentially this process is a fundamental process in about 10 to the 8 of these per second per core that's how rapid this process of excitation generation is and so that's how we know we can do these super fast calculations involving many many many walkers but the other is parallelization where we used hash functions before so essentially given the index of our determinant which can be several integers by the way but anyway that's a four or five integers you put it into a hash algorithm and what the hash algorithm does is basically return you a random number distributed between zero and one but it's a deterministic random number so for the same input I always get the same output which is what you need so then I can take this modulo the number of processors I've got to distribute our walkers among the processors so that's so in a parallelization scheme that's actually very important that every time a walker is spawned on a given determinant it's deterministically sent always to the same processor and then the annihilation is locally done on that processor so what you don't want to have two walkers on the same determinant but on different processors because that would then make the annihilation very difficult yeah so what we store is a ordered list of instantaneously occupied determinants together with a number of signed walkers on each determinant so now here we kind of move a little bit away from a strict walker representation so for each determinant we actually store the number of walkers on that determinant rather than just having a distributed list only of walkers so we here we have an ordered list and we have over here N of I1 is the number of walkers signed number of walkers on determinant I1 N of I2 signed number of walkers on I2 and so on so that is if you like what is instantaneously stored in our in our main memory list and so what happens is that we run over our list of parent determinants and for each parent determinant we loop over the number of walkers on the parent determinant and it attempts to and so each one of these typically makes multiple attempts at spawning walkers multiple if there are multiple walkers on it and we get a new list of newly spawned walkers which again once again we sort and then for each one of these so we now loop over this J list and we update our I list so to speak with the number of walkers that we have that we have introduced through the spawning event and so this basically takes care of the annihilation step so for example if there was a positive walker here so this was one and this happened to be minus one then then those two would annihilate each other and this would become zero and it's crucial that whenever a given determinant acquires zero occupancy it's then deleted from the list of instantaneously occupied determinants that's the crucial point of this algorithm because otherwise you end up with longer and longer longer lists with no occupancy in them so that this step if you like keeps the number of keeps the number of instantaneously occupied determinants small now if it happens that one of these J's is not present in this list then that means that we are about to introduce a new determinant into our list and then we would locate its position and insert it in there subject to the initiator criterion okay so merging an annihilation basically takes ns log nd steps so number of spawns and this is log nd is the number of occupied determinants that are instantaneously held so this process actually only takes a few percent of the overall cpu time it's very important to so because again naively speaking one would think that the annihilation it becomes a bottleneck in the algorithm but in practice it's a really minor component of the overall cost okay so let me just show you some examples these are from electron affinities so you calculate the energy of the atom and you calculate the energy of the anion you subtract the two and that gives you the electron affinity and so here you see we're doing it as a function of these basis sets augmented basis sets going from double zeta through to quintuple zeta so this is basically this just shows you here the size of the basis set so x is the cardinal number it basically tells you the highest angular momentum component and it basically goes as x cubed if I'm not mistaken if you multiply this out you get x cubed so when we talk about augmented vdz we're talking about 23 orbitals per atom first row when we're dealing with augmented vtz we're dealing with 46 augmented vqz we're dealing with 80 and so on so what's shown here is the basis set so we're really going to very large basis sets and the corresponding electron affinities and the good thing about electron affinities is that then they'd be measured to very very high accuracy so this is experiment where there's no doubt about the experiment in other words and also the non relativistic limit is probably good for these for these atoms but the point to see is that we converge from below the experimental value is here and this dotted line on either side is what the chemists call chemical accuracy is 1k cal per mole so that's the ideal accuracy you try to aim for you try to aim for you will try to be within these within these brackets and so you can see convergences towards experiment it's very interesting sodium for example sodium anion is not a beast that you would normally come with because you'd normally ionize a sodium rather than put an electron in sodium so you have this outer electron which is very loosely bound but you can see the extraordinary accuracy with which you can actually calculate this the biggest calculation in this plot here is the fluorine anion which has got 10 electrons in this V5Z which has got 127 orbitals so that's got a Hilbert space of 4 by 10 to the 15 and we needed about 8 million walkers to sample that problem in other words to get 50,000 walkers on the on the hot spot and so that is instantaneously if you like the fraction of the Hilbert space that's being instantaneously occupied is about a few in a billion so that's how sparse is the is the vector that we're actually sampling but it gives you a stable solution so this is another example now these are more complicated C2 molecule up to F2 molecule in these basis sets so these have got 108 spatial orbitals or 216 spin orbitals running from about 8 electrons to 14 electrons and the Hilbert space has now run from about 10 to the 11 through to 10 to the 19 so and these are the numbers of walkers you actually need to converge these calculations they run from a couple of million to about 100 million also in the case of NO now there's really interesting things about this one of them you can see for example the CO molecule we needed about 60 million walkers to converge and the fluorine molecule we needed about 50 million walkers to converge but the F2 molecule has got the Hilbert space has 5 orders of magnitude larger than the CO molecule 10 to the 19 as opposed to 10 to the 14 so why do we need more walkers in the case of the CO molecule and that's because intrinsically the wave function of the CO molecule turns out to be more complicated more multi-referenced less dominated by the Hartree-Fock determinant compared to the F2 molecule and that's kind of interesting that the algorithm senses the complexity of the wave function and that tells you okay I'm going to need more walkers to get to convergence and that's the kind of sort of intangible thing that you or I or essentially no one else has enough insight into these wave functions to be able to really talk about the complexity of them but the algorithm itself is sensitive to that and that's a very nice thing about it I don't know if we want to go through numbers but this plot is kind of interesting it shows a similar thing similar trend to the electron affinities again these are the dissociation energies so the energy that you need to fully dissociate a molecule so to do these calculations you do one calculation of the molecule at its equilibrium geometry and then one or possibly two calculations of the isolated atoms if there are these heteronuclear things and then you subtract the molecule from the sum of the two atoms and again as you increase your basis set you can see you approach the experimental numbers but now you can see chemical accuracy which is one k-cal per mole is on a finer scale compared to the results that we're getting so even when you're dealing with these quadruple zeta calculations so in the case of CO in principle that's a wave function that's been expanded at 10 to the 14 functions but that's still a couple of k-cals away from experimental value in other words you haven't hit the land of milk and honey which is chemical accuracy so it's sort of telling you just how difficult it is to achieve chemical accuracy for these systems but the slow convergence of these results is actually not so much a correlation problem at the full CI level but is a representation of that Coulomb cusp problem that you need to get to larger and larger basis sets to be able to describe the Coulomb cusp so typically in this cusp in this region from a chemical perspective there may be a couple of k-cals of correlation energy differential correlation energy buried there in other words if you do the same calculation on the atom and the molecule and the k-cals per mole because of this not capturing this region this short distance region and that is the error that is being seen here which I'm afraid is pretty depressing news because you know you do a huge amount of work and you're sure that you're going to hit this land of milk and honey and actually you don't you're actually short of it which we were within it but anyway that's a small system so what are we going to do in order to in order to accelerate basis set convergence and there are basically two techniques and they're both approximate and they involve either using extrapolations from earlier basis sets so that you can show that the correlation energy typically scales or the error in the correlation energy scales as 1 over x cubed third power of the basis set and that's a pretty general result actually it's essentially to do with the way you're filling space so when we're doing these TQ extrapolations we're essentially taking the correlation energies that you get from the triple zeta and from the quadruple zeta and then assuming this kind of an x cubed law extrapolating to the infinite correlation to the infinite basis set limit and then using those correlation energies to compute the dissociation energies and that's the value that you get for these dissociation energies they're actually that actually brings you within chemical accuracy which is a nice thing so in a way it's a confirmation if you like that the basis set error is the underlying source of problem here and the other thing which is a different technique is what we call F12 and that introduces what it does and what's a little bit about it tomorrow maybe is to treat the absence of a cusp as a perturbative correction in other words knowing the behavior of the wave function as two electrons come close to each other but are still somewhat distant you can then get an energy functional that well if you like you can extend this wave function down which is obviously an approximation at that point and work out in a perturbative manner what that would mean for the correlation energy and that's what these F12 corrections are now the good thing about the F12 corrections is that they're one-shot corrections so for a given basis set you can correct that number to get a certain value whereas the extrapolations obviously require two calculations and but anyway in both cases we now end up with a satisfactory description of the of these dissociation energies but this is the kind of thing that when you come to do quantum chemistry at least you have to start there is a lot of after you do these calculations there's some amount of work that you have to do if you like to polish things off so to speak it's somewhat tangential to at least the purposes of these that just okay I will okay so how much time 5-10 minutes okay so let me just finish off by looking at some of the structures of these wave functions that we get so what we've done is for a given molecule let's say let's take N2 for example is we've ordered the determinants in descending order of population so we have the leading determinant say that's we call that determinant 1 the second most populated determinant we call that determinant 2 and so on so for descending order we plot the corresponding populations so here you can see the leading determinant has got 50,000 walkers on it because that is the criteria that we used to to in our calculations and these are both on a log log scale so when you get to the 10th most important determinant you end up with a lower population so what I show here are two different geometries let's say for the C2 molecule one in which the C2 molecule is at equilibrium geometry and one in which the C2 molecule has been stretched in this case the 5 times the geometry and the feature that makes these problems difficult is that as you begin to stretch and you start to break these bonds covalent bonds the corresponding correlation problem becomes harder and harder and harder and the way that manifests itself is that you get determinants other than the Hartree-Fock determinant whose importance grows in magnitude and you see that very clearly so for example if you're dealing with the N2 molecule and at equilibrium geometry there's a determinant that's got 50,000 walkers on it and the next determinant which is this one it's got under 10,000 walkers and it's still pretty heavily compopulated but it's got under 10,000 and then as you go to the less and less important determinants you come down down down down until let's say the 100,000 most important determinant has about the initiator threshold so this distribution here is telling you about the rate of decay in the Hilbert space of important determinants so in this problem you would have something like 100,000 determinants that are definitely important now look what happens when you stretch the molecule when you stretch it you end up with a whole series I think of roughly 30 determinants which all end up to be equally important compared to the Hartree Fock so they've all got in this case about 50,000 determinants in them and then you get another series that has got another shell and then you get this precipitous jump down here and it decays off you typically end up with roughly the same number maybe within a factor of 10 or so of important determinants but the crucial point is that the algorithm somehow sorts out what is important and what's not important and distributes walkers accordingly now there's a very interesting effect that you notice that in the... these are the heteroatoms NO, CN and CO and these ones are the homonuclear molecules you notice that in the case of the heteronuclear atoms actually the stretched limit gets easier rather than more difficult so if you take CN initially you end up with more highly populated determinants so it's a more multi-reference problem but the decay in the Hilbert space is actually faster compared to the heteronuclear case so it's a very interesting thing that as you stretch out these two atoms if they're different this atom basically develops its own Hilbert space and this atom develops its own Hilbert space and the determinants essentially don't talk to each other apart from the low energy ones but that effect is not seen when you... when you're dealing with the homonuclear atoms where in the worst case scenario with CN2 you can see substantially more required a factor of 10 in the case of C2 again you start off being more multi-configurational among the important determinants but you end up more or less with roughly the same number of important determinants okay I think I'll probably stop there and if you have any questions I'll take your questions