From owner-chemistry@ccl.net Thu Mar 14 11:51:01 2024 From: "Susi Lehtola susi.lehtola###alumni.helsinki.fi" To: CCL Subject: CCL: Radial atomic DFT solver inconsistencies Message-Id: <-55110-240314053623-16690-RMDKsvbuYqWf6UVNznTi1Q]^[server.ccl.net> X-Original-From: Susi Lehtola Content-Language: en-US Content-Transfer-Encoding: 8bit Content-Type: text/plain; charset=UTF-8 Date: Thu, 14 Mar 2024 11:36:06 +0200 MIME-Version: 1.0 Sent to CCL by: Susi Lehtola [susi.lehtola=-=alumni.helsinki.fi] On 3/13/24 15:20, Szabolcs Goger szgoger^^gmail.com wrote: > > Sent to CCL by: "Szabolcs Goger" [szgoger\a/gmail.com] Dear Community, > > I am reaching out with hopes of gathering some thoughts on an issue that has > been troubling me for a bit. > > In particular, I've spent some time looking at free-atom solvers in different > quantum chemistry software, mainly because I am interested in Hirshfeld > partitioning, which uses free atomic densities as reference. As far as I can > tell, there are two approaches of treating free atoms: one can use the main > SCF solver without any major adjustment, or implement a different solver > specificall utilizing the fact that only a radial Schrdinger equation is to > be solved for a free atom. An example for the former philosophy is Q-Chem, an > example for the latter is FHI-Aims, with PySCF in an interesting spot since > either approach can be used relatively easily. Using normal SCF leads to incorrect results, since the atomic electron density breaks spherical symmetry for most atoms in the periodic table. The Q-Chem strategy also doesn't really work for many atoms in the periodic table, since the self-consistent field calculations can be very hard to converge without spherical symmetry. In my opinion, since neither the one-determinant Hartree-Fock (which is not really Hartree-Fock as it is defined in the literature!) or the Kohn-Sham density functional approximations lead to correct atomic multiplet structure (the resulting energy depends on the magnetic quantum number M), the best way to implement atomic solvers to set up molecular guesses or Hirshfeld analyses is to do the spherical averaging, which works nicely in DFT, and is indeed the standard way atomic calculations are done in the DFT community. However, if you have Hartree-Fock exchange, you will also introduce ghost interactions since the orbitals will experience full Coulomb interaction but only fractional exact exchange. I have recently discussed the theory and implementation within a finite-element program in S. Lehtola, Fully numerical calculations on atoms with fractional occupations and range-separated exchange functionals, Phys. Rev. A 101, 012516 (2020). doi:10.1103/PhysRevA.101.012516 arXiv:1908.02528 S. Lehtola, Meta-GGA density functional calculations on atoms with spherically symmetric densities in the finite element formalism, J. Chem. Theory Comput. 19, 2502 (2023). doi:10.1021/acs.jctc.3c00183 arXiv:2302.06284 The first of these works also presented the minimal-energy occupations for spherically symmetric atoms, for use in the superposition of atomic densities (SAD) guess; these wave functions can also be used to set up Hirshfeld decompositions. Having a general molecular Fock builder, it is a bit tedious but straightforward to solve the atomic radial-only problem; this is what is done in PySCF. I am also about to publish a reusable open source orbital optimization library called OpenOrbitalOptimizer, which can easily be used to implement spherically symmetric atomic calculations in any atomic-orbital program. https://github.com/susilehtola/OpenOrbitalOptimizer The library is able to handle an arbitrary number of types of particles in arbitrary symmetries. The library has already been interfaced with Psi4 in order to run spherically symmetric atomic SCF calculations for forming SAD guesses. The SAD calculations also allow an easy way to set up extended Hückel guesses, as I discussed in S. Lehtola, Assessment of initial guesses for self-consistent field calculations. Superposition of Atomic Potentials: simple yet efficient, J. Chem. Theory Comput. 15, 1593 (2019). doi:10.1021/acs.jctc.8b01089. arXiv:1810.11659 > The issue is the following: properties obtained with a DFT calculation are > different depending on which implementation is used. For example, in PySCF, > the PBE energy of a triplet carbon atom with aVTZ basis set is -37.7942 Eh > with the 3D solver and -37.74511 Eh with the radial one. Not only the energy > is different, but expectation values obtained from the wavefunction too. This > is curiously documented in the Q-Chem manual > https://manual.q-chem.com/latest/sect_TS.html on the example of the "volume" > of a free hydrogen atom, with codes designed with different philosophies > giving consistently different results. Both atomic volumes can also be > obtained by PySCF, depending whether the radial or the 3D solver is called. There are three things that affect this result. According to the manual, Q-Chem runs a *spin-unrestricted* calculation *without symmetry*, and the electron density is symmetrized afterwards, which formally yields fractional occupations of the orbitals. In contrast, FHI-aims and Quantum Espresso run *spin-restricted* calculations *with full atomic symmetry*, that is, the electron density is optimized for the fractional occupations. > There are some more pieces to the puzzle: for a hydrogen atom, the PBE > results are highly different with the two solvers, but with HF, there is no > difference for the hydrogen. This led me to believe that there is a > difference on how the XC functional should be defined for 3D and the quasi-1D > solver. For noble gas atoms, however, the difference between the solvers with > PBE is quite small, so the problem _might_ arise from a combination of > spin-orbitals with reduced-dimensional XC functionals. I would be tempted to > think of an issue regarding the spin indices, but I don't think multiple > codes would get it wrong in a reproducible way. What are you comparing here? For hydrogen, there is a big difference in total energy for doing a spin-restricted vs a spin-unrestricted PBE calculation; however, the radial 1s orbital turns out to be quite similar between the two calculations. For instance, in HelFEM I get the total energies -0.458928597 vs -0.499990369 Eh for the spin-restricted vs spin-unrestricted calculations of hydrogen, with a density maximum at 1.026772e+00 vs 9.827961e-01 bohr. For Hartree-Fock, you should observe a larger difference: -0.357709999 vs -0.500000000 Eh, and 1.138290e+00 vs 1.000000e+00 bohr. Note that PySCF only implements the spin-restricted variant; unrestricting the spin will yield you a lower energy. If you allow the radial orbitals to become spin-polarized, this will further lower the energy in the general case. And, if you allow the orbitals to break spatial symmetries, you arrive to the lowest possible energy. As I discussed in doi:10.1103/PhysRevA.101.012516 cited above, the F atom is already a good example of symmetry breaking. To get the lowest energy, a basis set of s and p functions is insufficient; you also need to introduce higher angular momenta like d, f, g, h etc to fully converge the total energy; this is one of the main reasons I think one should explicitly rely on the use of explicit spherical symmetry and fractional occupations. > However, this is where my skills end. This is why I am turning to the > community: is someone aware of this difference between 3D and radial atomic > solvers? If yes, would this cause any major issues down the line somewhere? I > know this discrepancy creates a headache for Hirshfeld reference volumes, but > maybe there is something else. > > For reference, a quick PySCF code to obtain the energy with the two solvers > is inserted at the end. I would like to draw attention here to our recent work on the (lack of) reproducibility of density functional approximations S. Lehtola and M. A. L. Marques, Reproducibility of density functional approximations: how new functionals should be reported, J. Chem. Phys. 159, 114116 (2023). doi:10.1063/5.0167763 arXiv:2307.07474 which may also cause differences between programs. The PBE functional, however, is relatively well standardized between programs, and you should get the same results between PySCF and Q-Chem, if you use the Libxc implementation in PySCF. As we discuss in the above article, by default, XCFun has a non-standard definition of PBE, simply because the PBE article two different values for the constant, while a third value was used in Burke's reference implementation, and this is also the value that is used in most implementations of PBE. Hope this helps, Susi -- ------------------------------------------------------------------------------ Mr. Susi Lehtola, PhD Academy of Finland research fellow susi.lehtola\a/helsinki.fi Associate professor, computational chemistry http://susilehtola.github.io/ University of Helsinki, Finland ------------------------------------------------------------------------------ Susi Lehtola, FT akatemiatutkija susi.lehtola\a/helsinki.fi dosentti, laskennallinen kemia http://susilehtola.github.io/ Helsingin yliopisto ------------------------------------------------------------------------------