分子模拟在化学工程中的应用

吕玲红，陆小华，刘维佳，朱育丹

（南京工业大学材料化学工程国家重点实验室，江苏 南京 210009）

近年来，随着实验仪器与劳动力成本的增加，人们做实验的花费也随之增长，据估计每过10年实验成本就要翻一番，而工业设计又离不开大量可靠的物性和热力学数据。分子模拟是利用分子理论与计算方法，对物理和化学过程中分子的微观行为进行仿真的一种技术。随着计算机硬件的不断发展，普通微机的性能不断提高，价格不断下降，采用计算机来模拟实际体系的成本也随之下降。研究者认为每一项新理论诞生都已经离不开分子模拟的工作[1]，因此，分子模拟技术愈来愈受到工业界和学术界的关注。多尺度的分子模拟涵盖了电子到分子团簇体系。近20年来，由于以量子和统计力学为基础的理论方法更加完善，分子之间力场的不断改进以及计算机运算速率和存储容量的不断增加，提供了更高效的算法，人们在分子层面上模拟物理和化学过程的能力得到了飞速的提高。除了可以提供基础认识以外，很多理论和模拟方法的可信度，使其足以应用到以下几个方面：预测非寻常条件下物质的性能，例如高温高压下的相态；研究在实验室条件下无法接触的剧毒化学和生物制剂；预测非理想条件下如纳米体系中物质的性质和过程；验证各种经验定律的局限性。对于验证各种经验定律的局限性而言，宏观定律在纳米级体系中的不适用是一个很好的例子，那么在纳米级尺寸下什么条件会使诸如菲克扩散定律、热力学第二定律，或者吉布斯表面热力学（表面张力、其他过剩表面性质、杨氏方程、拉普拉斯方程和开尔文方程等）出现不适用？要回答这些问题，采用分子模拟的方法要比直接实验的方法容易得多。在这里采用模拟方法还有很多现实原因，如实验方法耗时较长，花费较大，如果要达到令人满意的预测，可能需要巨大的经济投入。比较典型的一个例子就是近20年来稳定发展的多组分相平衡的预测，减少了大范围实验测定。

1 分子模拟的不同尺度

实验和模拟方法可以定义为以下4类：（1）电子尺度。这一尺度主要有基本粒子，如中子、质子和电子等构成；（2）原子尺度。物质主要包括原子，其行为遵循统计力学的定律；（3）介观尺度：主要包括原子团簇，每个团簇都还有大量的原子；（4）连续尺度。这一尺度中物质是连续的，许多熟知的宏观定律都可以适用，如连续性方程，能量守恒定律，热流傅里叶定律的连续定律等。其中分子模拟的方法主要涉及到前3类。

1.1 电子尺度

从头算方法是最精确而严格的，除了涉及到原子种类的信息，它不需要任何实验参数的输入[2,3]。当电子发生重排，或者其他基本粒子出现时，如在导体和半导体中、光激发装置中或者有化学键被破坏和产生的化学反应中，就需要这种计算原子尺度的方法。但是，电子结构的计算方法是计算的需要，会产生很多关于电子结构的细节，它的应用还只局限在原子核数量相对较少、时限比较短的计算。

1.2 原子尺度

原子级别的方法主要有蒙特卡洛（MC）模拟和分子动力学（MD）模拟[1,4]，它们共同之处在于两种方法都需要描述体系中分子之间相互作用的位能函数。在MC方法中，体系的微观构型通过从一定的几率分布取样得到，所需宏观性质通过系综平均获得。在MD方法中，所有的粒子先被赋予起始位置和动量，然后通过求解一组所有粒子的运动方程组，得到各粒子的运动轨迹，宏观性质通过时间平均获得。准爱高狄假设（Quasi-ergodichypothesis）假定两种方法是等价的，即沿着相轨道（沿时间）的积分均值等价于系综的均值。MC和MD方法可以进行上千甚至上百万的原子计算，时间为10 ns或者100 ns，这些都可以根据体系大小和复杂度建立。在利用这些方法时，会丢失电子方面的细节，但是对于大多数物理过程而言，电子的扰动非常小，这些就变得不重要了，大多数相转变和传递过程可以作为示例。

1.3 介观尺度模拟

在全原子经典处理中，分子在原子细节上得到模拟。这种处理虽然比较精确，但是对于大分子体系（聚合物和蛋白质等）或者时间超过1 ms的慢过程（例如自组装过程和慢扩散）却不再适用，对于这样的体系，各种介观尺度的方法得以发展，逐渐形成了一种“粗粒化”方法，但粗粒化会忽略掉原子尺度细节的问题。常用的粗粒化模拟方法有耗散粒子动力学模拟方法（dissipative particle dynamics，DPD）、格子玻耳兹曼方法（lattice Boltzmann，LB）、场论聚合物模拟方法（field-theoretic polymer simulation，FTPS）以及运用Martini力场的介观动力学等方法。此外，布朗分子动力学（Brownian dynamics），格子链蒙特卡罗（lattice Monte Carlo）和密度泛函（density functional theory）研究方法等，都属于介观尺度模拟的范畴。

1.4 反应力场方法

