☰
Navigation :

Kohn-Sham inversion

Thanks to the mapping between the electronic density and the Kohn-Sham potential it is possible to compute the exchange-correlation (xc) potential from the sole knowledge of the density. The usual SCF cycle of a ground-state calculation aims at finding the ground state density by iteratively updating the Hartree and xc potential. However, it is also possible to fix the density and ask the question: what is the xc potential that would produce this density, knowing the external potential acting on the electrons. The procedure that finds this xc potential is known as the Kohn-Sham inversion. This has been extensively used to learn about the exact xc potential for systems in which the exact density is known or where higher theory level can produce high-quality reference densities.

In this tutorial, we will explore how to perform the Kohn-Sham inversion for the simplest case: a one dipensional system with two electrons occupying the same orbital. In this case, the analytical expression of the xc potential is known. For more than two electrons, one needs to use iterative schemes, that we won’t discuss in this tutorial.

Input for the reference density

As usual, we will need to write an input file describing the system we want to calculate. In order to perform the Kohm-Sham inversion, we need first a density. Here, we will use Hartree-Fock theory to evaluate the ground-state of a 1D soft-Coulomb hydrogen atom.


CalculationMode = gs
FromScratch = yes

Dimensions = 1

BoxShape = parallelepiped
%Lsize
 15
%
Spacing = 0.05

%Species
'He1D' | species_soft_coulomb | softening | 1 | valence | 2
%

%Coordinates
'He1D' | 0.0
%

TheoryLevel = hartree_fock
AdaptivelyCompressedExchange = yes
ExperimentalFeatures = yes

%Output
 density | axis_x
%

ConvRelDens = 1e-8
EigensolverMaxIter = 50

Run this input file. This will produce the file ‘‘static/density.y=0,z=0’’. This is the density we will use in order to perform the Kohn-Sham inversion and find the local xc potential that produces it.

Rename the folder ‘‘static’’ to ‘‘static_ref’’, such that it does not get overwritten by the next run.

Input for the Kohn-Sham inversion

The Kohn-Shame inversion is done using the following input file.


CalculationMode = invert_ks
FromScratch = yes

ExperimentalFeatures = yes

Dimensions = 1

BoxShape = parallelepiped
%Lsize
 15
%
Spacing = 0.05

%Species
'He1D' | species_soft_coulomb | softening | 1 | valence | 2
%

%Coordinates
'He1D' | 0.0
%

XCFunctional = ks_inversion

InvertKSTargetDensity = "static_ref/density.y=0,z=0"
InvertKSmethod = two_particle

%Output
 potential | axis_x
 density | axis_x
%

We first note that here we change the CalculationMode to be ‘‘invert_ks’’ instead of gs. Other changes are mandatory. We need to set XCFunctional to ‘‘ks_inversion’’. Then, we need to instruct the code where to find the target density. This is specified thanks to the variable InvertKSTargetDensity. Importantly, whereas here we produced the target density using Octopus, with the same grid, there is no requirement on the grid in which the data is provided. If the grid is different, Octopus will perform internally an interpolation to obtain the density on the desired grid.

Finally, we need to choose the methodology used for the Kohn-Sham inversion. As in the present case we are considering the two-electron in one orbital case, we have set InvertKSmethod to ‘’two_particle’’.

Now run this input file and plot the potential. Does it looks like the expected result?

Comparison with the slater potential

In the first step of this tutorial, we performed a Hartree-Fock calculation in order to produce the target density. This is a nonlocal theory and in general the only way to obtain the corresponding local potential would be via the solution of the OEP equation. This can be done by Octopus. However, in our case, as we have a single orbital occupied, we know the analytical solution. This is the Slater potential.

We can therefore perform a Koh-Sham DFT calculation for our system using the following the input file.


CalculationMode = gs
FromScratch = yes

Dimensions = 1

BoxShape = parallelepiped
%Lsize
 15
%
Spacing = 0.05

%Species
'He1D' | species_soft_coulomb | softening | 1 | valence | 2
%

%Coordinates
'He1D' | 0.0   
%

XCFunctional = slater_x

%Output
 potential | axis_x
 density | axis_x
%

ConvRelDens = 1e-8
EigensolverMaxIter = 50

This calculation produces the same eigenvalue and density as the original HF calculation. We can finally compare the Slater potential and the one obtained from the Kohn-Sham inversion. There are few points to note:

These two points are related to the treatment of the asymptotic part of the potential and are limitations of the current implementation. To prevent numerical noise, the current implementation sets to zero the points where the density is too samll. As the xc potential can be defined up to a constant, this freedom is used to avoid discontinuity of the potential where we set it to zero.

By shifting the result of the Kohn-Sham inversion, we find that this matches perfectly the Slater potential. However, one can employ the fact that the interaction is a soft-Coulomb potential. This allows for an analytical treatment of the asymptotic region. This is done using the variable KSInversionAsymptotics to ‘‘xc_asymptotics_sc’’. If you rerun the Kohn-Sham inversion with this method, you should now properly reproduce perfectly the Slater potential.