Coverage for eminus/io/gth.py: 93.44%

1# SPDX-FileCopyrightText: 2021 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""GTH file handling.""" 

4 

5import importlib.resources 

6import pathlib 

7 

8import numpy as np 

9 

10from ..logger import log 

11 

12 

13def read_gth(atom, charge=None, psp_path="pbe"): 

14 """Read GTH files for a given atom. 

15 

16 Reference: Phys. Rev. B 54, 1703. 

17 

18 Args: 

19 atom: Atom name. 

20 

21 Keyword Args: 

22 charge: Valence charge. 

23 psp_path: Path to GTH pseudopotential files. Defaults to installation_path/psp/pbe. 

24 

25 Returns: 

26 GTH parameters. 

27 """ 

28 if psp_path in {"pade", "pbe"}: 

29 # The try-except block can be removed once eminus drops Python 3.7 support 

30 try: 

31 file_path = importlib.resources.files(f"eminus.psp.{psp_path}") 

32 except AttributeError: 

33 with importlib.resources.path("eminus.psp", psp_path) as p: 

34 file_path = p 

35 else: 

36 file_path = pathlib.Path(psp_path) 

37 

38 if charge is not None: 

39 f_psp = file_path.joinpath(f"{atom}-q{charge}") 

40 else: 

41 # When sorting normally, e.g., Ga-q13 would come before Ga-q3 

42 # Fix this by sorting over the charge integer numbers using the sorted key argument 

43 files = sorted(file_path.glob(f"{atom}-q*"), key=lambda s: int(str(s).split("-q")[-1])) 

44 try: 

45 f_psp = pathlib.Path(files[0]) 

46 except IndexError: 

47 if atom != "X": 

48 log.warning(f'There is no GTH pseudopotential in "{file_path}" for "{atom}".') 

49 return mock_gth() 

50 if len(files) > 1: 

51 log.info(f'Multiple pseudopotentials found for "{atom}". Continue with "{f_psp.name}".') 

52 

53 psp = {} 

54 cloc = np.zeros(4) 

55 rp = np.zeros(4) 

56 Nproj_l = np.zeros(4, dtype=int) 

57 h = np.zeros((4, 3, 3)) 

58 try: 

59 with open(f_psp, encoding="utf-8") as fh: 

60 atom = fh.readline() 

61 # Skip the first line, since we don't need the atom symbol here. If needed, use 

62 # psp["atom"] = atom.split()[0] # Atom symbol 

63 N_all = fh.readline().split() 

64 psp["Zion"] = sum(int(N) for N in N_all) # Ionic charge 

65 loc = fh.readline().split() 

66 psp["rloc"] = float(loc[0]) # Range of local Gaussian charge distribution 

67 # Skip the number of local coefficients, since we don't need it. If needed, use 

68 # psp["n_c_local"] = int(loc[1]) # Number of local coefficients 

69 for i, val in enumerate(loc[2:]): 

70 cloc[i] = float(val) 

71 psp["cloc"] = cloc # Coefficients for the local part 

72 lmax = int(fh.readline().split()[0]) 

73 psp["lmax"] = lmax # Maximal angular momentum in the non-local part 

74 for i in range(lmax): 

75 proj = fh.readline().split() 

76 rp[i], Nproj_l[i] = float(proj[0]), int(proj[1]) 

77 for k, val in enumerate(proj[2:]): 

78 h[i, 0, k] = float(val) 

79 for j in range(1, Nproj_l[i]): 

80 proj = fh.readline().split() 

81 for k, val in enumerate(proj): 

82 h[i, j, j + k] = float(val) 

83 # Copy upper triangle elements to the lower triangle 

84 for jtmp in range(3): 

85 for ktmp in range(i, 3): 

86 h[i, ktmp, jtmp] = h[i, jtmp, ktmp] 

87 psp["rp"] = rp # Projector radius for each angular momentum 

88 psp["Nproj_l"] = Nproj_l # Number of non-local projectors 

89 psp["h"] = h # Projector coupling coefficients per AM channel 

90 except FileNotFoundError: 

91 log.warning(f'There is no GTH pseudopotential for "{f_psp.name}".') 

92 return mock_gth() 

93 return psp 

94 

95 

96def mock_gth(): 

97 """Create a mock dictionary with all-zeros for atom species with no pseudopotential file. 

98 

99 Returns: 

100 GTH parameters (all zero). 

101 """ 

102 return { 

103 "Zion": 0, 

104 "rloc": 0, 

105 "cloc": np.zeros(4), 

106 "lmax": 0, 

107 "rp": np.zeros(4), 

108 "Nproj_l": np.zeros(4, dtype=int), 

109 "h": np.zeros((4, 3, 3)), 

110 }