反应力场分子动力学（ReaxFF-MD）[5]模拟方法是近年来兴起的一种新的模拟方法，与量子化学方法相比，它既可以模拟化学反应过程，也可以模拟结构和传递等物理性质[6]。反应力场是基于计算化学发展起来的，它能够填补量子化学（quantum chemical）和经典力场（empirical force field）之间的空缺。ReaxFF反应力场是以键级和键能的概念来处理化学反应过程中键的连接关系的变化，以这个键的变化作为化学反应的中心，同时伴随传统动力学过程中粒子的运动。由于在化学反应过程中，没有采用量子化学的计算方法，可以模拟较大体系，因此，其在空间和时间尺度上要比Car-Parrinellomolecular dynamics（CPMD）的方法更具优势。与传统MD模拟相比，ReaxFF-MD不仅能够得到流体的物理性质，而且能模拟体相流体与表面物质的化学反应，其模拟过程更接近于真实情况。

2 分子模拟的技术关键

计算机是计算化学的硬件基础，并行计算机的发展更加推动了计算化学的应用，使得计算的体系不断加大，计算的时间大幅加快。在上世纪70年代末，并行计算作为一个概念被广泛传播。在1980年，Seeger发表了开创性的论文《小型计算机用于并行计算：量子化学的有力工具》[7]。

现在，全世界的化学家使用万亿次和千万亿次的并行体系，小到十来个，大到几十万个处理器来模拟。计算化学软件也随着并行计算机的发展，不断地改进算法和代码，来充分利用硬件。和相对简单的从标量到向量的转变相比，并行体系要复杂得多，矢量处理器只要求有更长的循环，而并行计算平台要重写软件来处理数据和任务在处理器间的通信。

在过去几年里，硬件的发展如多核处理器已经改变了细粒度并行的概念，再一次促使计算化学家重新思考怎样利用一些硬件的潜力。最近出现的图形处理单元（GPU）在计算化学程序中的使用再一次改变了这种情形。这两个新的进展显示了计算化学家对计算能力的永不满足，不会止步于现在的标准的并行方案。

3 分子模拟在化工中的应用

以共聚高分子材料、功能膜、纳米碳管和层柱状多孔材料等为研究对象，经典分子动力学模拟与蒙特卡罗模拟在对其介观层次的微相结构及其演变过程、界面结构以及小分子和胶体颗粒的吸附和传递进行研究中起到了重要的作用。通过分子模拟手段建立分子结构、材料的介观微相结构与材料的宏观性能之间的定性和定量关系，可以为复杂材料的合成、加工和应用提供理论指导和基础数据。

介观多孔材料中流体相行为和传递性质的研究，致力于多孔氧化硅、活性炭纤维、介孔活性炭微球、碳纳米管、层柱材料、多金属复合钛酸盐功能陶瓷材料及大孔分子筛MCM-41等多孔材料的应用，研究氢、甲烷和二氧化碳等气体在孔道中的相行为、储存密度和输运过程，为筛选和制备实用的储能材料提供理论指导。

模拟和实验技术的结合是受限流体行为研究的重要技术发展方向，然而两者实际“结合”存在一定的困难，模拟研究被认为是“观察”受限下分子集群行为的有效工具，由于微观实验观测手段的缺乏或不够精密，实验研究者总希望模拟技术能把实验现象全部真实地重复出来，在揭示所有的机制后，对材料进行重新设计改良。与此同时，模拟研究者总期望实验能给出详尽的微观信息，再加以抽象进行模拟分析，这造成了模拟和实验的对话困难。因此，有必要在两者“定量”关联之前，进行模拟和实验的“定性互动”研究，即将分子模拟得到的“定性”认识，直接在实验中实践，优化实验方向。

3.1 纳-微尺度和表界面的量化计算及受限条件下流体的分子模拟

从相图的计算到过程的模拟，溶液热力学作为基础研究在传统工业发展过程中发挥了巨大的作用。随着科技的进步，生物、能源和材料等学科迅速发展，使得溶液热力学所涉及的体系对象越来越复杂。与单一的流体不同，现在所涉及的体系往往具有大比表面纳米结构的固体参与。对于流体分子在这样的受限空间下表现出的行为的认知，有助于更好地建立预测型的数理模型，同时也有助于作出前瞻性的决策。模拟和实验的结合是建立理论的发展方向，然而两者的结合并不是那么容易，所以在模拟和实验“定量关联”之前，进行“定性互动”的研究，从两头出发逼近问题是必要的。以下通过3个实例阐述量化计算与分子模拟研究如何和实验研究相互定性互动的。

3.1.1 TiO2表面化学吸附现象

TiO2的表面化学吸附性质非常复杂，不同晶型的表面上，同一种晶型的不同表面上的吸附行为都不同。虽然不同晶型的TiO2具有相同的化学组成，它们的晶体结构却不尽相同，造成了不同的表面化学活性和小分子在其上的不同吸附行为。小分子在TiO2表面的吸附行为受表面配位态和表面结构宽松度的影响，还与表面上的空穴或化学基团相关，空穴和表面基团能促进，也能抑制流体分子的吸附，进而改变催化剂的反应活性。通过“认知”反应物分子在TiO2表面的吸附行为可以更好地“认知”材料，并进一步对TiO2催化剂材料的表面进行调控，进而根据需求设计TiO2催化剂材料。

Vittadini等[8]通过第一性原理计算方法研究脱硝的反应机理，这种反应机理需要反应物NH3分子吸附在表面的Lewis酸性位上或者通过羟基吸附在表面，据此推测宽松一点的表面可以使更多的反应物NH3分子吸附在表面，所以可能更适应选择性催化还原（SCR）反应。Czekaj等[9,10]报道了在TiO2催化剂降解异氰酸的过程中，NO2分子可能会引起TiO2分子的失活，因为NO2分子吸附在TiO2表面形成NO3－，阻止了下一步反应物分子的吸附。催化剂表面反应物不同的吸附状态也会导致不同的反应路径。比如SCR反应中，NH3分子吸附在催化剂表面是最关键的一步，NH3分子吸附在催化剂表面Lewis酸性位和NH3分子吸附在催化剂表面Brönsted酸性的表面羟基－OH作用两种方式会引起不同的反应路径[11]。另外，催化剂表面羟基等基团也可能成为活性中心而改变反应物的吸附形态，并进一步影响催化反应的路径和性能。比如，Maira等[12]的研究结果表明甲苯吸附在表面孤立羟基和氢键羟基上能够引起羟基消失，使催化活性下降，且孤立羟基与甲苯作用生成苯乙醛，但氢键羟基与甲苯作用则生成CO2和H2O。

