杂原子对沥青质分子间相互作用影响的计算化学研究

叶正扬，龙 军，代振宇

(中国石化 石油化工科学研究院，北京 100083)

石油沥青质作为石油中的特征大分子表现出很强的聚集特性，以纳米聚集体(Nano-aggregation)或分形簇(Fractal cluster)的形式存在于原油、减压渣油及沥青中[1]。沥青质聚集会对石油采收、运输、加工各环节造成不利影响。在原油采收运输中，沥青质聚集会造成钻井口或管道堵塞，增加远程石油运输阻力；在石油炼制中，沥青质聚集将可能使分离或反应装置的管路堵塞，或造成催化剂中毒失效[2]。这些微观和介观非均相分布特征是由沥青质分子间相互作用所决定的。因此，对影响沥青质分子间相互作用定量研究的重要性日益凸显[3]。

沥青质中含有相当数量的杂原子，不同油源和分离方法[4]获得的沥青质杂原子含量和种类均有所不同，对沥青质的结构以及分子间的相互作用均有影响。定量地认识杂原子带来的影响，可以从元素组成的角度加深对沥青质分子间相互作用的基础认识[5]。

传统的实验手段很难实现对沥青质分子间相互作用的测定。富含沥青质的减压渣油或沥青体系所含的分子结构数以百万计，而每种分子的相对含量很低。因此，传统的分离手段只能利用分子的极性或溶解性差异分离出拟组分(Pseudo-component)[6]，例如四组分SARA(饱和分Saturates、芳香分Aromatics、胶质Resins、沥青质Asphaltenes)分离法，而很难分离出某种分子的纯净物进行相关的热力学量的测量。此外，对于沥青或减压渣油这样的复杂混合体系，某一种或几种分子无法决定整个体系的物理化学性质，因此，研究应当从寻找影响沥青质相互作用因素的规律入手[7]。许多研究表明，计算化学作为一种模拟实验的手段，对分子的结构、性质和相互作用能量的定量研究具有可靠性[8]，同时能够避免检测仪器和方法引入的系统误差。

1 研究方法与模型选择

1.1 沥青质分子结构特征

Mullins沥青质结构模型已被广为接受[9]。作为原油中相对分子质量最大，结构最复杂的分子种类，沥青质受到了大量的关注和研究。研究人员通过不同的实验手段对沥青质的相对分子质量、碳/氢比、不同杂化类型碳比、杂原子物种等有关沥青质的物理化学性质进行了深入的研究。根据荧光试验的结果，沥青质可抽象为椭球进行处理，长轴尺寸范围下限1.1～1.7 nm、上限2.5～2.8 nm[10-11]。Pomerantz等[12]通过两步激光解吸附/激光离子化质谱(Two-step laser desorption laser ionization mass spectrometry，L2MS)证明了岛式结构(又称手形结构)是沥青质的主要构型，平均相对分子质量在500～1500之间分布且峰值位于750[13-14]。

通过计算化学方法计算可能的沥青质结构的LUMO-HOMO能隙(这一数值与构成沥青质分子的稠环芳香核部分密切相关，受取代烷基侧链影响不大，故只考虑稠环芳香核的数量与排布)，对比沥青质样品的荧光发射光谱数据，Ruiz-Morales等[15-16]认为，绝大多数沥青质的稠环核由4～10个芳香环构成，其中最主要的是7环芳香烃结构。

1.2 沥青质模型分子的选择

为了研究杂原子对沥青质分子结构和分子间相互作用的影响，需要选择合适的模型分子作为研究的对象。沥青质分子建模的常用方法是按照沥青质样品的平均相对分子质量、平均碳/氢比、密度和平均杂原子含量等[7，17-18]平均结构参数建立一个平均分子模型。但是，这种方法并没有考虑沥青质分子的结构多样性，同时也很难找到一种符合沥青质样品整体平均性质的分子。由此，可以采用一种新的思路：针对沥青质中杂原子对沥青质结构及分子间相互作用的关系，选择结构相对简单、杂原子含量高的沥青质，尽量排除由相对分子质量、烷基侧链等其他因素带来的影响，通过这些结构合理且特点鲜明的模型分子，揭示杂原子对分子结构及其相互作用的影响。

