From owner-chemistry@ccl.net Mon Mar 18 16:03:01 2024 From: "=?UTF-8?Q?Szabolcs_G=C3=B3ger?= szgoger{:}gmail.com" To: CCL Subject: CCL: Radial atomic DFT solver inconsistencies Message-Id: <-55113-240318055228-19326-LEI0vKDHeOGSAfdpYC3dXQ[#]server.ccl.net> X-Original-From: =?UTF-8?Q?Szabolcs_G=C3=B3ger?= Content-Type: multipart/alternative; boundary="000000000000e6ab6b0613ec4fbc" Date: Mon, 18 Mar 2024 10:51:37 +0100 MIME-Version: 1.0 Sent to CCL by: =?UTF-8?Q?Szabolcs_G=C3=B3ger?= [szgoger]~[gmail.com] --000000000000e6ab6b0613ec4fbc Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable Dear Community, I must also say that this was an amazing response, and I have a lot of catching-up to do before being able to adequately address the issue. I want to add one thing, however: it was mentioned that "the radial 1s orbital turns out to be quite similar between the two [spin-restricted and spin-unrestricted PBE hydrogenic] calculations." My view is actually quite the opposite: depending on *what* expectation values we are interested, we might get significantly different results. For example, even though the density maxima (defined from \rho * r^2) occur at quite similar regions, integrals of the densities can be quite different. Apart from the -very different - expectation values of the energy, the matrix elements <0|r^3|0> also happen to not match, being 8.7 a.u. vs. 10.5 a.u. - this is exactly the difference that is mentioned in the Q-Chem manuals. This, to me, hints that the overall shape of the density is quite different. Apart from the fundamental questions of the proper definition of atomic solvers, I am also worried whether it is a problem that different definitions exist in different software. Some dispersion corrections relying on free atomic data already use empirical factors to correct for the mismatch of the different definitions, which is not an elegant solution at all, but at least the problem is addressed. I can't help to wonder whether there are other points of divergence that most people do not know about... Best regards, Szabolcs On Fri, Mar 15, 2024 at 7:09=E2=80=AFPM David Mannock dmannock-#-yahoo.com = < owner-chemistry(-)ccl.net> wrote: > 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! Thanks for providing such a thoughtful & detailed reply on > this subject! > > David Mannock > > On Thursday, March 14, 2024 at 10:27:00 a.m. MDT, Susi Lehtola > susi.lehtola###alumni.helsinki.fi wrote: > > > > Sent to CCL 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: "Szabolcs Goger" [szgoger\a/gmail.com] Dear Community, > > > > I am reaching out with hopes of gathering some thoughts on an issue tha= t > 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 solve= r > > specificall utilizing the fact that only a radial Schrdinger equation i= s > 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-Che= m > 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 witho= ut > spherical symmetry. > > In my opinion, since neither the one-determinant Hartree-Fock (which is n= ot > 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 analyse= s > 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 fo= r > 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 i= n > arbitrary symmetries. The library has already been interfaced with Psi4 i= n > 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=BCcke= l > 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 philosophie= s > > 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 i= s > 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 solver= s > 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 multipl= e > > 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; unrestrictin= g > 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 high= er > 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. 15= 9, > 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 fo= r > 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:_:helsinki.fi Associate professor, computational > chemistry > http://susilehtola.github.io/ University of Helsinki, Finland > > -------------------------------------------------------------------------= ----- > Susi Lehtola, FT akatemiatutkija > susi.lehtola:_:helsinki.fi dosentti, laskennallinen kemia > http://susilehtola.github.io/ Helsingin yliopisto > > -------------------------------------------------------------------------= ----- > > > > -=3D This is automatically added to each message by the mailing script = =3D-> the strange characters on the top line to the - sign. You can also> > E-mail to subscribers: CHEMISTRY-=C3=8Cl.net or use:> > E-mail to administrators: CHEMISTRY-REQUEST-=C3=8Cl.net or use> > > --000000000000e6ab6b0613ec4fbc Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable
Dear Community,

I must also = say that this was an amazing response, and I have a lot of catching-up to d= o before being able to adequately address the issue.

I want to add one thing, however: it was mentioned that "the radial= 1s orbital turns out to be quite similar between the two [spin-restricted = and spin-unrestricted PBE hydrogenic] calculations." My view is actual= ly quite the opposite: depending on what expectation values we are i= nterested, we might get significantly different results. For example, even = though the density maxima (defined from \rho * r^2) occur at quite similar = regions, integrals of the densities can be quite different. Apart from the = -very different - expectation values of the energy, the matrix elements <= ;0|r^3|0> also happen to not match, being 8.7 a.u. vs. 10.5 a.u. - this = is exactly the difference that is mentioned in the Q-Chem manuals. This, to= me, hints that the overall shape of the density is quite different.
<= div>
Apart from the fundamental questions of the proper defin= ition of atomic solvers, I am also worried whether it is a problem that dif= ferent definitions exist in different software. Some dispersion corrections= relying on free atomic data already use empirical factors to correct for t= he mismatch of the different definitions, which is not an elegant solution = at all, but at least the problem is addressed. I can't help to wonder w= hether there are other points of divergence that most people do not know ab= out...

Best regards,
Szabolcs
<= /div>
O= n Fri, Mar 15, 2024 at 7:09=E2=80=AFPM David Mannock dmannock-#-yahoo.com <owner-chemistry(-)ccl.net> wrote:
Susi,

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

David Mannock
=

=20
=20
On Thursday, March 14, 2024 at 10:27:00 a.m. MDT, Susi = Lehtola susi.lehtola###alumni.helsinki.fi <owner-chemistry-=C3=8Cl.net> wrote:



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: "Szabolcs=C2=A0 Goger" [szgoger= \a/gmail.com] Dear Commu= nity,
>
> I am r= eaching out with hopes of gathering some thoughts on an issue that has
<= /div>
> 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 Hirs= hfeld
> partitioning, which uses free atomic d= ensities as reference. As far as I can
> tell,= there are two approaches of treating free atoms: one can use the main
<= /div>
> SCF solver without any major adjustment, or impl= ement a different solver
> specificall utilizi= ng the fact that only a radial Schrdinger equation is to
> be solved for a free atom. An example for the former philosop= hy is Q-Chem, an
> example for the latter is F= HI-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 spherica= l symmetry for most atoms in the periodic table. The Q-Chem
strategy also doesn't really work for many atoms in the peri= odic table, since
the self-consistent field calcu= lations can be very hard to converge without
sphe= rical 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 approximat= ions lead to correct atomic multiplet structure
(= the resulting energy depends on the magnetic quantum number M), the best wa= y to
implement atomic solvers to set up molecular= guesses or Hirshfeld analyses is to
do the spher= ical 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 interactio= n but
only fractional exact exchange. I have rece= ntly discussed the theory and
implementation with= in a finite-element program in

S. Lehtola, Fully numerical calculations on atoms with fractional = occupations
and range-separated exchange function= als, Phys. Rev. A 101, 012516 (2020).
doi:10.1103= /PhysRevA.101.012516 arXiv:1908.02528

<= div dir=3D"ltr">S. Lehtola, Meta-GGA density functional calculations on ato= ms with spherically
symmetric densities in the fi= nite element formalism, J. Chem. Theory Comput. 19,
2502 (2023). doi:10.1021/acs.jctc.3c00183 arXiv:2302.06284

The first of these works also prese= nted the minimal-energy occupations for
spherical= ly symmetric atoms, for use in the superposition of atomic densities
(SAD) guess; these wave functions can also be used to s= et up Hirshfeld
decompositions.

Having a general molecular Fock builder= , it is a bit tedious but straightforward
to solv= e the atomic radial-only problem; this is what is done in PySCF. I am
also about to publish a reusable open source orbital o= ptimization 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 t= o handle an arbitrary number of types of particles in
arbitrary symmetries. The library has already been interfaced with Psi= 4 in order
to run spherically symmetric atomic SC= F calculations for forming SAD guesses.
The SAD c= alculations also allow an easy way to set up extended H=C3=BCckel guesses,<= br>
as I discussed in

<= /div>
S. Lehtola, Assessment of initial guesses for self-co= nsistent field
calculations. Superposition of Ato= mic Potentials: simple yet efficient, J. Chem.
Th= eory Comput. 15, 1593 (2019). doi:10.1021/acs.jctc.8b01089. arXiv:1810.1165= 9

> The issue is th= e following: properties obtained with a DFT calculation are
> different depending on which implementation is used. For ex= ample, in PySCF,
> the PBE energy of a triplet= carbon atom with aVTZ basis set is -37.7942 Eh
&= gt; with the 3D solver and -37.74511 Eh with the radial one. Not only the e= nergy
> is different, but expectation values o= btained from the wavefunction too. This
> is c= uriously documented in the Q-Chem manual
> ht= tps://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
> g= iving consistently different results. Both atomic volumes can also be
> obtained by PySCF, depending whether the radial o= r 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*, an= d the electron density is symmetrized afterwards, which formally
<= div dir=3D"ltr">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<= br>
fractional occupations.

> There are some more pieces to the puzzle: = for a hydrogen atom, the PBE
> results are hig= hly different with the two solvers, but with HF, there is no
> difference for the hydrogen. This led me to believe that t= here is a
> difference on how the XC functiona= l should be defined for 3D and the quasi-1D
> = solver. For noble gas atoms, however, the difference between the solvers wi= th
> 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&= #39;t think multiple
> codes would get it wron= g 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 cal= culation; however, the radial 1s
orbital turns ou= t to be quite similar between the two calculations. For
instance, in HelFEM I get the total energies -0.458928597 vs -0.49= 9990369 Eh for
the spin-restricted vs spin-unrest= ricted calculations of hydrogen, with a
density m= aximum at 1.026772e+00 vs 9.827961e-01 bohr.

=
For Hartree-Fock, you should observe a larger differ= ence: -0.357709999 vs
-0.500000000 Eh, and 1.1382= 90e+00 vs 1.000000e+00 bohr.

Note that PySCF only implements the spin-restricted variant; unres= tricting the
spin will yield you a lower energy. = If you allow the radial orbitals to become
spin-p= olarized, 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.11= 03/PhysRevA.101.012516 cited above, the F atom is
already a good example of symmetry breaking. To get the lowest energy, a b= asis
set of s and p functions is insufficient; yo= u also need to introduce higher
angular momenta l= ike 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 fract= ional 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 r= eference volumes, but
> maybe there is somethi= ng 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 r= ecent work on the (lack of)
reproducibility of de= nsity functional approximations

S. Lehtola and M. A. L. Marques, Reproducibility of density funct= ional
approximations: how new functionals should = be reported, J. Chem. Phys. 159,
114116 (2023). d= oi:10.1063/5.0167763 arXiv:2307.07474

<= div dir=3D"ltr">which may also cause differences between programs. The PBE = functional, however,
is relatively well standardi= zed between programs, and you should get the same
results between PySCF and Q-Chem, if you use the Libxc implementation in P= ySCF.
As we discuss in the above article, by defa= ult, 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 i= s used in most implementations of PBE.

=
Hope this helps,

Susi
--
--= ---------------------------------------------------------------------------= -
Mr. Susi Lehtola, PhD=C2=A0 =C2=A0 =C2=A0 =C2= =A0 =C2=A0 =C2=A0 Academy of Finland research fellow
susi.lehtola:_:helsin= ki.fi=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 Associate professor, computatio= nal chemistry
http://susilehtola.github.io/=C2=A0 =C2=A0 Un= iversity of Helsinki, Finland
-------------------= -----------------------------------------------------------
Susi Lehtola, FT=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0= =C2=A0 =C2=A0 akatemiatutkija
susi.lehtola:_:helsinki.fi=C2=A0 =C2=A0 = =C2=A0 =C2=A0 =C2=A0 dosentti, laskennallinen kemia
http://susil= ehtola.github.io/=C2=A0 =C2=A0 Helsingin yliopisto
------------------------------------------------------------------= ------------



-=3D This is automatically add= ed to each message by the mailing script =3D-
To = recover the email address of the author of the message, please change
the strange characters on the top line to the - sign. = You can also
look up the X-Original-From: line in= the mail header.

E-ma= il to subscribers: C= HEMISTRY-=C3=8Cl.net or use:

E-mail to administrators: CHEMISTRY-REQUE= ST-=C3=8Cl.net or use
=C2=A0 =C2=A0 =C2=A0 http://www.ccl.net/cgi-bin/ccl/send_ccl_message

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



If your mail bounces fro= m CCL with 5.7.1 error, check:
=C2=A0 =C2=A0 =C2= =A0 http://ww= w.ccl.net/spammers.txt



--000000000000e6ab6b0613ec4fbc--