单晶衍射仪系列3-衍射过程仿真

casjxm

单晶衍射仪主要部件包括X射线源，测角仪，探测器。X射线源一般是固定的，有时需要根据实际需求选择测角仪与探测器的合理配置，这时就需要对衍射过程进行仿真，输入测角仪，探测器以及样品的参数进行模拟测试，再看看测试数据是否满足需求。前面已经在相关帖子中说明了衍射原理以及晶体学的基本计算公式。


直接上代码（matlab）:
1. 计算面间距
参考：https://www.matstr.com/forum.php?mod=viewthread&tid=564&extra=
https://www.matstr.com/forum.php?mod=viewthread&tid=383&extra=

function result=hkldist(hkl, cellpara)
%hkl: row vector
        h=hkl(1);
        k=hkl(2);
        l=hkl(3);
        a=cellpara(1);
        b=cellpara(2);
        c=cellpara(3);
        d=cellpara(4)/180*pi;
        e=cellpara(5)/180*pi;
        f=cellpara(6)/180*pi;
                v=a*b*c*sqrt(1-cos(d)*cos(d)-cos(e)*cos(e)-cos(f)*cos(f)+2*cos(d)*cos(e)*cos(f));
                s11=b*b*c*c*sin(d)*sin(d);
                s22=a*a*c*c*sin(e)*sin(e);
                s33=a*a*b*b*sin(f)*sin(f);
                s12=a*b*c*c*(cos(d)*cos(e)-cos(f));
                s23=a*a*b*c*(cos(e)*cos(f)-cos(d));
                s13=a*b*b*c*(cos(f)*cos(d)-cos(e));
                dist=v/sqrt(s11*h*h+s22*k*k+s33*l*l+2*s12*h*k+2*s23*k*l+2*s13*l*h);
result=dist;

2. 测角仪（晶体）旋转
参考：https://www.matstr.com/forum.php?mod=viewthread&tid=491&extra=

function result=rot_all_mat(w,k,p,ok)
res=rot_phi_mat(w,k,p,ok)*rot_kap_mat(w,k,ok)*rot_omega_mat(w);
result=res;


3. 直空间与倒空间单胞参数的转换
参考：https://www.matstr.com/forum.php?mod=viewthread&tid=382&extra=

function result=realcellpara2recicellpara(realcellpara)
%realcellpara[6]
a=realcellpara(1);
b=realcellpara(2);
c=realcellpara(3);
alpha=realcellpara(4)/180*pi;
beta=realcellpara(5)/180*pi;
gama=realcellpara(6)/180*pi;
v=a*b*c*sqrt(1-cos(alpha)^2-cos(beta)^2-cos(gama)^2+2*cos(alpha)*cos(beta)*cos(gama));
a2=b*c*sin(alpha)/v;
b2=a*c*sin(beta)/v;
c2=a*b*sin(gama)/v;
alpha2=asin(v/(a*b*c*sin(beta)*sin(gama)))/pi*180;
beta2=asin(v/(a*b*c*sin(alpha)*sin(gama)))/pi*180;
gama2=asin(v/(a*b*c*sin(alpha)*sin(beta)))/pi*180;
v2=1/v;
res=[a2 b2 c2 alpha2 beta2 gama2 v2];

result=res;


4. 晶体学坐标系到直角坐标系转换
参考： https://www.matstr.com/forum.php?mod=viewthread&tid=381&extra=
function result=abc2xyz_mat(cellpara)

%x: row [x1,x2,x3]
%cellpara：row
%res: row [res1,res2,res3]
a=cellpara(1);
b=cellpara(2);
c=cellpara(3);
alpha=cellpara(4)/180*pi;
beta=cellpara(5)/180*pi;
gama=cellpara(6)/180*pi;

m=[a 0 0;
b*cos(gama) b*sin(gama) 0;
c*cos(beta) c*(cos(alpha)-cos(gama)*cos(beta))/sin(gama) c*sqrt(sin(gama)^2+2*cos(alpha)*cos(beta)*cos(gama)-cos(beta)^2-cos(alpha)^2)/sin(gama)];
% res=x*m
result=m;

5. 衍射点在直角坐标系中的位置

