Coverage for gpaw/utilities/ekin.py: 18%

1from math import pi 

2from typing import Tuple 

3import numpy as np 

4from ase.units import Ha 

5 

6from gpaw.typing import Array1D 

7from gpaw.setup import Setup 

8 

9 

10def ekin(dataset: Setup) -> Tuple[Array1D, Array1D, float]: 

11 """Calculate PAW kinetic energy contribution as a function of G.""" 

12 ds = dataset 

13 rgd = dataset.rgd 

14 de_j = ds.data.e_kin_jj.diagonal() 

15 phit_j = ds.pseudo_partial_waves_j 

16 e0 = -ds.Kc 

17 e_k: Array1D = 0.0 # type: ignore 

18 

19 for f, l, de, phit in zip(ds.f_j, ds.l_j, de_j, phit_j): 

20 if f == 0.0: 

21 continue 

22 phit_r = np.array([phit(r) for r in rgd.r_g]) 

23 G_k, phit_k = rgd.fft(phit_r * rgd.r_g**(l + 1), l) 

24 e_k += f * 0.5 * phit_k**2 * G_k**4 / (2 * pi)**3 

25 e0 -= f * de 

26 

27 return G_k, e_k, e0 

28 

29 

30def dekindecut(G: Array1D, de: Array1D, ecut: float) -> float: 

31 """Linear interpolation.""" 

32 dG = G[1] 

33 G0 = (2 * ecut)**0.5 

34 g = int(G0 / dG) 

35 dedG = np.polyval(np.polyfit(G[g:g + 2], de[g:g + 2], 1), G0) 

36 dedecut = dedG / G0 

37 return dedecut 

38 

39 

40if __name__ == '__main__': 

41 import argparse 

42 import matplotlib.pyplot as plt 

43 from gpaw.setup import create_setup 

44 

45 parser = argparse.ArgumentParser( 

46 description='Calculate approximation to the energy variation with ' 

47 'plane-wave cutoff energy. The approximation is to use the kinetic ' 

48 'energy from a PAW atom, which can be calculated efficiently on ' 

49 'a radial grid.') 

50 parser.add_argument('-d', '--derivative', type=float, metavar='ECUT', 

51 help='Calculate derivative of energy correction with ' 

52 'respect to plane-wave cutoff energy.') 

53 parser.add_argument('name', help='Name of PAW dataset.') 

54 args = parser.parse_args() 

55 

56 ds = create_setup(args.name) 

57 

58 G, de, e0 = ekin(ds) 

59 dG = G[1] 

60 

61 if args.derivative: 

62 dedecut = -dekindecut(G, de, args.derivative / Ha) 

63 print('de/decut({}, {} eV) = {:.6f}' 

64 .format(args.name, args.derivative, dedecut)) 

65 else: 

66 de = (np.add.accumulate(de) - 0.5 * de[0] - 0.5 * de) * dG 

67 

68 ecut = 0.5 * G**2 * Ha 

69 y = (de[-1] - de) * Ha 

70 ax = plt.figure().add_subplot(111) 

71 ax.plot(ecut, y) 

72 ax.set_xlim(300.0, 1000.0) 

73 ax.set_ylim(0.0, y[ecut > 300.0][0]) 

74 plt.show()