Schuler等[19]采用原子力显微镜观察沥青质单分子结构的结果显示，原油沥青质与煤沥青质在结构上有一定差异，原油沥青质所含的饱和烷基侧链使得成像相对困难。研究人员发现了一种杂原子相对含量很高的沥青质结构(见图1)，其相对分子质量在500左右，在沥青质相对分子质量分布中相对较小，9个芳香环以渺位缩合的形式稠合，且其中3个环含有杂原子，分子不含长支链。为了研究杂原子的影响，模型分子中的五元环分别采用噻吩、环戊二烯和吡咯结构，构成含硫沥青质、烃类沥青质和含氮沥青质的3种模型分子(分别为AS(Sulphuric asphaltene)、AC(Hydrocarbon asphaltene)和AN(Nitrogen asphaltene))，具体结构见图1。精确量子化学计算中，类似长支链的柔性结构会造成结构优化在40 kJ/mol的能量范围内的振荡，致使结构优化无法达到收敛标准，这一误差可能掩蔽由杂原子引入的相互作用差别。因此，笔者借鉴Mousavi等[3]的工作，在与烷基侧链关系不大的研究对象中采用去长侧链的沥青质模型。当然，这样处理会降低相互作用的绝对值，同时也会降低模型分子的相对分子质量，但在不影响研究目标的前提下是可以接受的。

图1 沥青质模型分子Fig.1 Model molecules of asphaltene(a)Sulphuric asphaltene (AS)observed by atomic force microscope[19]；(b)Hydrocarbon asphaltene (AC)；(c)Nitrogen asphaltene (AN)

Schuler等[20]的工作无法给出所观测到分子的相对含量，因此并不能说明观察到的分子结构在沥青质中具有代表性。但依然可以从其他角度来说明笔者采用的模型分子的可靠性。根据Clar模型，笔者采用的模型分子符合Y规则(Y-rule)所定义的稳定结构沥青质，其孤立双键碳原子和六隅共振芳香碳原子数量比(C-ratio)为0.333，与沥青质的荧光实验结构吻合良好，是一种理论上很稳定的沥青质结构[21]。模型分子的结构同时也符合地球化学的实验结论：噻吩是沥青质中硫的主要存在形式[22]；吡咯结构来源于干酪根成岩或熟成过程中引入有机物骨架的铵盐[23]，并被原子力显微镜实验大量发现[20]；环戊二烯则来源于生物质五碳糖代谢产物[23]。

1.3 计算化学方法和参数设置

笔者采用的计算化学方法为密度泛函理论(Density functional theory,DFT)和分子动力学(Molecular dynamics,MD)，前者主要针对分子及其二聚体的结构优化和分子间相互作用计算，后者可以给出多分子系统的径向分布函数(Radial distribution function,RDF)和空间定向相关函数(Spatial orientation correlation function,SOCF)，由此对多分子在系统中的分布进行描述。

各分子的结构优化与分子间相互作用的DFT计算使用GAUSSIAN 16软件实现。泛函采用一种长程矫正杂化密度泛函ωB97X-D，该泛函自带阻尼原子-原子色散矫正(Damped atom-atom dispersion corrections)[24]，对解决非键相互作用有很好的准确度[25-26]。基组采用3-zeta的Karlsruhe系基组def2-TZVP，这一基组能较好地平衡分子间相互作用的计算精度和计算成本[27]。手动添加最小弥散函数到上述基函数，这一最小扩增泛函称为ma-def2-TZVP(或简称ma-TZVP)[28]。计算两分子相互作用时引入的基组重叠误差(Basis set superposition error,BSSE)通过Counterpoise方式解决[29]。MD计算软件采用Material Studio 2017R2，力场为COMPASSⅡ，该力场在COMPASS的基础上加入对杂环和离子液体的优化，更符合笔者研究的模型分子体系[30]。

1.4 计算方案

GAUSSIAN 16中所有计算采用DFT方法，泛函选用ωB97X-D，基组选择ma-TZVP，基组输入文件参考文献[28]。AS、AC和AN 3种模型分子结构优化完成后，利用Merz-Kollman布居分析实现对模型分子内部的静电势(Electrostatic potential，ESP)分布进行计算。

