一种离子氢键交联超分子网络的分子模拟方法

技术领域

本发明属于高分子材料的计算模拟领域，尤其涉及一种离子氢键交联超分子网络的分子模拟方法。

背景技术

近十年来，基于非共价作用特别是氢键作用交联的聚合物材料已经成为人们研究的热点之一。当氢键作用强度足够强时，仅仅依靠它而构建的超分子聚合物网络便具有不错的力学性能；当它与其它非共价键协同作用时，它可以作为其中的牺牲键耗散能量来达到提高材料韧性，同时带来更好的自愈合能力。离子氢键正是其中一种特殊的氢键作用，它具有更强的相互作用，由此构成的超分子网络兼具较好地力学性能和自愈合能力，应用潜力非常好。然而，对于离子氢键作用在聚合物交联网络中的作用机制我们的认识还是有所不足的，大多数情况下研究人员都是通过间接的实验测试去解释材料的内在机理。分子动力学模拟经过多年的发展，现已成为研究物质微观结构的一种重要技术手段。在模拟中我们既可以轻松得到氢键的类型和数量信息，也可以观察它们的变化以及对材料内部结构的影响。由于一般的分子动力学模拟使用的都是不涉及化学反应的力场，于是在模拟含有离子氢键的体系时一定程度上忽视了基团之间的质子转移作用，造成结果的偏差。现今有关基于离子氢键交联的超分子材料模拟研究较少，如何以更为准确的方式在模拟中体现这种非共价作用对于这个领域的分子模拟研究具有重要意义。

发明内容

本发明旨在至少解决上述现有技术中存在的技术问题之一。为此，本发明提出一种离子氢键交联超分子网络的分子模拟方法。

在本发明的一些实施方式中，所述离子氢键交联超分子网络的分子模拟方法包括以下步骤：

S1：建立聚合物链原子结构模型，得到聚合物链的初始坐标文件；

S2：根据质子转移规律，修饰处理S1中所述聚合物链原子结构模型中可以形成离子氢键的基团，得到质子化的聚合物链结构模型；

S3：对S2所述质子化的聚合物链结构模型设置全原子力场，使其体系势能为键伸缩能、二面角扭曲能、范德华非键合势能以及库伦静电势能之和；

S4：选择原子电荷，对所述聚合物链进行片段分割，分别计算不同片段的原子电荷；

S5：确定分子链数量，建立模拟盒子，随机填充分子链，得到可用于模拟的低密度体系；

S6：对S5所述低密度体系进行平衡模拟流程，计算所述离子氢键基团之间的径向分布函数。

在本发明中，不使用计算量更大的反应力场的情况下，通过建立聚合物链中离子氢键基团的质子化模型并加以合理的原子电荷和模拟流程就可以取得不错的离子氢键交联超分子网络模拟效果。

在本发明的一些优选的实施方式中，使用Moltemplate程序构建聚合物链原子结构模型。

在本发明的一些优选的实施方式中，使用gromacs软件的pdb2gmx小程序处理聚合物链的初始结构文件得到了质子化的聚合物链结构模型和相应的拓扑关系文件。

在本发明的一些优选的实施方式中，上述S1中所述聚合物链的主链选自聚硅氧烷分子链、聚乙烯链、聚醚链中的任意一种。

在本发明的一些优选的实施方式中，S2中所述可以形成离子氢键的基团选自氨基以及不饱和羧基。

在本发明的一些优选的实施方式中，S2中所述修饰处理不包含离子氢键基团的重复单元结构保持不变，有机酸基团去掉氢原子而有机碱基团则得到一个氢原子的情况。

在本发明的一些优选的实施方式中，S3中所述全原子力场选自COMPASS力场、OPLS-AA力场、DREIDING力场、GAFF力场、AMBER力场；优选为OPLS-AA力场。

在本发明的一些优选的实施方式中，S4具体为：选择RESP原子电荷，所述片段用B97-3c方法在ORCA程序中进行几何结构优化，然后在B3LYP-D3/def2-TZVP级别的基组下计算优化后的不同片段的RESP电荷；在计算电荷时所述片段中间的重复结构单元两端用氢原子封闭，两个链端基则组合成一个完整分子；其中，规定不包含离子氢键基团的片段总电荷约束为0，存在一个有机酸基团的片段总电荷约束为-1，存在一个有机碱基团的片段总电荷约束为+1。

