Coverage for larch/xrd/cif2feff.py: 6%

185 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-10-16 21:04 +0000

1import os 

2from random import Random 

3from io import StringIO 

4 

5from xraydb import atomic_symbol, atomic_number, xray_edge 

6from larch.utils import fix_varname, strict_ascii, gformat 

7 

8from .amcsd_utils import PMG_CIF_OPTS, CifParser, Molecule, SpacegroupAnalyzer 

9 

10rng = Random() 

11 

12def get_atom_map(structure): 

13 """generalization of pymatgen atom map 

14 Returns: 

15 dict of ipots 

16 """ 

17 unique_pot_atoms = [] 

18 all_sites = [] 

19 for site in structure: 

20 for elem in site.species.elements: 

21 if elem.symbol not in unique_pot_atoms: 

22 unique_pot_atoms.append(elem.symbol) 

23 

24 atom_map = {} 

25 for i, atom in enumerate(unique_pot_atoms): 

26 atom_map[atom] = i + 1 

27 return atom_map 

28 

29 

30def read_cif_structure(ciftext): 

31 """read CIF text, return CIF Structure 

32 

33 Arguments 

34 --------- 

35 ciftext (string): text of CIF file or name of the CIF file. 

36 

37 Returns 

38 ------- 

39 pymatgen Structure object 

40 """ 

41 if CifParser is None: 

42 raise ValueError("CifParser from pymatgen not available. Try 'pip install pymatgen'.") 

43 try: 

44 cifstructs = CifParser(StringIO(ciftext), **PMG_CIF_OPTS) 

45 parse_ok = True 

46 file_found = True 

47 except: 

48 parse_ok = False 

49 file_found = False 

50 if os.path.exists(ciftext): 

51 file_found = True 

52 try: 

53 cifstructs = CifParser(ciftext, **PMG_CIF_OPTS) 

54 parse_ok = True 

55 except: 

56 parse_ok = False 

57 

58 try: 

59 cstruct = cifstructs.parse_structures()[0] 

60 except: 

61 raise ValueError('could not get structure from text of CIF file') 

62 

63 if not parse_ok: 

64 if not file_found: 

65 raise FileNotFoundError(f'file {ciftext:s} not found') 

66 else: 

67 raise ValueError('invalid text of CIF file') 

68 return cstruct 

69 

70def cif_sites(ciftext, absorber=None): 

71 "return list of sites for the structure" 

72 cstruct = read_cif_structure(ciftext) 

73 out = cstruct.sites 

74 if absorber is not None: 

75 abname = absorber.lower() 

76 out = [] 

77 for site in cstruct.sites: 

78 species = site.species_string.lower() 

79 if ',' in species and ':' in species: # multi-occupancy site 

80 for siteocc in species.split(','): 

81 sname, occ = siteocc.split(':') 

82 if sname.strip() == abname: 

83 out.append(site) 

84 elif species == abname: 

85 out.append(site) 

86 if len(out) == 0: 

87 out = cstruct.sites[0] 

88 return out 

89 

90 

91def cif2feffinp(ciftext, absorber, edge=None, cluster_size=8.0, absorber_site=1, 

92 site_index=None, extra_titles=None, with_h=False, 

93 version8=True, rng_seed=None): 

94 """convert CIF text to Feff8 or Feff6l input file 

95 

96 Arguments 

97 --------- 

98 ciftext (string): text of CIF file or name of the CIF file. 

99 absorber (string or int): atomic symbol or atomic number of absorbing element 

100 (see Note 1) 

101 edge (string or None): edge for calculation (see Note 2) [None] 

102 cluster_size (float): size of cluster, in Angstroms [8.0] 

103 absorber_site (int): index of site for absorber (see Note 3) [1] 

104 site_index (int or None): index of site for absorber (see Note 4) [None] 

105 extra_titles (list of str or None): extra title lines to include [None] 

106 with_h (bool): whether to include H atoms [False] 

107 version8 (bool): whether to write Feff8l input (see Note 5)[True] 

108 rng_seed (int or None): seed for RNG to get reproducible occupancy selections [None] 

109 Returns 

110 ------- 

111 text of Feff input file 

112 

113 Notes 

114 ----- 

115 1. absorber is the atomic symbol or number of the absorbing element, and 

116 must be an element in the CIF structure. 

117 2. If edge is a string, it must be one of 'K', 'L', 'M', or 'N' edges (note 

118 Feff6 supports only 'K', 'L3', 'L2', and 'L1' edges). If edge is None, 

119 it will be assigned to be 'K' for absorbers with Z < 58 (Ce, with an 

120 edge energy < 40 keV), and 'L3' for absorbers with Z >= 58. 

121 3. for structures with multiple sites for the absorbing atom, the site 

122 can be selected by the order in which they are listed in the sites 

123 list. This depends on the details of the CIF structure, which can be 

124 found with `cif_sites(ciftext)`, starting counting by 1. 

125 4. to explicitly state the index of the site in the sites list, use 

126 site_index (starting at 1!) 

127 5. if version8 is False, outputs will be written for Feff6l 

128 

129 """ 

