Coverage for eminus/potentials.py: 92.86%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""Collection of miscellaneous potentials.""" 

4 

5import numpy as np 

6from scipy.linalg import norm 

7 

8from .gth import init_gth_loc 

9from .logger import log 

10 

11 

12def harmonic(scf): 

13 """Harmonic potential. 

14 

15 Can be used for quantum dot calculations. 

16 

17 Args: 

18 scf: SCF object. 

19 

20 Returns: 

21 Harmonic potential in reciprocal space. 

22 """ 

23 atoms = scf.atoms 

24 freq = 2 

25 dr = norm(atoms.r - np.sum(atoms.a, axis=1) / 2, axis=1) 

26 Vharm = 0.5 * freq**2 * dr**2 

27 return atoms.Jdag(atoms.O(atoms.J(Vharm))) 

28 

29 

30def coulomb(scf): 

31 """All-electron Coulomb potential. 

32 

33 Reference: Bull. Lebedev Phys. Inst. 42, 329. 

34 

35 Args: 

36 scf: SCF object. 

37 

38 Returns: 

39 Coulomb potential in reciprocal space. 

40 """ 

41 atoms = scf.atoms 

42 Z = atoms.Z[0] # This potential should only be used for same species 

43 

44 # Ignore the division by zero for the first elements 

45 # One could do some proper indexing with [1:] but indexing is slow 

46 with np.errstate(divide="ignore", invalid="ignore"): 

47 Vcoul = -4 * np.pi * Z / atoms.G2 

48 Vcoul[0] = 0 

49 

50 Sf = np.sum(atoms.Sf, axis=0) 

51 return atoms.J(Vcoul * Sf) 

52 

53 

54def ge(scf): 

55 """Starkloff-Joannopoulos local pseudopotential for germanium. 

56 

57 Fourier-transformed by Tomas Arias. 

58 

59 Reference: Phys. Rev. B 16, 5212. 

60 

61 Args: 

62 scf: SCF object. 

63 

64 Returns: 

65 Germanium pseudopotential in reciprocal space. 

66 """ 

67 atoms = scf.atoms 

68 Z = 4 # This potential should only be used for germanium 

69 lamda = 18.5 

70 rc = 1.052 

71 Gm = np.sqrt(atoms.G2) 

72 Vps = np.empty_like(atoms.G2) 

73 

74 with np.errstate(divide="ignore", invalid="ignore"): 

75 Vps = ( 

76 -2 

77 * np.pi 

78 * np.exp(-np.pi * Gm / lamda) 

79 * np.cos(rc * Gm) 

80 * (Gm / lamda) 

81 / (1 - np.exp(-2 * np.pi * Gm / lamda)) 

82 ) 

83 for n in range(5): 

84 Vps = Vps + (-1) ** n * np.exp(-lamda * rc * n) / (1 + (n * lamda / Gm) ** 2) 

85 Vps = Vps * 4 * np.pi * Z / Gm**2 * (1 + np.exp(-lamda * rc)) - 4 * np.pi * Z / Gm**2 

86 

87 # Special case for G=(0,0,0) 

88 n = np.arange(1, 5) 

89 Vps[0] = ( 

90 4 

91 * np.pi 

92 * Z 

93 * (1 + np.exp(-lamda * rc)) 

94 * ( 

95 rc**2 / 2 

96 + 1 / lamda**2 * (np.pi**2 / 6 + np.sum((-1) ** n * np.exp(-lamda * rc * n) / n**2)) 

97 ) 

98 ) 

99 

100 Sf = np.sum(atoms.Sf, axis=0) 

101 return atoms.J(Vps * Sf) 

102 

103 

104def init_pot(scf): 

105 """Handle and initialize potentials. 

106 

107 Args: 

108 scf: SCF object. 

109 

110 Returns: 

111 Potential in reciprocal space. 

112 """ 

113 try: 

114 pot = IMPLEMENTED[scf.pot](scf) 

115 except KeyError: 

116 log.exception(f'No potential found for "{scf.pot}".') 

117 raise 

118 return pot 

119 

120 

121#: Map potential names with their respective implementation. 

122IMPLEMENTED = { 

123 "harmonic": harmonic, 

124 "coulomb": coulomb, 

125 "ge": ge, 

126 "gth": init_gth_loc, 

127}