# Which software is good with generally contracted basis sets?

Recently, I learnt that basis sets can be contracted as segmented (like def2 or 6-31G) or general (like cc-PVTZ, or ANO). The general contraction allows all the primitives to appear in all shells, whereas segmented does not do that. Most books on quantum chemistry also mention that QM programs usually cannot handle general contraction very well, and instead they duplicate the primitives for each shell, and computes the same integrals multiple times.

However, there is not much written about which QM programs can actually handle general contractions well, and which programs can't. The books mention a method called screening which prevents the same integral being calculated twice.

The manuals of the common softwares (Gaussian, GAMESS, orca etc.) don't mention clearly if they have any such method. (To add to the problem, main parts of the manuals are often old, and has not been updated as the program got updated). For example, GAMESS manual mentions an 'integral screening' but that seems to refer to ignoring 2e-integrals that contribute very little to the energy. Gaussian manual mentions that general contraction basis sets can be transformed to remove primitives that have $$\mathrm{<10^{-4}}$$ coefficients, but I am not sure if this is the same as screening. Orca mentions that segmented basis sets are faster, but doen't mention if it has anything to deal with general basis sets.

So, my questions are - 1) Which programs at present can properly use generally contracted basis sets (without duplicating primitives)? 2) If this can be done simply by checking if an integral has been calculated before, then why is it difficult to implement?

# General vs segmented contraction

Unless one is using a fully decontracted a.k.a "primitive" basis set, one uses contracted Gaussian-type orbitals (cGTOs), which are given as a linear combination of the primitive Gaussian-type orbitals (pGTOs):
$$\chi_i^\text{cGTO}({\bf r}) = \sum_{\alpha} d_{\alpha i} \chi_\alpha^\text{pGTO}({\bf r})$$.

Alternatively, you can write the cGTOs in terms of the pGTOs as a matrix equation: $${\bf X}^\text{cGTO} = {\bf X}^\text{pGTO} {\bf D}$$.

Let's assume our primitive exponents are given as a vector $$(\alpha_1, \alpha_2, \dots, \alpha_N)^{\rm T}$$ with $$\alpha_1 > \alpha_2 > \dots > \alpha_N$$. In a generally contracted basis set, $${\bf D}$$ typically looks like

$$\boldsymbol{D}=\left(\begin{array}{cccccc} x & \cdots & x & x & 0 & 0\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\ x & x & x & x & 0 & 0\\ x & x & x & x & 0 & 0\\ x & x & x & x & 1 & 0\\ x & x & x & x & 0 & 1 \end{array}\right)$$

where $$x$$ marks a non-zero entry. You have $$K$$ primitives (rows in D), which are contracted to $$N$$ atomic orbitals (columns filled with the $$x$$s), and then you have freed up the few most diffuse functions as primitive functions to be able to describe changes in the valence region.

Often the $$N$$ atomic orbitals are orthonormal to each other ($${\bf D}^{\rm T} {\bf SD}={\bf 1}$$ where $${\bf S}$$ is the overlap), but this is not always the case; IIRC the correlation consistent sets for the transition metals have significant overlap between the contracted functions, since the first functions come from atomic Hartree-Fock calculations, while the last come from configuration interaction calculations and are not orthogonal to the Hartree-Fock orbitals.

The thing about generally contracted basis sets is that there is some flexibility in the way you can pick the functions. All that matters in the end is the span of the basis, so you are free to make transformations to your basis as long as the span stays the same. This means you're free to substract columns from each other: for example, if $$x$$ and $$y$$ form a linearly independent basis, $$x$$ and $$x-y$$ also span the same basis.

Now, note that I have two free primitives in the basis, which means these functions can be fully represented by the basis. I can just drop these functions from the general contraction, since I can eliminate the coefficients by substracting the free primitives with the correct weight from the contracted functions. (This is called a Davidson rotation and is what "optimize general contractions" mean on the Basis Set Exchange; see Chem. Phys. Lett. 1996, 260, 514−518.) I now get $$\boldsymbol{D}=\left(\begin{array}{cccccc} x & \cdots & x & x & 0 & 0\\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\ x & x & x & x & 0 & 0\\ x & x & x & x & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 \end{array}\right)$$

