在计算材料学中极少的情况下,我们可能会需要随机产生一个模拟盒子内的原子结构,一般用于测试。
这样说,其实只要rand一个数组就可以了,但是我们又希望这个结构又能够大致满足一些物理上的限制,而不是纯粹数学上的随机,比如原子间距要比较合适。
这里提供的程序就是考虑了原子间距,而产生随机结构的。
function randStructure(numAtom,boxlength,dmin,dmax)
% Geanerate a structure randomly with the distance in specified range
%
% randStructure(numAtom,boxlength,dmin,dmax)
% numAtom: number of atoms in the structure.
% boxlength: length of the simulation box
% dmin: minimum of the distance between two atoms
% dmax: [optional] maximum of the distance between two atoms
%
% Recommendation:
%It's better to set dmax equaling to boxlength.
%
% Example:
% randStructure(10,4,2,5)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% log:
% 2012-12-10: Complete
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin == 3
dmax = boxlength;
end
%function randStructure()
%numAtom = input('Number of atoms:');
%boxlength = input('Scale factor:');
%dmin = input('minimum of the distance between two atoms:');
%dmax = input('maximum of the distance between two atoms:');
coord = rand(1,3)*boxlength; % the 1st atoms
n = 1; % number of atoms
while n < numAtom
tcoord = rand(1,3)*boxlength;
nlen = size(coord,1);
counter = 0;
for i=1:nlen
tdist = norm(coord(i,:)-tcoord);
if tdist >= dmin && tdist <= dmax
counter = counter + 1;
end
end
if counter == nlen
coord(end+1,:) = tcoord;
n = n + 1;
end
end
coord/boxlength
关于以上程序,有几点是需要特别说明的:
这里模拟的盒子为正交,并且各轴长度相等的。
如果第四个参数,即dmax不输入,默认赋boxlength。
最后一行将坐标转换成了分数坐标。如果需要绝对坐标,不除以boxlength就可以了。
在有些情况下,程序似乎不能结束。这是因为第一个原子的位置产生的不太好,使得在一个有限的盒子内,无法产生距离满足要求的其他原子。这时可以把程序终止掉,重新运行。如果这样做几次仍然不行,就要考虑是不是参数设置有问题,比如盒子太小,或者原子间距离太大。