Molecular wavefunctions in critic2

Introduction

Molecular wavefunctions and structures can be loaded in critic2 using the MOLECULE and LOAD keywords. There is no standard file format for the description of a molecular wavefunction and, even worse, wavefunction files whose format is supposedly the same are written by different programs in subtly different ways, leading to the erroneous calculation of molecular properties, and a great deal of confusion. The hope is that these notes serve to clarify some of these issues and also report how critic2 understands the wavefunction files written by various programs.

Relevant Equations for Gaussian Functions

Molecular spinorbitals (MO) are written as linear combination of atomic orbitals (AO): \begin{equation} \psi^\sigma_i({\bf r}) = \sum_j c^\sigma_j \chi_j({\bf r}) \end{equation} Each MO has an associated occupation $f_i^\sigma$, such that the electron density (for instance) is: \begin{equation} \rho({\bf r}) = \sum_\sigma \sum_i^{\rm occ} f_i^\sigma \lvert\psi_i^\sigma({\bf r})\rvert^2 \end{equation} The values of $f_i^\sigma$ go from zero to one but, in the case of a restricted or restricted-open wavefunction, we can have $f_i^\sigma = 2$ and the sum goes over spatial orbitals instead. Wavefunctions written in terms of natural orbitals (for instance, generated by a correlated wavefunction method) have fractional occupations. The molecular orbitals are normalized in such a way that the electron density integrates to the number of electrons.

The most natural way of reporting the molecular orbitals $c_j$ in wavefunction files is to have them combine with AOs that are normalized, such that the overlap matrix has a diagonal of ones. The AOs themselves are written as a linear combination of primitives, which are Gaussian-type functions (GTFs): \begin{equation} \chi_j({\bf r}) = N_{\rm AO}({\bf d},{\bf \alpha},l_x,l_y,l_z) \sum_k d_k^j g(\alpha_k,l_x,l_y,l_z;{\bf r}) \end{equation} where: \begin{equation} g(\alpha_k,l_x,l_y,l_z;{\bf r}) = N_{\rm GTF}(l_x,l_y,l_z,\alpha_k) x^{l_x}y^{l_y}z^{l_z} e^{-\alpha_k r^2} \end{equation} The $l_x$, $l_y$, and $l_z$ quantities are integers greater than or equal to zero, representing the angular momentum quantum number of the GTF. All the primitives that make up one AO share the same angular momentum numbers, but have different exponents. The $N_{\rm GTF}$ factor in front is the normalization factor of the GTF, and it depends on the angular momentum numbers and on the exponent.

We now calculate the normalization constant for a single GTF. The GTF can be written as the product of three one-dimensional Gaussian functions: \begin{equation} g(\alpha_k,l_x,l_y,l_z;{\bf r}) = g(\alpha_k,l_x;x)\times g(\alpha_k,l_y;y)\times g(\alpha_k,l_z;z) \end{equation} where: \begin{equation} g(\alpha,l_x;x) = N_1(\alpha,l_x) x^{l_x} e^{-\alpha x^2} \end{equation} The normalization factor can be derived easily using the integral: \begin{equation} \int_{-\infty}^{\infty} x^{2k} e^{-\alpha x^2} dx = \frac{(2k-1)!!}{(2\alpha)^k}\sqrt{\frac{\pi}{\alpha}} \end{equation} from where we have: \begin{equation} \langle g(\alpha_i,l_x;x) | g(\alpha_i,l_x;x) \rangle = N_1(\alpha,l_x)^2 \frac{(2l_x-1)!!}{(4\alpha)^{l_x}} \sqrt{\frac{\pi}{2\alpha}} \end{equation} \begin{equation} \langle g(\alpha_i,l_x;x) | g(\alpha_j,l_x;x) \rangle = N_1(\alpha_i,l_x)N_1(\alpha_j,l_x) \frac{(2l_x-1)!!}{(2(\alpha_i+\alpha_j))^{l_x}} \sqrt{\frac{\pi}{\alpha_i+\alpha_j}} \end{equation} Making the first expression equal to one, we have the normalization factor: \begin{equation} N_1(\alpha,l_x) = \frac{2^{l_x+1/4} \alpha^{l_x/2+1/4}}{\pi^{1/4} \sqrt{(2l_x-1)!!}} \end{equation} and the overlap between two Gaussian functions with different exponent but same angular momentum is: \begin{equation} \langle g(\alpha_i,l_x;x) | g(\alpha_j,l_x;x) \rangle = 2^{l_x+1/2} \frac{(\alpha_i\alpha_j)^{l_x/2+1/4}}{(\alpha_i+\alpha_j)^{l_x+1/2}} \end{equation}

