苯并咪唑类缓蚀剂的定量构效关系研究

摘 要采用 2D-QSAR 和 3D-QSAR 相结合的方法研究苯并咪唑类缓蚀剂分子结构与缓蚀性能之间的关系。首先，采用量子化学密度泛函理论( DFT ) 对 15 种苯并咪唑类缓蚀剂的量子化学参数进行计算，并利用遗传函数近似算法( GA ) 建立量化参数与缓蚀效率之间的二维定量构效关系( 2D-QSAR ) 模型；其次，采用比较分子力场分析法( CoMFA ) 建立缓蚀剂分子的最优三维定量构效关系 ( 3D-QSAR ) 模型；最后基于 3D-QSAR 结果设计出了 14 种新型苯并咪唑类缓蚀剂分子，并利用 2D-QSAR 模型进行验证。结果表明：苯并咪唑类缓蚀剂的缓蚀性能与其分子总能量 ET、最高占据轨道能量EHOMO 以及自旋极化 a( Spin polarization ) 正相关；立体场和静电场对缓蚀剂活性的贡献分别为 52.1%和 47.9%，在 R1 位置引入体积大小适中且具有一定电负性的基团可提高缓蚀效率；设计缓蚀剂分子均具有较高的理论缓蚀效率，最优 2D-QSAR 模型对所设计分子的活性预测结果与 3D-QSAR 模型预测结果相吻合。 关键词：苯并咪唑；缓蚀剂；定量构效关系；分子设计ABSTRACTIn this thesis, studied the relationship between the benzimidazole inhibitor structure and corrosion inhibition performance of combined 2D-QSAR and 3D-QSAR method. Firstly, quantum chemical density functional theory ( DFT ) was used to calculate the quantum chemical parameters of 15 kinds of benzimidazole corrosion inhibitor. Than the genetic function approximation algorithm ( GA ) was used to study of two dimension quantitative structure-activity relationship(2D-QSAR) between quantum chemical parameters and their inhibition efficiency. Comparative molecular field analysis ( CoMFA ) method was used to establish the optimal 3D-QSAR model. Finally, based on the deeply reach of relationship between structure and inhibition efficiency of inhibitor molecules, design 14 kinds of new inhibitor molecules. The 2D-QSAR results show that the corrosion performance of benzimidazole inhibitors positive correlate with molecular total energy ET 、the highest occupied orbital energies ( EHOMO ) and the Spin polarization of atomic energy; negative correlate with the first carbon atom nucleophilic attack index. 3D-QSAR study results show that the proportion of steric field and electrostatic field contribution to the inhibitor activity were 52.1% and 47.9% , introduce a functional group with moderate volume size, and certain electronegativity in the R1 position will improve the inhibition efficiency of corrosion inhibitor. The 3D-QSAR model prediction of designed molecular show high theoretical inhibition efficiency. The inhibitor molecules active forecast results by optimal 2D-QSAR model for the design molecules consistent with the 3D-QSAR model predictionsKeywords: Benzimidazole; Corrosion Inhibitor; Quantitative Structure Activity Relationship; Molecular Design中国石油大学(华东)毕业设计(论文)目 录第 1 章 绪 论.11.1 引言11.2 苯并咪唑类缓蚀剂.21.2.1 苯并咪唑类缓蚀剂简介.21.2.2 苯并咪唑类缓蚀剂研究现状.21.3 二维定量构效关系方法.51.3.1 二维定量构效关系方法简介.51.3.2 二维定量构效关系方法研究现状.61.4 三维定量构效关系方法.61.4.1 三维定量构效关系方法简介.61.4.2 三维定量构效关系方法研究现状.71.5 本文研究内容.7第 2 章 理论计算方法和软件简介.92.1 软件简介.92.1.1 Materials Studio 软件简介 92.1.2 SYBYL 软件简介 .92.2 理论计算方法.102.2.1 量子化学计算方法.102.2.2 遗传算法.102.2.3 PLS 算法 11第 3 章 苯并咪唑类缓蚀剂的二维定量构效关系研究.133.1 引言.133.2 研究对象.143.3 计算方法.143.3.1 量子化学计算.143.3.2 2D-QSAR 模型的建立 153.4 结果与讨论.163.4.1 2D-QSAR 模型 163.4.2 2D-QSAR 模型验证 18中国石油大学(华东)毕业设计(论文)3.5 本章小结.22第 4 章 苯并咪唑类缓蚀剂的三维定量构效关系研究.234.1 引言.234.2 计算方法.234.2.1 分子活性构象的确定.234.2.2 分子叠合.244.2.3 比较分子场分析法(CoMFA)模型构建 254.3 结果与讨论.254.3.1 CoMFA 模型 .254.3.2 CoMFA 等势图 .264.3.3 CoMFA 模型验证 .27表 4-3 3D-QSAR 模型模型预测的缓蚀效率值及其标准偏差 274.5 本章小结28第 5 章 新型苯并咪唑类缓蚀剂的分子设计.305.1 引言.305.2 基于 CoMFA 结果的分子设计 .305.3 基于 2D-QSAR 结果的活性预测 325.4 本章小结.32结 论.34致 谢.35参考文献.36附 录.42中国石油大学(华东)毕业设计(论文)0第 1 章 绪 论1.1 引言酸性溶液广泛应用于油气井酸化及酸洗过程中，而酸液对油井金属设备和油气管线的腐蚀非常严重 1-3，因腐蚀导致的管壁爆裂、管线失效和原油大量流失等问题给石油化工行业带来了巨大的经济损失，并造成了严重的环境污染甚至生态灾难。近年来，向环境介质中添加合理的缓蚀剂来防止金属及其合金发生腐蚀成为各大油田比较推崇的一种经济有效的防腐蚀措施 4。目前针对特定环境的高效缓蚀剂还比较缺乏，而传统缓蚀剂的研制大多建立在猜测和大量探索性实验的基础之上，具有极大的随机性和盲目性，并缺乏理论指导。随着各种理论计算方法和分子图形技术的发展，计算机辅助缓蚀剂分子设计，尤其是定量构效关系( Quantitative structure-activity relationship，QSAR ) 方法 5的出现，为油气田缓蚀剂的分子设计提供了新思路。相比定量构效关系方法在药物分子设计、环境化学和毒理学等领域取得的丰硕成果，其在缓蚀剂方面的研究还处于起步阶段。缓蚀剂的二维定量构效关系( 2D-QSAR ) 研究，可以探讨各种量子化学参量与缓蚀效率之间的定量关系，建立简单有效的回归模型 6-10，但其与从分子层面实现缓蚀剂的合理设计还有一定的距离。三维定量构效关系( 3D-QSAR ) 11方法考虑到了配体与受体作用时的空间取向，在不了解受体结构的情况下，将分子周围的力场分布与活性定量地联系起来，明确给出对现有分子进行修饰的建议，弥补了 2D-QSAR 在分子设计中的不足。利用 3D-QSAR 方法进行缓蚀剂分子设计相当于在一种受体未知情况下的小分子药物设计 12-15，将该方法引入缓蚀剂领域，对实现缓蚀剂的定向化设计具有重要的理论指导意义。苯并咪唑类化合物对碳钢在盐酸环境中的腐蚀具有较好的抑制作用，是一类优良的酸性缓蚀剂 3, 16-19。本文以苯并咪唑类缓蚀剂为研究对象，采用 DFT 理论计算了研究对象的各种量子化学参数，并利用遗传函数近似算法建立了苯并咪唑类缓蚀剂的2D-QSAR 模型。然后采用比较分子场分析( Comparative molecular field analysis，CoMFA ) 方法，对其分子结构与缓蚀性能的三维定量关系进行研究，建立了稳定且预测能力较强的模型，并以此为基础设计出一系列新型高效缓蚀剂分子，再利用最优 3D-QSAR 模型预测所设计分子的活性。为了更加深入的研究新型缓蚀剂分子的缓蚀性能，计算所设计分子的量子化学参数，用最优 2D-QSAR 模型预测其活性，为新中国石油大学(华东)毕业设计(论文)1型缓蚀剂的合成提供可靠的理论指导。1.2 苯并咪唑类缓蚀剂1.2.1 苯并咪唑类缓蚀剂简介苯并咪唑衍生物是一类含有两个 N 原子的苯并杂环化合物，由于其 1、3 位上的两个 N 原子的高活性而极易在金属表面附着 20。作为一种新型的缓蚀剂，苯并咪唑衍生物在酸性溶液中对于碳钢、铝、铜等金属的具有优良的缓蚀性 21，且其具有热稳定性、无特殊刺激性气味、毒性低等特点，是一种绿色的环境友好型缓蚀剂。文中研究苯并咪唑类缓蚀剂主要由以下四部分组成：苯并咪唑母环、苯并咪唑上含有官能团的 R1、R 2、R 3 以及 R4 取代基，如图 1-1 所示。苯并咪唑类缓蚀剂由于取代基官能团的在金属表面不同的吸附性质，从而具有不同的缓蚀性能 22。常见的取代基包括烷基、芳香基及含 N、O 和 S 等原子的官能团。NN R1R2R3R4图 1-1 研究对象基本结构苯并咪唑类缓蚀剂通过在金属表面吸附成膜而起到缓蚀保护作用的。其苯并咪唑基是配位能力比较强的配体，取代基中 N、O 和 S 原子中未成键的孤对电子，容易和金属中的空 d 轨道形成配位键，从而能够在金属表面形成化学吸附和物理吸附。研究表明，苯并咪唑衍生物在金属表面吸附规律遵循 Langmuir 吸附等温式 3, 18, 23-27或Frumkin 吸附等温式 23, 28。电化学研究指出，苯并咪唑衍生物为混合型缓蚀剂，由于缓蚀剂吸附在金属表面导致电极表面的部分闭塞效应和能量效应，金属表面氢离子浓度降低，氢过电位提高，从而同时抑制钢和铜在酸性溶液中的阳极和阴极过程 16, 29。 1.2.2 苯并咪唑类缓蚀剂研究现状近年来，国内外学者针对苯并咪唑类缓蚀剂，主要依靠传统的重量法和电化学方中国石油大学(华东)毕业设计(论文)2法研究其缓蚀性能及缓蚀规律，也有部分学者从苯并咪唑类衍生物本身结构性质出发，利用量子化学方法试图从原子层次揭示其缓蚀机理。Popova. A16等人在 2003 年采用交流阻抗法、极化曲线法以及重量法研究了 5 种苯并咪唑林类缓蚀剂( BI、 2-NH2-BI、2-SH-BI、1-BBI、1,2-DBBI ) 对于碳钢的缓蚀性能。在 1 mol/L 的 HCl 溶液中，1-BBI 在金属表面产生特殊形式的吸附，活活能随温度的升高而降低，缓蚀效率随之增大。2-SH-BI 在金属表面发生化学吸附，其余缓蚀剂发生静电吸附。缓蚀性能是由于缓蚀剂在金属表面的空间位阻效应和能量效应导致表面氢离子浓度的减小和氢过电位的增大而产生。Popova. A23等人在 2004 年采用重量法和极化曲线法研究了 8 种苯并咪唑类缓蚀剂在 1 mol/L 的 HCl 中的缓蚀性能，结果表明，缓蚀剂缓蚀效率随其溶解度的增大而增大，5(6)-COOH-BI 吸附等温式符合 Langmuir 吸附等温式，其余均符合 Frumkin 吸附等温式，物理吸附为该类缓蚀剂具有缓蚀性能的根本原因。Popova. A30等人在 2006 年报道的文献中同样采用交流阻抗法、极化曲线法研究了 5 中杂环化合物在 1 mol/L 的 HCl 溶液中对于碳钢的缓蚀性能，结果表明苯并咪唑( BI )缓蚀性能低于一唑( I ) 、三唑( BTA ) 和噻唑 ( BNS ) ，高于苯并噻二唑( BTD ) ，苯并噻二唑会加重碳钢的腐蚀。Khaled K. F. 3采用相同的方法研究了 6 种苯并咪唑类缓蚀剂在 1 mol/L 的 HCl 溶液中对于纯铁的缓蚀效率，结果表明，该类缓蚀剂的在浓度为 0.001-0.010 mol/L 之间，缓蚀效率随浓度的增大而增大，其吸附服从 Langmuir 吸附等温式，在铁表面的吸附同时抑制了铁腐蚀的阴极过程和阳极过程，为混合型缓蚀剂。通过分析 EIS 等效电路，发现缓蚀效率与最高占据轨道能量( EHOMO ) 有关，缓蚀效率随分子供电子能力的提高而增大。Damin A. Lpez31等人利用交流阻抗谱和线性极化法研究了苯并咪唑类缓蚀剂在CO2 饱和的 NaCl 溶液中对于调质钢和退火钢的缓蚀性能，研究表明，苯并咪唑类缓蚀剂的加入可以抑制退火钢的腐蚀，加重调质钢的腐蚀，碳钢的微观机构会对该类缓蚀剂的缓蚀性能具有影响。苯并咪唑类缓蚀剂的缓蚀作用与其对阴极质子化的抑制作用有关，缓蚀剂分子在金属表面的吸附阻滞了 FeCO3 的沉淀。Zhang Da-Quan32等人采用类似手段研究了 2-巯基苯并咪唑在 0.5 mol/L 的 HCl 溶液中对于铜的缓蚀性能。结果表明该缓蚀剂为阴极型缓蚀剂，量子化学计算结果显示该缓蚀剂具有较高的最高未占据轨道能量( EHOMO ) 、最低空轨道能量( ELUMO ) 及 电中国石油大学(华东)毕业设计(论文)3子密度。巯基的引入改变了电子的扰动和能量配置，从而提高了缓蚀性能。Subramanian R.24等人采用循环伏安法和交流阻抗法研究了 4 种苯并咪唑类缓蚀剂和苯并三唑以及四唑的吸附性质对于碱性溶液中铜表面氧化膜形成的影响。结果表明，吸附服从 Langmuir 吸附等温式，有 95%的缓蚀剂膜在加入缓蚀剂 800 s 之内形成。当唑环与电极表面平行时可增大电极表面覆盖范围，环上 电子与铜表面 d 轨道的相互作用可促进吸附，屏蔽溶液中的电极，从而发挥缓蚀性能。Aljorani .J18等人采用电化学方法研究了 2-巯基苯并咪唑、2-甲基苯并咪唑和苯并咪唑在 1 mol/L 的 HCl 溶液中对于碳钢的缓蚀性能。提出该类缓蚀剂在碳钢表面的吸附符合 Langmuir 吸附等温式，其在碳钢表面的吸附为一个自发、放热的物理吸附过程。Morales-Gil P25等人利用电化学测试手段以及 SEM 分析研究了 4 种杂环化合物对于钢管试样在 1 mol/L 的 H2SO4 溶液中的缓蚀作用，同样提出该类缓蚀剂属于混合型缓蚀剂，其吸附服从 Langmuir 吸附等温式。2- 巯基苯并咪唑在铁表面的吸附可以抑制的硫化物的形成。Abboud Y26等人合成了一种新型苯并咪唑类缓蚀剂 2, 2-双苯并咪唑，并利用失重法和动电位极化曲线法研究了其在 1 mol/L 的 HCl 溶液中对于碳钢的缓蚀性能。结果表明，该缓蚀剂为一种混合型缓蚀剂，但其对于阴极的抑制作用强于阳极。其缓蚀效率随着浓度的增大而增大，在 0.0001 mol/L 时达到最高。缓蚀剂分子的缓蚀性来源于其在碳钢表面的吸附，其吸附符合 Langmuir 吸附等温式。Jacinto Morales Roque33等人采用密度泛函理论和交流阻抗法研究了 3( 苯并咪唑-2-亚甲基 ) 胺对于碳钢的缓蚀性能。结果表明， N=C=N 离域环通过其最高占据轨道( HOMO ) 为金属表面的最低未占据轨道( LUMO ) 提供电子，从而在金属表面吸附成膜，抑制金属腐蚀。EIS 研究结果表面缓蚀剂的加入使得电极和缓蚀剂溶液界面的电容减小，充放电效应增强，使得缓蚀剂易于在金属表面吸附成膜。Christov M28等人采用失重法测量了 4 种苯并咪唑类缓蚀剂的缓蚀效率，并利用测量结果确定了碳钢/酸性溶液/缓蚀剂体系的 Frumkin 吸附等温式参数。并计算了最大吸附量、相互作用参数以及平衡常数，理论计算结果与实际吸附行为相吻合。Obot I. O. 34等人采用 DFT 理论对苯并咪唑啉及其衍生物的缓蚀性能进行了理论研究。结果表明巯基在杂环化合物中的引入促进了电子和轨道能量配置的多样性，从而促进了缓蚀作用，这与 Da-Quan Zhang32等人的研究结果相类似。甲基作为供电子中国石油大学(华东)毕业设计(论文)4体的引入导致电子密度的增加，同样可增强缓蚀作用。Mieczyslaw Scendo35等人研究了苯并咪唑和 2-巯基苯并咪唑在 Cu2+还原为 Cu1+过程中的抑制作用。结果表明，两种缓蚀剂均可在 Pt 电极表面形成可以抑制 Cu2+还原为 Cu1+的缓蚀剂膜，该膜还起到抑制极限电流波对流扩散的作用。论文指出缓蚀剂的缓蚀性能依赖于电场作用下的物质传输，采用 EQCN 技术可以探测到 Cu2+进入缓蚀剂膜，在缓蚀剂膜中被还原为 Cu1+之后滞留其中，多以 CuCl 的形式存在，在缓蚀性强的膜中，CuCl 的量较小。Benabdellah M. 27等人采用失重法、交流阻抗谱法、动电位极化法研究了 2-巯基苯并咪唑在 1 mol/L 的 HCl 溶液中温度和缓蚀剂浓度对其缓蚀性能的影响。结果表明，缓蚀效率随缓蚀剂浓度的增大而增大，在 0.001 mol/L 时达到最高，为 98%。该缓蚀剂为一种混合型缓蚀剂，且随着缓蚀剂浓度的升高，传输阻力增大。在 308K-353K 温度范围内，缓蚀剂在金属表面的吸附符合 Langumuir 吸附等温式，为一个自发的放热过程。张军 36等人采用量子化学计算和分子动力学模拟相结合的方法，对 4 种苯并咪唑类缓蚀剂抗 HCl 腐蚀的性能进行了理论评价，并对其缓蚀机理进行分析。全局活性指数的计算表明，四种分子中，2-巯基苯并咪唑分子具有最强的反应活性；对于其他三种分子，Fukui 指数和全电子密度分布指出，2-氨基苯并咪唑具有两个亲电攻击中心，可在金属表面形成双中心吸附，其缓蚀性能应优于 2-甲基苯并咪唑和苯并咪唑；缓蚀剂分子与三层铁原子表面相互作用的分子动力学模拟进一步确认 2-甲基苯并咪唑比苯并咪唑在金属表面吸附更稳定。相比于传统的重量法、极化曲线法和交流阻抗谱法等实验方法，对于苯并咪唑衍生物类缓蚀剂进行计算机分子模拟和量子化学计算研究的报道并不太多，其定量构效关系的研究更是少之又少，苯并咪唑类缓蚀剂缓蚀机理也并没有统一说法，所以有待于在这些方面开展更加深入的研究。1.3 二维定量构效关系方法1.3.1 二维定量构效关系方法简介所谓定量构效关系，就是定量的研究和描述有机物结构与活性之间的相互关系 5。传统的二维定量构效关系方法主要由有 Hansch 法、基团贡献法和分子连接性指数法等37。其中应用最广泛的是 Hansch 法。该方法假设同系列化合物某些生物活性的变化与中国石油大学(华东)毕业设计(论文)5其某些可测量的物理化学性质( 疏水性、电性质和空间立体性质等 ) 的变化相联系，并假定这些因子是彼此孤立的，采用多重自由能相关法，借助多元线性回归等统计方法得到定量构效关系模型。基团贡献法是 Free-Wilson 基于有机物亚结构信息和生物毒性的相关研究的一种方法。该模式认为有机物与受体间的毒性效应是该有机物特定位置上不同取代基团毒性贡献的加和。Free-Wilson 法仅适用于具有相同母体结构的有机物，常被用来对有机物进行毒性初评。分子连接性指数法( Molecular connective index，MCI ) 是由 Kier 和 Hal 提出的。该方法根据分子中各个骨架原子排列或相连接的方式来描述其结构性质。MCI 是一种拓扑学参数，有零阶项、一阶项、二阶项、等等，可以根据分子的结构式和原子的点价计算得到，与有机物的毒性数据有较好的相关性。MCI 能较强地反映分子的立体结构，但反映分子电子结构的能力较弱，因此缺乏明确的物理意义，但由于其具有方便、简单且不依赖于实验等优点，近年来得到广泛应用和发展 38, 39。1.3.2 二维定量构效关系方法研究现状定量构效关系的研究起于 19 世纪中叶，最初考察的结构信息主要以化合物的化学性质为主，1863 年，Cros 首次观察到酒精浓度对哺乳动物的酒精毒性有着重要的影响。1868 年 Crum - Brown 和 Fraser 指出可用化学结构的函数来表示化合物的生物活性，由此发表了第一个 QSAR 方程 40。1990 年前后，Meyer 和 Overton 等提出了类脂学说，指出脂- 水分配系数 41对活性有重要影响，直到现在脂-水分配系数依然是描述生物活性的重要参量。尽管构效关系思想开展较早，但是直到 20 世纪 60 年代以物理化学参数为主的分子结参量才开始在 QSAR 研究中出现，例如分子的理化性质、量子化学参数和分子拓扑指数等，这类参量大多仅能反映分子的二维结构信息，因此所使用的方法被称为二维定量构效关系。传统的二维定量构效关系方法( 2D-QSAR ) 主要有 Hansch 法 42、Free-Wilson 法 43和分子链接性指数法( MCI，Molecular Connectivity Index ) 44等，其中以 Hansch 和 Fujita 等提出来的 Hansch 法应用最为广泛。由于 2D-QSAR 中使用的结构参量大多具有明确的物理意义，因此目前 2D-QSAR 方法依然在药物分子设计及其他定量构效关系研究中发挥着重要的作用。1.4 三维定量构效关系方法1.4.1 三维定量构效关系方法简介中国石油大学(华东)毕业设计(论文)6虽然 2D-QSAR 方法取得了一定的成功，但是传统 2D-QSAR 方法有许多局限性，其中最主要的就是在研究新化合物或新取代基的定量构效关系时，很难找到相对应的量子化学参数；以及在设计合成由 2D-QSAR 预测的高活性的化合物时，必须具有相当的物理有机化学、量子化学、计算化学等方面的知识；2D-QSAR 结果很难以图像的方式直观显示出来。与 2D-QSAR 相比，3D-QSAR 方法引入了与生物活性分子三维结构信息有关的因素作为变量，物理化学意义更加明确，能间接反应分子与受体之间的非键相互作用特征，更能深刻地揭示分子与受体的相互作用机理，并且分析结果以不同颜色势能面显示，便于直观。目前 3D-QSAR 研究主要是在药物化学、环境化学和生态保护等领域，并取得了巨大的成功。关于缓蚀剂及油田化学药剂方面的 3D-QSAR 研究尚未见报道，因此，借鉴 3D-QSAR 方法开展缓蚀剂分子设计及油田化学药剂研究方面的工作具有一定的创新性。1.4.2 三维定量构效关系方法研究现状自 20 世纪 80 年代诞生三维定量构效关系方法( 3D-QSAR ) 以来，其在药物分子设计、环境化学和毒理学等领域取得了丰硕的成果，已经取代 2D-QSAR 成为计算机辅助药物分子设计的主要方法，目前每年都有大量关于药物分子 3D-QSAR 研究的论文发表，其理论基础已经非常完善 45-47。3D-QSAR 认为分子周围的力场分布对活性有着重要的影响，因此以化合物的生物活性与分子周围不同位置力场分布情况之间的定量关系为基础，不仅可以分析分子与受体结构间的微观作用机理，而且还可以对新型同类分子的设计提供有力的建议。此外，3D-QSAR 结果以不同颜色势能面反映分子中各个基团对活性的贡献，可明确给出对现有分子进行修饰的建议。但是目前 3D-QSAR 方法在缓蚀剂领域还没有得到应用，将这种思想引入缓蚀剂分子设计能够显著提高该领域的理论深度。3D-QSAR 方法是基于配体的计算机辅助药物分子设计中应用最为广泛的一种，主要以受体结构未知的小分子为研究对象。而油气田常用的缓蚀剂分子多为有机物小分子，如苯胺类衍生物、咪唑啉及其衍生物、硫胺及其衍生物等。利用 3D-QSAR 方法进行缓蚀剂分子设计相当于受体未知情况下的小分子药物设计。如果将计算机分子模拟技术与 3D-QSAR 两种方法有机结合起来 48, 49，充分发挥二者的优势，扬长避短，一定能够更深入地理解分子三维结构对活性的影响机制，为新型缓蚀剂的研发提供有力的理论支撑。1.5 本文研究内容中国石油大学(华东)毕业设计(论文)7采用 Materials studio 软件中的 DMol3 模块对训练集缓蚀剂分子的量子化学参数进行计算，通过遗传函数近似算法( Genetic function approximation Algorithm ) 分析建立最优 2D-QSAR 模型，并利用测试集对模型预测能力进行验证。利用 SYBYL_X 1.2 软件 QSAR with CoMFA 模块建立训练集分子结构与缓蚀效率之间的最优 3D-QSAR 模型，分析分子周围立体场和静电场分布对缓蚀性能的影响，并在此基础上设计出一系列新型缓蚀剂分子结构。将 2D-QSAR 与 3D-QSAR 方法相结合，利用最优 2D-QSAR 模型对设计出的新型缓蚀剂分子活性进行预测，验证两种方法相结合的可能，为新型缓蚀剂的设计合成提供可靠的理论指导。中国石油大学(华东)毕业设计(论文)8第 2 章 理论计算方法和软件简介2.1 软件简介2.1.1 Materials Studio 软件简介Materials Studio 是 Accelrys 公司专门为材料科学领域研究者开发的一款可运行在PC 上的模拟软件，它可以解决当今化学、材料工业中的一系列重要问题。研究者可以通过 Materials Studio 操作平台方便地建立三维结构模型，并对各种晶体、无定型结构以及高分子材料的性质及相关过程进行深入的研究。模拟的内容包括了催化剂、聚合物、固体及表面、晶体与衍射、化学反应等材料和化学研究领域的主要课题 50。本论文主要采用的 MS 软件包中模块主要包括：( 1 ) Materials Visualizer：用于缓蚀剂分子结构搭建和编辑以及分子数据库的建立；( 2 ) DMol3：用于分子的优化和量子化学参数的计算；( 3 ) QSAR：用于苯并咪唑类缓蚀剂实验数据分析和 2D-QSAR 模型构建。2.1.2 SYBYL 软件简介SYBYL 软件是美国 Tripos 公司旗下的一款用于研究分子结构与性质、进行计算机辅助分子设计的一款专业化分子模拟软件，目前是国内外药物化学及环境化学领域应用最为广泛的软件，涵盖了分子结构搭建与优化、结构生物学、基于配体的分子设计、基于受体的分子设计、化学组合库、化学信息学、同源建模和虚拟筛选等内容。前期 SYBYL 软件只能用于 Linux 系统，其操作性和界面友好性都受到了一定的限制，近几年开发出的 SYBYL_X 软件可以分别适用于 Windows 和 Linux 系统，为使用者提供了友好的用户界面和更加完善的性能，目前该版本软件已经得到了广泛的应用。本论文选取 SYBYL_X 1.2 for Windows 版本，其中用到的模块或功能主要包括：( 1 ) SYBYL/Base：用于缓蚀剂分子结构的搭建和编辑以及分子数据库的建立；( 2 ) Dynamics：分子力学/分子动力学模块，用于缓蚀剂分子结构优化及最优构象的搜寻，软件中可以选择的力场有 Tripos、AMBER、MMFF94 等；( 3 ) Database Alignment：对数据库中的分子进行叠合，分析不同分子的结构差异；中国石油大学(华东)毕业设计(论文)9( 4 ) MSS( Molecular Spreadsheet ) ：用于训练集及测试集数据库文件的显示，分子周围各种分子场的计算和数据图谱分析等功能；( 5 ) QSAR with CoMFA：用于三维网格结构大小、扫描步长、探针原子的定义及缓蚀剂分子结构与缓蚀性能间 CoMFA 模型的构建；( 6 ) View QSAR：以不同颜色势能面的形式显示 CoMFA 模型的分析结果，为分子设计提供借鉴性指导。2.2 理论计算方法2.2.1 量子化学计算方法MS 中的量子化学计算方法主要包括从头算( ab-innitio ) 方法、半经验( Semi-empirical ) 方法以及密度泛函理论 ( DFT ) 方法。DFT 方法在 Hartree-Fock( HF ) 理论中，通过对电子动能和势能的平均化处理，借助变分法或数值方法，可以得到Schrdinger 方程的近似解。密度泛函理论中，体系状态的基本表征是电子密度而不是波函数。解密度泛函方程同样可以得到相当准确的体系能量以及其它性质。对大的体系，DFT 耗时比传统的超 HF 从头算方法要少可观的 1-2 个数量级。它也可以处理有机、无机、金属、非金属体系，几乎可以囊括周期表中的所有元素的化合物 51。近年来DFT 方法在缓蚀机理研究和缓蚀剂分子设计方面展现出较大的优势，并逐渐成为该领域强有力的研究工具。采用密度泛函理论可以处理较大分子，进行诸如电子转移与分布、轨道能量、电离能、偶极矩和 Fukui 函数等方面的分析和计算，这些研究提供了缓蚀过程中的基本信息，能够使人们更深刻地认识缓蚀剂作用机理。在本课题主要采用密度泛函理论计算苯并咪唑类缓蚀剂的量子化学参数。2.2.2 遗传算法遗传算法( GA ) 是一类借鉴生物界的进化规律( 适者生存，优胜劣汰遗传机制 ) 演化而来的随机化搜索方法 52。它是由美国的 J. Holland 教授 1975 年首先提出，其主要特点是直接对结构对象进行操作，不存在求导和函数连续性的限定；具有内在的隐并行性和更好的全局寻优能力；采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应地调整搜索方向，不需要确定的规则。遗传算法的这些性质，已被人们广泛地应用于组合优化、机器学习、信号处理、自适应控制和人工生命等领域。它中国石油大学(华东)毕业设计(论文)10是现代有关智能计算中的关键技术 53。 遗传算法的基本运算过程如下( 图 2-1 ) ： ( 1 ) 初始化：设置进化代数计数器 t=0，设置最大进化代数 T，随机生成 M 个个体作为初始群体 P( 0 )。 YNY执行选择算子随机产生初始种群执行变异算子执行交叉算子收敛准则满足？Random0,1PcRandom0,1Pm输出搜索结果YNN图 2-1 遗传函数算法流程图( 2 ) 个体评价：计算群体 P( t )中各个个体的适应度。 ( 3 ) 选择运算:将选择算子作用于群体。选择的目的是把优化的个体直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。选择操作是建立在群体中个体的适应度评估基础上的。 ( 4 ) 交叉运算；将交叉算子作用于群体。所谓交叉是指把两个父代个体的部分结构加以替换重组而生成新个体的操作。遗传算法中起核心作用的就是交叉算子。 ( 5 ) 变异运算：将变异算子作用于群体。即是对群体中的个体串的某些基因座上的基因值作变动。群体 P(t)经过选择、交叉、变异运算之后得到下一代群体 P(t 1)。 ( 6 ) 终止条件判断:若 t=T,则以进化过程中所得到的具有最大适应度个体作为最优中国石油大学(华东)毕业设计(论文)11解输出，终止计算。2.2.3 PLS 算法偏 最 小 二 乘 法 是 一 种 新 型 的 多 元 统 计 数 据 分 析 方 法 ， 它 于 1983 年 由 伍 德 ( S.Wold ) 和 阿 巴 诺 ( C.Albano ) 等 人 首 次 提 出 54。 近 几 十 年 来 ， 它 在 理 论 、 方 法和 应 用 方 面 都 得 到 了 迅 速 的 发 展 。 长 期 以 来 ， 模 型 式 的 方 法 和 认 识 性 的 方 法 之 间 的界 限 分 得 十 分 清 楚 。 而 偏 最 小 二 乘 法 则 把 它 们 有 机 的 结 合 起 来 了 ， 在 一 个 算 法 下 ，可 以 同 时 实 现 回 归 建 模 ( 多 元 线 性 回 归 )、 数 据 结 构 简 化 ( 主 成 分 分 析 ) 以 及 两组 变 量 之 间 的 相 关 性 分 析 ( 典 型 相 关 分 析 )。 这 是 多 元 统 计 数 据 分 析 中 的 一 个 飞 跃55。 偏 最 小 二 乘 回 归 可 以 解 决 主 成 分 分 析 法 漏 掉 相 关 性 较 小 变 量 的 问 题 。 它 采 用 对变 量 X 和 Y 都 进 行 分 解 的 方 法 ， 从 变 量 X 和 Y 中 同 时 提 取 成 分 ( 通 常 称 为 因 子 )，再 将 因 子 按 照 它 们 之 间 的 相 关 性 从 大 到 小 排 列 56。 现 在 ， 我 们 要 建 立 一 个 模 型 ，我 们 只 要 决 定 选 择 几 个 因 子 参 与 建 模 就 可 以 了 。作 为 一 个 多 元 线 性 回 归 方 法 ， 偏 最 小 二 乘 回 归 的 主 要 目 的 是 要 建 立 一 个 线 性 模型 ： Y=XB+E， 其 中 Y 是 具 有 m 个 变 量 、 n 个 样 本 点 的 响 应 矩 阵 ， X 是 具 有 p 个变 量 、 n 个 样 本 点 的 预 测 矩 阵 ， B 是 回 归 系 数 矩 阵 ， E 为 噪 音 校 正 模 型 ， 与 Y 具有 相 同 的 维 数 。 在 通 常 情 况 下 ， 变 量 X 和 Y 被 标 准 化 后 再 用 于 计 算 ， 即 减 去 它 们的 平 均 值 并 除 以 标 准 偏 差 。 偏 最 小 二 乘 回 归 和 主 成 分 回 归 一 样 ， 都 采 用 得 分 因 子 作 为 原 始 预 测 变 量 线 性 组合 的 依 据 ， 所 以 用 于 建 立 预 测 模 型 的 得 分 因 子 之 间 必 须 线 性 无 关 。 例 如 ： 假 如 我 们现 在 有 一 组 响 应 变 量 Y( 矩 阵 形 式 )和 大 量 的 预 测 变 量 X( 矩 阵 形 式 )， 其 中 有 些变 量 严 重 线 性 相 关 ， 我 们 使 用 提 取 因 子 的 方 法 从 这 组 数 据 中 提 取 因 子 ， 用 于 计 算 得分 因 子 矩 阵 ： T=XW， 最 后 再 求 出 合 适 的 权 重 矩 阵 W， 并 建 立 线 性 回 归 模 型 ：Y=TQ+E， 其 中 Q 是 矩 阵 T 的 回 归 系 数 矩 阵 ， E 为 误 差 矩 阵 。 一 旦 Q 计 算 出 来 后 ，前 面 的 方 程 就 等 价 于 Y=XB+E， 其 中 B=WQ， 它 可 直 接 作 为 预 测 回 归 模 型 。 在 建 模 当 中 ， 偏 最 小 二 乘 回 归 产 生 了 pxc 的 权 重 矩 阵 W， 矩 阵 W 的 列 向 量 用于 计 算 变 量 X 的 列 向 量 的 nxc 的 得 分 矩 阵 T。 不 断 的 计 算 这 些 权 重 使 得 响 应 与 其相 应 的 得 分 因 子 之 间 的 协 方 差 达 到 最 大 。 普 通 最 小 二 乘 回 归 在 计 算 Y 在 T 上 的 回归 时 产 生 矩 阵 Q， 即 矩 阵 Y 的 载 荷 因 子 ( 或 称 权 重 )， 用 于 建 立 回 归 方 程 ：中国石油大学(华东)毕业设计(论文)12Y=TQ+E。 一 旦 计 算 出 Q， 则 可 以 得 出 方 程 ： Y=XB+E， 其 中 B=WQ， 从 而 建 立 最终 模 型 。中国石油大学(华东)毕业设计(论文)13第 3 章 苯并咪唑类缓蚀剂的二维定量构效关系研究3.1 引言有机缓蚀剂的分子结构对其缓蚀性能有决定性的影响。近年来，随着计算机技术的迅猛发展，缓蚀剂的定量构效关系( QSAR ) 理论逐渐成为研究热点。定量构效关系研究的实质是从已知的样本集中提取有关化合物结构、活性的相关信息，建立某种数学模型，然后运用该模型去预测未知化合物的活性，从而为新分子的设计、评价提供理论依据。表 3-1 18 种苯并咪唑衍生物分子结构Compound Structure IE(%) Compound Structure IE(%)1NH52.760 10NH68.9502NHCl45.060 11N2CH378.0303NHO58.190 12NS3CH367.5454NH276.930 13 NH2O269.5805NHS58.125 14 NH2C3397.9856NHC2l52.800 15 NH2C397.7857NHC2O58.755 16 NHO250.7708NHC268.970 17NHO257.4909NHC2S53.985 18NHO2280.700中国石油大学(华东)毕业设计(论文)14因此，建立缓蚀剂 QSAR 模型( 或方程 ) 的首要条件是要保障有足够量的和可靠的样本集。样本集的数量与可靠性是构建 QSAR 模型的前提。通过量子化学计算，可以全面获得有关缓蚀剂分子的电子结构和几何结构的信息，如分子前线轨道能量、原子电荷、极化率、偶极矩，软硬度、Fukui 指数等。与传统的经验参数相比，量子化学参数对化合物结构的描述更加全面，理论意义更加明确，将量子化学参数应用于缓蚀剂的定量构效关系研究中已获得不少成功的模型 57。本章采用量子化学密度泛函理论( DFT ) 方法，计算苯并咪唑类缓蚀剂分子最稳定构型下的结构参数，建立所选 15 个缓蚀剂分子的各种量子化学参数与其抑制碳钢在 HCl 腐蚀环境中的缓蚀效率之间的定量构效关系模型，并对模型的科学性、稳定性及预测能力进行了评价，从量子化学角度揭示了该类缓蚀剂分子的活性差异。3.2 研究对象本文选取具有相同骨架结构的 18 种苯并咪唑缓蚀剂作为研究对象，分子结构如表3-1 所示，其缓蚀效率值取自实验文献 58，分别为在 1 mol/L 的 HCl 腐蚀环境中，添加50 mg/L 的不同类型缓蚀剂时对 Q235 钢的缓蚀效率，实验温度为 298K，腐蚀时间 10 h，通过失重法获得腐蚀速率，然后按照公式 3-1 计算得出缓蚀效率值：(3-1)00()/1%tIEW其中，W 0 和 Wt 分别代表加入缓蚀剂前后钢片的质量损失，单位为 g。3.3 计算方法3.3.1 量子化学计算本文所有分子的初始构型均采用 Accelrys 公司的 Materials Studio 软件包中的Visualizer 模块构建，并利用该软件包中的 DMol3 模块 59-61在不施加任何对称性限制的条件下计算所有研究对象，运用密度泛函理论广义梯度近似 GGA-PW9162交换关联函数，在双数值轨道基组加轨道极化函数( DNP ) 水平上进行全电子计算 62，热拖尾效应和实空间的截断半径分别采用 0.005 Ha 与 4.4 ，用以改进计算，所有计算均考虑自旋非限制性，并在该理论基础上对咪唑类衍生物中各分子的最高占据轨道能量( EHOMO ) 和最低未占轨道能量( ELUMO ) 、分子中原子的 Mulliken 电荷分布、各原子对前线轨道的贡献及 Fukui 指数等参量进行计算，并对其频率进行分析，确保计算后的分子结构均为势能面上的极小点。中国石油大学(华东)毕业设计(论文)153.3.2 2D-QSAR 模型的建立本论文选取 18 个苯并咪唑缓蚀剂分子中的前 15 个分子作为预测集分子构建QSAR 模型，剩余 3 个分子作为测试集。首先对前 15 个分子的实验测得的缓蚀效率进行单变量分析，经验证，缓蚀效率值呈现正态分布，保证了文中回归模型的合理性。前线轨道理论认为，发生化学反应时，电子的跃迁是由反应物分子前线轨道之间的相互作用引起的，前线轨道( 如最高占有轨道 HOMO、最低空轨道 LUMO ) 对化学反应起决定性作用。因此，论文选用了最高占有轨道能量 EHOMO、最低空轨道能量 ELUMO 及两者轨道能隙 ( 即分离能 ) E 作为表征分子间相互作用的参数。此外，根据传统化学理论，化学反应中存在着电子转移，而分子内的原子电荷是引起电子相互作用的主要原因，因此选取了分子中共有基团苯并咪唑环上非氢原子净电荷( qC( 1-6、8 ) 以及 qN( 7、9 )，苯并咪唑环原子编号如图 3-1 所示 )、咪唑环上非氢原子净电荷之和( qsum ) 作为表征咪唑啉类缓蚀剂分子电子特性的参数。Fukui 指数可以精确地定位缓蚀剂分子中每个原子的活性强弱及分子的亲核、亲电反应位，是确定分子反应活性中心及其强弱的有效方法。在外加势场 确定的条件rv下，电子密度 对电子数 N 的一阶导数即为 Fukui 函数 63：r rf(3-2)vrfrNFukui 函数 在有限差分近似下可表示为：rf(3-3)1iiifrq(3-4)iiiN式中， 表示分子为中性时原子 i 具有的电量、 qi( N+1 )和 qi( N-1 )则分别代表Nqi分子是阴离子和阳离子时，原子 i 所具有的电量； 和 则表示亲核 Fukui 指数rfifi与亲电 Fukui 指数，表征原子 i 得(失)电子能力的强弱，其值越大，得 (失)电子的能力越强 64。所以，论文还选取了苯并环上非氢原子的亲核( f+C( 1-6 )、f +C( 8 )、f +N( 7 )、f+N( 9 ) )亲电 ( f-C( 1-6 )、f -C( 8 )、f -N( 7 )、f -N( 9 ) )指数及其和作为苯并咪唑环在不同官能团取代结构的亲核亲电性用于构建 2D-QSAR 模型。另外，本论文还选取了不同方向的偶极矩 ( x )、( y )、( z )、总偶极矩 以及总能量 ET 和自旋极化率 a( Spin polarization ) 作为构建 2D-QSAR 模型的量子化学参数。利用 Accelrys 公司的 Materials Studio 软件包中的 QSAR 模块对建模用量子化学参中国石油大学(华东)毕业设计(论文)16数进行相关性分析，建立 2D-QSAR 模型。在模型建立过程中，论文所选取的遗传函数算法可以在数据回归过程中自动剔除掉相关性较大的参数，从而避免产生过拟合。本论文建模过程中将最高回归代数( Maximum generations ) 设为 1000，初始方程长度( Initial equation length，方程长度只方程中的自变量数 ) 设为 3，最大方程长度( Maximum equation length ) 设为 4，其余参数取默认值，然后通过遗传函数近似得到量子化学参量与缓蚀剂性能之间的关系式，即 2D-QSAR 模型。图 3-1 苯并咪唑环上非氢原子编号3.4 结果与讨论3.4.1 2D-QSAR 模型表 3-2 实验数据单变量分析统计参数Number of sample points Range Maximum Minimum Mean Median15 52.925 97.985 45.06 67.03 67.545Variance Standard deviationMean absolute deviationSkewness Kurtosis228.147 15.6347 11.938 0.706627 -0.54694对研究对象实测缓蚀效率单变量分析的结果如表 3-2，表中数据基本呈现正态分布，因而不再进行数据转换。正态分布最重要的参量是偏斜度( Skewness ) 和峰度( Kurtosis )。偏斜度是数据分布的三阶矩阵，可以表征数据分布的对称性。偏斜度若为正值，则数据分布向正方向偏移，对称分布的偏斜度为零。峰度是数据分布的四阶矩阵，表征与正态分布相关的峰态系数。单变量数值分析采用费舍尔峰度自动从以上定义中扣去3.0，正态分布峰度为 0.0；如果峰度是正值，则该组数据的分布要比正态分布陡峭，反之则较为平缓。中国石油大学(华东)毕业设计(论文)17利用相关系数矩阵对论文选取的量子化学参量进行相关性分析，结果如表 3-3( 见附录 )所示。表中显示了所有数据两两之间的相关系数，并将相关系数在 0.70000-0.90000 之间的系数用黄色标记，相关系数大于 0.90000 的用红色标记。分析表中数据可知，与缓蚀效率相关性较大的参数有 x 方向偶极矩 ( x )、最高占据轨道能量EHOMO、苯并环 7、9 号 N 原子所带电荷 qN( 7 )、qN( 9 ) 以及自旋极化 a。进一步分析这 5 个参量两两之间的相关系数可知，x 方向偶极矩与最高占据轨道能量 EHOMO 之间呈正相关，相关系数为 0.77908；与 7 号 N 原子所带电荷 qN( 7 ) 之间呈负相关，相关系数为-0.81957 ，与 9 号氮原子所带电荷 qN( 9 ) 及自旋极化 a 之间的相关系数较小，分别为-0.50615 、0.61858。最高占据轨道能量 EHOMO 与 7 号 N 原子所带电荷 qN( 7 ) 呈强负相关，相关系数为 -0.92056；与 9 号氮原子所带电荷 qN( 9 ) 之间呈较弱的负相关，相关系数为-0.58700，与自旋极化 a 之间呈较强正相关，相关系数为 0.70342；7 号 N原子所带电荷 qN( 7 ) 与 9 号氮原子所带电荷 qN( 9 ) 之间的相关性较弱，相关系数为0.42927；与自旋极化 a 之间呈较弱的负相关，相关系数为-0.48688。9 号氮原子所带电荷 qN( 9 ) 与自旋极化 a 之间呈较强负相关，相关系数为-0.76238。由于回归方程会自动避免相关系数高于 0.9000 的数据同时出现，以避免数据冗余。可以判定，在合理的 2D-QSAR 模型中，不会出现最高占据轨道能量 EHOMO 与 7 号 N原子所带电荷 qN( 7 ) 共存的现象。完成对相关量化参量相关性分析之后，采用遗传函数近似方法构建 10 个 2D-QSAR 模型。关系式分别如下： -HO+MTIE1=0.2386*E28.41739*48.256*fsum547fC() E6(3-5)TLUOI .90.51.907fC208263*f21743(3-6)HOMIE3 7.1950E5.9680*fC120.89*a (3-7)HOI4 6.238 * 47.25 f6 3.541a 15947(3-8)中国石油大学(华东)毕业设计(论文)18THO