Source code for gpaw.utilities.ekin

from math import pi
from typing import Tuple
import numpy as np
from ase.units import Ha

from gpaw.typing import Array1D
from gpaw.setup import Setup


[docs]def ekin(dataset: Setup) -> Tuple[Array1D, Array1D, float]: """Calculate PAW kinetic energy contribution as a function of G.""" ds = dataset rgd = dataset.rgd de_j = ds.data.e_kin_jj.diagonal() phit_j = ds.pseudo_partial_waves_j e0 = -ds.Kc e_k: Array1D = 0.0 # type: ignore for f, l, de, phit in zip(ds.f_j, ds.l_j, de_j, phit_j): if f == 0.0: continue phit_r = np.array([phit(r) for r in rgd.r_g]) G_k, phit_k = rgd.fft(phit_r * rgd.r_g**(l + 1), l) e_k += f * 0.5 * phit_k**2 * G_k**4 / (2 * pi)**3 e0 -= f * de return G_k, e_k, e0
def dekindecut(G: Array1D, de: Array1D, ecut: float) -> float: """Linear interpolation.""" dG = G[1] G0 = (2 * ecut)**0.5 g = int(G0 / dG) dedG = np.polyval(np.polyfit(G[g:g + 2], de[g:g + 2], 1), G0) dedecut = dedG / G0 return dedecut if __name__ == '__main__': import argparse import matplotlib.pyplot as plt from gpaw.setup import create_setup parser = argparse.ArgumentParser( description='Calculate approximation to the energy variation with ' 'plane-wave cutoff energy. The approximation is to use the kinetic ' 'energy from a PAW atom, which can be calculated efficiently on ' 'a radial grid.') parser.add_argument('-d', '--derivative', type=float, metavar='ECUT', help='Calculate derivative of energy correction with ' 'respect to plane-wave cutoff energy.') parser.add_argument('name', help='Name of PAW dataset.') args = parser.parse_args() ds = create_setup(args.name) G, de, e0 = ekin(ds) dG = G[1] if args.derivative: dedecut = -dekindecut(G, de, args.derivative / Ha) print('de/decut({}, {} eV) = {:.6f}' .format(args.name, args.derivative, dedecut)) else: de = (np.add.accumulate(de) - 0.5 * de[0] - 0.5 * de) * dG ecut = 0.5 * G**2 * Ha y = (de[-1] - de) * Ha plt.plot(ecut, y) plt.xlim(300, 1000) plt.ylim(0, y[ecut > 300][0]) plt.show()