Guo等[13]通过第一性原理计算的方法考察了TiO2表面配位性以及表面结构宽松性对TiO2-B（100）表面吸附行为的影响。NH3分子在TiO2-B（100）表面的吸附能要大于其在TiO2-B（001）表面的吸附能，这意味着相对于NH3分子，TiO2-B（100）表面要比TiO2-B（001）表面更活泼。TiO2-B（100）表面比TiO2-B（001）表面能够接受NH3分子更多的电子造成了TiO2-B（100）表面更强的酸性，表面配位的不同带来了表面不同的电子分配，并进而影响了小分子在TiO2表面的吸附行为，随着TiO2-B（100）表面NH3覆盖率的增加，分子吸附方式的平均吸附能单调下降，与覆盖率成较好的线性关系，其函数关系式是y=111.00-36.29x ，但是二聚体吸附方式的平均吸附能并没有太大的变化。由于TiO2-B（100）表面特殊的宽松结构，NH3在表面上平均吸附能的下降率要远小于rutile（011）表面。尽管差分结合能随着NH3覆盖率的增加而显著下降，但二聚体形成能随着NH3覆盖率的增加只有少量的下降。在1/3 ML覆盖率条件下的差分结合能大约是1 ML覆盖率下的2倍。二聚体形成能的最大值与最小值之间的差值为4.82 kJ/mol。在5/6覆盖率的条件下，NH3分子吸附在TiO2-B（100）表面Ti5c位与NH3分子和已吸附在TiO2-B（100）表面的NH3形成二聚体的概率差不多相等。TiO2表面结构的不同会带来表面吸附分子之间作用力的不同，进而影响分子在高覆盖率下的吸附，带来TiO2表面吸附能力的不同，可能会影响SCR反应的效率。

3.1.2 碳管内壁修饰基团对受限水分子结构的影响

纳米受限下水分子的行为（比如，生物分子的空腔中[14]，细胞膜的水通道中[15]，以及纳米多孔催化剂的孔道中[16]）受到了很大的关注。理论和实验研究[17-19]表明，在受限条件下水分子往往表现出异于体相的性质，这些特殊的性质会影响到生物大分子[20]和纳米材料的功能[21]。许多分子模拟研究采用碳纳米管（CNT）作为圆柱状孔道模型来研究受限作用对于水分子结构的影响，一些特殊的受限水分子结构相继被报道[22-24]。

然而，CNT内壁具有均一的疏水特征，而真实的孔道往往呈现化学非均一性，不能完全反映出真实孔道的一些特征。比如，在蛋白质孔道表面有许多极性和非极性残基分布[25]；许多纳米多孔氧化物由于壁面被羟基化后呈现部分的亲水特征[26]；在疏水的多孔碳表面也分布有亲水性的氧化活性点[27]。这些受限界面上的官能团被认为是导致纳米尺度下许多奇异现象的原因，许多理论工作开始研究这些基团对水分子性质的影响[28,29]，发现水通道Aquaporin-1的孔道结构的中部比较狭窄，而在细胞膜表面具有比较大的开口，水链的取向会受到靠近中部狭窄区域处的残基影响。

为了能给设计仿生的纳米器件提供指导，将-COOH基团修饰在CNT的中部区域来模拟Aquaporin-1生物水通道的形状。利用MD的模拟方法研究修饰基团对水分子静态结构的影响。尽管改性后CNT的尺寸比Aquaporin-1生物水通道大，但是由于改性CNT中部修饰了-COOH基团，导致中部实际孔道空间变小，所以形状上类似Aquaporin-1生物水通道。考察了管内和管口区域的受限水分子取向和空间分布[30]，结果表明，在CNT两侧未修饰的区域（S1和S2区域），水分子具有相反且稳定的择优取向，偶极指向CNT中部的-COOH基团。尽管管口与修饰基团有一定的距离，但是管口的水分子的偶极也受到影响，呈现择优取向。这种诱导作用随着-COOH修饰层数的增加而增大，管径的减小而增强。此外，中部内壁的-COOH对管口水分子的分布影响不大，但对管内水分子的分布影响较大，-COOH基团与管口的距离会引起不同的轴向分布。

综上所述，在CNT中部修饰亲水性基团使得CNT内以及管口处的水分子氢键结构重新发生了调整，促使基团两端的水分子团簇形成较为规整的双偶极的“赝相”结构，在认识受限水分子结构的基础上，进一步通过深入的动态性质的分析，认识结构与水分子动态传输性质之间的关系，这有助于认识生物通道高通量形成的原因，同时也为CNT膜材料以及仿生的纳米流体控制器件的设计提供一些思路。

3.1.3 碳覆盖二氧化钛狭缝模型中水分子的行为

对于工业应用而言，TiO2材料往往需要大比表面积。提高材料的比表面积的方法有两种：一是材料的纳米化；二是材料具有丰富的孔道结构。后者因其具有纳米/微米二级结构越来越受到关注。He等合成了一种具有介孔结构的TiO2材料[31]，并将其应用于光催化[32]、电化学[33]、药物担载[34]和工业催化[35]等领域。