function result=hkl2xyz(hkl,omega,kappa,phi,phiok,reci_orien_mat)
%hkl: row vector
%res: row vector
%reci_orien_mat=[da*/dx db*/dx dc*/dx; da*/dy db*/dy dc*/dy;da*/dz db*/dz dc*/dz;]
rotmat=rot_all_mat(omega,kappa,phi,phiok);
res=hkl*reci_orien_mat'*rotmat';
result=res;

6. 衍射点直角坐标对应的球坐标角度
参考：https://blog.csdn.net/qq_35097289/article/details/94003524

function result=xyz2phisita(xyz)
r=norm(xyz);
phi=acos(xyz(3)/r);
sita=acos(xyz(1)/r/sin(phi));
res=[phi/pi*180,sita/pi*180];
result=res;

7. 生成单张衍射照片中的hkl

function result=creatframehkl(reci_orien_mat,peakposierr,lamda,resol,cellpara,omega,kappa,phi,phiok,phiSitaRange)
%reci_orien_mat=[da*/dx db*/dx dc*/dx; da*/dy db*/dy dc*/dy;da*/dz db*/dz dc*/dz;]
hkl=creatspherehkl(resol,cellpara);
rotmat=rot_all_mat(omega/180*pi,kappa/180*pi,phi/180*pi,phiok/180*pi);
n=size(hkl);
xyz=hkl*reci_orien_mat'*rotmat';