对3种模型分子的二聚体(AS-AS、AC-AC和AN-AN)进行结构优化后计算最低能量结构相互作用，利用Counterpoise方法对相互作用能量进行BSSE矫正。以分子间质心距离和分子取向夹角作为变量，在沥青质二聚体最低能量结构附近分别对2个变量进行手动调节。在研究二聚体质心距和分子间相互作用的关系时，保持分子取向夹角等于最低能量结构，手动调节质心距为最低能量结构的0.90、0.95、1.05、1.10、1.20和1.40倍。而在研究二聚体分子取向和分子间相互作用的关系时同理，保持二聚体质心间距等于最低能量结构，手动调节分子主轴夹角相较最低能量结构顺时针旋转30°、60°、90°、120°、150°和180°。调整后二聚体的相互作用能量计算方法同最低能量结构的分子间相互作用。

分子动力学计算在Material Studio 2017R中完成，计算力场均采用COMPASSⅡ，相互作用截断半径为1.55 nm。构建动力学计算盒子(Amorphous cell)并放入100种相同的模型分子，3种模型分子(AS、AC和AN)则对应3个动力学盒子，每个盒子均经过相同的优化及动力学模拟过程。分子动力学模拟在Forcite模块进行，经历4个阶段：(1)盒子结构优化和能量最小化；(2)退火，从200 K到500 K经历5次循环；(3)NPT系综，Nosé-Hoover恒温器(温度分别为250、293和353 K)和Andersen恒压器(压力为1 GPa)相应的条件下对盒子进行体积收缩，整个过程500 ps，步长0.5 fs；再对每个盒子分别在相同温度正则系综下进行动力学计算，这一过程500 ps，步长0.5 fs；(4)采用Nosé-Hoover恒温器和Andersen恒压器对每个盒子在阶段(3)中的温度和1×10-4GPa下，NPT系综进行运算，以模拟真实体系，总模拟时间大于1 ns，步长0.5 fs。

2 结果与讨论

2.1 模型分子能量最低结构及静电势分布

表1是3种模型分子AS、AC和AN的体积、偶极矩和四极矩的结果。由表1可以看出：AS的偶极矩略大于AN，而AN远大于AC。由于3种沥青质均为结构相对复杂的平面分子，偶极矩并不能很好地表征其内部的静电势分布。考察各分子的无迹四极矩，其绝对值|Θ|大小如表1所示。AC的无迹四极矩远小于另外2种模型分子；和偶极矩的结果不同，AN的|Θ|大于AS。3种沥青质模型分子均属于稠环芳香化合物，单从主对角元zz分量来比较，3种分子在垂直分子平面方向上的电子延展程度均远大于分子平面，但AC在该方向上的延展性要明显小于另外两者。

表1 沥青质模型分子的体积和静电势分布相关参数Table 1 Volumes and electrostatic potential distribution parameters of asphaltene model molecules

3种沥青质模型分子具体的电荷分布如图2所示。由图2可以看出：3种模型分子的电荷分布不均匀程度有所不同。从原子来看，AN分子的N原子静电势为-0.4～-0.55 e，而AS分子的S原子仅为-0.07～-2.7 e。但如果以芳香环作为分析单元则情况有所不同，AS分子中的并联噻吩环整体静电势为-0.547 e，远大于AN分子中并联吡咯环部分；后者只有-0.268 e，和AC分子同位置的并联异戊二烯烃环的-0.216 e接近。由此可见，不同杂原子的静电势大小和其所构成杂环整体的静电势大小变化趋势并不完全一致。

分子体积参考0.002 a.u.电子密度等值面来计算其凝聚态相的范德华体积。3种分子的体积由大到小的顺序依次为AS、AC、AN，这与各自所含杂原子的原子半径大小排列相似。

The numbers on the atoms represent the electrostatic potential.图2 模型分子静电势分布图Fig.2 Electrostatic potential distribution of model molecules(a)Sulphuric asphaltene (AS);(b)Hydrocarbon asphaltene (AC)；(c)Nitrogen asphaltene (AN)

2.2 模型分子二聚体结构优化及相互作用

对3种模型分子及其二聚体在最低能量下结构参数的计算结果如表2所示。所取的模型分子均为近似平面分子，分子内扭矩可以忽略，因此对两分子间最低能量的结构可以用分子质心距离来描述分子间距，用分子长轴夹角来描述分子取向[31-32]。为了可以和SOCF的结果更好地比较，分子选择长轴方向的无方向向量来定义分子取向。

