今天建模时要做个粒子无法闯过的球面，我采用的方面是在球面上放置overlap的粒子，这些粒子和其他粒子的相互作用采用LJ势，这样就可以保证其他粒子进不了球内，呵呵。那么如何产生均匀的球面呢？

我使用的办法是：
1.先计算表面密度是1.0时大概需要的粒子数N；
2.用MATLAB程序把粒子随机分布到球面上，并加上原点处的粒子写如data数据，这样就有N+1个粒子；
3.所有球面粒子和原点处的组成 N个bonds
4.fixed 原点粒子，采用bond_style ：harmonic  bond_coeff    1 5000.0 【数值设置尽量大】 5.0  保证粒子不远离球面
5.粒子间相互作用采用pair_style           soft 1.0291   pair_coeff      * * 0.0 10.0，利用 NVE 系综run, run的过程中并不断调整截断距离使pair相互作用处于0和非0边缘。这样所有粒子间距离就可以通过对势来限制为相同。

问题的关键点是不断调整截断距离，先小后慢慢增大，不断速度置零，找到温度和对势的边缘（0），这样就可以使粒子间距离几乎相等，从而实现分布均匀。

MATLAB代码：
% 320 atoms of spherical face
clc;clear;
%%%%%%%%%%%%%%%%%%%%%
%  parameters
R = 5;
Nf =320;
%%%%%%%%%%%%%%%%%%%%%
sita = pi*rand(1,Nf);
fai = 2*pi*rand(1,Nf);
x = R * (sin(sita).*cos(fai));
y = R * (sin(sita).*sin(fai));
z = R * cos(sita);
XYZ =[ x' y' z'];
plot3(x,y,z,'r*')

部分LAMMPS代码：我run这个代码N次，# 处是第一次run使用的

bond_style           harmonic
special_bonds        0.0 1.0 1.0
read_restart         face1.equ

bond_coeff           1 5000.0  5.0
group                org  id 321
group                face id <= 320
pair_style           soft 1.0291
pair_coeff           * * 0.0 10.0
thermo               500
velocity             face  set 0 0 0 units box
#velocity             face  create 2.0 4928459  units box
fix                  1 face nve
#fix                  2 face langevin 2.0 0.25 10.0 9887415
fix                  3 org  nve/noforce
run                  10000
#unfix                2
velocity             face  set 0 0 0 units box
pair_coeff           * * 10.0 10.0
run                  10000
velocity             face  set 0 0 0 units box
dump                 1  all custom 10  atom.dat  tag  mol  type  q  x  y  z
run                  20


[ 本帖最后由 redream 于 2008-7-4 14:19 编辑 ]

真的不错啊

恩 很有启发

嗯，楼主是利用lammps作的，如果用XMD好像就简单的了！XMD中有一个constrain的命令。好像是：CONSTRAIN CAVITY SPHERE spring xcenter ycenter zcenter radius
我也是才接触MD，楼主见笑了！

   
很有意思

有意思，不过有点不明白，为什么刚开始你不让粒子在球面上比较均匀分布，而让它随机在球面上分布？

[ 本帖最后由 eyra 于 2008-7-9 10:48 编辑 ]

问的好，你可以搞个均匀分布的程序，看看是不是相当的麻烦，如何把sita,fai角合适的分配才能达到均匀分布，这个问题值得想想

平均分布我觉得可以采用图形变换来做，
不必一定采用极坐标系。
因为球面是一个约束条件，
而球面的引入导致了球面上的点虽然采用三维表示，
但是实际上是二维的点，
所以可以考虑降维。
因此要均匀分布，做两个变换就可以，
首先是降维的变换，
即第一个变换A把球面变换到平面曲边正方形，
然后是坐标变换，
即第二个变换B把曲边正方形变换为标准的单位正方形，
然后是原子坐标分布，
在这个单位正方形里平均分布原子坐标C，
最后是逆变换还原三维坐标，
用A^{-1}B^{-1}C就可以求出球面上平均分布的原子坐标。
基本上都是计算机图形学的内容，
我想可能用opengl或者d3d的贴图功能也能直接实现。

[ 本帖最后由 realasking 于 2008-7-12 12:48 编辑 ]

真的不错啊

  有创意

一个小建议：　用四元数来处理球面上的点是非常方便的

用四元数来处理球面上的点是非常方便的？？？

说说四元函数怎么处理啊

可以参考书目The Art of Molecular Dynamics Simulation@2004Cambridge。

它可以做到在球面上均匀随机的分布点，且更方便的对这个点进行任意你想要的旋转位移。只是用它来做你提出事情，有点大材小用的意思。

很有意思啊