随着我国风电的不断发展,风电装机容量不断增加,风电并网比例也极大提升。图1为2010—2017年我国的风电装机容量。图2为我国“十三五”电力系统规划中非化石发电量的占比。风电并网将会给电力系统带来两方面的影响,一是风电机组将会替代一部分传统化石能源机组出力,减少化石能源的消耗以及污染物的排放,符合我国“十三五”电力系统规划的要求;二是为解决风电并网后所带来的风电消纳问题,系统将配置更多的旋转备用去消纳风电,不仅会导致系统成本上升,而且在某些极端情况下,风电解列还会造成大量的功率缺额,使整个电力系统面临风险。因此,对风电备用问题的研究具有极为重要的现实意义和工程意义。

由于固液界面和纳米孔道的协同作用，TiO2纳米孔内流体的行为值得我们注意。TiO2表面研究和纳米孔道内流体行为的研究都各自展开，并十分成熟。Gelb等[17]认为纳米孔道上化学性质的变化会造成流体分子吸附的异常。Romero-Vargas Castrillón等[36]在水分子与化学非均一表面相互作用方面进行了长期的研究，他们通过表面极性对水分子行为影响的研究得出，表面的极性将控制着水分子与表面的接触过程。

Wei等[37]分别建立了金红石型二氧化钛（110）面（rutile TiO2（110））和石墨型碳（0001）面（graphite carbon（0001））组成的狭缝模型，用MD的方法对比研究了水分子在其中的结构性质、吸附和流动性等行为，发现表面的化学性质对水分子的吸附和流动的影响显著。具体来说，水分子在碳狭缝中的流动性要远远高于二氧化钛的狭缝，而水分子在二氧化钛表面的吸附要强于碳表面。同时，二氧化钛表面的水分子强吸附，通过水分子的氢键结构影响到水分子的流动性。在某些工业应用中，水分子在TiO2表面的强吸附往往会起到负面作用，如将TiO2材料用于电池电极[33]。

Wei等[38]建立了具有不同覆盖率情况下碳覆盖的TiO2的狭缝模型，并沿用水分子作为流体模型，研究其在该种具有非均一表面的狭缝中表现的行为（图1）。结果表明，碳的覆盖仅仅是阻碍了水分子在覆盖区域的吸附，它并不能对水分子与TiO2暴露区域的相互作用、强吸附和低流动性产生影响。但是碳覆盖确实增加了水分子在这些狭缝中的流动性。随着碳覆盖率的上升，流动性最多增加了200%。造成这一结果的原因是碳的覆盖打破了水分子与TiO2间形成的氢键网络[38]。

图1 碳覆盖的TiO2的狭缝模型[38]Fig.1 TiO2 slit model with covering of graphite slice[38]

上述结论对于碳覆盖TiO2材料的实验研究有一定的指导，但研究中的碳覆盖是理想化的，与实验相比，其分散性极高，实验中往往难以达到这样的分散度。如果实验中想避免水分子与TiO2的接触，最好能够尽量多的覆盖碳层，仅少部分暴露的TiO2就会通过氢键网络，明显降低水分子的流动性。另外，孔道大小对于水分子的流动性也存在影响，尽管它并不是决定因素。但是如果孔道过于狭窄，必定会阻碍水分子的流动。为了让流体能够有更多交换的机会，孔道的大小至少是流体分子大小的3倍以上。

通过上述研究实例，可以看出模拟和实验的“定性互动”研究可以在一定程度上帮助指导材料的制备和应用，然而，这对于建立受限流体传递理论的建立是远远不够的，需要在“定性互动”研究的基础之上，进行“定量关联”。利用密度泛函理论方法实现模拟和实验结果的定量关联，进而建立数理模型，是受限流体未来发展的重要研究手段。

3.2 介观结构受限条件下流体粗粒化模拟

复杂材料涉及非常广泛的应用领域，研究其介观结构的形成及受限于其中的流体介质的行为，有重要的实际意义。对于介观结构的形成及其动力学问题的研究，大多由胶体与界面化学和高分子学者完成，而计算机模拟和理论研究成果则多是由统计力学和凝聚态物理学家进行的，这是一个具有巨大挑战性的多学科交叉领域。在化学化工领域，主要以共聚高分子材料、功能膜和多孔材料等为研究重点，对其介观层次的微相分离结构及其演变过程、受限环境的影响、界面结构以及小分子和胶体颗粒的吸附和传递进行研究，建立分子结构、材料的微相分离结构与材料的宏观性能之间的定性和定量关系，为复杂材料的合成、加工和应用提供理论指导和基础数据。目前，对各个层次的规律的研究已经取得了长足的进步，但对不同尺度间关系的研究还不够，特别是不同尺度性质的相互推算还有待解决。不同领域的研究者普遍认为这是一种多尺度的研究，必须采用多尺度研究方法。除了实验研究外，计算机分子模拟、统计力学理论、密度泛函理论和粗粒化方法等的相互结合是取得成功的关键。目前比较现实的做法是先研究相邻尺度的联系以及预测方法，在取得相邻尺度结构的关系后，可以进一步构建多尺度的理论和模拟方法。

3.2.1 介观结构形成和演变的模拟