由表2可以看出，分子间距按AS、AC、AN的顺序依次递减，AC与AN的分子取向相近，而AS与另外两者的分子取向相差略大。此外，分子的堆叠方式也有所不同：AS和AN两分子空间位置关系近似于沿垂直分子平面的法向量平移旋转，两分子主轴夹角分别为137°和153°；AC2个分子的堆叠位置关系则不同，其中一分子可以视作另一分子“翻面”后沿分子平面法向量进行平移旋转。

表2 模型分子二聚体最低能量下的结构参数Table 2 Structural parameters of model molecular dimers at the lowest energy

最低能量结构下，AS分子间的相互作用最弱，而AN之间最强，能量相差接近45 kJ/mol，杂原子对沥青质稠环芳核的分子间相互作用影响显著。

2.3 分子间相互作用与分子质心距离及分子取向的关系

在2.2节已经得到了最低能量构型下的分子间相互作用，在0 K真空环境下，所有的分子均会保持这一构型和相互作用状态。但在沥青的真实工作条件环境(250～350 K、1×10-4GPa)下，分子的振动和转动分子被解放，分子之间的相互作用和构型并不局限于最低能量态。对各模型分子偏离最低能量状态的相互作用和构型的关系进行研究，有助于理解杂原子对沥青质分子在真实环境下分子间相互作用的影响。

研究二聚体构型偏离最低能量结构的分子间相互作用可以参考单分子能量和构型关系的研究思路，采用势能面扫描的方法。在单分子构型与能量关系的研究中，以分子结构的单一结构参数作为变量(例如键长、键角或扭矩)，在最低能量结构附近的区间内依照一个极小步长进行扫描，由此获得该变量和分子能量之间的关系曲线。根据这条曲线可以估计分子获得能量后偏离最低能量结构的难易程度。对于分子间相互作用的计算而言，体系远大于单分子，计算量较大且需要的计算方法相对昂贵，因此，完全按势能面扫描方法精确扫描二聚体构型-相互作用能量曲线的耗时和计算资源的占用难以想象。笔者以各模型分子二聚体的分子质心间距和分子主轴夹角取向2个变量作为扫描结构参数，分别按照1.4节的方案对偏离最低能量构型的二聚体能量进行计算，从而绘制出结构参数-分子间相互作用的关系曲线，通过能量随结构参数变化的大小和对结构参数的敏感程度来判断不同模型分子二聚体偏离最低能量结构的难易程度，从而研究分子的有序性。

笔者定义相对质心距离为质心距离与最低能量状态下质心距离的比值，相对相互作用为分子间相互作用与最低能量状态下分子相互作用的比值。分子间相互作用及相对分子间相互作用与分子质心距离的关系见图3。不同沥青质分子间的相对质心距离与分子间相互作用能量呈现类似的关系，尽管相互作用能量的绝对值有所差异。这一结果在相对相互作用与相对质心距离的关系(见图3(b))中更加明显，在不同分子间这一关系差别不大。

AS—Sulphuric asphaltene；AC—Hydrocarbon asphaltene；AN—Nitrogen asphaltene图3 分子间相互作用与分子质心距离的关系Fig.3 Relationship between mass center distance and intermolecular interaction(a)Intermolecular interaction；(b)Relative intermolecular interaction

对比3种模型分子沥青质分子间相互作用和构型的关系曲线可以发现，最低能量态下的质心距离和分子的体积呈正相关；从3种模型分子相互作用-质心距离之间的关系来看，AN分子的相互作用对质心距离最不敏感；3种模型分子二聚体质心距偏离最低能量构型程度相同时，AN分子间的相互作用变化最小，而AC和AS尽管分子极性和体积差别较大，分子间相互作用和质心距的关系却相似。由此认为杂原子对质心距和分子间相互作用关系的影响仍有待进一步的研究。

