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.
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
> 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:_:
helsinki.fi
Associate professor, computational
chemistry
------------------------------------------------------------------------------
Susi Lehtola, FT
akatemiatutkija
susi.lehtola:_:
helsinki.fi
dosentti, laskennallinen kemia
------------------------------------------------------------------------------
-= This is automatically added to each message by the mailing script
=-
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.
If your mail bounces from CCL with 5.7.1
error, check: