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

197 statements  

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

1import os 

2from random import Random 

3 

4from xraydb import atomic_symbol, atomic_number, xray_edge 

5from larch.utils.strutils import fix_varname, strict_ascii 

6 

7from .amcsd_utils import (SpacegroupAnalyzer, Molecule, IMolecule, IStructure) 

8 

9rng = Random() 

10 

11def get_atom_map(structure): 

12 """generalization of pymatgen atom map 

13 Returns: 

14 dict of ipots 

15 """ 

16 unique_pot_atoms = [] 

17 all_sites = [] 

18 for site in structure: 

19 for elem in site.species.elements: 

20 if elem.symbol not in unique_pot_atoms: 

21 unique_pot_atoms.append(elem.symbol) 

22 

23 atom_map = {} 

24 for i, atom in enumerate(unique_pot_atoms): 

25 atom_map[atom] = i + 1 

26 return atom_map 

27 

28 

29def read_structure(structure_text, fmt="cif"): 

30 """read structure from text 

31 

32 Arguments 

33 --------- 

34 structure_text (string): text of structure file 

35 fmt (string): format of structure file (cif, poscar, etc) 

36 

37 Returns 

38 ------- 

39 pymatgen Structure object or Molecule object 

40 """ 

41 if Molecule is None: 

42 raise ImportError("pymatgen required. Try 'pip install pymatgen'.") 

43 try: 

44 if fmt.lower() in ('cif', 'poscar', 'contcar', 'chgcar', 'locpot', 'cssr', 'vasprun.xml'): 

45 struct = IStructure.from_str(structure_text, fmt, merge_tol=5.e-4) 

46 else: 

47 struct = IMolecule.from_str(structure_text, fmt) 

48 parse_ok = True 

49 file_found = True 

50 

51 except: 

52 parse_ok = False 

53 file_found = False 

54 if os.path.exists(structure_text): 

55 file_found = True 

56 fmt = os.path.splitext(structure_text)[-1].lower() 

57 try: 

58 if fmt.lower() in ('cif', 'poscar', 'contcar', 'chgcar', 'locpot', 'cssr', 'vasprun.xml'): 

59 struct = IStructure.from_file(structure_text, merge_tol=5.e-4) 

60 else: 

61 struct = IMolecule.from_file(structure_text) 

62 parse_ok = True 

63 except: 

64 parse_ok = False 

65 

66 if not parse_ok: 

67 if not file_found: 

68 raise FileNotFoundError(f'file {structure_text:s} not found') 

69 else: 

70 raise ValueError('invalid text of structure file') 

71 return struct 

72 

73def structure_sites(structure_text, absorber=None, fmt='cif'): 

74 "return list of sites for the structure" 

75 struct = read_structure(structure_text, fmt=fmt) 

76 out = struct.sites 

77 if absorber is not None: 

78 abname = absorber.lower() 

79 out = [] 

80 for site in struct.sites: 

81 species = site.species_string.lower() 

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

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

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

85 if sname.strip() == abname: 

86 out.append(site) 

87 elif species == abname: 

88 out.append(site) 

89 if len(out) == 0: 

90 out = struct.sites[0] 

91 return out 

92 

93def parse_structure(structure_text, fmt='cif', fname="default.filename"): 

94 try: 

95 struct = read_structure(structure_text, fmt=fmt) 

96 except ValueError: 

97 return '# could not read structure file' 

98 

99 return {'formula': struct.composition.reduced_formula, 'sites': struct.sites, 'structure_text': structure_text, 'fmt': fmt, 'fname': fname} 

100 

101 

102def structure2feffinp(structure_text, absorber, edge=None, cluster_size=8.0, 

103 absorber_site=1, site_index=None, extra_titles=None, 

104 with_h=False, version8=True, fmt='cif', rng_seed=None): 

105 """convert structure text to Feff8 or Feff6l input file 

106 

107 Arguments 

108 --------- 

109 structure_text (string): text of CIF file or name of the CIF file. 

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

111 (see Note 1) 

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

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

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

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

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

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

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

119 fmt (string): format of structure file (cif, poscar, etc) [cif] 

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

121 Returns 

122 ------- 

123 text of Feff input file 

124 

125 Notes 

126 ----- 

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

128 must be an element in the CIF structure. 

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

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

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

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

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

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

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

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

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

138 site_index (starting at 1!) 

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

140 

141 """ 

