润湿-编写in文件

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

值得注意的是，纳米尺度下的润湿，线张力的作用不可忽视，会使平衡接触角偏离宏观观测，因此建议本文建议如果润湿不涉及研究线张力的影响时，优先采用二维液滴，即圆柱状液滴。

本文推荐两篇参考文献，一篇是中文学位论文，里面给出了润湿的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或加入QQ群994359511