1 前言

亲水抗污材料因其出色的抗蛋白质吸附特性被广泛应用于船壳涂料^[1]、医疗器械^[2]等方面。该材料的抗蛋白质吸附主要来源于其表面水合层所形成的屏障^{[3, 4]}。根据其亲水机制，亲水抗污材料主要可以分为两大类：氢键型水合材料和离子型水合材料。

聚乙二醇(Poly(ethylene glycol)，PEG)是典型的氢键型水合的亲水抗污材料。它是由环氧乙烷开环聚合形成的链状聚合物，结构式为HO(CH₂CH₂O)_nH。聚乙二醇通常呈电中性，分子呈线性规整型螺旋结构或反式构型^[5]，亲水的-O-分布于聚合物链使得聚合物很容易与水结合。

大量对聚乙二醇的研究表明，聚乙二醇在水溶液中有着优良的低摩擦特性^[6]。Heeb等^[7]发现，金表面经过PEG5000修饰后，摩擦系数可以从1.00降低到0.03左右。Spencer课题组在对聚乙二醇水润滑方面的研究做了大量工作，主要研究了梳状聚合物聚L-赖氨酸接枝聚乙二醇(poly(L-lysine)-g-poly(ethylene glycol)，PLL-g-PEG)在各种表面的接枝技术和在各类液体中的摩擦特性^[8~10]。他们发现在聚二甲基硅氧烷表面覆盖PLL-g-PEG后，材料的边界润滑性能有显著提升，摩擦系数从1.0左右下降到0.02左右，润湿角也降至30°左右^[11]。另一方面，他们在表面覆盖有金属氧化物的摩擦计中加入纯缓冲液与含有PLL-g-PEG的缓冲液进行测量对比，发现与纯缓冲液相比，添加PLL-g-PEG后测得的摩擦系数从0.2下降到0.05^[12]。

亲水抗污材料的低摩擦特性与表面形成的水合层密切相关。Raviv等人对光滑硅表面在离子溶液中的摩擦研究表明，通过表面吸附离子、离子吸附水分子，进而在表面上形成具有较强的流动性的水层，可以起到润滑的作用^[13]。Leng等^[14]随后对相同体系进行的模拟也验证了硅表面间纳米级受限水层的流动性。Klein在后续对带电聚合物^[15]和聚合物电解质^[16]的水润滑摩擦的研究中延用了该机理模型对实验结果进行解释，均与实验结果符合得很好。Kyomoto等人通过将各类亲水材料(Oligo(ethylene glycol) methacrylate(OEGMA)，2-Methacryloyloxyethyl phosphorylcholine(MPC)等)接枝到交联聚乙烯表面^[17]，研究了它们在水溶液中的摩擦行为和耗损状态，从实验的角度证明了亲水材料表面的水合层与其在水溶液中的低摩擦特性之间的密切联系。然而，目前对于电中性亲水材料低摩擦特性的理论研究，尤其是从微观角度对亲水材料水润滑摩擦行为的分析还较少。传统实验法采用的原子力显微镜等分析仪器无法获取摩擦过程中水分子的动态行为，因此无法直接考察原子级的微观摩擦机理。分子动力学模拟方法可以直接考察微观原子的动态相互作用和摩擦过程，因此，它扮演着传统实验方法不可替代的角色。本文采用的分子动力学模拟方法的可靠性已经在之前的研究中得到了成功验证。以磷脂单分子膜在盐溶液中的摩擦为例，He等人^[18]采用分子模拟研究预测钠盐浓度升高会导致表面摩擦系数增大，这一结论随后得到了Klein等人^[19]的实验证实。

