Coverage for larch/xrd/struct2xas.py: 0%
615 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
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
4"""
5Struct2XAS: convert CIFs and XYZs files to FDMNES and FEFF inputs
6"""
7# main imports
8import os
9import json
10import time
11import tempfile
12import numpy as np
14# pymatgen
15from pymatgen.core import Structure, Element, Lattice
16from pymatgen.io.xyz import XYZ
17from pymatgen.io.cif import CifParser
18from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
20import larch.utils.logging as logging
21from larch.utils import mkdir, unixpath
22from larch.utils.strutils import fix_filename, unique_name, strict_ascii
23from larch.site_config import user_larchdir
24from larch.io import read_ascii
25from larch.math.convolution1D import lin_gamma, conv
27try:
28 import pandas as pd
29 from pandas.io.formats.style import Styler
31 HAS_PANDAS = True
32except ImportError:
33 HAS_PANDAS = False
35try:
36 import py3Dmol
38 HAS_PY3DMOL = True
39except ImportError:
40 HAS_PY3DMOL = False
43__author__ = ["Beatriz G. Foschiani", "Mauro Rovezzi"]
44__email__ = ["beatrizgfoschiani@gmail.com", "mauro.rovezzi@esrf.fr"]
45__credits__ = ["Jade Chongsathapornpong", "Marius Retegan"]
46__version__ = "2023.3.0"
49# initialize the logger
50logger = logging.getLogger("struct2xas", level="INFO")
53def _get_timestamp() -> str:
54 """return a custom time stamp string: YYY-MM-DD_HHMM"""
55 return "{0:04d}-{1:02d}-{2:02d}_{3:02d}{4:02d}".format(*time.localtime())
58def _pprint(matrix):
59 """pretty print a list of lists (in case Pandas is not available)
60 from: https://stackoverflow.com/questions/13214809/pretty-print-2d-list
61 an alternative could be: https://pypi.org/project/tabulate/
62 """
63 s = [[str(e) for e in row] for row in matrix]
64 lens = [max(map(len, col)) for col in zip(*s)]
65 fmt = "\t".join("{{:{}}}".format(x) for x in lens)
66 table = [fmt.format(*row) for row in s]
67 print("\n".join(table))
70def xyz2struct(molecule):
71 """Convert pymatgen molecule to dummy pymatgen structure"""
73 alat, blat, clat = 1, 1, 1
75 # Set the lattice dimensions in each direction
76 for i in range(len(molecule) - 1):
77 if molecule.cart_coords[i][0] > molecule.cart_coords[i + 1][0]:
78 alat = molecule.cart_coords[i][0]
79 if molecule.cart_coords[i][1] > molecule.cart_coords[i + 1][1]:
80 blat = molecule.cart_coords[i][1]
81 if molecule.cart_coords[i][2] > molecule.cart_coords[i + 1][2]:
82 clat = molecule.cart_coords[i][2]
84 # Set the lattice dimensions in each direction
85 lattice = Lattice.from_parameters(
86 a=alat, b=blat, c=clat, alpha=90, beta=90, gamma=90
87 )
89 # Create a list of species
90 species = [Element(sym) for sym in molecule.species]
92 # Create a list of coordinates
93 coords = molecule.cart_coords
95 # Create the Structure object
96 struct = Structure(lattice, species, coords, coords_are_cartesian=True)
97 return struct
100def structure_folders():
101 """
102 return folders under USER_LARCHDIR for FDMNES, Feff, and structure files,
103 making sure each exists
104 """
105 folders = {}
106 for name in ("feff", "fdmnes", "mp_structs"):
107 folders[name] = unixpath(os.path.join(user_larchdir, name))
108 mkdir(folders[name])
109 return folders
112class Struct2XAS:
113 """Class to convert data from CIF and XYZ files to FDMNES and FEFF inputs"""
115 def __init__(self, file, abs_atom) -> None:
116 """
118 Arguments
119 ---------
120 file : str
121 full path string to cif or xyz file
122 abs_atom : str
123 Absorber element in the structure, e.g.: "Fe", "Ni".
125 Returns
126 -------
127 None
129 ..note::
131 --> IMPORTANT: <--
133 for xyz files:
134 Structures from xyz files are always considered non-symmetric for
135 the lack of lattice information. For creating the object from an
136 XYZ file, a non-symmetric `structure` object from pymatgen is
137 generated (spacegroup: P1) from a `molecule` one via
138 :func:`xyz2struct`. The lattice parameters chosen for this
139 structure are arbitrary and are based on the size of the molecule,
140 as are the fractional coordinates. Therefore, the analysis of this
141 structure is limited to the central atoms and is not valid for
142 atoms at the border of the molecule.
144 for cif files:
145 For creating the object from cif file, a pymatgen structure is
146 generated with symmetry information from cif file.
147 """
149 self.file = file
150 self.abs_atom = abs_atom
151 self.frame = 0
152 self.abs_site = 0
153 self.is_cif = False
154 self.is_xyz = False
155 self.read_structure(file)
156 self.nabs_sites = len(self.get_abs_sites())
157 self.elems = self._get_elems()
158 self.species = self._get_species()
159 self.parent_path = user_larchdir
160 self.folders = structure_folders()
161 logger.info(
162 f"Frames: {self.nframes}, Absorbing sites: {self.nabs_sites}. (Indexes for frames and abs_sites start at 0)"
163 )
165 def set_frame(self, frame=0):
166 """For multiframe xyz file, set the frame index site."""
167 self.frame = frame
168 if self.molecules is not None:
169 self.struct = xyz2struct(self.molecules[self.frame])
170 self.mol = self.molecules[self.frame]
171 else:
172 logger.error("single frame structure")
173 logger.debug(f"frame idx: {frame} ")
175 def get_frame(self):
176 """For multiframe xyz file, get the frame index site."""
177 return self.frame
179 def set_abs_site(self, abs_site=0):
180 """Set the crystallographic index site for the absorbing atom."""
181 try:
182 self.get_abs_sites()[abs_site]
183 except IndexError:
184 logger.error("absorber site index out of range")
185 raise IndexError("absorber site index out of range")
186 logger.debug(f"abs_site idx: {abs_site} ")
187 self.abs_site = abs_site
189 def get_abs_site(self):
190 """Get the crystallographic index site for the absorbing atom."""
191 return self.abs_site
193 def _get_elems(self):
194 """Get elements present in the structure"""
195 elems_list = []
196 for _, elem in enumerate(self.struct):
197 # if self.is_cif:
198 for e in elem.species.elements:
199 if e.name not in elems_list:
200 elems_list.append(e.name)
202 # if self.is_xyz:
203 # if elem.species_string not in elems_list:
204 # elems_list.append(elem.species_string)
206 return elems_list
208 def _get_species(self):
209 """Get elements present in the structure"""
210 species_list = []
211 species_list_occu = []
212 for _, elem in enumerate(self.struct):
213 if not self.full_occupancy:
214 if str(elem.species) not in species_list_occu:
215 species_list_occu.append(str(elem.species))
216 self.atoms_occu = species_list_occu
218 for e in elem.species.elements:
219 if str(e) not in species_list:
220 species_list.append(str(e))
222 return species_list
224 def read_structure(self, file):
225 """Read and initialize the structure/molecule from the input file"""
226 # Split the file name and extension
227 self.file = file
228 if os.path.isfile(self.file):
229 # file_dirname = os.path.dirname(file)
230 file = os.path.basename(self.file)
231 split_file = os.path.splitext(file)
232 self.file_name = split_file[0]
233 self.file_ext = split_file[1]
234 else:
235 self.file_name, self.file_ext = None, None
236 errmsg = f"{self.file} not found"
237 logger.error(errmsg)
238 raise FileNotFoundError(errmsg)
240 if self.file_ext == ".cif":
241 self.is_cif = True
242 self.xyz = None
243 self.molecules = None
244 self.mol = None
245 self.nframes = 1
246 self.cif = CifParser(self.file)
247 self.struct = Structure.from_file(self.file)
248 # self.struct = self.cif.get_structures()[0] #: NOT WORKING!
249 logger.debug("structure created from a CIF file")
250 elif self.file_ext == ".xyz":
251 self.is_xyz = True
252 self.cif = None
253 self.xyz = XYZ.from_file(self.file)
254 self.molecules = self.xyz.all_molecules
255 self.mol = self.molecules[self.frame]
256 self.nframes = len(self.molecules)
257 self.struct = xyz2struct(self.mol)
258 logger.debug("structure created from a XYZ file")
259 elif self.file_ext == ".mpjson":
260 self.is_xyz = False
261 self.is_cif = True
262 self.molecules = None
263 self.mol = None
264 self.nframes = 1
265 self.struct = Structure.from_dict(json.load(open(self.file, "r")))
266 logger.debug("structure created from JSON file")
268 else:
269 errmsg = "only CIF and XYZ files are currently supported"
270 logger.error(errmsg)
271 raise NotImplementedError(errmsg)
273 def _get_atom_sites(self):
274 """get atomic positions from the cif file
276 Returns
277 -------
278 atom_sites : list of list
279 same output as self.get_abs_sites()
280 """
281 if not self.is_cif:
282 errmsg = "not a CIF file!"
283 logger.error(errmsg)
284 raise NotImplementedError(errmsg)
285 cf = self.cif.as_dict()
286 cf = cf[list(cf.keys())[0]]
288 atom_lists = []
289 try:
290 labels = cf["_atom_site_label"]
291 natoms = len(labels)
292 try:
293 type_symbols = cf["_atom_site_type_symbol"]
294 except KeyError:
295 type_symbols = None
296 pass
297 if type_symbols is None:
298 atom_lists.append(labels)
299 else:
300 atom_lists.append(type_symbols)
301 atom_lists.append(cf["_atom_site_fract_x"])
302 atom_lists.append(cf["_atom_site_fract_y"])
303 atom_lists.append(cf["_atom_site_fract_z"])
304 atom_lists.append(cf["_atom_site_occupancy"])
305 except KeyError:
306 errmsg = (
307 "the CIF file does not contain all required `_atom_site_*` information"
308 )
309 logger.error(errmsg)
310 raise KeyError(errmsg)
311 try:
312 atom_lists.append(cf["_atom_site_symmetry_multiplicity"])
313 except KeyError:
314 atom_lists.append(["?"] * natoms)
316 try:
317 atom_lists.append(cf["_atom_site_Wyckoff_symbol"])
318 except KeyError:
319 atom_lists.append(["?"] * natoms)
321 atom_sites = []
322 for ikey, key in enumerate(zip(*atom_lists)):
323 atom_row = [
324 ikey,
325 key[0],
326 np.array(
327 [key[1].split("(")[0], key[2].split("(")[0], key[3].split("(")[0]],
328 dtype=float,
329 ),
330 key[5] + key[6],
331 np.array([None, None, None]),
332 float(key[4]),
333 None,
334 ]
335 atom_sites.append(atom_row)
337 # atom_site_keys = [key for key in cf.keys() if '_atom_site' in key]
338 # atom_lists = [cf[key] for key in atom_site_keys]
339 # atom_sites = []
340 # for ikey, key in enumerate(zip(*atom_lists)):
341 # atom_row = [ikey, key[1], np.array([key[4], key[5], key[6]], dtype=float), key[2] + key[3], np.array([None, None, None]), float(key[8]), None]
342 # atom_sites.append(atom_row)
343 return atom_sites
345 def _get_idx_struct(self, atom_coords):
346 """get the index of the pymatgen Structure corresponding to the given atomic coordinates"""
347 for idx, atom in enumerate(self.struct):
348 if np.allclose(atom.coords, atom_coords, atol=0.001) is True:
349 return idx
350 errmsg = f"atomic coordinates {atom_coords} not found in self.struct"
351 logger.error(errmsg)
352 # raise IndexError(errmsg)
353 return None
355 def get_abs_sites(self):
356 """
357 Get information about the possible absorbing sites present of the
358 structure.
360 ..note:: If the structure has a readable symmetry given by a cif file,
361 the method will return just equivalent sites. If the structure does not
362 have symmetry or the symmetry is not explicit in the files, this method
363 will return all possible sites for absorber atoms.
365 Returns
366 -------
367 abs_sites : list of lists
368 The lists inside the list contain the following respective
369 information:
370 [
371 abs_idx, # index identifier of the absorber site
372 # used by self.set_abs_site().
373 species_string, # The specie for absorber sites.
374 frac_coords, # Fractional coordinate position
375 # If the structure was created using xyz file,
376 # the frac. coords. are arbitrary and the lattice
377 # parameters are based on the molecule size.
378 wyckoff_symbols, # Wyckoff site for absorber sites.
379 # For structures created from xyz
380 # files, Wyckoff sites are always equal to 1a.
381 # (No symmetry)
382 cart_coords, # Cartesian coordinate position
383 occupancy, # Occupancy for absorber sites. For structures
384 # created from xyz files,
385 # occupancy are always equal to 1.
386 structure index # Original index in the pymatgen structure
387 # (private usage)
388 ]
389 """
391 abs_sites = []
392 if self.is_cif:
393 sym_struct = SpacegroupAnalyzer(self.struct).get_symmetrized_structure()
395 # Get multiples sites for absorber atom
396 for idx, sites in enumerate(sym_struct.equivalent_sites):
397 sites = sorted(
398 sites, key=lambda s: tuple(abs(x) for x in s.frac_coords)
399 )
400 site = sites[0]
401 abs_row = [idx, site.species_string]
402 abs_row.append([j for j in np.round(site.frac_coords, 4)])
403 abs_row.append(sym_struct.wyckoff_symbols[idx])
404 abs_row.append(np.array([j for j in np.round(site.coords, 4)]))
405 if self.abs_atom in abs_row[1]:
406 try:
407 ats_occ = abs_row[1].split(",")
408 at_occ = [at for at in ats_occ if self.abs_atom in at][0]
409 occupancy = float(at_occ.split(":")[1])
410 self.full_occupancy = False
411 except Exception:
412 occupancy = 1
413 self.full_occupancy = True
414 abs_row.append(occupancy)
415 abs_row.append(self._get_idx_struct(abs_row[4]))
416 abs_sites.append(abs_row)
418 if self.is_xyz:
419 self.full_occupancy = True
420 k = 0
421 for idx, elem in enumerate(self.struct):
422 if f"{self.abs_atom}" in elem:
423 abs_sites += [
424 (
425 int(f"{k}"),
426 (str(elem.specie)) + f"{k}",
427 elem.coords.round(4),
428 "1a",
429 elem.coords.round(4),
430 1,
431 idx,
432 )
433 ]
434 k += 1
435 if len(abs_sites) == 0:
436 _errmsg = f"---- Absorber {self.abs_atom} not found ----"
437 logger.error(_errmsg)
438 raise AttributeError(_errmsg)
439 return abs_sites
441 def get_abs_sites_info(self):
442 """pretty print for self.get_abs_sites()"""
443 header = [
444 "idx_abs",
445 "specie",
446 "frac_coords",
447 "wyckoff_site",
448 "cart_coords",
449 "occupancy",
450 "idx_in_struct",
451 ]
452 abs_sites = self.get_abs_sites()
453 if HAS_PANDAS:
454 df = pd.DataFrame(
455 abs_sites,
456 columns=header,
457 )
458 df = Styler(df).hide(axis="index")
459 return df
460 else:
461 matrix = [header]
462 matrix.extend(abs_sites)
463 _pprint(matrix)
465 def get_atoms_from_abs(self, radius):
466 """Get atoms in sphere from absorbing atom with certain radius"""
467 abs_sites = self.get_abs_sites()
469 if self.is_cif:
470 nei_list = self.struct.get_sites_in_sphere(
471 abs_sites[self.abs_site][4], float(radius)
472 ) # list os neighbors atoms in sphere
473 sites = []
475 for i in range(len(nei_list)):
476 nei_list[i].cart_coords = (
477 nei_list[i].coords - abs_sites[self.abs_site][4]
478 )
480 for i, _ in enumerate(nei_list):
481 if np.allclose(nei_list[i].cart_coords, [0, 0, 0], atol=0.01) is True:
482 sites.append(
483 [
484 nei_list[i],
485 f"{self.abs_atom}(abs)",
486 0.000,
487 {f"{self.abs_atom}(abs)": 1},
488 ]
489 )
490 nei_list.remove(nei_list[i])
492 for i in range(len(nei_list)):
493 occu_dict = dict(nei_list[i].as_dict()["species"])
494 sites += [
495 [
496 nei_list[i],
497 str(nei_list[i].species.elements[0]),
498 round(np.linalg.norm(nei_list[i].cart_coords - [0, 0, 0]), 5),
499 occu_dict,
500 ]
501 ]
502 if self.is_xyz:
503 nei_list = self.mol.get_sites_in_sphere(
504 abs_sites[self.abs_site][4], float(radius)
505 ) # list os neighbors atoms in sphere
506 sites = []
508 for i in range(len(nei_list)):
509 nei_list[i].cart_coords = (
510 nei_list[i].coords - abs_sites[self.abs_site][4]
511 )
513 for i, _ in enumerate(nei_list):
514 if np.allclose(nei_list[i].cart_coords, [0, 0, 0], atol=0.01) is True:
515 sites.append([nei_list[i], f"{self.abs_atom}(abs)", 0.000])
516 nei_list.remove(nei_list[i])
518 sites += [
519 [
520 nei_list[i],
521 nei_list[i].species_string,
522 round(np.linalg.norm(nei_list[i].cart_coords - [0, 0, 0]), 5),
523 ]
524 for i in range(len(nei_list))
525 ]
527 sites = sorted(sites, key=lambda x: x[2])
528 return sites
530 def get_coord_envs(self):
531 """
532 For structures from cif files, this method will try to find the
533 coordination environment type and return the elements and the
534 coordination env. symbol from the first using the classes from pymatgen
535 as LocalGeometryFinder(), BVAnalyzer(), MultiWeightsChemenvStrategy()
536 and LightStructureEnvironments() .
538 > coordination env. symbol.
539 - `S:4` - Square Plane
540 - `T:4` - Tetrahedral
541 - `T:5` - Trigonal bipyramid
542 - `S:5` - Square pyramidal
543 - `O:6` - Octahedral
544 - `T:6` - Trigonal prism
545 > ce_fraction:
546 probability for given coordination env. (between 0 and
547 1)
549 > CSM:
550 a measure of the degree of symmetry in the coordination
551 environment. It is based on the idea that symmetric
552 environments are more stable than asymmetric ones, and
553 is calculated using a formula that takes into account
554 the distances and angles between the coordinating
555 atoms. The CSM can be understood as a distance to a
556 shape and can take values between 0.0 (if a given
557 environment is perfect) and 100.0 (if a given
558 environment is very distorted). The environment of the
559 atom is then the model polyhedron for which the
560 similarity is the highest, that is, for which the CSM
561 is the lowest.
563 > permutation:
564 possible permutation of atoms surrounding the central
565 atom. This is a list that indicates the order in which
566 the neighboring atoms are arranged around the central
567 atom in the coordination environment. The numbering
568 starts from 0, and the list indicates the indices of
569 the neighboring atoms in this order. For example, in
570 the second entry of the list above, the permutation [0,
571 2, 3, 1, 4] means that the first neighboring atom is in
572 position 0, the second is in position 2, the third is
573 in position 3, the fourth is in position 1, and the
574 fifth is in position 4. The permutation is used to
575 calculate the csm value.
577 > site:
578 element in the coordination environment and its
579 coordinates (cartesian and fractional).
581 > site_index:
582 structure index for the coordinated atom.
585 For structures from the xyz file the methods will try to return the
586 elements (but not the coord. env. symbol) for the first coordination
587 env. shell using the the class CrystalNN from pymatgen, which gives
588 better results than BrunnerNN_real and CutOffDictNN (as previously
589 tested)
591 List of lists:
593 [0]: Info about which site is being analyzed.
595 [1]: Coord. env as dictionary.
597 [2]: Info about coord. env.
598 > site:
599 element in the coordination environment and its coordinates
600 (cartesian and fractional).
602 > image:
603 image is defined as displacement from original site in
604 structure to a given site. i.e. if structure has a site at
605 (-0.1, 1.0, 0.3), then (0.9, 0, 2.3) -> jimage = (1, -1,
606 2). Note that this method takes O(number of sites) due to
607 searching an original site.
609 > weight:
610 quantifies the significance or contribution of each
611 coordinated site to the central site's coordination.
613 > site_index:
614 structure index for the coerdinated atom.
617 """
618 abs_sites = self.get_abs_sites()
619 idx_abs_site = abs_sites[self.abs_site][-1]
621 if self.is_cif:
622 from pymatgen.analysis.bond_valence import BVAnalyzer
623 from pymatgen.analysis.chemenv.coordination_environments.coordination_geometry_finder import (
624 LocalGeometryFinder,
625 )
626 from pymatgen.analysis.chemenv.coordination_environments.structure_environments import (
627 LightStructureEnvironments,
628 )
629 from pymatgen.analysis.chemenv.coordination_environments.chemenv_strategies import (
630 # SimplestChemenvStrategy,
631 MultiWeightsChemenvStrategy,
632 # WeightedNbSetChemenvStrategy,
633 )
635 lgf = LocalGeometryFinder()
636 lgf.setup_structure(self.struct)
638 bva = BVAnalyzer() # Bond Valence Analyzer
639 try:
640 valences = bva.get_valences(structure=self.struct)
641 except ValueError:
642 valences = "undefined"
644 coord_env_list = []
645 se = lgf.compute_structure_environments(
646 max_cn=6,
647 valences=valences,
648 only_indices=[idx_abs_site],
649 only_symbols=["S:4", "T:4", "T:5", "S:5", "O:6", "T:6"],
650 )
652 dist_1st_shell = se.voronoi.neighbors_distances[idx_abs_site][0]["max"]
653 logger.debug(dist_1st_shell)
654 # atom_coord = lgf.compute_coordination_environments(self.struct)
655 strategy = MultiWeightsChemenvStrategy.stats_article_weights_parameters()
656 # strategy = SimplestChemenvStrategy(distance_cutoff=1.1, angle_cutoff=0.3)
657 lse = LightStructureEnvironments.from_structure_environments(
658 strategy=strategy, structure_environments=se
659 )
660 coord_env_ce = lse.coordination_environments[idx_abs_site]
661 ngbs_sites = lse._all_nbs_sites
662 coord_env_list.append(
663 [
664 f"Coord. Env. for Site {abs_sites[self.abs_site][0]}",
665 coord_env_ce,
666 ngbs_sites,
667 ]
668 )
670 if self.is_xyz:
671 from pymatgen.analysis.local_env import CrystalNN
673 obj = CrystalNN()
674 coord_env_list = []
675 coord_env = obj.get_nn_info(self.struct, idx_abs_site)
676 for site in coord_env:
677 site["site"].cart_coords = self.struct[site["site_index"]].coords
678 coord_dict = obj.get_cn_dict(self.struct, idx_abs_site)
679 coord_env_list.append(
680 [
681 f"Coord. Env. for Site {abs_sites[self.abs_site][0]}",
682 {"ce_symbol": f"Elements Dict = {coord_dict}"},
683 coord_env,
684 ]
685 )
687 return coord_env_list
689 def get_coord_envs_info(self):
690 """
691 Class with summarized and more readable information from get_coord_envs() method
692 """
694 coord_env = self.get_coord_envs()[0]
695 abs_site_coord = self.get_abs_sites()[self.abs_site][4]
697 elems_dist = []
698 for site in coord_env[2]:
699 if self.is_cif:
700 coord_sym = [
701 coord_env[1][i]["ce_symbol"] for i in range(len(coord_env[1]))
702 ]
703 elems_dist.append(
704 (
705 site["site"].species,
706 round(np.linalg.norm(site["site"].coords - abs_site_coord), 5),
707 )
708 )
709 if self.is_xyz:
710 coord_sym = coord_env[1]["ce_symbol"]
711 elems_dist.append(
712 (
713 site["site"].species,
714 round(
715 np.linalg.norm(site["site"].cart_coords - abs_site_coord), 5
716 ),
717 )
718 )
719 elems_dist = sorted(elems_dist, key=lambda x: x[1])
720 print(
721 f"Coord. Env. from absorber atom: {self.abs_atom} at site {self.get_abs_site()}"
722 )
723 print(coord_sym)
724 header = ["Element", "Distance"]
725 if HAS_PANDAS:
726 df = pd.DataFrame(data=elems_dist, columns=header)
727 return df
728 else:
729 matrix = [header]
730 matrix.extend(elems_dist)
731 _pprint(matrix)
733 def make_cluster(self, radius):
734 """Create a cluster with absorber atom site at the center.
736 Arguments
737 ---------
738 radius :float
739 cluster radius [Angstrom]
741 Returns
742 -------
743 atoms : list
744 species and coords for the new cluster structure
745 """
747 selected_site = self.get_abs_sites()[self.abs_site]
748 cluster = self.mol.get_sites_in_sphere(selected_site[-3], radius)
750 # showing and storing cartesian coords and species
751 atoms = []
753 # abs_atom at the cluster center
754 for i in range((len(cluster))):
755 try:
756 species = round(Element((cluster[i].specie).element).Z)
757 except AttributeError:
758 species = round(Element(cluster[i].specie).Z)
760 # getting species, after atomic number() and rounding
761 coords = (
762 cluster[i].coords - selected_site[2]
763 ) # cartesial coords and ""frac_coords"" are the same for molecule structure (a = b = c = 1)
764 coords = np.around(coords, 5)
765 dist = round(np.linalg.norm(coords - [0, 0, 0]), 5)
766 atoms.append((species, coords, dist))
767 atoms = sorted(atoms, key=lambda x: x[2])
768 return atoms
770 def make_input_fdmnes(
771 self, radius=7, parent_path=None, template=None, green=True, **kwargs
772 ):
773 """
774 Create a fdmnes input from a template.
776 Arguments
777 ---------
778 radius : float, [7]
779 radius for fdmnes calculation in Angstrom
780 parent_path : str, [None]
781 path to the parent directory where the input files are stored
782 if None it will create a temporary directory
783 template : str, [None]
784 full path to the template file
785 green : boolean [True]
786 True: use `Green` (muffin-tin potentials, faster)
787 False: use finite-difference method (slower)
789 ..note:: SCF is always used
790 ..note:: for further information about FDMNES keywords, search for "FDMNES users guide"
792 Returns
793 -------
794 None -> writes FDMNES input to disk
795 directory structure: {parent_path}/fdmnes/{self.file_name}/{self.abs_atom}/frame{self.frame}/site{self.abs_site}/
797 """
798 replacements = {}
799 replacements.update(**kwargs)
800 replacements["version"] = __version__
802 if template is None:
803 template = os.path.join(
804 os.path.dirname(os.path.realpath(__file__)), "templates", "fdmnes.tmpl"
805 )
807 if parent_path is None:
808 parent_path = self.folders["fdmnes"]
810 self.outdir = os.path.join(
811 parent_path,
812 self.file_name,
813 self.abs_atom,
814 f"frame{self.frame}",
815 f"site{self.abs_site}",
816 )
818 method = "green" if green else ""
819 absorber = ""
820 crystal = ""
821 occupancy = ""
822 comment = f"input structure: {self.file_name}{self.file_ext}\ncreation date: {_get_timestamp()}"
824 if self.is_cif:
825 try:
826 selected_site = self.get_abs_sites()[self.abs_site]
827 except IndexError:
828 logger.error("IndexError: check if abs_atom is correct")
830 if not selected_site[-2] == 1:
831 logger.warning("the selected site does not have full occupancy!")
833 # SpacegroupAnalyzer to get symmetric structure
834 analyzer = SpacegroupAnalyzer(self.struct)
835 structure = analyzer.get_refined_structure()
837 symmetry_data = analyzer.get_symmetry_dataset()
838 group_number = symmetry_data["number"]
839 group_choice = symmetry_data["choice"]
841 # FDMNES doesn't recognize 2 as a space group.
842 if group_number == 2:
843 group_number = "P-1"
845 crystal = f"{crystal}"
846 replacements["crystal"] = "crystal"
848 group = f"{group_number}"
849 if group_choice:
850 group += f":{group_choice}"
851 replacements["group"] = f"spgroup\n{group}"
853 unique_sites = []
854 for sites in analyzer.get_symmetrized_structure().equivalent_sites:
855 sites = sorted(
856 sites, key=lambda s: tuple(abs(x) for x in s.frac_coords)
857 )
858 unique_sites.append((sites[0], len(sites)))
859 sites = str()
860 if self.full_occupancy:
861 for i, (site, _) in enumerate(unique_sites):
862 try:
863 e = site.specie
864 except AttributeError:
865 e = Element(site.species_string.split(":")[0])
866 sites += "\n" + (
867 f"{e.Z:>2d} {site.a:12.8f} {site.b:12.8f} {site.c:12.8f}"
868 f" {site.species_string:>4s}"
869 )
870 else:
871 occupancy = "occupancy"
872 for site, _ in unique_sites:
873 for i, e in enumerate(site.species.elements):
874 occu = site.as_dict()["species"][i]["occu"]
875 sites += "\n" + (
876 f"{e.Z:>2d} {site.a:12.8f} {site.b:12.8f} {site.c:12.8f} {occu}"
877 f" {str(e):>4s}"
878 )
879 lat = structure.lattice
880 replacements["lattice"] = (
881 f"{lat.a:<12.8f} {lat.b:12.8f} {lat.c:12.8f} "
882 f"{lat.alpha:12.8f} {lat.beta:12.8f} {lat.gamma:12.8f}"
883 )
885 absorber = f"{absorber}"
886 for i in range(len(unique_sites)):
887 if (
888 np.allclose(unique_sites[i][0].coords, selected_site[4], atol=0.01)
889 is True
890 ):
891 replacements["absorber"] = f"absorber\n{i+1}"
893 # absorber = f"{absorber}"
894 # replacements["absorber"] = f"Z_absorber\n{round(Element(elem).Z)}"
896 if self.is_xyz:
897 replacements["crystal"] = "molecule"
899 atoms = self.make_cluster(radius=radius)
900 sites = str()
901 for i in range(len(atoms)):
902 e = atoms[i][0]
903 c = atoms[i][1]
904 sites += "\n" + (
905 f"{e:>2d} {c[0]:12.8f} {c[1]:12.8f} {c[2]:12.8f}"
906 f" {Element.from_Z(e).name}"
907 )
909 absorber = f"{absorber}"
910 for i in range(len(atoms)):
911 if np.allclose(atoms[i][1], [0, 0, 0], atol=0.01) is True:
912 replacements["absorber"] = f"absorber\n{i+1}"
914 replacements["group"] = ""
916 lat = self.struct.lattice
917 replacements["lattice"] = (
918 f"{lat.a:<12.8f} {lat.b:12.8f} {lat.c:12.8f} "
919 f"{lat.alpha:12.8f} {lat.beta:12.8f} {lat.gamma:12.8f}"
920 )
922 replacements["sites"] = sites
923 replacements["radius"] = radius
924 replacements["method"] = method
925 replacements["comment"] = comment
926 replacements["occupancy"] = occupancy
928 try:
929 os.makedirs(self.outdir, mode=0o755)
930 except FileExistsError:
931 pass
933 # Write the input file.
934 fnout = os.path.join(self.outdir, "job_inp.txt")
935 with open(fnout, "w") as fp, open(template) as tp:
936 inp = tp.read().format(**replacements)
937 fp.write(inp)
939 # Write the fdmfile.txt.
940 with open(os.path.join(self.outdir, "fdmfile.txt"), "w") as fp:
941 fp.write("1\njob_inp.txt")
943 logger.info(f"written FDMNES input -> {fnout}")
945 def make_input_feff(
946 self,
947 radius=7,
948 parent_path=None,
949 template=None,
950 feff_comment="*",
951 edge="K",
952 sig2=None,
953 debye=None,
954 **kwargs,
955 ):
956 """
957 Create a FEFF input from a template.
959 Arguments
960 ---------
961 radius : float, [7]
962 radius for feff calculation [Angstrom].
963 parent_path : str, [None]
964 path to the parent directory where the input files are stored
965 if None it will create a temporary directory
966 template : str, [None]
967 full path to the template file
968 feff_coment : str, ["*"]
969 comment character used in the input file
970 sig2 : float or None, [None]
971 SIG2 keywork, if None it will be commented
972 debye : list of two floats or None, [None]
973 DEBYE keyword, if None it will be commented, otherwise:
974 debye=[temperature, debye_temperature]
975 temperatue : float
976 temperature at which the Debye-Waller factors are calculated [Kelvin].
977 debye_temperature : float
978 Debye Temperature of the material [Kelvin].
980 ..note:: refer to [FEFF documentation](https://feff.phys.washington.edu/feffproject-feff-documentation.html)
982 Returns
983 -------
984 None -> writes FEFF input to disk
985 directory structure: {parent_path}/feff/{self.file_name}/{self.abs_atom}/frame{self.frame}/site{self.abs_site}/
986 """
987 replacements = {}
988 replacements.update(**kwargs)
989 replacements["version"] = __version__
991 if parent_path is None:
992 parent_path = self.folders["feff"]
993 self.outdir = os.path.join(
994 parent_path,
995 self.file_name,
996 self.abs_atom,
997 f"frame{self.frame}",
998 f"site{self.abs_site}",
999 )
1001 if template is None:
1002 template = os.path.join(
1003 os.path.dirname(os.path.realpath(__file__)),
1004 "templates",
1005 "feff_exafs.tmpl",
1006 )
1008 if sig2 is None:
1009 use_sig2 = "*"
1010 sig2 = 0.005
1011 else:
1012 use_sig2 = ""
1014 if debye is None:
1015 use_debye = "*"
1016 temperature, debye_temperature = 0, 0
1017 else:
1018 use_debye = ""
1019 temperature, debye_temperature = debye[0], debye[1]
1021 feff_comment = f"{feff_comment}"
1022 edge = f"{edge}"
1023 radius = f"{radius}"
1024 use_sig2 = f"{use_sig2}"
1025 sig2 = f"{sig2}"
1026 use_debye = f"{use_debye}"
1027 temperature = f"{temperature}"
1028 debye_temperature = f"{debye_temperature}"
1030 if self.is_cif:
1031 sites = self.get_atoms_from_abs(radius)
1032 ipot_list = [(0, Element(f"{self.abs_atom}").Z, sites[0][1])]
1033 ipot = {f"{self.abs_atom}(abs)": 0}
1034 elems = self.species
1036 for i, elem in enumerate(elems):
1037 ipot[elem] = i + 1
1038 for i in range(1, len(sites)):
1039 for j, _ in enumerate(sites[i][0].species.elements):
1040 if str(sites[i][0].species.elements[j]) in elems:
1041 ipot_list.append(
1042 (
1043 ipot[str(sites[i][0].species.elements[j])],
1044 sites[i][0].species.elements[j].Z,
1045 str(sites[i][0].species.elements[j]),
1046 )
1047 )
1048 pot = list(dict.fromkeys(ipot_list))
1049 potentials = str("* ipot Z tag [lmax1 lmax2 xnatph sphinph]")
1050 for i, _ in enumerate(pot):
1051 potentials += "\n" + (f"{pot[i][0]:>5} {pot[i][1]:>3} {pot[i][2]:>5}")
1053 atoms_list = [
1054 (sites[0][0].cart_coords, 0, sites[0][1], sites[0][2], sites[0][3])
1055 ]
1056 for i in range(1, len(sites)):
1057 atoms_list.append(
1058 (
1059 sites[i][0].cart_coords,
1060 ipot[str(sites[i][0].species.elements[0])],
1061 sites[i][1],
1062 sites[i][2],
1063 sites[i][3],
1064 ) # cart, ipot, tag, dist
1065 )
1067 if self.is_xyz:
1068 sites = self.get_atoms_from_abs(radius=radius)
1069 ipot_list = [(0, Element(f"{self.abs_atom}").Z, sites[0][1])]
1070 elems = []
1071 ipot = {}
1072 for i in range(1, len(sites)):
1073 if sites[i][0].species_string not in elems:
1074 elems.append((sites[i][0].species_string))
1075 for i, elem in enumerate(elems):
1076 ipot[elem] = i + 1
1077 for i in range(1, len(sites)):
1078 if sites[i][0].species_string in elems:
1079 ipot_list.append(
1080 (
1081 ipot[sites[i][0].species_string],
1082 sites[i][0].specie.Z,
1083 sites[i][0].species_string,
1084 )
1085 )
1086 pot = list(dict.fromkeys(ipot_list))
1087 potentials = str("* ipot Z tag [lmax1 lmax2 xnatph sphinph]")
1088 for i in range(len(pot)):
1089 potentials += "\n" + (f"{pot[i][0]:>5} {pot[i][1]:>3} {pot[i][2]:>5}")
1090 atoms_list = []
1091 for i in range(len(sites)):
1092 atoms_list.append(
1093 (
1094 sites[i][0].cart_coords,
1095 ipot_list[i][0],
1096 sites[i][1], # tag
1097 sites[i][2],
1098 )
1099 )
1100 atoms = str(
1101 "* x y z ipot tag distance occupancy"
1102 )
1103 at = atoms_list
1104 for i in range(len(at)):
1105 if self.full_occupancy:
1106 atoms += "\n" + (
1107 f"{at[i][0][0]:10.6f} {at[i][0][1]:10.6f} {at[i][0][2]:10.6f} { int(at[i][1])} {at[i][2]:>5} {at[i][3]:10.5f} *1 "
1108 )
1109 else:
1110 choice = np.random.choice(
1111 list(at[i][4].keys()), p=list(at[i][4].values())
1112 )
1113 atoms += "\n" + (
1114 f"{at[i][0][0]:10.6f} {at[i][0][1]:10.6f} {at[i][0][2]:10.6f} {ipot[choice]} {choice:>5} {at[i][3]:10.5f} *{at[i][4]}"
1115 )
1117 title = f"TITLE {self.file_name}{self.file_ext}\nTITLE {_get_timestamp()}\nTITLE site {self.abs_site}"
1119 replacements["feff_comment"] = feff_comment
1120 replacements["edge"] = edge
1121 replacements["radius"] = radius
1122 replacements["use_sig2"] = use_sig2
1123 replacements["sig2"] = sig2
1124 replacements["use_debye"] = use_debye
1125 replacements["temperature"] = temperature
1126 replacements["debye_temperature"] = debye_temperature
1127 replacements["potentials"] = potentials
1128 replacements["atoms"] = atoms
1129 replacements["title"] = title
1130 # replacements[""] =
1132 try:
1133 os.makedirs(self.outdir, mode=0o755)
1134 except FileExistsError:
1135 pass
1137 # Write the input file.
1138 fnout = os.path.join(self.outdir, "feff.inp")
1139 with open(fnout, "w") as fp, open(template) as tp:
1140 inp = tp.read().format(**replacements)
1141 fp.write(inp)
1143 logger.info(f"written FEFF input -> {fnout}")
1145 def _get_xyz_and_elements(self, radius):
1146 """
1147 Get information about cartesian coords and elements surrounding the central atom given
1148 a radius.
1150 Args:
1151 > radius(float):
1152 radius from the central atom [Angstrom].
1154 return list of elements with coords and list of elements, both lists of strings.
1156 """
1157 sites = self.get_atoms_from_abs(radius)
1158 coords = []
1159 elements = []
1160 for _, site in enumerate(sites):
1161 try:
1162 coords.append(
1163 (str((site[0].species).elements[0].name), site[0].cart_coords)
1164 )
1165 except AttributeError:
1166 coords.append((str(site[0].specie), site[0].cart_coords))
1168 output_str = str(len(coords)) + "\n\n"
1169 for element, coords in coords:
1170 coords_str = " ".join([f"{c:.6f}" for c in coords])
1171 output_str += f"{element} {coords_str}\n"
1172 if element not in elements:
1173 elements.append(element)
1174 elements = sorted(elements)
1175 return output_str, elements
1177 def _round_up(self, x):
1178 rounded_x = np.ceil(x * 100) / 100
1179 return rounded_x
1181 def visualize(self, radius=2.5, unitcell=False):
1182 """
1183 Display a 3D visualization for material local structure.
1185 Args:
1186 > radius (float):
1187 radius visualization from the central atom.
1189 > unitcell (boolean):
1190 if True, allows the visualization of the structure unit cell.
1192 return 3D structure visualization from py3Dmol.
1193 """
1194 if HAS_PY3DMOL is False:
1195 logger.error("py3Dmol not installed! -> run `pip install py3Dmol`")
1196 return
1198 radius = self._round_up(radius)
1200 xyz, elems = self._get_xyz_and_elements(radius)
1202 a = self.struct.lattice.a
1203 b = self.struct.lattice.b
1204 c = self.struct.lattice.c
1205 alpha = self.struct.lattice.alpha
1206 beta = self.struct.lattice.beta
1207 gamma = self.struct.lattice.gamma
1208 xyzview = py3Dmol.view(
1209 width=600, height=600
1210 ) # http://3dmol.org/doc/GLViewer.html#setStyle
1211 xyzview.addModel(xyz, "xyz")
1213 if unitcell is True:
1214 m = xyzview.getModel()
1215 m.setCrystData(a, b, c, alpha, beta, gamma)
1216 xyzview.addUnitCell()
1218 colors = [
1219 "red",
1220 "green",
1221 "blue",
1222 "orange",
1223 "yellow",
1224 "white",
1225 "purple",
1226 "pink",
1227 "brown",
1228 "black",
1229 "gray",
1230 "cyan",
1231 "magenta",
1232 "olive",
1233 "navy",
1234 "teal",
1235 "maroon",
1236 "turquoise",
1237 "indigo",
1238 "salmon",
1239 ]
1240 color_elems = {}
1241 for idx, elem in enumerate(self.elems):
1242 color_elems[f"{elem}"] = colors[idx]
1244 for idx, elem in enumerate(elems):
1245 color = color_elems[f"{elem}"]
1246 xyzview.setStyle(
1247 {"elem": f"{elem}"},
1248 {
1249 "stick": {
1250 "radius": 0.1,
1251 "opacity": 1,
1252 "hidden": False,
1253 "color": f"{color}",
1254 },
1255 "sphere": {"color": f"{color}", "radius": 0.4, "opacity": 1},
1256 },
1257 )
1258 xyzview.addLabel(
1259 "Abs",
1260 {
1261 "fontColor": "black",
1262 "fontSize": 14,
1263 "backgroundColor": "white",
1264 "backgroundOpacity": 0.8,
1265 "showBackground": True,
1266 },
1267 {"index": 0},
1268 )
1270 xyzview.zoomTo()
1271 xyzview.show()
1273 logger.info(color_elems)
1274 if not self.full_occupancy:
1275 logger.warning("3D displayed image does not consider partial occupancy")
1276 logger.info("check atoms occupancy here:", self.atoms_occu)
1277 # color_elems = {}
1278 # for idx, elem in enumerate(self.species_occu):
1279 # color_elems[f"{list(elem.values())[0]}"] = colors[idx]
1280 # print("Label:\n", color_elems)
1283def get_fdmnes_info(file, labels=("energy", "mu")):
1284 """Get info from the fdmnes output file such as edge energy, atomic number Z,
1285 and fermi level energy, and returns a group with the storage information
1287 Parameters:
1289 file (str): path to the fdmnes output file.
1290 Obs: The INPUT file must have the "Header" keyword to use this function in the OUTPUT file
1292 """
1293 group = read_ascii(file, labels=labels)
1295 with open(group.path) as f:
1296 line = f.readlines()[3]
1297 header = line.split()
1298 (
1299 e_edge,
1300 Z,
1301 e_fermi,
1302 ) = (float(header[0]), float(header[1]), float(header[6]))
1303 print(
1304 f"Calculated Fermi level: {e_fermi}\nAtomic_number: {Z}\nEnergy_edge: {e_edge}"
1305 )
1307 group.e_edge = e_edge
1308 group.Z = Z
1309 group.e_fermi = e_fermi
1311 return group
1314def convolve_data(
1315 energy, mu, group, fwhm=1, linbroad=[1.5, 0, 50], kernel="gaussian", efermi=None
1316):
1317 """
1318 Function for manual convolution using Convolution1D from larch and returning a group
1320 Generic discrete convolution
1322 Description
1323 -----------
1325 This is a manual (not optimized!) implementation of discrete 1D
1326 convolution intended for spectroscopy analysis. The difference with
1327 commonly used methods is the possibility to adapt the convolution
1328 kernel for each convolution point, e.g. change the FWHM of the
1329 Gaussian kernel as a function of the energy scale.
1331 Resources
1332 ---------
1334 .. [WPconv] <http://en.wikipedia.org/wiki/Convolution#Discrete_convolution>
1335 .. [Fisher] <http://homepages.inf.ed.ac.uk/rbf/HIPR2/convolve.htm>
1336 .. [GP1202] <http://glowingpython.blogspot.fr/2012/02/convolution-with-numpy.html>
1338 """
1340 gamma_e = lin_gamma(energy, fwhm=fwhm, linbroad=linbroad)
1341 mu_conv = conv(energy, mu, kernel=kernel, fwhm_e=gamma_e, efermi=efermi)
1342 group.conv = mu_conv
1343 return group
1346def save_cif_from_mp(api_key, material_id, parent_path=None):
1347 """Collect a CIF file from the Materials Project Database, given the material id
1349 Parameters
1350 ----------
1351 api_key : str
1352 api-key from Materials Project
1353 material id : str
1354 material id (format mp-xxxx) from Materials Project
1355 parent_path : str
1356 path where to store the CIF files
1357 if None, a temporary one is created
1359 Returns
1360 -------
1361 [str, str] : parent_path, CIF file name
1363 """
1364 if parent_path is None:
1365 parent_path = structure_folders()["mp_structs"]
1367 from pymatgen.ext.matproj import _MPResterLegacy
1369 cif = _MPResterLegacy(api_key).get_data(material_id, prop="cif")
1370 pf = _MPResterLegacy(api_key).get_data(material_id, prop="pretty_formula")[0][
1371 "pretty_formula"
1372 ]
1373 outfn = f"{pf}_{material_id}.cif"
1374 with open(file=os.path.join(parent_path, outfn), mode="w") as f:
1375 f.write(cif[0]["cif"])
1376 logger.info(f"{material_id} -> {outfn}")
1377 return [parent_path, outfn]
1380def save_mp_structure(api_key, material_id, parent_path=None):
1381 """Save structure from Materials Project Database as json, given the material id
1383 Parameters
1384 ----------
1385 api_key : str
1386 api-key from Materials Project
1387 material id : str
1388 material id (format mp-xxxx) from Materials Project
1389 parent_path : str
1390 path where to store the Structure files
1391 if None, user_larchdir + 'mp_structs' is used
1393 Returns
1394 -------
1395 name of structure file, which will have an 'mpjson' extension
1397 Notes
1398 ------
1399 The structure is saved as json that can be loaded with
1400 from pymatgen.core import Structure
1401 import json
1402 struct = Structure.from_dict(json.load(open(filename, 'r')))
1404 """
1406 if parent_path is None:
1407 parent_path = structure_folders()["mp_structs"]
1409 try:
1410 from mp_api.client import MPRester
1411 except ImportError:
1412 print("need to install mp_api: pip install mp_api")
1414 mpr = MPRester(api_key)
1415 results = mpr.summary.search(
1416 material_ids=[material_id], fields=["structure", "formula_pretty"]
1417 )
1418 formula = results[0].formula_pretty
1419 structure = results[0].structure
1421 outfile = os.path.join(parent_path, f"{formula}_{material_id}.mpjson")
1422 with open(outfile, "w") as fh:
1423 fh.write(structure.to_json())
1424 logger.info(f"saved {material_id} to {outfile}")
1426 return outfile