Obtaining all-electron wave functions and electrostatic potential

Wave functions

To get from the pseudo (PS) wave function to the all-electron (AE) wave function, a PAW correction needs to be added. The correction will add the cusp and all the wiggles necessary for the wave function to be orthogonal to all the frozen core states. This tutorial show how to use the PS2AE class for interpolating the PS wave function to a fine grid where the PAW corrections can be represented:

class gpaw.utilities.ps2ae.PS2AE(calc: gpaw.calculator.GPAW, grid_spacing: float = 0.05, n: int = 2, h=None)[source]

Transform PS to AE wave functions.

Interpolates PS wave functions to a fine grid and adds PAW corrections in order to obtain true AE wave functions.

Create transformation object.

calc: GPAW calculator object

The calcalator that has the wave functions.

grid_spacing: float

Desired grid-spacing in Angstrom.

n: int

Force number of points to be a mulitiple of n.

Here is the code for plotting some AE wave functions for a HLi dimer using a PAW dataset for Li with a frozen 1s orbital (get_wave_function()):

import matplotlib.pyplot as plt
from ase import Atoms
from ase.units import Bohr
from gpaw.utilities.ps2ae import PS2AE
from gpaw import GPAW

hli = Atoms('HLi', positions=[[0, 0, 0], [0, 0, 1.6]])
hli.calc = GPAW(txt='hli.txt', mode='fd')

# Transformer:
t = PS2AE(hli.calc, grid_spacing=0.05)

for n, color in enumerate(['green', 'red']):
    ps = t.get_wave_function(n, ae=False)
    ae = t.get_wave_function(n)
    norm = t.gd.integrate(ae**2) * Bohr**3
    print('Norm:', norm)
    assert abs(norm - 1) < 1e-2
    i = ps.shape[0] // 2
    x = t.gd.coords(2) * Bohr

    # Interpolated PS and AE wfs:
    plt.plot(x, ps[i, i], '--', color=color,
    plt.plot(x, ae[i, i], '-', color=color,

    # Raw PS wfs:
    ps0 = hli.calc.get_pseudo_wave_function(n, pad=True)
    gd = hli.calc.wfs.gd
    i = ps0.shape[0] // 2
    X = gd.coords(2) * Bohr
    plt.plot(X, ps0[i, i], 'o', color=color)

plt.plot(x, 0 * x, 'k')
plt.xlabel('z [Ang]')
plt.ylabel('wave functions [Ang$^{-3/2}$]')
PS2AE.get_wave_function(n: int, k: int = 0, s: int = 0, ae: bool = True)gpaw.hints.ArrayND[source]

Interpolate wave function.

n: int

Band index.

k: int

K-point index.

s: int

Spin index.

ae: bool

Add PAW correction to get an all-electron wave function.

Electrostatic potential

The relevant formulas can be found here: Note on electrostatic potential.

Here is how to extract the AE potential from a gpw-file using the get_electrostatic_potential() method:

import matplotlib.pyplot as plt
from ase.units import Bohr
from gpaw.utilities.ps2ae import PS2AE
from gpaw import GPAW

calc = GPAW('hli.gpw', txt=None)

# Avarage PS potentials:
vh, vli = calc.get_atomic_electrostatic_potentials()
zh, zli = calc.atoms.positions[:, 2]
rh, rli = [max(setup.rcut_j) * Bohr for setup in calc.setups]
plt.plot([zh - rh, zh + rh], [vh, vh], label=r'$\tilde v$(H)')
plt.plot([zli - rli, zli + rli], [vli, vli], label=r'$\tilde v$(Li)')

# Transformer:
t = PS2AE(calc, grid_spacing=0.05)

# Interpolated PS and AE potentials:
ps = t.get_electrostatic_potential(ae=False)
ae = t.get_electrostatic_potential()
i = ps.shape[0] // 2
z = t.gd.coords(2) * Bohr

plt.plot(z, ps[i, i], '--', label=r'$\tilde v$')
plt.plot(z, ae[i, i], '-', label=r'$v$')

# Raw PS potential:
ps0 = calc.get_electrostatic_potential()
gd = calc.hamiltonian.finegd
i = ps0.shape[0] // 2
Z = gd.coords(2) * Bohr
plt.plot(Z, ps0[i, i], 'o')

plt.plot(z, 0 * z, 'k')
plt.xlabel('z [Ang]')
plt.ylabel('potential [eV]')
plt.ylim(bottom=-100, top=10)

The figure also shows the avarage PS potentials at the atomic sites calculated with the get_atomic_electrostatic_potentials() method.

PS2AE.get_electrostatic_potential(ae: bool = True, rcgauss: float = 0.02)gpaw.hints.ArrayND[source]

Interpolate electrostatic potential.

Return value in eV.

ae: bool

Add PAW correction to get the all-electron potential.

rcgauss: float

Width of gaussian (in Angstrom) used to represent the nuclear charge.

Pseudo density


PS2AE.get_pseudo_density(add_compensation_charges: bool = True)gpaw.hints.ArrayND[source]

Interpolate pseudo density.