介孔材料孔径介于微孔与大孔之间，具有大表面积和三维孔道结构。有序介孔碳材料是近年来迅速发展起来的一类新型介孔材料。其不仅具有很好的疏水性、耐酸、耐碱、耐热性能和化学稳定性、较高的机械强度和良好的导电性，还可在微米尺度内保持高度的孔道有序性。在催化、吸附、分离、生物和电导材料等领域有巨大的应用潜力。Wang等[39]采用耗散粒子动力学模拟方法对有序介孔碳制备机理进行了探索。借助DPD方法研究了亲水端和疏水端组成不同的PEO-PPO-PEO型三嵌段聚合物在乙醇溶液中总体的自组装规律。研究了温度、剪切速率和体系各物质间不同体积比对F127在乙醇溶液中自组装结构的影响，并结合实验考察了F127和酚醛树脂在乙醇溶液中的自组装行为，从模拟的微观角度简单探讨了有序介孔碳的制备机理，为软模板法制备介孔碳提供了介观尺度上的理解和认识，研究了F127/酚醛树脂/乙醇体系下有序介孔碳致孔机制和制备机理。

图2 小纳米颗粒体系形成的不同自组装结构的快照[41]Fig.2 Snapshots of the self-assembled structures from the system[41](a)-hollow sphere; (b)-porous structure; (c)-perforated lamella with star polymers through nanospheres; (d)-hexagonally packed cylinders of nanospheres. The nanospheres are colored red, and the star polymers are colored green

3.2.2 纳微受限条件下流体性质的粗粒化模拟

界面粗糙度是影响纳微尺度催化反应的重要因素。其对于催化反应活性的影响主要由于其界面粗糙度不同，引起其上担载的分子和团簇的粘附力和稳定性不同。对于传递的影响则在于粗糙度对孔道内壁流动流体产生的摩擦力不同，从而影响流体的传递过程。现有的MD方法可以考察壁面化学非均一性对受限流体作用的影响，但用于研究表面粗糙度时，由于体系大小的限制显得捉襟见肘（如模型化纳米尺度凹凸表面），采用粗粒化模拟方法是解决该问题的可行途径。

已有的粗粒化模拟研究中对于模拟力场的构建存在多种方式。应用最广的是将一组原子或联合原子基团映射成一个粗粒化珠子[42,43]。利用全原子或者联合原子模拟结果或实验数据拟合得到粗粒化珠子之间的相互作用。这种粗粒化策略得到的特定粗粒化模型是基于关键的化学结构信息的，所以用于对特定体系的定量（或半定量）预测。原则上来说，粗粒化力场与相互作用势函数都要能够再现在联合原子模拟中观察到的结构性质，但有时也要求能够再现动力学与热力学性质[44]。由于粗粒化力场参数化的过程中只考虑了少量的结构信息，因此标准Martini力场参数不具有很好的普适性。针对特定体系，需要以实验数据和全原子MD模拟结果作为参考基准对粗粒化力场中针对催化剂表界面结构行为的模拟参数进行调试，得到更具针对性的力场对实验设计提供理论预测与指导。

对于受限条件下小分子的粗粒化主要是把小分子当成一个珠子（bead），且不考虑分子的极性。Mashayak等[45]采用粗粒化分子动力学模拟（CG-MD）和基于半连续理论的经验势（EQT）预测受限于石墨烯狭缝状管道中水的结构，水分子用单点各向同性粗粒化势能描述，势能参数参考全原子SPC/E水MD（AA-MD）模拟的结果得到。由受限水密度曲线的结果来看，CG-MD和EQT的结果与AA-MD模拟一致。

Maibaum[46]提出受限于疏水性纳米管的水的格子模型。利用分析方法和计算机模拟，发现管充满和全空条件在平衡下达到能量相等的条件。进一步发现水分子流以一定比例以简单随机规则通过充满管，该现象能用水脉冲传导解释。这两个结果与Hummer等[22]原子模拟的结果一致。最后分析了充满和全空管的转化，发现管入口处密度涨落在过程中起决定速率的作用。

薄液膜润滑对先进技术装置摩擦可靠性很关键，如微电子机械系统和硬盘驱动。然而由于膜结构和分子尺度的运动认识有限，薄膜的粘附和摩擦性能以及隐藏的机理不清。Fukuda[47]研究纳米厚度受限于两固体表面间perfluoropolyether（PFPE）膜的粘附行为，用粗粒化分子动力学模拟研究膜厚作为参数。与实验结果一致，模拟结构显示由PFPE膜产生的粘附力随着固体-固体空间的增加，先到达一个极大值，然后下降，小于4 nm厚的膜粘附力极大值增加非常显著，当呈现粘附力极大值时，膜在固体-固体空间轻度张开，比原始厚度略厚，因而密度比平衡密度稍低。

对于流体等软物质来说，已有比较完备的粗粒化方法和力场利用。为了可描述纳米微米介尺度结构的粗糙度，必须对固体的表面也进行粗粒化。固体粗粒化模型的建立目前还处于初步阶段，Steele[48,49]建立了气体小分子在点阵平面上的力场形式，这实际是对平整表面力场势能函数的一种粗粒化。研究者采用该力场进行了粗粒化模拟[50]。Zhao[51]在此基础上更考虑了石墨片层体系的偶极和四极的影响。上述处理只是针对平滑的无限延伸表面体系，因此，无法反映介尺度结构的粗糙度。排列棒类狭缝（RSP）模型被用于研究介孔碳的吸附[52,53]孔碳处理为同轴石墨片层的圆柱，再将这些圆柱排列为一定宽度的狭缝。吸附气体分子与该模型间的作用力场，则以气体分子与碳柱作用力之和来得到。虽然这种粗粒处理方法可以准确描述介孔碳对气体分子的吸附并与实验结果一致，但这种高度模型化的处理仍然不能用来描述界面粗糙度。如何将真实粗糙表面进行单元划分，建立介尺度结构粗粒化模型，是今后需要考虑的研究方向之一。

