Source code for matador.fingerprints.pdf

# coding: utf-8
# Distributed under the terms of the MIT License.

""" This submodule defines classes for computing, combining
and convolving pair distribution functions.


from itertools import combinations_with_replacement
import itertools
import copy
from math import ceil, copysign
import time

import numpy as np
import numba

from matador.utils.cell_utils import frac2cart, cart2volume
from matador.utils.cell_utils import standardize_doc_cell
from matador.utils.chem_utils import get_stoich
from matador.fingerprints.fingerprint import Fingerprint, FingerprintFactory

[docs]class PDF(Fingerprint): """This class implements the calculation and comparison of pair distribution functions. Attributes: r_space (ndarray) : 1-D array containing real space grid gr (ndarray): 1-D array containing total PDF dr (float): real-space grid spacing in Å rmax (float): extent of real-space grid in Å label (str): structure label elem_gr (dict): dict with pairs of element symbol keys, containing 1-D arrays of projected PDFs (if calculated) number_density (float): number density for renormalisation and comparison with other PDFs kwargs (dict): arguments used to create PDF """ def __init__(self, doc, lazy=False, **kwargs): """Initialise parameters and run PDF (unless lazy=True). Parameters: doc (dict) : matador document to calculate PDF of Keyword Arguments: dr (float) : bin width for PDF (Angstrom) (DEFAULT: 0.01) gaussian_width (float) : width of Gaussian smearing (Angstrom) (DEFAULT: 0.01) num_images (int/str) : number of unit cell images include in PDF calculation (DEFAULT: 'auto') max_num_images (int) : cutoff number of unit cells before crashing (DEFAULT: 50) rmax (float) : maximum distance cutoff for PDF (Angstrom) (DEFAULT: 15) projected (bool) : optionally calculate the element-projected PDF standardize (bool) : standardize cell before calculating PDF lazy (bool) : if True, calculator is not called when initializing PDF object timing (bool) : if True, print the total time taken to calculate the PDF """ prop_defaults = { "dr": 0.01, "gaussian_width": 0.1, "rmax": 15, "num_images": "auto", "style": "smear", "debug": False, "timing": False, "low_mem": False, "projected": True, "max_num_images": 50, "standardize": True, } # read and store kwargs self.kwargs = prop_defaults self.kwargs.update( {key: kwargs[key] for key in kwargs if kwargs[key] is not None} ) # useful data for labelling structure = copy.deepcopy(doc) self.spg = structure.get("space_group", "") if self.kwargs.get("standardize"): structure = standardize_doc_cell(structure) self.spg = structure["space_group"] self.stoichiometry = structure.get( "stoichiometry", get_stoich(structure["atom_types"]) ) # private variables self._num_images = self.kwargs.get("num_images") self._lattice = np.asarray(structure["lattice_cart"]) self._poscart = np.asarray( frac2cart(structure["lattice_cart"], structure["positions_frac"]) ).reshape(-1, 3) self._types = structure["atom_types"] self._num_atoms = len(self._poscart) self._volume = cart2volume(self._lattice) self._image_vec = None # public variables self.rmax = self.kwargs.get("rmax") self.number_density = self._num_atoms / self._volume self.dr = self.kwargs.get("dr") self.r_space = None = None self.elem_gr = None self.label = None if self.kwargs.get("label"): self.label = self.kwargs["label"] elif "text_id" in structure: self.label = " ".join(structure["text_id"]) if not lazy: if self.kwargs.get("timing"): start = time.time() self.calc_pdf() if self.kwargs.get("timing"): end = time.time() print("PDF calculated in {:.3f} s".format(end - start))
[docs] def calc_pdf(self): """Wrapper to calculate PDF with current settings.""" if self._image_vec is None: if self.kwargs.get("debug"): start = time.time() self._set_image_trans_vectors() if self.kwargs.get("debug"): end = time.time() print("Image vectors length = {}".format(len(self._image_vec))) print("Set image trans vectors in {} s".format(end - start)) if self.kwargs.get("projected"): if self.kwargs.get("debug"): start = time.time() if self.elem_gr is None: self._calc_projected_pdf() if self.kwargs.get("debug"): end = time.time() print("Calculated projected PDF {} s".format(end - start)) if is None: self._calc_unprojected_pdf()
[docs] def calculate(self): """Alias for `self.calc_pdf`.""" self.calc_pdf()
def _calc_distances(self, poscart, poscart_b=None): """Calculate PBC distances with cdist. Parameters: poscart (numpy.ndarray): array of absolute atomic coordinates. Keyword Arguments: poscart_b (numpy.ndarray): absolute positions of a second type of atoms, where only A-B distances will be calculated. Returns: numpy.ndarray: pair d_ij matrix with values > rmax < 1e-12 removed. """ from matador.utils.cell_utils import calc_pairwise_distances_pbc return calc_pairwise_distances_pbc( poscart, self._image_vec, self._lattice, self.rmax, filter_zero=True, compress=True, poscart_b=poscart_b, debug=self.kwargs.get("debug"), ) def _calc_unprojected_pdf(self): """Wrapper function to calculate distances and output a broadened and normalised PDF. Sets and self.r_space to G(r) and r respectively. """ if self.elem_gr is not None: self._calc_unprojected_pdf_from_projected() else: distances = self._calc_distances(self._poscart) self.r_space = np.arange(0, self.rmax + self.dr, self.dr) = self._get_broadened_normalised_pdf( distances, style=self.kwargs.get("style"), gaussian_width=self.kwargs.get("gaussian_width"), ) def _calc_projected_pdf(self): """Calculate broadened and normalised element-projected PDF of a matador document. Sets self.elem_gr of e.g. Li2Zn3 to { ('Li', 'Li'): G_{Li-Li}(r), ('Li', 'Zn'): G_{Li-Zn}(r), ('Zn', 'Zn'): G_{Zn-Zn}(r) } """ # initalise dict of element pairs with correct keys style = self.kwargs.get("style") gw = self.kwargs.get("gaussian_width") self.r_space = np.arange(0, self.rmax + self.dr, self.dr) elem_gr = dict() for comb in combinations_with_replacement(set(self._types), 2): elem_gr[tuple(set(comb))] = np.zeros_like(self.r_space) for elem_type in elem_gr: poscart = [ self._poscart[i] for i in range(len(self._poscart)) if self._types[i] == elem_type[0] ] poscart_b = ( [ self._poscart[i] for i in range(len(self._poscart)) if self._types[i] == elem_type[1] ] if len(elem_type) == 2 else None ) distances = self._calc_distances(poscart, poscart_b=poscart_b) elem_gr[elem_type] = len(elem_type) * self._get_broadened_normalised_pdf( distances, style=style, gaussian_width=gw ) self.elem_gr = elem_gr def _calc_unprojected_pdf_from_projected(self): """ " Reconstruct full PDF from projected.""" = np.zeros_like(self.r_space) for key in self.elem_gr: += self.elem_gr[key] @staticmethod @numba.njit def _normalize_gr(gr, r_space, dr, num_atoms, number_density): """Normalise a broadened PDF, ignoring the Gaussian magnitude.""" norm = 4 * np.pi * (r_space + dr) ** 2 * dr * num_atoms * number_density return np.divide(gr, norm) @staticmethod @numba.njit def _dist_hist(distances, r_space, dr): """Bin the pair-wise distances according to the radial grid. Parameters: distances (numpy.ndarray): array of pair-wise distances. r_space (numpy.ndarray): radial grid dr (float): bin width. """ hist = np.zeros_like(r_space) for dij in distances: hist[ceil(dij / dr)] += 1 return hist def _get_broadened_normalised_pdf( self, distances, style="smear", gaussian_width=0.1 ): """Broaden the values provided as distances and return G(r) and r_space of the normalised PDF. Parameters: distances (numpy.ndarray): distances used to calculate PDF Keyword arguments: style (str): either 'smear' or 'histogram' gaussian_width (float): smearing width in Angstrom^1/2 Returns: gr (np.ndarray): G(r), the PDF of supplied distances """ if style == "histogram" or gaussian_width == 0: gr = self._dist_hist(distances, self.r_space, self.dr) else: # otherwise do normal smearing hist = self._dist_hist(distances, self.r_space, self.dr) gr = self._broadening_unrolled(hist, self.r_space, gaussian_width) gr = self._normalize_gr( gr, self.r_space, self.dr, self._num_atoms, self.number_density ) return gr @staticmethod def _get_image_trans_vectors_auto(lattice, rmax, dr, max_num_images=50): """Finds all "images" (integer 3-tuples, supercells) that have atoms within rmax + dr + longest LV of the parent lattice. Parameters: lattice (list): list of lattice vectors. rmax (float): maximum radial distance. dr (float): PDF bin width. Keyword arguments: max_num_images (int): the greatest integer multiple of LVs to search out to. Returns: list: list of int 3-tuples of cells, up to rmax or if max_num_images is exceeded, just up to 1 cell away. """ image_vec = set() # find longest combination of single LV's max_trans = 0 _lattice = np.asarray(lattice) products = list(itertools.product(range(-1, 2), repeat=3)) for prod in products: trans = prod @ _lattice length = np.sqrt(np.sum(trans**2)) if length > max_trans: max_trans = length unit_vector_lengths = np.sqrt(np.sum(_lattice**2, axis=1)) limits = [ int((dr + rmax + max_trans) / length) for length in unit_vector_lengths ] for ind, limit in enumerate(limits): if abs(limit) > max_num_images: limits[ind] = int(copysign(max_num_images, limit)) products = itertools.product(*(range(-lim, lim + 1) for lim in limits)) for prod in products: trans = prod @ _lattice length = np.sqrt(np.sum(trans**2)) if length <= rmax + dr + max_trans: image_vec.add(prod) return image_vec def _set_image_trans_vectors(self): """Sets self._image_vec to a list/generator of image translation vectors, based on self._num_images. If self._num_images is an integer, create all 3-member integer combinations up to the value. If self._num_images is 'auto', create all translation vectors up to length self.rmax. e.g. self._image_vec = [[1, 0, 1], [0, 1, 1], [1, 1, 1]]. """ if self._num_images == "auto": self._image_vec = self._get_image_trans_vectors_auto( self._lattice, self.rmax, self.dr, max_num_images=self.kwargs.get("max_num_images"), ) else: self._image_vec = list( itertools.product( range(-self._num_images, self._num_images + 1), repeat=3 ) )
[docs] def get_sim_distance(self, pdf_b, projected=False): """Return the similarity between two PDFs.""" return PDFOverlap(self, pdf_b, projected=projected).similarity_distance
[docs] def pdf(self): """Return G(r) and the r_space for easy plotting.""" try: return (self.r_space, except AttributeError: return (None, None)
[docs] def plot_projected_pdf(self, **kwargs): """Plot projected PDFs. Keyword arguments: keys (list): plot only a subset of projections, e.g. [('K', )]. other_pdfs (list of PDF): other PDFs to plot. """ from matador.plotting.pdf_plotting import plot_projected_pdf plot_projected_pdf(self, **kwargs)
[docs] def plot(self, **kwargs): """Plot PDFs. Keyword arguments: other_pdfs (list of PDF): other PDFs to add to the plot. """ from matador.plotting.pdf_plotting import plot_pdf plot_pdf(self, **kwargs)
[docs]class PDFFactory(FingerprintFactory): """This class computes PDF objects from a list of structures, as concurrently as possible. The PDFs are stored under the `pdf` key inside each structure dict. Attributes: nprocs (int): number of concurrent processes. """ default_key = "pdf" fingerprint = PDF
[docs]class PDFOverlap: """Calculate the PDFOverlap between two PDF objects, pdf_a and pdf_b, with number density rescaling. Attributes: pdf_a (PDF): first PDF to compare. pdf_b (PDF): second PDF to compare. fine_dr (float): fine grid scale on which to compare. similarity_distance (float): "distance" between PDFs. overlap_int (float): the value of the overlap integral. """ def __init__(self, pdf_a, pdf_b, projected=False): """Perform the overlap and similarity distance calculations. Parameters: pdf_a (PDF): first PDF to compare. pdf_b (PDF): second PDF to compare. Keyword arguments: projected : if True, attempt to use projected PDFs. """ self.pdf_a = pdf_a self.pdf_b = pdf_b self.fine_dr = self.pdf_a.dr / 2.0 # initialise with large number self.similarity_distance = 1e10 self.overlap_int = 0 if projected: if isinstance(pdf_a.elem_gr, dict) and isinstance(pdf_b.elem_gr, dict): self.projected_pdf_overlap() else: print("Projected PDFs missing, continuing with total.") self.pdf_overlap() else: self.pdf_overlap()
[docs] def pdf_overlap(self): """Calculate the overlap of two PDFs via a simple meshed sum of their difference. """ self.overlap_int = 0 self.similarity_distance = 1e10 self.fine_space = np.arange(0, self.pdf_a.rmax, self.fine_dr) self.fine_gr_a = np.interp(self.fine_space, self.pdf_a.r_space, self.fine_gr_b = np.interp(self.fine_space, self.pdf_b.r_space, # scaling factor here is normalising to number density density_rescaling_factor = pow( self.pdf_b.number_density / (self.pdf_a.number_density), 1 / 3 ) rescale_factor = density_rescaling_factor self.fine_gr_a = np.interp( self.fine_space, rescale_factor * self.fine_space, self.fine_gr_a ) self.fine_gr_a = self.fine_gr_a[: int(len(self.fine_space) * 0.75)] self.fine_gr_b = self.fine_gr_b[: int(len(self.fine_space) * 0.75)] self.fine_space = self.fine_space[: int(len(self.fine_space) * 0.75)] overlap_fn = self.fine_gr_a - self.fine_gr_b worst_case_overlap_int = np.trapz( np.abs(self.fine_gr_a), dx=self.pdf_a.dr / 2.0 ) + np.trapz(np.abs(self.fine_gr_b), dx=self.pdf_b.dr / 2.0) self.overlap_int = np.trapz(np.abs(overlap_fn), dx=self.pdf_a.dr / 2.0) self.similarity_distance = self.overlap_int / worst_case_overlap_int self.overlap_fn = overlap_fn
[docs] def projected_pdf_overlap(self): """Calculate the overlap of two projected PDFs via a simple meshed sum of their difference. """ self.fine_space = np.arange(0, self.pdf_a.rmax, self.fine_dr) self.overlap_int = 0 self.similarity_distance = 1e10 elems = set(key for key in self.pdf_a.elem_gr) if elems != set(key for key in self.pdf_b.elem_gr): for key in self.pdf_b.elem_gr: elems.add(key) # pad out missing elements with zero PDFs for key in elems: if key not in self.pdf_a.elem_gr: self.pdf_a.elem_gr[key] = np.zeros_like(self.pdf_a.r_space) if key not in self.pdf_b.elem_gr: self.pdf_b.elem_gr[key] = np.zeros_like(self.pdf_b.r_space) self.fine_elem_gr_a, self.fine_elem_gr_b = dict(), dict() for key in elems: self.fine_elem_gr_a[key] = np.interp( self.fine_space, self.pdf_a.r_space, self.pdf_a.elem_gr[key] ) self.fine_elem_gr_b[key] = np.interp( self.fine_space, self.pdf_b.r_space, self.pdf_b.elem_gr[key] ) # scaling factor here is normalising to number density density_rescaling_factor = pow( (self.pdf_b.number_density) / (self.pdf_a.number_density), 1 / 3 ) rescale_factor = density_rescaling_factor for key in elems: self.fine_elem_gr_a[key] = np.interp( self.fine_space, rescale_factor * self.fine_space, self.fine_elem_gr_a[key], ) for key in elems: self.fine_elem_gr_a[key] = self.fine_elem_gr_a[key][ : int(len(self.fine_space) * 0.75) ] self.fine_elem_gr_b[key] = self.fine_elem_gr_b[key][ : int(len(self.fine_space) * 0.75) ] self.fine_space = self.fine_space[: int(len(self.fine_space) * 0.75)] for key in elems: overlap_fn = self.fine_elem_gr_a[key] - self.fine_elem_gr_b[key] worst_case_a = np.trapz( np.abs(self.fine_elem_gr_a[key]), dx=self.pdf_a.dr / 2.0 ) worst_case_b = np.trapz( np.abs(self.fine_elem_gr_b[key]), dx=self.pdf_b.dr / 2.0 ) worst_case_overlap = worst_case_a + worst_case_b overlap = np.trapz(np.abs(overlap_fn), dx=self.pdf_a.dr / 2.0) self.overlap_int += overlap / worst_case_overlap self.similarity_distance = self.overlap_int / len(elems)
[docs] def plot_diff(self): """Simple plot for comparing two PDFs.""" from matador.plotting.pdf_plotting import plot_diff_overlap plot_diff_overlap(self)
[docs] def plot_projected_diff(self): """Simple plot for comparing two projected PDFs.""" from matador.plotting.pdf_plotting import plot_projected_diff_overlap plot_projected_diff_overlap(self)
[docs]class CombinedProjectedPDF: """Take some computed PDFs and add them together.""" def __init__(self, pdf_cursor): """Create CombinedPDF object from list of PDFs. Parameters: pdf_cursor (:obj:`list` of :obj:`PDF`): list of PDF objects to combine. """ self.dr = min([pdf.dr for pdf in pdf_cursor]) self.rmax = min([pdf.rmax for pdf in pdf_cursor]) self.r_space = np.arange(0, self.rmax + self.dr, self.dr) self.label = "Combined PDF" if any([not pdf.elem_gr for pdf in pdf_cursor]): raise RuntimeError("Projected PDFs not found.") keys = {key for pdf in pdf_cursor for key in pdf.elem_gr} self.elem_gr = {key: np.zeros_like(self.r_space) for key in keys} for pdf in pdf_cursor: for key in pdf.elem_gr: self.elem_gr[key] += np.interp( self.r_space, pdf.r_space, pdf.elem_gr[key] )
[docs] def plot_projected_pdf(self): """Plot the combined PDF.""" from matador.plotting.pdf_plotting import plot_projected_pdf plot_projected_pdf(self)