142 try: 

143 struct = read_structure(structure_text, fmt=fmt) 

144 except ValueError: 

145 return '# could not read structure file' 

146 

147 global rng 

148 if rng_seed is not None: 

149 rng.seed(rng_seed) 

150 

151 is_molecule = False 

152 

153 if isinstance(struct, IStructure): 

154 sgroup = SpacegroupAnalyzer(struct).get_symmetry_dataset() 

155 space_group = sgroup["international"] 

156 else: 

157 space_group = 'Molecule' 

158 is_molecule = True 

159 

160 

161 if isinstance(absorber, int): 

162 absorber = atomic_symbol(absorber_z) 

163 absorber_z = atomic_number(absorber) 

164 

165 if edge is None: 

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

167 

168 edge_energy = xray_edge(absorber, edge).energy 

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

170 

171 unique_pot_atoms = [] 

172 for site in struct: 

173 for elem in site.species.elements: 

174 if elem.symbol not in unique_pot_atoms: 

175 unique_pot_atoms.append(elem.symbol) 

176 

177 atoms_map = {} 

178 for i, atom in enumerate(unique_pot_atoms): 

179 atoms_map[atom] = i + 1 

180 

181 if absorber not in atoms_map: 

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

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

184 

185 

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

187 site_tags = {} 

188 absorber_count = 0 

189 for sindex, site in enumerate(struct.sites): 

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

191 if len(site_species) > 1: 

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

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

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

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

196 else: 

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

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

199 if absorber in site_species: 

200 absorber_count += 1 

201 if absorber_count == absorber_site: 

202 absorber_index = sindex 

203 

204 if site_index is not None: 

205 absorber_index = site_index - 1 

206 

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

208 

209 center = struct[absorber_index].coords 

210 sphere = struct.get_neighbors(struct[absorber_index], cluster_size) 

211 symbols = [absorber] 

212 coords = [[0, 0, 0]] 

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

214 

215 for i, site_dist in enumerate(sphere): 

216 s_index = site_dist[0].index 

217 

218 site_symbol = site_atoms[s_index].pop() 

219 tags.append(site_tags[s_index]) 

220 symbols.append(site_symbol) 

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

222 cluster = Molecule(symbols, coords) 

223 

224 out_text = ['*** feff input generated by xraylarch structure2feff using pymatgen ***'] 

225 

226 if extra_titles is not None: 

227 for etitle in extra_titles[:]: 

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

229 etitle = 'TITLE ' + etitle 

230 out_text.append(etitle) 

231 

232 out_text.append(f'TITLE Formula: {struct.composition.reduced_formula:s}') 

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

234 out_text.append(f'TITLE # sites: {struct.num_sites}') 

235 

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

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

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

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

240 

241 for i, site in enumerate(struct): 

242 # The method of obtaining the cooridanates depends on whether the structure is a molecule or not 

243 if is_molecule: 

244 fc = site.coords 

245 else: 

246 fc = site.frac_coords 

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

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

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

250 

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

252 

253 if version8: 

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

255 out_text.append('S02 1.0') 

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

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

258 out_text.append('EXAFS 20.0') 

259 out_text.append('NLEG 6') 

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

261 out_text.append('*SCF 5.0') 

262 

263 else: 

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

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

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

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

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

269 

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

271 '* POLARIZATION 0 0 0', '', 

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

273 

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

275 # (maybe fractional occupation) is not included 

276 

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

278 ipot_map = {} 

279 next_ipot = 0 

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

281 sym = site.species_string 

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

283 continue 

284 if sym in ipot_map: 

285 ipot = ipot_map[sym] 

286 else: 

287 next_ipot += 1 

288 ipot_map[sym] = ipot = next_ipot 

289 

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

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

292 

293 

294 ipot, z = 0, absorber_z 

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

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

297 z = atomic_number(sym) 

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

299 

300 out_text.append('') 

301 out_text.append('ATOMS') 

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

303 

304 acount = 0 

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

306 acount += 1 

307 if acount > 500: 

308 break 

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

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

311 

312 out_text.append('') 

313 out_text.append('* END') 

314 out_text.append('') 

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