分子取向与分子间相互作用的关系见图4。图4(a)是二聚体的分子取向夹角与分子间相互作用的关系，图4(b)是取向夹角与相对相互作用的关系。杂原子对分子间取向有显而易见的影响，但各模型分子取向与相互作用能量之间并没有明显的规律。AS分子间相互作用变化不大，最大能垒在150°处约为75 kJ/mol；AC分子间相互作用变化幅度略小于AS。最特殊的是AN分子，偏离最稳定结构后二聚体能量急剧上升，甚至相互作用在某些取向上为正值，这意味着分子间产生的巨大排斥作用已经超过了吸引作用，生成一个极为不利的构形；2个能垒分别位于30°和150°，最高相互作用高出最稳定构型316 kJ/mol，这种构形显然几乎不可能出现。

根据Hunter等[33]的研究，出现上述现象并不意外。对于取平行排列的高对称性分子对，π-π相互作用主要有以下原则：(1)分子错位排列的原因是静电相互作用规避电子之间的斥力，范德华相互作用在两分子全重叠排列下最强，但不足以补偿此构型下π电子排斥造成能量升高；(2)范德华相互作用主要与分子的大小有关，分子之间呈错位排列，吸引作用主要来自于分子骨架的σ和π电子云的相互吸引；(3)对于弱极性的π共轭体系，高电负性杂原子有利于增强吸引作用，并减小π-π排斥；(4)甲基有利于增强π相互作用，但不至于使共轭稠环芳核骨架结构发生形变；更长的饱和烷基侧链会导致π-π排斥增强，不利于错位平行堆叠的稳定。

笔者所取的模型分子属于含高电负性杂原子的积极取代弱极性稠环芳核，对称性低于Hunter采取的模型，但依然可以观察到2个分子之间错位排列的现象。由于3种模型分子除X位点原子外基本骨架完全一样，故在分子处于最低能量取向的构形时，仅改变分子质心距离，3个模型分子的能量随相对质心距离变化的规律基本一致的；而在取向发生变化时，π-π排斥急剧上升，从而出现较大能垒；对于分子间距较小的AN分子，取向改变所造成的影响最为明显，π-π排斥甚至完全抵消了范德华相互作用和σ-π吸引，导致二聚体能量高于2倍的单分子能量。

AS—Sulphuric asphaltene；AC—Hydrocarbon asphaltene；AN—Nitrogen asphaltene图4 分子取向与分子间相互作用的关系Fig.4 Relationship between molecular orientation angle and intermolecular interaction(a)Intermolecular interaction；(b)Relative intermolecular interaction

具体从结果来看，杂原子对分子取向的影响是由杂原子形成的负电中心造成的。氮原子作为AN分子的负电中心，在最稳定构型下交错排列以避免严重的π-π排斥；一旦分子取向发生变化，这些负电中心逐渐靠拢，产生巨大的排斥，这种排斥甚至在呈某些角度时超过了分子间相互吸引的作用。而对于AS分子，硫原子作为负电中心较氮原子而言负电势绝对值较小，同时受到硫原子体积影响分子间距较AC、AN2种分子明显更大，故而分子取向偏离引起相互作用绝对值减小的程度远不如AN。由此可知，决定分子取向和相互作用关系的关键因素应当是分子内原子静电势的分布。对比其他描述分子内部静电势分布的方案：偶极矩难以描述复杂分子的静电分布，AS有最大的偶极矩，但分子内部的杂原子中心负电势绝对值并不高。无迹四极矩显示AN的极性略大于AS，定性解释AN和AS分子取向-相互作用能量关系之间的巨大差异。如果以官能团(稠环)作为静电势分布的基本单位，AS中的并联噻吩环静电势最低，但这个结果无法解释AS分子取向相对自由的现象。可见，以稠合芳环这样的官能团作为基本单位分析分子内静电势分布依旧过于粗略，无法分析分子取向变化和分子间相互作用之间的关系。

2.4 多分子模型在常压和一定温度范围内的分布统计

在实际体系中，大量分子的分布和取向会偏离最稳定构型，具有一定的随机性，同时也受到温度和压强的影响。利用分子动力学方法，在压力为1×10-4GPa下，选择温度分别为250、293和353 K的环境模拟低温、室温和高温下的沥青质分子分布有关的RDF和SOCF。