本文将通过对低聚乙二醇自组装单分子膜(Oligo (ethylene glycol)，OEG-SAM)在纯水中摩擦行为的分子动力学模拟，分析表面水合层的结构和动态行为对水润滑效果的影响。聚乙二醇是典型的氢键型亲水抗污材料，对它的研究可以加深对这一类材料水润滑特性的理解。选用自组装单分子膜则是因其易于在实验条件下制备，结构稳定。本文主要通过改变两表面的间距，借此改变两个水合层的相互作用，进而引起水合层结构和动态行为的变化，考察其对摩擦特性的影响。

2 模型建立和模拟方法

本模拟选择链末端带有四个聚乙二醇重复单元的OEG-SAM作为摩擦表面。表面的单链结构式为S(CH₂)₁₀(OCH₂CH₂)₄OH。以Au(111)作为基底构建的摩擦模型如图 1所示，部分水分子被限制在摩擦表面间作为摩擦介质，两侧各有水盒子与表面间的水相连，层间水可以与水盒子内的水分子进行自由交换，同时，水盒子中水分子采用恒压过程(1 atm)保持主体水(Bulk Water)的特征。上下表面沿x轴方向相对运动。由于模拟采用周期性边界条件，所以表面在x方向无限大。

模型使用CHARMM(version c30b1)软件(CHARMM软件是采用CHARMM力场的分子动力学模拟软件，适用于聚合物、溶液等体系的模拟)建立。模型建立步骤如下：首先构建摩擦表面。单个表面含有48个OEG聚合物链(x方向6个，y方向8个)，链的位置按照Au(111)的晶体结构排列，晶体结构为($\sqrt 3 \times \sqrt 3 $) R30°。聚合物链末端的S原子与基底的Au原子成键，S与S的相邻距离为0.496 nm，单个链所占基底面积为21.6Ų。链与表面法向量即z轴正方向夹角设为30°，链的聚乙二醇部分呈螺旋状构型，此时整个系统的能量最低。之后用相同方法构建第二个表面，并将两块表面如图 1所示平行放置，聚合物链端的羟基相对，两表面的中间位置为z轴零点。根据Raviv等人的实验研究^[13]报道，PEG水润滑时两个表面的间距可以低至8Å，因此表面间距(末端羟基的距离)分别设为8Å、10Å、15Å、25Å和50Å。将建成的结构以共轭梯度算法进行5000步的能量最小化，过程中硫原子的位置固定不变。在得到的表面结构中添加水分子，水分子模型为三点电荷模型(TIP3P)。水分子作为一个长方体盒子添加进系统，密度为1 g·cm^-3，x和y方向的长度分别为26Å和97.5Å，z方向的长度随表面间距变化。最终形成了由两层相向的SAM表面、层间水和两侧水盒子构成的水摩擦体系，再以相同方式进行能量最小化。将完成的结构升温至300K，并进行50 ps(1 ps=10^-12s)的平衡操作。完成所有上述步骤后，用于分子动力学模拟的初始结构建立完成，将初始结构导入LAMMPS软件中进行模拟计算。

本模拟采用的力场与Zheng^[20]等人在对OEG-SAM的研究中采用的力场参数相同。势能函数式包含键伸缩势能、Urey-Bradley角弯曲势能，键角弯曲势能、二面角扭曲势能、离平面振动势能和非键相互作用，包括范德华和库伦相互作用。完整的势能函数式表示为：

$\begin{array}{l} U = \sum\nolimits_{{\rm{bonds}}} {{k_{\rm{b}}}{{(b - {b_\mathit{\boldsymbol{0}}})}^2} + \sum\nolimits_{{\rm{UB}}} {{k_{{\rm{UB}}}}{{(r - {r_\mathit{\boldsymbol{0}}})}^2} + \sum\nolimits_{{\rm{angles}}} {{k_{\rm{ \mathsf{ θ} }}}{{(\theta - {\theta _\mathit{\boldsymbol{0}}})}^2}} } } \\ \quad + \sum\nolimits_{{\rm{dihedrals}}} {{k_\phi }[1 + \cos (n\phi - \delta )]} \\ \quad + \sum\nolimits_{{\rm{impropers}}} {{k_\chi }{{(\chi - {\chi _0})}^2} + \sum\nolimits_{{\rm{vDW}}} {{\varepsilon _{ij}}[{{(\frac{{{R_{{\rm{mi}}{{\rm{n}}_{ij}}}}}}{{{r_{ij}}}})}^{12}} - 2(} } \frac{{{R_{{\rm{mi}}{{\rm{n}}_{ij}}}}}}{{{r_{ij}}}}{)^6}] + \sum\nolimits_{{\rm{Coulombic}}} {\frac{{{q_i}{q_j}}}{{4\pi {\varepsilon _0}{r_{ij}}}}} \end{array}$