r2=(1/lamda)^2;
ewalda=1/lamda;
xyz2=xyz(:,1)-ewalda;
xyz3=[xyz2,xyz(:,2:3)];
xyz4=xyz3.^2;
xyz5=abs(sum(xyz4')'-r2);

xyz6=[];
hkl2=[];
for i=1:n(1)
        if(xyz5(i)<peakposierr)
                xyz6=[xyz6;xyz(i,:)];
                hkl2=[hkl2;hkl(i,:)];
        end
end

m=size(hkl2);
hkl3=[];
for i=1:m(1)
        phiSita=xyz2phisita(xyz6(i,:));
        if((phiSita(1)>phiSitaRange(1))&&(phiSita(1)<phiSitaRange(2))&&(phiSita(2)>phiSitaRange(3))&&(phiSita(2)<phiSitaRange(4)))
        hkl3=[hkl3;hkl2(i,:)];
        end
end        

result=hkl3;


8. 生成全ewald球衍射点

function result=creatspherehkl(resol,cellpara)
        res=[];
        for i=1:100
                for j=1:100
                        for k=1:100
                        hkl=[i,j,k];
                        if(hkldist(hkl,cellpara)>resol)
                                res=[res;hkl];
                        end
                        end
                end                        
        end

result=res;


9. 倒空间与直空间取向矩阵的变换
参考：https://www.matstr.com/forum.php?mod=viewthread&tid=563&extra=

function result=reciorienmat2realorienmat(reci_orien_mat)
%reci_orien_mat=[da*/dx db*/dx dc*/dx; da*/dy db*/dy dc*/dy;da*/dz db*/dz dc*/dz;]
recicellpara=orienmat2cellpara(reci_orien_mat);
va2=reci_orien_mat(:,1)';
vb2=reci_orien_mat(:,2)';
vc2=reci_orien_mat(:,3)';
va=cross(vb2,vc2)/recicellpara(7);
vb=cross(vc2,va2)/recicellpara(7);
vc=cross(va2,vb2)/recicellpara(7);
real_orien_mat=[va' vb' vc'];

result=real_orien_mat;


10. 取向矩阵与单胞参数的转换

function result=orienmat2cellpara(orien_mat)
%reci_orien_mat=[da*/dx db*/dx dc*/dx; da*/dy db*/dy dc*/dy;da*/dz db*/dz dc*/dz;]
%real_orien_mat=[da/dx db/dx dc/dx; da/dy db/dy dc/dy;da/dz db/dz dc/dz;]

va=orien_mat(:,1)';
vb=orien_mat(:,2)';
vc=orien_mat(:,3)';
a=norm(va);
b=norm(vb);
c=norm(vc);
alpha=acos(dot(vb,vc)/b/c)/pi*180;
beta=acos(dot(vc,va)/c/a)/pi*180;
gama=acos(dot(va,vb)/a/b)/pi*180;
v=a*b*c*sqrt(1-cos(alpha/180*pi)^2-cos(beta/180*pi)^2-cos(gama/180*pi)^2+2*cos(alpha/180*pi)*cos(beta/180*pi)*cos(gama/180*pi));

res=[a b c alpha beta gama v];
result=res;


11. 探测器测试范围

function result=detectorrange(width,height,sita2,distance,centery,centerz)
s_phi1=90-atan((height/2+centerz)/distance)/pi*180;
s_phi2=90+atan((height/2-centerz)/distance)/pi*180;
s_sita1=180-atan((width/2+centery)/distance)/pi*180-sita2;
s_sita2=180+atan((width/2-centery)/distance)/pi*180-sita2;
res=[s_phi1,s_phi2,s_sita1,s_sita2];
result=res;

12. omega扫描测试
function result=scan_omega(reci_orien_mat,peakposierr,lamda,resol,cellpara,stepangle,anglerange,phiok)
%anglerange=omegamin,omegamax,kappa,phi

%Eiger4M: 155mmX163mm
%pilatus2M: 254mmX289mm（HxV）
width=254;
height=289;
sita2=0;
distance=70;
centery=0;
centerz=0;

phiSitaRange=detectorrange(width,height,sita2,distance,centery,centerz);
m=size(anglerange);
nrun=m(1);
m=(anglerange(:,2)-anglerange(:,1))/stepangle;
nframe=floor(m)+1;
hkl=[];
multiple=[];
for i=1:nrun
        runhkl=[];
        for j=1:nframe(i)
                omega=anglerange(i,1)+(j-1)*stepangle;
                kappa=anglerange(i,3);
                phi=anglerange(i,4);
                framehkl=creatframehkl(reci_orien_mat,peakposierr,lamda,resol,cellpara,omega,kappa,phi,phiok,phiSitaRange);
                n=size(framehkl);
                if(n(1)>=1)
                        plot3(framehkl(:,1),framehkl(:,2),framehkl(:,3),'b*');
                end
                pause(2);
                runhkl=[runhkl;framehkl];
        end
        [uniqueRows, ~, rowIndices] = unique(runhkl, 'rows');
        rowCounts = accumarray(rowIndices, 1);
        hkl=[hkl;uniqueRows];
        multiple=[multiple;rowCounts];
                        n=size(hkl);
                if(n(1)>=1)
                         plot3(hkl(:,1),hkl(:,2),hkl(:,3),'r*');
                 end
                 pause(2);
        
end
        
        [uniqueRows, ~, rowIndices] = unique(hkl, 'rows');
        rowCounts = accumarray(rowIndices, 1);
        m=size(hkl);
        n=size(uniqueRows);
        for i=1:n(1)
                for j=1:m(1)
                        if(rowIndices(j)==i)
                                rowCounts(i)=rowCounts(i)+multiple(j)-1;
                        end
                end
        end

        hklm3=[];
        for i=1:n(1)
                if(rowCounts(i)>=3)
                        hklm3=[hklm3;uniqueRows(i)];
                end
        end
        n3=size(hklm3);
        
        spherehkl=creatspherehkl(resol,cellpara);
        p=size(spherehkl);
        complete=n(1)/p(1);
        complete3=n3(1)/p(1);
res=[complete,complete3];

result=res;


13. 初始化并进行测试
function result=initiallzation()
realcellpara=[6 7 8 90 106.5 90 1000];
stepangle=0.5;
%anglerange=omegamin,omegamax,kappa,phi
anglerange=[20 50 45 60;30 60 45 60;70 90 45 80];
phiok=45;
resol=0.83;
lamda=0.71073;
peakposierr=3;
recicellpara=realcellpara2recicellpara(realcellpara);
%reci_orien_mat=[da*/dx db*/dx dc*/dx; da*/dy db*/dy dc*/dy;da*/dz db*/dz dc*/dz;]
%real_orien_mat=[da/dx db/dx dc/dx; da/dy db/dy dc/dy;da/dz db/dz dc/dz;]
abctoxyzmat=abc2xyz_mat(realcellpara);
real_orien_mat=abctoxyzmat';
reci_orien_mat=reciorienmat2realorienmat(real_orien_mat);

complete=scan_omega(reci_orien_mat,peakposierr,lamda,resol,realcellpara,stepangle,anglerange,phiok)
result=complete;

%%%%%%%%%%%%%%%%%%%%%%%%
测试过程，在matlab中运行initiallzation：
complete =
    0.8655    0.8480

搁浅

有意思