# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
from __future__ import division, unicode_literals
"""
This module defines classes for point defects.
"""
__author__ = "Bharat Medasani, Nils E. R. Zimmermann"
__copyright__ = "Copyright 2011, The Materials Project"
__version__ = "1.0"
__maintainer__ = "Bharat Medasani, Nils E. R. Zimmermann"
__email__ = "mbkumar@gmail.com, n.zimmermann@tuhh.de"
__status__ = "Production"
__date__ = "Nov 28, 2016"
import abc
import itertools
import numpy as np
import logging
import six
import pandas as pd
from collections import defaultdict
from scipy.spatial import Voronoi
from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import linkage, fcluster
from pymatgen.core.periodic_table import Specie, Element, get_el_sp
from pymatgen.core.sites import PeriodicSite
from pymatgen.core.structure import Structure
from pymatgen.io.vasp.outputs import Chgcar
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.io.zeopp import get_voronoi_nodes, get_void_volume_surfarea, \
get_high_accuracy_voronoi_nodes
from pymatgen.command_line.gulp_caller import get_energy_buckingham, \
get_energy_relax_structure_buckingham
from pymatgen.analysis.local_env import LocalStructOrderParams, \
MinimumDistanceNN, VoronoiNN, ValenceIonicRadiusEvaluator , \
cn_opt_params
from pymatgen.analysis.structure_analyzer import RelaxationAnalyzer
from pymatgen.analysis.structure_matcher import StructureMatcher
from pymatgen.analysis.phase_diagram import get_facets
from pymatgen.analysis.bond_valence import BVAnalyzer
from pymatgen.util.coord import pbc_diff
from pymatgen.vis.structure_vtk import StructureVis
from monty.dev import requires
try:
from skimage.feature import peak_local_max
peak_local_max_found = True
except ImportError:
peak_local_max_found = False
from six.moves import filter
from six.moves import map
from copy import deepcopy
logger = logging.getLogger(__name__)
motif_cn_op = {}
for cn, di in cn_opt_params.items():
for mot, li in di.items():
motif_cn_op[mot] = {'cn': int(cn), 'optype': li[0]}
motif_cn_op[mot]['params'] = deepcopy(li[1]) if len(li) > 1 else None
[docs]class Defect(six.with_metaclass(abc.ABCMeta, object)):
"""
Abstract class for point defects
"""
[docs] @abc.abstractmethod
def enumerate_defectsites(self):
"""
Enumerates all the symmetrically distinct defects.
"""
raise NotImplementedError()
@property
def structure(self):
"""
Returns the structure without any defects
Useful for Mott-Littleton calculations.
"""
return self._structure
@property
def struct_radii(self):
"""
Radii of elements in the structure
"""
return self._rad_dict
@property
def struct_valences(self):
"""
Valence of elements in the structure
"""
return self._valence_dict
[docs] def defectsite_count(self):
"""
Returns the number of symmetrically distinct defect sites
"""
return len(self._defect_sites)
[docs] def get_defectsite(self, n):
"""
Returns the defect site at the index.
"""
return self._defect_sites[n]
[docs] def get_defectsite_multiplicity(self, n):
"""
Returns the symmtric multiplicity of the defect site at the index.
"""
return self._defect_site_multiplicity[n]
[docs] def get_defectsite_coordination_number(self, n):
"""
Coordination number of interstitial site.
Args:
n: Index of interstitial list
"""
return self._defectsite_coord_no[n]
[docs] def get_coordinated_sites(self, n):
"""
The sites in structure surrounding the defect site.
Args:
n: Index of defects list
"""
return self._defect_coord_sites[n]
[docs] def get_coordinated_elements(self, n):
"""
Elements of sites in structure surrounding the defect site.
Args:
n: Index of defect list
"""
coordinated_species = []
for site in self._defect_coord_sites[n]:
coordinated_species.append(site.specie.symbol)
return list(set(coordinated_species))
[docs] @abc.abstractmethod
def make_supercells_with_defects(self, scaling_matrix):
"""
Generate the supercell with input multipliers and create the defect.
First supercell has no defects.
To create unit cell with defect pass unit matrix.
"""
raise NotImplementedError()
[docs]class Vacancy(Defect):
"""
Subclass of Defect to generate vacancies and their analysis.
Args:
structure: pymatgen.core.structure.Structure
valences: valences of elements as a dictionary
radii: Radii of elements as a dictionary
"""
def __init__(self, structure, valences, radii):
self._structure = structure
self._valence_dict = valences
self._rad_dict = radii
# Store symmetrically distinct sites, their coordination numbers
# coordinated_sites, effective charge
symm_finder = SpacegroupAnalyzer(self._structure)
symm_structure = symm_finder.get_symmetrized_structure()
equiv_site_seq = symm_structure.equivalent_sites
self._defect_sites = []
self._defect_site_multiplicity = []
for equiv_sites in equiv_site_seq:
self._defect_sites.append(equiv_sites[0])
self._defect_site_multiplicity.append(len(equiv_sites))
self._vac_site_indices = []
for site in self._defect_sites:
for i in range(len(self._structure.sites)):
if site == self._structure[i]:
self._vac_site_indices.append(i)
vnn = VoronoiNN()
self._defectsite_coord_no = []
self._defect_coord_sites = []
for i in self._vac_site_indices:
self._defectsite_coord_no.append(
vnn.get_cn(self._structure, i, use_weights=True))
self._defect_coord_sites.append(
vnn.get_nn(self._structure, i))
# Store the ionic radii for the elements in the structure
# (Used to computing the surface are and volume)
# Computed based on valence of each element
self._vac_eff_charges = None
self._vol = None
self._sa = None
# @property
# def valence_dict(self):
# return self._valence_dict
[docs] def enumerate_defectsites(self):
"""
Returns symmetrically distinct vacancy sites
"""
return self._defect_sites
[docs] def get_defectsite_structure_indices(self):
"""
Returns indices of symmetrically distinct vacancy sites
"""
return self._vac_site_indices
[docs] def get_defectsite_structure_index(self, n):
"""
index of the vacacy site in the structure.sites list
Args:
n:
Index of vacancy list
"""
return self._vac_site_indices[n]
[docs] def get_defectsite_effective_charge(self, n):
"""
Effective charge (In Kroger-Vink notation, cation vacancy has
effectively -ve charge and anion vacancy has +ve charge.)
Args:
n: Index of vacancy list
Returns:
Effective charnge of defect site
"""
# Effective charge (In Kroger-Vink notation, cation vacancy has
# effectively -ve charge and anion vacancy has +ve charge.) Inverse
# the BVAnalyzer.get_valences result.
el = self.get_defectsite(n).species_string
return -self._valence_dict[el]
# if not self._vac_eff_charges:
# self._vac_eff_charges = []
# for site in self.enumerate_defectsites():
# specie = site.specie.symbol
# self._vac_eff_charges.append(-self._valence_dict[specie])
# return self._vac_eff_charges[n]
[docs] def get_coordsites_min_max_charge(self, n):
"""
Minimum and maximum charge of sites surrounding the vacancy site.
Args:
n: Index of vacancy list
"""
bv = BVAnalyzer()
struct_valences = bv.get_valences(self._structure)
coordinated_site_valences = []
def _get_index(site):
for i in range(len(self._structure.sites)):
if site.is_periodic_image(self._structure.sites[i]):
return i
raise ValueError("Site not found")
for site in self._defect_coord_sites[n]:
ind = _get_index(site)
coordinated_site_valences.append(struct_valences[ind])
coordinated_site_valences.sort()
return coordinated_site_valences[0], coordinated_site_valences[-1]
# deprecated
[docs] def get_volume(self, n):
"""
Volume of the nth vacancy
Args:
n: Index of symmetrically distinct vacancies list
Returns:
floating number representing volume of vacancy
"""
if not self._vol:
self._vol = []
self._sa = []
um = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
sc = self.make_supercells_with_defects(um)[1:]
rad_dict = self.struct_radii
for i in range(len(sc)):
vol, sa = get_void_volume_surfarea(sc[i], rad_dict)
self._vol.append(vol)
self._sa.append(sa)
return self._vol[n]
# deprecated
[docs] def get_surface_area(self, n):
"""
Surface area of the nth vacancy
Args:
n: Index of symmetrically distinct vacancies list
Returns:
floating number representing volume of vacancy
"""
if not self._sa:
self._vol = []
self._sa = []
um = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
supercells = self.make_supercells_with_defects(um)[1:]
rad_dict = self.struct_radii
for sc in supercells:
vol, sa = get_void_volume_surfarea(sc, rad_dict)
self._vol.append(vol)
self._sa.append(sa)
return self._sa[n]
def _supercell_with_defect(self, scaling_matrix, defect_site):
sc = self._structure.copy()
sc.make_supercell(scaling_matrix)
oldf_coords = defect_site.frac_coords
coords = defect_site.lattice.get_cartesian_coords(oldf_coords)
newf_coords = sc.lattice.get_fractional_coords(coords)
sc_defect_site = PeriodicSite(defect_site.species_and_occu, newf_coords,
sc.lattice,
properties=defect_site.properties)
for i in range(len(sc.sites)):
# if sc_defect_site == sc.sites[i]:
if sc_defect_site.distance(sc.sites[i]) < 1e-3:
del sc[i]
return sc
raise ValueError('Something wrong if reached here')
[docs] def make_supercells_with_defects(self, scaling_matrix, species=None,
limit_return_structures=False):
"""
Generate sequence of supercells in pymatgen.core.structure.Structure
format, with each supercell containing one vacancy.
Args:
scaling_matrix: super cell scale parameters in matrix forms
species: Species in list format only for which vacancy supercells
are required. If not specified all the species are considered.
limit_return_structures: Boolean or positive number
If number, only that many structures are returned.
Returns:
Supercells with vacancies. First supercell has no defects.
"""
sc_with_vac = []
sc = self._structure.copy()
sc.make_supercell(scaling_matrix)
sc_with_vac.append(sc)
if not species:
species = sc.symbol_set
if not limit_return_structures:
limit_return_structures = self.defectsite_count()
for defect_site in self.enumerate_defectsites():
if len(sc_with_vac) <= limit_return_structures:
if isinstance(defect_site.specie, Specie):
site_specie = defect_site.specie.element.symbol
elif isinstance(defect_site.specie, Element):
site_specie = defect_site.specie.symbol
else:
raise TypeError("site specie is neither Specie nor Element")
if site_specie in species:
sc_with_vac.append(self._supercell_with_defect(
scaling_matrix, defect_site))
return sc_with_vac
[docs]class Interstitial(Defect):
"""
Subclass of Defect to generate interstitial sites
"""
def __init__(self, structure, valences, radii, site_type='voronoi_vertex',
accuracy='Normal', symmetry_flag=True, oxi_state=False):
"""
Given a structure, generate symmetrically distinct interstitial sites.
For a non-ionic structure, use oxi_state=True and give atomic radii.
Args:
structure: pymatgen.core.structure.Structure
valences: Dictionary of oxidation states of elements in
{el:valence} form
radii: Radii of elemnts in the structure
site_type: "voronoi_vertex" uses voronoi nodes
"voronoi_edgecenter" uses voronoi polyhedra edge centers
"voronoi_facecenter" uses voronoi polyhedra face centers
"all" combines vertices, edgecenters and facecenters.
Default is "voronoi_vertex"
accuracy: Flag denoting whether to use high accuracy version
of Zeo++. Options are "Normal" and "High". Default is normal.
symmetry_flag: If True, only returns symmetrically distinct sites
oxi_state: If False, input structure is considered devoid of
oxidation-state decoration. And oxi-state for each site is
determined. Use True, if input structure is oxi-state
decorated. This option is useful when the structure is
not electro-neutral after deleting/adding sites. In that
case oxi-decorate the structure before deleting/adding the
sites.
"""
if not oxi_state:
self._structure = ValenceIonicRadiusEvaluator(structure).structure
else:
self._structure = structure
self._valence_dict = valences
self._rad_dict = radii
"""
Use Zeo++ to obtain the voronoi nodes. Apply symmetry reduction
and the symmetry reduced voronoi nodes are possible candidates
for interstitial sites.
"""
if accuracy == "Normal":
high_accuracy_flag = False
elif accuracy == "High":
high_accuracy_flag = True
else:
raise NotImplementedError("Accuracy setting not implemented.")
if accuracy == "High":
if site_type in ('voronoi_facecenter', 'voronoi_edgecenter', 'all'):
raise NotImplementedError(
"Site type not implemented for the accuracy setting")
vor_node_sites, vor_edgecenter_sites, vor_facecenter_sites = \
symmetry_reduced_voronoi_nodes(self._structure, self._rad_dict,
high_accuracy_flag, symmetry_flag)
if site_type == 'voronoi_vertex':
possible_interstitial_sites = vor_node_sites
elif site_type == 'voronoi_facecenter':
possible_interstitial_sites = vor_facecenter_sites
elif site_type == 'voronoi_edgecenter':
possible_interstitial_sites = vor_edgecenter_sites
elif site_type == "all":
possible_interstitial_sites = vor_node_sites + \
vor_facecenter_sites + vor_edgecenter_sites
else:
raise ValueError("Input site type not implemented")
# Do futher processing on possibleInterstitialSites to obtain
# interstitial sites
self._defect_sites = possible_interstitial_sites
self._defectsite_coord_no = []
self._defect_coord_sites = []
self._defect_coord_charge = []
self._radii = []
for site in self._defect_sites:
coord_no, coord_sites, chrg = self._get_coord_no_sites_chrg(site)
self._defectsite_coord_no.append(coord_no)
self._defect_coord_sites.append(coord_sites)
self._defect_coord_charge.append(chrg)
for site in self._defect_sites:
vor_radius = site.properties.get('voronoi_radius', None)
if vor_radius:
vor_radius = float(vor_radius)
self._radii.append(vor_radius)
def _get_coord_no_sites_chrg(self, site):
"""
Compute the coordination number and coordination charge
Args:
site:
pymatgen.core.sites.Site
"""
struct = self._structure.copy()
struct.append(site.specie.symbol, site.frac_coords)
vnn = VoronoiNN()
coord_no = vnn.get_cn(struct, -1, use_weights=True)
coord_sites = vnn.get_nn(struct, -1)
# In some cases coordination sites to interstitials include
# interstitials also. Filtering them.
def no_inter(site):
return not site.specie.symbol == 'X'
coord_sites = filter(no_inter, coord_sites)
coord_chrg = 0
if self._valence_dict:
for site, weight in vnn.get_voronoi_polyhedra(struct, -1).items():
if not site.specie.symbol == 'X':
coord_chrg += weight * self._valence_dict[
site.species_string]
return coord_no, coord_sites, coord_chrg
[docs] def enumerate_defectsites(self):
"""
Enumerate all the symmetrically distinct interstitial sites.
The defect site has "X" as occupied specie.
"""
return self._defect_sites
[docs] def append_defectsite(self, site):
"""
Append a site to list of possible interstitials
Args:
site: pymatgen.core.sites.Site
"""
raise NotImplementedError()
[docs] def delete_defectsite(self, n):
"""
Remove a symmetrically distinct interstitial site
Args:
n: Index of interstitial site
"""
del self._defect_sites[n]
[docs] def get_coordsites_charge_sum(self, n):
"""
Total charge of the interstitial coordinated sites.
Args:
n: Index of interstitial list
"""
return self._defect_coord_charge[n]
[docs] def get_coordsites_min_max_charge(self, n):
"""
Minimum and maximum charge of sites surrounding the interstitial site.
Args:
n: Index of symmetrical distinct interstitial site
"""
coord_site_valences = []
for site in self._defect_coord_sites[n]:
coord_site_valences.append(self._valence_dict[site.specie.symbol])
coord_site_valences.sort()
return coord_site_valences[0], coord_site_valences[-1]
[docs] def get_radius(self, n):
"""
Volume of the nth interstitial
Args:
n: Index of symmetrically distinct vacancies list
Returns:
floating number representing radius of interstitial sphere
"""
return self._radii[n]
[docs] def get_radii(self):
return self._radii
[docs] def reduce_defectsites(self):
"""
If multiple defect sites have same voronoi radius, only one is kept.
Useful if the symmetry based reduction of initial sites returned
from Zeo++ is not working properly due to deviation in ideal lattice
coordinates.
"""
distinct_radii = list(set(self._radii))
for rad in distinct_radii:
ind = self._radii.index(rad) # Index of first site with 'rad'
for i in reversed(list(range(ind + 1, len(self._radii)))):
# Backward search for remaining sites so index is not changed
if self._radii[i] == rad:
self._defect_sites.pop(i)
self._defectsite_coord_no.pop(i)
self._defect_coord_sites.pop(i)
self._radii.pop(i)
[docs] def radius_prune_defectsites(self, radius):
"""
Remove all the defect sites with voronoi radius less than input radius
"""
for i in reversed(list(range(len(self._radii)))):
if self._radii[i] < radius:
self._defect_sites.pop(i)
self._defectsite_coord_no.pop(i)
self._defect_coord_sites.pop(i)
self._radii.pop(i)
[docs] def prune_defectsites(self, el="C", oxi_state=4, dlta=0.1):
"""
Prune all the defect sites which can't acoomodate the input elment
with the input oxidation state.
"""
rad = float(Specie(el, oxi_state).ionic_radius) - dlta
self.radius_prune_defectsites(rad)
[docs] def prune_close_defectsites(self, dist=0.2):
"""
Prune the sites that are very close.
"""
# print self.defectsite_count()
ind = 0
while ind < self.defectsite_count():
# i = ind + 1
# while i < self.defectsite_count():
i = self.defectsite_count() - 1
# print ind, i
while i > ind:
d = self._defect_sites[ind].distance(self._defect_sites[i])
# print d, dist
if d < dist:
self._defect_sites.pop(i)
# self._defectsite_coord_no.pop(i)
# self._defect_coord_sites.pop(i)
# self._radii.pop(i)
# i += 1
i -= 1
ind += 1
# print self.defectsite_count()
def _supercell_with_defect(self, scaling_matrix, defect_site, element):
sc = self._structure.copy()
sc.make_supercell(scaling_matrix)
oldf_coords = defect_site.frac_coords
coords = defect_site.lattice.get_cartesian_coords(oldf_coords)
# print coords
newf_coords = sc.lattice.get_fractional_coords(coords)
for i in range(3):
coord = newf_coords[i]
if coord < 0:
while (coord < 0):
coord = coord + 1
newf_coords[i] = coord
elif coord > 1:
while (coord > 1):
coord = coord - 1
newf_coords[i] = coord
# print newf_coords
# sc_defect_site = PeriodicSite(element, newf_coords,
# sc.lattice)
try:
sc.append(element, newf_coords, coords_are_cartesian=False,
validate_proximity=True)
except ValueError:
sc = None
finally:
return sc
[docs] def make_supercells_with_defects(self, scaling_matrix, element):
"""
Returns sequence of supercells in pymatgen.core.structure.Structure
format, with each supercell containing an interstitial.
First supercell has no defects.
"""
sc_list_with_interstitial = []
sc = self._structure.copy()
sc.make_supercell(scaling_matrix)
sc_list_with_interstitial.append(sc)
for defect_site in self.enumerate_defectsites():
sc_with_inter = self._supercell_with_defect(
scaling_matrix, defect_site, element
)
if sc_with_inter:
sc_list_with_interstitial.append(sc_with_inter)
return sc_list_with_interstitial
[docs]class InterstitialAnalyzer(object):
"""
Use GULP to compute the interstitial formation energy, relaxed structures.
Works only for metal oxides due to the use of Buckingham Potentials.
Args:
inter: pymatgen.defects.point_defects.Interstitial
el: Element name in short hand notation ("El")
oxi_state: Oxidtation state
scd: Super cell dimension as number. The scaling is equal along xyz.
"""
def __init__(self, inter, el, oxi_state, scd=2):
self._inter = inter
self._el = el
self._oxi_state = oxi_state
self._scd = scd
self._relax_energies = []
self._norelax_energies = []
self._relax_struct = []
[docs] def get_energy(self, n, relax=True):
"""
Formation Energy for nth symmetrically distinct interstitial.
"""
if relax and not self._relax_energies:
self._relax_analysis()
if not relax and not self._norelax_energies:
no_inter = self._inter.defectsite_count()
inter_gulp_kw = ('qok',)
val_dict = self._inter.struct_valences
val_dict[self._el] = self._oxi_state # If element not in structure
scd = self._scd
scale_mat = [[scd, 0, 0], [0, scd, 0], [0, 0, scd]]
sc = self._inter.make_supercells_with_defects(scale_mat, self._el)
blk_energy = get_energy_buckingham(sc[0])
for i in range(1, no_inter + 1):
inter_energy = get_energy_buckingham(
sc[i], keywords=inter_gulp_kw, valence_dict=val_dict
)
form_energy = inter_energy - blk_energy
self._norelax_energies.append(form_energy)
if relax:
return self._relax_energies[n]
else:
return self._norelax_energies[n]
def _relax_analysis(self):
"""
Optimize interstitial structures
"""
no_inter = self._inter.defectsite_count()
inter_gulp_kw = ('optimise', 'conp', 'qok')
val_dict = self._inter.struct_valences
scd = self._scd
scale_mat = [[scd, 0, 0], [0, scd, 0], [0, 0, scd]]
sc = self._inter.make_supercells_with_defects(scale_mat, self._el)
blk_energy, rlx_struct = get_energy_relax_structure_buckingham(sc[0])
self._relax_struct.append(rlx_struct)
val_dict[self._el] = self._oxi_state # If element not in structure
for i in range(1, no_inter + 1):
energy, rlx_struct = get_energy_relax_structure_buckingham(
sc[i], keywords=inter_gulp_kw, valence_dict=val_dict
)
form_energy = energy - blk_energy
self._relax_energies.append(form_energy)
self._relax_struct.append(rlx_struct)
[docs] def get_relaxed_structure(self, n):
"""
Optimized interstitial structure
Args:
n: Symmetrically distinct interstitial index
.. note::
To get relaxed bulk structure pass -1.
-ve index will not work as expected.
"""
if not self._relax_struct:
self._relax_analysis()
return self._relax_struct[n + 1]
[docs] def get_percentage_volume_change(self, n):
"""
Volume change after the introduction of interstitial
Args:
n: Symmetrically distinct interstitial index
"""
if not self._relax_struct:
self._relax_analysis()
blk_struct = self._relax_struct[0]
def_struct = self._relax_struct[n + 1:n + 2][0]
del def_struct.sites[-1]
rv = RelaxationAnalyzer(blk_struct, def_struct)
return rv.get_percentage_volume_change()
[docs] def get_percentage_lattice_parameter_change(self, n):
"""
Lattice parameter change after the introduction of interstitial
Args:
n: Symmetrically distinct interstitial index
"""
if not self._relax_struct:
self._relax_analysis()
blk_struct = self._relax_struct[0]
def_struct = self._relax_struct[n + 1:n + 2][0]
del def_struct.sites[-1]
rv = RelaxationAnalyzer(blk_struct, def_struct)
return rv.get_percentage_lattice_parameter_changes()
[docs] def get_percentage_bond_distance_change(self, n):
"""
Bond distance change after the introduction of interstitial
Args:
n: Symmetrically distinct interstitial index
"""
if not self._relax_struct:
self._relax_analysis()
blk_struct = self._relax_struct[0]
def_struct = self._relax_struct[n + 1:n + 2][0]
del def_struct.sites[-1]
# print def_struct
rv = RelaxationAnalyzer(blk_struct, def_struct)
return rv.get_percentage_bond_dist_changes()
[docs] def relaxed_structure_match(self, i, j):
"""
Check if the relaxed structures of two interstitials match
Args:
i: Symmetrically distinct interstitial index
j: Symmetrically distinct interstitial index
.. note::
To use relaxed bulk structure pass -1.
-ve index will not work as expected
"""
if not self._relax_struct:
self._relax_analysis()
sm = StructureMatcher()
struct1 = self._relax_struct[i + 1]
struct2 = self._relax_struct[j + 1]
return sm.fit(struct1, struct2)
[docs]class StructureRelaxer(object):
def __init__(self, structure):
self._unrelax_struct = structure
self.relax()
[docs] def relax(self):
energy, rlx_struct = get_energy_relax_structure_buckingham(
self._unrelax_struct)
self._relax_struct = rlx_struct
[docs] def get_relaxed_structure(self):
return self._relax_struct
[docs]class InterstitialStructureRelaxer(object):
"""
Performs structural relaxation for each interstitial supercell.
Args:
interstitial: Unrelaxed interstitial
el: Species string in short notation
oxi_state: Oxidation state of the element
supercell_dim: Dimensions of super cell
"""
def __init__(self, interstitial, el, oxi_state, supercell_dim=2):
self._inter = interstitial
self._scd = supercell_dim
self._el = el
self._oxi_state = oxi_state
self._relax_structs = []
self._relax_energies = []
[docs] def relax(self):
"""
Optimize interstitial structures
"""
no_inter = self._inter.defectsite_count()
inter_gulp_kw = ('optimise', 'conp', 'qok')
val_dict = self._inter.struct_valences
scd = self._scd
scale_mat = [[scd, 0, 0], [0, scd, 0], [0, 0, scd]]
sc = self._inter.make_supercells_with_defects(scale_mat, self._el)
blk_energy, rlx_struct = get_energy_relax_structure_buckingham(sc[0])
self._relax_structs.append(rlx_struct)
self._relax_energies.append(blk_energy)
val_dict[self._el] = self._oxi_state # If element not in structure
for i in range(1, no_inter + 1):
try:
energy, rlx_struct = get_energy_relax_structure_buckingham(
sc[i], keywords=inter_gulp_kw, valence_dict=val_dict
)
self._relax_energies.append(energy)
self._relax_structs.append(rlx_struct)
except:
self._relax_energies.append(None)
self._relax_structs.append(None)
def is_empty(lst):
for value in lst:
if value:
return False
return True
if is_empty(self._relax_energies):
raise IOError('Relaxation failed')
[docs] def relaxed_structure_match(self, i, j):
"""
Check if the relaxed structures of two interstitials match
Args:
i: Symmetrically distinct interstitial index
j: Symmetrically distinct interstitial index
.. note::
Index 0 corresponds to bulk.
"""
if not self._relax_structs:
self.relax()
sm = StructureMatcher()
struct1 = self._relax_structs[i]
struct2 = self._relax_structs[j]
return sm.fit(struct1, struct2)
[docs] def relaxed_energy_match(self, i, j):
"""
Check if the relaxed energies of two interstitials match
Args:
i: Symmetrically distinct interstitial index
j: Symmetrically distinct interstitial index
.. note::
Index 0 corresponds to bulk.
"""
if not self._relax_energies:
self.relax()
energy1 = self._relax_energies[i]
energy2 = self._relax_energies[j]
return energy1 == energy2
[docs] def get_relaxed_structure(self, n):
"""
Get the relaxed structure of nth symmetrically distinct interstitial.
Args:
n: Symmetrically distinct interstitial index
.. note::
0 corresponds to relaxed bulk structure
"""
if not self._relax_structs:
self.relax()
return self._relax_structs[n]
[docs] def get_relaxed_energy(self, n):
"""
Get the relaxed structure of nth symmetrically distinct interstitial.
Args:
n: Symmetrically distinct interstitial index
.. note::
0 corresponds to relaxed bulk structure
"""
if not self._relax_energies:
self.relax()
return self._relax_energies[n]
[docs] def get_relaxed_interstitial(self):
"""
Get the relaxed structure of nth symmetrically distinct interstitial.
Args:
n: Symmetrically distinct interstitial index
"""
if not self._relax_energies:
self.relax()
energies = self._relax_energies[:]
structs = self._relax_structs[:]
distinct_energy_set = set(energies[1:]) # only interstitial energies
if None in distinct_energy_set:
distinct_energy_set.remove(None)
distinct_structs = [structs[0]] # bulk
distinct_energies = [energies[0]]
for energy in distinct_energy_set:
ind = energies.index(energy)
distinct_structs.append(structs[ind])
distinct_energies.append(energies[ind])
return RelaxedInterstitial(
distinct_structs, distinct_energies, self._inter.struct_valences
)
[docs]class RelaxedInterstitial(object):
"""
Stores the relaxed supercell structures for each interstitial
Used to compute formation energies, displacement of atoms near the
the interstitial.
Args:
struct_list: List of structures(supercells). The first structure should
represent relaxed bulk structure and the subsequent ones
interstitial structures (with the extra interstitial site
appended at the end).
energy_list: List of energies for relaxed interstitial structures.
The first energy should correspond to bulk structure
valence_dict: Valences of elements in dictionary form
"""
def __init__(self, struct_list, energy_list, valence_dict):
self._blk_struct = struct_list[0]
struct_list.pop(0)
self._structs = struct_list
self._blk_energy = energy_list[0]
energy_list.pop(0)
self._energies = energy_list
self._valence_dict = valence_dict
self._coord_no = []
self._coord_sites = []
self._coord_charge_no = []
[docs] def get_percentage_volume_change(self, n):
"""
Volume change after the introduction of interstitial
Args:
n: index of interstitials
"""
def_struct = self._structs[n:n + 1][0]
del def_struct.sites[-1]
rv = RelaxationAnalyzer(self._blk_struct, def_struct)
return rv.get_percentage_volume_change()
[docs] def get_percentage_lattice_parameter_change(self, n):
"""
Lattice parameter change after the introduction of interstitial
Args:
n: index of interstitials
"""
def_struct = self._structs[n:n + 1][0] # copy
del def_struct.sites[-1]
rv = RelaxationAnalyzer(self._blk_struct, def_struct)
return rv.get_percentage_lattice_parameter_changes()
[docs] def get_percentage_bond_distance_change(self, n):
"""
Bond distance change after the introduction of interstitial.
Args:
n: index of interstitials
"""
def_struct = self._structs[n:n + 1][0] # copy
del def_struct.sites[-1]
rv = RelaxationAnalyzer(self._blk_struct, def_struct)
return rv.get_percentage_bond_dist_changes()
[docs] def get_bulk_structure(self):
"""
Return relaxed bulk structure
"""
return self._blk_struct
[docs] def get_interstitial_structure(self, n):
"""
Return relaxed bulk structure
"""
return self._structs[n]
[docs] def defect_count(self):
"""
Returns the number of distinct interstitials
"""
return len(self._structs)
[docs] def get_defectsite(self, n):
"""
Returns the defect site of nth interstitial.
Args:
n: Index of interstitial
"""
return self._structs[n][-1]
[docs] def get_coordination_number(self, n):
"""
Coordination number for nth interstitial.
Args:
n: Index of interstitials
"""
if not self._coord_no:
self._coord_find()
return self._coord_no[n]
[docs] def get_charge_coordination_number(self, n):
"""
Charge coordination number for nth interstitial.
Args:
n: Index of interstitials
"""
if not self._coord_charge_no:
self._coord_find()
return self._coord_charge_no[n]
[docs] def get_coordinated_sites(self, n):
"""
Coordinated sites for nth interstitial.
Args:
n: Index of interstitials
"""
if not self._coord_sites:
self._coord_find()
return self._coord_sites[n]
[docs] def get_coordinated_bulk_sites(self, n):
"""
Bulk sites corresponding to the coordinated sites for nth interstitial.
Args:
n: Index of interstitials
"""
blk_sites = []
for site in self.get_coordinated_sites(n):
site_index = self._structs[n].sites.index(site)
blk_sites.append(self._blk_struct[site_index])
return blk_sites
[docs] def get_coordinated_site_displacement(self, n):
"""
Compute the total displacement of coordinated sites from the
interstitial sites during the relaxation
Args:
n: Index of defect site
"""
coord_sites = self.get_coordinated_sites(n)
coord_blk_sites = self.get_coordinated_bulk_sites(n)
dist_sum = 0
for i in range(len(coord_sites)):
dist_sum += coord_sites[i].distance_from_point(coord_blk_sites[i])
# How to compute the average?
return dist_sum
def _coord_find(self):
"""
calls VoronoiNN to compute the coordination number,
coordination charge
"""
for i in range(self.defect_count()):
struct = self._structs[i].copy()
vnn = VoronoiNN()
self._coord_no.append(vnn.get_cn(struct, -1, use_weights=True))
self._coord_sites.append(vnn.get_nn(struct, -1))
coord_chrg = 0
for site, weight in vnn.get_voronoi_polyhedra(struct, -1).items():
coord_chrg += weight * self._valence_dict[site.species_string]
self._coord_charge_no.append(coord_chrg)
[docs]def symmetry_reduced_voronoi_nodes(
structure, rad_dict, high_accuracy_flag=False, symm_flag=True):
"""
Obtain symmetry reduced voronoi nodes using Zeo++ and
pymatgen.symmetry.finder.SpacegroupAnalyzer
Args:
strucutre: pymatgen Structure object
rad_dict: Dictionary containing radii of spcies in the structure
high_accuracy_flag: Flag denotting whether to use high accuracy version of Zeo++
symm_flag: Flag denoting whether to return symmetrically distinct sites only
Returns:
Symmetrically distinct voronoi nodes as pymatgen Strucutre
"""
def add_closest_equiv_site(dist_sites, equiv_sites):
if not dist_sites:
dist_sites.append(equiv_sites[0])
else:
avg_dists = []
for site in equiv_sites:
dists = [site.distance(dst_site, jimage=[0, 0, 0])
for dst_site in dist_sites]
avg_dist = sum(dists) / len(dist_sites)
avg_dists.append(avg_dist)
min_avg_dist = min(avg_dists)
ind = avg_dists.index(min_avg_dist)
dist_sites.append(equiv_sites[ind])
def cmp_memoize_last_site(f): # Compares and stores last site
def not_duplicates(site1, site2):
if site1.distance(site2) < 1e-5:
return False
else:
return True
cmp_memoize_last_site.cache = None
def helper(x):
if not cmp_memoize_last_site.cache:
cmp_memoize_last_site.cache = f(x)
return True
y = f(x)
if not_duplicates(cmp_memoize_last_site.cache, y):
cmp_memoize_last_site.cache = y
return True
else:
return False
return helper
@cmp_memoize_last_site
def check_not_duplicates(site):
return site
if not symm_flag:
if not high_accuracy_flag:
vor_node_struct, vor_edgecenter_struct, vor_facecenter_struct = \
get_voronoi_nodes(structure, rad_dict)
return vor_node_struct.sites, vor_edgecenter_struct.sites, \
vor_facecenter_struct.sites
else:
# Only the nodes are from high accuracy voronoi decomposition
vor_node_struct = \
get_high_accuracy_voronoi_nodes(structure, rad_dict)
# Before getting the symmetry, remove the duplicates
vor_node_struct.sites.sort(key=lambda site: site.voronoi_radius)
# print type(vor_node_struct.sites[0])
dist_sites = filter(check_not_duplicates, vor_node_struct.sites)
return dist_sites, None, None
if not high_accuracy_flag:
def get_dist_sites(vor_struct):
SpgA = SpacegroupAnalyzer
try:
symmetry_finder = SpgA(vor_struct, symprec=1e-1)
symm_struct = symmetry_finder.get_symmetrized_structure()
except:
vor_struct.merge_sites(0.1, 'delete')
symmetry_finder = SpgA(vor_struct, symprec=1e-1)
symm_struct = symmetry_finder.get_symmetrized_structure()
equiv_sites_list = symm_struct.equivalent_sites
if not equiv_sites_list:
dist_sites = vor_struct.sites
else:
dist_sites = []
for equiv_sites in equiv_sites_list:
add_closest_equiv_site(dist_sites, equiv_sites)
return dist_sites
vor_node_struct, vor_edgecenter_struct, vor_facecenter_struct = \
get_voronoi_nodes(structure, rad_dict)
node_dist_sites = get_dist_sites(vor_node_struct)
edgecenter_dist_sites = get_dist_sites(vor_edgecenter_struct)
facecenter_dist_sites = get_dist_sites(vor_facecenter_struct)
return node_dist_sites, edgecenter_dist_sites, facecenter_dist_sites
else:
# Only the nodes are from high accuracy voronoi decomposition
vor_node_struct = \
get_high_accuracy_voronoi_nodes(structure, rad_dict)
# Before getting the symmetry, remove the duplicates
vor_node_struct.sites.sort(key=lambda site: site.voronoi_radius)
# print type(vor_node_struct.sites[0])
dist_sites = list(filter(check_not_duplicates, vor_node_struct.sites))
# Ignore symmetry from ha voronoi nodes
# Increase the symmetry precision to 0.25
# spg = SpacegroupAnalyzer(structure,symprec=1e-1).get_spacegroup()
# Remove symmetrically equivalent sites
# i = 0
# while (i < len(dist_sites)-1):
# sites1 = [dist_sites[i]]
# sites2 = [dist_sites[i+1]]
# if spg.are_symmetrically_equivalent(sites1,sites2):
# del dist_sites[i+1]
# else:
# i = i+1
node_dist_sites = dist_sites
return (node_dist_sites, None, None)
# vor_edge_symmetry_finder = SpacegroupAnalyzer(
# vor_edgecenter_struct, symprec=1e-1)
# vor_edge_symm_struct = vor_edge_symmetry_finder.get_symmetrized_structure()
# edgecenter_equiv_sites_list = vor_edge_symm_struct.equivalent_sites
# edgecenter_dist_sites = []
# for equiv_sites in edgecenter_equiv_sites_list:
# add_closest_equiv_site(edgecenter_dist_sites, equiv_sites)
# if not edgecenter_equiv_sites_list:
# edgecenter_dist_sites = vor_edgecenter_struct.sites
# vor_fc_symmetry_finder = SpacegroupAnalyzer(
# vor_facecenter_struct, symprec=1e-1)
# vor_fc_symm_struct = vor_fc_symmetry_finder.get_symmetrized_structure()
# facecenter_equiv_sites_list = vor_fc_symm_struct.equivalent_sites
# facecenter_dist_sites = []
# for equiv_sites in facecenter_equiv_sites_list:
# add_closest_equiv_site(facecenter_dist_sites, equiv_sites)
# if not facecenter_equiv_sites_list:
# facecenter_dist_sites = vor_facecenter_struct.sites
# return node_dist_sites, edgecenter_dist_sites, facecenter_dist_sites
[docs]class StructureMotifInterstitial(Defect):
"""
Generate interstitial sites at positions
where the interstitialcy is coordinated by nearest neighbors
in a way that resembles basic structure motifs
(e.g., tetrahedra, octahedra). The algorithm is called InFiT
(Interstitialcy Finding Tool), it was introducted by
Nils E. R. Zimmermann, Matthew K. Horton, Anubhav Jain,
and Maciej Haranczyk (Front. Mater., 4, 34, 2017),
and it is used by the Python Charged Defect Toolkit
(PyCDT: D. Broberg et al., Comput. Phys. Commun., in press, 2018).
"""
def __init__(self, struct, inter_elem,
motif_types=("tetrahedral", "octahedral"),
op_threshs=(0.3, 0.5),
dl=0.2, doverlap=1, facmaxdl=1.01, verbose=False):
"""
Generates symmetrically distinct interstitial sites at positions
where the interstitial is coordinated by nearest neighbors
in a pattern that resembles a supported structure motif
(e.g., tetrahedra, octahedra).
Args:
struct (Structure): input structure for which symmetrically
distinct interstitial sites are to be found.
inter_elem (string): element symbol of desired interstitial.
motif_types ([string]): list of structure motif types that are
to be considered. Permissible types are:
tet (tetrahedron), oct (octahedron).
op_threshs ([float]): threshold values for the underlying order
parameters to still recognize a given structural motif
(i.e., for an OP value >= threshold the coordination pattern
match is positive, for OP < threshold the match is
negative.
dl (float): grid fineness in Angstrom. The input
structure is divided into a grid of dimension
a/dl x b/dl x c/dl along the three crystallographic
directions, with a, b, and c being the lengths of
the three lattice vectors of the input unit cell.
doverlap (float): distance that is considered
to flag an overlap between any trial interstitial site
and a host atom.
facmaxdl (float): factor to be multiplied with the maximum grid
width that is then used as a cutoff distance for the
clustering prune step.
verbose (bool): flag indicating whether (True) or not (False;
default) to print additional information to screen.
"""
# Initialize interstitial finding.
self._structure = struct.copy()
self._motif_types = motif_types[:]
if len(self._motif_types) == 0:
raise RuntimeError("no motif types provided.")
self._op_threshs = op_threshs[:]
self.cn_motif_lostop = {}
self.target_cns = []
for imotif, motif in enumerate(self._motif_types):
if motif not in list(motif_cn_op.keys()):
raise RuntimeError("unsupported motif type: {}.".format(
motif))
cn = int(motif_cn_op[motif]['cn'])
if cn not in self.target_cns:
self.target_cns.append(cn)
if cn not in list(self.cn_motif_lostop.keys()):
self.cn_motif_lostop[cn] = {}
tmp_optype = motif_cn_op[motif]['optype']
if tmp_optype == 'tet_max':
tmp_optype = 'tet'
if tmp_optype == 'oct_max':
tmp_optype = 'oct'
self.cn_motif_lostop[cn][motif] = LocalStructOrderParams(
[tmp_optype],
parameters=[motif_cn_op[motif]['params']], cutoff=-10.0)
self._dl = dl
self._defect_sites = []
self._defect_types = []
self._defect_site_multiplicity = []
self._defect_cns = []
self._defect_opvals = []
rots, trans = SpacegroupAnalyzer(
struct)._get_symmetry()
nbins = [int(struct.lattice.a / dl), \
int(struct.lattice.b / dl), \
int(struct.lattice.c / dl)]
dls = [struct.lattice.a / float(nbins[0]), \
struct.lattice.b / float(nbins[1]), \
struct.lattice.c / float(nbins[2])]
maxdl = max(dls)
if verbose:
print("Grid size: {} {} {}".format(nbins[0], nbins[1], nbins[2]))
print("dls: {} {} {}".format(dls[0], dls[1], dls[2]))
struct_w_inter = struct.copy()
struct_w_inter.append(inter_elem, [0, 0, 0])
natoms = len(list(struct_w_inter.sites))
trialsites = []
# Loop over trial positions that are based on a regular
# grid in fractional coordinate space
# within the unit cell.
for ia in range(nbins[0]):
a = (float(ia) + 0.5) / float(nbins[0])
for ib in range(nbins[1]):
b = (float(ib) + 0.5) / float(nbins[1])
for ic in range(nbins[2]):
c = (float(ic) + 0.5) / float(nbins[2])
struct_w_inter.replace(
natoms - 1, inter_elem, coords=[a, b, c],
coords_are_cartesian=False)
if len(struct_w_inter.get_sites_in_sphere(
struct_w_inter.sites[natoms - 1].coords,
doverlap)) == 1:
neighs_images_weigths = MinimumDistanceNN(
tol=0.8, cutoff=6).get_nn_info(
struct_w_inter, natoms - 1)
neighs_images_weigths_sorted = sorted(
neighs_images_weigths, key=lambda x: x['weight'],
reverse=True)
for nsite in range(1, len(neighs_images_weigths_sorted)+1):
if nsite not in self.target_cns:
continue
allsites = [neighs_images_weigths_sorted[i]['site'] for i in range(nsite)]
indices_neighs = [i for i in range(len(allsites))]
allsites.append(struct_w_inter.sites[natoms - 1])
for mot, ops in self.cn_motif_lostop[nsite].items():
opvals = ops.get_order_parameters(
allsites, len(allsites) - 1,
indices_neighs=indices_neighs)
if opvals[0] > op_threshs[motif_types.index(mot)]:
cns = {}
for isite in range(nsite):
site = neighs_images_weigths_sorted[isite]['site']
if isinstance(site.specie, Element):
elem = site.specie.symbol
else:
elem = site.specie.element.symbol
if elem in list(cns.keys()):
cns[elem] = cns[elem] + 1
else:
cns[elem] = 1
trialsites.append({
"mtype": mot,
"opval": opvals[0],
"coords": struct_w_inter.sites[
natoms - 1].coords[:],
"fracs": np.array([a, b, c]),
"cns": dict(cns)})
break
# Prune list of trial sites by clustering and find the site
# with the largest order parameter value in each cluster.
nintersites = len(trialsites)
unique_motifs = []
for ts in trialsites:
if ts["mtype"] not in unique_motifs:
unique_motifs.append(ts["mtype"])
labels = {}
connected = []
for i in range(nintersites):
connected.append([])
for j in range(nintersites):
dist, image = struct_w_inter.lattice.get_distance_and_image(
trialsites[i]["fracs"],
trialsites[j]["fracs"])
connected[i].append(
True if dist < (maxdl * facmaxdl) else False)
include = []
for motif in unique_motifs:
labels[motif] = []
for i, ts in enumerate(trialsites):
labels[motif].append(i if ts["mtype"] == motif else -1)
change = True
while change:
change = False
for i in range(nintersites - 1):
if change:
break
if labels[motif][i] == -1:
continue
for j in range(i + 1, nintersites):
if labels[motif][j] == -1:
continue
if connected[i][j] and labels[motif][i] != \
labels[motif][j]:
if labels[motif][i] < labels[motif][j]:
labels[motif][j] = labels[motif][i]
else:
labels[motif][i] = labels[motif][j]
change = True
break
unique_ids = []
for l in labels[motif]:
if l != -1 and l not in unique_ids:
unique_ids.append(l)
if verbose:
print("unique_ids {} {}".format(motif, unique_ids))
for uid in unique_ids:
maxq = 0.0
imaxq = -1
for i in range(nintersites):
if labels[motif][i] == uid:
if imaxq < 0 or trialsites[i]["opval"] > maxq:
imaxq = i
maxq = trialsites[i]["opval"]
include.append(imaxq)
# Prune by symmetry.
multiplicity = {}
discard = []
for motif in unique_motifs:
discard_motif = []
for indi, i in enumerate(include):
if trialsites[i]["mtype"] != motif or \
i in discard_motif:
continue
multiplicity[i] = 1
symposlist = [trialsites[i]["fracs"].dot(
np.array(m, dtype=float)) for m in rots]
for t in trans:
symposlist.append(trialsites[i]["fracs"] + np.array(t))
for indj in range(indi + 1, len(include)):
j = include[indj]
if trialsites[j]["mtype"] != motif or \
j in discard_motif:
continue
for sympos in symposlist:
dist, image = struct.lattice.get_distance_and_image(
sympos, trialsites[j]["fracs"])
if dist < maxdl * facmaxdl:
discard_motif.append(j)
multiplicity[i] += 1
break
for i in discard_motif:
if i not in discard:
discard.append(i)
if verbose:
print("Initial trial sites: {}\nAfter clustering: {}\n"
"After symmetry pruning: {}".format(
len(trialsites), len(include),
len(include) - len(discard)))
c = 0
for i in include:
if i not in discard:
self._defect_sites.append(
PeriodicSite(
Element(inter_elem),
trialsites[i]["fracs"],
self._structure.lattice,
to_unit_cell=False,
coords_are_cartesian=False,
properties=None))
self._defect_types.append(trialsites[i]["mtype"])
self._defect_cns.append(trialsites[i]["cns"])
self._defect_site_multiplicity.append(multiplicity[i])
self._defect_opvals.append(trialsites[i]["opval"])
[docs] def enumerate_defectsites(self):
"""
Get all defect sites.
Returns:
defect_sites ([PeriodicSite]): list of periodic sites
representing the interstitials.
"""
return self._defect_sites
[docs] def get_motif_type(self, i):
"""
Get the motif type of defect with index i (e.g., "tet").
Returns:
motif (string): motif type.
"""
return self._defect_types[i]
[docs] def get_coordinating_elements_cns(self, i):
"""
Get element-specific coordination numbers of defect with index i.
Returns:
elem_cn (dict): dictionary storing the coordination numbers (int)
with string representation of elements as keys.
(i.e., {elem1 (string): cn1 (int), ...}).
"""
return self._defect_cns[i]
[docs] def get_op_value(self, i):
"""
Get order-parameter value of defect with index i.
Returns:
opval (float): OP value.
"""
return self._defect_opvals[i]
[docs] def make_supercells_with_defects(self, scaling_matrix):
"""
Generate a sequence of supercells
in which each supercell contains a single interstitial,
except for the first supercell in the sequence
which is a copy of the defect-free input structure.
Args:
scaling_matrix (3x3 integer array): scaling matrix
to transform the lattice vectors.
Returns:
scs ([Structure]): sequence of supercells.
"""
scs = []
sc = self._structure.copy()
sc.make_supercell(scaling_matrix)
scs.append(sc)
for ids, defect_site in enumerate(self._defect_sites):
sc_with_inter = sc.copy()
sc_with_inter.append(defect_site.species_string,
defect_site.frac_coords,
coords_are_cartesian=False,
validate_proximity=False,
properties=None)
if not sc_with_inter:
raise RuntimeError("could not generate supercell with"
" interstitial {}".format(ids + 1))
scs.append(sc_with_inter.copy())
return scs
[docs]class TopographyAnalyzer(object):
"""
This is a generalized module to perform topological analyses of a crystal
structure using Voronoi tessellations. It can be used for finding potential
interstitial sites. Applications including using these sites for
inserting additional atoms or for analyzing diffusion pathways.
Note that you typically want to do some preliminary postprocessing after
the initial construction. The initial construction will create a lot of
points, especially for determining potential insertion sites. Some helper
methods are available to perform aggregation and elimination of nodes. A
typical use is something like::
a = TopographyAnalyzer(structure, ["O"], ["P"])
a.cluster_nodes()
a.remove_collisions()
"""
def __init__(self, structure, framework_ions, cations, tol=0.0001,
max_cell_range=1, check_volume=True,
constrained_c_frac=0.5, thickness=0.5):
"""
Init.
Args:
structure (Structure): An initial structure.
framework_ions ([str]): A list of ions to be considered as a
framework. Typically, this would be all anion species. E.g.,
["O", "S"].
cations ([str]): A list of ions to be considered as non-migrating
cations. E.g., if you are looking at Li3PS4 as a Li
conductor, Li is a mobile species. Your cations should be [
"P"]. The cations are used to exclude polyhedra from
diffusion analysis since those polyhedra are already occupied.
tol (float): A tolerance distance for the analysis, used to
determine if something are actually periodic boundary images of
each other. Default is usually fine.
max_cell_range (int): This is the range of periodic images to
construct the Voronoi tesselation. A value of 1 means that we
include all points from (x +- 1, y +- 1, z+- 1) in the
voronoi construction. This is because the Voronoi poly
extends beyond the standard unit cell because of PBC.
Typically, the default value of 1 works fine for most
structures and is fast. But for really small unit
cells with high symmetry, you may need to increase this to 2
or higher.
check_volume (bool): Set False when ValueError always happen after
tuning tolerance.
constrained_c_frac (float): Constraint the region where users want
to do Topology analysis the default value is 0.5, which is the
fractional coordinate of the cell
thickness (float): Along with constrained_c_frac, limit the
thickness of the regions where we want to explore. Default is
0.5, which is mapping all the site of the unit cell.
"""
self.structure = structure
self.framework_ions = set([get_el_sp(sp) for sp in framework_ions])
self.cations = set([get_el_sp(sp) for sp in cations])
# Let us first map all sites to the standard unit cell, i.e.,
# 0 ≤ coordinates < 1.
# structure = Structure.from_sites(structure, to_unit_cell=True)
# lattice = structure.lattice
# We could constrain the region where we want to dope/explore by setting
# the value of constrained_c_frac and thickness. The default mode is
# mapping all sites to the standard unit cell
s = structure.copy()
constrained_sites = []
for i, site in enumerate(s):
if site.frac_coords[2] >= constrained_c_frac - thickness and \
site.frac_coords[
2] <= constrained_c_frac + thickness:
constrained_sites.append(site)
structure = Structure.from_sites(sites=constrained_sites)
lattice = structure.lattice
# Divide the sites into framework and non-framework sites.
framework = []
non_framework = []
for site in structure:
if self.framework_ions.intersection(site.species_and_occu.keys()):
framework.append(site)
else:
non_framework.append(site)
# We construct a supercell series of coords. This is because the
# Voronoi polyhedra can extend beyond the standard unit cell. Using a
# range of -2, -1, 0, 1 should be fine.
coords = []
cell_range = list(range(-max_cell_range, max_cell_range + 1))
for shift in itertools.product(cell_range, cell_range, cell_range):
for site in framework:
shifted = site.frac_coords + shift
coords.append(lattice.get_cartesian_coords(shifted))
# Perform the voronoi tessellation.
voro = Voronoi(coords)
# Store a mapping of each voronoi node to a set of points.
node_points_map = defaultdict(set)
for pts, vs in voro.ridge_dict.items():
for v in vs:
node_points_map[v].update(pts)
logger.debug("%d total Voronoi vertices" % len(voro.vertices))
# Vnodes store all the valid voronoi polyhedra. Cation vnodes store
# the voronoi polyhedra that are already occupied by existing cations.
vnodes = []
cation_vnodes = []
def get_mapping(poly):
"""
Helper function to check if a vornoi poly is a periodic image
of one of the existing voronoi polys.
"""
for v in vnodes:
if v.is_image(poly, tol):
return v
return None
# Filter all the voronoi polyhedra so that we only consider those
# which are within the unit cell.
for i, vertex in enumerate(voro.vertices):
if i == 0:
continue
fcoord = lattice.get_fractional_coords(vertex)
poly = VoronoiPolyhedron(lattice, fcoord, node_points_map[i],
coords, i)
if np.all([-tol <= c < 1 + tol for c in fcoord]):
if len(vnodes) == 0:
vnodes.append(poly)
else:
ref = get_mapping(poly)
if ref is None:
vnodes.append(poly)
logger.debug("%d voronoi vertices in cell." % len(vnodes))
# Eliminate all voronoi nodes which are closest to existing cations.
if len(cations) > 0:
cation_coords = [
site.frac_coords for site in non_framework
if self.cations.intersection(site.species_and_occu.keys())]
vertex_fcoords = [v.frac_coords for v in vnodes]
dist_matrix = lattice.get_all_distances(cation_coords,
vertex_fcoords)
indices = np.where(dist_matrix == np.min(dist_matrix, axis=1)[
:, None])[1]
cation_vnodes = [v for i, v in enumerate(vnodes) if i in indices]
vnodes = [v for i, v in enumerate(vnodes) if i not in indices]
logger.debug("%d vertices in cell not with cation." % len(
vnodes))
self.coords = coords
self.vnodes = vnodes
self.cation_vnodes = cation_vnodes
self.framework = framework
self.non_framework = non_framework
if check_volume:
self.check_volume()
[docs] def check_volume(self):
# Basic check for volume of all voronoi poly sum to unit cell volume
# Note that this does not apply after poly combination.
vol = sum((v.volume for v in self.vnodes)) + sum((v.volume for v in
self.cation_vnodes))
if abs(vol - self.structure.volume) > 1e-8:
raise ValueError(
"Sum of voronoi volumes is not equal to original volume of "
"structure! This may lead to inaccurate results. You need to "
"tweak the tolerance and max_cell_range until you get a "
"correct mapping.")
[docs] def cluster_nodes(self, tol=0.2):
"""
Cluster nodes that are too close together using a tol.
Args:
tol (float): A distance tolerance. PBC is taken into account.
"""
lattice = self.structure.lattice
vfcoords = [v.frac_coords for v in self.vnodes]
# Manually generate the distance matrix (which needs to take into
# account PBC.
dist_matrix = np.array(lattice.get_all_distances(vfcoords, vfcoords))
dist_matrix = (dist_matrix + dist_matrix.T) / 2
for i in range(len(dist_matrix)):
dist_matrix[i, i] = 0
condensed_m = squareform(dist_matrix)
z = linkage(condensed_m)
cn = fcluster(z, tol, criterion="distance")
merged_vnodes = []
for n in set(cn):
poly_indices = set()
frac_coords = []
for i, j in enumerate(np.where(cn == n)[0]):
poly_indices.update(self.vnodes[j].polyhedron_indices)
if i == 0:
frac_coords.append(self.vnodes[j].frac_coords)
else:
fcoords = self.vnodes[j].frac_coords
# We need the image to combine the frac_coords properly.
d, image = lattice.get_distance_and_image(frac_coords[0],
fcoords)
frac_coords.append(fcoords + image)
merged_vnodes.append(VoronoiPolyhedron(
lattice,
np.average(frac_coords, axis=0),
poly_indices,
self.coords
))
self.vnodes = merged_vnodes
logger.debug("%d vertices after combination." %
len(self.vnodes))
[docs] def remove_collisions(self, min_dist=0.5):
"""
Remove vnodes that are too close to existing atoms in the structure
Args:
min_dist(float): The minimum distance that a vertex needs to be
from existing atoms.
"""
vfcoords = [v.frac_coords for v in self.vnodes]
sfcoords = self.structure.frac_coords
dist_matrix = self.structure.lattice.get_all_distances(vfcoords,
sfcoords)
all_dist = np.min(dist_matrix, axis=1)
new_vnodes = []
for i, v in enumerate(self.vnodes):
if all_dist[i] > min_dist:
new_vnodes.append(v)
self.vnodes = new_vnodes
[docs] def get_structure_with_nodes(self):
"""
Get the modified structure with the voronoi nodes inserted. The
species is set as a DummySpecie X.
"""
new_s = Structure.from_sites(self.structure)
for v in self.vnodes:
new_s.append("X", v.frac_coords)
return new_s
[docs] def print_stats(self):
latt = self.structure.lattice
def get_min_dist(fcoords):
n = len(fcoords)
dist = latt.get_all_distances(fcoords, fcoords)
all_dist = [dist[i, j] for i in range(n) for j in range(i + 1, n)]
return min(all_dist)
voro = [s[1] for s in self.vertices]
print("Min dist between voronoi vertices centers = %.4f" %
get_min_dist(voro))
def get_non_framework_dist(fcoords):
cations = [site.frac_coords for site in self.non_framework]
dist_matrix = latt.get_all_distances(cations, fcoords)
min_dist = np.min(dist_matrix, axis=1)
assert len(cations) == len(min_dist)
return np.linalg.norm(min_dist), min(min_dist), max(min_dist)
print(len(self.non_framework))
print("MSE dist voro = %s" % str(get_non_framework_dist(voro)))
[docs] def write_topology(self, fname="Topo.cif"):
new_s = Structure.from_sites(self.structure)
for v in self.vnodes:
new_s.append("Mg", v.frac_coords)
new_s.to(filename=fname)
[docs] def analyze_symmetry(self, tol):
s = Structure.from_sites(self.framework)
site_to_vindex = {}
for i, v in enumerate(self.vnodes):
s.append("Li", v.frac_coords)
site_to_vindex[s[-1]] = i
print(len(s))
finder = SpacegroupAnalyzer(s, tol)
print(finder.get_space_group_operations())
symm_structure = finder.get_symmetrized_structure()
print(len(symm_structure.equivalent_sites))
return [[site_to_vindex[site] for site in sites]
for sites in symm_structure.equivalent_sites if sites[
0].specie.symbol == "Li"]
[docs] def vtk(self):
if StructureVis is None:
raise NotImplementedError("vtk must be present to view.")
lattice = self.structure.lattice
vis = StructureVis()
vis.set_structure(Structure.from_sites(self.structure))
for v in self.vnodes:
vis.add_site(PeriodicSite("K", v.frac_coords, lattice))
vis.add_polyhedron(
[PeriodicSite("S", c, lattice, coords_are_cartesian=True)
for c in v.polyhedron_coords],
PeriodicSite("Na", v.frac_coords, lattice),
color="element",
draw_edges=True,
edges_color=(0, 0, 0))
vis.show()
[docs]class VoronoiPolyhedron(object):
"""
Convenience container for a voronoi point in PBC and its associated polyhedron.
"""
def __init__(self, lattice, frac_coords, polyhedron_indices, all_coords,
name=None):
self.lattice = lattice
self.frac_coords = frac_coords
self.polyhedron_indices = polyhedron_indices
self.polyhedron_coords = np.array(all_coords)[list(polyhedron_indices),
:]
self.name = name
[docs] def is_image(self, poly, tol):
frac_diff = pbc_diff(poly.frac_coords, self.frac_coords)
if not np.allclose(frac_diff, [0, 0, 0], atol=tol):
return False
to_frac = lambda c: self.lattice.get_fractional_coords(c)
for c1 in self.polyhedron_coords:
found = False
for c2 in poly.polyhedron_coords:
d = pbc_diff(to_frac(c1), to_frac(c2))
if not np.allclose(d, [0, 0, 0], atol=tol):
found = True
break
if not found:
return False
return True
@property
def coordination(self):
return len(self.polyhedron_indices)
@property
def volume(self):
return calculate_vol(self.polyhedron_coords)
def __str__(self):
return "Voronoi polyhedron %s" % self.name
[docs]class ChargeDensityAnalyzer(object):
"""
Analyzer to find potential interstitial sites based on charge density. The
`total` charge density is used.
"""
def __init__(self, chgcar):
"""
Initialization.
Args:
chgcar (pmg.Chgcar): input Chgcar object.
"""
self.chgcar = chgcar
self.structure = chgcar.structure
self.extrema_coords = [] # list of frac_coords of local extrema
self.extrema_type = None # "local maxima" or "local minima"
self._extrema_df = None # extrema frac_coords - chg density table
self._charge_distribution_df = None # frac_coords - chg density table
[docs] @classmethod
def from_file(cls, chgcar_filename):
chgcar = Chgcar.from_file(chgcar_filename)
return cls(chgcar=chgcar)
@property
def charge_distribution_df(self):
if self._charge_distribution_df is None:
return self._get_charge_distribution_df()
else:
return self._charge_distribution_df
@property
def extrema_df(self):
if self.extrema_type is None:
logger.warning("Please run ChargeDensityAnalyzer.get_local_extrema first!")
return self._extrema_df
def _get_charge_distribution_df(self):
"""
Return a complete table of fractional coordinates - charge density.
"""
# Fraction coordinates and corresponding indices
axis_grid = np.array(
[np.array(self.chgcar.get_axis_grid(i)) /
self.structure.lattice.abc[i] for i in range(3)])
axis_index = np.array([range(len(axis_grid[i])) for i in range(3)])
data = {}
for index in itertools.product(*axis_index):
a, b, c = index
f_coords = (axis_grid[0][a], axis_grid[1][b], axis_grid[2][c])
data[f_coords] = self.chgcar.data["total"][a][b][c]
# Fraction coordinates - charge density table
df = pd.Series(data).reset_index()
df.columns = ['a', 'b', 'c', 'Charge Density']
self._charge_distribution_df = df
return df
def _update_extrema(self, f_coords, extrema_type,
threshold_frac=None, threshold_abs=None):
"""Update _extrema_df, extrema_type and extrema_coords"""
if threshold_frac is not None:
if threshold_abs is not None:
logger.warning( # Exit if both filter are set
"Filter can be either threshold_frac or threshold_abs!")
return
assert 0 <= threshold_frac <= 1, "threshold_frac range is [0, 1]!"
# Return empty result if coords list is empty
if len(f_coords) == 0:
df = pd.DataFrame({}, columns=['A', 'B', 'C', "Chgcar"])
self._extrema_df = df
self.extrema_coords = []
logger.info("Find {} {}.".format(len(df), extrema_type))
return
data = {}
unit = 1 / np.array(self.chgcar.dim) # pixel along a, b, c
for fc in f_coords:
a, b, c = tuple(map(int, fc / unit))
data[tuple(fc)] = self.chgcar.data["total"][a][b][c]
df = pd.Series(data).reset_index()
df.columns = ['a', 'b', 'c', 'Charge Density']
ascending = (extrema_type == "local minima")
if threshold_abs is None:
threshold_frac = threshold_frac \
if threshold_frac is not None else 1.0
num_extrema = int(threshold_frac * len(f_coords))
df = df.sort_values(by="Charge Density",
ascending=ascending)[0: num_extrema]
df.reset_index(drop=True, inplace=True) # reset major index
else: # threshold_abs is set
df = df.sort_values(by="Charge Density", ascending=ascending)
df = df[df["Charge Density"] <= threshold_abs] if ascending \
else df[df["Charge Density"] >= threshold_abs]
extrema_coords = []
for row in df.iterrows():
fc = np.array(row[1]["a": "c"])
extrema_coords.append(fc)
self._extrema_df = df
self.extrema_type = extrema_type
self.extrema_coords = extrema_coords
logger.info("Find {} {}.".format(len(df), extrema_type))
[docs] @requires(peak_local_max_found,
"get_local_extrema requires skimage.feature.peak_local_max module"
" to be installed. Please confirm your skimage installation.")
def get_local_extrema(self, find_min=True,
threshold_frac=None, threshold_abs=None):
"""
Get all local extrema fractional coordinates in charge density,
searching for local minimum by default. Note that sites are NOT grouped
symmetrically.
Args:
find_min (bool): True to find local minimum else maximum, otherwise
find local maximum.
threshold_frac (float): optional fraction of extrema shown, which
returns `threshold_frac * tot_num_extrema` extrema fractional
coordinates based on highest/lowest intensity.
E.g. set 0.2 to show the extrema with 20% highest or lowest
intensity. Value range: 0 <= threshold_frac <= 1
Note that threshold_abs and threshold_frac should not set in the
same time.
threshold_abs (float): optional filter. When searching for local
minima, intensity <= threshold_abs returns; when searching for
local maxima, intensity >= threshold_abs returns.
Note that threshold_abs and threshold_frac should not set in the
same time.
Returns:
extrema_coords (list): list of fractional coordinates corresponding
to local extrema.
"""
sign, extrema_type = 1, "local maxima"
if find_min:
sign, extrema_type = -1, "local minima"
# Make 3x3x3 supercell
# This is a trick to resolve the periodical boundary issue.
total_chg = sign * self.chgcar.data["total"]
total_chg = np.tile(total_chg, reps=(3, 3, 3))
coordinates = peak_local_max(total_chg, min_distance=1)
# Remove duplicated sites introduced by supercell.
f_coords = [coord / total_chg.shape * 3 for coord in coordinates]
f_coords = [f - 1 for f in f_coords if
all(np.array(f) < 2) and all(np.array(f) >= 1)]
# Update information
self._update_extrema(f_coords, extrema_type,
threshold_frac=threshold_frac,
threshold_abs=threshold_abs)
return self.extrema_coords
[docs] def cluster_nodes(self, tol=0.2):
"""
Cluster nodes that are too close together using a tol.
Args:
tol (float): A distance tolerance. PBC is taken into account.
"""
lattice = self.structure.lattice
vf_coords = self.extrema_coords
if len(vf_coords) == 0:
if self.extrema_type is None:
logger.warning(
"Please run ChargeDensityAnalyzer.get_local_extrema first!")
return
new_f_coords = []
self._update_extrema(new_f_coords, self.extrema_type)
return new_f_coords
# Manually generate the distance matrix (which needs to take into
# account PBC.
dist_matrix = np.array(lattice.get_all_distances(vf_coords, vf_coords))
dist_matrix = (dist_matrix + dist_matrix.T) / 2
for i in range(len(dist_matrix)):
dist_matrix[i, i] = 0
condensed_m = squareform(dist_matrix)
z = linkage(condensed_m)
cn = fcluster(z, tol, criterion="distance")
merged_fcoords = []
for n in set(cn):
frac_coords = []
for i, j in enumerate(np.where(cn == n)[0]):
if i == 0:
frac_coords.append(self.extrema_coords[j])
else:
f_coords = self.extrema_coords[j]
# We need the image to combine the frac_coords properly.
d, image = lattice.get_distance_and_image(frac_coords[0],
f_coords)
frac_coords.append(f_coords + image)
merged_fcoords.append(np.average(frac_coords, axis=0))
self._update_extrema(merged_fcoords, extrema_type=self.extrema_type)
logger.debug(
"{} vertices after combination.".format(len(self.extrema_coords)))
[docs] def remove_collisions(self, min_dist=0.5):
"""
Remove predicted sites that are too close to existing atoms in the
structure.
Args:
min_dist (float): The minimum distance (in Angstrom) that
a predicted site needs to be from existing atoms. A min_dist
with value <= 0 returns all sites without distance checking.
"""
s_f_coords = self.structure.frac_coords
f_coords = self.extrema_coords
if len(f_coords) == 0:
if self.extrema_type is None:
logger.warning(
"Please run ChargeDensityAnalyzer.get_local_extrema first!")
return
new_f_coords = []
self._update_extrema(new_f_coords, self.extrema_type)
return new_f_coords
dist_matrix = self.structure.lattice.get_all_distances(f_coords,
s_f_coords)
all_dist = np.min(dist_matrix, axis=1)
new_f_coords = []
for i, f in enumerate(f_coords):
if all_dist[i] > min_dist:
new_f_coords.append(f)
self._update_extrema(new_f_coords, self.extrema_type)
return new_f_coords
[docs] def get_structure_with_nodes(self, find_min=True, min_dist=0.5, tol=0.2,
threshold_frac=None, threshold_abs=None):
"""
Get the modified structure with the possible interstitial sites added.
The species is set as a DummySpecie X.
Args:
find_min (bool): True to find local minimum else maximum, otherwise
find local maximum.
min_dist (float): The minimum distance (in Angstrom) that
a predicted site needs to be from existing atoms. A min_dist
with value <= 0 returns all sites without distance checking.
tol (float): A distance tolerance of nodes clustering that sites too
closed to other predicted sites will be merged. PBC is taken
into account.
threshold_frac (float): optional fraction of extrema, which returns
`threshold_frac * tot_num_extrema` extrema fractional
coordinates based on highest/lowest intensity.
E.g. set 0.2 to insert DummySpecie atom at the extrema with 20%
highest or lowest intensity.
Value range: 0 <= threshold_frac <= 1
Note that threshold_abs and threshold_frac should not set in the
same time.
threshold_abs (float): optional filter. When searching for local
minima, intensity <= threshold_abs returns; when searching for
local maxima, intensity >= threshold_abs returns.
Note that threshold_abs and threshold_frac should not set in the
same time.
Returns:
structure (Structure)
"""
structure = self.structure.copy()
self.get_local_extrema(find_min=find_min,
threshold_frac=threshold_frac,
threshold_abs=threshold_abs)
self.remove_collisions(min_dist)
self.cluster_nodes(tol=tol)
for fc in self.extrema_coords:
structure.append("X", fc)
return structure
[docs]def calculate_vol(coords):
if len(coords) == 4:
coords_affine = np.ones((4, 4))
coords_affine[:, 0:3] = np.array(coords)
return abs(np.linalg.det(coords_affine)) / 6
else:
simplices = get_facets(coords, joggle=True)
center = np.average(coords, axis=0)
vol = 0
for s in simplices:
c = list(s.coords)
c.append(center)
vol += calculate_vol(c)
return vol