If all pair-wise distances are cast into a histogram, the scattering equation can be written: $$ I \left( Q \right) = \left| f(Q) \right|^2 \left[ N + 2 \sum_{i=1}^{N_{bins}} n_i {sin \left(Q r_{i} \right) \over {Q r_{i}}} \right] $$ where $n_i$ is the number of atomic pairs within the bin corresponding to a distance $r_i$.
Load external libraries.
import math
import numpy as np
import numexpr as ne
import numba as nb
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist
import os
n_cpu = os.cpu_count()
try:
ref_intensity = np.loadtxt("ref_intensity.txt")
except :
print("No reference intensity file found. An error will occur when evaluating the relative error in the computed intensities.")
Now we load the atomic coordinates and the coefficients required for the calculation of the atomic scattering factor. The first 11 elements of the coeff array are the coefficients given by Waasmaier & Kirfel for each element of the periodic table. The last element is the absorption/anomalous dispersion coefficient as given in the Brennan-Cowan database.
# Atomic coordinates within cluster
coords = np.loadtxt("Au.xyz")
# Atomic scattering factor coefficients
coeff = np.array([16.777390,19.317156,32.979683,5.595453,10.576854,-6.279078,0.122737,
8.621570,1.256902,38.008820,0.000601, -4.014837862516359+7.717066375953113j])
# Define 2theta computation range
tth = np.arange(10, 180, 0.2)*np.pi/180
# X-ray wavelength
wl = 1.5406
# Step (in Angstroms) to compute the distance histogram
step=0.0005
Global variables computed from the above
N, dim = np.shape(coords) #Extract the number of atoms and the number of dimensions
Q = 4*np.pi*np.sin(tth/2) / wl #Compute the scatttering vector range
intensity_pt = np.loadtxt("ref_intensity.txt")
%load_ext Cython
%%cython --compile-args=-fopenmp --link-args=-fopenmp
import numpy as np
from cython.parallel import prange
import cython
cimport numpy as cnp
from libc.math cimport sqrt
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def r_ij_cython(double[:,::1] coords):
cdef:
double[::1] r
double tmp
cnp.int64_t i,j,k,l
cnp.int64_t N = coords.shape[0]
cnp.int64_t dim = coords.shape[1]
r = np.zeros(int((N*N-N)/2), dtype=np.float64)
for i in prange(N, nogil=True):
for j in range(i+1,N):
l = int(i * (N - 1) - i * (i + 1) / 2 + j - 1)
tmp = 0.0
for k in range(dim):
tmp = tmp + (coords[i,k]-coords[j,k])**2
r[l] = sqrt(tmp)
return r
%timeit r_ij_cython(coords)
r_ij_ct = np.asarray(r_ij_cython(coords))
bins=np.arange(r_ij_ct.min(), r_ij_ct.max()+step, step)
Numpy has a built-in function to evaluate a histogram from a given data set. The result of this evaluation will be used to check the validity of our calculations.
%timeit np.histogram(r_ij_ct, bins)[0]
Nr_sp = np.histogram(r_ij_ct, bins)[0]
Compute the histogram from the pairwise distance distribution: serial algorithm
%%cython --compile-args=-fopenmp --link-args=-fopenmp
import numpy as np
import cython
cimport numpy as cnp
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def hist_cython(double[::1] data, double[::1] bins):
cdef:
double[::1] hist
double max_bins, min_bins
cnp.int64_t i, data_size, hist_size, index
data_size = data.size
hist_size = bins.size
max_bins = max(bins)
min_bins = min(bins)
hist = np.zeros(hist_size-1)
step = (max_bins-min_bins)/(hist_size-1)
for i in range(data_size):
index = int((data[i]-min_bins)/step)
hist[index] = hist[index] + 1
return hist
Nr_ct = hist_cython(r_ij_ct, bins)
%timeit hist_cython(r_ij_ct, bins)
print("Is the result correct:",np.allclose(Nr_ct, Nr_sp))
Parallel algorithm: in order to avoid the different threads to write at the same moment in the same memory block, each thread owns a copy of the histogram. All copies are summed up in the end.
%%cython --compile-args=-fopenmp --link-args=-fopenmp
import numpy as np
from cython.parallel import prange, threadid
import cython
cimport numpy as cnp
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def hist_cython_p(double[::1] data, double[::1] bins, int threads):
cdef:
double[::1] hist
double[:,::1] hist_ar
double max_bins, min_bins
cnp.int64_t i, j, data_size, hist_size, index
data_size = data.size
hist_size = bins.size-1
max_bins = max(bins)
min_bins = min(bins)
hist = np.zeros(hist_size)
hist_ar = np.zeros((threads,hist_size))
step = (max_bins-min_bins)/(hist_size)
for i in prange(data_size, nogil=True, num_threads = threads):
index = int((data[i]-min_bins)/step)
hist_ar[threadid(), index] = hist_ar[threadid(),index] + 1
for i in range(threads):
for j in range(hist_size):
hist[j] = hist[j] + hist_ar[i,j]
return hist
Nr_ctp = hist_cython_p(r_ij_ct, bins, n_cpu)
%timeit hist_cython_p(r_ij_ct, bins, n_cpu)
print("Is the result correct:",np.allclose(Nr_ctp, Nr_sp))
First compute the atomic scattering factor
def asf(th,wl,coeff):
n = 0
sum = coeff[5]
while n <= 4:
sum = sum + coeff[n] * np.exp ( -coeff[n+6] * ( np.sin(th) / wl )**2)
n = n+1
return sum + coeff[-1]
f_at = asf(tth/2, wl, coeff)
Finally, compute the intensity using the histogram.
%%cython --compile-args=-fopenmp --link-args=-fopenmp
import numpy as np
import math
from cython.parallel import prange
import cython
cdef extern from "math.h" nogil:
double sin(double x)
cdef extern from "complex.h" nogil:
double cabs(double complex)
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.cdivision(True)
def Debye_cython_hist(double[::1] Q, double[::1] r, double[::1] w,
long N, double complex[::1] f_at):
cdef:
double tmp
long i_r, i_Q, size_r, size_Q
double[::1] res
size_r = r.size
size_Q = Q.size
res = np.zeros((size_Q))
for i_Q in prange(size_Q, nogil=True):
tmp = 0.0
for i_r in range(size_r):
tmp = tmp + w[i_r]*sin(Q[i_Q]*r[i_r])/(Q[i_Q]*r[i_r])
res[i_Q] = (N + 2*tmp)*cabs(f_at[i_Q])**2
return res
bins=bins[:-1:]
%timeit Debye_cython_hist(Q, bins, Nr_ctp, N,f_at)
intensity_ct_hist = Debye_cython_hist(Q, bins,Nr_ctp, N,f_at)
print("Error:", abs(intensity_ct_hist-ref_intensity).max()/ref_intensity.max())
%matplotlib notebook
plt.semilogy(tth*180/np.pi, intensity_ct_hist)
plt.show()