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
« 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
5from xraydb import atomic_symbol, atomic_number, xray_edge
6from larch.utils import fix_varname, strict_ascii, gformat
8from .amcsd_utils import PMG_CIF_OPTS, CifParser, Molecule, SpacegroupAnalyzer
10rng = Random()
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)
24 atom_map = {}
25 for i, atom in enumerate(unique_pot_atoms):
26 atom_map[atom] = i + 1
27 return atom_map
30def read_cif_structure(ciftext):
31 """read CIF text, return CIF Structure
33 Arguments
34 ---------
35 ciftext (string): text of CIF file or name of the CIF file.
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
58 try:
59 cstruct = cifstructs.parse_structures()[0]
60 except:
61 raise ValueError('could not get structure from text of CIF file')
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
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
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
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
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
129 """
130 try:
131 cstruct = read_cif_structure(ciftext)
132 except ValueError:
133 return '# could not read CIF file'
135 global rng
136 if rng_seed is not None:
137 rng.seed(rng_seed)
139 sgroup = SpacegroupAnalyzer(cstruct).get_symmetry_dataset()
140 space_group = sgroup["international"]
142 if isinstance(absorber, int):
143 absorber = atomic_symbol(absorber_z)
144 absorber_z = atomic_number(absorber)
146 if edge is None:
147 edge = 'K' if absorber_z < 58 else 'L3'
149 edge_energy = xray_edge(absorber, edge).energy
150 edge_comment = f'{absorber:s} {edge:s} edge, around {edge_energy:.0f} eV'
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)
158 atoms_map = {}
159 for i, atom in enumerate(unique_pot_atoms):
160 atoms_map[atom] = i + 1
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})')
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
184 if site_index is not None:
185 absorber_index = site_index - 1
187 # print("Got sites ", len(cstruct.sites), len(site_atoms), len(site_tags))
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}']
195 for i, site_dist in enumerate(sphere):
196 s_index = site_dist[0].index
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)
204 out_text = ['*** feff input generated by xraylarch cif2feff using pymatgen ***']
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)
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}')
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}')
226 out_text.extend(['* ', '', ''])
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')
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}')
245 out_text.extend(['', 'EXCHANGE 0', '',
246 '* POLARIZATION 0 0 0', '',
247 'POTENTIALS', '* IPOT Z Tag'])
249 # loop to find atoms actually in cluster, in case some atom
250 # (maybe fractional occupation) is not included
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
265 dist = cluster.get_distance(0, i+1)
266 at_lines.append((dist, site.x, site.y, site.z, ipot, sym, tags[i+1]))
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}')
275 out_text.append('')
276 out_text.append('ATOMS')
277 out_text.append(f'* x y z ipot tag distance site_info')
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}')
287 out_text.append('')
288 out_text.append('* END')
289 out_text.append('')
290 return strict_ascii('\n'.join(out_text))