3.3 反应力场模拟

传统的MD模拟可以考察受限条件下流体的性质，也可以用来研究诸如表面亲疏水性造成水的停留时间等表面现象，得到不同表面的流体传递相关的系数，为实验设计提供可靠的依据。量子化学计算能够考虑不同分子在孔表面的反应活性，这种反应活性在一定程度上也会影响到流体的传递。MD模拟虽然能够反应孔内流体的传递性质，但是由于其采用的力场形式仅仅局限于粒子之间的物理作用，它无法体现粒子之间的化学相互作用，无法模拟化学反应的进行。量子化学计算由于计算量耗费大，使得量子化学无法应用于大体系中（如超过100个原子的体系）。以TiO2作为催化剂活性载体为研究对象，研究其催化反应中流体的传递性质和化学反应，利用MD模拟和量子化学计算只能考虑过程中的一部分性质，因此无法仅凭MD模拟或仅凭量子化学计算来描述TiO2作为催化剂活性载体的体系的整个传递和反应过程，势必要找到一种兼具MD模拟和量子化学计算特性的模拟方法。反应力场方法为研究此类问题提供了可能。

2001年，van Duin等[5]为执行大尺度反应体系（1 000个以上原子）的分子动力学模拟，开发了反应力场（ReaxFF）方法。它一方面将键长与键级联系起来，另一方面将键能和键级也联系起来，并由此模拟原子间键的分解。除非键作用外，分子内各部分能量均通过键级来表达。此外，ReaxFF包含库仑力和Morse势（范德华力），可以描述原子间非键作用。体系总能量由下式表示：

首个可以应用的ReaxFF力场是针对烃类开发的，研究的结果表明，计算结果的精度可与半经验量化计算方法-参数化模型法3（PM3）得到的结果媲美，且计算速率提高了100倍。经过10多年的发展，ReaxFF已经可以应用于各种催化和反应过程的模拟[54,55]，同时对表面性质也进行了研究[56,57]，对于反应中物理性质的关注还有赖于分析方法的完善和体系的进一步扩大。目前ReaxFF的主要问题是力场的普适性和完备性，几乎每研究一个新的体系，都需要对该体系的力场先进行拟合[58,59]，而这一过程的工作量很大。

许多研究者利用这一新的模拟技术进行了研究。Raju等[60]考察了水在300 K时，不同TiO2晶型的不同表面上的吸附和解离，结果表明，ReaxFF方法预测的分子解离吸附构型与实验结果一致。此外，分别分析了吸附水分子和吸附层外水分子氢键的强度，并尝试将ReaxFF与GCMC方法进行结合。Senftle等[61]采用基于ReaxFF混合GCMC/MD的方法对于氧化物热力学稳定性进行了研究，得到了Pd935团簇在温度300～1 300 K，氧气压力1.013×10-9～1.013×105Pa下的氧化理论相图。

反应力场模拟方法是对现有模拟方法的重要革新，在此之前的模拟方法，是对化学反应的性质和物理性质在不同尺度上加以分割，分别进行研究的，而ReaxFF则打破了这一限制，实现了物理性质和化学性质研究的统一。虽然现有的模拟体系对于真实体系的应用范围还较小，要得到体系结构传递性质，计算的速率还不能令人满意，但随着计算机计算速率进一步提高，并行技术的发展和方法的完备，反应力场模拟方法必将在化工应用中占据越来越重要的位置。

4 展 望

分子模拟技术发展至今，已取得了很大的进步。并行技术的提升使得可研究的体系越来越大，与实验结果的对比度也越来越高。多尺度模拟的展开也拓展了分子模拟技术可仿真的体系范围。当前的模拟结果已由定性向定量转化，并且可以弥补实验检测手段的不足，从分子的角度诠释物理化学现象的本质。因此，对化工过程中的关键问题进行多尺度的模拟是未来大力发展的方向，分子模拟必将在化工过程和设计中发挥越来越大的作用。

[1] Frenkel D, Smit B. Understanding molecular simulation: from algorithms to applications[M]. San Diego: Academic Press, 2002: 1-6,23-105.

[2] Cramer C J. Essentials of computational chemistry: theories and models[M]. Chichester: Wiley, 2004: 165-166.

[3] Jensen F. Introduction to computational chemistry[M]. Chichester: Wiley, 2007: 14-21.

[4] Allen M P, Tildesley D J. Computer simulation of liquids[M]. Oxford: Clarendon Press, 1987: 71-139.

[5] Van Duin A C T, Dasgupta S, Lorant F, et al. ReaxFF: a reactive force field for hydrocarbons[J]. The Journal of Physical Chemistry A,2001, 105(41): 9396-9409.

[6] Zhang L Z, Van Duin A C T, Zybin S V, et al. Thermal decomposition of hydrazines from reactive dynamics using the ReaxFF reactive force field[J]. The Journal of Physical Chemistry B, 2009, 113: 10770-10778.

[7] Seeger R. Parallel processing on minicomputers: a powerful tool for quantum chemistry[J]. Journal of Computational Chemistry, 1981,2(2): 168-176.

[8] Vittadini A, Casarin M, Selloni A. First principles studies of vanadia-titania monolayer catalysts: mechanisms of NO selective reduction[J]. The Journal of Physical Chemistry B, 2005, 109(5): 1652-1655.

[9] Czekaj I, Piazzesi G, Kröcher O, et al. DFT modeling of the hydrolysis of isocyanic acid over the TiO2anatase (001) surface: adsorption of HNCO species[J]. Surface Science, 2006, 600(24): 5158-5167.