(1)

其中，等式右边的前五项是键的相互作用，k_b、k_UB、k_θ、k_φ和k_χ分别是各项键的相互作用常数，而b₀、r₀、θ₀、δ和χ₀则分别是各项的平衡值，后两项代表非键相互作用，范德华作用描述为12-6 Lennard-Jones(LJ)势能，其中ε_ij是势阱深度，R_min是LJ势能的最小值，交叉LJ势能参数ε_ij=(ε_iε_j)^1/2，R_minij=(R_mini+R_minj)/2是原子部分电荷，q_ij是原子i与j之间的距离。在本研究中，短程范德华作用开关函数的截断距离范围设置为10Å和12Å，长程静电作用采用particle-mesh Ewald^[21]方法计算。

在接下来的分子动力学模拟中，x轴和y轴方向采用周期性边界条件，z轴方向在上下表面末端硫原子所在平面处设置反射边界条件。使上下表面分别沿x轴的正负方向相对运动，运动速度为0.0001Å/fs(10 m·s^-1)，运动步长设为1.0 fs(1 fs=10^-15s)。采用Verlet算法对牛顿运动方程求解。为提高模拟效率，所有涉及到氢原子的键都通过RATTLE算法进行限制。整个系统采用NPT系综控温控压。模拟温度通过Nose-Hoover热浴法设定在300K，耦合时间常数为0.1 ps。压强通过Nose-Hoover算法保持为1 atm，压力在y方向进行调节。系统运行3 ns(1 ns=10^-9 s)后再运行1 ns记录数据用于分析。每隔1.0 ps记录一次系统状态，同时，每隔1fs记录表面x方向受力与表面z方向受力。所有的模拟均在Linux计算机集群上运行，分子动力学模拟软件为LAMMPS。

3 模拟结果与分析

本文采用的表面摩擦系数μ计算方法为：

$\mu = \frac{{{F_{\rm{x}}}}}{{{F_{\rm{z}}}}}$

(2)

其中F_x表示x方向的表面受力即表面摩擦力，F_z表示z方向的表面受力即表面法相压力。计算获得的不同表面间距下的OEG-SAM表面摩擦系数如图 2所示。

由图 2可见，OEG-SAM表面摩擦系数随表面间距的增大而减小。间距15Å时，摩擦系数为0.68。这与同属亲水抗污材料的磷酸胆碱自组装单分子膜的模拟摩擦系数相近^[18]。摩擦系数在表面间距15Å附近有较大变化，表面相距小于15Å时，摩擦系数随表面间距增大而下降的速率明显快于表面间距大于15Å时的下降速率。表面间距大于15Å后，表面摩擦系数受表面间距的影响显著减小。造成摩擦系数在15Å处存在拐点的主要原因是两表面水合层的相对位置关系。当表面间距低于15Å时，两表面的水合层直接接触，摩擦系数较高。当表面间距等于或大于15Å时，两表面的水合层之间存在类似主体水层中的水，摩擦系数低。此结论可以从亲水材料的水润滑摩擦与材料表面水合层中水分子的结构和动态行为关系推得，因此下面将进一步分析不同间距下表面间水分子的结构和动态行为。

