wmaee.scopes.structural_analysis

 1import numpy as np
 2import pandas as pd
 3from itertools import product
 4
 5def rdf(struct, rmin=1, rmax=10, nbins=90):
 6    """
 7    Calculate radial distribution function (RDF) for a given structure
 8    :param struct: (pymatgen.Structure) the structure
 9    :param rmin: (float) Rmin of the radial mesh (default: 1 (Angstroem))
10    :param rmax: (float) Rmax of the radial mesh (default: 10 (Angstroem))
11    :param nbins: (int) number of bins (points) along the radial mesh (default: 90)
12    :return: (r:numpy.array, rdf:numpy.array) radial mesh, RDF
13    """   
14
15    # define radial mesh for RDF
16    bin_edges = np.linspace(rmin, rmax, nbins+1)
17    r = [0.5*(bin_edges[i+1]+bin_edges[i]) for i in range(len(bin_edges)-1)]
18    dr = [(bin_edges[i+1]+bin_edges[i]) for i in range(len(bin_edges)-1)]
19
20    sc = [int(np.ceil(rmax/length)) for length in struct.lattice.abc]
21    struct.make_supercell(sc)
22    dists = struct.distance_matrix.reshape(-1)
23    hist, _ = np.histogram(dists, bins=bin_edges)
24    return np.array(r), np.array([hist[i]/(4*np.pi*r[i]**2*dr[i]) for i in range(len(hist))])
25
26def rdf_partial(struct, rmin=1, rmax=10, nbins=90):
27    """
28    Calculate radial distribution function (RDF) including all partial RDFs for a given structure
29    :param struct: (pymatgen.Structure) the structure
30    :param rmin: (float) Rmin of the radial mesh (default: 1 (Angstroem))
31    :param rmax: (float) Rmax of the radial mesh (default: 10 (Angstroem))
32    :param nbins: (int) number of bins (points) along the radial mesh (default: 90)
33    :return: (pandas.DataFrame) [R, RDF_tot, RDF_X-X, RDF_X-Y, ...]
34    """
35
36    bin_width = (rmax-rmin)/nbins
37
38    # get elements for partial RDF
39    species = sorted(list(struct.symbol_set))
40    cols = tuple(species[i]+'-'+species[j] for i in range(len(species)) for j in range(i, len(species)))
41
42    # create empty RDF
43    RDF = pd.DataFrame(columns=('r', 'tot') + cols)
44    row = {col: 0 for col in cols}
45    row['tot'] = 0
46    for i in range(nbins):
47        row['r'] = rmin+(i+0.5)*bin_width
48        RDF = RDF.append(row, ignore_index=True)
49        
50    # get SC sizes necessary for covering the desired range
51    SC = [int(np.ceil(rmax/length)) for length in struct.lattice.abc]
52    shifts = product(np.arange(-SC[0], SC[0]+1), np.arange(-SC[1], SC[1]+1), np.arange(-SC[2], SC[2]+1))
53    # list of translational vertors to be applied during RDF evaluation to go beyond
54    # period boundaries
55    shifts = [np.dot(shift, struct.lattice.matrix) for shift in shifts]
56
57    # for s in struct.sites:
58    for s in struct.sites:
59        for t in struct.sites:
60            # which bond
61            col = '-'.join(sorted([s.species_string, t.species_string]))
62            # get all SC distances
63            dist = [np.linalg.norm(s.coords-t.coords+shift) for shift in shifts]
64            iR = np.floor((np.array(dist)-rmin)/bin_width).astype(int)
65            for i in iR:
66                if i>=0 and i<nbins:
67                    RDF['tot'][i] += 1
68                    RDF[col][i] += 1
69
70    # normalize per area
71    RDF['tot'] = RDF['tot']/(4*np.pi*bin_width*RDF['r']**2)
72    for col in cols:
73        RDF[col] = RDF[col]/(4*np.pi*bin_width*RDF['r']**2)
74    
75    return RDF