我的目的就是要模拟两个原子的拉伸状态，其中一个原子固定，另一个原子被施加位移边解条件，就是强迫他走一定的距离。原子势用的是lenard-jones potential. (指数是12 和 6的那个）。但是我用附件的程序运行的结果是无穷大。为什么？

程序中，yupxilong是1， sigma是1，截断距离是2.5，即是说原子间距离超过2.5就没有力了。两个原子初始状态距离是1，末了状态距离是2。




因为这个问题很简单，我希望大家把它给讨论明白了，好谢谢。

[ 本帖最后由 redream 于 2008-7-23 23:40 编辑 ]

yupxilong = epsion

Administrator@TIANWD ~
$ lmps<in.try
LAMMPS (21 May 2008)
dimension       2
units             real
atom_style          atomic
boundary        p p p
neighbor              0.5  bin
neigh_modify      every 1

#read in coordinate information
read_data       data.try
Reading data file ...
  orthogonal box = (0 0 -0.5) to (1 1 0.5)
  1 by 1 by 1 processor grid
  2 atoms

pair_style          lj/cut 2.5
pair_coeff      *  *  1   2.5   #这个相互作用是排斥的啊，2.5*powe(2,1/6)= 2.75>2.5(r_cut)，你再想想初态 应该给什么样的

group           upper id  2
1 atoms in group upper
group           lower id  1
1 atoms in group lower
velocity        all    create 0.01 887723

displace_atoms  upper move 0 0.02 0 units box # 即使移动这么个距离，排斥还是很大，你原子都不知道跑哪了
System init for displace_atoms ...
Displacing atoms ...
fix             hold lower planeforce 0 1 0
fix             3  all nve

# run
timestep              0.005
thermo                10

dump                    1 all xyz 10 dump.try
run                       5000
Setting up run ...
Memory usage per processor = 1.33614 Mbytes
Step Temp E_pair E_mol TotEng Press
       0         0.01 5.8207661e+25            0 5.8207661e+25 2.3947242e+31【看看你的作用势，你怎么不看之前的建议改正呢】
      10          nan    482206.16            0          nan          nan
      20          nan    482206.16            0          nan          nan
      30          nan    482206.16            0          nan          nan
      40          nan    482206.16            0          nan          nan
      50          nan    482206.16            0          nan          nan  
Loop time of 0.512 on 1 procs for 5000 steps with 2 atoms

Pair  time (%) = 0.0130051 (2.54006)
Neigh time (%) = 0 (0)
Comm  time (%) = 0.026009 (5.07989)
Outpt time (%) = 0.461 (90.0392)
Other time (%) = 0.0119854 (2.3409)

Nlocal:    2 ave 2 max 2 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost:    118 ave 118 max 118 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs:    53 ave 53 max 53 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 53
Ave neighs/atom = 26.5
Neighbor list builds = 0
Dangerous builds = 0

这种错误是应该的，你的排斥相互作用很大啊

[ 本帖最后由 redream 于 2008-7-23 09:56 编辑 ]

还有
velocity        all    create 0.01 887723  不对，只能给第二个原子速度，你这个命令都给了
用velocity  upper set 命令

fix             3  all nve不对，第一个不动,用fix  1 lower nve/noforce or fix setforce 把力值为零
第二个fix nve


一句话，你不是很理解LAMMPS的命令，想当然的用一些你认为能实现你目的的命令实现你的目的，其实LAMMPS执行的结果根本不是你想像的

建议：静心好好看看常用的LAMMPS命令；看看别人是如何选择势函数和相应的参数的

[ 本帖最后由 redream 于 2008-7-23 09:54 编辑 ]

我有一些不同的理解：
我用的是
Pair Coeffs
1  1  1   2.5

而不是pair_coeff，所以1（# of atom types)   1(epslon)   1(sigma)     2.5(cut distance)
请参考lammps精简本手册第9页。
所以，在力为零的时候r=sigma*powe(2,1/6)=powe(2,1/6)=1.122462, 当r<1.122462是斥力，当
r>1.122462是吸引力。对不对？

按你的势设置，你在code的第一行加上echo screen 输出到屏幕，把运行的结果让我看看。

我还是建议你仔细琢磨一下LAMMPS的命令。你把atom的dump文件用vmd看看第一个原子是不是运动的

[ 本帖最后由 redream 于 2008-7-23 15:06 编辑 ]

引用:

原帖由 sunysb 于 2008-7-23 06:50 发表
我的目的就是要模拟两个原子的拉伸状态，其中一个原子固定，另一个原子被施加位移边解条件，就是强迫他走一定的距离。原子势用的是lenard-jones potential. (指数是12 和 6的那个）。但是我用附件的程序运行的结果是 ...

建议你做拉伸之前将系统平衡一下，检查一下系统平衡后的状态（粒子速度，能量，原子间距，两原子间相互作用，等等）。
另外，你的边界条件是p p p ，这表示你的体系里有无数个这样的镜像。

[ 本帖最后由 lwan 于 2008-7-23 19:34 编辑 ]

请问Iwan，什么命令可以将系统平衡呢？

不是什么命令使系统平衡，是 先不加外力或其他让原子拉伸，运行一段时间，所谓的驰豫时间，后达到平衡态，然后加力或其他做非平衡拉伸模拟

也就是说，你初态给的是随机的，这个并不是系统的平衡状态。

模拟就是像做试验一样，实际试验的初始条件是什么样的，这就是初态。我们模拟时给的初态不是试验上的初态，要先让给的初态达到真正的试验初态才可以

说的不好，希望你能理解，试验只是打比方

你可以发帖求一个拉伸方面的简单例子，先学学人家怎么做的。
论坛QQ群不知道你加入么 18632934

redream您说的很好，但是我有一点不明白，比如说，我给了一个原子在某一位置，而让另一原子固定，这时两个原子靠lennard-jones势互相吸引或排斥，所以当我松开手之后他们两个就像在作简谐振动一样，振动个不停不是吗？他们怎么
会到平衡位置后停下来呢？





这次程序我按照您说的改了一下，可以运行了，就是一直在作振动，可以看到有原子在动。

# 3D CNT simulation of uniaxial tension

echo                   screen
dimension           2
units                    real
atom_style          atomic
boundary            p p p
timestep              0.005

pair_style           lj/cut   2.5
neighbor            0.3     bin
neigh_modify     every   1

#read in coordinate information
read_data          data.try
group upper       id    2
group lower        id    1

# 预平衡
fix                       1   lower nve/noforce
fix                       2   upper  nve
fix                       3   upper  temp/rescale 10  0.1  0.0  0.02  0.5
run                     5000
fix                       3   upper  temp/rescale 10  0.0  0.0  0.02  0.5
run                     5000
unfix                   1
unfix                   2
unfix                   3

# 如果不放心，把二的速度置零，因为1不动，不用管
velocity              upper  set  0  0  0  units box

#  拉伸模拟，lower 不动，upper移动一定距离，驰豫
displace_atoms upper move 0   1.0   0   units box
fix                       1   lower nve/noforce
fix                       2   upper  nve

timestep             0.0005
thermo               10
dump                 1  all   xyz   10   dump.try
run                     50000

=============================================
下面的代码你想做什么，没看懂
displace_atoms upper move 0 1.0 0 units box 移动上面的原子
fix             hold lower planeforce 0 1 0             且不停的拉下面的原子
#fix              hold lower nve/noforce                 你不是保持下面的不动么
fix             3  upper nve

fix             hold lower planeforce 0 1 0  

    * ID, group-ID are documented in fix command
    * lineforce = style name of this fix command
    * x y z = 3-vector that is normal to the plane

Adjust the forces on each atom in the group so that it's motion will be in the plane specified by the normal vector (x,y,z). This is done by subtracting out components of force perpendicular to the plane.

ID的名字是hold
Group-ID的名字是lower
planeforce是fix命令的style
这个命令是要把lower这组原子的运动限制在（0，1，0）这个平面上，无论发生什么
这个原子总在XZ平面内运动，然而upper这个原子是沿着y轴向上运动，这样一来就是一端固定，另一端作往复运动了。

我运行了redream的程序，可是结果是
ERROR: Computed temperature for fix temp/berendsen cannot be 0.0
可是我们用的是fix temp/rescale 并没有用fix temp/berendsen  。
这是怎么回事？

void FixTempRescale::end_of_step()
{
  double t_current = temperature->compute_scalar();
  if (t_current == 0.0)
    error->all("Computed temperature for fix temp/berendsen cannot be 0.0");
  double delta = update->ntimestep - update->beginstep;
  delta /= update->endstep - update->beginstep;
  double t_target = t_start + delta * (t_stop-t_start);
  // rescale velocity of appropriate atoms if outside window

代码中提示错误，估计steve还没有注意到，应该改为temp/rescale的