3.1 表面间水的结构分析

对表面间水的结构分析包括分析表面间水的密度分布、计算水分子在表面上的径向分布函数以及统计表面与水分子形成的氢键数目。

(1) 表面间水的密度分布

本文采用的分析水的密度分布的方法是：将表面间的水在垂直于z轴方向分成若干1Å厚的水层，分别求得各层水的密度，并以模型中z轴的坐标为自变量，作出水分子在z轴方向的密度分布曲线。

如图 3所示，各个间距下表面间的水的密度分布曲线的线型都相似。水的密度随z轴坐标的分布主要分为三个区域：(1)表面内部，(2)表面附近，(3)两表面之间的主体部分。由于OEG具有强亲水特性，因此部分水分子会进入OEG表面区域，由图 3可见，这部分水的密度与表面间距没有明显关系，各个间距的表面附近都存在一层密度约为1.2 g·cm^-3的水合层。两表面中间的主体水密度约为1.0 g·cm^-3，这与实验条件下主体水的密度相近。

表面之间水在表面间距较小时有明显的水合层。在初始表面间距下，表面间均可观察到有两层较强的水合层。这些层的存在主要是由水合特征决定的，可以从下面的径向分布函数分析得到其存在原因。表面通过氢键结合表面邻近水分子形成水合层，水合层通过与邻近水分子作用，最终影响更多水分子的密度分布。在本研究中，单层水分子厚度约为3Å。链的运动会导致链的倾斜，降低链的高度，增大了表面间距，因此能够给水的分层提供更多的空间。

表面间水的密度分布曲线展现了表面附近水合层的存在以及表面间水的分层行为。这与Cummings等^[22]在对硅表面间受限水的密度分布分析结果相似。不同的是，在他们的研究中，单个硅表面上有明显的两层高密度水分子。这是因为其研究的硅表面为刚性表面，表面对近表面水分子结构的影响强于本文采用的柔性表面。此外，在其研究中表面是静止的，而在本文的研究中表面是运动的，这也会削弱表面对水分子结构的影响。

(2) 水分子在表面上的径向分布函数

为了分析水分子在表面附近的结构，模拟计算了水分子在表面上的径向分布函数(Radial Distribution Function, RDF)。其函数式为：

${g_{ij}} = \frac{{\left\langle {{N_{ij}}(r)/V(r)} \right\rangle }}{{{\rho _{j,{\rm{bulk}}}}}}$

(3)

函数描述的是以i粒子为参考粒子，计算j粒子的密度随i与j之间的距离变化而变化的情况。其中N_ij(r)是距离i粒子r距离的体积为V(r)的球壳内j粒子数的统计平均值。ρ_{j, bulk}是j粒子的体相密度。本文计算径向分布函数选择的两部分粒子分别为表面末端羟基的氧原子和表面间水分子的氧原子，截断距离为10Å。

从图 4中可以看出，各个间距表面间的水分子在表面上的径向分布函数曲线几乎重合，这表明表面间距对表面附近水合层结构的影响不明显。RDF曲线在3Å附近有一个峰值，表明在距离表面末端羟基3Å处水分子密度最高，这与水分子的密度分布函数的结果相一致。由此可知，表面附近第一层水合层在表面之上3-4Å内。曲线在第一个波峰之后再没有峰值较大的波峰。这与Zheng^[20]等人在对OEG-SAM抗蛋白质吸附模拟中获取的RDF曲线有一定差异。在该模拟中，水分子在OEG-SAM末端羟基氧原子周围的径向分布函数曲线有较为明显的第二个峰，位于5Å附近，g_ij值在1.2左右。本文中的RDF曲线在相同位置也有波峰，但峰值不够突出。这同样是由于本文中的OEG-SAM表面是运动的，在结构上较难形成稳定的第二层水合层。图 4还表明径向分布函数在7.5Å以后便趋于稳定，大于此值，水就表现出主体水的特征，因此，这也从一个侧面解释了摩擦系数在15Å附近存在拐点的原因。

