24小时热门版块排行榜

南方科技大学公共卫生及应急管理学院2025级博士研究生招生报考通知

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

txcokokok

木虫 (小有名气)

应助: 1 (幼儿园)
金币: 26.8
散金: 1578
红花: 6
帖子: 162
在线: 29.1小时
虫号: 2071953
注册: 2012-10-19
性别: GG
专业: 金属材料的力学行为

[求助] lammps 计算粘度 lj单位转化为 metal单位出错已有3人参与

最近做粘度的计算，修改lammps里面提供的例子，Muller-Plane的方法，但是计算的结果对不上，所以一步步排查，最后发现单位转换这里都没搞定。
原例子都是用lj单位写的，我想把它改成metal单位，就是按照manual里面lj单位的说明改，改了的地方都下划线标出了，想请大家帮忙看看是不是哪里错了。确实很繁琐，我把过程都写出来了，还请大家多多帮忙，也好帮助其他遇到这个问题的人。最后附上全文的word版，红色标出比较明显，还有lj单位和metal单位的in文件和log文件。

http://tieba.baidu.com/p/2483520076
按照上面这个帖子的方法，确定参数如下：
sigma：3.405e-10 m
epsilon： 1.67e-21 J
玻尔兹曼常数Kb：1.38065E-23 J/K
摩尔质量：39.948 gram/mol
单个原子质量：6.63352E-26 kg
(用摩尔质量/阿伏伽德罗常数)

具体代码做了如下修改，lj部分为源代码，metal部分为我自己修改后的，源代码见lammps安装包：lammps/Examples/VISCOSITY/in.mp.2d
lj:
variable x equal 20 #x方向长度
variable y equal 20 #y方向长度
metal:
保持不变

lj:
variable rho equal 0.6 #lj单位中的约化密度（相当于metal单位中的晶格常数）
variable  t equal 1.0 #温度
variable rc equal 2.5 #截断半径
metal:
variable rho equal 6.217 #晶格常数
(由rho*=rho•sigma^dim，和rho=N/V；
二维模拟，dim=2，采用sq2结构所含原子数N=2，V=a^2，a为晶格常数；
得0.6=2•(sigma/a) ^2，所以a=6.217 埃)
variable  t equal 120.96 #温度
(由T*=T•Kb/epsilon ，1.0*1.67e-21/(1.38065E-23)=120.96 K)
variable rc equal 8.5125 #截断半径
(由x*=x/sigma，2.5*3.405e-10 =8.5125埃)

lj:
units lj
dimension 2
atom_style atomic
neigh_modify delay 0 every 1
lattice                   sq2 ${rho}
region                simbox block 0 $x 0 $y -0.1 0.1
create_box          1 simbox
create_atoms       1 box
metal:
保持不变

lj:
pair_style    lj/cut ${rc}
pair_coeff    * * 1 1
metal:
pair_style    lj/cut ${rc}
pair_coeff    * * 1.042332128634570e-2  3.405
(势阱常数1.042332128634570e-2 eV，零势能距离 3.405 埃)

lj:
mass          * 1.0
timestep       0.005
velocity       all create $t 97287
metal:
mass          * 39.948
(质量改为相对摩尔质量)
timestep          0.01
（lj单位默认时间步长0.005，由t* = t (epsilon / (m•sigma^2)^1/2，得lj下默认时间为10fs，即0.01ps)
velocity       all create $t 97287

lj:
fix          1 all nve
fix          2 all langevin $t $t 0.1 498094
fix          3 all enforce2d
metal:
fix          1 all nve
fix          2 all langevin $t $t 0.2 498094
(由damp = damping parameter (time units))
fix          3 all enforce2d

lj:
# equilibration run
thermo          1000
run                5000
unfix 2
# turn on Muller-Plathe driving force and equilibrate some more
velocity all scale $t
fix 4 all viscosity 100 x y 20
fix 5 all ave/spatial 20 50 1000 y center 0.05 vx &
      units reduced file profile.mp.2d
# equilibration run
variable dVx equal f_5[11][3]-f_5[1][3]
thermo          1000
thermo_style custom step temp epair etotal press f_4 v_dVx
run 20000
# data gathering run
# reset fix viscosity to zero flux accumulation
unfix 4
fix 4 all viscosity 100 x y 20
#100步交换一次动量，交换x方向的动量分量，另一个分量在y方向，共20层
metal:
保持不变

lj:
variable visc equal -(f_4/(2*(step*0.005-125)*lx+1.0e-10))/(v_dVx/(ly/2))
#http://lammps.sandia.gov/threads/msg43465.html
#http://www.52souji.net/lammps-command-fix-viscosity/
#f_4输出总的动量交换量，除以时间和截面积即为动量通量；
#(step*0.005-125)为时间，0.005为lj单位默认的时间步长，125为到目前为止（20000+5000=25000）步的时间：25000*0.005=125，为防止25000步计算时分母为0,加上1.0e-10；
#lx为截面积（二维，z方向为1）；
#v_dvx为速度分量，除以ly/2为速度梯度
variable visc equal -((f_4/(2*(step*0.005-125)*lx+1.0e-10))/(v_dVx/(ly/2)))*9.07812e-05
metal:
variable visc equal -((f_4/(2*(step*0.01-250)*lx+1.0e-10))/(v_dVx/(ly/2)))*1.66e-5

lj:
fix vave all ave/time 1000 1 1000 v_visc ave running start 26000
#ave running 每1000步输出一次到目前为止的v_visc的累积平均值，作为vave的值,从26000步开始
thermo_style custom step temp f_4 v_dVx v_visc f_vave
run                50000
metal:
保持不变

回复此楼