如何快速找到布拉格反射与矮胖？
对于x射线衍射，需要找到所谓的劳厄方程的解。
G_hkl - k_in + |k_in|*(sin(theta) cos(phi)，sin(theta) sin(phi)，cos(theta))=0
当G_hkl是给定的三维向量时，k_in可以选择为(0,0,1)，θ和phi是满足方程的自由参数。在一个典型的实验中，G_hkl绕x轴旋转，旋转的每一步都需要找到方程的解。在给定的旋转条件下，方程不能有多个解。
我编写了这个python脚本来找到这些解决方案，但是对于我的应用程序来说，它还不够快。
import numpy as np
import time
# Just initialization of variables. This is fast enough
res_list = []
N_rot = 100
Ghkl = np.array([0,0.7,0.7])
Ghkl_norm = np.linalg.norm(Ghkl)
kin = np.array([0,0,1])
kin_norm = np.linalg.norm(kin)
alpha_step = 2*np.pi/(N_rot-1)
rot_mat = np.array([[1,0,0],[0,np.cos(alpha_step),-np.sin(alpha_step)],[0,np.sin(alpha_step),np.cos(alpha_step)]])
# You can deduce theta from the norm of the vector equation
theta = np.arccos(1-(Ghkl_norm**2/(2*kin_norm**2)))
sint = np.sin(theta)
cost = np.cos(theta)
# This leaves only phi as paramter to find
# I find phi by introducing a finite test vector
# and testing whether the norm of the vector equation is close
# to zero for any of those test phis
phi_test = np.linspace(0,2*np.pi,200)
kout = kin_norm * np.array([sint * np.cos(phi_test), sint * np.sin(phi_test), cost + 0*phi_test]).T
##############
start_time = time.time()
for j in range(100): # just to make it longer to measure the time
res_list = []
for i in range(N_rot): # This loop is too slow
# Here the norm of the vector equation is calculated for all phi_test
norm_vec = np.linalg.norm(Ghkl[np.newaxis, :] - kin[np.newaxis, :] + kout, axis=1)
if (norm_vec < 0.01 * kin_norm).any(): # check whether any fulfills the criterion
minarg = np.argmin(norm_vec)
res_list.append([theta, phi_test[minarg]])
Ghkl = np.dot(rot_mat,Ghkl)
print('Time was {0:1.2f} s'.format( (time.time()-start_time)))
# On my machine it takes 0.3s and res_list should be
# [[1.0356115365192968, 1.578689775673263]]你知道用一种更快的方法来计算这个问题吗?从概念上讲，要么通过求解完全不同的方程，要么用我现有的方法使它更快？
回答 2
Stack Overflow用户
发布于 2016-09-27 14:23:33
在每次迭代中都会更新Ghkl并在下一个迭代中重新使用，这就存在依赖关系。相应的封闭形式可能很难找到。因此，我将专注于改进最内部循环中的其余代码的性能。
现在，我们正在计算norm_vec，我认为可以用下面列出的两种方法来加速计算。
Appproach #1使用Scipy's cdist -
from scipy.spatial.distance import cdist
norm_vec = cdist(kout,(kin-Ghkl)[None]) # Rest of code stays the sameAppproach #2使用np.einsum -
sums = (Ghkl-kin)+kout
norm_vec_sq = np.einsum('ij,ij->i',sums,sums)
if (norm_vec_sq < (0.01 * kin_norm)**2 ).any():
minarg = np.argmin(norm_vec_sq) # Rest of code stays the same运行时测试-
使用问题中列出的输入，我们得到了以下结果：
In [91]: %timeit np.linalg.norm(Ghkl- kin + kout, axis=1)
10000 loops, best of 3: 31.1 µs per loop
In [92]: sums = (Ghkl-kin)+kout
...: norm_vec_sq = np.einsum('ij,ij->i',sums,sums)
...:
In [93]: %timeit (Ghkl-kin)+kout # Approach2 - step1
100000 loops, best of 3: 7.09 µs per loop
In [94]: %timeit np.einsum('ij,ij->i',sums,sums) # Approach2 - step2
100000 loops, best of 3: 3.82 µs per loop
In [95]: %timeit cdist(kout,(kin-Ghkl)[None]) # Approach1
10000 loops, best of 3: 44.1 µs per loop因此，方法1没有显示这些输入大小有任何改进，但是方法2正在绕过3x加速计算norm_vec。
关于为什么np.einsum 在这里是有益的原因的简短说明:好吧，einsum对于元素乘法和和约化来说是很棒的。我们正在利用这一问题的本质。np.linalg.norm给我们在输入的平方版本上沿着第二轴进行的求和。因此，einsum对应的方法就是将相同的输入提供给它的两个输入的两倍，从而处理平方，然后用它的字符串表示法来表示它的和约化，失去第二个轴。这两种方法都是一次完成的，这可能就是为什么这里速度如此之快的原因。
Stack Overflow用户
发布于 2022-11-03 21:34:28
只要对原子位置进行傅里叶变换，就可以得到任意原子平面的衍射。傅里叶空间分辨率是由你的整体晶格尺寸来定义的。
有关该实现的更多细节可以在这里找到：https://iopscience.iop.org/article/10.1088/1361-651X/ac860c
我将用存储库更新这个答案，一旦我将它公之于众，它就会计算XRD。
https://stackoverflow.com/questions/39726139
复制相似问题
腾讯云开发者