[10] Piazzesi G, Elsener M, Kröcher O, et al. Influence of NO2on the hydrolysis of isocyanic acid over TiO2[J]. Applied Catalysis B:Environmental, 2006. 65(3/4): 169-174.

[11] Yuan R M, Fu G, Xu X, et al. Bronsted-NH4+mechanism versus nitrite mechanism: new insight into the selective catalytic reduction of NO by NH3[J]. Physical Chemistry Chemical Physics, 2011, 13(2): 453-460.

[12] Maira A J, Coronado J M, Augugliaro V, et al. Fourier transform infrared study of the performance of nanostructured TiO2particles for the photocatalytic oxidation of gaseous toluene[J]. Journal of Catalysis, 2001, 202(2): 413-420.

[13] Guo X J, Liu W J, Fang W, et al. DFT study of coverage-depended adsorption of NH3on TiO2-B(100) surface[J]. Physical Chemistry Chemical Physics, 2012, 14(48): 16618-16625.

[14] Pan Y S, Birkedal H, Pattison P, et al. Molecular dynamics study of tryptophylglycine: a dipeptide nanotube with confined water[J]. The Journal of Physical Chemistry B, 2004, 108(20): 6458-6466.

[15] Borgnia M, Nielsen S, Engel A, et al. Cellular and molecular biology of the aquaporin water channels[J]. Annual Review of Biochemistry,1999, 68: 425-458.

[16] Desbiens N, Demachy I, Fuchs A H, et al. Water condensation in hydrophobic nanopores[J]. Angewandte Chemie-International Edition,2005, 44(33): 5310-5313.

[17] Gelb L D, Gubbins K E, Radhakrishnan R, et al. Phase separation in confined systems[J]. Reports on Progress in Physics, 1999, 62(12):1573-1659.

[18] Mamontov E, Burnham C J, Chen S H, et al. Dynamics of water confined in single- and double-wall carbon nanotubes[J]. The Journal of Chemical Physics, 2006, 124(19): 194703.

[19] Byl O, Liu J C, Wang Y, et al. Unusual hydrogen bonding in water-filled carbon nanotubes[J]. Journal of The American Chemical Society,2006, 128(37): 12090-12097.

[20] Pratt L R, Pohorille A. Hydrophobic effects and modeling of biophysical aqueous solution interfaces[J]. Chem Rev, 2002, 102(8):2671-2691.

[21] Ohba T, Kanoh H, Kaneko K. Affinity transformation from hydrophilicity to hydrophobicity of water molecules on the basis of adsorption of water in graphitic nanopores[J]. Journal of the American Chemical Society, 2004, 126(5): 1560-1562.

[22] Hummer G, Rasaiah J C, Noworyta J P. Water conduction through the hydrophobic channel of a carbon nanotube[J]. Nature, 2001,414(6860): 188-190.

[23] Wang J, Zhu Y, Zhou J, et al. Diameter and helicity effects on static properties of water molecules confined in carbon nanotubes[J].Physical Chemistry Chemical Physics, 2004, 6(4): 829-835.

[24] Shiomi J, Kimura T, Maruyama S. Molecular dynamics of ice-nanotube formation inside carbon nanotubes[J]. The Journal of Physical Chemistry C, 2007, 111(33): 12188-12193.

[25] Zhou R H, Huang X H, Margulis C J, et al. Hydrophobic collapse in multidomain protein folding[J]. Science, 2004, 305(5690):1605-1609.

[26] Giovambattista N, Debenedetti P G, Rossky P J. Hydration behavior under confinement by nanoscale surfaces with patterned hydrophobicity and hydrophilicity[J]. The Journal of Physical Chemistry C, 2007, 111(3): 1323-1332.

[27] Kimura T, Kanoh H, Kanda T, et al. Cluster-associated filling of water in hydrophobic carbon micropores[J]. The Journal of Physical Chemistry B, 2004, 108(37): 14043-14048.

[28] Murata K, Mitsuoka K, Hirai T, et al. Structural determinants of water permeation through aquaporin-1[J]. Nature, 2000, 407(6804):599-605.

[29] Tajkhorshid E, Nollert P, Jensen M O, et al. Control of the selectivity of the aquaporin water channel family by global orientational tuning[J]. Science, 2002, 296(5567): 525-530.

[30] Zhu Y D, Wei M J, Shao Q, et al. Molecular dynamics study of pore inner wall modification effect in structure of water molecules confined in single-walled carbon nanotubes[J]. The Journal of Physical Chemistry C, 2009, 113(3): 882-889.

[31] He M, Lu X H, Feng X, et al. A simple approach to mesoporous fibrous titania from potassium dititanate[J]. Chemical Communications,2004, (19): 2202-2203.

[32] Li W, Liu C, Zhou Y X, et al. Enhanced photocatalytic activity in anatase/TiO2(B) core-shell nanofiber[J]. The Journal of Physical Chemistry C, 2008, 112(51): 20539-20545.

[33] Lu L H, Zhu Y D, Li F J, et al. Carbon titania mesoporous composite whisker as stable supercapacitor electrode material[J]. Journal of Materials Chemistry, 2010, 20(36): 7645-7651.

[34] Li Q, Wang X, Lu X, et al. The incorporation of daunorubicin in cancer cells through the use of titanium dioxide whiskers[J].Biomaterials, 2009, 30(27): 4708-4715.

[35] Zhu Y D, Li W, Zhou Y X, et al. Low-temperature CO oxidation of gold catalysts loaded on mesoporous TiO2whisker derived from potassium dititanate[J]. Catalysis Letters, 2009, 127(3): 406-410.

