How to calculate FCI matrix elements from an FCIDUMP file

I'm interested in writing a simple program (as an exercise) that takes an FCIDUMP file (e.g. from a Molpro calculation) and can determine FCI matrix elements (e.g. to construct the entire FCI matrix to diagonalise, or to use in FCIQMC). I would like to know how to do it so that I can program it myself.

Reading in the FCIDUMP file itself is just an exercise in programming (and can be done in one line with pySCF: pyscf.tools.fcidump.read("fcidump_file", molpro_orbsym=True)), so let's start with all the integrals given $$(\psi_i\psi_j|\psi_k\psi_l)$$ and $$h_{ij}$$, as well as the extra parameters in the first block (which I am less sure how to use). How does one calculate the (FCI) Hamiltonian matrix element $$H_{ij}$$ given these values?

I'm aware of the Slater-Condon rules and know they will be needed here, but I'm not sure how to generate a general matrix element from these values, or how the parameters that come before the integrals are useful (presumably ORBSYM can reduce the size of the matrix).

• I figured that was a small component of the question that could even be answered with a comment, but if you think it warrants an entirely new question, I can do that. I am more interested in how to construct it myself, though. – tmph Mar 15 at 22:17
• +1, but the second question about existing software that can do this, should be asked as a separate question with the software-recommendations tag. You're right that ORBSYM tells you information that can help to reduce the overall cost of things. So you have the Hamiltonian (i.e. integrals), and you want to calculate the matrix elements $H_{ij} = \langle \Phi_i | H | \Phi_j \rangle$ where $|\Phi\rangle$ are Slater determinants with any number of excitations? Have you tried to see what happens when the operator $H$ acts on a wavefunction $| \Phi_j \rangle$? – Nike Dattani Mar 15 at 22:20
• Assuming $|\Psi\rangle$ is a determinant, then the single-particle part of the Hamiltonian is easy (you have an eigenket), but the $r^{-1}_{12}$ part I don't know how to handle beyond what Szabo/Ostlund does in chapter 2.3.4. FWIW, I'm also not sure how to enumerate all possible determinants given only an FCIDUMP file, much less symmetry-restrict them, but yes ultimately what I want is the matrix element between two determinants with any number of excitations. – tmph Mar 15 at 23:38
• The number of possible determinants before considering spatial symmetry is (N choose a)(N choose b) where N is the number of spatial orbitals, a is the number of alpha electrons, and b is the number of beta electrons. How to reduce this number using the point-group (symmetry) ought to be the subject of another, separate question. As for the rest, maybe start with the 2nd quantized H: docs.microsoft.com/en-us/azure/quantum/user-guide/libraries/…, $a$ and $a^\dagger$ will change # of excitations in deterimants. – Nike Dattani Mar 16 at 0:07
• Thank you, I will experiment with a (small) matrix without symmetry adaption for now, since this is not meant to be a production-ready program. You are right that the problem is a lot easier in second quantisation, didn't occur to me somehow. – tmph Mar 16 at 2:33

Let's go back to the original paper where the format was first proposed: A determinant based full configuration interaction program by Peter Knowles (lead author of the MOLPRO program, along with Hans-Joachim Werner) and Nick Handy in 1989. This algorithm is still used by a lot of people today, and is considered one of the most efficient ways to do FCI. Some people call this a string-based approach.

In a "string-based" approach, Slater determinants are represented by a binary string, such as:

$$\tag{1} |\Phi\rangle =|11111111110010000000000000\rangle.$$

Here the 1s, 2s, 2px, 2py, and 2pz orbitals are all doubly occupied (the first 10 ones), and the 3s orbital is completely empty, then the 3px orbital has one electron, then everything else in the basis set is empty. This is probably a double-excitation of the HF determinant which would look like:

$$\tag{2} |\textrm{HF}\rangle = |11111111111000000000000000\rangle.$$

You can now enumerate all possible "strings" that you wish to allow, based on your molecule's spatial symmetry and spin symmetry (for example), and label these determinants as $$|\Phi_i\rangle$$ where $$i$$ goes from 1 (for the reference determinant, which is usually the Hartree-Fock determinant) to the number of determinants in your list.

Now, the question you ask is how to get the matrix elements:

$$\tag{3} H_{ij} = \langle \Phi _i|\hat{H} |\Phi_j \rangle,$$

for the following Hamiltonian operator:

$$\tag{4}\hat{H} = \sum_{pq} h_{pq} \hat{a}_p^\dagger \hat{a}_q + \frac12 \sum_{pqrs} h_{pqrs} \hat{a}_p^\dagger \hat{a}_q^\dagger \hat{a}_r \hat{a}_s + \hat{H}_\textrm{nuc}.$$

You have all the one-electron integrals $$H_{pq}$$ and $$H_{pqrs}$$ from the FCIDUMP file already, so now you just have to apply the creation and annihilation operators to the bit-strings, remembering that single-excitations will remove an electron (i.e. a 1 in the bit-string) from its present orbital and place it into one of the unoccupied orbitals (i.e. a 0 in the bit-string), and a double-excitation will do this twice.

• I would not recommend doing this for a large basis set or a large number of electrons, because your number of matrix elements will explode very quickly.
• I would also recommend that while writing a program to do this can be a good exercise for a beginner, it would be best to stick to a well-established software if you want to do any serious calculations, rather than re-inventing the wheel.

In your question you said that you might try to implement FCIQMC. If so, please do not store the whole $$H_{ij}$$. The whole point of FCIQMC is that you do not sore the whole FCI matrix, but you effectively only use parts of it as they are needed. There's an open source FCIQMC code called NECI (N-electron Configuration Interaction) which can be downloaded from GitHub here and there's also a paper about it for which I'm a co-author. There's a lot of work that went into this software, so again, be careful to balance your education journey with practicality in trying to accomplish research goals.

• Wow, this is an amazing answer, thank you so much!! In this answer, did you address the use of the FCIDUMP parameters MS2,ORBSYM,ISYM,OCC,CLOSED somehow? I guess NORB and NELEC determine (half) the length of the string and the number of ones respectively. Also, I actually use NECI already. I am primarily trying to understand it by making minimal examples, which is exactly why I asked about matrix elements and not just the FCI matrix. I also appreciate your advice. :) – tmph Mar 16 at 2:31
• Thanks for the kind words! I didn't address those FCIDUMP parameters in my answer: I guess I assumed that the reader knows what MS2,ORBSYM,ISYM,OCC,CLOSED,NORB,NELEC mean. I suppose a new question with the title "What do the parameters at the top of an FCIDUMP file mean?" would be valuable to a lot of future users as well. Also, are you using the GitHub version of NECI, BitBucket version, or the version that comes with MOLPRO? – Nike Dattani Mar 16 at 2:41
• I have an idea of what they mean, but don't really know how to use them. For this string representation, would it be easier to use two strings, one for beta and one for alpha electrons? Seems simpler to take care of spin symmetries and matching the indices in the FCIDUMP file with the positions in the string that way, e.g. $|1100\rangle=|10\rangle|10\rangle$ (one of each spin in orbital 1, empty orbital 2). As for NECI, I use the BitBucket version. – tmph Mar 16 at 15:00