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

1# SPDX-FileCopyrightText: 2023 The eminus developers 

2# SPDX-License-Identifier: Apache-2.0 

3"""GTH file handling.""" 

4 

5import importlib.resources 

6import pathlib 

7 

8from eminus import backend as xp 

9from eminus.logger import log 

10 

11 

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

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

14 

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

16 

17 Args: 

18 atom: Atom name. 

19 

20 Keyword Args: 

21 charge: Valence charge. 

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

23 

24 Returns: 

25 GTH parameters. 

26 """ 

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

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

29 else: 

30 file_path = pathlib.Path(psp_path) 

31 

32 if charge is not None: 

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

34 else: 

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

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

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

38 try: 

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

40 except IndexError: 

41 if atom != "X": 

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

43 return mock_gth() 

44 if len(files) > 1: 

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

46 

47 psp = {} 

48 cloc = xp.zeros(4) 

49 rp = xp.zeros(4) 

50 Nproj_l = xp.zeros(4, dtype=int) 

51 h = xp.zeros((4, 3, 3)) 

52 try: 

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

54 atom = fh.readline() 

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

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

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

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

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

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

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

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

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

64 cloc[i] = float(val) 

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

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

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

68 for i in range(lmax): 

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

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

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

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

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

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

75 for k, val in enumerate(proj): 

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

77 # Copy upper triangle elements to the lower triangle 

78 for jtmp in range(3): 

79 for ktmp in range(i, 3): 

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

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

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

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

84 except FileNotFoundError: 

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

86 return mock_gth() 

87 return psp 

88 

89 

90def mock_gth(): 

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

92 

93 Returns: 

94 GTH parameters (all zero). 

95 """ 

96 return { 

97 "Zion": 0, 

98 "rloc": 0, 

99 "cloc": xp.zeros(4), 

100 "lmax": 0, 

101 "rp": xp.zeros(4), 

102 "Nproj_l": xp.zeros(4, dtype=int), 

103 "h": xp.zeros((4, 3, 3)), 

104 }