But, it does not stop here: I can eliminate entries in the fully contracted block, too, by adding and substracting vectors. Let's say the contracted block looks like this: $$\boldsymbol{D}=\left(\begin{array}{cccc} x & x & x & x\\ x & x & x & x\\ x & x & x & x\\ x & x & x & x\\ x & x & x & x\\ x & x & x & x\\ x & x & x & x\\ x & x & x & x \end{array}\right)$$

I can take the first column from the left and substract it from the other ones to make the coefficient of the tightest exponent vanish. Next, I take the second column and make the coefficient of the second-tightest function go away.. and keep going until I've made it to the right-most column.

I can also do the same thing in reverse: I take the last column and make the most diffuse coefficient vanish in the other columns, and repeat until I am back on the right. Then I get $$\boldsymbol{D}=\left(\begin{array}{cccc} x & 0 & 0 & 0\\ x & x & 0 & 0\\ x & x & x & 0\\ x & x & x & x\\ x & x & x & x\\ 0 & x & x & x\\ 0 & 0 & x & x\\ 0 & 0 & 0 & x \end{array}\right)$$

Suddenly, you see that the contraction matrix has developed quite a few zeros. If the contraction length is not very long, $$N\approx K$$, the columns are mostly zeros (the extreme case is $$N=K$$ where you can choose $${\bf D}={\bf 1}$$ i.e. to just use uncontracted primitives).

Moreover, one can get even more sparsity by discarding small elements $$x$$ in the columns, and (possibly) reoptimizing the contraction coefficients and exponents for a pre-specified sparsity pattern, e.g.

$$\boldsymbol{D}=\left(\begin{array}{cccc} x & 0 & 0 & 0\\ x & 0 & 0 & 0\\ x & x & 0 & 0\\ 0 & x & 0 & 0\\ 0 & x & x & 0\\ 0 & 0 & x & x\\ 0 & 0 & x & x\\ 0 & 0 & 0 & x \end{array}\right)$$

Frank Jensen has written a very nice article on this a few years ago, see J. Chem. Theory Comput. 2014, 10, 3, 1074–1085.

# Integral evaluation

What does this have to do with integral evaluation? The two-electron integrals $$(pq|rs)$$ can be evaluated analytically for Gaussian basis functions. However, the analytical integrals are known only for the primitives, which means that to get the integrals for the contracted functions one has to calculate $$(pq|rs) = \sum_{\alpha \beta \gamma \delta} d_{\alpha p} d_{\beta q} d_{\gamma r} d_{\delta s} (\alpha \beta|\gamma \delta)$$, instead, where $$\alpha \beta \gamma \delta$$ stand for the primitive functions.

If you have a segmented basis set, where there are a lot of zeros in the contraction matrix, it is more efficient to evaluate the integrals one shell quartet at a time; that is, you recompute all the primitive integrals for each $$(pq|rs)$$ you're aiming for, because most primitives only contribute to a single contracted basis function. This means that you have simpler code, but also has the additional benefit that you can easily perform the contractions within the two-electron integral evaluation algorithm, which makes some steps cheaper in the computation, depending on the used integrals evaluation algorithm. Indeed, most codes assume this approach; e.g. this is how Gaussian, Psi4, Q-Chem, and Orca operate.

However, if you have a generally contracted basis set (especially when $$K \gg N$$ and $$N \gg 1$$), the programs that rely internally on segmented basis sets run into troble, as instead of computing the $$K^4$$ primitive shell quartets, you basically need to compute the $$K^4$$ quartets $$N^4$$ times. Screening for combinations of primitives that have small contributions helps a bit, but this is still extremely slow.

As far as I understand, the generally contracted programs compute all primitive shell quartets once, and then form the $$N^4$$ contractions with the $${\bf D}$$'s per the equation above to get the target integrals $$(pq|rs)$$. (Disclaimer: I have never written a general contraction code.)

Programs that handle generally contracted basis sets well are far fewer, but they do exist. For instance, OpenMolcas handles general contractions. PySCF is also very efficient with generally contracted basis set. Both programs are open source and freely available. NWChem is also open source and freely available, and claims to support general contractions (see the opening paragraph of their documentation on basis sets), but it's so slow that I'm not sure whether this is really the case; see below.

