相对角距离方法的Matlab实现

之前过冷水在推文中三维空间分布函数绘制实例中和大家分享了对分布函数g(r)的程序实现方法。只要你认真学习专研总有新的发现，这不过冷水就接触到了一种叫做相对角距离的方法，应用该方法可以得到一个完整的峰值函数，了解液态结构的应该知道称之为第一配位球层对分布函数。图像如下：

相对角较算法：该方法给出了判断以i原子为中心,j原子在其配位层内的条件。对于原子i，如果对所有原子k满足关系式，则认为原子j是在i的配位层内。

 

将上述的相对角公式带入到我们之前定义的对分布函数公式

我们就可以得到复合相对角方法对分布函数

至此相对角方法介绍完毕，公式就是这么简洁，有问题的是需要如何编程实现？

在这里我们再重新讲一下对分布函数g(r)的编程思路

（1）采用循环的方法统计所有原子i和原子j的距离，将所有距离划入到不同的具体梯度内，统计在对应梯度的个数，统计不同梯度的距离所占的百分比    

、、

（2）因为要求平均所以需要对百分比除以一个平均密度，进行归一化。

（3）根据不同距离梯度生成合适的梯度距离数据，对距离数据和概率密度数据进行插值拟合，就得到的对分布函数图像。

因为提出了相对角距离判断公式所以需要更改第一步的统计对应梯度的个数

比如说说以前ij距离为7梯度间距是为0.1那么就在第70个梯度区间，梯度区间计数+1。

更改后的为ij距离为7，ijk之间的关系不满足相对距离判断条件。所以所以梯度计数保持不变。剩余的处理流程和之前一致，所以代码为：

clear
tic
height=11.1087999344000004;width=11.1087999344000004;long=11.1087999344000004;
load('data.mat')
N=52;
n=60;  %划分区间个数
rc = sqrt(width^2+height^2+long^2); %搜索圆的最大半径
dr = rc/n;%确定半径区间
num = round(rc/dr);
[m,n]=size(a);
%读数据晶型区间划分，每32个坐标数据为一个结构，储存在元胞中
for i=1:length(a)/N
   b= a(1+(N*(i-1)):N*i,:);
   Centers{i,1}=b;
end
%划分距离区间
gr=zeros(num,1);
for n=1:length(Centers)
    centers=Centers{n,1};
    particle_num=length(centers);
    [row,col] = size(centers);
    for i=1:(row-1)
        for j = i+1:row
            %这两个语句的目的是剔除掉原子i,j的坐标
            centersij=centers;
            centersij([i,j],:)=[];
            %增加的判据语句RAD
            a =sqrt(sum((centers(i,:)-centers(j,:)).^2));
            b=sqrt((centers(j,1)-centersij(:,1)).^2+(centers(j,2)-centersij(:,2)).^2+(centers(j,3)-centersij(:,3)).^2);
            c=sqrt((centers(i,1)-centersij(:,1)).^2+(centers(i,2)-centersij(:,2)).^2+(centers(i,3)-centersij(:,3)).^2);
            %b=sqrt(arrayfun(@(a)((centers(j,1)-a).^2),centersij(:,1))+arrayfun(@(a)((centers(j,2)-a).^2),centersij(:,2))+arrayfun(@(a)((centers(j,3)-a).^2),centersij(:,3)));
            %c=sqrt(arrayfun(@(a)((centers(i,1)-a).^2),centersij(:,1))+arrayfun(@(a)((centers(i,2)-a).^2),centersij(:,2))+arrayfun(@(a)((centers(i,3)-a).^2),centersij(:,3)));
            cosa=(b.^2+c.^2-a.^2)./(2*b.*c);
            RAD=1/a.^2-cosa./b.^2;
            distance =sqrt(sum((centers(i,:)-centers(j,:)).^2));
            if distance <= rc%计算
                lane = round(distance/dr);%将粒子划分不同梯度内
                %进行计数前加一个判据条件
                if  RAD>=0 
                    gr(lane) = gr(lane)+1;%做个数累计
                else
                    continue;
                end
            end
        end
    end
end
%绘制径向分布函数
[row,col] = size(gr);
percent = zeros(row,1);
gobal_rho = particle_num / (3/4*rc^3);%全局数密度
for i=1:row
    temp = gr(i);%不同半径下的单位原内的原子个数
    temp = temp / (particle_num*length(Centers));
    temp = temp / (4/3*pi*((i*dr)^3-((i-1)*dr)^3));%某一个半径梯度下的局部密度
    percent(i) = temp / gobal_rho;
end
 
index=nonzeros(linspace(0,rc,num+1))';
percent = reshape(percent,1,row);
xlim([0, 20]);
ylim([0, 20]);
values = spcrv([[index(1) index index(end)];[percent(1) percent percent(end)]],3)';
r=values(:,1);g=values(:,2);
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(index,percent,'LineWidth',3,'Color',[1 0 0]);
plot(r,g,'LineWidth',3,'Color',[0 0 1]);
plot(index,percent,'MarkerFaceColor',[0 1 0],'MarkerSize',8,'Marker','o','LineStyle','none','Color',[0 1 0]);
box(axes1,'on');
hold(axes1,'off');
set(axes1,'FontName','Times New Roman','FontSize',15,'LineWidth',3);

运行代码就可以得到对分布函数图像。我们求第一配位层的对分布函数图像主要就是为了求配位数，过冷水一并给出配位数的计算公式。所谓的配位数就是4πr²ρ₀g(r)与Y轴包围的面积。

在被积函数4πr²ρ₀g(r)相当复杂时，就只能采取数值积分的求法。选择采用大数定理 ：在（a,b）区域内均匀随机抽样得到N个点x₁、x₂、x₃、...、x_N;求这些点上被积函数的值f(x₁)、f(x₂)、f(x₃)......f(x_N)、f(x₁)、于是f(x)在[a,b]区域的平均值：

%% 配位数计算
load('data.mat');
fRAD= fit(r,g, 'smoothingspline');
rgRAD=linspace(r(1),r(end),1000)';
intgtotal_cal=rho*4*pi*rgRAD.^2.*fRAD(rgRAD);
z=(max(rgRAD)-min(rgRAD)).*(sum(intgtotal_cal))./1000;

过冷水作为动力学的初学者也是在不断摸索中。最开始第一次自己编出对分布函数图像，虽然现在看来好简单，当时就对我可难了，能做出来特别有成就感，学习能够让我感到快乐，希望大家阅读过冷水的推文也能够感觉到快乐。本期内容就这么多，欢迎大家留言讨论。

往期回顾>>>>>>

大数定理的小陷阱

sym、sum妙用案例

 一文带你轻松搞定matlab各种数据读写

过冷水带您走进matlab数据转换新世界

十万个matlab编程问题征集，欢迎来问

MATLAB编程爱好者Q群正式升级为2000人群，欢迎加入！