LAMMPS分子动力学基于周期扰动法的黏度计算

关键词：黏度,周期扰动法，SPC/E水分子，分子动力学，lammps

目前分子动力学计算黏度主要有以下方法：（1）基于 Green - Kubo 关系的方法。从微观角度出发，利用压力张量自相关函数积分计算黏度。理论基础强，能考虑复杂微观因素，但计算量极大，对计算机性能和时间要求高，积分上限选择需谨慎。（2）爱因斯坦关系法。通过分析粒子扩散行为间接求黏度，依据爱因斯坦关系，由粒子扩散系数计算。计算相对简单，只需粒子运动轨迹信息。不过对于粒子间作用强、扩散不明显的复杂体系误差较大，粒子半径选择也会影响结果。（3）非平衡分子动力学方法。对系统施加外力场或速度梯度使其处于非平衡态，根据非平衡条件下的响应，如应力和速度梯度，依牛顿黏性定律计算。可模拟复杂流场中黏度，但施加外力和处理边界条件要注意，理论较复杂。

本文将介绍非平衡分子动力学计算黏度的方法之一：周期扰动法。周期扰动法计算黏度的原理主要基于对流体系统施加周期性的外力扰动，然后通过分析系统的响应来计算黏度。具体如下：

(1)施加周期扰动：在流体系统中引入一个周期性变化的外力场，比如在模拟体系的某个方向上按照一定的频率和振幅施加外力。这个外力会使流体中的粒子或分子产生相应的运动和位移。

(2)系统响应分析：流体系统在受到周期外力扰动后会产生响应。由于流体具有黏性，粒子或分子的运动不会立即跟上外力的变化，而是会有一定的滞后。这种滞后表现为流体内部产生剪切应力，并且剪切应力会随着时间周期性地变化。

(3)计算黏度：通过测量流体系统在周期扰动下的应力响应和外力的关系，结合相关的理论公式来计算黏度。在周期扰动法中，通常会测量系统在多个周期内的响应，并对数据进行统计分析，以获得准确的黏度值。

接下来我们以计算水的黏度为例子，用lammps中的周期扰动法计算其黏度。

图1.左图为水分子系统的初始结构，右图为经过1ps的时间扰动之后的结构

图2.log.lammps文件输出结果局部截图

图2中，最后一列为黏度的倒数。经过几百ps时间之后，平均值会趋于一定定值。图3展示的就是黏度的倒数平均值随时间的变化。我们发现经过400ps之后，平均值趋于一个定值，约为1480/Pa/s。求倒数之后，得到黏度大约为0.68mPa.s(计算温度为300K)。

图3。黏度的倒数平均值随时间的变化。

最后附上计算代码