Getting the all-electron density

The variational quantity of the PAW formalism is the pseudo-density \(\tilde{n}\). This is also the density returned by the get_pseudo_density() method of the GPAW calculator. Sometimes it is desirable to work with the true all-electron density. The PAW formalism offers a recipe for reconstructing the all-electron density from the pseudo-density, and in GPAW, this can be reached by the method get_all_electron_density() of the GPAW class:

get_all_electron_density(spin=None, gridrefinement=2, pad=True)

Return reconstructed all-electron density array.

The get_all_electron_density() method is used in the same way as you would normally use the get_pseudo_density() method, i.e.:

from gpaw import GPAW
from ase.build import molecule

calc = GPAW(mode='fd')
mol = molecule('C6H6', calculator=calc)
mol.center(vacuum=5)
E = mol.get_potential_energy()
nt = calc.get_pseudo_density()
n_ae = calc.get_all_electron_density()

would give you the pseudo-density in nt and the all-electron density in n_ae.

As the all-electron density has more structure than the pseudo-density, it is necessary to refine the density grid used to represent the pseudo-density. This can be done using the gridrefinement keyword of the get_all_electron_density method for n_ae_fine:

n_ae_fine = calc.get_all_electron_density(gridrefinement=2)

Current only the values 1, 2, and 4 are supported (2 is default).

The all-electron density will always integrate to the total number of electrons of the considered system (independent of the grid resolution), while the pseudo density will integrate to some more or less arbitrary number. This fact is illustrated in the following example.

See also

Bader

Example: NaCl

As an example of application, consider the three systems Na, Cl, and NaCl. The pseudo- and all-electron densities of these three systems can be calculated with the script NaCl.py:

# web-page: all_electron.csv
import numpy as np

from ase.build import molecule
from gpaw import GPAW
from ase.parallel import paropen

unitcell = np.array([6.5, 6.6, 9.])
gridrefinement = 2

f = paropen('all_electron.csv', 'w')

for formula in ('Na', 'Cl', 'NaCl',):
    if formula in ['Na', 'Cl']:
        hund = True
    else:
        hund = False
    calc = GPAW(mode='fd',
                xc='PBE',
                h=0.18,
                hund=hund,
                convergence={'eigenstates': 1e-8},
                txt=formula + '.txt')

    sys = molecule(formula, cell=unitcell, calculator=calc)
    sys.center()
    sys.get_potential_energy()

    # Get densities
    nt = calc.get_pseudo_density()
    n = calc.get_all_electron_density(gridrefinement=gridrefinement)

    # Get integrated values
    dv = sys.get_volume() / calc.get_number_of_grid_points().prod()
    It = nt.sum() * dv
    I = n.sum() * dv / gridrefinement**3
    print('%-4s,%4.2f,%5.2f' % (formula, It, I), file=f)

f.close()

The result for the integrated pseudo- and all-electron densities of the three systems is:

formula	ñ	n
Na	1.88	11.00
Cl	7.50	17.00
NaCl	9.36	28.00

From which we see that the all-electron densities integrate to the total number of electrons in the system, as expected.