在本发明的一些优选的实施方式中，S5中所述分子链的数量范围为40～60。

在本发明的一些优选的实施方式中，S5中所述建立模拟盒子和所述随机填充分子链是通过Packmol程序完成的。

在本发明的一些优选的实施方式中，S5中所述模拟盒子为立方体，尺寸为10nm～15nm。

在本发明的一些优选的实施方式中，S5中不同的所述分子链的原子间距≥0.2nm。

在本发明的一些优选的实施方式中，S6所述平衡模拟流程依次为能量最小化、体系压缩、高温弛豫、降温、等温等压下的NPT平衡。

在本发明的一些优选的实施方式中，S6所述平衡模拟采用gromacs软件进行；其中所述能量最小化使用steep方法，它可以调整体系中不合理结构，能有效防止初次分子动力学模拟出现崩溃的情况。

在本发明的一些优选的实施方式中，S6中，在完成体系压缩后使用NVT系综，让模拟体系在600K下充分弛豫，模拟采用LINCS约束所有具有氢原子参与连接的化学键，步长为2fs，该过程的模拟总时间为200ns；弛豫过程结束后进行线性降温，温度从600K降到300K，同时放开键约束，步长改为1fs，时间为10ns；最后的平衡模拟采用NPT系综，在300K和1.0bar(接近一个大气压)下进行100ns的模拟；所有的平衡模拟的温度控制采用V-rescale方法，耦合时间常数为0.2ps，压力控制采用Berendsen方法，耦合时间常数为0.5ps。

在本发明的一些优选的实施方式中，采用密度和总势能作为体系平衡态判定的标准。

在本发明的一些优选的实施方式中，所述平衡模拟中平衡状态的判定标准为体系的密度和势能变化＜5％。

本发明的有益效果为：

(1)本发明有助于解决涉及离子氢键作用的分子动力学模拟问题，合理运用质子化的初始结构模型可以在模拟中更好地体现出离子氢键的强相互作用。

(2)本发明采用片段分割法计算大分子链的原子电荷，保证准确性的同时大大减少了计算资源，为大分子量的聚合物单链原子电荷计算提供了思路。

附图说明

下面结合附图和实施例对本发明做进一步的说明，其中：

图1为本发明实施例1中氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷的化学结构式。

图2为本发明质子化的模型和现有技术非质子化模型中氨基和不饱和羧基质心之间的径向分布函数图。

具体实施方式

以下将结合实施例对本发明的构思及产生的技术效果进行清楚、完整地描述，以充分地理解本发明的目的、特征和效果。显然，所描述的实施例只是本发明的一部分实施例，而不是全部实施例，基于本发明的实施例，本领域的技术人员在不付出创造性劳动的前提下所获得的其他实施例，均属于本发明保护的范围。

实施例1

本实施例以氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷为例，实现一种离子氢键交联超分子网络的分子模拟方法，具体过程为：

(1)建立氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷两种聚合物链的全原子结构模型，其化学结构式如图1所示，两者的聚合度都为160，接枝率都为5％，其中的氨基和不饱和羧基都是随机接枝在聚硅氧烷主链上的；

(2)根据质子转移规律，对其中可以形成离子氢键的基团进行修饰处理，在这个模拟体系中，氨基和不饱和羧基在形成离子氢键的过程中存在质子转移，不饱和羧基脱去一个氢原子，转变为去质子化状态；而氨基则得到一个氢原子，变为质子化状态，依此对(1)中得到的全原子结构模型进行修改。

(3)对氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷非共价交联体系选择OPLS-AA力场，这种力场应用范围广，原子类型众多，对一些复杂的有机结构也有相应的力场参数描述。该力场的总势能为键伸缩能、键角弯曲能、二面角扭曲能、范德华非键合势能以及库伦静电势能之和。依据文献，Si的非键合范德华作用参数σ和ε分别定为0.4nm和0.4184kJ/mol，以Si为中心的二面角扭曲项参数设置为0。