[36] Romero-Vargas Castrillón S, Giovambattista N S, Aksay I A, et al. Effect of surface polarity on the structure and dynamics of water in nanoscale confinement[J]. The Journal of Physical Chemistry B, 2009, 113(5): 1438-1446.

[37] Wei M J, Zhou J, Lu X H, et al. Diffusion of water molecules confined in slits of rutile TiO2(110) and graphite(0001)[J]. Fluid Phase Equilibria, 2011, 302(1-2): 316-320.

[38] Wei M J, Zhang L Z, Lu L H, et al. Molecular behavior of water in TiO2nano-slits with varying coverages of carbon: a molecular dynamics simulation study[J]. Physical Chemistry Chemical Physics, 2012, 14(48): 16536-16543.

[39] Wang K L, Zhao J G, Fu M L, et al. Experimental and molecular simulation study for the preparation of ordered mesoporous carbons[J].Aerosol and Air Quality Research, 2013, 13(3): 1034-1044.

[40] Whitesides G M, Grzybowski B. Self-assembly at all scales[J]. Science, 2002, 295(5564): 2418-2421.

[41] Cheng L S, Cao D P. Self-assembly of star-polymer-attached nanospheres for polymer nanocomposites[J]. The Journal of Physical Chemistry C, 2010, 114(13): 5732-5740.

[42] Baschnagel J, Binder K, Doruker P, et al. Advances in polymer science:viscoelasticity, atomistic models[J]. Statistical Chemistry, 2000,152: 41-156.

[43] Fukunaga H, Takimoto J, Doi M. A coarse-graining procedure for flexible polymer chains with bonded and nonbonded interactions[J].The Journal of Chemical Physics 2002, 116(18): 8183-8190.

[44] Muller-Plathe F. Coarse-graining in polymer simulation: from the atomistic to the mesoscopic scale and back[J]. Chem Phys Chem, 2002,3(9): 754-769.

[45] Mashayak S Y, Aluru N R. Thermodynamic state-dependent structure-based coarse-graining of confined water[J]. The Journal of Chemical Physics, 2012, 137(21): 214707.

[46] Maibaum L, Chandler D. A coarse-grained model of water confined in a hydrophobic tube[J]. The Journal of Physical Chemistry B, 2003,107(5): 1189-1193.

[47] Fukuda M, Zhang H D, Ishiguro T, et al. Adhesion properties of nanometer-thick perfluoropolyether films confined between solid surfaces: a coarse-grained molecular dynamics study[J]. Tribology Letters, 2013, 51(3): 479-487.

[48] Steele W A. The physical interaction of gases with crystalline solids I: gas-solid energies and properties of isolated adsorbed atoms[J].Surface Science, 1973, 36(1): 317-352.

[49] Steele W A. The interaction of rare gas atoms with graphitized carbon black[J]. The Journal of Physical Chemistry, 1978, 82(7): 817-821.

[50] Stecki J. Steele (10-4-3) potential due to a solid wall[J]. Langmuir, 1997, 13(3): 597-598.

[51] Zhao X C, Johnson J K. An effective potential for adsorption of polar molecules on graphite[J]. Molecular Simulation, 2005, 31(1): 1-10

[52] Takahiro O, Jin M, Katsumi K, et al. Adsorption properties of templated mesoporous carbon (CMK-1) for nitrogen and supercritical methane-experiment and GCMC simulation[J]. The Journal of Physical Chemistry B, 2002, 106(25): 6523-6528.

[53] Peng X, Cao D P, Wang W C. Heterogeneity characterization of ordered mesoporous carbon adsorbent CMK-1 for methane and hydrogen storage: GCMC simulation and comparison with experiment[J]. The Journal of Physical Chemistry C, 2008, 112(33): 13024-13036.

[54] Castro-Marcano F, van Duin A C T. Comparison of thermal and catalytic cracking of 1-heptene from ReaxFF reactive molecular dynamics simulations[J]. Combustion and Flame, 2013, 160(4): 766-775.

[55] Odegard G M, Jensen B D, Gowtham S, et al. Predicting mechanical response of crosslinked epoxy using ReaxFF[J]. Chemical Physics Letters, 2014, 591(20): 175-178.

[56] Kulkarni A D, Truhlar D G, Srinivasan S G, et al. Oxygen interactions with silica surfaces: coupled cluster and density functional investigation and the development of a new ReaxFF potential[J]. The Journal of Physical Chemistry C, 2013, 117(1): 258-269.

[57] Mayernick A D, Batzill M, van Duin A C T, et al. A reactive force-field (ReaxFF) monte carlo study of surface enrichment and step structure on yttria-stabilized zirconia[J]. Surface Science, 2010, 604(17): 1438-1444.

[58] Larsson H R, van Duin A C T, Hartke B. Global optimization of parameters in the reactive force field ReaxFF for SiOH[J]. Journal of Computational Chemistry, 2013, 34(25): 2178-2189.

[59] Jarvi T T, Kuronen A, Hakala M, et al. Development of a ReaxFF description for gold[J]. European Physical Journal B, 2008, 66(1):75-79.

[60] Raju M, Kim S Y, van Duin A C T, et al. ReaxFF reactive force field study of the dissociation of water on titania surfaces[J]. The Journal of Physical Chemistry C, 2013, 117(20): 10558-10572.

[61] Senftle T P, Meyer R J, Janik M J, et al. Development of a ReaxFF potential for Pd/O and application to palladium oxide formation[J].The Journal of Chemical Physics, 2013, 139(4): 044109.