On the commercial side, IIRC Molpro has excellent support for general contractions, and is faster than OpenMolcas.

# Benchmarks

To see whether your program supports general contractions is quite simple: just try running e.g. a transition metal system with a correlation consistent basis set. For instance, in cc-pVDZ the first-row transition metals have a [20s16p8d2f|6s5p3d1f] contraction pattern, so $$K=20$$ and $$N=6$$ for s orbitals, $$K=16$$ and $$N=5$$ for p orbitals, and $$K=8$$ and $$N=3$$ for d orbitals, which make them quite heavy.

Although I'm usually against mindless benchmarks, here one makes sense. I'll take the Zn atom and Zn2 dimer (R=2.3 Å) as examples, since they're nice closed-shell systems.

I performed all benchmarks on my laptop with 4 cores, using Fedora 32 packages: psi4-1.3.2-2.fc32.x86_64, python3-pyscf-1.7.5-1.fc32.x86_64 employing the optimized integrals library qcint-4.0.5-1.fc32.x86_64, nwchem-openmpi-7.0.2-1.fc32.x86_64, and OpenMolcas-20.10-1.fc32.x86_64.

## Psi4

As I said above, Psi4 is an example of a program that assumes segmented basis sets.

Setting scf_type = pk to enforce exact integrals instead of the default use of density fitting (which would make the calculation a LOT faster), the HF calculation on the zinc atom took about 19 seconds wall time. Final energy -1777.84665520791896

Running the dimer took considerably longer: 216 seconds of wall time. Final energy -3555.64477390092725

## PySCF

The PySCF calculation on the Zn atom runs in about 0.7 seconds of wall time. Final energy -1777.846655207936010

The Zn2 dimer took about 1.4 seconds wall time. Final energy -3555.644773900906785

## OpenMolcas

Running the Zn atom took about 2 seconds. Final energy -1777.846655208

The Zn2 dimer took about 5 seconds. Final energy -3555.644773901

## NWChem

After running for 5 minutes on four cores for the Zn atom, without seeing even one iteration of SCF, I cancelled the execution...

• Thanks for the great answer! One question, did you run the benchmarks with direct scf or disk storage? – Shoubhik R Maiti Nov 25 '20 at 12:00
• @ShoubhikRMaiti neither: these are all in-core calculations; there are so few basis functions that everything fits into memory. – Susi Lehtola Nov 25 '20 at 18:55

Integral Screening

It is possible to place upper bounds on integrals based on the Cauchy-Schwartz inequality: $$\left( a a | b b \right) \ge \left( a c | b d \right)$$ If the first integral is computed and deemed insignificant, the latter is as well. One can further combine this with the appropriate density matrix element and decide that the resulting contribution to the Fock/Kohn-Sham matrix is insignificant and can be discarded. Obviously, this applies to the case of a "direct" SCF, in which most or all integrals are recalculated on every iteration.

I have not yet heard of other uses for the word screening in relation to integrals. It is not entirely unreasonable to discard a primitive from a contracted basis function if its coefficient is too small, but I think this is a decision to be taken by the user, that is, by using custom input of the basis set.

Why recalculate integrals (If this can be done simply by checking if an integral has been calculated before, then why is it difficult to implement?)

There are too many integrals to keep them in main memory. Some programs may write some to disk, but there is a non-zero writing/reading time cost associated with that. Some programs try to be smart and put the most expensive integrals on disk.

Which programs at present can properly use generally contracted basis sets (without duplicating primitives)? This is hard to judge. Some programs may provide the references to the integral algorithm used (In case of Gaussian, I found: M. Head-Gordon and J. A. Pople, J. Chem. Phys., 89 (1988) 5777), where it may take an expert to decided whether the algorithm is geared towards general or segmented.

The reason why some algorithms work well with one but not the other is that in case of general contraction, one must keep different intermediates in memory to be efficient. If the basis set is not constructed accordingly, one ends up not benefitting from keeping the intermediates for which other savings/more efficient ways were discarded. I found the PRISM algorithm paper illuminating, together with the papers that cite it, it may be a decent entry point into the subject: P. M. W. Gill, J. A. Pople, Int. J. of Quantum Chem., Vol. 40, (1991) 753.