For the three-dimensional GTFs, we have that the normalization factor is the product of the three one-dimensional $N_{\rm GTF}$: \begin{equation} N_{\rm GTF}(l_x,l_y,l_z,\alpha_k) = \frac{2^{l+3/4} \alpha^{l/2+3/4}}{\pi^{3/4} \sqrt{(2l_x-1)!!(2l_y-1)!!(2l_z-1)!!}} \end{equation} where $l = l_x + l_y + l_z$. The overlap between two GTFs with different exponents but same $l$ and center is: \begin{equation} \langle \alpha_i | \alpha_j \rangle = 2^{l+3/2} \frac{(\alpha_i\alpha_j)^{l/2+3/4}}{(\alpha_i+\alpha_j)^{l+3/2}} \end{equation} With this result, we can now find the normalization factor for a single atomic orbital: \begin{equation} \langle\chi|\chi\rangle = N_{\rm AO}({\bf d},{\bf \alpha},l_x,l_y,l_z)^2 2^{l+3/2} \sum_{ij} d_i d_j \frac{(\alpha_i\alpha_j)^{l/2+3/4}}{(\alpha_i+\alpha_j)^{l+3/2}} = 1 \end{equation} Solving for the normalization constant: \begin{equation} N_{\rm AO}({\bf d},{\bf \alpha},l_x,l_y,l_z) = \left[2^{l+3/2} \sum_{ij} d_i d_j \frac{(\alpha_i\alpha_j)^{l/2+3/4}}{(\alpha_i+\alpha_j)^{l+3/2}} \right]^{-1/2} \end{equation} Alternatively, we can insert the normalization constants of both GTFs in the sum to give: \begin{equation} N_{\rm AO}({\bf d},{\bf \alpha},l_x,l_y,l_z) = \left[ \frac{(2l_x-1)!!(2l_y-1)!!(2l_z-1)!! \pi^{3/2}}{2^l} \sum_{ij} \frac{d_i d_j N_{\rm GTF}^i N_{\rm GTF}^j}{(\alpha_i+\alpha_j)^{l+3/2}} \right]^{-1/2} \end{equation}

Wavefunction Data Structure in Critic2

For the purpose of calculating the electron density and other fields at given points in space, critic2 does not need the atomic orbital information. Instead, the MO coefficients are multiplied by the contraction coefficients and a list of coefficient/exponent pairs is generated. Specifically, critic2 writes a given molecular orbital as a sum over primitive GTFs: \begin{equation} \psi^\sigma_i({\bf r}) = \sum_k C_{ki}^\sigma x^{l^k_x}y^{l^k_y}z^{l^k_z} e^{-\alpha_k r^2} \end{equation} where the $C_{ki}^\sigma$ are a combination of the MO coefficients $c$, the contraction coefficients $d$ and the normalization constants: \begin{equation} C^\sigma_{ki} = c^\sigma_{ki} \times d_k \times N_{\rm AO}({\bf d},{\bf \alpha},l^k_x,l^k_y,l^k_z) \times N_{\rm GTF}(l^k_x,l^k_y,l^k_z,\alpha_k) \end{equation} Note that the same $c^\sigma_{ki}$ is repeated in different $k$s, since the same MO coefficient multiplies all primitive functions in the same AO. Also, note that the normalization constant for the AO depends on the contraction coefficients and exponents of the whole shell, not just primitive $k$.

In practice, the following is the basic information critic2 uses to calculate the molecular scalar fields:

