# flake8: noqa
"""Definition of the XrDebye class.
This module defines the XrDebye class for calculation
of X-ray scattering properties from atomic cluster
using Debye formula.
Also contains routine for calculation of atomic form factors and
X-ray wavelength dict.
"""
from math import acos, cos, exp, pi, sin, sqrt
import numpy as np
from ase.data import atomic_numbers
# Table (1) of
# D. WAASMAIER AND A. KIRFEL, Acta Cryst. (1995). A51, 416-431
waasmaier = {
# a1 b1 a2 b2 a3 b3 a4 b4 a5 b5 c
'C': [2.657506, 14.780758, 1.078079, 0.776775, 1.490909, 42.086843, -4.241070, -0.000294, 0.713791, 0.239535, 4.297983],
'N': [11.893780, 0.000158, 3.277479, 10.232723, 1.858092, 30.344690, 0.858927, 0.656065, 0.912985, 0.217287, -11.804902],
'O': [2.960427, 14.182259, 2.5088111, 5.936858, 0.637053, 0.112726, 0.722838, 34.958481, 1.142756, 0.390240, 0.027014],
'P': [1.950541, 0.908139, 4.146930, 27.044953, 1.494560, 0.071280, 1.522042, 67.520190, 5.729711, 1.981173, 0.155233],
'S': [6.372157, 1.514347, 5.154568, 22.092528, 1.473732, 0.061373, 1.635073, 55.445176, 1.209372, 0.646925, 0.154722],
'Cl': [1.446071, 0.052357, 6.870609, 1.193165, 6.151801, 18.343416, 1.750347, 46.398394, 0.634168, 0.401005, 0.146773],
'Ni': [13.521865, 4.077277, 6.947285, 0.286763, 3.866028, 14.622634, 2.135900, 71.966078, 4.284731, 0.004437, -2.762697],
'Cu': [14.014192, 3.738280, 4.784577, 0.003744, 5.056806, 13.034982, 1.457971, 72.554793, 6.932996, 0.265666, -3.774477],
'Pd': [6.121511, 0.062549, 4.784063, 0.784031, 16.631683, 8.751391, 4.318258, 34.489983, 13.246773, 0.784031, 0.883099],
'Ag': [6.073874, 0.055333, 17.155437, 7.896512, 4.173344, 28.443739, 0.852238, 110.376108, 17.988685, 0.716809, 0.756603],
'Pt': [31.273891, 1.316992, 18.445441, 8.797154, 17.063745, 0.124741, 5.555933, 40.177994, 1.575270, 1.316997, 4.050394],
'Au': [16.777389, 0.122737, 19.317156, 8.621570, 32.979682, 1.256902, 5.595453, 38.008821, 10.576854, 0.000601, -6.279078],
}
wavelengths = {
'CuKa1': 1.5405981,
'CuKa2': 1.54443,
'CuKb1': 1.39225,
'WLa1': 1.47642,
'WLa2': 1.48748
}
[docs]
class XrDebye:
"""
Class for calculation of XRD or SAXS patterns.
"""
def __init__(self, atoms, wavelength, damping=0.04,
method='Iwasa', alpha=1.01, warn=True):
"""
Initilize the calculation of X-ray diffraction patterns
Parameters:
atoms: ase.Atoms
atoms object for which calculation will be performed.
wavelength: float, Angstrom
X-ray wavelength in Angstrom. Used for XRD and to setup dumpings.
damping : float, Angstrom**2
thermal damping factor parameter (B-factor).
method: {'Iwasa'}
method of calculation (damping and atomic factors affected).
If set to 'Iwasa' than angular damping and q-dependence of
atomic factors are used.
For any other string there will be only thermal damping
and constant atomic factors (`f_a(q) = Z_a`).
alpha: float
parameter for angular damping of scattering intensity.
Close to 1.0 for unplorized beam.
warn: boolean
flag to show warning if atomic factor can't be calculated
"""
self.wavelength = wavelength
self.damping = damping
self.mode = ''
self.method = method
self.alpha = alpha
self.warn = warn
self.twotheta_list = []
self.q_list = []
self.intensity_list = []
self.atoms = atoms
# TODO: setup atomic form factors if method != 'Iwasa'
[docs]
def set_damping(self, damping):
""" set B-factor for thermal damping """
self.damping = damping
[docs]
def get(self, s):
r"""Get the powder x-ray (XRD) scattering intensity
using the Debye-Formula at single point.
Parameters:
s: float, in inverse Angstrom
scattering vector value (`s = q / 2\pi`).
Returns:
Intensity at given scattering vector `s`.
"""
pre = exp(-self.damping * s**2 / 2)
if self.method == 'Iwasa':
sinth = self.wavelength * s / 2.
positive = 1. - sinth**2
if positive < 0:
positive = 0
costh = sqrt(positive)
cos2th = cos(2. * acos(costh))
pre *= costh / (1. + self.alpha * cos2th**2)
f = {}
def atomic(symbol):
"""
get atomic factor, using cache.
"""
if symbol not in f:
if self.method == 'Iwasa':
f[symbol] = self.get_waasmaier(symbol, s)
else:
f[symbol] = atomic_numbers[symbol]
return f[symbol]
I = 0.
fa = [] # atomic factors list
for a in self.atoms:
fa.append(atomic(a.symbol))
pos = self.atoms.get_positions() # positions of atoms
fa = np.array(fa) # atomic factors array
for i in range(len(self.atoms)):
vr = pos - pos[i]
I += np.sum(fa[i] * fa * np.sinc(2 * s *
np.sqrt(np.sum(vr * vr, axis=1))))
return pre * I
[docs]
def get_waasmaier(self, symbol, s):
r"""Scattering factor for free atoms.
Parameters:
symbol: string
atom element symbol.
s: float, in inverse Angstrom
scattering vector value (`s = q / 2\pi`).
Returns:
Intensity at given scattering vector `s`.
Note:
for hydrogen will be returned zero value."""
if symbol == 'H':
# XXXX implement analytical H
return 0
elif symbol in waasmaier:
abc = waasmaier[symbol]
f = abc[10]
s2 = s * s
for i in range(5):
f += abc[2 * i] * exp(-abc[2 * i + 1] * s2)
return f
if self.warn:
print('<xrdebye::get_atomic> Element', symbol, 'not available')
return 0
[docs]
def calc_pattern(self, x=None, mode='XRD', verbose=False):
r"""
Calculate X-ray diffraction pattern or
small angle X-ray scattering pattern.
Parameters:
x: float array
points where intensity will be calculated.
XRD - 2theta values, in degrees;
SAXS - q values in 1/A
(`q = 2 \pi \cdot s = 4 \pi \sin( \theta) / \lambda`).
If ``x`` is ``None`` then default values will be used.
mode: {'XRD', 'SAXS'}
the mode of calculation: X-ray diffraction (XRD) or
small-angle scattering (SAXS).
Returns:
list of intensities calculated for values given in ``x``.
"""
self.mode = mode.upper()
assert (mode in ['XRD', 'SAXS'])
result = []
if mode == 'XRD':
if x is None:
self.twotheta_list = np.linspace(15, 55, 100)
else:
self.twotheta_list = x
self.q_list = []
if verbose:
print('#2theta\tIntensity')
for twotheta in self.twotheta_list:
s = 2 * sin(twotheta * pi / 180 / 2.0) / self.wavelength
result.append(self.get(s))
if verbose:
print(f'{twotheta:.3f}\t{result[-1]:f}')
elif mode == 'SAXS':
if x is None:
self.twotheta_list = np.logspace(-3, -0.3, 100)
else:
self.q_list = x
self.twotheta_list = []
if verbose:
print('#q\tIntensity')
for q in self.q_list:
s = q / (2 * pi)
result.append(self.get(s))
if verbose:
print(f'{q:.4f}\t{result[-1]:f}')
self.intensity_list = np.array(result)
return self.intensity_list
[docs]
def write_pattern(self, filename):
""" Save calculated data to file specified by ``filename`` string."""
with open(filename, 'w') as fd:
self._write_pattern(fd)
def _write_pattern(self, fd):
fd.write('# Wavelength = %f\n' % self.wavelength)
if self.mode == 'XRD':
x, y = self.twotheta_list, self.intensity_list
fd.write('# 2theta \t Intesity\n')
elif self.mode == 'SAXS':
x, y = self.q_list, self.intensity_list
fd.write('# q(1/A)\tIntesity\n')
else:
raise Exception('No data available, call calc_pattern() first.')
for i in range(len(x)):
fd.write(f' {x[i]:f}\t{y[i]:f}\n')
[docs]
def plot_pattern(self, filename=None, show=False, ax=None):
""" Plot XRD or SAXS depending on filled data
Uses Matplotlib to plot pattern. Use *show=True* to
show the figure and *filename='abc.png'* or
*filename='abc.eps'* to save the figure to a file.
Returns:
``matplotlib.axes.Axes`` object."""
import matplotlib.pyplot as plt
if ax is None:
plt.clf() # clear figure
ax = plt.gca()
if self.mode == 'XRD':
x, y = np.array(self.twotheta_list), np.array(self.intensity_list)
ax.plot(x, y / np.max(y), '.-')
ax.set_xlabel('2$\\theta$')
ax.set_ylabel('Intensity')
elif self.mode == 'SAXS':
x, y = np.array(self.q_list), np.array(self.intensity_list)
ax.loglog(x, y / np.max(y), '.-')
ax.set_xlabel('q, 1/Angstr.')
ax.set_ylabel('Intensity')
else:
raise Exception('No data available, call calc_pattern() first')
if show:
plt.show()
if filename is not None:
fig = ax.get_figure()
fig.savefig(filename)
return ax