Equation Generator for Equation-of-Motion Coupled Cluster Assisted by CAS

Hello everyone, my name is Raul Quintero and today I will talk about an algorithm that generates automatically equations that are important in quantum chemistry. So this algorithm has been made in Mathematica. So the outline of this person. First I will start with a quantum chemistry background. So I will introduce some quantum mechanic concepts and also quantum chemistry concepts. And then with this background we will discuss about the main subject that is the couple cluster equation and equation motion couple cluster. So I will explain why these equations are important and why do we need to generate them automatically. And then I will present the algorithm in, in an interactive way. I mean we will change from, we will change from Latex to Mathematica to see how the algorithm works. So starting with the quantum chemistry background. So what is Quantum chemistry? According ChatGPT is a branch of chemistry that uses the principle of quantum mechanics to understand and predict the behavior of atom and molecules. And according Wikipedia, is the branch of physical chemistry focused on the application of quantum mechanics to chemical system. So yeah, it's basically quantum mechanics applied to chemical problems. So which kind of chemical problems? So here I have three examples so we can understand chemical reaction. So in this example, in the, in the Y axis we have the energy and in the x axis we have the reaction progress. So at the beginning of the reaction we will have two separated products, A and B. On the addition of both will be a given energy. So then when they start interacting they will form a transition state higher in energy and finally they will B they will bind. So those compo. Those. Those elements are, are bounded will from form a compound that is stable with respect to transition state. Another important use of quantum chemistry is a drug design. Whenever a drug is designed, we, we need to target protein, an enzyme or something in our organism. So in this example I'm showing in this orange cycle is the enzyme and this blue dot is a drop. So we can model the interaction energy between both. So it actually helps to optimize the drug candidates. So if we have many candidates, we can select the one with the minimum energy because they will interact, they will interact more strongly and the drop will be more efficient. And the third application is spectroscopy. So spectroscopy is basically the study of the interaction of light and matter. So regarding our subject, this is more related. So in this example I'm showing an emission expectorant. Each of this line represents a transition energy from a higher energy state to a lower energy state. So with respect these three examples we see that we have something in common is that we are, we're measuring the energy somehow. So how do we obtain those energy. So we do it through the Schrodinger equation. Sorry for that. Yeah, we do it through the electronic Schrodinger equation. First of all, what, what electronic. Because atom and molecules are made out of electrons. So the electrons define to find their behavior. So the electronic shooting equation describes the behavior of electrons under wave functions. So in this equation on the right hand side, first we have the Hamiltonian that is composed by F that is also called fog operator. And B is the electron Electron operator. F will contain the interaction between electron and the nucleus. Also the kinetic energy of the electrons and also the electron within mean field. And B will contain just interaction between electrons. F is also called two body, one body term and B two body term. Psi is the wave function is what we want to obtain. And E is the energy of the system. So this is an eigenvalue equation. So when we solve it, we obtain this wave function on the energy. So but we, so we start with a, with a Hamiltonian is what we propose in this equation. And then make the Hamilton acting on the wave function to solve this, this equation. Another key concept in quantum mechanics is the quantization of energy. So here I'm showing again the emission spectra. So for this emission spectra, imagine we have a gas phase atom. So if we apply some, some current, the atoms will emit light. So here we are, we are measuring this light emitted for the atom. And each one of these line correspond the transition from a higher, from a photon from a higher energy state to the lower one. So, so yeah, those, so each one of these energy state will correspond to a different wave function. But as I mentioned before, we are working with electrons. So what are the rules that electrons must follow? So they, they, they follow something that we call Paul crucial principle that says that no two electrons in the same atom can have identical values for all four of their quantum numbers. So yeah, in the right hand side we have this image with two electrons with the same energy but different spin. So if we have two electrons of the same with the same energy but with the same spin is not possible. And for simplicity we'll represent the two electrons with different spin just with dots. But since we have just two spins, we cannot have more than two electrons for a given energy. So another important concept in, in quantum chemistry is the electronic configurations. So basically the electronic configuration is the distribution of electron in different orbitals. In the right hand side, I'm just depicting this is not how atom look like. I'm just depicting an atom and apart two electrons are occupying each of the Internal. More internal orbitals. So in this case we have five orbitals and six electrons. So a simpler way to depict this image is by just using these orbital energies and electrons. So we have a simple image. And in this, this is the ground state configuration, and it is also called the reference configuration, where all the electrons are in the, in the lower orbital energies. So some notation. The the electrons that are occupying the lower energy levels are called occupated states. And the occupated, the unoccupied orbital energies are called virtual states. So I will mention that. So it's important to have that in mind. So, but we can also represent different states. So again, if we start from a ground state configuration, we can perform jumps of these electrons. For example, we can make this electron in the level three jump to level four and from the level one to level five. So when we are doing this single jump, we say that we are obtaining single state configurations. And again, if we start again with our ground state configurations and we perform two jumps simultaneously, we are now talking about double excited state configurations. So here I'm representing some of them. Of course, we have all the possibilities of two jumps. But also if we perform T electric, if we perform a jump of three elections simultaneously, we obtain triple excited state configurations. So how do we represent those downstate configurations mathematically? So consider again the ground state. So if we number the positions in the orbital energies as follows, 1, 2, 3, 4, 5, 6, we can say that we are creating in an empty space. It's not empty space, it's physical vacuum. That is just the orbital energies that we are creating electrons. So for example, here I'm filling this, this space first with electron 1 and electron 2 that are occupying the lowest orbital energy. Then 3 and 4, there is the second one, and then 5 and 6. So in this way we are representing this configuration. So this way of representing this is called second quantization. But it's also important to notice that each one of these this number is an operator. And whenever the operators contain a dagger, that means that they are creating electrons in the vacuum. So really useful shortcore notation for this ground state is just by depicting a zero. So in this way we are also saying that we are in the reference state. So, so just to so be aware, because it's easy to confuse this, this reference state with the physical vacuum. This physical vacuum again is the absence of electron and just orbital energies. But whenever there is a zero, we have electrons all in the ground state. So in general, we can write a ground state configuration for a given number of electron and orbital energies as follow. Here we are depicting a product of these creation operators upon the physical backing and we can also represent excited state configurations. So for example, if we start again from the ground state configuration, if, if we want to represent from it excited state configuration, I need to remove the electron 5 and create electron 7. So I am depicting that here by putting a 5 without the dagger. That means that it's an annulation operator that this five without the dagger will annihilate the five with the dagger. And we are creating an electron in the seventh position. So we are using this seven operator. So after this simple operation, we obtain this excited state configuration in second quantization. So we need to annulate one electron and create another in the reference state. So that's the idea behind creating excited state configuration. So again, this formulation is named second quantization. In general, if we want to write this single excited configuration in general, we can do it by using a pair of creation and relation operators acting upon the reference state. So here I is annulating I. So in this way we obtain this general single state configuration. So another. Another notation note here is that the indexes A, B and C creates particle in the virtual orbital spaces. And also the indexes I, J and K remove electrons from occupied orbital energy. So they create something that we call holes. And if they can act as both, because we have examples where they can act as both, we use the indexes pqr. So again, IBC ABC we will call them particles. IJNK will call them holds is another shortcut notation that is very useful. But how can we relate this configuration that I have described with the Schrodinger equation? So for that we need to use a concept within quantum mechanics, that is the quantum superposition that refer to the capability for a given state to exist in multiple states concurrently. So for that we will use the the Schrodinger cuts. This mental experiment consists in a cut inside a box that is dead and alive at the same time. So if we want to, if we want to build the Schillinger equation of that, we need to sum the alive gut state and the dead cut state. So also we can see that there is a number multiplying each one of the states. So it numbers just refers how much present is a given state in our wave function. But okay, in our case we are. We want to represent atom and molecules. So instead of using this cut, we will use the electronic configuration that I already described. So in this example I'm representing a genetic wave function. So these numbers that I Will call coefficients are multiplying each one, each one of these configuration states. So actually when we obtain the Schrodinger equation, what we obtain are those, those coefficients. So we know how to solve, we know how to obtain the shooting, the, the wave function that is by the Schrodinger equation. So let's come, let's consider it again. So if we solve this eigenvalue problem, the eigenvectors of this solution will be the different wave functions and the eigenvalues will be the energy. So, so in this example I'm showing that when I'm solving this, the first eigenvector will be the ground state wave function, that its eigenvalue will be the ground state energy. So the first excited wave function, it's over here. And the energy, the energy associated with, with this wave function will be the second second energy, and so on and so forth. So actually there is a method to obtain all this wave function and all of this, all of these energies. And it's called full CI, but it's really expensive. So actually we can, it, it can just be used for like small molecules or small atoms. But of course we have, let's say more economic way to also obtain all of the wave functions and the energies. This method is called coupled cluster. An equation of motion coupled cluster. So now we can talk about them with the backgrounds. So the textbook definition of coupled cluster is that it's a wave function that is built as an exponential parameterization of the reference state. So here the argument of the exponential is the T and T is a mathematical object that operates in the reference reference state configuration to generate excited configurations. So this T is also expanded in single excitations T1 double excitations T2, triple excitations T3. So let's take an example. What, what does the T2, what does the T2 do to the, to the reference state? So what it does is it creates all the all the possible double excited state configurations. So we are creating a zoom here of all possible second state configurations. So here we have, remember that we have a part of creation relation. So we are performing two jumps, but we also have something different with respectful CI, that is the amplitudes. So let's say that this amp. So for one, for one given jump we will have an amplitude. So if we apply the Copper cluster to the Schrodinger equation, we will obtain an expression like this. So here I also multiply by the left this exponential minus T. So whenever this Hamiltonian is enclosed, so here the exponential minus T is multiplied this Hamiltonian. And then we have the exponential T. So if we have this Hamiltonian between these two exponential, we call it similarity transform Hamiltonian. So I will talk about it a little bit later. So how do we solve these equations computationally? So first we, we need we truncate the the exponential up to fourth order and then we truncate the expansion. So in this example I'm truncating T equal to T1. So I'm obtaining couple cluster singles. The singles come from single excitation. But if, if I have T equal to T1 plus T2 will be coupled cluster singles and doubles and so on and so forth. And we can also have T equal to T2, that will be coupler cluster doubles. And then after, after this approximation we need to project the equations here. I'm projecting the equations in the reference state and also in the single state configuration. So if we have this set of equations, what we're doing with the first equation is obtain the energy. And in the second equations we will obtain the amplitudes of the single excited state configurations. So this amplitude will correspond to the single excited state. So in the end this amplitude will contribute to the, to this energy. So let's say that we will have a correlated energy by by these T1s. So we can obtain all the wave functions with this method or for with this way of solving them. But it's very expensive. So instead what we do is just obtain the minimum energy that correspond to the ground state. So this way of doing it is achievable even for big molecules. So we obtain the ground states. And the ground state is very important in quantum chemistry for a lot of reasons. So yeah, with copper cluster we obtain just the ground state energy that is depicted here, but not excited states. So here you see that I change the C from the T's that are the amplitude. So when we saw we obtain these amplitudes and when we apply the Hamiltonian to, to this wave function, we will obtain this correlated energy by this amplitude by these amplitudes. But how do we obtain with couple cluster the excited state wave functions? So we do it with equation of motion couple cluster. So after obtaining the solution of coupled cluster the amplitudes of the wave functions, we need to define another operator that is our operator. So we also need to truncate this operator. So the order of this operator will be devices in which we represent the similarity transform Hamiltonian. For example, if I truncate R, if R is equal to R1 and R2, I will just truncate this expansion to two. So we will have this operator in a basis of S and this. So we will just have 4 by 42 by 2 matrix. So we just need to solve this this equation to obtain the eigen energies that are the transition state and the eigenvectors. Now with eigenvectors we can build the excited state wave functions. So yeah, to summarize this, what we need to obtain first is the orbit and energies. We do it with a method that we call Hardy from. And then that we have our reference that is formed by these orbital energies. We can do copper cluster in coupled cluster. We select these excitations. We truncate this T to obtain the similar transform Hamiltonian the ground state energy and also the ground state wave function. And then we do the equation of motion coupler cluster when we define R to obtain the transition energies and also the excited state energies. But okay, I haven't answered yet. Why do we need to use a program to generate these equations? So I prepare a small example. So this, this element here will correspond to the 2 comma 1 element of the equation of motion matrix. So 2, the 2 will correspond to d and the 1 will correspond to S. So to simplify everything, we will just consider t equal to t1. So again, following the recipe, if we do the Taylor expansion, we can truncate it here up to third order, because fourth order is zero and zero order is also zero. So we have this expression. But okay, let's take again this the simplest example of this expansion. So then we need to express this block W2 comma 1 in second quantized form. So when we do it, we obtain first, these double excited configurations will be this with the string, the strings of second quantization operators and the H, if you remember, it will be splitted in the F term and in the B term the one electron and two electron operators. And then we have the T1 that it will span all the possible single excited state configurations. And then S that is just a single excited state. So how do we evaluate this element to obtain for example, results of atom and molecules? So what we use is combinatorics. So for example, here I take the example of the two electron terms. So we need to use this formula. Actually this formula fulfills the Pauli exclusion principle. So what we do with this formula, it's to generate all the possible combinations between the different sets. Since we have 12 operators, if you can count them in four curly brackets, we will obtain combinations. So that's the reason what we are using a computer to generate these equations. But this is even for a Computer, I am just giving a really simple term. We have a much more complicated term, so this will become really difficult even for a computer. So we have a smarter way to do it. So yeah, this is another way to define this second energy term. So what I'm doing here is just form this all possible bar combinations just between the first set of second quantization operators and the last set. If we form all the possible combinations, we obtain just eight terms and we can notice that the first four terms are really similar up to exchanging of indexes. So indeed we can just define one. So in this case I'm defining one called ia. So this is, this will be my definition of this first term. So you see, if I exchange the term I will, I exchange the I from the J, I will obtain this term, the A for B, this term, the I and the A for J and B, I will obtain this term. So I will need to represent this expression and the same case for the blue and the green terms, we just need to represent them once together with these equations. So we can say that we have just three unique many body term or intermediates. So those terms are super important. We base our work in this in those terms. So for, for the next section, section I want to, I want to share Mathematica notebook to see practical examples. So I will start with the general algorithm. So the input will be the cluster operator and equation of motion operator are those that I already described. So in the next step we transform this cluster operation or equation motion cluster operators into second quantized form of this equation of motion matrix that I already showed you. And I also obtain the amplitudes from the cluster operators. In the next step we use as an input this equation of motion matrix in second quantization to obtain the many body term that we, that we saw in the last part. And then we use this intermediate and second quantization. And in this step is actually the harder one of the algorithm. We do all the possible combinations that I, that I mentioned before, the power combinations. And finally we need to operate this delta clinical in the summation to obtain everything in terms of the amplitudes and the integrals, the 1 and 2 electron integrals that I described at the beginning. So yeah, this is the input. So here I'm using this particle connotation. So if I'm using this, it exactly means that I having T1 and T2 and in the equation of motion operator R1 and R2. So it's a bit. So yeah, what we are doing here is just saying, okay, we want to work with this, with this matrix two by two, because the equation of motion operator contain R1 and R2. Okay, so then we'll generate the amplitude and string. So if I. If I take this d the two, one block that I already mentioned in the example. So in this step we are transforming the operators into second quantized form of this matrix. So for example, if we. So we store everything in this variable called strings. So for the one body operator and for. For the block 2 comma 1, we obtain this string in the second quantile form. You can notice that the indexes are different because we are interested in obtaining like higher orders, like up to kintuples quadruples. So if I use I, J, K, they start to mess between each other. So it's better. This notation is more systematic, but it's more systematical. But it's the same idea. And then what I do is I obtain the many body terms. So just need to run that. And yeah, what I'm doing here is. This expression in term, in terms of these many body terms will be the following. And it's a total because I have everything stored in a list. And we also obtain the many body strings that are those many, many body terms, but in second quantized form. And I will close that. And then we need to generate all possible word combination. So yeah, this one is taking a bit more because I took a long string just to show. But yeah, here is. Here I have already the output. So you can see that whenever we have like a false entry, we already put a new to evaluate it again. So yeah, here I'm showing all the possible part combinations that we obtain. So we also have a sign associated to these part combinations. And yeah, we obtain a lot. And yeah, the final steps, we generate these many body terms in terms of the integrals and the amplitudes. So yeah, we obtain these equations. So yeah, the algorithm is like this. It's. If you are not really familiar with the notation, could be a little bit complicated. For the first time I understand that. So I would run. So yeah, I have examples of the runs, some examples that I already run, I haven't yet. But yeah, let's. Let's obtain the couple cluster singles and doubles. So we need to define singles and doubles that correspond T1 and T2. And the equation of motion operator is not important here. So we define it as 00. So now we are running the calc. No, it's taking a lot because I have a lot of things open here. Wait a bit until we obtain the result. Okay, we can wait for the result. And yeah, for the. For the equation of motion coupled cluster, we need to define the two operators T1, T2, R1, R2. Here I already run. Run the calculation. So we already obtained the expression. So first we obtain the many body term. This is the execution time. I don't know why. The other example is taking a lot. It's already finished. So yeah, we obtained the couple cluster equations here. So first we obtain the energy. First we obtain the expression when we project in the single excited state. And then with the double excited state. Yeah, for the equation of motion coupled cluster, we are using these many body terms and then the expression of the equation of motion. So we see that here we have this many body terms already here. So we need to just define them and define these equations and define. Go together to obtain the equations. I have actually examples how they should be programmed and I also have numerical examples. I'm sorry, I just wrote the numerical example. So if we want to calculate the. The excited state and the ground state of helium, we can use this, this model that is. That is equation of motion numeric. So we have two electrons, two orbitals. The atom is centered here. I'm. I'm using cartesian coordinates and we need to. We need to provide with some integrals and also with the. With the methodology that we're using. So yeah, first as we first obtain the orbital energies. So here are the. We call the molecular orbital energies. Then it's doing the copper cluster singles and double equations. So you can see that is iterating until it obtains the minima, the minimal energy of the ground state. And then it start with coupled cluster. Equation of motion. Coupled cluster. Sorry. And this is the matrix. I'm just printing it. This example is very small. So for. I have another example that is bigger. But yeah, the eigenvalues that we obtain are the transition energies. So yeah, like actually so I have a more complicated example that is water molecule. So here it's a little bit bigger the input because I need to define the coordinates. So again I obtain with hydrofog the molecular orbital energies. And then it needs. It takes more time to. For this iteration to obtain the ground state energy of couple cluster singles and doubles. And look, this is the matrix of equation of motion. This time is huge. So yeah, we need to diagonalize this matrix. It's very big. So in the end we obtain all the transitional energies that mean we obtain the spectra and from. From. From my side. That's it. Thanks a lot for tuning. So. Yep, thanks.