Atomic identifier for the primitives’ centers (icenter).
The “type” of primitive, a combined integer index that represent the three values $l_x$, $l_y$, and $l_z$ (itype).
Primitive exponents (e).
Molecular orbital occupations (occ).
The $C^\sigma_{ki}$ coefficients (cmo).

A lot of additional information like the spins, wavefunction type, etc. is also written down by critic2, but these are the basic variables of the wavefunction type.

Order of Primitive Types Within the Same Shell

Wavefunction file formats that give the AO exponents, AO coefficients, and the MO coefficients do not strictly contain all the information necessary to reproduce the wavefunction. The reason is that the AO exponents and coefficients are given per shell, but there is no information regarding how the AOs within that shell are ordered.

For example, the hydrogen atom in the cc-pVTZ basis set has a d shell with a single primitive. This primitive’s exponent is 1.057 and the coefficient is 1.0. In a given calculation using Cartesian primitives, there are six MO coefficients associated with this shell, corresponding to the AOs with angular parts xx, xy, xz, yy, yz, and zz. To which of this six angular parts do those six coefficients correspond? The wavefunction file gives no indication how this question is answered, and different programs have different orders for their primitives. If the source of the wavefunction file is a closed-source program, this order has to be deduced using model molecules made specifically for this purpose.

The Cartesian primitive order in critic2 is:

p shell: x, y, z.
d shell: xx, yy, zz, xy, xz, yz.
f shell: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz.
g shell: xxxx, yyyy, zzzz, xxxy, xxxz, xyyy, yyyz, xzzz, yzzz, xxyy, xxzz, yyzz, xxyz, xyyz, xyzz.

Most programs do not follow this order. Translating the MO coefficients from the source order to the critic2 order is done using a permutation vector, typtrans. This vector transforms the integer type of the original primitive into the value it has in critic2.

Spherical and Cartesian primitives

The primitive GTFs used in the description of the MOs can also be given as spherical GTFs: \begin{equation} g(\alpha,l,m;{\bf r}) = S_{lm}({\bf r}) e^{-\alpha_k r^2} \end{equation} where $S_{lm}$ is a real solid harmonic function: \begin{equation} S_{lm} = \sqrt{2} {\rm Im}\left[ r^l Y_l^{|m|}\right] \quad {\rm if\ } m < 0
\end{equation} \begin{equation} S_{lm} = Y_l^0 \quad {\rm if\ } m = 0
\end{equation} \begin{equation} S_{lm} = \sqrt{2} {\rm Re}\left[ r^l Y_l^m\right] \quad {\rm if\ } m > 0
\end{equation} where the $Y_l^m$ are spherical harmonics.

The transformation between spherical and Cartesian primitives can be carried out by replacing the $2l+1$ MO coefficients corresponding to the spherical primitives with the ${\binom {3+l-1} l}$ MO coefficients for the Cartesian primitives. Since the $S_{lm}$ functions can be written as a linear combination of Cartesian monomials of degree $l$, the transformation between spherical and Cartesian MO coefficients involves a simple multiplication of the coefficients for a given shell, arranged in vector form, with the transformation matrix between spherical and Cartesian GTFs. This transformation matrix is given in Schlegel and Frisch, Int. J. Quantum Chem. 54 (1995) 83. In particular, the coefficients can be extracted from Table 1 (note the table gives $r^l Y_l^m$ instead of $S_{lm}$). They can also be calculated numerically or analytically from equation 15. The numerical calculation of these coefficients has been implemented in this script.

SP shells

SP-shells are s and p shells that share the same number of primitives and contraction exponents but have different coefficients. They are used in Pople’s basis sets (e.g. 6-31G*). Some wavefunction file formats can give a shell type as SP. When that happens, two sets of contraction coefficients are given that correspond to the s and p shells. If an sp shell is found, critic2 simply unfolds it into the corresponding s and p shells, which are then independent of each other. This is slightly inefficient, but it is by far the simplest solution.

Wavefunction File Formats in Critic2

Current Implementation Status

