# -*- coding: utf-8 -*-
"""
Created on Tue Jul 12 20:18:04 2016
Derived from:
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
https://github.com/lammps/lammps/tree/lammps-icms/src/USER-DIFFRACTION
This package contains the commands neeed to calculate x-ray and
electron diffraction intensities based on kinematic diffraction
theory. Detailed discription of the computation can be found in the
following works:
Coleman, S.P., Spearot, D.E., Capolungo, L. (2013) Virtual
diffraction analysis of Ni [010] symmetric tilt grain boundaries,
Modelling and Simulation in Materials Science and Engineering, 21
055020. doi:10.1088/0965-0393/21/5/055020
Coleman, S.P., Sichani, M.M., Spearot, D.E. (2014) A computational
algorithm to produce virtual x-ray and electron diffraction patterns
from atomistic simulations, JOM, 66 (3), 408-416.
doi:10.1007/s11837-013-0829-3
Coleman, S.P., Pamidighantam, S. Van Moer, M., Wang, Y., Koesterke, L.
Spearot D.E (2014) Performance improvement and workflow development
of virtual diffraction calculations, XSEDE14,
doi:10.1145/2616498.2616552
@author: chris sewell
"""
#TODO Averaging over (thermal) phase space
#TODO Analysis of triclinic cells
#TODO Calculation of structure factor coefficients from atom charge & type (rather than pre-defining ionic state)
#TODO fast fourier transform?
#TODO parallelization of fourier summation (http://ipyparallel.readthedocs.io/en/latest/intro.html#getting-started)
import math
import numpy as np
import pandas as pd
from ..shared import get_data_path
from . import data
from . import basic
from .. import plotting
def _set_thetas(min2theta=1.,max2theta=179.):
""" set min and max angles to assess
Properties
----------
min2theta : float
minimum 2 theta range to explore (degrees)
max2theta : float
maximum 2 theta range to explore (degrees)
Returns
-------
min_theta : float
minimum theta range to explore (radians)
max_theta : float
maximum theta range to explore (radians)
"""
# Process angles
min_theta = math.radians(min2theta) / 2.
max_theta = math.radians(max2theta) / 2.
return min_theta, max_theta
def _original_compute_rmesh(sim_abc, wlambda, min_theta, max_theta,
rspace=[1,1,1], periodic=[True,True,True], manual=False):
"""Compute full reciprocal lattice mesh, here for prosperity
Properties
----------
sim_abc : numpy.array((3,3))
a,b,c cell vectors (length units)
wlambda : float
radiation wavelength (length units)
x-rays usually in the range 0.1 to 100 Angstroms
min_theta : float
minimum theta range to explore (radians)
max_theta : float
maximum theta range to explore (radians)
rspace : list of floats
parameters to multiply the spacing of the reciprocal lattice nodes
in the h, k, and l directions respectively
periodic : list of bools
whether periodic boundary in the h, k, and l directions respectively
manual : bool
use manual spacing of reciprocal lattice points based on the values of the rspac parameters
(good for comparing diffraction results from multiple simulations, but small rspace required).
Returns
-------
rmesh : np.array((N,3))
mesh of k points defining reciprocal lattice
"""
# get cell parameters
a,b,c,alpha,beta,gamma = basic.lattparams_bb(sim_abc)
cell_lenghts = (a,b,c)
if not np.allclose([alpha,beta,gamma],[90,90,90]):
raise ValueError("Compute XRD does not work with triclinic structures")
# maximum reciprocal lattice vector |K|,
# calculated from Bragg's law
Kmax = 2 * math.sin(max_theta) / wlambda
# Calculate spacing between reciprocal lattice points
# Using distance based on periodic repeating distance
if manual:
inv_cell_lengths = [1.0,1.0,1.0]
else:
if not np.any(periodic):
raise ValueError("Compute XRD must have at least one periodic boundary unless manual spacing specified")
# calc inverse dimension for periodic directions
inv_cell_lengths = [1./cell_lenghts[i] if periodic[i] else np.nan for i in range(3)]
ave_inv = np.nanmean(inv_cell_lengths)
# Use the average inverse dimensions for non-periodic directions
inv_cell_lengths = [icl if periodic[i] else ave_inv for i,icl in enumerate(inv_cell_lengths)]
# resolution (i.e. spacing) of reciprocal space points
dK = [inv_cell_lengths[i] * rspace[i] for i in range(3)]
# maximum integer value for K points in each dimension
Knmax = [math.ceil(Kmax / dK[i]) for i in range(3)]
# create the full reprocal lattice indices grid
rmesh = np.mgrid[-Knmax[0]:Knmax[0]+1:1,
-Knmax[1]:Knmax[1]+1:1,
-Knmax[1]:Knmax[1]+1:1].reshape(3,-1).T
# resize reciprocal mesh to correct spacing
for i in range(3):
rmesh[:,i] *= dK[i]
return rmesh, imesh
def _compute_rmesh_triclinic(sim_abc, wlambda, min_theta, max_theta,
rspace=[1.,1.,1.],manual=False,periodic=[True,True,True]):
"""Compute full reciprocal lattice mesh
Properties
----------
sim_abc : numpy.array((3,3))
a,b,c cell vectors (length units)
wlambda : float
radiation wavelength (length units)
x-rays usually in the range 0.1 to 100 Angstroms
min_theta : float
minimum theta range to explore (radians)
max_theta : float
maximum theta range to explore (radians)
rspace : list of floats
parameters to multiply the spacing of the reciprocal lattice nodes
in the h, k, and l directions respectively
manual : bool
use manual spacing of reciprocal lattice points based on the values of the rspace parameters
(good for comparing diffraction results from multiple simulations, but small rspace required).
periodic : list of bools
whether periodic boundary in the h, k, and l directions respectively
Returns
-------
rmesh : np.array((N,3))
mesh of k points defining reciprocal lattice
"""
if not np.any(periodic):
raise ValueError('at least one direction must be periodic')
# maximum reciprocal lattice vector |K|,
# calculated from Bragg's law
Kmax = 2 * math.sin(max_theta) / wlambda
# Calculate the rimitive reciprocal lattice vectors
a,b,c = sim_abc
a_recip = np.cross(b,c)/(np.dot(a,np.cross(b,c)))
b_recip = np.cross(c,a)/(np.dot(a,np.cross(b,c)))
c_recip = np.cross(a,b)/(np.dot(a,np.cross(b,c)))
recip_lengths = map(np.linalg.norm, [a_recip,b_recip,c_recip])
# get mean length of periodic directions
mean_length = np.mean(np.array(recip_lengths)[np.array(periodic)])
# set non-periodic directions as the mean length of periodic ones
recip_lengths = [r if p else mean_length for r,p in zip(recip_lengths,periodic)]
# maximum integer value for K points in each dimension
Knmax = [math.ceil(Kmax / (recip_lengths[i]*rspace[i])) for i in range(3)]
# create the full reciprocal lattice indices grid
imesh = np.mgrid[-Knmax[0]:Knmax[0]+1:1,
-Knmax[1]:Knmax[1]+1:1,
-Knmax[1]:Knmax[1]+1:1].reshape(3,-1).T
if not manual:
# calulate reciprocal lattice
rmesh = (np.einsum('j,i->ji',imesh[:,0],a_recip)*rspace[0] +
np.einsum('j,i->ji',imesh[:,1],b_recip)*rspace[1] +
np.einsum('j,i->ji',imesh[:,2],c_recip)*rspace[2])
else:
# calulate reciprocal lattice
rmesh = (imesh[:,0]*rspace[0] +
imesh[:,1]*rspace[1] +
imesh[:,2]*rspace[2])
return rmesh, imesh
def _restrict_rmesh(rmesh, imesh, wlambda, min_theta, max_theta):
"""filter mesh points to only those in Eswald's sphere
(and angular limits)
Parameters
----------
rmesh : np.array((N,3))
mesh of k points defining reciprocal lattice
wlambda : float
radiation wavelength (length units)
x-rays usually in the range 0.1 to 100 Angstroms
min_theta : float
minimum theta range to explore (radians)
max_theta : float
maximum theta range to explore (radians)
Returns
-------
rmesh_sphere : np.array((N,3))
mesh of k points defining reciprocal lattice,
retricted to Eswald's shere (and angular limits)
k_mods : np.array((N,1))
modulus for each k-point
thetas : np.array((N,1))
angles for each k-point (radians)
"""
# calculate the length (squared) of each mesh vector
K_sqr = np.sum(np.square(rmesh),axis=1)
# select only mesh points within the Eswald sphere radius
radius_mask = K_sqr * wlambda**2 <= 2**2 # i.e. (2sin(pi))^2
# calculate the angle of each remaining mesh vector
K = np.sqrt(K_sqr[radius_mask])
theta = np.arcsin(wlambda * np.sqrt(K_sqr[radius_mask]) * 0.5)
# select only mesh points within the angular limits
angle_mask = np.logical_and(theta <= max_theta, theta >= min_theta)
# return remaining mesh points
return rmesh[radius_mask][angle_mask],imesh[radius_mask][angle_mask], K[angle_mask], theta[angle_mask]
[docs]def get_sf_coeffs():
datapath = get_data_path('xray_scattering_factors_coefficients.csv',
module=data)
return pd.read_csv(datapath, index_col=0,comment='#')
def _calc_struct_factors(atoms_df,rmesh_sphere,k_mods):
""" calculate atomic scattering factors, fj,
for each atom at each reciprocal lattice point
Parameters
----------
atoms_df : pandas.DataFrame
a dataframe of info for each atom, including columns; x,y,z,type
rmesh_sphere : np.array((N,3))
mesh of k points defining reciprocal lattice,
retricted to Eswald's shere (and angular limits)
k_mods : np.array((N,1))
modulus for each k-point, only required for calclating Lorentz-polarization factor
Returns
-------
intensities : dict(np.array((N,1)))
structure factor for each k-point (values), for each atom type (keys)
"""
# calculate (|K|/2)^2
K_2_sqr = (0.5*k_mods)**2
# get the structure factor coefficients
sf_coeffs_df = get_sf_coeffs()
struct_factors = {}
for atype in atoms_df.type.unique():
sfs = sf_coeffs_df.loc[atype]
struct_factors[atype] = 0
struct_factors[atype] += sfs.A1 * np.exp(-sfs.B1*K_2_sqr)
struct_factors[atype] += sfs.A2 * np.exp(-sfs.B2*K_2_sqr)
struct_factors[atype] += sfs.A3 * np.exp(-sfs.B3*K_2_sqr)
struct_factors[atype] += sfs.A4 * np.exp(-sfs.B4*K_2_sqr)
struct_factors[atype] += sfs.C
return struct_factors
def _calc_intensities(atoms_df, rmesh_sphere, wlambda, struct_factors,
thetas=None,k_mods=None,use_Lp=True):
""" calculate diffraction intensities for each atom at each reciprocal lattice point
Parameters
----------
atoms_df : pandas.DataFrame
a dataframe of info for each atom, including columns; x,y,z,type
rmesh_sphere : np.array((N,3))
mesh of k points defining reciprocal lattice,
retricted to Eswald's shere (and angular limits)
wlambda : float
radiation wavelength (length units)
x-rays usually in the range 0.1 to 100 Angstroms
k_mods : np.array((N,1))
modulus for each k-point, only required for calclating Lorentz-polarization factor
thetas : np.array((N,1))
angles for each k-point (radians), only required for calclating Lorentz-polarization factor
use_Lp : bool
switch to apply Lorentz-polarization factor
Returns
-------
2thetas : np.array((N,1))
2-theta (degrees) for each k-point
intensities : np.array((N,1))
intensity for each k-point
hkl : np.array((N,3))
h,k,l for each k-point,
N.B.: if analysing a crystal then they will need to be divided by the number of unit cells (in each direction)
"""
# compute F(K)
F = np.zeros(rmesh_sphere.shape[0]) + 0*1j
for xyz,atype in zip(atoms_df[['x','y','z']].values, atoms_df.type):
inner_dot = 2 * np.pi * np.dot(xyz,rmesh_sphere.T)
F += struct_factors[atype] * (np.cos(inner_dot) + 1j*np.sin(inner_dot))
# compute Lp(theta)
if use_Lp:
sin_thetas = 0.5*k_mods*wlambda
Lp = (1+np.cos(2*thetas)**2)/(np.cos(thetas)*sin_thetas**2)
else:
Lp = 1.
# calculate intensities
return Lp*F*np.conjugate(F)/float(atoms_df.shape[0])
[docs]def xrd_compute(atoms_df, meta_data,wlambda, min2theta=1.,max2theta=179., lp=True,
rspace=[1,1,1], manual=False,periodic=[True,True,True]):
r"""Compute predicted x-ray diffraction intensities for a given wavelength
Properties
----------
atoms_df : pandas.DataFrame
a dataframe of info for each atom, including columns; x,y,z,type
meta_data : pandas.Series
data of a,b,c crystal vectors (as tuples, e.g. meta_data.a = (0,0,1))
wlambda : float
radiation wavelength (length units)
typical values are Cu Ka = 1.54, Mo Ka = 0.71 Angstroms
min2theta : float
minimum 2 theta range to explore (degrees)
max2theta : float
maximum 2 theta range to explore (degrees)
lp : bool
switch to apply Lorentz-polarization factor
use_triclinic : bool
use_triclinic
rspace : list of floats
parameters to adjust the spacing of the reciprocal lattice nodes
in the h, k, and l directions respectively
manual : bool
use manual spacing of reciprocal lattice points based on the values of the c parameters
(good for comparing diffraction results from multiple simulations, but small c required).
periodic : list of bools
whether periodic boundary in the h, k, and l directions respectively
Returns
-------
df : pandas.DataFrame
theta (in degrees), I, h, k and l for each k-point
Notes
-----
This is an implementation of the virtual x-ray diffraction pattern algorithm
by Coleman *et al*. [ref1]_
The algorithm proceeds in the following manner:
1. Define a crystal structure by position (x,y,z) and atom/ion type.
2. Define the x-ray wavelength to use
3. Compute the full reciprocal lattice mesh
4. Filter reciprocal lattice points by those in the Eswald's sphere
5. Compute the structure factor at each reciprocal lattice point, for each atom type
6. Compute the x-ray diffraction intensity at each reciprocal lattice point
7. Group and sum intensities by angle
reciprocal points of the lattice are computed such that:
.. math::
\mathbf{K} = {m_{1}}\cdot \mathbf{b}_{1}+{m_{2}}\cdot \mathbf{b}_{2}+{m_{3}}\cdot \mathbf{b}_{3}
where,
.. math::
\begin{aligned}
\mathbf {b_{1}} &= {\frac {\mathbf {a_{2}} \times \mathbf {a_{3}} }{\mathbf {a_{1}} \cdot (\mathbf {a_{2}} \times \mathbf {a_{3}} )}}\\
\mathbf {b_{2}} &= {\frac {\mathbf {a_{3}} \times \mathbf {a_{1}} }{\mathbf {a_{2}} \cdot (\mathbf {a_{3}} \times \mathbf {a_{1}} )}}\\
\mathbf {b_{3}} &= {\frac {\mathbf {a_{1}} \times \mathbf {a_{2}} }{\mathbf {a_{3}} \cdot (\mathbf {a_{1}} \times \mathbf {a_{2}} )}}
\end{aligned}
The reciprocal k-point modulii of the x-ray is calculated from Bragg's law:
.. math::
\left| {\mathbf{K}_{\lambda}} \right| = \frac{1}{{d_{\text{hkl}} }} = \frac{2\sin \left( \theta \right)}{\lambda }
This is used to construct an Eswald's sphere,
and only reciprocal lattice ponts within are retained, as illustrated:
.. image:: ../images/xrd_mesh.jpg
The atomic scattering factors, fj, accounts for the reduction in
diffraction intensity due to Compton scattering, with coefficients based on
the electron density around the atomic nuclei.
.. math::
f_j \left(\frac{\sin \theta}{\lambda}\right)
= \left[ \sum\limits^4_i a_i \exp \left( -b_i \frac{\sin^2 \theta}{\lambda^2} \right)\right] + c
= \left[ \sum\limits^4_i a_i \exp \left( -b_i \left(\frac{\left| {\mathbf{K}} \right|}{2}\right)^2 \right)\right] + c
The relative diffraction intensity from x-rays is computed at each
reciprocal lattice point through:
.. math::
I_x(\mathbf{K}) = Lp(\theta) \frac{F(\mathbf{K})F^*(\mathbf{K})}{N}
such that:
.. math::
F ({\mathbf{K}} )=
\sum\limits_{j = 1}^{N} {f_{j}.e^{\left( {2\pi i \, {\mathbf{K}} \cdot {\mathbf{r}}_{j} } \right)}}
= \sum\limits_{j = 1}^{N} {f_j.\left[ \cos \left( 2\pi \mathbf{K} \cdot \mathbf{r}_j \right) + i \sin \left( 2\pi \mathbf{K} \cdot \mathbf{r}_j \right) \right]}
and the Lorentz-polarization factor is:
.. math::
Lp(\theta) = \frac{1+\cos^2 (2\theta)}{\cos(\theta)\sin^2(\theta)}
References
----------
.. [ref1] 1.Coleman, S. P., Sichani, M. M. & Spearot, D. E.
A Computational Algorithm to Produce Virtual X-ray and Electron Diffraction
Patterns from Atomistic Simulations. JOM 66, 408–416 (2014).
"""
sim_abc = np.asarray([meta_data.a,meta_data.b,meta_data.c])
min_theta, max_theta = _set_thetas(min2theta,max2theta)
rmesh, imesh = _compute_rmesh_triclinic(sim_abc,wlambda,min_theta, max_theta,rspace, manual, periodic)
rmesh_sphere,imesh_sphere, k_mods, thetas = _restrict_rmesh(rmesh,imesh,wlambda,min_theta, max_theta)
struct_factors = _calc_struct_factors(atoms_df,rmesh_sphere,k_mods)
I = _calc_intensities(atoms_df,rmesh_sphere,wlambda,struct_factors,thetas,k_mods,use_Lp=lp)
return (pd.DataFrame({'theta':np.real(np.degrees(2*thetas)),'I':np.real(I),
'h':imesh_sphere[:,0].astype(int),'k':imesh_sphere[:,1].astype(int),'l':imesh_sphere[:,2].astype(int)}))
##TODO add other symmetry-equivalent reflections
def _stack_hkl(df, sym='none'):
hkl = df[['h','k','l']].values
if sym=='mmm':
hkl = np.abs(hkl)
elif sym=='m-3m':
hkl = np.sort(np.abs(hkl))[:,::-1]
return np.vstack({tuple(row) for row in hkl})
[docs]def xrd_group_i(df,tstep=None,sym_equiv_hkl='none'):
""" group xrd intensities by theta
Parameters
----------
df : pandas.DataFrame
containing columns; ['theta','I'] and optional ['h','k','l']
tstep: None or float
if not None, group thetas within ranges i to i+tstep
sym_equiv_hkl : str
group hkl by symmetry-equivalent refelections; ['none','mmm','m-3m']
Returns
-------
df : pandas.DataFrame
Notes
-----
if grouping by theta step, then the theta value for each group will be the intensity weighted average
Crystal System | Laue Class | Symmetry-Equivalent Reflections | Multiplicity
-------------- | ---------- | ------------------------------- | ------------
Triclinic | -1 | hkl ≡ -h-k-l | 2
Monoclinic | 2/m | hkl ≡ -hk-l ≡ -h-k-l ≡ h-kl | 4
Orthorhombic | mmm | hkl ≡ h-k-l ≡ -hk-l ≡ -h-kl ≡ -h-k-l ≡ -hkl ≡ h-kl ≡ hk-l | 8
Tetragonal | 4/m | hkl ≡ -khl ≡ -h-kl ≡ k-hl ≡ -h-k-l ≡ k-h-l ≡ hk-l ≡ -kh-l | 8
| 4/mmm | hkl ≡ -khl ≡ -h-kl ≡ k-hl ≡ -h-k-l ≡ k-h-l ≡ hk-l ≡ -kh-l |
| | ≡ khl ≡ -hkl ≡ -k-hl ≡ h-kl ≡ -k-h-l ≡ h-k-l ≡ kh-l ≡ -hk-l | 16
Cubic | m-3 | hkl ≡ -hkl ≡ h-kl ≡ hk-l ≡ -h-k-l ≡ h-k-l ≡ -hk-l ≡ -h-kl |
| | ≡ klh ≡ -klh ≡ k-lh ≡ kl-h ≡ -k-l-h ≡ k-l-h ≡ -kl-h ≡ -k-lh |
| | ≡ lhk ≡ -lhk ≡ l-hk ≡ lh-k ≡ -l-h-k ≡ l-h-k ≡ -lh-k ≡ -l-hk | 24
| m-3m | hkl ≡ -hkl ≡ h-kl ≡ hk-l ≡ -h-k-l ≡ h-k-l ≡ -hk-l ≡ -h-kl |
| | ≡ klh ≡ -klh ≡ k-lh ≡ kl-h ≡ -k-l-h ≡ k-l-h ≡ -kl-h ≡ -k-lh |
| | ≡ lhk ≡ -lhk ≡ l-hk ≡ lh-k ≡ -l-h-k ≡ l-h-k ≡ -lh-k ≡ -l-hk |
| | ≡ khl ≡ -khl ≡ k-hl ≡ kh-l ≡ -k-h-l ≡ k-h-l ≡ -kh-l ≡ -k-hl |
| | ≡ lkh ≡ -lkh ≡ l-kh ≡ lk-h ≡ -l-k-h ≡ l-k-h ≡ -lk-h ≡ -l-kh |
| | ≡ hlk ≡ -hlk ≡ h-lk ≡ hl-k ≡ -h-l-k ≡ h-l-k ≡ -hl-k ≡ -h-lk | 48
"""
if tstep is None:
g=df.groupby('theta')
g_df = pd.DataFrame({'I':g.I.sum(),'multiplicity':g.I.count()}).reset_index()
else:
df1 = df.copy()
df1['theta_ind'] = np.digitize(df.theta,np.arange(0,180,tstep))
df1['wt_theta'] = df.theta * df.I
g = df1.groupby('theta_ind')
g_df = pd.DataFrame({'theta':g.wt_theta.sum()/g.I.sum(),'I':g.I.sum(),'multiplicity':g.I.count()})
g_df['I_norm'] = g_df.I/g_df.I.max()
if set(['h','k','l']).issubset(df.columns):
g_df['hkl'] = g.apply(_stack_hkl,sym_equiv_hkl).tolist()
return g_df
[docs]def xrd_plot(df, icol='I_norm', imin=0.01, barwidth=1.,
hkl_num=0, hkl_trange=[0.,180.], incl_multi=False,
label_inline=True, label_trange=[0.,180.],
ax=None, **kwargs):
""" create plot of xrd pattern
Properties
----------
df : pd.DataFrame
containing columns ['theta',icol] and optional ['hkl','multiplicity']
icol : str
column name for intensity data
imin : float
minimum intensity to display
barwidth : float or None
if None then the barwidth will be the bin width
hkl_num : int
number of k-point values to label
hkl_trange : [float,float]
theta range within which to label peaks with k-point values
label_inline : bool
place k-point labels inline or at top of plot
label_trange : [float,float]
if not inline, theta range within which to fit k-point labels
incl_multi : bool
add multiplicity to k-point label
kwargs : optional
additional arguments for bar plot (e.g. label, color, alpha)
Returns
-------
plot : ipymd.plotting.Plotting
a plot object
"""
if ax is None:
plot = plotting.Plotter()
ax = plot.axes
else:
plot=None
df1 = df[df[icol]>=imin]
ax.bar(df1.theta-0.5*barwidth,df1[icol],barwidth,**kwargs)
ax.set_xlabel(r'Scatteting Angle ($2 \theta$)')
ax.set_ylabel('Relative Intensity')
ax.set_yticklabels([])
ax.grid(True)
if 'label' in kwargs.keys():
ax.legend(framealpha=0.5)
if hkl_num > 0 and 'hkl' in df.columns:
l,u = hkl_trange
if label_inline:
df2 = df1[(df1.theta>=l)&(df1.theta<=u)].sort_values(icol,ascending=True).iloc[-hkl_num:]
for idx, row in df2.iterrows():
hkl_string = '\n'.join([str(hkl) for hkl in row.hkl])
if incl_multi and 'multiplicity' in df.columns:
hkl_string += '\nm='+str(row.multiplicity)
hkl_string = hkl_string.replace('. ',' ')
hkl_string = hkl_string.replace('[ ','[')
hkl_string = hkl_string.replace('.]',']')
hkl_string = hkl_string.replace(' ',' ')
bbox=dict(facecolor='white',edgecolor='none',alpha=0.8,pad=0)
xydata = (row.theta+barwidth*0.5,row[icol])
ax.annotate(hkl_string, xy=xydata, bbox=bbox)
else:
df2 = df1[(df1.theta>=l)&(df1.theta<=u)].sort_values(icol,ascending=True).iloc[-hkl_num:]
df2.sort_values('theta',ascending=True,inplace=True)
spacing = np.linspace(label_trange[0],label_trange[1],hkl_num,endpoint=False).tolist()
for idx, row in df2.iterrows():
hkl_string = '\n'.join([str(hkl) for hkl in row.hkl])
if incl_multi and 'multiplicity' in df.columns:
hkl_string += '\nm='+str(row.multiplicity)
hkl_string = hkl_string.replace('. ',' ')
hkl_string = hkl_string.replace('[ ','[')
hkl_string = hkl_string.replace('.]',']')
hkl_string = hkl_string.replace(' ',' ')
bbox=dict(facecolor='white',edgecolor='none',alpha=0.8,pad=0)
xydata = (row.theta,row[icol])
xytext = (spacing.pop(0),df[icol].max()*1.2)
if xytext[0]<xydata[0]:
cstyle="arc3,rad=.1"
else:
cstyle="arc3,rad=-.1"
ax.annotate(hkl_string, xy=xydata,xytext=xytext, bbox=bbox,
arrowprops=dict(arrowstyle="->",
linestyle='solid',connectionstyle=cstyle))
return plot
##TODO identification and classification of peaks