(4)选择RESP原子电荷用于该体系的模拟，为了减少计算量，将氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷分割成如下片段：二甲基硅氧烷单元、氨基接枝单元、不饱和羧基接枝单元以及首尾端基，处于中间的重复单元用氢原子封闭，首尾端基则是组合在一起，然后分别进行原子电荷计算。上述片段都是用B97-3c方法在ORCA程序中进行几何结构优化，然后在B3LYP-D3/def2-TZVP级别的基组下计算优化后的片段的RESP电荷；电荷计算中设置必要的约束，其中二甲基硅氧烷片段总电荷约束为0，而带有一个质子化氨基的片段总电荷约束为+1，去质子化的羧基部分则约束为-1。

(5)模拟体系中氨基接枝聚硅氧烷链和不饱和羧基接枝聚硅氧烷链的数量都设置为20，这样可以在保证模拟结果的有效性下不必耗费过多的计算资源；使用Packmol程序完成建立模拟盒子和填充分子链的任务，这40条聚合物链随机堆积在一个大模拟盒子中，构成了可用于模拟的低密度体系，即建立了质子化初始结构模型。

(6)在建立好质子化初始结构模型的基础上，加上合适的模拟流程才能得到更好的模拟效果。模拟采用gromacs软件完成，其中平衡模拟流程主要由体系能量最小化、体系压缩、高温弛豫、降温、常温常压下的NPT平衡组成；能量最小化使用的是steep方法，它可以调整体系中不合理结构，能有效防止初次分子动力学模拟出现崩溃的情况。在完成体系压缩后使用NVT系综，让模拟体系在600K下充分弛豫，模拟采用LINCS约束全部有氢原子参与连接的化学键，步长为2fs，该过程的模拟总时间为200ns；弛豫过程结束后体系会经历一个线性降温的阶段，温度从600K降到300K，同时放开键约束，步长改为1fs，时间为10ns。最后的平衡模拟采用NPT系综，在300K和1.0bar(接近一个大气压)下进行100ns的模拟；所有的平衡模拟的温度控制和压力控制分别采用V-rescale方法和Berendsen方法，耦合时间常数分别为0.2ps和0.5ps。采用密度和总势能作为体系平衡态判定的标准，当其波动值处于5％的范围内时则认为体系基本处于平衡态。

对比例1

本对比例与实施例1的区别在于本对比例建立的模型为非质子化模型，其余步骤与实施例1相同。具体过程为：建立了聚合度为160，接枝率为5％的氨基接枝聚硅氧烷和不饱和羧基接枝聚硅氧烷两种聚合物链结构模型；采用片段分割法计算RESP原子电荷，每个片段的总电荷都约束0，计算使用的基组和流程都与实施例1保持一致；构建了20条氨基接枝聚硅氧烷和20条不饱和羧基接枝聚硅氧烷组成的模拟体系，平衡模拟使用了与实施例1相同的OPLS-AA力场参数、模拟流程和参数。

试验例

通过分析平衡模拟中氨基和不饱和羧基之间的径向分布函数以及氢键数量来表征实施例1中质子化结构模型和对比例1中非质子化结构模型在模拟中对离子氢键作用的描述效果。如图2所示，两种结构模型中的氨基和不饱和羧基质心之间的径向分布函数表现出明显的差异，质子化模型的峰值更高，说明了其中的离子氢键基团彼此之间的聚集程度更高。表1显示的是两种结构模型在平衡态中主要基团之间形成的氢键数量平均值。

表1两种结构模型在平衡态中主要基团之间形成的氢键数量平均值

由此可以知道在质子化结构模型中氢键数量更多且氨基和不饱和羧基之间形成的氢键占据主导地位。从对比的结果来看，相较于一般的非质子化模型，本申请经过修改处理后的质子化模型在模拟中能更好地体现离子氢键这种强相互作用，可以为这一类体系的模拟提供指引思路。

上面对本发明实施例作了详细说明，但是本发明不限于上述实施例，在所属技术领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化，皆属于本发明的范围之内。