Program	Format	Works with at least these versions…
Gaussian	wfn	g03(E.01), g09(A.02,B.01,C.01,D.01,E.01), g16(A.03,B.01,C.01)
	wfx	g09(B.01,C.01,D.01,E.01), g16(A.03,B.01,C.01)
	fchk	g03(E.01), g09(A.02,B.01,C.01,D.01,E.01), g16(A.03,B.01,C.01)
psi4	molden	1.1, 1.4, 1.9.1
orca	wfn	4.0
	wfx	4.2
	molden	4.0, 4.2, 5.0.3, 6.1.0
xtb	molden	6.7.1 (GFN1-xTB and GFN2-xTB)
PySCF	molden	2.14.0

The limitations for each of these formats are listed below, along with notes on their particularities. Primitives with angular momentum of h and above are not supported at present, but can be implemented if necessary. Other programs may work with critic2 if they write these formats, but you should test that the results are correct (e.g. by integrating the molecular charge using MOLCALC "$1").

Gaussian Wfn and Wfx File Formats (wfn, wfx)

The Gaussian wfn and wfx file formats contain exactly the information described above, so reading the wavefunction information from them is straightforward (versioning and formatting problems notwithstanding).

Orca Wfn File

The structure of these files seems to be the same as in Gaussian, but early versions of orca are a little finnicky. In particular, version 4.0.2 has the following issues:

There can be lines such as CENTRE ASSIGNMENT that have no information associated with them, and this can throw off the parser in critic2.
g orbitals are not written to orca wfn files.
Unrestricted wavefunctions are not written.

Some or all of these issues seem to have been fixed in version 4.2.

Orca Wfx File

These are written by orca 4.2 but not by orca 4.0.2. They seem to follow the same format as Gaussian and no problems reading these files have been found. Orca seems to refuse writing a wfx file if the calculation contains ECPs.

Gaussian Formatted Checkpoint File (fchk)

The primitive order is:
- p shell: x, y, z.
- d shell: xx, yy, zz, xy, xz, yz.
- f shell: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz.
- g shell: zzzz, yzzz, yyzz, yyyz, yyyy, xzzz, xyzz, xyyz, xyyy, xxzz, xxyz, xxyy, xxxz, xxxy, xxxx.
The fchk file may contain SP shells. These are indicated as shell type -1.
Spherical and Cartesian shells are possible in a fchk file. Cartesian shells are indicated by a positive type (e.g. +2 is a Cartesian d shell) and spherical shells by a negative type (-2 is a spherical d shell).
The contraction coefficients are not multiplied by any normalization constant.

Molden files (molden)

The molden website, where the specifcation of this format was given, seems to be down at present. The following lists what is implemented in critic2 at the moment, and I have no way of checking whether these follow the “official” molden format or not. But each program seems to have its own individual opinion on how to write these files anyway. In the following, the particular version of the molden file written by a program is called its “dialect”. Critic2 implements six dialects: generic, psi4, orca, xtb, PySCF, and ADF (STOs).

