From owner-chemistry@ccl.net Wed Mar 13 10:31:00 2024
From: "Szabolcs Goger szgoger^^gmail.com" <owner-chemistry[-]server.ccl.net>
To: CCL
Subject: CCL: Radial atomic DFT solver inconsistencies
Message-Id: <-55109-240313092047-1507-wsnC1oRmKSkmT5b0vLTfnQ[-]server.ccl.net>
X-Original-From: "Szabolcs  Goger" <szgoger]_[gmail.com>
Date: Wed, 13 Mar 2024 09:20:45 -0400


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.

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 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.

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.

Best regards,
Szabolcs Goger

import pyscf
> from pyscf.scf import atom_hf, atom_ks
import numpy as np

mol = pyscf.M(
    atom = 'H 0 0 0',
    basis = 'augccpvtz',
    symmetry = True,
spin = 1,
verbose=6
)

myhf = mol.KS()
myhf.xc = 'PBE'
myhf.kernel()

sphericalHF = atom_ks.get_atm_nrks(mol, xc='PBE')