详细内容

LAMMPS讲解57-LAMMPS的data文件导出与in文件编写

MS中导出LAMMPSdata文件

在导出模型之前先做两件事,第一是把模型boxZ方向加大一点,可以避免Z方向上由于周期性而导致非物理的作用和现象发生,也就是说避免基底在液滴上方对水分子有作用,第二是给模型赋力场。

原来的模型Z方向box53.72Å,现在扩大为200 Å,如下图所示,注意中间的图中不要勾选。

image.png 

接下来是给模型赋力场,这是后面使用lammps自带的程序将模型转为data文件的关键一步,否则程序无法识别原子。注意按照下面图片中的步骤来,改勾选的地方勾选,不该勾选的地方不勾。

image.png 

image.png 

赋完力场之后就可以导出转成data文件了,导出格式为car,然后会在文件夹中得到两个文件Layer.carLayer.mdf

image.png 

lammps中自带的msi2lmp.execvff.frc复制到存放Layer.carLayer.mdf的文件夹,打开cmd,输入命令即可转化为data文件,如下图所示:

msi2lmp.exe Layer -class 1 -f \Users\Administrator\Desktop\wetting_on_cu/cvff  -i>info.data

image.png 

注意使用cvff力场文件的绝对路径,否则会报错。用OVITO打开生成的data文件,如下图所示,生成的data文件box上下边界与MS中不一致,可以在data文件中直接修改即可。

image.png 

接下来是修改势函数,直接从力场文件中导出的参数有时候会有错漏或者需要人为进行修改以适应研究目的。这里假设基底原子和水分子之间的相互作用由12-6 LJ势函数描述,水分子之间的相互作用由12-6 LJ势函数和库仑力来描述。并且润湿时基底固定不动,改变基底原子和水分子之间的势函数参数来调节基底的润湿性。因此基底原子与基底原子之间的相互作用参数可以随意设置,水分子这里假设使用spce模型,具体参数可从论文中获得。势函数参数可以写在data文件中,也可以写在in文件中。当势函数函数参数写在data文件里面时,在in文件中指定所需的不同原子间的相互作用的混合模式即可,在in文件中则可以直接指定不同原子间的相互作用,因此本文中的势函数参数均在in文件中。

编写in文件

编写in文件需要注意为模拟完成后求接触角提供便利,目前论文中求接触角主要有两种方式,第一是直接在轨迹文件中拟合液滴,第二是求得气液界面,然后用球或者圆拟合气液界面,得到接触角。无疑第二种方法是更适合的,求取气液界面需要获得水分子密度的分布,在lammps中,可以进行对区域进行圆柱形划分,也可以进行长方体划分。本例中采用的液滴是球形的,使用圆柱形划分会更省时间,但圆柱形划分要确定圆心的位置,如果圆心位置稍有差别,结果可能会截然不同,因此本文推荐使用更稳健的长方体划分。本文将data文件命名为wetting.datadata文件中有原子质量信息以及拓扑信息,如下图所示,势函数参数在in文件中。

 

in文件入所示:

# ----------------- Init Section -----------------

dimension               3

units                   real

boundary                p p f

atom_style              full

bond_style              harmonic

angle_style             harmonic

neighbor                2.0 bin

neigh_modify            every 1 delay 0 check yes

pair_style              lj/cut/coul/long 10 12

kspace_style            pppm 1e-4

kspace_modify           slab 3.0

timestep                1.0

# ----------------- Atom Definition Section -------

read_data            system.data

# ----------------- group define ------------------

group           water type 2 3

group           Cu type 1

# ----------------- Settings Section --------------

pair_coeff      1 1 0.0      0.0          

pair_coeff      1 2 0.449    2.752         

pair_coeff      1 3 0.0      0.0           

pair_coeff      2 2 0.1553   3.166         

pair_coeff      2 3 0.0      0.0           

pair_coeff      3 3 0.0      0.0           