The dialect can be given explicitly by the user, either with a keyword after the file name (LOAD file.molden PSI4, ORCA, XTB, or PYSCF) or before it (LOAD MOLDEN_PSI4 file, MOLDEN_ORCA, MOLDEN_XTB, or MOLDEN_PYSCF). This is always the safest option.
If no dialect is given, critic2 guesses it from the contents of the file (not from its extension: .molden and .molden.input are treated the same way). The checks are made in this order:
- ADF (STOs): the file contains an [STO] block.
- orca: the line after the [Title] tag contains the words “created by orca_2mkl”.
- xtb: the [Title] block is empty (the line after [Title] is already the next tag).
- PySCF: the line right after [Molden Format] contains “made by pyscf”.
- psi4: the line right after [Molden Format] is [Atoms] (AU) (psi4 writes neither a title nor a program line).
- generic: none of the above. In this case, critic2 prints a warning suggesting that the dialect be given explicitly.
The psi4 and xtb checks are heuristic and only reflect the files written by the versions in the table above. A file written by another program that happens to look like one of these may be read with the wrong dialect, which is why giving the dialect explicitly is recommended.
The generic dialect follows the usual molden convention: the contraction coefficients are not multiplied by any normalization constant, and the MO coefficients of Cartesian shells refer to Cartesian functions that are normalized individually. The other GTO dialects are deviations from this convention, described below.
The primitive order is:
- p shell: x, y, z.
- d shell: xx, yy, zz, xy, xz, yz.
- f shell: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz.
- g shell: xxxx, yyyy, zzzz, xxxy, xxxz, xyyy, yyyz, xzzz, yzzz, xxyy, xxzz, yyzz, xxyz, xyyz, xyzz.
The spherical functions in a shell are written in the order $m = 0, +1, -1, +2, -2, +3, -3, +4, -4$.
Spherical primitives are possible in a molden file. They are indicated by the “5d”, “7f”, “5d7f”, “5d10f”, and “9g” keywords, which is somewhat confusing. If none of these keywords are present, all primitives are assumed Cartesian.
- “5d” and “5d7f” are equivalent, and they indicate that d and f shells are spherical.
- “5d10f”: the d shells are spherical but the f shells are Cartesian.
- “7f”: the f shells are spherical but the d shells are Cartesian.
- “9g”: the g shells are spherical.
If several of these keywords appear in the same file, the last one that refers to the d and f shells is the one that counts, and “9g” applies regardless of the order. This is the literal reading of the keywords, and it is used in the generic, psi4, and xtb dialects. However, some programs (orca, PySCF) write “5d”, “7f”, and “9g” on separate lines to indicate that d, f, and g shells are all spherical, which in the literal reading would give Cartesian d shells. Orca files are always read as spherical, and in the PySCF dialect the keywords add up instead (see below and issue #87).
ROHF and fractional-occupation wavefunctions are currently not supported by critic2 with molden files.
SP shells are not implemented.
If the file contains several wavefunctions (several [Molden Format] tags), only the first one is read and a warning is printed.

Psi4 Molden File

Critic2 will use the psi4 dialect if no indication is given by the user and the [Atoms] (AU) tag comes right after [Molden Format].
Psi4 does not write SP shells. The SP shells are already unpacked into the corresponding s and p shells in the molden file.
The contraction coefficients are not multiplied by any normalization constant.
The MO coefficients of the Cartesian shells refer to Cartesian functions that all have the normalization constant of the $x^l$, $y^l$, and $z^l$ functions in the same shell, instead of being normalized individually. Critic2 multiplies the MO coefficient of the $x^iy^jz^k$ function ($l = i+j+k$) by: \begin{equation} \sqrt{\frac{(2i-1)!!(2j-1)!!(2k-1)!!}{(2l-1)!!}} \end{equation} For instance, the coefficients of the xy functions in a d shell are divided by $\sqrt{3}$, and those of the xx, yy, and zz functions are unchanged. Spherical shells are not affected. Before this correction, critic2 gave the wrong electron density for psi4 molden files with Cartesian basis sets.
The molecular coordinates in the molden file are different from the input file (the molecule is, at least, translated).
Psi4 does not write h primitives to the molden file (but the file is still generated!).
Psi4 fchk files with Cartesian g shells give an electron density that does not agree with the density calculated by psi4 itself (the corresponding molden file does). The Cartesian d and f fchk files and all spherical fchk files are fine.

Orca Molden File

Critic2 will use the orca molden dialect if no indication is given by the user and the words “created by orca_2mkl” appear under the [TITLE] tag. Orca writes files with the .molden.input extension, but the dialect is detected from the contents, not from the extension.
Orca does not write SP shells. The SP shells are already unpacked into the corresponding s and p shells in the molden file.
Orca writes spherical primitives only. In a file written by orca_2mkl (with the title above), the “5d”, etc. keywords are ignored and all shells are read as spherical. This is also true if the orca dialect is given explicitly by the user. If the orca dialect is forced on a file that was not written by orca_2mkl, the keywords are used as in the generic dialect.
Some MO coefficients of the f and g shells have flipped sign. The spherical functions in a shell are written in the order $m = 0, +1, -1, +2, -2, +3, -3, +4, -4$, and the coefficients with flipped sign are those for $\vert m\vert = 3$ in the f shells and $\vert m\vert \geq 3$ in the g shells. All other coefficients, including those for the f and g functions with $\vert m\vert \leq 2$, have the usual sign. This was verified with orca 5.0.3 and 6.1.0 by comparing with the wfx files generated by orca_2aim.
The contraction coefficients are multiplied by a constant that is different from the normalization constant of a GTF. In particular, if the shell has angular momentum $l$, all contraction coefficients in the molden file have been multiplied by: \begin{equation} N(l,\alpha_k) = \frac{2^{l+3/4} \alpha_k^{l/2+3/4}}{\pi^{3/4} \sqrt{(2l-1)!!}} \end{equation} where $\alpha_k$ are their respective exponents. This is the normalization factor for the AOs with angular terms $x^l$, $y^l$, and $z^l$ within that shell.

Xtb Molden File

Critic2 will use the xtb dialect if no indication is given by the user and the [Title] block is empty. Xtb writes the molden file as molden.input (with the --molden option), but the dialect is detected from the contents, so the file can be renamed.
The contraction coefficients are multiplied by the same constant as in orca files ($N(l,\alpha_k)$, see above).
Xtb writes no “5d”, etc. keywords and its shells are Cartesian. Unlike in orca files, the shells are not forced to be spherical, and the MO coefficients are used as given (there are no sign flips).
Xtb uses a minimal valence basis set, so integrating the electron density gives the number of valence electrons (e.g. 26 for SOCl2 and 58 for Fe(CO)5). The density and the dipole moment calculated by critic2 agree with the values printed by xtb for GFN1-xTB and GFN2-xTB.

PySCF Molden File

Critic2 will use the PySCF dialect if no indication is given by the user and the line after [Molden Format] contains the words “made by pyscf” (PySCF writes its version there). The following was determined from the molden writer in PySCF 2.14.0 (pyscf/tools/molden.py) and verified by comparing the electron density calculated by critic2 with the density calculated by PySCF for spherical and Cartesian basis sets up to g shells, RHF, RKS, and UHF wavefunctions, and a basis set with ECPs.
For spherical basis sets, PySCF writes the “5d”, “7f”, and “9g” keywords on separate lines, which means that the d, f, and g shells are all spherical. In the PySCF dialect, the keywords add up: the d shells are spherical if “5d”, “5d7f”, or “5d10f” appear, the f shells are spherical if “5d”, “5d7f”, or “7f” appear, and the g shells are spherical if “9g” appears. In the literal reading of the other dialects, “7f” after “5d” makes the d shells Cartesian, which is wrong for these files.
For Cartesian basis sets, PySCF writes the “6d”, “10f”, and “15g” keywords, which critic2 ignores (Cartesian is the default). PySCF normalizes each Cartesian function individually before writing the file, as in the generic dialect.
The contraction coefficients are not multiplied by any normalization constant, and the spherical functions are in the usual order ($m = 0, +1, -1, \ldots$). Apart from the treatment of the spherical keywords, the PySCF dialect is the same as the generic dialect.
PySCF does not write h and higher shells by default (they are removed from the basis set before writing the file), so the density from a molden file with such a basis set is approximate.
If the basis set has ECPs, PySCF writes a [core] block with the number of core electrons on each atom. Critic2 ignores this block, so the density is that of the valence electrons only, the same as in PySCF.
ROHF wavefunctions are written by PySCF with occupations 1 and 2, and they are not supported (see above).

ADF Molden File (with STOs)

Critic2 will use the ADF molden dialect if the [STO] tag appears in it. Critic2 can handle the STO wavefunction and calculate all the usual properties. Naturally, the equations above do not apply.
Since the atomic orbitals are simple exponentials, there are no contractions and the angular momentum of each AO is indicated explicitly. Therefore, reading this input file is a lot more straightforward than the GTF counterparts.