13. Self-Consistent Field Methods#
👩🏽💻 Assignment: Complete the Below Notebook#
To do this, complete the code blocks indicated by ### START YOUR CODE HERE and ### END YOUR CODE HERE .
The detailed steps you’ll follow are:
Read the theoretical background on the self-consistent field (SCF) method.
Implement the “core Hamiltonian” guess.
Implement the addition of an electric field to the Hamiltonian and compute the dipole moment and polarizability.
Motivation of Mean-Field (Average Potential) Methods#
In atomic units, the electronic Hamiltonian for an \(N\)-electron molecule, with \(P\) atoms with atomic numbers \(Z_A\) at positions \(\mathbf{R}_A\), is given by
\[\hat{H} = -\sum_{i=1}^N \frac{1}{2} \nabla_i^2 - \sum_{i=1}^N \sum_{A=1}^P \frac{Z_A}{|\mathbf{r}_i - \mathbf{R}_A|} + \sum_{i<j} \frac{1}{|\mathbf{r}_i - \mathbf{r}_j|}\]
This Hamiltonian cannot be solved by separation of variables because the electrons are coupled through the electron-electron repulsion term. Mean-field, or average potential, approximations decouple electrons by treating the effect of all other electrons on a given electron as an average potential (operator), \(\hat{w}(\mathbf{r})\). (I.e., individual electrons feel only the effects of the “cloud” of other electrons.) With this revision, the electronic Hamiltonian becomes,
\[\hat{H}_{\text{avg}} = -\sum_{i=1}^N \frac{1}{2} \nabla_i^2 + \sum_{i=1}^N v(\mathbf{r}_i) + \sum_{i=1}^N \hat{w}(\mathbf{r}_i)$
\]
where the external potential is defined as
\[v(\mathbf{r}) = -\sum_{A=1}^P \frac{Z_A}{|\mathbf{r} - \mathbf{R}_A|}\]
The electronic Schrödinger equation for the mean-field Hamiltonian is
\[\hat{H}_{\text{avg}} \Phi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N) = E \Phi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)\]
which can be solved by separation of variables, yielding a set of one-electron equations,
\[\left[-\frac{1}{2} \nabla^2 + v(\mathbf{r}) + \hat{w}(\mathbf{r})\right] \phi_p(\mathbf{r}) = \epsilon_p \phi_p(\mathbf{r})\]
where \(\phi_p(\mathbf{r})\) is the \(p\)-th molecular orbital. However, the average potential \(\hat{w}\) depends on the wavefunction \(\Phi\), which is a Slater determinant of the molecular orbitals \(\phi_p\). Making this dependence explicit, we have,
\[\left[-\frac{1}{2} \nabla^2 + v(\mathbf{r}) + \hat{w}[\{\phi_q\}](\mathbf{r})\right] \phi_p(\mathbf{r}) = \epsilon_p \phi_p(\mathbf{r})\]
This is a nonlinear Schrödinger equation because the operator \(\hat{w}\) depends on the orbitals \(\{\phi_q\}\), which are the solutions to the equation. Since we cannot write the Hamiltonian without knowing the orbitals, and cannot know the orbitals without writing the Hamiltonian, we have a circular problem. The simplest self-consistent-field (SCF) method solves this problem by starting from an initial guess for the orbitals, constructing the Hamiltonian, solving for the orbitals, and then updating the Hamiltonian using those orbitals. This procedure is iterated until the orbitals used to construct the Hamiltonian are the same as the orbitals obtained from solving the Schrödinger equation, at which point we have a self-consistent solution.
Average Potential Models#
Hartree Model#
The Hartree model is the simplest average potential model. It approximates the average potential as the electrostatic potential generated by the total electron density, \(\rho(\mathbf{r})\). Therefore, the Hartree average potential is given by $\(w_{\text{H}}(\mathbf{r}) = \int \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}'\)\( where \)\rho(\mathbf{r})$ is the total electron density. The Hartree model does not account for the fact that an electron does not interact with itself, leading to an overestimation of the electron-electron repulsion energy. This is known as the self-interaction error.
Fermi-Amaldi Model#
The Fermi-Amaldi model is the simplest way to impose a self-interaction correction. It accounts for the fact that an \(\alpha\)-spin electron interacts with all of the \(\beta\)-spin electrons but only \(N_\alpha - 1\) of the \(\alpha\)-spin electrons, and vice versa for \(\beta\)-spin electrons. Therefore, the Fermi-Amaldi average potential is given by $\( \begin{align*} w_{\text{FA};\alpha}(\mathbf{r}) &= \frac{N_\alpha-1}{N_\alpha} \int \frac{\rho_\alpha(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' + \int \frac{\rho_\beta(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}'\\ w_{\text{FA};\beta}(\mathbf{r}) &= \frac{N_\beta-1}{N_\beta} \int \frac{\rho_\beta(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' + \int \frac{\rho_\alpha(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} d\mathbf{r}' \end{align*} \)\( where \)\rho_\alpha(\mathbf{r})\( and \)\rho_\beta(\mathbf{r})\( are the electron densities for \)\alpha\( and \)\beta$ spins, respectively.
Kohn-Sham Density Functional Theory#
Kohn-Sham density functional theory (KS-DFT) is a more sophisticated average potential model that includes an exchange-correlation potential, \(v_{\text{xc}}(\mathbf{r})\), to account for the effects of electron exchange and correlation. The KS-DFT average potential is given by $\(\hat{w}_{\text{KS}}(\mathbf{r}) = w_{\text{H}}(\mathbf{r}) + v_{\text{xc}}(\mathbf{r})\)\( where \)v_{\text{xc}}(\mathbf{r})\( is the exchange-correlation potential, which is a functional of the electron (spin) densities \)\rho_\alpha(\mathbf{r})\( and \)\rho_\beta(\mathbf{r})\(. The exact form of \)v_{\text{xc}}$ is unknown, and various approximations have been developed to model it, such as the local density approximation (LDA) and the generalized gradient approximation (GGA).
The simplest—and one of the first—exchange-correlation functionals is the local spin density approximation (LSDA), which is a simple scaling of the exchange energy from the homogeneous electron gas. If only exchange effects are included, then one has $\(v_{\text{xc};\alpha}(\mathbf{r}) = -\left(\frac{6}{\pi}\right)^{1/3} \rho_\alpha(\mathbf{r})^{1/3}\)\( \)\(v_{\text{xc};\beta}(\mathbf{r}) = -\left(\frac{6}{\pi}\right)^{1/3} \rho_\beta(\mathbf{r})^{1/3}\)$
Hartree-Fock Theory#
Hartree-Fock (HF) theory is an average potential model that includes an exchange potential, \(\hat{K}\), to account for the effects of electron exchange. The HF average potential is given by $\(\hat{w}_{\text{HF}}(\mathbf{r}) = w_{\text{H}}(\mathbf{r}) + \hat{k}(\mathbf{r})\)\( where \)\(\hat{k}(\mathbf{r}) \phi_p(\mathbf{r}) = -\sum_{q=1}^N \delta_{\sigma_q \sigma_p}\int \frac{\phi_q^*(\mathbf{r}') \phi_p(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}d\mathbf{r}' \phi_q(\mathbf{r}) \)$ The exchange operator, accounting for the Pauli exclusion principle, only acts on orbitals with the same spin.
Initial Guess#
In the separated-atom limit, a molecule’s density is the sum of the isolated atomic densities. Even at (near) equilibrium geometries, this promolecular approximation is often a good starting point. When we do this, we say we are using the Harris functional as an initial guess. The Harris functional can be used with any pure density functional (including the Hartree, Fermi-Amaldi, and LSDA).
An easier alternative is the core Hamiltonian guess, which is the solution to the one-electron Schrödinger equation with only the kinetic energy and nuclear attraction terms. This guess isn’t as accurate, but it is easy to compute and for systems with weakly bound electrons, it sometimes performs better than the (raw) Harris functional guess. (An even better guess, however, is often the Harris functional guess but using the electron density of the molecular cation.)
External Packages Needed#
We will use several packages to construct the program.
AtomDBis a database of atomic properties. It can be used to generate the promolecular density for the initial guess.IODatais a package for reading and writing data. We will use it to read the molecular geometry from an XYZ file.GBasisis a package for evaluating Gaussian basis functions and integrals related to them.QC-Gridis a package for evaluating numerical integrals (and derivatives).QC-BFitis a package for fitting functions (like the electron density) to Gaussian basis sets. We use it to build fast (approximate) models for the Coulomb potential.basis_set_exchangedownloads Gaussian basis sets.
## Load Packages
!pip install qc-atomdb
!pip install qc-iodata
!pip install qc-gbasis
!pip install qc-grid
!pip install basis_set_exchange
!pip install qc-bfit
Requirement already satisfied: qc-atomdb in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (0.0.2.post5)
Requirement already satisfied: numpy>=1.16 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (2.4.4)
Requirement already satisfied: scipy>=1.4 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (1.17.1)
Requirement already satisfied: msgpack>=1.0.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (1.1.2)
Requirement already satisfied: msgpack-numpy>=0.4.8 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (0.4.8)
Requirement already satisfied: h5py>=3.6.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (3.16.0)
Requirement already satisfied: importlib-resources>=3.0.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (6.5.2)
Requirement already satisfied: pooch>=1.8.1 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-atomdb) (1.9.0)
Requirement already satisfied: platformdirs>=2.5.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from pooch>=1.8.1->qc-atomdb) (4.9.4)
Requirement already satisfied: packaging>=20.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from pooch>=1.8.1->qc-atomdb) (26.0)
Requirement already satisfied: requests>=2.19.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from pooch>=1.8.1->qc-atomdb) (2.33.1)
Requirement already satisfied: charset_normalizer<4,>=2 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from requests>=2.19.0->pooch>=1.8.1->qc-atomdb) (3.4.7)
Requirement already satisfied: idna<4,>=2.5 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from requests>=2.19.0->pooch>=1.8.1->qc-atomdb) (3.11)
Requirement already satisfied: urllib3<3,>=1.26 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from requests>=2.19.0->pooch>=1.8.1->qc-atomdb) (2.6.3)
Requirement already satisfied: certifi>=2023.5.7 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from requests>=2.19.0->pooch>=1.8.1->qc-atomdb) (2026.2.25)
Requirement already satisfied: qc-iodata in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (1.0.1)
Requirement already satisfied: numpy>=1.26.4 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-iodata) (2.4.4)
Requirement already satisfied: scipy>=1.13.1 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-iodata) (1.17.1)
Requirement already satisfied: attrs>=21.3.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-iodata) (26.1.0)
Collecting qc-gbasis
Downloading qc_gbasis-0.1.0-py3-none-any.whl.metadata (7.4 kB)
Requirement already satisfied: scipy>=1.11.1 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-gbasis) (1.17.1)
Requirement already satisfied: importlib-resources in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-gbasis) (6.5.2)
Requirement already satisfied: sympy in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-gbasis) (1.14.0)
Requirement already satisfied: numpy>=1.22 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-gbasis) (2.4.4)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from sympy->qc-gbasis) (1.3.0)
Downloading qc_gbasis-0.1.0-py3-none-any.whl (97 kB)
Installing collected packages: qc-gbasis
Successfully installed qc-gbasis-0.1.0
Collecting qc-grid
Downloading qc_grid-0.0.9.post1-py3-none-any.whl.metadata (4.5 kB)
Requirement already satisfied: numpy>=1.16 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-grid) (2.4.4)
Collecting pytest>=8.0.0 (from qc-grid)
Downloading pytest-9.0.2-py3-none-any.whl.metadata (7.6 kB)
Requirement already satisfied: scipy>=1.4 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-grid) (1.17.1)
Requirement already satisfied: importlib_resources in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-grid) (6.5.2)
Requirement already satisfied: sympy in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from qc-grid) (1.14.0)
Collecting iniconfig>=1.0.1 (from pytest>=8.0.0->qc-grid)
Downloading iniconfig-2.3.0-py3-none-any.whl.metadata (2.5 kB)
Requirement already satisfied: packaging>=22 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from pytest>=8.0.0->qc-grid) (26.0)
Collecting pluggy<2,>=1.5 (from pytest>=8.0.0->qc-grid)
Downloading pluggy-1.6.0-py3-none-any.whl.metadata (4.8 kB)
Requirement already satisfied: pygments>=2.7.2 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from pytest>=8.0.0->qc-grid) (2.20.0)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from sympy->qc-grid) (1.3.0)
Downloading qc_grid-0.0.9.post1-py3-none-any.whl (66.1 MB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/66.1 MB ? eta -:--:--
━━━━━━━━╺━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 13.9/66.1 MB 76.2 MB/s eta 0:00:01
━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━ 23.1/66.1 MB 60.1 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━╸━━━━━━━━━━━━━━━━━━━ 34.6/66.1 MB 64.1 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━ 48.2/66.1 MB 64.5 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━━ 54.5/66.1 MB 56.7 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━ 61.9/66.1 MB 51.9 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸ 66.1/66.1 MB 48.7 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 66.1/66.1 MB 46.3 MB/s 0:00:01
?25hDownloading pytest-9.0.2-py3-none-any.whl (374 kB)
Downloading pluggy-1.6.0-py3-none-any.whl (20 kB)
Downloading iniconfig-2.3.0-py3-none-any.whl (7.5 kB)
Installing collected packages: pluggy, iniconfig, pytest, qc-grid
?25l
━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━ 2/4 [pytest]
━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━ 2/4 [pytest]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [qc-grid]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [qc-grid]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [qc-grid]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4/4 [qc-grid]
?25h
Successfully installed iniconfig-2.3.0 pluggy-1.6.0 pytest-9.0.2 qc-grid-0.0.9.post1
Collecting basis_set_exchange
Downloading basis_set_exchange-0.12-py3-none-any.whl.metadata (6.7 kB)
Requirement already satisfied: jsonschema in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from basis_set_exchange) (4.26.0)
Collecting argcomplete (from basis_set_exchange)
Downloading argcomplete-3.6.3-py3-none-any.whl.metadata (16 kB)
Collecting regex (from basis_set_exchange)
Downloading regex-2026.4.4-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl.metadata (40 kB)
Collecting unidecode (from basis_set_exchange)
Downloading Unidecode-1.4.0-py3-none-any.whl.metadata (13 kB)
Requirement already satisfied: attrs>=22.2.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from jsonschema->basis_set_exchange) (26.1.0)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from jsonschema->basis_set_exchange) (2025.9.1)
Requirement already satisfied: referencing>=0.28.4 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from jsonschema->basis_set_exchange) (0.37.0)
Requirement already satisfied: rpds-py>=0.25.0 in /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages (from jsonschema->basis_set_exchange) (0.30.0)
Downloading basis_set_exchange-0.12-py3-none-any.whl (40.1 MB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/40.1 MB ? eta -:--:--
━━━━━━━╸━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.9/40.1 MB 42.2 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━ 28.3/40.1 MB 88.2 MB/s eta 0:00:01
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 40.1/40.1 MB 81.4 MB/s 0:00:00
?25hDownloading argcomplete-3.6.3-py3-none-any.whl (43 kB)
Downloading regex-2026.4.4-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl (800 kB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/800.6 kB ? eta -:--:--
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 800.6/800.6 kB 176.8 MB/s 0:00:00
?25hDownloading Unidecode-1.4.0-py3-none-any.whl (235 kB)
Installing collected packages: unidecode, regex, argcomplete, basis_set_exchange
?25l
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0/4 [unidecode]
━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━━━━━━━━━━━ 2/4 [argcomplete]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╺━━━━━━━━━ 3/4 [basis_set_exchange]
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4/4 [basis_set_exchange]
?25h
Successfully installed argcomplete-3.6.3 basis_set_exchange-0.12 regex-2026.4.4 unidecode-1.4.0
Collecting qc-bfit
Downloading qc_BFit-0.1.2-py3-none-any.whl.metadata (49 kB)
Downloading qc_BFit-0.1.2-py3-none-any.whl (561 kB)
?25l ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/561.1 kB ? eta -:--:--
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 561.1/561.1 kB 13.1 MB/s 0:00:00
?25h
Installing collected packages: qc-bfit
Successfully installed qc-bfit-0.1.2
# Load the Water molecule geometry from an XYZ file (using IOData)
from urllib.request import urlretrieve
from iodata import load_one
# Load the Water molecule geometry from an XYZ file (using IOData)
from urllib.request import urlretrieve
from iodata import load_one
url = "https://raw.githubusercontent.com/theochem/iodata/main/iodata/test/data/water.xyz"
urlretrieve(url, "water.xyz")
h2o = load_one("water.xyz")
print(h2o.atnums)
print(h2o.atcoords)
[1 8 1]
[[ 1.48123726 -0.93019123 -0. ]
[-0. 0.11720081 -0. ]
[-1.48123726 -0.93019123 -0. ]]
# Select a basis set (using basis_set_exchange)
import iodata.periodic
import iodata.utils
from gbasis.parsers import parse_gbs, make_contractions
import basis_set_exchange as bse
import tempfile, iodata
def get_ao_basis_from_bse(atnums, atcoords, basis_name, coord_types="spherical"):
r"""Fetch the atomic orbital basis set from the online repository.
The basis will be downloaded from the basis set exchange repository and returned in the gaussian94
format. The basis set exchange repository is available at:
https://www.basissetexchange.org/
Parameters
----------
atnums : list of int
Atomic numbers of the atoms in the molecule.
atcoords : numpy.ndarray[float, shape=(n_atoms, 3)]
Coordinates of the atoms in the molecule.
basis_name : str
Name of the atomic orbital basis set to fetch.
coord_types : {"cartesian"/"c", list/tuple of "cartesian"/"c" or "spherical"/"p", "spherical"/"p"}
Types of the coordinate system for the contractions.
If "cartesian" or "c", then all of the contractions are treated as Cartesian contractions.
If "spherical" or "p", then all of the contractions are treated as spherical contractions.
If list/tuple, then each entry must be a "cartesian" (or "c") or "spherical" (or "p") to
specify the coordinate type of each `GeneralizedContractionShell` instance.
Default value is "spherical".
"""
local_atnums = [int(i) for i in atnums]
atoms = [iodata.periodic.num2sym[i] for i in local_atnums]
gbs_basis =bse.api.get_basis(name = basis_name, elements=local_atnums, fmt="gaussian94")
# save basis to tmp file
with tempfile.NamedTemporaryFile(delete=False, mode='w', suffix='.gbs') as tmp_data:
tmp_data.write(gbs_basis)
fname = tmp_data.name
basis_dict = parse_gbs(fname)
ao_basis = make_contractions(basis_dict, atoms, atcoords, coord_types)
return ao_basis
bse.get_all_basis_names()
['2ZaP',
'2ZaPa-NR',
'2ZaPa-NR-CV',
'3-21G',
'3ZaP',
'3ZaPa-NR',
'3ZaPa-NR-CV',
'4-31G',
'4ZaP',
'4ZaPa-NR',
'4ZaPa-NR-CV',
'5-21G',
'5ZaP',
'5ZaPa-NR',
'5ZaPa-NR-CV',
'6-21G',
'6-31++G',
'6-31++G*',
'6-31++G**',
'6-31++G**-J',
'6-31+G',
'6-31+G*',
'6-31+G**',
'6-31+G*-J',
'6-311++G',
'6-311++G(2d,2p)',
'6-311++G(3df,3pd)',
'6-311++G*',
'6-311++G**',
'6-311++G**-J',
'6-311+G',
'6-311+G(2d,p)',
'6-311+G*',
'6-311+G**',
'6-311+G*-J',
'6-311G',
'6-311G(2df,2pd)',
'6-311G(d,p)',
'6-311G*',
'6-311G**',
'6-311G**-RIFIT',
'6-311G-J',
'6-311xxG(d,p)',
'6-31G',
'6-31G(2df,p)',
'6-31G(3df,3pd)',
'6-31G(d,p)',
'6-31G*',
'6-31G**',
'6-31G**-RIFIT',
'6-31G*-Blaudeau',
'6-31G-Blaudeau',
'6-31G-J',
'6ZaP',
'6ZaPa-NR',
'7ZaPa-NR',
'AHGBS-5',
'AHGBS-7',
'AHGBS-9',
'AHGBSP1-5',
'AHGBSP1-7',
'AHGBSP1-9',
'AHGBSP2-5',
'AHGBSP2-7',
'AHGBSP2-9',
'AHGBSP3-5',
'AHGBSP3-7',
'AHGBSP3-9',
'ANO-DK3',
'ANO-R',
'ANO-R0',
'ANO-R1',
'ANO-R2',
'ANO-R3',
'ANO-RCC',
'ANO-RCC-MB',
'ANO-RCC-VDZ',
'ANO-RCC-VDZP',
'ANO-RCC-VQZP',
'ANO-RCC-VTZ',
'ANO-RCC-VTZP',
'ANO-VT-DZ',
'ANO-VT-QZ',
'ANO-VT-TZ',
'ATZP-ZORA',
'Ahlrichs TZV',
'Ahlrichs VDZ',
'Ahlrichs VTZ',
'Ahlrichs pVDZ',
'CADPAC-TZ2P',
'CRENBL',
'CRENBL ECP',
'CRENBS',
'CRENBS ECP',
'Cologne DKH2',
'DFO+-NRLMOL',
'DFO-1',
'DFO-1-BHS',
'DFO-2',
'DFO-NRLMOL',
'DZ (Dunning-Hay)',
'DZ + Double Rydberg (Dunning-Hay)',
'DZ + Rydberg (Dunning-Hay)',
'DZP (Dunning-Hay)',
'DZP + Diffuse (Dunning-Hay)',
'DZP + Rydberg (Dunning-Hay)',
'FANO-5Z',
'FANO-6Z',
'FANO-DZ',
'FANO-QZ',
'FANO-TZ',
'Grimme vDZP',
'HGBS-5',
'HGBS-7',
'HGBS-9',
'HGBSP1-5',
'HGBSP1-7',
'HGBSP1-9',
'HGBSP2-5',
'HGBSP2-7',
'HGBSP2-9',
'HGBSP3-5',
'HGBSP3-7',
'HGBSP3-9',
'IGLO-II',
'IGLO-III',
'Koga unpolarized',
'LANL08',
'LANL08(d)',
'LANL08(f)',
'LANL08+',
'LANL2DZ',
'LANL2DZ ECP',
'LANL2DZdp',
'LANL2TZ',
'LANL2TZ(f)',
'LANL2TZ+',
'MIDI',
'MIDI!',
'MIDIX',
'MINI',
'NASA Ames ANO',
'NASA Ames ANO2',
'NASA Ames cc-pCV5Z',
'NASA Ames cc-pCVQZ',
'NASA Ames cc-pCVTZ',
'NASA Ames cc-pV5Z',
'NASA Ames cc-pVQZ',
'NASA Ames cc-pVTZ',
'NLO-V',
'NMR-DKH (TZ2P)',
'ORP',
'PAW-L05',
'PAW-L1',
'PAW-L1-contracted',
'PAW-L2',
'PAW-L2-contracted',
'PB4-D',
'PB4-F1',
'PB4-F2',
'PB5-D',
'PB5-F',
'PB5-G',
'PB6-D',
'PB6-F',
'PB6-G',
'PB6-H',
'Partridge Uncontracted 1',
'Partridge Uncontracted 2',
'Partridge Uncontracted 3',
'Partridge Uncontracted 4',
'PsX-DZ',
'PsX-QZ',
'PsX-TZ',
'Pt - mDZP',
'Roos Augmented Double Zeta ANO',
'Roos Augmented Triple Zeta ANO',
'SARC-DKH2',
'SARC-ZORA',
'SARC2-QZV-DKH2',
'SARC2-QZV-DKH2-JKFIT',
'SARC2-QZV-ZORA',
'SARC2-QZV-ZORA-JKFIT',
'SARC2-QZVP-DKH2',
'SARC2-QZVP-DKH2-JKFIT',
'SARC2-QZVP-ZORA',
'SARC2-QZVP-ZORA-JKFIT',
'SBKJC Polarized (p,2d) - LFK',
'SBKJC-ECP',
'SBKJC-VDZ',
'SBO4-DZ(d)-3G',
'SBO4-DZ(d,p)-3G',
'SBO4-SZ-3G',
'STO-2G',
'STO-3G',
'STO-3G*',
'STO-4G',
'STO-5G',
'STO-6G',
'SV (Dunning-Hay)',
'SV + Double Rydberg (Dunning-Hay)',
'SV + Rydberg (Dunning-Hay)',
'SVP (Dunning-Hay)',
'SVP + Diffuse (Dunning-Hay)',
'SVP + Diffuse + Rydberg (Dunning-Hay)',
'SVP + Rydberg (Dunning-Hay)',
'Sadlej pVTZ',
'Sadlej+',
'Sapporo-DKH3-DZP',
'Sapporo-DKH3-DZP-2012',
'Sapporo-DKH3-DZP-2012-diffuse',
'Sapporo-DKH3-DZP-diffuse',
'Sapporo-DKH3-QZP',
'Sapporo-DKH3-QZP-2012',
'Sapporo-DKH3-QZP-2012-diffuse',
'Sapporo-DKH3-QZP-diffuse',
'Sapporo-DKH3-TZP',
'Sapporo-DKH3-TZP-2012',
'Sapporo-DKH3-TZP-2012-diffuse',
'Sapporo-DKH3-TZP-diffuse',
'Sapporo-DZP',
'Sapporo-DZP-2012',
'Sapporo-DZP-2012-diffuse',
'Sapporo-DZP-diffuse',
'Sapporo-QZP',
'Sapporo-QZP-2012',
'Sapporo-QZP-2012-diffuse',
'Sapporo-QZP-diffuse',
'Sapporo-TZP',
'Sapporo-TZP-2012',
'Sapporo-TZP-2012-diffuse',
'Sapporo-TZP-diffuse',
'Scaled MINI',
'Stuttgart RLC',
'Stuttgart RLC ECP',
'Stuttgart RSC 1997',
'Stuttgart RSC 1997 ECP',
'Stuttgart RSC ANO',
'Stuttgart RSC Segmented + ECP',
'TZ (Dunning-Hay)',
'TZP-ZORA',
'UGBS',
'WTBS',
'Wachters+f',
'acv2z-J',
'acv3z-J',
'acv4z-J',
'admm-1',
'admm-2',
'admm-3',
'ano-pV5Z',
'ano-pVDZ',
'ano-pVQZ',
'ano-pVTZ',
'apr-cc-pV(Q+d)Z',
'asigmaDZ',
'asigmaQZ',
'asigmaTZ',
'aug-admm-1',
'aug-admm-2',
'aug-admm-3',
'aug-ano-pV5Z',
'aug-ano-pVDZ',
'aug-ano-pVQZ',
'aug-ano-pVTZ',
'aug-cc-pCV5Z',
'aug-cc-pCVDZ',
'aug-cc-pCVDZ-DK',
'aug-cc-pCVQZ',
'aug-cc-pCVQZ-DK',
'aug-cc-pCVTZ',
'aug-cc-pCVTZ-DK',
'aug-cc-pV(5+d)Z',
'aug-cc-pV(D+d)Z',
'aug-cc-pV(Q+d)Z',
'aug-cc-pV(T+d)Z',
'aug-cc-pV5Z',
'aug-cc-pV5Z-DK',
'aug-cc-pV5Z-OPTRI',
'aug-cc-pV5Z-PP',
'aug-cc-pV5Z-PP-OPTRI',
'aug-cc-pV5Z-PP-RIFIT',
'aug-cc-pV5Z-RIFIT',
'aug-cc-pV6Z',
'aug-cc-pV6Z-RIFIT',
'aug-cc-pV7Z',
'aug-cc-pVDZ',
'aug-cc-pVDZ-DK',
'aug-cc-pVDZ-DK3',
'aug-cc-pVDZ-OPTRI',
'aug-cc-pVDZ-PP',
'aug-cc-pVDZ-PP-OPTRI',
'aug-cc-pVDZ-PP-RIFIT',
'aug-cc-pVDZ-RIFIT',
'aug-cc-pVDZ-X2C',
'aug-cc-pVQZ',
'aug-cc-pVQZ-DK',
'aug-cc-pVQZ-DK3',
'aug-cc-pVQZ-OPTRI',
'aug-cc-pVQZ-PP',
'aug-cc-pVQZ-PP-OPTRI',
'aug-cc-pVQZ-PP-RIFIT',
'aug-cc-pVQZ-RIFIT',
'aug-cc-pVQZ-X2C',
'aug-cc-pVTZ',
'aug-cc-pVTZ-DK',
'aug-cc-pVTZ-DK3',
'aug-cc-pVTZ-J',
'aug-cc-pVTZ-OPTRI',
'aug-cc-pVTZ-PP',
'aug-cc-pVTZ-PP-OPTRI',
'aug-cc-pVTZ-PP-RIFIT',
'aug-cc-pVTZ-RIFIT',
'aug-cc-pVTZ-X2C',
'aug-cc-pwCV5Z',
'aug-cc-pwCV5Z-DK',
'aug-cc-pwCV5Z-PP',
'aug-cc-pwCV5Z-PP-OPTRI',
'aug-cc-pwCV5Z-PP-RIFIT',
'aug-cc-pwCV5Z-RIFIT',
'aug-cc-pwCVDZ',
'aug-cc-pwCVDZ-DK3',
'aug-cc-pwCVDZ-PP',
'aug-cc-pwCVDZ-PP-OPTRI',
'aug-cc-pwCVDZ-PP-RIFIT',
'aug-cc-pwCVDZ-RIFIT',
'aug-cc-pwCVDZ-X2C',
'aug-cc-pwCVQZ',
'aug-cc-pwCVQZ-DK',
'aug-cc-pwCVQZ-DK3',
'aug-cc-pwCVQZ-PP',
'aug-cc-pwCVQZ-PP-OPTRI',
'aug-cc-pwCVQZ-PP-RIFIT',
'aug-cc-pwCVQZ-RIFIT',
'aug-cc-pwCVQZ-X2C',
'aug-cc-pwCVTZ',
'aug-cc-pwCVTZ-DK',
'aug-cc-pwCVTZ-DK3',
'aug-cc-pwCVTZ-PP',
'aug-cc-pwCVTZ-PP-OPTRI',
'aug-cc-pwCVTZ-PP-RIFIT',
'aug-cc-pwCVTZ-RIFIT',
'aug-cc-pwCVTZ-X2C',
'aug-ccX-5Z',
'aug-ccX-DZ',
'aug-ccX-QZ',
'aug-ccX-TZ',
'aug-mcc-pV5Z',
'aug-mcc-pV6Z',
'aug-mcc-pV7Z',
'aug-mcc-pV8Z',
'aug-mcc-pVQZ',
'aug-mcc-pVTZ',
'aug-pV7Z',
'aug-pc-0',
'aug-pc-1',
'aug-pc-2',
'aug-pc-3',
'aug-pc-4',
'aug-pcH-1',
'aug-pcH-2',
'aug-pcH-3',
'aug-pcH-4',
'aug-pcJ-0',
'aug-pcJ-0_2006',
'aug-pcJ-1',
'aug-pcJ-1_2006',
'aug-pcJ-2',
'aug-pcJ-2_2006',
'aug-pcJ-3',
'aug-pcJ-3_2006',
'aug-pcJ-4',
'aug-pcJ-4_2006',
'aug-pcS-0',
'aug-pcS-1',
'aug-pcS-2',
'aug-pcS-3',
'aug-pcS-4',
'aug-pcSseg-0',
'aug-pcSseg-1',
'aug-pcSseg-2',
'aug-pcSseg-3',
'aug-pcSseg-4',
'aug-pcX-1',
'aug-pcX-2',
'aug-pcX-3',
'aug-pcX-4',
'aug-pcseg-0',
'aug-pcseg-1',
'aug-pcseg-2',
'aug-pcseg-3',
'aug-pcseg-4',
'aug-seg-cc-pV5Z-PP',
'aug-seg-cc-pVDZ-PP',
'aug-seg-cc-pVQZ-PP',
'aug-seg-cc-pVTZ-PP',
'aug-seg-cc-pwCV5Z-PP',
'aug-seg-cc-pwCVDZ-PP',
'aug-seg-cc-pwCVQZ-PP',
'aug-seg-cc-pwCVTZ-PP',
'binning 641',
'binning 641(d)',
'binning 641(df)',
'binning 641+',
'binning 641+(d)',
'binning 641+(df)',
'binning 962',
'binning 962(d)',
'binning 962(df)',
'binning 962+',
'binning 962+(d)',
'binning 962+(df)',
'cc-pCV5Z',
'cc-pCVDZ',
'cc-pCVDZ-DK',
'cc-pCVDZ-F12',
'cc-pCVDZ-F12-OPTRI',
'cc-pCVDZ-F12-RIFIT',
'cc-pCVQZ',
'cc-pCVQZ-DK',
'cc-pCVQZ-F12',
'cc-pCVQZ-F12-OPTRI',
'cc-pCVQZ-F12-RIFIT',
'cc-pCVTZ',
'cc-pCVTZ-DK',
'cc-pCVTZ-F12',
'cc-pCVTZ-F12-OPTRI',
'cc-pCVTZ-F12-RIFIT',
'cc-pV(5+d)Z',
'cc-pV(D+d)Z',
'cc-pV(Q+d)Z',
'cc-pV(T+d)Z',
'cc-pV5Z',
'cc-pV5Z(fi/sf/fw)',
'cc-pV5Z(fi/sf/lc)',
'cc-pV5Z(fi/sf/sc)',
'cc-pV5Z(pt/sf/fw)',
'cc-pV5Z(pt/sf/lc)',
'cc-pV5Z(pt/sf/sc)',
'cc-pV5Z-DK',
'cc-pV5Z-F12',
'cc-pV5Z-F12(rev2)',
'cc-pV5Z-JKFIT',
'cc-pV5Z-PP',
'cc-pV5Z-PP-RIFIT',
'cc-pV5Z-RIFIT',
'cc-pV6Z',
'cc-pV6Z-RIFIT',
'cc-pV8Z',
'cc-pV9Z',
'cc-pVDZ',
'cc-pVDZ(fi/sf/fw)',
'cc-pVDZ(fi/sf/lc)',
'cc-pVDZ(fi/sf/sc)',
'cc-pVDZ(pt/sf/fw)',
'cc-pVDZ(pt/sf/lc)',
'cc-pVDZ(pt/sf/sc)',
'cc-pVDZ(seg-opt)',
'cc-pVDZ-DK',
'cc-pVDZ-DK3',
'cc-pVDZ-F12',
'cc-pVDZ-F12(rev2)',
'cc-pVDZ-F12-OPTRI',
'cc-pVDZ-F12-OPTRI+',
'cc-pVDZ-PP',
'cc-pVDZ-PP-RIFIT',
'cc-pVDZ-RIFIT',
'cc-pVDZ-X2C',
'cc-pVQZ',
'cc-pVQZ(fi/sf/fw)',
'cc-pVQZ(fi/sf/lc)',
'cc-pVQZ(fi/sf/sc)',
'cc-pVQZ(pt/sf/fw)',
'cc-pVQZ(pt/sf/lc)',
'cc-pVQZ(pt/sf/sc)',
'cc-pVQZ(seg-opt)',
'cc-pVQZ-DK',
'cc-pVQZ-DK3',
'cc-pVQZ-F12',
'cc-pVQZ-F12(rev2)',
'cc-pVQZ-F12-OPTRI',
'cc-pVQZ-F12-OPTRI+',
'cc-pVQZ-JKFIT',
'cc-pVQZ-PP',
'cc-pVQZ-PP-RIFIT',
'cc-pVQZ-RIFIT',
'cc-pVQZ-X2C',
'cc-pVTZ',
'cc-pVTZ(fi/sf/fw)',
'cc-pVTZ(fi/sf/lc)',
'cc-pVTZ(fi/sf/sc)',
'cc-pVTZ(pt/sf/fw)',
'cc-pVTZ(pt/sf/lc)',
'cc-pVTZ(pt/sf/sc)',
'cc-pVTZ(seg-opt)',
'cc-pVTZ-DK',
'cc-pVTZ-DK3',
'cc-pVTZ-F12',
'cc-pVTZ-F12(rev2)',
'cc-pVTZ-F12-OPTRI',
'cc-pVTZ-F12-OPTRI+',
'cc-pVTZ-JKFIT',
'cc-pVTZ-PP',
'cc-pVTZ-PP-RIFIT',
'cc-pVTZ-RIFIT',
'cc-pVTZ-X2C',
'cc-pwCV5Z',
'cc-pwCV5Z-DK',
'cc-pwCV5Z-PP',
'cc-pwCV5Z-PP-RIFIT',
'cc-pwCV5Z-RIFIT',
'cc-pwCVDZ',
'cc-pwCVDZ-DK3',
'cc-pwCVDZ-PP',
'cc-pwCVDZ-PP-RIFIT',
'cc-pwCVDZ-RIFIT',
'cc-pwCVDZ-X2C',
'cc-pwCVQZ',
'cc-pwCVQZ-DK',
'cc-pwCVQZ-DK3',
'cc-pwCVQZ-PP',
'cc-pwCVQZ-PP-RIFIT',
'cc-pwCVQZ-RIFIT',
'cc-pwCVQZ-X2C',
'cc-pwCVTZ',
'cc-pwCVTZ-DK',
'cc-pwCVTZ-DK3',
'cc-pwCVTZ-PP',
'cc-pwCVTZ-PP-RIFIT',
'cc-pwCVTZ-RIFIT',
'cc-pwCVTZ-X2C',
'ccJ-pV5Z',
'ccJ-pVDZ',
'ccJ-pVQZ',
'ccJ-pVTZ',
'ccX-5Z',
'ccX-DZ',
'ccX-QZ',
'ccX-TZ',
'ccemd-2',
'ccemd-3',
'coemd-2',
'coemd-3',
'coemd-4',
'coemd-ref',
'd-aug-cc-pV5Z',
'd-aug-cc-pV6Z',
'd-aug-cc-pVDZ',
'd-aug-cc-pVQZ',
'd-aug-cc-pVTZ',
'deMon2k-DZVP-GGA',
'def2-ECP',
'def2-QZVP',
'def2-QZVP-RIFIT',
'def2-QZVPD',
'def2-QZVPP',
'def2-QZVPP-RIFIT',
'def2-QZVPPD',
'def2-QZVPPD-RIFIT',
'def2-SV(P)',
'def2-SV(P)-JKFIT',
'def2-SV(P)-RIFIT',
'def2-SVP',
'def2-SVP-RIFIT',
'def2-SVPD',
'def2-SVPD-RIFIT',
'def2-TZVP',
'def2-TZVP-RIFIT',
'def2-TZVPD',
'def2-TZVPD-RIFIT',
'def2-TZVPP',
'def2-TZVPP-RIFIT',
'def2-TZVPPD',
'def2-TZVPPD-RIFIT',
'def2-mTZVP',
'def2-mTZVPP',
'def2-mTZVPP-RIJ',
'def2-universal-JFIT',
'def2-universal-JKFIT',
'dgauss-a1-dftjfit',
'dgauss-a1-dftxfit',
'dgauss-a2-dftjfit',
'dgauss-a2-dftxfit',
'dgauss-dzvp',
'dgauss-dzvp2',
'dgauss-tzvp',
'dhf-ECP',
'dhf-QZVP',
'dhf-QZVPP',
'dhf-SV(P)',
'dhf-SVP',
'dhf-TZVP',
'dhf-TZVPP',
'dyall-aae2z',
'dyall-aae3z',
'dyall-aae4z',
'dyall-aae5z',
'dyall-acv2z',
'dyall-acv3z',
'dyall-acv4z',
'dyall-acv5z',
'dyall-ae2z',
'dyall-ae3z',
'dyall-ae4z',
'dyall-ae5z',
'dyall-av2z',
'dyall-av3z',
'dyall-av4z',
'dyall-av5z',
'dyall-cv2z',
'dyall-cv3z',
'dyall-cv4z',
'dyall-cv5z',
'dyall-v2z',
'dyall-v3z',
'dyall-v4z',
'dyall-v5z',
'epc-10s10p10d10f',
'epc-8s8p8d',
'jgauss-dzp',
'jgauss-qz2p',
'jgauss-qzp',
'jgauss-tzp1',
'jgauss-tzp2',
'jorge-5ZP',
'jorge-5ZP-DKH',
'jorge-5ZP-ZORA',
'jorge-6ZP',
'jorge-6ZP-DKH',
'jorge-6ZP-ZORA',
'jorge-A5ZP',
'jorge-A6ZP',
'jorge-ADZP',
'jorge-AQZP',
'jorge-ATZP',
'jorge-ATZP-ZORA',
'jorge-DZP',
'jorge-DZP-DKH',
'jorge-DZP-ZORA',
'jorge-QZP',
'jorge-QZP-DKH',
'jorge-QZP-ZORA',
'jorge-TZP',
'jorge-TZP-DKH',
'jorge-TZP-ZORA',
'jul-cc-pV(D+d)Z',
'jul-cc-pV(Q+d)Z',
'jul-cc-pV(T+d)Z',
'jun-cc-pV(D+d)Z',
'jun-cc-pV(Q+d)Z',
'jun-cc-pV(T+d)Z',
'lcecp-0-QZVP',
'lcecp-0-QZVPP',
'lcecp-0-SV(P)',
'lcecp-0-SVP',
'lcecp-0-TZVP',
'lcecp-0-TZVPP',
'lcecp-1-QZVP',
'lcecp-1-QZVPP',
'lcecp-1-SV(P)',
'lcecp-1-SVP',
'lcecp-1-TZVP',
'lcecp-1-TZVPP',
'lcecp-2-QZVP',
'lcecp-2-QZVPP',
'lcecp-2-SV(P)',
'lcecp-2-SVP',
'lcecp-2-TZVP',
'lcecp-2-TZVPP',
'm6-31G',
'm6-31G*',
'maug-cc-pV(D+d)Z',
'maug-cc-pV(Q+d)Z',
'maug-cc-pV(T+d)Z',
'may-cc-pV(Q+d)Z',
'may-cc-pV(T+d)Z',
'modified-LANL2DZ',
'pSBKJC',
'pV6Z',
'pV7Z',
'pc-0',
'pc-1',
'pc-2',
'pc-3',
'pc-4',
'pcH-1',
'pcH-2',
'pcH-3',
'pcH-4',
'pcJ-0',
'pcJ-0_2006',
'pcJ-1',
'pcJ-1_2006',
'pcJ-2',
'pcJ-2_2006',
'pcJ-3',
'pcJ-3_2006',
'pcJ-4',
'pcJ-4_2006',
'pcS-0',
'pcS-1',
'pcS-2',
'pcS-3',
'pcS-4',
'pcSseg-0',
'pcSseg-1',
'pcSseg-2',
'pcSseg-3',
'pcSseg-4',
'pcX-1',
'pcX-2',
'pcX-3',
'pcX-4',
'pcemd-2',
'pcemd-3',
'pcemd-4',
'pcseg-0',
'pcseg-1',
'pcseg-2',
'pcseg-3',
'pcseg-4',
'pecJ-1',
'pecJ-2',
'pob-DZVP-rev2',
'pob-TZVP',
'pob-TZVP-rev2',
's3-21G',
's3-21G*',
's6-31G',
's6-31G*',
'sap_grasp_large',
'sap_grasp_small',
'sap_helfem_large',
'sap_helfem_small',
'saug-ano-pV5Z',
'saug-ano-pVDZ',
'saug-ano-pVQZ',
'saug-ano-pVTZ',
'seg-cc-pV5Z-PP',
'seg-cc-pVDZ-PP',
'seg-cc-pVQZ-PP',
'seg-cc-pVTZ-PP',
'seg-cc-pwCV5Z-PP',
'seg-cc-pwCVDZ-PP',
'seg-cc-pwCVQZ-PP',
'seg-cc-pwCVTZ-PP',
'sigmaDZ',
'sigmaDZHF',
'sigmaQZ',
'sigmaSZHF',
'sigmaTZ',
'sigmaTZHF',
'un-ccemd-ref',
'un-pcemd-ref',
'x2c-JFIT',
'x2c-JFIT-universal',
'x2c-QZVPPall',
'x2c-QZVPPall-2c',
'x2c-QZVPPall-2c-s',
'x2c-QZVPPall-s',
'x2c-QZVPall',
'x2c-QZVPall-2c',
'x2c-QZVPall-2c-s',
'x2c-QZVPall-s',
'x2c-SV(P)all',
'x2c-SV(P)all-2c',
'x2c-SV(P)all-s',
'x2c-SVPall',
'x2c-SVPall-2c',
'x2c-SVPall-s',
'x2c-TZVPPall',
'x2c-TZVPPall-2c',
'x2c-TZVPPall-s',
'x2c-TZVPall',
'x2c-TZVPall-2c',
'x2c-TZVPall-s']
# Load basis set
basis = get_ao_basis_from_bse(h2o.atnums, atcoords=h2o.atcoords, basis_name="STO-6G", coord_types = "spherical")
Secular Equations#
Given a set of basis functions \(\{\chi_\lambda\}\), we can expand the molecular orbitals as linear combinations of the basis functions,
\[\phi_p(\mathbf{r}) = \sum_\lambda c_{\lambda p} \chi_\lambda(\mathbf{r})\]
where \(c_{\lambda p}\) is the coefficient of the \(\lambda\)-th basis function in the expansion of the \(p\)-th molecular orbital. Substituting this expansion into the one-electron equations and projecting onto a basis function \(\chi_\kappa\), we obtain the secular equations,
\[\sum_\lambda \left[ \int \chi_\kappa^*(\mathbf{r}) \left(-\frac{1}{2} \nabla^2 + v(\mathbf{r}) + \hat{w}[\{\phi_q\}](\mathbf{r})\right) \chi_\lambda(\mathbf{r}) d\mathbf{r} \right] c_{\lambda p} = \epsilon_p \sum_\lambda \left[ \int \chi_\kappa^*(\mathbf{r}) \chi_\lambda(\mathbf{r}) d\mathbf{r} \right] c_{\lambda p}\]
Defining the Hamiltonian matrix elements as
\[F_{\kappa \lambda} = \int \chi_\kappa^*(\mathbf{r}) \left(-\frac{1}{2} \nabla^2 + v(\mathbf{r}) + \hat{w}[\{\phi_q\}](\mathbf{r})\right) \chi_\lambda(\mathbf{r}) d\mathbf{r}\]
and the overlap matrix elements as
\[S_{\kappa \lambda} = \int \chi_\kappa^*(\mathbf{r}) \chi_\lambda(\mathbf{r}) d\mathbf{r}\]
the secular equations can be written in matrix form as
\[\mathbf{F} \mathbf{c}_p = \epsilon_p \mathbf{S} \mathbf{c}_p\]
where \(\mathbf{F}\) is the Fock matrix, \(\mathbf{S}\) is the overlap matrix, and \(\mathbf{c}_p\) is the vector of coefficients for the \(p\)-th molecular orbital. This is a generalized eigenvalue problem that can be solved to obtain the molecular orbitals and their corresponding energies. Note that one needs to (usually) solve for the \(\alpha\)- and \(\beta\)-spin orbitals separately, which leads to separate Fock matrices for each spin.
The Fock matrix is a sum of three parts, the kinetic energy matrix, the nuclear attraction matrix, and the average potential matrix,
\[F_{\kappa \lambda} = T_{\kappa \lambda} + V_{\kappa \lambda} + W_{\kappa \lambda}\]
where
\[T_{\kappa \lambda} = \int \chi_\kappa^*(\mathbf{r}) \left(-\frac{1}{2} \nabla^2\right) \chi_\lambda(\mathbf{r}) d\mathbf{r}\]
\[V_{\kappa \lambda} = \int \chi_\kappa^*(\mathbf{r}) v(\mathbf{r}) \chi_\lambda(\mathbf{r}) d\mathbf{r}\]
\[W_{\kappa \lambda} = \int \chi_\kappa^*(\mathbf{r}) \hat{w}[\{\phi_q\}](\mathbf{r}) \chi_\lambda(\mathbf{r}) d\mathbf{r}\]
Density Matrices#
These formulae are most simply expressed in the language of density matrices. The one-electron reduced density matrix is defined as
\[
\gamma(\mathbf{r}, \mathbf{r}') = \sum_{i=1}^K n_p \phi_p(\mathbf{r}) \phi_p^*(\mathbf{r}')
\]
where \(n_p\) is the occupation number of the \(p\)-th molecular orbital. The electron density can be obtained from the density matrix by evaluating it at \(\mathbf{r} = \mathbf{r}'\),
\[\rho(\mathbf{r}) = \gamma(\mathbf{r}, \mathbf{r}) = \sum_{i=1}^K n_p |\phi_p(\mathbf{r})|^2\]
In terms of atomic orbitals, the density matrix can be expressed as
\[\gamma(\mathbf{r}, \mathbf{r}') = \sum_{p=1}^K n_p c_{\kappa p}^* c_{\lambda p}^* \chi_{\kappa}^*(\mathbf{r}') \chi_{\lambda}(\mathbf{r})\]
The quantity
\[
D_{\kappa \lambda} = \sum_{p=1}^K n_p c_{\kappa p}^* c_{\lambda p}
\]
is known as the density matrix (in the atomic orbital basis). The density matrix is normalized by its relation to the total electron number,
\[
\sum_{\kappa \lambda} D_{\kappa \lambda}S_{\kappa \lambda} = \iint \gamma(\mathbf{r}, \mathbf{r}') \delta(\mathbf{r} - \mathbf{r}') d\mathbf{r} d\mathbf{r}' = \sum_{\kappa,\lambda}\sum_{p=1}^K n_p c_{\kappa p}^* c_{\lambda p} S_{\kappa \lambda} = N
\]
import numpy as np
import scipy
from scipy.special import erf
import json
import atomdb
from iodata import load_one
from grid.onedgrid import GaussLegendre
from grid.rtransform import BeckeRTransform
from grid.becke import BeckeWeights
from grid.molgrid import MolGrid
from gbasis.integrals.overlap import overlap_integral
from gbasis.integrals.kinetic_energy import kinetic_energy_integral
from gbasis.integrals.nuclear_electron_attraction import nuclear_electron_attraction_integral
from gbasis.integrals.electron_repulsion import electron_repulsion_integral
from gbasis.evals.eval import evaluate_basis
from gbasis.evals.density import evaluate_density
# Make a decent(ish) grid.
def make_grid(atnums, atcoords, rad_points=100, preset="ultrafine"):
'''Make a grid object from a molecule
Args:
atnums (list of int): Atomic numbers of the atoms in the molecule
atcoords (numpy.ndarray): Coordinates of the atoms in the molecule
rad_points (int): number of radial points
preset (str): preset grid type
"coarse", "medium", "fine", "veryfine", "ultrafine", "insane"
Returns:
mol_grid (Grid): grid object
'''
# create a radial grid with rad_points points
oned_grid = GaussLegendre(rad_points)
radial_grid = BeckeRTransform(0.0, R=1.5).transform_1d_grid(oned_grid)
# create a pruned molecular grid with the preset angular convention
mol_grid = MolGrid.from_preset(atnums, atcoords, preset, rgrid=radial_grid, store=True)
return mol_grid
# Compute overlap integrals
def compute_overlap(ao_basis):
r"""Compute the overlap matrix in the AO basis.
Parameters
----------
ao_basis : list of `GeneralizedContractionShell` instances
Atomic orbital basis set for the molecule.
Returns
-------
S : numpy.ndarray[float, shape=(n_ao, n_ao)]
Overlap matrix in the AO basis.
"""
return overlap_integral(ao_basis)
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In[6], line 7
3 from scipy.special import erf
4 import json
5 import atomdb
6 from iodata import load_one
----> 7 from grid.onedgrid import GaussLegendre
8 from grid.rtransform import BeckeRTransform
9 from grid.becke import BeckeWeights
10 from grid.molgrid import MolGrid
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/__init__.py:25
21 """Grid Module."""
22 # flake8: noqa
---> 25 from grid.angular import *
26 from grid.atomgrid import *
27 from grid.basegrid import *
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/angular.py:75
72 import numpy as np
73 from importlib_resources import files
---> 75 from grid.basegrid import Grid
77 # Lebedev dictionary for grid's number of points (keys) and degrees (values)
78 LEBEDEV_NPOINTS = {
79 6: 3,
80 18: 5,
(...) 110 5810: 131,
111 }
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/basegrid.py:25
22 import numpy as np
23 from scipy.spatial import cKDTree
---> 25 from grid.utils import (
26 convert_cart_to_sph,
27 generate_orders_horton_order,
28 solid_harmonics,
29 )
32 class Grid:
33 """Basic Grid class for grid information storage."""
File /opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/grid/utils.py:23
20 """Utility Module."""
22 import numpy as np
---> 23 from scipy.special import sph_harm
25 _bragg = np.array(
26 [
27 np.nan, # index 0, place holder
(...) 114 ]
115 )
117 _cambridge = np.array(
118 [
119 np.nan, # index 0, place holder
(...) 206 ]
207 )
ImportError: cannot import name 'sph_harm' from 'scipy.special' (/opt/hostedtoolcache/Python/3.14.3/x64/lib/python3.14/site-packages/scipy/special/__init__.py)
# Evaluate the Coulomb potential of a Gaussian and r^2 times a Gaussian (on a grid)
def coulomb_gaussian_s(r, a):
"""
Coulomb potential of exp(-a r^2), for scalar or array r.
Returns V(r) = (pi/a)^(3/2) * erf(sqrt(a) r) / r
with the correct finite limit at r=0.
"""
r = np.asarray(r, dtype=float)
out = np.empty_like(r)
mask0 = np.isclose(r, 0.0)
mask = ~mask0
out[mask] = (np.pi / a) ** 1.5 * erf(np.sqrt(a) * r[mask]) / r[mask]
out[mask0] = 2.0 * np.pi / a
return out
def coulomb_gaussian_r2(r, b):
"""
Coulomb potential of r^2 exp(-b r^2), for scalar or array r.
Returns
V(r) = (3 pi^(3/2)/(2 b^(5/2))) * erf(sqrt(b) r)/r
+ (2 pi / b^2) * exp(-b r^2)
with the correct finite limit at r=0.
"""
r = np.asarray(r, dtype=float)
out = np.empty_like(r)
mask0 = np.isclose(r, 0.0)
mask = ~mask0
out[mask] = (
(3.0 * np.pi**1.5 / (2.0 * b**2.5)) * erf(np.sqrt(b) * r[mask]) / r[mask]
+ (2.0 * np.pi / b**2) * np.exp(-b * r[mask] ** 2)
)
out[mask0] = 5.0 * np.pi / b**2
return out
def coulomb_potential(r, c, a, d=None, b=None):
"""
Coulomb potential of
rho(r) = sum_i c_i exp(-a_i r^2) + sum_j d_j r^2 exp(-b_j r^2)
Parameters
----------
r : float or ndarray
Radial coordinate(s).
c, a : 1D arrays
Coefficients and exponents for the pure Gaussians.
d, b : 1D arrays, optional
Coefficients and exponents for the r^2 Gaussians.
Returns
-------
V : ndarray
Coulomb potential evaluated at r.
"""
r = np.asarray(r, dtype=float)
V = np.zeros_like(r)
c = np.asarray(c, dtype=float)
a = np.asarray(a, dtype=float)
for ci, ai in zip(c, a):
V += ci * coulomb_gaussian_s(r, ai)
if d is not None and b is not None:
d = np.asarray(d, dtype=float)
b = np.asarray(b, dtype=float)
for dj, bj in zip(d, b):
V += dj * coulomb_gaussian_r2(r, bj)
return V
def compute_harris_lsda(atnums, atcoords, ao_basis, n_electrons):
r"""Compute the LSDA Harris-functional guess for the exchange-correlation potential matrix.
The LSDA Harris-functional guess is computed by evaluating the Fock matrix using the promolecular density
as the input density. The exchange-correlation potential is computed using the local spin-density approximation (LSDA).
Parameters
----------
atnums : list of int
Atomic numbers of the atoms in the molecule.
atcoords : numpy.ndarray[float, shape=(n_atoms, 3)]
Coordinates of the atoms in the molecule.
ao_basis : list of `GeneralizedContractionShell` instances
Atomic orbital basis set for the molecule.
n_electrons : numpy.ndarray[int, shape=(2,)]
Number of electrons in the molecule (of each spin channel).
Returns
-------
harris_guess : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
LSDA exchange-correlation potential matrix elements. The first component of the last dimension
corresponds to the alpha spin, and the second component corresponds to the beta spin.
"""
# We use `AtomDB` to compute the promolecular density of the molecule. We assume
# that the spin-densities are proportional to the total density.
# Make a grid. No need to be super accurate here.
mol_grid = make_grid(atnums, atcoords, rad_points=50, preset="medium")
# Compute the promolecular density on the grid.
# evaluate the promolecular density on the grid
promol = atomdb.make_promolecule(atnums, atcoords, dataset="gaussian")
promol_dens = promol.density(mol_grid.points, log=True)
# Get spin-densities.
n_a, n_b = n_electrons
#The promolecular density is always neutral, so the number of electrons *expected* is the
#sum of the atomic numbers.
n_total = np.sum(atnums)
# Components of the density are n_alpha/n_total and n_beta/n_total times the total density.
promol_dens_a = promol_dens * n_a / (n_total)
promol_dens_b = promol_dens * n_b / (n_total)
# Compute the LSDA exchange-correlation potential on the grid.
vxc_a = (6/np.pi)**(1/3) * promol_dens_a ** (1/3) # LSDA exchange potential for alpha spin
vxc_b = (6/np.pi)**(1/3) * promol_dens_b ** (1/3) # LSDA exchange potential for beta spin
# Evaluate the exchange-correlation potential matrix elements in the AO basis.
phi = evaluate_basis(basis, mol_grid.points)
# AO matrices for alpha and beta
Fxc = np.empty((phi.shape[0], phi.shape[0], 2))
Fxc[:, :, 0] = np.einsum("ip,p,jp->ij", phi, mol_grid.weights * vxc_a, phi)
Fxc[:, :, 1] = np.einsum("ip,p,jp->ij", phi, mol_grid.weights * vxc_b, phi)
return Fxc
def compute_harris_hartree(atnums, atcoords, ao_basis, n_electrons):
r"""Compute the Hartree Harris-functional guess for the Fock matrix.
The Hartree Harris-functional guess is computed by evaluating the Fock matrix using the promolecular density
as the input density.
Parameters
----------
atnums : list of int
Atomic numbers of the atoms in the molecule.
atcoords : numpy.ndarray[float, shape=(n_atoms, 3)]
Coordinates of the atoms in the molecule.
ao_basis : list of `GeneralizedContractionShell` instances
Atomic orbital basis set for the molecule.
n_electrons : numpy.ndarray[int, shape=(2,)]
Number of electrons in the molecule (of each spin channel).
Returns
-------
hartree_pot : numpy.ndarray[float, shape=(n_ao, n_ao)]
Hartree potential matrix elements.
"""
n_a, n_b = n_electrons
n_total = sum(atnums)
# We use `BFit` to obtain an approximate density with the same number of electrons as the promolecular density,
# expanded in s-type Gaussians. The Coulomb potential is then easy to evaluate on a grid.
# Make a grid. No need to be super accurate here.
mol_grid = make_grid(atnums, atcoords, rad_points=50, preset="medium")
# Load atomic fit data from BFit.
url = "https://raw.githubusercontent.com/theochem/BFit/master/bfit/data/kl_slsqp_results.json"
urlretrieve(url, "kl_slsqp_results.json")
with open("kl_slsqp_results.json") as file:
data = json.load(file)
# For each atom, we compute the Coulomb potential on the grid and sum them up to get the total Hartree potential.
# Loop over the atoms
hartree_pot = np.zeros((mol_grid.points.shape[0],))
for atnum, atcoord in zip(atnums, atcoords):
# Get the BFit model for the atom
# Convert the atomic number to the element symbol
element = iodata.periodic.num2sym[atnum]
# BFit uses lowercase element symbols, so we convert to lowercase.
element = element.lower()
num_s = data[element]["num_s"]
num_p = data[element]["num_p"]
coefficients = np.array(data[element]["coeffs"])
exponents = np.array(data[element]["exps"])
# For every grid point, compute its distance to the atom in question.
dist = np.linalg.norm(mol_grid.points - atcoord, axis=1)
# Evaluate the Coulomb potential from the atom at each grid point using the BFit model.
# Break the coefficients and exponents into s-type and p-type contributions.
# Our code uses unnormalized Gaussians but BFit produces normalized Gaussians, so we need to
# multiply the coefficients by the normalization factors (cf. eqs. (22) and (23)).
c_s = coefficients[:num_s] * (exponents[:num_s] / np.pi )**(3/2)
a_s = exponents[:num_s]
c_p = coefficients[num_s:num_s+num_p] * 2/3 *exponents[num_s:num_s+num_p]**(5/2) / np.pi**(3/2)
a_p = exponents[num_s:num_s+num_p]
hartree_pot += coulomb_potential(dist, c_s, a_s, d=c_p, b=a_p)
# Evaluate the Hartree potential matrix elements in the AO basis.
phi = evaluate_basis(ao_basis, mol_grid.points)
vj = hartree_pot * (n_a +n_b) / n_total # Scale the Hartree potential if the molecule is not neutral.
# AO matrices for alpha and beta
fj = np.empty((phi.shape[0], phi.shape[0]))
fj[:, :] = np.einsum("ip,p,jp->ij", phi, mol_grid.weights * vj, phi)
return fj
def compute_harris_guess(atnums, atcoords, ao_basis, n_electrons, guess_type='fermi-amaldi'):
r"""Compute the Harris-functional guess for the Fock matrix.
The Harris-functional guess is computed by evaluating the Fock matrix using the promolecular density
as the input density.
Parameters
----------
atnums : list of int
Atomic numbers of the atoms in the molecule.
atcoords : numpy.ndarray[float, shape=(n_atoms, 3)]
Coordinates of the atoms in the molecule.
ao_basis : list of `GeneralizedContractionShell` instances
Atomic orbital basis set for the molecule.
n_electrons : numpy.ndarray[int, shape=(2,)]
Number of electrons in the molecule (of each spin channel).
guess_type : {"fermi-amaldi", "lsda", "hartree", "core"}
Type of Harris-functional guess to compute.
Choices are core, LSDA (exchange-only), Fermi-Amaldi, and Hartree.
Returns
-------
harris_guess : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
Harris-functional guess for the Fock matrix. The first component of the last dimension corresponds to the
alpha spin, and the second component corresponds to the beta spin.
"""
n_a, n_b = n_electrons
kin_energy = kinetic_energy_integral(ao_basis)
nuc_attraction = nuclear_electron_attraction_integral(ao_basis, atcoords, atnums)
# Assign the core Hamiltonian for a and b.
fock = np.empty((kin_energy.shape[0], kin_energy.shape[1], 2))
fock[:, :, 0] = kin_energy + nuc_attraction
fock[:, :, 1] = fock[:, :, 0]
# Assign the hartree potential using compute_harris_hartree.
if guess_type != "core":
vj = compute_harris_hartree(atnums, atcoords, ao_basis, n_electrons)
fock[:, :, 0] += vj
fock[:, :, 1] += vj
if guess_type == "lsda":
fock += compute_harris_lsda(atnums, atcoords, ao_basis, n_electrons)
elif guess_type == "fermi-amaldi":
fock[:, :, 0] -= (1/n_a) * vj
fock[:, :, 1] -= (1/n_b) * vj
#make sure there is no weird case where the user choose a forbidden option
elif (guess_type != "core" and guess_type != "hartree"):
raise ValueError(f"Invalid guess type: {guess_type}. Must be one of 'core', 'fermi-amaldi', 'lsda', or 'hartree'.")
return fock
# Solve for the eigenvalues and eigenvectors of the Harris guess Fock matrix to get the initial guess density matrix.
import scipy
from sympy.abc import phi
s = compute_overlap(basis)
n_electrons = np.array([5, 5]) # Number of electrons in the molecule (of each spin channel).
fock = compute_harris_guess(h2o.atnums, h2o.atcoords, basis, n_electrons, guess_type='core')
eps_a, phi_a = scipy.linalg.eigh(fock[:, :, 0], s) # Solve the generalized eigenvalue problem for the alpha spin channel
eps_b, phi_b = scipy.linalg.eigh(fock[:, :, 1], s) # Solve the generalized eigenvalue problem for the beta spin channel
# Put the eigenvalues and MOs into a 2-d array for easier indexing.
eps = np.empty((eps_a.shape[0], 2))
eps[:, 0] = eps_a
eps[:, 1] = eps_b
phi = np.empty((phi_a.shape[0], phi_a.shape[1], 2))
phi[:, :, 0] = phi_a
phi[:, :, 1] = phi_b
print(phi.shape)
print(eps)
(7, 7, 2)
[[-33.03926628 -33.03926628]
[ -8.34758288 -8.34758288]
[ -7.77418972 -7.77418972]
[ -7.48618707 -7.48618707]
[ -7.47336021 -7.47336021]
[ -4.22880819 -4.22880819]
[ -4.2186479 -4.2186479 ]]
def compute_density_matrix(eps, phi, n_electrons):
r"""Compute the density matrix from the eigenvalues and eigenvectors of the Fock matrix.
Parameters
----------
eps : numpy.ndarray[float, shape=(n_ao, 2)]
Eigenvalues of the Fock matrix.
phi : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
Eigenvectors of the Fock matrix.
n_electrons : numpy.ndarray[int, shape=(2,)]
Number of electrons in the molecule (of a given spin channel).
Returns
-------
density_matrix : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
Density matrix in the AO basis.
"""
n_a, n_b = n_electrons
# Compute the density matrix as D = sum_i |phi_i><phi_i| where i runs over the occupied orbitals.
density_matrix = np.zeros((phi.shape[0], phi.shape[0], 2))
density_matrix[:, :, 0] = np.einsum("pi,qi->pq", phi[:, :n_a, 0], phi[:, :n_a, 0])
density_matrix[:, :, 1] = np.einsum("pi,qi->pq", phi[:, :n_b, 1], phi[:, :n_b, 1])
return density_matrix
def check_normalization(density_matrix, n_electrons, overlap):
r"""Check the normalization of the density matrix.
Parameters
----------
density_matrix : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
Density matrix in the AO basis.
n_electrons : numpy.ndarray[int, shape=(2,)]
Number of electrons in the molecule (of each spin channel).
overlap : numpy.ndarray[float, shape=(n_ao, n_ao)]
Overlap matrix in the AO basis.
Returns
-------
is_normalized : bool
True if the density matrix is normalized, False otherwise.
"""
n_a, n_b = n_electrons
# The number of electrons can be computed as N = Tr(D S) where D is the density matrix and S is the overlap matrix.
n_a_computed = np.trace(density_matrix[:, :, 0] @ overlap)
n_b_computed = np.trace(density_matrix[:, :, 1] @ overlap)
print(n_a_computed, n_b_computed, "number of electrons in each spin channel computed from the density matrix")
return np.isclose(n_a_computed, n_a) and np.isclose(n_b_computed, n_b)
print(check_normalization(compute_density_matrix(eps, phi, n_electrons), n_electrons, s ))
4.999999999999999 4.999999999999999 number of electrons in each spin channel computed from the density matrix
True
def get_eri_ao(basis):
"""
Return AO two-electron integrals in chemist notation:
eri[mu, nu, lam, sig] = (mu nu | lam sig)
"""
return electron_repulsion_integral(basis, notation="chemist")
vee = get_eri_ao(basis)
def build_j_ao(D, eri):
"""
Build the AO Coulomb matrix
J_mn = sum_ls D_ls (mn|ls)
Parameters
----------
D : (nao, nao) ndarray
AO density matrix.
eri : (nao, nao, nao, nao) ndarray
AO ERI tensor in chemist notation.
Returns
-------
J : (nao, nao) ndarray
AO Coulomb matrix.
"""
D = np.asarray(D)
return np.einsum("ls,mnls->mn", D, eri, optimize=True)
def build_k_ao(D, eri):
"""
Build the AO exchange matrix
K_mn = sum_ls D_ls (ml|ns)
Parameters
----------
D : (nao, nao) ndarray
AO density matrix.
eri : (nao, nao, nao, nao) ndarray
AO ERI tensor in chemist notation.
Returns
-------
K : (nao, nao) ndarray
AO exchange matrix.
"""
D = np.asarray(D)
return np.einsum("ls,mlns->mn", D, eri, optimize=True)
def build_vxc_lsda(atnums, atcoords, ao_basis, D):
r"""Compute the LSDA contribution to the Fock matrix for the exchange-correlation potential matrix.
The LSDA is computed by evaluating the Fock matrix using the electron density
as computed from the density matrix.
The exchange-correlation potential is computed using the
local spin-density approximation (LSDA).
Parameters
----------
atnums : list of int
Atomic numbers of the atoms in the molecule.
atcoords : numpy.ndarray[float, shape=(n_atoms, 3)]
Coordinates of the atoms in the molecule.
ao_basis : list of `GeneralizedContractionShell` instances
Atomic orbital basis set for the molecule.
D : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
Density matrix in the AO basis
Returns
-------
F_lsda : numpy.ndarray[float, shape=(n_ao, n_ao, 2)]
LSDA exchange-correlation potential matrix elements.
The first component of the last dimension corresponds to the alpha
spin, and the second component corresponds to the beta spin.
"""
# We use `AtomDB` to compute the promolecular density of the molecule. We assume
# that the spin-densities are proportional to the total density.
# Make a grid. No need to be super accurate here.
mol_grid = make_grid(atnums, atcoords, rad_points=100, preset="veryfine")
# Compute the density on the grid.
rdm_a = D[:, :, 0]
rdm_b = D[:, :, 1]
rho_a = evaluate_density(rdm_a, ao_basis, mol_grid.points)
rho_b = evaluate_density(rdm_b, ao_basis, mol_grid.points)
#Check normalization
print("number of alpha-spin electrons", mol_grid.integrate(rho_a))
print("number of beta-spin electrons", mol_grid.integrate(rho_b))
# Compute the LSDA exchange-correlation potential on the grid.
vxc_a = (6/np.pi)**(1/3) * rho_a ** (1/3) # LSDA exchange potential for alpha spin
vxc_b = (6/np.pi)**(1/3) * rho_b ** (1/3) # LSDA exchange potential for beta spin
# Evaluate the exchange-correlation potential matrix elements in the AO basis.
phi = evaluate_basis(basis, mol_grid.points)
# AO matrices for alpha and beta
Fxc = np.empty((phi.shape[0], phi.shape[0], 2))
Fxc[:, :, 0] = np.einsum("ip,p,jp->ij", phi, mol_grid.weights * vxc_a, phi)
Fxc[:, :, 1] = np.einsum("ip,p,jp->ij", phi, mol_grid.weights * vxc_b, phi)
return Fxc
def build_fock_ao(h_core, density_matrix, eri, overlap, n_electrons, atnums, atcoords, ao_basis, calc_type = 'fermi-amaldi'):
"""
Build the AO Fock matrix for a given density matrix and ERI tensor.
Parameters
----------
h_core : (nao, nao) ndarray
AO core Hamiltonian matrix.
density_matrix : (nao, nao, 2) ndarray
AO density matrix for both spin channels.
eri : (nao, nao, nao, nao) ndarray
AO ERI tensor in chemist notation.
overlap : (nao, nao) ndarray
AO overlap matrix.
n_electrons : array-like of shape (2,)
Number of electrons in the molecule (of each spin channel).
calc_type : str
Type of calculation to perform.
Choices are 'fermi-amaldi', 'lsda', 'hartree', and 'hartree-fock'.
Returns
-------
F_a : (nao, nao) ndarray
AO Fock matrix for the alpha spin channel.
F_b : (nao, nao) ndarray
AO Fock matrix for the beta spin channel.
"""
D_a = density_matrix[:, :, 0]
D_b = density_matrix[:, :, 1]
J_a = build_j_ao(D_a, eri)
J_b = build_j_ao(D_b, eri)
K_a = build_k_ao(D_a, eri)
K_b = build_k_ao(D_b, eri)
# For consistency check that the number of electrons targetted is consistent with the density matrix given.
check_normalization(density_matrix, n_electrons, overlap)
if calc_type == 'fermi-amaldi':
F_a = h_core + J_a + J_b - (1/n_electrons[0]) * J_a
F_b = h_core + J_a + J_b - (1/n_electrons[1]) * J_b
elif calc_type == 'hartree':
F_a = h_core + J_a + J_b
F_b = h_core + J_a + J_b
elif calc_type == 'lsda':
F_xc = build_vxc_lsda(atnums, atcoords, ao_basis, density_matrix)
F_a = h_core + J_a + J_b + F_xc[:, :, 0]
F_b = h_core + J_a + J_b + F_xc[:, :, 1]
elif calc_type == 'hartree-fock':
F_a = h_core + J_a + J_b - K_a
F_b = h_core + J_a + J_b - K_b
else:
raise ValueError(f"Invalid calculation type: {calc_type}. Must be one of 'fermi-amaldi', 'lsda', 'hartree', or 'hartree-fock'.")
return F_a, F_b
def scf_cycle(h_core, density_matrix, eri, s, n_electrons, atnums, atcoords, ao_basis, calc_type='hartree', max_iter=100, tol=1e-4):
"""
Perform the SCF cycle to self-consistency.
Parameters
----------
h_core : (nao, nao) ndarray
AO core Hamiltonian matrix.
density_matrix : (nao, nao, 2) ndarray
Initial AO density matrix for both spin channels.
eri : (nao, nao, nao, nao) ndarray
AO ERI tensor in chemist notation.
s : (nao, nao) ndarray
AO overlap matrix.
n_electrons : array-like of shape (2,)
Number of electrons in the molecule (of each spin channel).
calc_type : str
Type of calculation to perform. Choices are 'fermi-amaldi', 'lsda', and 'hartree'.
max_iter : int
Maximum number of SCF iterations to perform.
tol : float
Convergence tolerance for the density matrix.
Returns
-------
eps : (n_ao, 2) ndarray
Eigenvalues of the final Fock matrix.
phi : (n_ao, n_ao, 2) ndarray
Eigenvectors of the final Fock matrix.
"""
nao = h_core.shape[0]
for iteration in range(max_iter):
F_a, F_b = build_fock_ao(h_core, density_matrix, eri, s, n_electrons, atnums, atcoords, ao_basis, calc_type)
eps_a, phi_a = scipy.linalg.eigh(F_a, s)
eps_b, phi_b = scipy.linalg.eigh(F_b, s)
# Put the eigenvalues and MOs into a 2-d array for easier indexing.
eps = np.empty((eps_a.shape[0], 2))
eps[:, 0] = eps_a
eps[:, 1] = eps_b
phi = np.empty((phi_a.shape[0], phi_a.shape[1], 2))
phi[:, :, 0] = phi_a
phi[:, :, 1] = phi_b
# Compute the new density matrix.
density_matrix_new = compute_density_matrix(eps, phi, n_electrons)
# Check for convergence
delta_dm = np.linalg.norm(density_matrix - density_matrix_new)
print(f"Iteration {iteration+1}: delta_dm = {delta_dm:.6e}")
if delta_dm < tol:
print(f"SCF converged in {iteration+1} iterations.")
return eps, phi
density_matrix = density_matrix_new *0.1 + density_matrix *0.9 # Simple damping to help convergence.
raise RuntimeError("SCF did not converge within the maximum number of iterations.")
h_core = kinetic_energy_integral(basis) + nuclear_electron_attraction_integral(basis, h2o.atcoords, h2o.atnums)
eri = get_eri_ao(basis)
eps, phi = scf_cycle(h_core, compute_density_matrix(eps, phi, n_electrons), eri, s, n_electrons, h2o.atnums, h2o.atcoords, basis, calc_type='lsda', max_iter=200)
print("Final eigenvalues (orbital energies):", eps)
4.999999999999999 4.999999999999999 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 4.999999937316277
number of beta-spin electrons 4.999999937316277
Iteration 1: delta_dm = 3.574831e+00
4.999999999999999 4.999999999999999 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 4.999999959181298
number of beta-spin electrons 4.999999959181298
Iteration 2: delta_dm = 3.158304e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 4.999999978598807
number of beta-spin electrons 4.999999978598807
Iteration 3: delta_dm = 2.797596e+00
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 4.999999995599962
number of beta-spin electrons 4.999999995599962
Iteration 4: delta_dm = 2.479467e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000010060662
number of beta-spin electrons 5.000000010060662
Iteration 5: delta_dm = 2.202876e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000021835491
number of beta-spin electrons 5.000000021835491
Iteration 6: delta_dm = 1.967207e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000030849837
number of beta-spin electrons 5.000000030849837
Iteration 7: delta_dm = 1.769761e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000371803734
number of beta-spin electrons 5.0000000371803734
Iteration 8: delta_dm = 1.604899e+00
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000041084401
number of beta-spin electrons 5.000000041084401
Iteration 9: delta_dm = 1.465153e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000042961406
number of beta-spin electrons 5.000000042961406
Iteration 10: delta_dm = 1.343284e+00
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000043273542
number of beta-spin electrons 5.000000043273542
Iteration 11: delta_dm = 1.233747e+00
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000042468313
number of beta-spin electrons 5.000000042468313
Iteration 12: delta_dm = 1.132994e+00
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000040930003
number of beta-spin electrons 5.000000040930003
Iteration 13: delta_dm = 1.039061e+00
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000038962035
number of beta-spin electrons 5.000000038962035
Iteration 14: delta_dm = 9.509891e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000036789702
number of beta-spin electrons 5.000000036789702
Iteration 15: delta_dm = 8.683707e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000034572069
number of beta-spin electrons 5.000000034572069
Iteration 16: delta_dm = 7.910514e-01
4.999999999999999 4.999999999999999 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000032416173
number of beta-spin electrons 5.000000032416173
Iteration 17: delta_dm = 7.189670e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000030390001
number of beta-spin electrons 5.000000030390001
Iteration 18: delta_dm = 6.520575e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000028533069
number of beta-spin electrons 5.000000028533069
Iteration 19: delta_dm = 5.902294e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000026864776
number of beta-spin electrons 5.000000026864776
Iteration 20: delta_dm = 5.333422e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000253904995
number of beta-spin electrons 5.0000000253904995
Iteration 21: delta_dm = 4.812099e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000002410629
number of beta-spin electrons 5.00000002410629
Iteration 22: delta_dm = 4.336074e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000023002175
number of beta-spin electrons 5.000000023002175
Iteration 23: delta_dm = 3.902807e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000022064635
number of beta-spin electrons 5.000000022064635
Iteration 24: delta_dm = 3.509569e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000212782565
number of beta-spin electrons 5.0000000212782565
Iteration 25: delta_dm = 3.153531e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000020626946
number of beta-spin electrons 5.000000020626946
Iteration 26: delta_dm = 2.831849e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000200948
number of beta-spin electrons 5.0000000200948
Iteration 27: delta_dm = 2.541721e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019666637
number of beta-spin electrons 5.000000019666637
Iteration 28: delta_dm = 2.280439e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019328234
number of beta-spin electrons 5.000000019328234
Iteration 29: delta_dm = 2.045425e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001906659
number of beta-spin electrons 5.00000001906659
Iteration 30: delta_dm = 1.834250e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018870018
number of beta-spin electrons 5.000000018870018
Iteration 31: delta_dm = 1.644648e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018728094
number of beta-spin electrons 5.000000018728094
Iteration 32: delta_dm = 1.474525e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000186316385
number of beta-spin electrons 5.0000000186316385
Iteration 33: delta_dm = 1.321955e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018572646
number of beta-spin electrons 5.000000018572646
Iteration 34: delta_dm = 1.185177e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018544168
number of beta-spin electrons 5.000000018544168
Iteration 35: delta_dm = 1.062589e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018540272
number of beta-spin electrons 5.000000018540272
Iteration 36: delta_dm = 9.527401e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001855588
number of beta-spin electrons 5.00000001855588
Iteration 37: delta_dm = 8.543163e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018586697
number of beta-spin electrons 5.000000018586697
Iteration 38: delta_dm = 7.661338e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018629102
number of beta-spin electrons 5.000000018629102
Iteration 39: delta_dm = 6.871272e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018680075
number of beta-spin electrons 5.000000018680075
Iteration 40: delta_dm = 6.163391e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018737148
number of beta-spin electrons 5.000000018737148
Iteration 41: delta_dm = 5.529101e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018798243
number of beta-spin electrons 5.000000018798243
Iteration 42: delta_dm = 4.960700e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018861698
number of beta-spin electrons 5.000000018861698
Iteration 43: delta_dm = 4.451288e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018926152
number of beta-spin electrons 5.000000018926152
Iteration 44: delta_dm = 3.994685e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000018990566
number of beta-spin electrons 5.000000018990566
Iteration 45: delta_dm = 3.585362e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019054071
number of beta-spin electrons 5.000000019054071
Iteration 46: delta_dm = 3.218372e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019116024
number of beta-spin electrons 5.000000019116024
Iteration 47: delta_dm = 2.889289e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019175939
number of beta-spin electrons 5.000000019175939
Iteration 48: delta_dm = 2.594155e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019233458
number of beta-spin electrons 5.000000019233458
Iteration 49: delta_dm = 2.329428e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019288338
number of beta-spin electrons 5.000000019288338
Iteration 50: delta_dm = 2.091941e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019340431
number of beta-spin electrons 5.000000019340431
Iteration 51: delta_dm = 1.878860e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000193896215
number of beta-spin electrons 5.0000000193896215
Iteration 52: delta_dm = 1.687651e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019435909
number of beta-spin electrons 5.000000019435909
Iteration 53: delta_dm = 1.516045e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019479309
number of beta-spin electrons 5.000000019479309
Iteration 54: delta_dm = 1.362013e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019519861
number of beta-spin electrons 5.000000019519861
Iteration 55: delta_dm = 1.223736e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019557655
number of beta-spin electrons 5.000000019557655
Iteration 56: delta_dm = 1.099589e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019592784
number of beta-spin electrons 5.000000019592784
Iteration 57: delta_dm = 9.881140e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001962535
number of beta-spin electrons 5.00000001962535
Iteration 58: delta_dm = 8.880067e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001965549
number of beta-spin electrons 5.00000001965549
Iteration 59: delta_dm = 7.980980e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019683349
number of beta-spin electrons 5.000000019683349
Iteration 60: delta_dm = 7.173408e-03
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000197090175
number of beta-spin electrons 5.0000000197090175
Iteration 61: delta_dm = 6.447964e-03
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019732669
number of beta-spin electrons 5.000000019732669
Iteration 62: delta_dm = 5.796236e-03
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019754407
number of beta-spin electrons 5.000000019754407
Iteration 63: delta_dm = 5.210681e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019774344
number of beta-spin electrons 5.000000019774344
Iteration 64: delta_dm = 4.684536e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019792648
number of beta-spin electrons 5.000000019792648
Iteration 65: delta_dm = 4.211735e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000198094146
number of beta-spin electrons 5.0000000198094146
Iteration 66: delta_dm = 3.786839e-03
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000198247525
number of beta-spin electrons 5.0000000198247525
Iteration 67: delta_dm = 3.404966e-03
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019838779
number of beta-spin electrons 5.000000019838779
Iteration 68: delta_dm = 3.061736e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019851568
number of beta-spin electrons 5.000000019851568
Iteration 69: delta_dm = 2.753220e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019863272
number of beta-spin electrons 5.000000019863272
Iteration 70: delta_dm = 2.475888e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019873907
number of beta-spin electrons 5.000000019873907
Iteration 71: delta_dm = 2.226575e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001988363
number of beta-spin electrons 5.00000001988363
Iteration 72: delta_dm = 2.002438e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019892458
number of beta-spin electrons 5.000000019892458
Iteration 73: delta_dm = 1.800924e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019900491
number of beta-spin electrons 5.000000019900491
Iteration 74: delta_dm = 1.619741e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019907808
number of beta-spin electrons 5.000000019907808
Iteration 75: delta_dm = 1.456829e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019914452
number of beta-spin electrons 5.000000019914452
Iteration 76: delta_dm = 1.310340e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019920495
number of beta-spin electrons 5.000000019920495
Iteration 77: delta_dm = 1.178613e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019925953
number of beta-spin electrons 5.000000019925953
Iteration 78: delta_dm = 1.060156e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019930946
number of beta-spin electrons 5.000000019930946
Iteration 79: delta_dm = 9.536274e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019935456
number of beta-spin electrons 5.000000019935456
Iteration 80: delta_dm = 8.578230e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019939547
number of beta-spin electrons 5.000000019939547
Iteration 81: delta_dm = 7.716603e-04
5.0 5.0 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000199432435
number of beta-spin electrons 5.0000000199432435
Iteration 82: delta_dm = 6.941665e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019946605
number of beta-spin electrons 5.000000019946605
Iteration 83: delta_dm = 6.244673e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019949654
number of beta-spin electrons 5.000000019949654
Iteration 84: delta_dm = 5.617769e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.0000000199524
number of beta-spin electrons 5.0000000199524
Iteration 85: delta_dm = 5.053890e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019954887
number of beta-spin electrons 5.000000019954887
Iteration 86: delta_dm = 4.546687e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019957148
number of beta-spin electrons 5.000000019957148
Iteration 87: delta_dm = 4.090452e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019959184
number of beta-spin electrons 5.000000019959184
Iteration 88: delta_dm = 3.680055e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019961023
number of beta-spin electrons 5.000000019961023
Iteration 89: delta_dm = 3.310880e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019962692
number of beta-spin electrons 5.000000019962692
Iteration 90: delta_dm = 2.978782e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019964195
number of beta-spin electrons 5.000000019964195
Iteration 91: delta_dm = 2.680030e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019965559
number of beta-spin electrons 5.000000019965559
Iteration 92: delta_dm = 2.411272e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001996679
number of beta-spin electrons 5.00000001996679
Iteration 93: delta_dm = 2.169491e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019967903
number of beta-spin electrons 5.000000019967903
Iteration 94: delta_dm = 1.951975e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019968902
number of beta-spin electrons 5.000000019968902
Iteration 95: delta_dm = 1.756288e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001996981
number of beta-spin electrons 5.00000001996981
Iteration 96: delta_dm = 1.580234e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.00000001997063
number of beta-spin electrons 5.00000001997063
Iteration 97: delta_dm = 1.421843e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019971373
number of beta-spin electrons 5.000000019971373
Iteration 98: delta_dm = 1.279339e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019972037
number of beta-spin electrons 5.000000019972037
Iteration 99: delta_dm = 1.151129e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019972637
number of beta-spin electrons 5.000000019972637
Iteration 100: delta_dm = 1.035776e-04
5.000000000000003 5.000000000000003 number of electrons in each spin channel computed from the density matrix
C:\Users\ayers\AppData\Local\Programs\Python\Python312\Lib\site-packages\grid\atomgrid.py:895: UserWarning: Lebedev weights are negative which can introduce round-off errors.
sphere_grid = AngularGrid(degree=deg_i, method=method)
number of alpha-spin electrons 5.000000019973171
number of beta-spin electrons 5.000000019973171
Iteration 101: delta_dm = 9.319901e-05
SCF converged in 101 iterations.
Final eigenvalues (orbital energies): [[-13.42794989 -13.42794989]
[ 0.32196893 0.32196893]
[ 0.65586839 0.65586839]
[ 1.12855171 1.12855171]
[ 1.18939881 1.18939881]
[ 1.39706767 1.39706767]
[ 1.5593528 1.5593528 ]]
Energy Computation#
After the SCF procedure has converged, we can compute the total electronic energy. If the electrons moved independently in the average potential, then the total energy would simply be the sum of the orbital energies. However, this double counts the electron-electron repulsion energy: each electron interacts with \(N-1\) other electrons, so summing over all the electrons counts each electron-electron interaction twice. We can correct for this by subtracting the spurious interaction terms from the sum of orbital energies or just directly compute the energy from the relevant terms in the Hamiltonian. I.e.,
\[\begin{split}
\begin{split}
E_{\text{Hartree}} &= \operatorname{Tr}[\hat{h}_{\text{core}}(D_\alpha + D_\beta)] + \tfrac{1}{2} \operatorname{Tr}[(D_\alpha + D_\beta)(J_\alpha + J_\beta)]\\
E_{\text{Fermi-Amaldi}} &= E_{\text{Hartree}} - \left( \tfrac{1}{N_{\alpha}} \operatorname{Tr}[D_\alpha J_\alpha] - \tfrac{1}{N_{\beta}} \operatorname{Tr}[D_\beta J_\beta] \right)\\
E_{\text{LSDA}} &= E_{\text{Hartree}} + \tfrac{4}{3}\int \rho(\mathbf{r}) v_{\text{xc}}(\mathbf{r}) d\mathbf{r}\\
E_{\text{HF}} &= E_{\text{Hartree}} - \frac{1}{2} \operatorname{Tr}[D_\alpha K_\alpha + D_\beta K_\beta]
\end{split}
\end{split}\]
def compute_total_energy(s, h_core, eri, density_matrix, atnums, atcoords, ao_basis, calc_type='hartree'):
"""Compute the total energy of the system given the density matrix and integrals.
Parameters
----------
s : (nao, nao) ndarray
AO overlap matrix.
h_core : (nao, nao) ndarray
AO core Hamiltonian matrix.
eri : (nao, nao, nao, nao) ndarray
Two-electron integrals in the AO basis.
density_matrix : (nao, nao) ndarray
Density matrix.
atnums : (natoms,) ndarray
Atomic numbers.
atcoords : (natoms, 3) ndarray
Atomic coordinates.
ao_basis : list of str
Basis set for each atom.
calc_type : str
Type of energy expression to evaluate.
Choices are 'fermi-amaldi', 'lsda', 'hartree', and 'hartree-fock'.
Returns
-------
float
Total energy of the system.
"""
D_a = density_matrix[:, :, 0]
D_b = density_matrix[:, :, 1]
J_a = build_j_ao(D_a, eri)
J_b = build_j_ao(D_b, eri)
# Compute the Hartree energy
energy = np.sum((D_a + D_b) * h_core) + 0.5 * np.sum((D_a + D_b) * (J_a + J_b))
if calc_type == 'hartree':
pass # The Hartree energy is just the sum of the core Hamiltonian and the Coulomb energy, which we have already computed.
elif calc_type == 'fermi-amaldi':
# Compute the number of electrons of each spin.
n_electrons_a = np.trace(D_a @ s)
n_electrons_b = np.trace(D_b @ s)
energy -= 1.0/n_electrons_a * np.sum(D_a * J_a) + 1.0/n_electrons_b * np.sum(D_b * J_b)
elif calc_type == 'hartree-fock':
K_a = build_k_ao(D_a, eri)
K_b = build_k_ao(D_b, eri)
energy -= 0.5 * np.sum(D_a * K_a) + 0.5 * np.sum(D_b * K_b)
elif calc_type == 'lsda':
F_xc = build_vxc_lsda(atnums, atcoords, ao_basis, density_matrix)
energy += np.sum((D_a * F_xc[:, :, 0] + D_b * F_xc[:, :, 1]))*4/3
else:
raise ValueError(f"Invalid calculation type: {calc_type}. Must be one of 'fermi-amaldi', 'lsda', 'hartree', or 'hartree-fock'.")
return energy
Electron Correlation#
We can approximate the effects of electron correlation by using either the variational principle or perturbation theory. In the variational approach, we use a linear combination of the Slater determinants that diagonalize the average-potential approximation. In perturbation theory, we consider the difference between the exact Hamiltonian and the average-potential Hamiltonian as a perturbation and compute corrections to the energy and wavefunction order by order in perturbation theory. The simplest perturbative correction is the Møller-Plesset second-order (MP2) correction, which accounts for the missing correlation energy in Hartree-Fock theory.
\[
E_{\mathrm{MP2}}^{\mathrm{UHF}}
=
\frac{1}{2}
\sum_{iajb}
\frac{(ia|jb)\left[(ia|jb)-(ib|ja)\right]}
{\epsilon_i^\alpha+\epsilon_j^\alpha-\epsilon_a^\alpha-\epsilon_b^\alpha}
+
\frac{1}{2}
\sum_{IAJB}
\frac{(IA|JB)\left[(IA|JB)-(IB|JA)\right]}
{\epsilon_I^\beta+\epsilon_J^\beta-\epsilon_A^\beta-\epsilon_B^\beta}
+
\sum_{iaJB}
\frac{(ia|JB)^2}
{\epsilon_i^\alpha+\epsilon_J^\beta-\epsilon_a^\alpha-\epsilon_B^\beta}.
\]
where the integrals are written in the molecular-orbital basis.
def ump2_energy(eps_a, eps_b, c_a, c_b, eri_ao, nalpha, nbeta):
"""
Compute the unrestricted MP2 correlation energy from AO-basis ERIs.
Parameters
----------
eps_a : (nmo_a,) ndarray
Alpha-spin orbital energies.
eps_b : (nmo_b,) ndarray
Beta-spin orbital energies.
c_a : (nao, nmo_a) ndarray
Alpha-spin MO coefficients.
c_b : (nao, nmo_b) ndarray
Beta-spin MO coefficients.
eri_ao : (nao, nao, nao, nao) ndarray
AO-basis ERIs in chemist's notation:
eri_ao[mu, nu, lam, sig] = (mu nu | lam sig)
nalpha : int
Number of occupied alpha orbitals.
nbeta : int
Number of occupied beta orbitals.
Returns
-------
emp2 : float
Total UMP2 correlation energy.
e_aa : float
Same-spin alpha-alpha contribution.
e_bb : float
Same-spin beta-beta contribution.
e_ab : float
Opposite-spin alpha-beta contribution.
Notes
-----
Uses the formulas
E_aa = 1/2 sum_{iajb} (ia|jb)[(ia|jb) - (ib|ja)] / D_ijab
E_bb = 1/2 sum_{IAJB} (IA|JB)[(IA|JB) - (IB|JA)] / D_IJAB
E_ab = sum_{iaJB} (ia|JB)^2 / D_iJaB
where
D_ijab = eps_i^a + eps_j^a - eps_a^a - eps_b^a
D_IJAB = eps_I^b + eps_J^b - eps_A^b - eps_B^b
D_iJaB = eps_i^a + eps_J^b - eps_a^a - eps_B^b
"""
eps_a = np.asarray(eps_a)
eps_b = np.asarray(eps_b)
c_a = np.asarray(c_a)
c_b = np.asarray(c_b)
eri_ao = np.asarray(eri_ao)
# Occupied / virtual partitions
c_a_occ = c_a[:, :nalpha]
c_a_vir = c_a[:, nalpha:]
c_b_occ = c_b[:, :nbeta]
c_b_vir = c_b[:, nbeta:]
eps_a_occ = eps_a[:nalpha]
eps_a_vir = eps_a[nalpha:]
eps_b_occ = eps_b[:nbeta]
eps_b_vir = eps_b[nbeta:]
# -------- alpha-alpha block: g_aa[i,a,j,b] = (i a | j b) --------
g_aa = np.einsum(
"mi,na,lj,sb,mnls->iajb",
c_a_occ, c_a_vir, c_a_occ, c_a_vir, eri_ao,
optimize=True
)
denom_aa = (
eps_a_occ[:, None, None, None]
+ eps_a_occ[None, None, :, None]
- eps_a_vir[None, :, None, None]
- eps_a_vir[None, None, None, :]
)
exch_aa = g_aa.transpose(0, 3, 2, 1) # (i b | j a)
e_aa = 0.5 * np.sum(g_aa * (g_aa - exch_aa) / denom_aa)
# -------- beta-beta block: g_bb[I,A,J,B] = (I A | J B) --------
g_bb = np.einsum(
"mI,nA,lJ,sB,mnls->IAJB",
c_b_occ, c_b_vir, c_b_occ, c_b_vir, eri_ao,
optimize=True
)
denom_bb = (
eps_b_occ[:, None, None, None]
+ eps_b_occ[None, None, :, None]
- eps_b_vir[None, :, None, None]
- eps_b_vir[None, None, None, :]
)
exch_bb = g_bb.transpose(0, 3, 2, 1) # (I B | J A)
e_bb = 0.5 * np.sum(g_bb * (g_bb - exch_bb) / denom_bb)
# -------- alpha-beta block: g_ab[i,a,J,B] = (i a | J B) --------
g_ab = np.einsum(
"mi,na,lJ,sB,mnls->iaJB",
c_a_occ, c_a_vir, c_b_occ, c_b_vir, eri_ao,
optimize=True
)
denom_ab = (
eps_a_occ[:, None, None, None]
+ eps_b_occ[None, None, :, None]
- eps_a_vir[None, :, None, None]
- eps_b_vir[None, None, None, :]
)
e_ab = np.sum(g_ab * g_ab / denom_ab)
emp2 = e_aa + e_bb + e_ab
return emp2
density_matrix = compute_density_matrix(eps, phi, n_electrons)
eps, phi = scf_cycle(h_core, density_matrix, eri, s, n_electrons, h2o.atnums, h2o.atcoords, basis, calc_type='hartree-fock', max_iter=200, tol=1e-4)
eps_a, eps_b = eps[:, 0], eps[:, 1]
c_a, c_b = phi[:, :, 0], phi[:, :, 1]
density_matrix = compute_density_matrix(eps, phi, n_electrons)
energy = compute_total_energy(s, h_core, eri, density_matrix, h2o.atnums, h2o.atcoords, basis, calc_type='hartree-fock')
print("Hartree-Fock total energy:", energy)
emp2 = ump2_energy(eps_a=eps_a,eps_b=eps_b,c_a=c_a,c_b=c_b,eri_ao=eri, nalpha=n_electrons[0], nbeta=n_electrons[1])
print("UMP2 correlation energy:", emp2)
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 1: delta_dm = 7.207171e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 2: delta_dm = 6.436001e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 3: delta_dm = 5.748130e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 4: delta_dm = 5.134604e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 5: delta_dm = 4.587424e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 6: delta_dm = 4.099444e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 7: delta_dm = 3.664282e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 8: delta_dm = 3.276240e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 9: delta_dm = 2.930228e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 10: delta_dm = 2.621702e-01
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 11: delta_dm = 2.346605e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 12: delta_dm = 2.101316e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 13: delta_dm = 1.882602e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 14: delta_dm = 1.687578e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 15: delta_dm = 1.513669e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 16: delta_dm = 1.358578e-01
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 17: delta_dm = 1.220254e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 18: delta_dm = 1.096870e-01
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 19: delta_dm = 9.867934e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 20: delta_dm = 8.885691e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 21: delta_dm = 8.008994e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 22: delta_dm = 7.226271e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 23: delta_dm = 6.527203e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 24: delta_dm = 5.902593e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 25: delta_dm = 5.344243e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 26: delta_dm = 4.844850e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 27: delta_dm = 4.397907e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 28: delta_dm = 3.997622e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 29: delta_dm = 3.638837e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 30: delta_dm = 3.316964e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 31: delta_dm = 3.027921e-02
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 32: delta_dm = 2.768081e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 33: delta_dm = 2.534221e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 34: delta_dm = 2.323479e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 35: delta_dm = 2.133315e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 36: delta_dm = 1.961476e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 37: delta_dm = 1.805965e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 38: delta_dm = 1.665013e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 39: delta_dm = 1.537054e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 40: delta_dm = 1.420699e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 41: delta_dm = 1.314722e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 42: delta_dm = 1.218036e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 43: delta_dm = 1.129679e-02
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 44: delta_dm = 1.048802e-02
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 45: delta_dm = 9.746505e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 46: delta_dm = 9.065576e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 47: delta_dm = 8.439322e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 48: delta_dm = 7.862495e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 49: delta_dm = 7.330438e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 50: delta_dm = 6.839008e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 51: delta_dm = 6.384516e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 52: delta_dm = 5.963674e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 53: delta_dm = 5.573542e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 54: delta_dm = 5.211491e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 55: delta_dm = 4.875161e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 56: delta_dm = 4.562434e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 57: delta_dm = 4.271401e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 58: delta_dm = 4.000339e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 59: delta_dm = 3.747691e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 60: delta_dm = 3.512046e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 61: delta_dm = 3.292120e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 62: delta_dm = 3.086747e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 63: delta_dm = 2.894861e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 64: delta_dm = 2.715491e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 65: delta_dm = 2.547745e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 66: delta_dm = 2.390806e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 67: delta_dm = 2.243922e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 68: delta_dm = 2.106402e-03
5.000000000000003 5.000000000000003 number of electrons in each spin channel computed from the density matrix
Iteration 69: delta_dm = 1.977609e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 70: delta_dm = 1.856953e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 71: delta_dm = 1.743891e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 72: delta_dm = 1.637919e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 73: delta_dm = 1.538569e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 74: delta_dm = 1.445408e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 75: delta_dm = 1.358034e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 76: delta_dm = 1.276072e-03
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 77: delta_dm = 1.199175e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 78: delta_dm = 1.127017e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 79: delta_dm = 1.059297e-03
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 80: delta_dm = 9.957336e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 81: delta_dm = 9.360632e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 82: delta_dm = 8.800406e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 83: delta_dm = 8.274367e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 84: delta_dm = 7.780373e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 85: delta_dm = 7.316423e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 86: delta_dm = 6.880645e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 87: delta_dm = 6.471287e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 88: delta_dm = 6.086712e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 89: delta_dm = 5.725386e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 90: delta_dm = 5.385873e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 91: delta_dm = 5.066829e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 92: delta_dm = 4.766995e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 93: delta_dm = 4.485191e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 94: delta_dm = 4.220311e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 95: delta_dm = 3.971320e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 96: delta_dm = 3.737247e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 97: delta_dm = 3.517181e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 98: delta_dm = 3.310269e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 99: delta_dm = 3.115710e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 100: delta_dm = 2.932754e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 101: delta_dm = 2.760696e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 102: delta_dm = 2.598877e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 103: delta_dm = 2.446677e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 104: delta_dm = 2.303514e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 105: delta_dm = 2.168842e-04
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 106: delta_dm = 2.042150e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 107: delta_dm = 1.922958e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 108: delta_dm = 1.810813e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 109: delta_dm = 1.705294e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 110: delta_dm = 1.606001e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 111: delta_dm = 1.512563e-04
5.000000000000002 5.000000000000002 number of electrons in each spin channel computed from the density matrix
Iteration 112: delta_dm = 1.424628e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 113: delta_dm = 1.341867e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 114: delta_dm = 1.263972e-04
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 115: delta_dm = 1.190653e-04
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 116: delta_dm = 1.121635e-04
5.0 5.0 number of electrons in each spin channel computed from the density matrix
Iteration 117: delta_dm = 1.056664e-04
5.000000000000001 5.000000000000001 number of electrons in each spin channel computed from the density matrix
Iteration 118: delta_dm = 9.954993e-05
SCF converged in 118 iterations.
Hartree-Fock total energy: -84.8335824085104
UMP2 correlation energy: -0.03539372403748492