Coverage for gpaw/elf.py: 86%

1"""This module defines an ELF function.""" 

2from __future__ import annotations 

3 

4import sys 

5 

6import numpy as np 

7from gpaw.core import UGArray 

8from gpaw.fd_operators import Gradient 

9from gpaw.new.ase_interface import GPAW, ASECalculator 

10from gpaw.new.calculation import DFTCalculation 

11 

12 

13def elf(nt_sR: np.ndarray, 

14 nt_grad2_sR: np.ndarray, 

15 taut_sR: np.ndarray, 

16 ncut: float | None = None) -> np.ndarray: 

17 """Pseudo electron localisation function (ELF). 

18 

19 See: 

20 

21 Becke and Edgecombe, J. Chem. Phys., vol 92 (1990) 5397 

22 

23 More comprehensive definition in 

24 M. Kohout and A. Savin, Int. J. Quantum Chem., vol 60 (1996) 875-882 

25 

26 Parameters 

27 ========== 

28 nt_sR: 

29 Pseudo valence density. 

30 nt_grad2_sR: 

31 Squared norm of the density gradient. 

32 taut_sR: 

33 Kinetic energy density. 

34 ncut: 

35 Minimum density cutoff parameter. 

36 

37 Returns 

38 ======= 

39 np.ndarray: 

40 Array of ELF values. 

41 """ 

42 

43 # Fermi constant 

44 cF = 3.0 / 10 * (3 * np.pi**2)**(2 / 3) 

45 

46 eps = 1e-11 

47 nt_sR = nt_sR.copy() 

48 nt_sR[nt_sR < eps] = eps 

49 

50 if nt_sR.shape[0] == 2: 

51 # Kouhut eq. (9) 

52 D0 = 2**(2 / 3) * cF * (nt_sR[0]**(5 / 3) + 

53 nt_sR[1]**(5 / 3)) 

54 

55 taut = taut_sR.sum(axis=0) 

56 D = taut - (nt_grad2_sR[0] / nt_sR[0] + nt_grad2_sR[1] / nt_sR[1]) / 8 

57 else: 

58 # Kouhut eq. (7) 

59 D0 = cF * nt_sR[0]**(5 / 3) 

60 taut = taut_sR[0] 

61 D = taut - nt_grad2_sR[0] / nt_sR[0] / 8 

62 

63 elf_R = 1.0 / (1.0 + (D / D0)**2) 

64 

65 if ncut is not None: 

66 nt = nt_sR.sum(axis=0) 

67 elf_R[nt < ncut] = 0.0 

68 

69 return elf_R 

70 

71 

72def elf_from_dft_calculation(dft: DFTCalculation | ASECalculator, 

73 ncut: float = 1e-6) -> UGArray: 

74 """Calculate the electronic localization function. 

75 

76 Parameters 

77 ========== 

78 dft: 

79 DFT-calculation object. 

80 ncut: 

81 Density cutoff below which the ELF is zero. 

82 

83 Returns 

84 ======= 

85 UGArray: 

86 ELF values. 

87 """ 

88 if isinstance(dft, ASECalculator): 

89 dft = dft.dft 

90 density = dft.density 

91 density.update_ked(dft.ibzwfs) 

92 taut_sR = density.taut_sR 

93 assert taut_sR is not None 

94 nt_sR = density.nt_sR 

95 grad_v = [Gradient(nt_sR.desc._gd, v, n=2) for v in range(3)] 

96 gradnt2_sR = nt_sR.new(zeroed=True) 

97 for gradnt2_R, nt_R in zip(gradnt2_sR, nt_sR): 

98 for grad in grad_v: 

99 gradnt_R = grad(nt_R) 

100 gradnt2_R.data += gradnt_R.data**2 

101 elf_R = nt_sR.desc.empty() 

102 elf_R.data[:] = elf( 

103 nt_sR.data, gradnt2_sR.data, taut_sR.data, ncut) 

104 return elf_R 

105 

106 

107if __name__ == '__main__': 

108 e_R = elf_from_dft_calculation(GPAW(sys.argv[1]).dft, 0.001) 

109 e_R.isosurface(isomin=0.8, isomax=0.8)