图5是3种模型分子在不同温度下的径向分布函数图。从图5可以看出：3种模型分子均呈现短程有序，长程随机分布的状态；从峰形和分布来看，随着温度升高，3种模型分子的第一分布峰与参考原点之间的距离增加；温度对3种分子RDF第一峰位置的影响程度由大到小的顺序为AS、AC、AN，且AS的分布峰随着温度的上升明显变宽，可见硫原子的存在使分子更少地约束在最低能量结构的质心间距上，分子间距分布的随机性更强。由图3(a)可知，当不同模型分子的二聚体质心距偏离最低能量构型的程度相同时，AS分子能量上升最少，这一现象可以用来解释AS分子径向分布相对无序的现象。反观AN分子径向分布峰位置与峰高随温度变化不大，峰形相对另两种模型分子也更加尖锐，可见氮原子使得分子二聚体结构更加稳定，分布随机性更小。

SOCF可以观察到更明显的现象来支撑上述结论。图6是3种模型分子在不同温度下的空间定向相关函数图。随着温度的升高，分子间取向倾向于随机分布，低温(250 K)情况下3种模型分子均相对有序排列，而在较高温度下不同分子排列取向的差异逐步增大。从各模型分子二聚体分子间相互作用与分子取向的关系(见图4)可以看出，AS分子主轴夹角偏离最低能量构型时，其二聚体分子间相互作用的变化最不明显；由图6可知，在整个MD计算温度(250～353 K)范围内，AS分子的主轴取向分布相比其他2种分子更加随机，这和量子化学计算的结果相一致。具体来说，在温度为293 K和353 K 下，AS分子的空间定向相关函数波动无规律，且在某些距离上出现近乎垂直于参考向量的情况。而AN分子的空间定向相关函数相比之下则更有规律，即使在353 K的高温下取向也相对趋于一致，受温度的影响更小，这与2.3节展示的AN分子的极高转动能垒对应，显示了2种方法在定性上具有较好的对应性。

AS—Sulphuric asphaltene；AC—Hydrocarbon asphaltene；AN—Nitrogen asphaltene图5 各模型分子在不同温度下的径向分布函数Fig.5 Radial distribution functions of model molecules at different temperatures(a)250 K；(b)293 K；(c)353 K

AS—Sulphuric asphaltene；AC—Hydrocarbon asphaltene；AN—Nitrogen asphaltene图6 在不同温度下各模型分子的空间定向相关函数Fig.6 Spatial orientation correlation functions of model molecules at different temperatures(a)250 K；(b)293 K；(c)353 K

3 结 论

(1)杂原子的存在使分子内部的静电势分布区别明显，杂原子作为分子的负电中心存在。偶极矩对沥青质这样的复杂分子极性描述有局限性，不能很好地解释分子间相互作用和分子取向间的关系；而无迹四极矩能够定性描述这些复杂分子的极性以及所带来的相互作用特点变化。

(2)杂原子的存在通过影响π-π相互作用影响沥青质分子间相互作用的绝对值。在沥青质的稠环芳核中，对比不含杂原子的沥青质，噻吩结构的硫原子会减小分子间相互作用，而吡咯结构的氮原子则会增加分子间相互作用。

(3)最低能量构型下的分子间距和沥青质分子体积呈正相关，而结构近似的沥青质分子体积由其所含的杂原子决定。含硫沥青质、烃类沥青质和含氮沥青质3种模型分子的质心间距和相互作用大小的关系相近，其中的差别决定因素依旧需要进一步探索。

(4)杂原子能够影响沥青质分子二聚体的稳定性，主要通过影响分子间取向和分子间相互作用的关系来实现。当分子取向改变时，稳定构型的错位堆叠结构发生变化，π-π排斥不同程度增加；而含有不同杂原子的沥青质分子电荷分布不均匀，尤其是杂原子会成为整个分子的负电中心，故而氮原子的存在会大大增加取向改变时的能垒，使得分子的取向趋于固定。硫原子尽管也是一个较弱的负电中心，但由于其体积较大，使二聚体最低能量构型质心距较远，综合因素的结果反而会略微减小分子取向改变的能垒，使分子的取向更加自由。

(5)在不同温度下，含氮原子沥青质无论在分子分布和分子取向上保持着很好的稳定性，整体呈现较强的有序分布。而含硫沥青质则分布相对自由，尤其在高温状态下无序随机分布特点更加明显。