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
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1import os
2from random import Random
4from xraydb import atomic_symbol, atomic_number, xray_edge
5from larch.utils.strutils import fix_varname, strict_ascii
7from .amcsd_utils import (SpacegroupAnalyzer, Molecule, IMolecule, IStructure)
9rng = Random()
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)
23 atom_map = {}
24 for i, atom in enumerate(unique_pot_atoms):
25 atom_map[atom] = i + 1
26 return atom_map
29def read_structure(structure_text, fmt="cif"):
30 """read structure from text
32 Arguments
33 ---------
34 structure_text (string): text of structure file
35 fmt (string): format of structure file (cif, poscar, etc)
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
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
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
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
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'
99 return {'formula': struct.composition.reduced_formula, 'sites': struct.sites, 'structure_text': structure_text, 'fmt': fmt, 'fname': fname}
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
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
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
141 """
142 try:
143 struct = read_structure(structure_text, fmt=fmt)
144 except ValueError:
145 return '# could not read structure file'
147 global rng
148 if rng_seed is not None:
149 rng.seed(rng_seed)
151 is_molecule = False
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
161 if isinstance(absorber, int):
162 absorber = atomic_symbol(absorber_z)
163 absorber_z = atomic_number(absorber)
165 if edge is None:
166 edge = 'K' if absorber_z < 58 else 'L3'
168 edge_energy = xray_edge(absorber, edge).energy
169 edge_comment = f'{absorber:s} {edge:s} edge, around {edge_energy:.0f} eV'
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)
177 atoms_map = {}
178 for i, atom in enumerate(unique_pot_atoms):
179 atoms_map[atom] = i + 1
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})')
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
204 if site_index is not None:
205 absorber_index = site_index - 1
207 # print("Got sites ", len(cstruct.sites), len(site_atoms), len(site_tags))
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}']
215 for i, site_dist in enumerate(sphere):
216 s_index = site_dist[0].index
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)
224 out_text = ['*** feff input generated by xraylarch structure2feff using pymatgen ***']
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)
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}')
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')
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}')
251 out_text.extend(['* ', '', ''])
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')
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}')
270 out_text.extend(['', 'EXCHANGE 0', '',
271 '* POLARIZATION 0 0 0', '',
272 'POTENTIALS', '* IPOT Z Tag'])
274 # loop to find atoms actually in cluster, in case some atom
275 # (maybe fractional occupation) is not included
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
290 dist = cluster.get_distance(0, i+1)
291 at_lines.append((dist, site.x, site.y, site.z, ipot, sym, tags[i+1]))
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}')
300 out_text.append('')
301 out_text.append('ATOMS')
302 out_text.append(f'* x y z ipot tag distance site_info')
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}')
312 out_text.append('')
313 out_text.append('* END')
314 out_text.append('')
315 return strict_ascii('\n'.join(out_text))