分子动力学的径向分布函数绘制

3年前浏览4804

这是过冷水之前在工作中的一张图片。这张图片的实际含义是以坐标原点为参考点，距离原点距离x处的小球个数f(x)分布图。这个问题理解起来好容易，统计距离和小球个数的关系，but实际问题是过冷水要处理的是这样的问题

请问请告诉我随意找一个小球为原点，距离其x处的小球个数怎么求？随意找一个小球为原点。距离x处的小球平均个数怎么求？这样的问题在物理中会经常出现，过冷水就带领大家讲复杂问题简单化，1：首先要读取图片，将图片上的圆数据化。2：以半径r进行搜索圆的个数，将整个图片上的原按照r距离划分到不同距离区间中(dr=r/N);3:小球之间统计距离，划分到对应区间中。

代码如下

clc;
clear
image = imread('D:\桌面\小球分布.png');
[height,width,channel] = size(image);%确定输入突变的长宽
rc = sqrt(width^2 height^2)/2; %搜索圆的最大半径
n=60;  %划分区间个数
dr = rc/n;%确定半径区间
%检测输入图片上的圆
gray = rgb2gray(image);  %灰度化
bw = imbinarize(gray);   %二值化
bw_inverse = imcomplement(bw); %二值图像反转，即0变1,1变0
[centers,radius] = imfindcircles(bw_inverse,[5,20]); %检测图像上半径在5~20之间的圆，也可以是其它，看实际情况需要
particle_num = length(radius(:)); %统计圆的个数
%计算每个单位圆内的原子个数
[row,col] = size(centers);
num = round(rc/dr);
gr=zeros(num,1);%任意一个粒子为中心，不同半径梯度下的原子个数
for i=1:(row-1)
    for j = i 1:row
        distance = sqrt((centers(i,1)-centers(j,1))^2   (centers(i,2)-centers(j,2))^2);%计算每两个原子间的距离
        if distance <= rc%计算
            lane = round(distance/dr);%将粒子划分不同梯度内
            gr(lane) = gr(lane) 1;%做个数累计
        end
    end
End
r=nonzeros(linspace(0,sqrt(width^2 height^2),row 1))'

计算得到的r~g(r)就是我们想要的距离和小球个数的统计结果画图可得：

这幅图看上去比较怪不太容易看出来规律，这只是第一问的结果，小球个数也是以每个小球为中心累计统计叠加出来的总的总的结果，那么请问以一个小球问中心其余小球分布到底应该怎么算？这就需要求平均了，怎么说呢？过冷在思考这个问题的时候为了形象理解抽象出了这么一个问题：

2.2个小孩看起来像个笑话。但这就是数学，不知读者能否理解这个2.2，如果这个回答不对，那么准确的回答应该怎么说？在该问题的基础上，过冷水又构思了这么一个问题：

读者只有理解了上述问题才能理解过冷水接下来的操作，在求得的总f(r)基础上除以n(小球个数)

这样我们就得到任意一个小球为圆心距离r处圆内的小球的个数了。在统计中我们更关心的是概率，在r处的找到小球的概率，而不是r处圆内的小球个数。（半径梯度对应的是圆环面积）

程序处理为：

[row,col] = size(gr);
percent = zeros(row,1);
for i=1:col
    temp = gr(1,i);%不同半径下的单位原内的原子个数
    temp = temp / particle_num;
    percent(1,i)= temp / (pi*((i*dr)^2-((i-1)*dr)^2));%某一个半径梯度下的局部密度
end
figure1 = figure;
axes1 = axes('Parent',figure1,'Position',[0.131090512540894 0.11 0.775 0.815]);
hold(axes1,'on');
plot(r,percent,'MarkerFaceColor',[0 0.447058826684952 0.74117648601532],'MarkerSize',10,'Marker','o','LineWidth',1,'LineStyle','--');
ylabel('\rho');
xlabel('r');
box(axes1,'on');
hold(axes1,'off');
set(axes1,'FontSize',14,'LineWidth',3);[object Object]