pyCADMium: Chemical Atoms in Diatomic Molecules. A prolate spheroidal Python module for embedding calculations

Most practical calculations use a finite basis set of fixed functional form to represent operators and quantities like potentials, orbitals, and densities (Hill, 2013; Szabo & Ostlund, 2012). Such approaches require analytic expressions for the matrix representation of various operators and often exhibit poor basis set convergence. The basis set incompleteness error can be minimized by increasing the number of basis functions used, but it cannot be entirely eliminated in practice because of the lack of parameters (for instance, exponents, contraction coefficients, etc.) for large basis sets. On the other hand, grid-based methods intrinsically allow for an accurate representation of operators and permit the use of arbitrarily large basis sets (Andrade et al., 2015). Nevertheless, the number of points required to achieve a significant accuracy can quickly become unmanageable. pyCADMium presented in this work, uses a prolate spheroidal (PS) grid to circumvent these issues for atoms and diatomic molecules.


Introduction
Diatomic molecules are among the most useful systems to test new ideas in quantum chemistry: • They highlight the present and ubiquitous problems of modern approximations (Cohen et al., 2008). • By defining two fragments with an atom each, diatomics are ideal to implement and test quantum embedding methods (Chávez & Wasserman, 2020).
Most practical calculations use a finite basis set of fixed functional form to represent operators and quantities like potentials, orbitals, and densities (Hill, 2013;Szabo & Ostlund, 2012). Such approaches require analytic expressions for the matrix representation of various operators and often exhibit poor basis set convergence. The basis set incompleteness error can be minimized by increasing the number of basis functions used, but it cannot be entirely eliminated in practice because of the lack of parameters (for instance, exponents, contraction coefficients, etc.) for large basis sets. On the other hand, grid-based methods intrinsically allow for an accurate representation of operators and permit the use of arbitrarily large basis sets (Andrade et al., 2015). Nevertheless, the number of points required to achieve a significant accuracy can quickly become unmanageable. pyCADMium presented in this work, uses a prolate spheroidal (PS) grid to circumvent these issues for atoms and diatomic molecules.

Summary
In this work, we introduce pyCADMium, a Python module that uses a PS coordinate grid to accurately perform computational chemistry calculations on systems with cylindrical symmetry. The name is an acronym for Chemical Atoms in Diatomic Molecules. pyCADMium originated in a proprietary programming language but has been rewritten from the ground up as an open-source alternative. The code has been the main driving force behind the development of "Partition Density Functional Theory," (Chávez & Wasserman, 2020;Elliott et al., 2010;Jiang et al., 2018;Nafziger et al., 2017) a method that uses quantum embedding to lower the cost of a calculation and fixes problems inherent to density functional approximations (Cohen et al., 2008).
The PS grid is a coordinate system formed by revolving an elliptic coordinate plane around the line intersecting the two foci. These planes are formed by ellipses and hyperbolae that share the same focus (Arfken & Weber, 1999). Atoms and diatomics are ideally suited for this coordinate system given that each foci can be used to allocate an atom. Additionally, the PS grid is denser around the foci in its cartesian representation, where we expect wave functions of molecular systems to change rapidly (Ryabinkin et al., 2017).

Functionality
Consider a PS grid with foci located at (− /2, 0) and ( /2, 0), where represents the internuclear separation in a diatomic of interest. Place one (A) or two atoms (A, B) at each of these foci and specify its charge, number of electrons, and number of atomic and/or molecular orbitals. Continue by specifying the symmetry of orbitals to calculate.
Once the set of fragments and/or molecule is defined pyCADMium can perform: 1. Density functional theory (DFT) calculation. Choose a density functional approximation up to the generalized gradient approximation (GGA) rung (Tao et al., 2003) from the library of exchange correlation functionals, libxc (Marques et al., 2012), and perform a self-consistent Kohn-Sham DFT calculation to find the energies, orbitals and density.
2. Density-to-Potential inversion. Given any density on the PS grid, perform a numerical inversion (Shi et al., 2022) to find the multiplicative external potential that reproduces the input density (Wu & Yang, 2003). This problem is ill-posed when the potential is expressed on a basis-set, but it is well-posed when done on a grid (Jensen & Wasserman, 2018).
3. Partition-DFT. Given a molecule of interest, fragment its external potential and find the set of densities associated with them. The sum of these densities, for most cases, will not be equal to the density of the full system, but we want them to be. Perform a numerical inversion to find the non-additive kinetic potential, that when added to the exact nonadditive external-hartree-exchange-correlation potential becomes the unique embedding potential required for the sum of isolated fragments to reproduce the molecular density as well as minimize the sum of fragment energies (Elliott et al., 2010).

Statement of Need
PS coordinates have proven to be accurate in calculations using atoms and diatomic molecules (Becke, 1982). Despite these coordinates being used thoroughly in literature, the options for freely-available modules that focus on embedding applications are almost non-existent. One example worth mentioning is DARSEC (Makmal et al., 2009), a FORTRAN-based code that uses the same prolate spheroidal methodology in the context of the optimized effective potential equation. Their use of the coordinates and basic operators, such as derivatives and integration, is analogous to our approach. A key difference is our use of libxc to obtain the functional derivatives needed to obtain exchange-correlation potentials. Moreover, their codebase is not openly available.
Our code was designed from the ground up with embedding calculations in mind. As described in Nafziger & Wasserman (2014), it can be used to perform and quickly develop other embedding calculations. Last but not least, a repository of Jupyter Notebooks is available on GitHub that includes various examples of the functionalities available in the code.