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, h=0.05, n=2)[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.
h: 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.center(vacuum=2.5)
hli.calc = GPAW(txt='hli.txt', mode='fd')
hli.get_potential_energy()

# Transformer:
t = PS2AE(hli.calc, h=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,
             label=r'$\tilde\psi_{}$'.format(n))
    plt.plot(x, ae[i, i], '-', color=color,
             label=r'$\psi_{}$'.format(n))

    # 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}$]')
plt.legend()
plt.savefig('hli-wfs.png')
hli.calc.write('hli.gpw')
../../_images/hli-wfs.png
PS2AE.get_wave_function(n, k=0, s=0, ae=True)[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)

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

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

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

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

plt.plot(x, 0 * x, 'k')
plt.xlabel('z [Ang]')
plt.ylabel('potential [eV]')
plt.ylim(ymin=-100, ymax=10)
plt.legend()
plt.savefig('hli-pot.png')
../../_images/hli-pot.png
PS2AE.get_electrostatic_potential(ae=True, rcgauss=0.02)[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.