【lammps模拟】水分子热导模拟代码分析

最新推荐文章于 2024-04-26 00:00:03 发布

元素魔方科研服务

最新推荐文章于 2024-04-26 00:00:03 发布

阅读量2k

点赞数

原文链接：https://mp.weixin.qq.com/s?__biz=Mzg5Njc1MzgwMg==&mid=2247486589&idx=1&sn=13f091fee2a2ea551bb63d88319314c1&chksm=c07d0598f70a8c8e39a7972fb32cc739116d22f0c8d228f81fad89cbb391a89ccda9a793d216&scene=126&&sessionid=0

版权

【9月开学季】

关注“元素魔方”小程序

更多优惠等你来领！

背景

分子动力学模拟是一种研究分子体系结构与性质的重要方法，已被广泛应用于生、化、环、材等学科领域。

LAMMPS分子动力学程序是一款免费的软件，可模拟气态、液态、固态以及混合态体系，其计算速度快，计算能力强，采用不同的力场合边界条件来模拟全原子，聚合物，生物，金属，粒状和粗粒化等。

水

水是研究最为广泛的一种流体，目前已有多种分子模型可供选择。水分子是极性分子，原子带有明显的极化电荷，这可以在Lennard－Jones势的基础上引入静电作用来描述。不过，由于水分子间可以形成氢键，加大了势函数描述氢键方向性和饱和性的难度。在不修改势函数表达式的基础上，可以引入只有电荷没有质量的虚原子来描述水分子的电荷分布，从而更好地描述氢键，这样不会使势函数过于复杂化。

lammps NEMD算法模拟水分子热导率

此案例输入的热量方式：fit ehex命令

这个命令可以直接设置热流大小，不需要再画热量曲线求斜率，因此，相对简单。

模拟的原理和石墨烯热导一样，在传热方向设置一个热源和一个冷源，在热源持续输入热量，在冷源持续抽走热量，由此产生温度梯度，模拟完成将热流和温度梯度等参数带入公式即可算出热导率。

Green–Kubo热导模拟代码详解

所用的模型如下图所示，沿z方向进行热传导：

in文件关键语句已经注释，供参考！

# 定义热流 variable F           equal 0.075# 步长             variable dt             equal 3.0  #总的模拟时间 variable tprod        equal 1000000#根据模拟时间换算成模拟步数 variable Nprod        equal floor(${tprod}/${dt})# SPC/E 参数 variable epsOO         equal 0.15535 variable sigOO         equal 3.166 variable theta         equal 109.47#定义截断半径 variable cutC          equal (xhi-xlo)/2. variable cutLJ    equal 11# 模型基本设置 units        real dimension    3 boundary     p p p atom_style   full #读取模型文件 read_data    "data.spce"# 原子分组 group    O type 2 group    H type 1 group    water type 1 2#力场参数设置 pair_style      lj/cut/coul/long ${cutLJ} ${cutC} pair_coeff      2 2 ${epsOO} ${sigOO} pair_coeff      1 2 0 0 pair_coeff      1 1 0 0 kspace_style   ewald 1.e-5 bond_style      harmonic               angle_style     harmonic               bond_coeff      1 1000.00 1.000 angle_coeff     1 100.0 ${theta}#近邻参数定义，因时间不长较大，因此设置一个较大的参数 neighbor  4.5 bin pair_modify tail yes#定义模拟结果采样频率参数 variable    Nsamp    equal 10 variable    Nrepeat  equal floor(${Nprod}/${Nsamp}) variable    Nevery   equal ${Nsamp}*${Nrepeat} # 计算盒子质心速度 (vcmx, vcmy, vcmz) variable vcmx equal "vcm(all,x)" variable vcmy equal "vcm(all,y)" variable vcmz equal "vcm(all,z)" variable vcm2 equal v_vcmx*v_vcmx+v_vcmy*v_vcmy+v_vcmz*v_vcmz#计算温度、压力、能量等 compute   cT  all temp compute   cP  all pressure thermo_temp compute   cPe all pe compute   cKe all ke variable  vE  equal c_cKe+c_cPe#设置分区坐标点，沿z方向热传导，只计算z方向关键坐标点 variable Lz          equal zhi-zlo variable delta       equal 4  variable dz          equal ${Lz}/60 variable zlo_Thi     equal -${Lz}/4.-${delta}/2. variable zhi_Thi     equal ${zlo_Thi}+${delta} variable zlo_Tlo     equal ${Lz}/4.-${delta}/2. variable zhi_Tlo     equal ${zlo_Tlo}+${delta}# 设置冷源和热源区域 region   Thi_region     block INF INF INF INF ${zlo_Thi} ${zhi_Thi} region   Tlo_region     block INF INF INF INF ${zlo_Tlo} ${zhi_Tlo}#计算冷源和热源温度 compute  cTlo   water temp/region Tlo_region compute  cThi   water temp/region Thi_region#校正温度 variable  Tlo_act    equal c_cTlo/2*3 variable  Thi_act    equal c_cThi/2*3#在热源持续输入热量，在冷源抽取相同的热量，以产生温差 fix fHi all ehex 1  ${F}   region Thi_region com constrain fix fLo all ehex 1 -${F}   region Tlo_region com constrain# nve系综 fix fNVE all nve#根据动能计算温度 variable  kB  equal  0.001987204 compute   ke  water  ke/atom variable  T   atom   c_ke/${kB}#固定键角 fix    fRattle all rattle 1e-10 400 0 b 1 a 1 reset_timestep  0 timestep        ${dt}#输出温度、能量等计算结果 thermo_style    custom step temp v_Tlo_act v_Thi_act ke pe etotal press v_vcm2 thermo_modify   flush yes thermo          1000#分块计算温度并输出，根据此数据可得温度梯度 compute         cchT  all  chunk/atom bin/1d z lower ${dz} fix             fchT  all  ave/chunk  ${Nsamp} ${Nrepeat} ${Nevery} cchT v_T file out.Tspce_ehex#输出能量 fix             fE    all  ave/time 10 500 5000 v_vE file out.Espce_ehex run             ${Nprod}