# How do Quantum Software Packages Work?

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 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 takes 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 a simple math like 2+2 for 6 million times it takes way longer than what these packages do with such complicated integrals!! Anyone knows what their trick is? Do they have 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 the symmetric situations.

• I think many of these packages use highly optimized numerical libraries for linear algebra like LAPACK and BLAS. They may also implement some parallel programming techniques (eg, multi-threading or MPI) to take advantage of multi-core processors. – commanderdata Sep 23 '14 at 1:33
• >The only unanswered part is an algorithm for finding the symmetric situations || actually, gamess puts this load onto user. Gaussian somehow can guess symmetry without it explicitly mentioning in input. It is not very hard, though. D2h and friends is easily detectable for linear molecules, and all elements of symmetry must go through geometrical center of the molecule, tranforming one atom into anouther. Since all symmetry operations are representable in matrix form, it is possible for each to points to concstuct such a matrix and try to see if it is indeed and element of symmetry. – permeakra Sep 28 '14 at 10:48
• I also suggest to google for article "what every programmer should know about memory". It gives interesting insights to optimizations. – permeakra Sep 28 '14 at 10:49
• Are you sure you C++ code was that efficient? Non-CS guys in scientific computing are notorious about making totally false claims about their codes. – Greg Sep 28 '14 at 14:19
• @Aug You should look into the source code of free software computing packages. Psi4 is written mostly in C++ and Python. MPQC is mainly C++. There may be others available as well. There may also be some helpful information in the Q-Chem manual – LordStryker Sep 29 '14 at 12:17

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 use of multi-core processing. But these are actually minor optimizations.

UPDATE

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

• Yes I absolutely agree with the efficacy of rejection of negligible integrals especially in larger molecules as there will be more zero overlaps. I have used Obara-Saika for Overlap Integrals and Kinetic energy. I found that Rys iteration is way faster for e-e integrals . But I couldn't find a definite formula for selection-rejection of the integrals (I only know "Schwarz inequality" but it is not that effective). Do you know any effective selection-rejection criteria? I think the answer of this question may be such criteria. – Aug Sep 23 '14 at 3:09
• This may not be the best forum for your question. I might try ccl.net where there are more method development people. – Geoff Hutchison Sep 23 '14 at 3:28
• I'd add that many times you can exploit molecular point group symmetry. Even just having one symmetry element can reduce the storage by 2 and the work by 4. D2h, which has order eight, can reduce the work by 64! – jjgoings Sep 23 '14 at 19:55
• @jjgoings I think this symmetry sounds really promising but it is so difficult to find symmetric situations ( I mean by coding ) because 4 integrals can come from any different atoms but be symmetric yet. By any chance you may know a reliable algorithm to apply this symmetry? – Aug Sep 28 '14 at 4:37
• >> Yes, there are numeric tricks (using BLAS, LAPACK, etc.) and use of multi-core processing. But these are actually minor optimizations. || I.. wouldn't say so. use of multi-core and parallel runs may reduce job time by orders of magniture, though parallelization involves significant overhead. Also, some packages include GPGPU support. – permeakra Sep 28 '14 at 8:43

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

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, like 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 understanding integral evaluation. He has also written several exercises and tutorials, which you may find useful.

• Thank you for the good suggestions. I wrote some "Homemade codes" to find the symmetry ( it is not very efficient as it can not find signed symmetries/ plane symmetries and also is rather slow) but I cannot say how faster the whole calculation is now ( maybe 10 times faster )! I got surprised and now I know that the main magic is the symmetry and negligible integral rejection. I'm gonna dig into the suggested sources and will post the results here. – Aug Sep 29 '14 at 19:52

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 term as a product of single particle orbital function thus greatly reduces 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 ist 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 types 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 frequency should be positive and for 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 condition. It enables us to solve the problem for minimum non repititive unit cell ( Irreducible Kpoints). Finding the symmetry is also an crucial part of this types of problem.

Cut off energy is used to reduce the number of basis sets considered. And also optimization techniques also has great effect on computational cost, you may use 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 on online), 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 cost dramatically. One another potential method is to use GPU, but it is still in it's primary stage.