(3) 表面与水分子形成的氢键数

OEG的亲水机制在于OEG链中重复单元的O原子和末端羟基中的O原子会与水分子形成氢键。因此，表面与水分子形成氢键的数目与表面对水合层的作用强度有着紧密的联系。

本研究对于氢键的形成采用结构定义：受体与供体间的距离小于0.35 nm且夹角小于60°。模拟中受体为OEG链中重复单元的氧原子和末端羟基中的氧原子，氢键的供体为水分子的O-H。计算获得的各个间距表面与水分子形成的氢键数目如表 1所示。结果表明，随着表面间距的变化，OEG-SAM表面与水分子形成的氢键数目基本保持不变，这说明在本研究所考察的表面间距范围内，从氢键角度考察的表面的亲水作用未受到显著影响。

3.2 表面间水分子的动态行为分析

表面间水分子的动态行为可通过水分子的平动行为和水分子的转动行为(偶极取向变化)进行分析，分析可采用停留时间分布函数和偶极自相关函数。由于表面间的水在结构上会出现分层，层与层之间的水分子的动态行为可能存在差异，所以对水分子的平动行为和转动行为的分析以分层进行考察。分层以表面末端羟基所处平面为起始点，3Å为厚度单位划分。

(1) 水分子的平动行为

摩擦过程中，水分子作为润滑介质会有各方向的平动行为，其中与摩擦密切相关的是水分子在z方向上的平动(垂直于表面方向)。通过z方向的平动，水分子得以接近或远离摩擦表面，由此带来表面结构与能量的变化，影响摩擦特性。为了分析水分子在z方向的平动行为，选择对OEG-SAM表面间的水分子在z方向做停留时间分布(C_r(t))计算。(C_r(t)表示的是水分子在t时刻出现在给定区域内的概率，其函数表达式^[23]为：

${C_{\rm{r}}}(t) = \frac{1}{N}\sum\nolimits_{{\rm{j = 1}}}^{\rm{N}} {\frac{{\left\langle {{P_{\rm{R}}}_j(0){P_{\rm{R}}}_j(t)} \right\rangle }}{{{{\left\langle {{P_{\rm{R}}}_j(0)} \right\rangle }^2}}}} $

(4)

其中，P_R为二元函数，N为体系内水分子个数，j原子在设定区域内时P_Rj=1，否则P_Rj=0。以初始时刻设定区域内的粒子为样本，一旦粒子移出区域，后续的函数值计算就不再考虑该粒子。初值C_r(0)=1，随着时间的变化，移出区域的粒子数增多，P_R=0的情况越多，C_r(t)逐渐减小。

本文首先分析的是OEG-SAM表面上的水分子在z方向平动行为强度的分层状况。如图 5所示，在表面间距50Å时，表面上的水分子按z方向的停留时间分布可以分为三层。前三层水分子的停留时间差异较大，而第四层与第三层水分子的停留时间分布曲线几乎重合。由此可以推断，在表面间距50Å时，运动表面上两层水分子的平动行为与主体水层中的水分子的平动行为相异。距离表面最近的一层水受到表面较强的影响，因此有较长的停留时间。介于第一层水层和主体水层之间有第二层水层，目前无法判定其主要受表面的控制还是仅作为水合层与主体水层间的过渡存在。各层中的水分子的平动行为强度大小依次为：第一层水层＞第二层水层＞主体水层。为了了解其他表面间距下的OEG-SAM表面之间的水是否也有同样的平动行为分层，尤其是在表面间距8Å和10Å时，这样的第二层水层是否还存在，本文计算了所有表面间距下表面间水分子的分层停留时间分布。

从图 6中可以看出，表面上的水分子在z方向平动行为方面的分层在各个表面间距下都存在，且分层的结果相似。OEG-SAM表面上的有第一层水层和第二层水层各一层。各层水分子的停留时间都符合：第一层水层＞第二层水层＞主体水层。在表面间距8Å和10Å时，尽管表面之间只有3层和4层水分子，还是形成了第二层水层，并因此缺少了主体水层。所以从受力方面考虑，第一层和第二层水层本质上都属于水合层，平动行为受到表面的影响。这与上述Zheng^[20]等人在对静止的OEG-SAM抗蛋白质吸附的模拟中发现OEG-SAM表面形成两层水合层相吻合。综上，在水分子的z方向的平动行为分布上，间距8Å和间距10Å的表面间的所有水层在平动行为方面都属于受表面作用的水合层，两者之间的差异则在于水合层的层数。而间距15Å以上的表面间才有自由的主体水层。

(2) 水分子的转动行为

摩擦过程中，摩擦介质粒子除了平动外还伴随着自身的转动，摩擦介质粒子的转动行为会带来能量的变化，同样也会影响表面的摩擦特性。分子的转动行为可以用偶极自相关函数表征。偶极自相关函数C_μ的函数表达式^[24]为：

${C_{\rm{ \mathit{ μ} }}}(t) = \left\langle {\overrightarrow {\mu (0)} \cdot \overrightarrow {\mu (t)} } \right\rangle = \left\langle {{\rm{cos}}\theta (t)} \right\rangle $

(5)

函数描述的是系统内目标粒子在t时刻的偶极$\overrightarrow {\mu (t)} $与初始时刻的偶极$\overrightarrow {\mu (0)} $的夹角余弦值θ(t)的均值。在使用中将偶极$\overrightarrow {\mu (t)} $化为单位向量，则函数初值为1。随着运动时间的增长，粒子偶极与初始时刻的偶极的夹角越来越大，导致总余弦值均值C_μ(t)越来越小。因此，偶极自相关函数随时间变化下降得越快，代表水分子在层内的转动行为越剧烈。对水分子偶极自相关函数的计算依然是先对间距50Å的表面之间的水分子分层进行。

由图 7可知，表面附近的两层水层中，其水分子的偶极自相关函数随时间的下降速率低于主体水层。这表明，与平动行为类似，表面上水分子的转动行为也有两层与主体水层不同。对水分子的平动行为进行的分层讨论在转动行为上依然有效。因此，同样可以将表面间的水分子在转动行为上也划分为第一层水层、第二层水层和主体水层，各层水分子转动行为强度的大小顺序是：第一层水层＞第二层水层＞主体水层。接下来对各个间距表面间的水分子的转动行为进行分层分析。

由图 8可知，水分子的转动行为在各个间距的表面间也有明显的分层，并且分层结构与平动行为的分层结构类似。单表面上第一层水层和第二层水层都只有一层。各个间距表面之间的各层水分子转动行为的强度排列依然是：第一层水层＞第二层水层＞主体水层。在表面间距8Å和间距10Å时，水分子的转动行为仅有第一层水层和第二层水层，而没有主体水层；在间距大于15Å时，出现了主体水层。因此在转动行为方面，第一层和第二层水层本质上都属于水合层，受到表面的作用。与水分子的平动行为稍有不同的是，水分子在第二层水层和主体水层中的转动行为与表面间距的对应关系比较明显。

4 结论

对自组装单分子膜表面(OEG-SAM)水润滑摩擦的分子动力学模拟表明，其表面间的摩擦主要取决于水合层的相互作用。OEG-SAM表面对近表面两层水分子的动态行为影响较为明显。表面间距小于15Å时，两表面水合层间距离过近并存在部分重叠，因此表面摩擦系数较高。表面间距大于15Å时，两表面水合层没有直接接触，表面之间除受表面影响的水分子层外还有不受表面影响的主体水层，该层中的水分子流动性强，可以起到较好的润滑效果，因此表面摩擦系数较低。