130 try: 

131 cstruct = read_cif_structure(ciftext) 

132 except ValueError: 

133 return '# could not read CIF file' 

134 

135 global rng 

136 if rng_seed is not None: 

137 rng.seed(rng_seed) 

138 

139 sgroup = SpacegroupAnalyzer(cstruct).get_symmetry_dataset() 

140 space_group = sgroup["international"] 

141 

142 if isinstance(absorber, int): 

143 absorber = atomic_symbol(absorber_z) 

144 absorber_z = atomic_number(absorber) 

145 

146 if edge is None: 

147 edge = 'K' if absorber_z < 58 else 'L3' 

148 

149 edge_energy = xray_edge(absorber, edge).energy 

150 edge_comment = f'{absorber:s} {edge:s} edge, around {edge_energy:.0f} eV' 

151 

152 unique_pot_atoms = [] 

153 for site in cstruct: 

154 for elem in site.species.elements: 

155 if elem.symbol not in unique_pot_atoms: 

156 unique_pot_atoms.append(elem.symbol) 

157 

158 atoms_map = {} 

159 for i, atom in enumerate(unique_pot_atoms): 

160 atoms_map[atom] = i + 1 

161 

162 if absorber not in atoms_map: 

163 atlist = ', '.join(atoms_map.keys()) 

164 raise ValueError(f'atomic symbol {absorber:s} not listed in CIF data: ({atlist})') 

165 

166 site_atoms = {} # map xtal site with list of atoms occupying that site 

167 site_tags = {} 

168 absorber_count = 0 

169 for sindex, site in enumerate(cstruct.sites): 

170 site_species = [e.symbol for e in site.species] 

171 if len(site_species) > 1: 

172 s_els = [s.symbol for s in site.species.keys()] 

173 s_wts = [s for s in site.species.values()] 

174 site_atoms[sindex] = rng.choices(s_els, weights=s_wts, k=1000) 

175 site_tags[sindex] = f'({site.species_string:s})_{1+sindex:d}' 

176 else: 

177 site_atoms[sindex] = [site_species[0]] * 1000 

178 site_tags[sindex] = f'{site.species_string:s}_{1+sindex:d}' 

179 if absorber in site_species: 

180 absorber_count += 1 

181 if absorber_count == absorber_site: 

182 absorber_index = sindex 

183 

184 if site_index is not None: 

185 absorber_index = site_index - 1 

186 

187 # print("Got sites ", len(cstruct.sites), len(site_atoms), len(site_tags)) 

188 

189 center = cstruct[absorber_index].coords 

190 sphere = cstruct.get_neighbors(cstruct[absorber_index], cluster_size) 

191 symbols = [absorber] 

192 coords = [[0, 0, 0]] 

193 tags = [f'{absorber:s}_{1+absorber_index:d}'] 

194 

195 for i, site_dist in enumerate(sphere): 

196 s_index = site_dist[0].index 

197 

198 site_symbol = site_atoms[s_index].pop() 

199 tags.append(site_tags[s_index]) 

200 symbols.append(site_symbol) 

201 coords.append(site_dist[0].coords - center) 

202 cluster = Molecule(symbols, coords) 

203 

204 out_text = ['*** feff input generated by xraylarch cif2feff using pymatgen ***'] 

205 

206 if extra_titles is not None: 