• Thanks! I had one more question. There are too many integrals to keep them in main memory. As I understand, the integrals are calculated and stored as matrices. Could it be possible to take the indices of the primitives that are duplicate... – Shoubhik R Maiti Nov 23 '20 at 15:02
• ...because the basis set stays the same throughout the calculation. So, it should be possible to identify which combinations of primitives are the same, even before the calculation actually starts, and then skip/replace them by index? – Shoubhik R Maiti Nov 23 '20 at 15:04
• Basically, you are probably right. However, you will run into the same storage problem, I think. – TAR86 Nov 23 '20 at 15:07
• The answer is a bit besides the point. Also, I'm not sure if memory is really an issue for storage in main memory here, since the question is really about computing individual shell quartets that share the same primitive exponents. Although the number of quartets grows like N^4, where N is the number of shells, the present question can be limited to single types of shell quartets, i.e. by fixing the indices of the four centers. – Susi Lehtola Nov 23 '20 at 20:55

## MOLPRO

Continuing Susi's excellent benchmarking for $$\ce{Zn}$$ and $$\ce{Zn_2}$$, MOLPRO 2012 gives for $$\ce{Zn}$$:

RHF-SCF energy:    -1777.84665521
CPU time for INT:      0.42 SEC
CPU time for RHF:      0.01 SEC
Disk used for INT:    12.87 MB


However the RHF energy was converged already at iteration 1, so there wasn't much work to do:

 ITERATION    DDIFF          GRAD             ENERGY        2-EL.EN.            DIPOLE MOMENTS         DIIS   ORB.
1      0.000D+00      0.000D+00     -1777.84665521   1412.514729    0.00000    0.00000    0.00000    0    start
2      0.000D+00      0.328D-06     -1777.84665521   1412.514530    0.00000    0.00000    0.00000    0    orth

Final occupancy:   6   2   2   1   2   1   1   0

!RHF STATE 1.1 Energy              -1777.846655206719


For $$\ce{Zn_2}$$:

RHF-SCF energy:    -3555.64477389
CPU time for INT:      1.85 SEC
CPU time for RHF:      0.03 SEC
Disk used for INT:    28.64 MB


This time the program actually did something for the RHF iterations:

 ITERATION    DDIFF          GRAD             ENERGY        2-EL.EN.            DIPOLE MOMENTS         DIIS   ORB.
1      0.000D+00      0.000D+00     -3555.64052105   3233.546012    0.00000    0.00000    0.00000    0    start
2      0.000D+00      0.302D-02     -3555.64412736   3238.277667    0.00000    0.00000    0.00000    1    diag
3      0.131D-01      0.237D-02     -3555.64463460   3236.375359    0.00000    0.00000    0.00000    2    diag
4      0.373D-02      0.709D-03     -3555.64472506   3236.809227    0.00000    0.00000    0.00000    3    diag
5      0.210D-02      0.294D-03     -3555.64476998   3236.644635    0.00000    0.00000    0.00000    4    diag
6      0.215D-02      0.737D-04     -3555.64477361   3236.670571    0.00000    0.00000    0.00000    5    diag
7      0.548D-03      0.192D-04     -3555.64477389   3236.660433    0.00000    0.00000    0.00000    6    diag
8      0.178D-03      0.250D-05     -3555.64477389   3236.659036    0.00000    0.00000    0.00000    7    diag
9      0.189D-04      0.443D-06     -3555.64477389   3236.659009    0.00000    0.00000    0.00000    8    diag
10      0.373D-05      0.512D-07     -3555.64477389   3236.659037    0.00000    0.00000    0.00000    0    orth

Final occupancy:   8   3   3   1   8   3   3   1

!RHF STATE 1.1 Energy              -3555.644773894570


This was all with only 1 core but to compare it to Susi's benchmarks, we would need to be using the same hardware. I have put my output files in the folder for this question on GitHub, in case anyone wants to try on their hardware (but when the runtimes are so short, comparisons are probably not so valuable: perhaps we need to do this on a bigger basis set so that the programs actually have to take a substantial amount of time calculating things). Also if anyone wants to put their files in the Git repo, I can run them on the same hardware I used to do the above calculations, and can share the results here.