bond_coeff      1 540.6336   1.0  

angle_coeff     1 50         109.47

# ----------------- Run Section -----------------

velocity        Cu set 0 0 0

fix             fix_Cu Cu setforce 0 0 0

#能量优化

dump            1 all atom 100 mini.xyz

minimize        1e-5 1e-7 1000 10000

#undump         1

#初始温度

velocity        water create 298 98763 rot yes dist gaussian units box

fix             fix_water water shake 0.0001 10 0 b 1 a 1

#计算温度

compute         newTemp water temp

#控温

fix             1 water nvt temp 298 298 100

fix_modify      1 temp newTemp

#计算Z方向密度分布

compute          z water chunk/atom bin/1d z 3.61 1 units box

fix              2 water ave/chunk 100 150 15000 z density/mass norm sample file density_Z.profile

dump             2 all custom 5000 dump.lammpstrj id mol type x y z

run              30000

 

data文件生成和in文件编写的思路大体如上所示,遇到问题可加QQ670911765讨论。值得注意的是,纳米尺度下的润湿,线张力的作用不可忽视,会使平衡接触角偏离宏观观测,因此建议本文建议如果润湿不涉及研究线张力的影响时,优先采用二维液滴,即圆柱状液滴。本文推荐两篇参考文献,一篇是中文学位论文,里面给出了润湿的in文件,可以作为参考,一篇是采用二维液滴研究基底与液滴相互作用对润湿性的影响,可为研究润湿的同学提供一些理论指导,推荐论文如下:

Effect of Solid−Liquid Interactions on Substrate Wettability and Dynamic Spreading of Nanodroplets: A Molecular Dynamics Study

https://sci-hub.se/10.1021/acs.jpcc.0c07919

固着纳米水滴润湿行为的分子动力学模拟

https://xueshu.baidu.com/usercenter/paper/show?paperid=652d7a7b1261e10f631ccaa7ca776d69&site=xueshu_se

使用下面MATLAB代码可以方便测量接触角

clear

rho=imread('C:\Users\baolu\Desktop\3.png');

imshow(rho)

set(gcf,'outerposition',get(0,'screensize'));

hold on

[p1,p2] = ginput(3);

x=[p1(1),p2(1)];

y=[p1(2),p2(2)];

linex=[x(1),y(1)];

liney=[x(2),y(2)];

plot(linex,liney,'r-','Linewidth',3);

hold on

z=[p1(3),p2(3)];

A=[x(1)-y(1),x(2)-y(2);z(1)-y(1),z(2)-y(2)];

B=[x(1)^2-y(1)^2+x(2)^2-y(2)^2;z(1)^2-y(1)^2+z(2)^2-y(2)^2];

 

ab=A\B;

a=ab(1)/2;

b=ab(2)/2;

circleCenter = [a,b];

 

c2 = (x(1)-a)^2+(x(2)-b)^2;

radius = sqrt(c2);

theta = (0:pi/360:2*pi)';

Result = zeros(size(theta,1),4);

for i = 1: size(theta,1)

    Result(i,1) = i;

    Result(i,2) = theta(i);

    Result(i,3) = circleCenter(1) + radius*cos(theta(i));

    Result(i,4) = circleCenter(2) + radius*sin(theta(i));

end

plot(Result(:,3),Result(:,4),'r-','Linewidth',3);

hold on

ca_left=90-asin((circleCenter(2)-x(2))/radius)/3.1415926*180;

ca_right=90-asin((circleCenter(2)-y(2))/radius)/3.1415926*180;

 

 

感谢鲍路瑶老师的分享,内容来自于鲍老师分享出来的资料

如有需要添加微信:lmp_zhushou  进入微信群,帮助他人,共建社区

获取完整版lammps讲义可以加微信lmp_zhushou或加入QQ994359511

该部分由武汉大学阿湖宝博士编写


最新评论
请先登录才能进行回复登录
技术支持: CLOUD | 管理登录
seo seo