Susi,

Wow!= This is an amazing answer covering many aspects of a difficult problem. As= a non-expert, I was very impressed and clearly have some reading to do! Th= anks for providing such a thoughtful & detailed reply on this subject!<= /div>

David Mannock

=20

On Thursday, March 14, 2024 at 10:27:00 a.m. MDT, Susi = Lehtola susi.lehtola###alumni.helsinki.fi <owner-chemistry * ccl.net> w= rote:

Sent to CC= L by: Susi Lehtola [susi.lehtola=3D-=3Dalumni.helsinki.fi]

On 3/13/24 15:20, Szabolcs Goger szgoger^^gmail.com wrote:

> Sent to CCL by: "S= zabolcs 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-a= tom solvers in different

> quantum chemistry s= oftware, mainly because I am interested in Hirshfeld

> partitioning, which uses free atomic densities as reference. As fa= r as I can

> tell, there are two approaches of= treating free atoms: one can use the main

> S= CF solver without any major adjustment, or implement a different solver
=

> specificall utilizing the fact that only a radi= al 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 int= eresting spot since

> either approach can be u= sed relatively easily.

Using normal SCF leads to incorrect results, since the atomic electron den= sity

breaks spherical symmetry for most atoms in = the periodic table. The Q-Chem

strategy also does= n't really work for many atoms in the periodic table, since

the self-consistent field calculations can be very hard to conve= rge without

spherical symmetry.

In my opinion, since neither the one-de= terminant Hartree-Fock (which is not

really Hartr= ee-Fock as it is defined in the literature!) or the Kohn-Sham

density functional approximations lead to correct atomic multi= plet structure

(the resulting energy depends on t= he magnetic quantum number M), the best way to

im= plement atomic solvers to set up molecular guesses or Hirshfeld analyses is= to

do the spherical averaging, which works nicel= y in DFT, and is indeed the

standard way atomic c= alculations are done in the DFT community.

However, if you have Hartree-Fock exchange, you will a= lso introduce ghost

interactions since the orbita= ls will experience full Coulomb interaction but

o= nly fractional exact exchange. I have recently discussed the theory and
=

implementation within a finite-element program in

S. Lehtola, Fully numeri= cal calculations on atoms with fractional occupations

and range-separated exchange functionals, Phys. Rev. A 101, 012516 (20= 20).

doi:10.1103/PhysRevA.101.012516 arXiv:1908.0= 2528

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.j= ctc.3c00183 arXiv:2302.06284

The first of these works also presented the minimal-energy occupat= ions 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 strai= ghtforward

to solve the atomic radial-only proble= m; this is what is done in PySCF. I am

also about= to publish a reusable open source orbital optimization library called
<= /div>

OpenOrbitalOptimizer, which can easily be used to imp= lement spherically

symmetric atomic calculations = in any atomic-orbital program.

https://git= hub.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=C3=BCckel guesses,

as I = discussed in

S. Lehtol= a, Assessment of initial guesses for self-consistent field

calculations. Superposition of Atomic Potentials: simple yet effi= cient, J. Chem.

Theory Comput. 15, 1593 (2019). d= oi:10.1021/acs.jctc.8b01089. arXiv:1810.11659

> The issue is the following: properties obtaine= d with a DFT calculation are

> different depen= ding on which implementation is used. For example, in PySCF,

> the PBE energy of a triplet carbon atom with aVTZ basis se= t 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 t= oo. This

> is curiously documented in the Q-Ch= em manual

> https://manual.q-chem.com/latest/= sect_TS.html on the example of the "volume"

&= gt; of a free hydrogen atom, with codes designed with different philosophie= s

> giving consistently different results. Bot= h atomic volumes can also be

> obtained by PyS= CF, depending whether the radial or the 3D solver is called.

There are three things that affect t= his result.

According = to the manual, Q-Chem runs a *spin-unrestricted* calculation *without

symmetry*, and the electron density is symmetrized aft= erwards, which formally

yields fractional occupat= ions of the orbitals.

= In contrast, FHI-aims and Quantum Espresso run *spin-restricted* calculatio= ns

*with full atomic symmetry*, that is, the elec= tron density is optimized for the

fractional occu= pations.

> There ar= e some more pieces to the puzzle: for a hydrogen atom, the PBE

> results are highly different with the two solvers, but w= ith HF, there is no

> difference for the hydro= gen. 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 rega= rding 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 differ= ence 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 calcula= tions. For

instance, in HelFEM I get the total en= ergies -0.458928597 vs -0.499990369 Eh for

the sp= in-restricted vs spin-unrestricted calculations of hydrogen, with a

density maximum at 1.026772e+00 vs 9.827961e-01 bohr.

For Hartree-Fock, you sh= ould 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 spi= n-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 i= s

already a good example of symmetry breaking. To= get the lowest energy, a basis

set of s and p fu= nctions is insufficient; you also need to introduce higher

angular momenta like d, f, g, h etc to fully converge the total e= nergy; this is

one of the main reasons I think on= e should explicitly rely on the use of

explicit s= pherical symmetry and fractional occupations.

> However, this is where my skills end. This is = why I am turning to the

> community: is someon= e aware of this difference between 3D and radial atomic

> solvers? If yes, would this cause any major issues down the l= ine somewhere? I

> know this discrepancy creat= es a headache for Hirshfeld reference volumes, but

> maybe there is something else.

> For reference, a quick PySCF code to obtain the e= nergy with the two solvers

> is inserted at th= e 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, Re= producibility of density functional

approximation= s: how new functionals should be reported, J. Chem. Phys. 159,

114116 (2023). doi:10.1063/5.0167763 arXiv:2307.07474

which may also cause differen= ces 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

thi= s is also the value that is used in most implementations of PBE.

<= div dir=3D"ltr">

Hope this helps,

Susi

---------------------------------------------------= ---------------------------

Mr. Susi Lehtola, PhD= Academy of Finland research fell= ow

susi.lehtola:_:helsinki.fi = Associate professor, computational chemistry

http= ://susilehtola.github.io/ University of Helsinki, Finland=

------------------------------------------------= ------------------------------

Susi Lehtola, FT&n= bsp; akatemiatutkij= a

susi.lehtola:_:helsinki.fi = dosentti, laskennallinen kemia

http://susilehtola.= github.io/ Helsingin yliopisto

= ---------------------------------------------------------------------------= ---

-=3D This is automatically added to eac= h message by the mailing script =3D-

To recover t= he email address of the author of the message, please change

the strange characters on the top line to the * sign. You can a= lso

look up the X-Original-From: line in the mail= header.

E-mail to sub= scribers: CHEMISTRY * ccl.net or use:

&nb= sp; http://www.ccl.net/cgi-bin/ccl/send_ccl_message

E-mail to administrators: CHEMISTRY-REQUEST * ccl.net or use

&n= bsp; http://www.ccl.net/cgi-bin/ccl/send_ccl_message

http://www.ccl.net/chemistr= y/sub_unsub.shtml

= Before posting, check wait time at: http://www.ccl.net

Job: htt= p://www.ccl.net/jobs

Conferences: http://server.ccl.net/chemistry/announcements/conferences/

Search Messages: http://= www.ccl.net/chemistry/searchccl/index.shtml

<= br>

If your mail bounces from CCL with 5.7.1 error, c= heck:

http://www.ccl.net/spammers.txt=

RTFI: http://w= ww.ccl.net/chemistry/aboutccl/instructions/

<= br>