How do quantum software packages Work?

Question

If one wants to calculate a moderate size Alkane (with say 10-15 carbons, assuming 100 electrons, with Restricted Hartree Fock based methods) we can simply say that the electron-electron part will be $100^4 = 100$ million integrals. If we remove the non-unique ones it will become $\frac{100M}{8}= 12.5 M$ Integrals. If we do a pre-screening it may become around 6 million. Even if we go for DFT it will be around $\frac{100^3}{8} =125,000$ .

I have calculated a small system like Methane with just 10 electrons (around 400 unique/pre-screened integrals) with an efficient C++ code (based on Rys Integration which is among the most efficient methods) and it took around 1 hour. When 400 cycles take 1 hour I can assume that 6 million cycles may take years! But software packages do that just in seconds!

My question is if I just want to do simple math like 2+2 for 6 million times it takes way longer than what these packages do with such complicated integrals!! Does anyone know what their trick is? Do they have a pre-calculated database or some approximation methods?

UPDATE: there are many good hints in this question and thanks to all for that. A combination of symmetry, finding negligible integrals, and using optimized methods work fine. The only unanswered part is an algorithm for finding symmetric situations.

Answer 1

There are many tricks. Since I'm not a developer, I'm not privy to all of them, but there are some simple things to start.

• Some integrals are predictably close to zero: Even in fully first-principals calculations, there are parts of the code that prune out integrals that will clearly be negligible. This cuts out a lot of work. Consider, for example, overlap or Coulomb integrals between atoms that are very far apart.
• Gaussian Orbitals: Using Gaussian-type orbitals are incredibly fast because most quantum integrals can be derived as simple recurrence relationships. That is, you almost never need to do the integration. See, for example, Obara and Saika and the recent PRISM algorithm from Gill and Pople.

I believe there are some actual quantum developers on here and perhaps they'd be willing to share other tricks. But these are big - eliminate any integral you don't need, and don't actually take the integrals in the first place. ;-)

Yes, there are numeric tricks (using BLAS, LAPACK, etc.) and the use of multi-core processing. But these are actually minor optimizations.

UPDATE

I realized that I forgot to talk about the convergence of the self-consistent field which is also optimized. There are various convergence strategies, most notably DIIS, which help to extrapolate the final orbital coefficients and electron density.

Answer 2

First things first, I would check your code: there is no reason 400 unique integrals should take an hour, even with the most naive implementation. You are computing at a rate of an integral every ten seconds. While I agree with all the comments pointing to increasing parallelization, I think efficient scalar code is more important. Why scale up inefficient code? Considering GPUs is a whole 'nother beast, and from the problems, my coworkers have expressed trying to program them I'd recommend avoiding them for the time being.

The best thing you can do is to exploit molecular symmetry. Heads up though: this won't fix everything, because many interesting molecules lack symmetry. But for methane, as you tried, this should help a ton. Of course, this isn't as simple as just looking up some algorithm, because depending on your program you need to rotate the nuclear orientations to be symmetric, or develop routines to detect point group symmetry, etc. That being said, there are --- as far as I know --- two primary ways of exploiting point group symmetry for integrals. The first is to create symmetry-adapted basis functions, so your basis set is symmetric, which means your wave function/integrals/Fock matrix/whatever will be symmetric as well. Ernie Davidson did this in the seventies. A slightly more modern reference by Helgaker and Taylor is given here. Daniel Crawford's electronic structure tutorial explores this as well (it's really well written --- a great introduction).

The other way you can do it, and I believe this is more popular, is to create symmetry unique integrals and then build a partial Fock matrix only to take care of symmetrization later. If you see something like "generating petite list" in the output of your favorite electronic structure program, this is probably what it's doing. This method allows for integral screening more efficiently, I believe. Dupuis and King did this in the late seventies, based on earlier work by Dacre and Elder. These are the standard references you'll find. There may be some more recent work that I'm unaware of. The problem of point group symmetry had been effectively solved by 1980.

Of course, many interesting molecules don't have symmetry to take advantage of, so symmetry won't fix all your problems. The Schwarz equality is well addressed in Molecular Electronic Structure Theory (actually, most electronic structure theory is well addressed in that lovely pink book). If you haven't checked that out, please do. I think it will answer many of your questions about how modern codes work.

I don't want to repeat the other answers (Geoff gave you the papers I'd have suggested as well), so I'll finish with two lesser-considered references. First, there is an awesome (but tiny!) book called Parallel Computing in Quantum Chemistry that lays out the basics for efficient code including integral evaluation. It's not too hard to read, and I found the work on local MP2 to be pretty interesting. I'm often surprised more people haven't heard of it. The other place I would check out is Ed Valeev's LibInt package. He's a pretty talented programmer, and it's used in many modern packages like Psi and MPQC and ORCA, not to mention free for anyone to download and use. Some of my coworkers have played around with it for their own "homemade" electronic structure code, and speak highly of it. The user's manual is actually quite helpful in understanding integral evaluation. He has also written several exercises and tutorials, which you may find useful.

Answer 3

To further reduce storage requirements, integrals may be calculated on the fly, employing different caching strategies. For some methods storing all the integrals is plainly impossible.

Answer 4

There are two things here:

1. Quantum chemistry problem
2. Optimization problem

Quantum chemistry or the Schrodinger equation is simplified using DFT and some other techniques:

a) ECP or effective core potentials and pseudopotentials: Instead of solving for all electrons just the outer shell electrons are solved and core electrons are considered using pre-calculated integrals
b) ri approximation: Resolution of identity, it treats two electron terms as a product of single-particle orbital function thus greatly reducing computational cost

And as mentioned above there are some other good techniques are used to deal with quantum chemistry part.

Optimization problem: In optimization usually second derivative (Hessian matrix) is used to find minima or maxima but in computational chemistry code only the first derivative or jacobian is used, thus it is highly sensitive to initial guess and computational setup.

Calculating Hessian is the most costly operation of computational chemistry code, either approximate hessian is used or no hessian is used in this type of optimization problem. If you want to make sure that you got the right geometry you can calculate Hessian with your optimized geometry and for geometry optimization, all the frequencies should be positive and for the transition state, only one frequency should be negative but you will find that to meet this requirement is highly demanding.

Moreover to reduce computational cost we can apply periodic boundary conditions. It enables us to solve the problem for minimum nonrepetitive unit cell ( Irreducible Kpoints). Finding symmetry is also a crucial part of this type of problem.

Cut-off energy is used to reduce the number of basis sets considered. And also optimization techniques also have a great effect on computational cost, you may use a fast but aggressive algorithm but that may not be able to find the minima or maxima.

You can read VASP or TURBOMOLE manual (available freely online), and there you will find some good techniques for reducing the computational cost. And above all, as it is mentioned above efficient implementation of parallelization can help to reduce computational costs dramatically. One another potential method is to use GPU, but it is still in its primary stage.