207 for etitle in extra_titles[:]: 

208 if not etitle.startswith('TITLE '): 

209 etitle = 'TITLE ' + etitle 

210 out_text.append(etitle) 

211 

212 out_text.append(f'TITLE Formula: {cstruct.composition.reduced_formula:s}') 

213 out_text.append(f'TITLE SpaceGroup: {space_group:s}') 

214 out_text.append(f'TITLE # sites: {cstruct.num_sites}') 

215 

216 out_text.append('* crystallographics sites: note that these sites may not be unique!') 

217 out_text.append(f'* using absorber at site {1+absorber_index:d} in the list below') 

218 out_text.append(f'* selected as absorber="{absorber:s}", absorber_site={absorber_site:d}') 

219 out_text.append('* index X Y Z species') 

220 for i, site in enumerate(cstruct): 

221 fc = site.frac_coords 

222 species_string = fix_varname(site.species_string.strip()) 

223 marker = ' <- absorber' if (i == absorber_index) else '' 

224 out_text.append(f'* {i+1:3d} {fc[0]:.6f} {fc[1]:.6f} {fc[2]:.6f} {species_string:s} {marker:s}') 

225 

226 out_text.extend(['* ', '', '']) 

227 

228 if version8: 

229 out_text.append(f'EDGE {edge:s}') 

230 out_text.append('S02 1.0') 

231 out_text.append('CONTROL 1 1 1 1 1 1') 

232 out_text.append('PRINT 1 0 0 0 0 3') 

233 out_text.append('EXAFS 20.0') 

234 out_text.append('NLEG 6') 

235 out_text.append(f'RPATH {cluster_size:.2f}') 

236 out_text.append('*SCF 5.0') 

237 

238 else: 

239 edge_index = {'K': 1, 'L1': 2, 'L2': 3, 'L3': 4}[edge] 

240 out_text.append(f'HOLE {edge_index:d} 1.0 * {edge_comment:s} (2nd number is S02)') 

241 out_text.append('CONTROL 1 1 1 0 * phase, paths, feff, chi') 

242 out_text.append('PRINT 1 0 0 0') 

243 out_text.append(f'RMAX {cluster_size:.2f}') 

244 

245 out_text.extend(['', 'EXCHANGE 0', '', 

246 '* POLARIZATION 0 0 0', '', 

247 'POTENTIALS', '* IPOT Z Tag']) 

248 

249 # loop to find atoms actually in cluster, in case some atom 

250 # (maybe fractional occupation) is not included 

251 

252 at_lines = [(0, cluster[0].x, cluster[0].y, cluster[0].z, 0, absorber, tags[0])] 

253 ipot_map = {} 

254 next_ipot = 0 

255 for i, site in enumerate(cluster[1:]): 

256 sym = site.species_string 

257 if sym == 'H' and not with_h: 

258 continue 

259 if sym in ipot_map: 

260 ipot = ipot_map[sym] 

261 else: 

262 next_ipot += 1 

263 ipot_map[sym] = ipot = next_ipot 

264 

265 dist = cluster.get_distance(0, i+1) 

266 at_lines.append((dist, site.x, site.y, site.z, ipot, sym, tags[i+1])) 

267 

268 

269 ipot, z = 0, absorber_z 

270 out_text.append(f' {ipot:4d} {z:4d} {absorber:s}') 

271 for sym, ipot in ipot_map.items(): 

272 z = atomic_number(sym) 

273 out_text.append(f' {ipot:4d} {z:4d} {sym:s}') 

274 

275 out_text.append('') 

276 out_text.append('ATOMS') 

277 out_text.append(f'* x y z ipot tag distance site_info') 

278 

279 acount = 0 

280 for dist, x, y, z, ipot, sym, tag in sorted(at_lines, key=lambda x: x[0]): 

281 acount += 1 

282 if acount > 500: 

283 break 

284 sym = (sym + ' ')[:2] 

285 out_text.append(f' {x: .5f} {y: .5f} {z: .5f} {ipot:4d} {sym:s} {dist:.5f} * {tag:s}') 

286 

287 out_text.append('') 

288 out_text.append('* END') 

289 out_text.append('') 

290 return strict_ascii('\n'.join(out_text))