一种基于分子模拟的油析热力学相图的预测方法

1.本发明属于化学工程工业结晶技术领域，具体涉及一种基于分子模拟的油析热力学相图的预测方法。


背景技术：

2.在溶析或者冷却结晶过程中，溶液达到过饱和状态后，在出晶前会出现新的液相，这一现象通常被称作油析现象，有时也被称作液液分层或者叫作液液相分离。油析现象发生后一般会对产品质量造成影响，大部分的研究都表明油析会使晶体包藏杂质，降低最终产品的纯度，还会造成产品聚结。此外，在大规模工业生产中，因油析增加生产成本和安全隐患也同样是不可忽视的。黏附在搅拌装置上的油相会增加搅拌阻力，进而增加设备能耗。同时，黏附在结晶器器壁的油相增加了污垢热阻，扰乱了正常的传热过程。
3.对于有机小分子物质油析现象的报道目前较多着重于油析结晶过程的热力学性质研究，很多研究人员主要是通过测定油析结晶过程中的相图来研究结晶过程的油析现象及油析现象的控制。通过热力学相图研究油析区域形成的条件及影响因素，可为实际工艺设计提供直接的优化方案，对改善操作条件、优化结晶工艺、提高产品质量具有重要意义。然而，由于缺乏高效便捷的测量方法，油析体系热力学相图的测定需要耗费大量的时间和人力成本。更严重的是，测定过程通常使用大量有害生物和昂贵的有机溶剂，对环境极易造成污染。
4.目前油析相图的预测方法多数集中于热力学模型拟合方法，例如nrtl模型、pc-saft模型等。其中有研究者根据二羟基二苯砜的溶解度数据利用pc-saft模型成功预测了在二羟基二苯砜-水-丙酮和二羟基二苯砜-水-异丙醇三元体系中的油析现象。然而，以上方法均属于宏观热力学研究的范畴，精度和适用范围有限，从分子水平上准确高效地预测油析热力学相图仍然是一个难题。
5.因此，探索更廉价、更有效、更快速且普适性好的油析热力学相图预测方法仍是一项重大而紧迫的挑战。


技术实现要素：

6.为了克服现有基于热力学模型的油析热力学相图预测方法精度低、适用范围窄的缺陷，
7.本发明提供了一种基于分子模拟的油析热力学相图的预测方法，通过基于密度泛函理论的量子化学计算、波函数分析以及分子动力学模拟的方法，从分子间相互作用角度出发，对其油析现象发生的溶剂体系和油析曲线进行了预测，可以解决上述背景技术中提到的问题。
8.本发明的技术方案如下：
9.(1)确定待预测物质的化学结构。
10.(2)对该物质进行气相下的构型搜索，进而得到气相下稳定的二聚体构型。
11.(3)对气相下稳定的二聚体构型进行独立梯度模型分析，判断在二聚体构型中溶质分子间的相互作用类型。同时对气相下稳定的二聚体构型进行分子静电势表面分析，判断二聚体构型此时能够提供氢键供体和受体的能力。
12.(4)以分子静电势表面分析的结果为基础计算二聚体构型的分子极性指数值，结合独立梯度模型和分子静电势表面分析的结果对溶剂范围进行划定。
13.(5)选择在筛选范围内的溶剂，对该物质进行溶剂相下的构型搜索。如果依旧能得到稳定的二聚体构型，则把该溶剂视为可以发生油析现象的溶剂体系；如果二聚体构型被破坏，则需要更换溶剂后重新进行构型搜索。
14.(6)在获得溶剂体系后，构建特定溶质浓度的无定形单元，通过分子动力学模拟来计算特定温度下油析体系的总能量，利用所有温度的模拟结果即可绘制总能量-温度曲线。
15.(7)以特定体系下的总能量-温度曲线为基础，结合剪应力自相关函数分析、径向分布函数分析和表面能计算的结果，可推测出油析发生的温度。
16.(8)将不同浓度下的油析温度点相连接即可得到该体系的油析曲线。
17.(9)发生油析的溶剂体系和油析曲线即可组成油析热力学相图。
18.步骤(1)中待测物质分别为吡唑醚菌酯、香兰素、乙基香兰素、柠檬酸和苯甲酸。
19.步骤(2)、(3)和(5)中二聚体构型搜索、几何优化以及能量计算均是通过orca软件先在b97-3c的方法下进行预优化，再在不低于b3lyp-d3(bj)/ma-def2-tzvp的精度下完成的。
20.步骤(3)和(4)中独立梯度模型和分子静电势表面的分析是利用multiwfn软件，通过读取由orca软件计算产生的molden文件中的能量信息和波函数信息计算得来的。
21.步骤(4)中溶剂体系分别为水、甲醇、乙醇、1-丙醇、2-丙醇、1-丁醇、2-丁醇、异丁醇、1-戊醇、2-戊醇、异戊醇、1,2-丙二醇、丙三醇、异丙醚、n-己烷、环己烷、正辛烷、异辛烷、二氯甲烷、乙酸甲酯、乙酸乙酯、乙酸正丙酯、乙酸异丙酯、乙酸正丁酯、丙酮、乙腈、n-甲基吡咯烷酮、n,n-二甲基甲酰胺、二甲基亚砜、1,4-二氧六环。
22.步骤(6)和(7)中剪应力自相关函数分析、径向分布函数分析和表面能计算是利用materials studio软件的forcite模块在compass力场下实现的，在分子动力学和分子力学模拟中，andersen恒温器被用来控制温度，berendsen恒压器被用来控制压力。
23.上述方法中具有以下有益效果：
24.a)本方法实现了对吡唑醚菌酯等五种物质油析相图的预测，亦可适用于其他有机化合物，更利于推广应用。
25.b)本方法应用高精度的分子模拟手段和波函数分析，通过多角度考察分子团簇之间及与溶剂分子之间的相互作用强度和类型，确立最佳的预测指标，增加预测的准确性，提升预测能力。
26.c)本方法输入文件仅需要化合物与溶剂分子的分子结构，不需要大量的热力学数据，避免了测定数据耗费的大量时间、人力成本以及有害生物和昂贵的有机溶剂使用。
附图说明
27.图1：本发明中油析热力学相图预测方法的流程图；
28.图2：本发明中吡唑醚菌酯二聚体构型独立梯度模型分析结果；
29.图3：本发明中香兰素、乙基香兰素、柠檬酸和苯甲酸在选定溶剂体系下油析现象的辅助实验照片；其中，其中：(a)为香兰素在水中的油析现象的辅助实验照片，(b)为乙基香兰素在水中的油析现象的辅助实验照片，(c)为柠檬酸在水中的油析现象的辅助实验照片，(d)为苯甲酸在水 异丙醇中的油析现象的辅助实验照片。
30.图4：本发明中吡唑醚菌酯特定浓度下的总能量—温度曲线；
31.图5：本发明中吡唑醚菌酯在环己烷和异丙醇混合溶剂下油析相图的预测和实验结果对比图；
32.图6：本发明中香兰素、乙基香兰素、柠檬酸和苯甲酸在各自能够发生油析现象的溶剂体系下的油析相图预测和实验结果对比图；其中：(a)为香兰素在水中的油析相图预测和实验结果对比图，(b)为乙基香兰素在水中的油析相图预测和实验结果对比图，(c)为柠檬酸在水中的油析相图预测和实验结果对比图，(d)为苯甲酸在水 异丙醇中的油析相图预测和实验结果对比图。
具体实施方式
33.下面结合附图，对本发明的技术方案进行具体说明。
34.应该指出，以下详细说明都是示例性的，旨在对本技术提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本技术所属技术领域的普通技术人员通常理解的相同含义。
35.一种基于分子模拟的油析热力学相图的预测方法，包括以下步骤：
36.(1)按照图1所示的流程，首先进行物质分子的建模。吡唑醚菌酯、香兰素、乙基香兰素、柠檬酸和苯甲酸作为待预测物质，同时选择水、甲醇、乙醇、1-丙醇、2-丙醇、1-丁醇、2-丁醇、异丁醇、1-戊醇、2-戊醇、异戊醇、1,2-丙二醇、丙三醇、异丙醚、n-己烷、环己烷、正辛烷、异辛烷、二氯甲烷、乙酸甲酯、乙酸乙酯、乙酸正丙酯、乙酸异丙酯、乙酸正丁酯、丙酮、乙腈、n-甲基吡咯烷酮、n,n-二甲基甲酰胺、二甲基亚砜、1,4-二氧六环作为候选溶剂体系，并对分子进行建模。
37.(2)按照图1所示的流程，接着对该物质进行气相下的构型搜索。通过orca软件将这一批不同的二聚体构型先在b97-3c的方法下进行预优化，再在不低于b3lyp-d3(bj)/ma-def2-tzvp的精度下计算每个二聚体构型的单点能，并进行排序，最终得到气相下稳定的一批二聚体构型。
38.(3)按照图1所示的流程，对气相下稳定的二聚体构型进行表面静电势分析和独立梯度模型分析。基于multiwfn软件，通过读取由orca软件计算产生的molden文件中的能量信息和波函数信息，对气相下稳定的二聚体构型进行独立梯度模型分析，判断在二聚体构型中溶质分子间的相互作用类型。非共价相互作用的梯度等值面所包围的区域表示存在弱相互作用的区域，对应散点图所映射的颜色表示交互的类型(蓝色代表氢键，绿色代表色散，红色代表位阻)。同时对气相下稳定的二聚体构型进行分子静电势表面分析，判断二聚体构型相对原本单独的分子其能够提供氢键供体和受体的能力是否发生改变，供受体能力的变化会导致与特定溶剂结合的能力发生改变。
39.(4)按照图1所示的流程，以分子静电势表面分析的结果为基础计算二聚体构型的分子极性指数值。分子极性指数的大小反映了二聚体构型极性的大小，结合独立梯度模型
和分子静电势表面分析的结果对溶剂范围进行划定，例如，吡唑醚菌酯二聚体构型中溶质分子之间均存在大面积绿色区域，体现相应弱相互作用区域电子密度很小，因此为色散作用主导的π
…
π堆积相互作用(如图2)，同时分子静电势表面分析结果表明，其正极值位于甲氧基-och3位置，负极值位于酯羰基c＝o位置，且相比单独吡唑醚菌酯分子，二聚体构型所具有的氢键供体能力大幅降低，但其仍然保持较高水平的氢键受体能力。因此，所选择的溶剂体系只有同时满足具备较弱的氢键供体能力和与二聚体构型的极性存在较大差异两个条件，才能维持原有的二聚体构型。因此，可以判断出分子极性指数较小且氢键供体能力较弱的溶剂(正己烷、异辛烷、环己烷、正辛烷、异丙醚、2-戊醇、异戊醇和异丁醇)不会对溶质的二聚体构型产生影响，进而可能发生油析现象。
40.(5)按照图1所示的流程，选择在筛选范围内的溶剂，对该物质进行溶剂相下的构型搜索。通过orca软件将这一批不同的二聚体构型先在b97-3c的方法以及溶剂化模型下进行预优化，再在不低于b3lyp-d3(bj)/ma-def2-tzvp的精度同时在溶剂化模型下计算每个二聚体构型的单点能，并进行排序，最终得到特定溶剂环境下稳定的一批二聚体构型。如果依旧能得到稳定的二聚体构型，则把该溶剂视为可以发生油析现象的溶剂体系；如果二聚体构型被破坏，则需要更换溶剂后重新进行构型搜索。例如，分别在正己烷、异辛烷、环己烷、正辛烷、异丙醚、2-戊醇、异戊醇和异丁醇体系下进行构型搜索，结果表明吡唑醚菌酯依旧以二聚体构型存在，因此，上述溶剂体系被认为是可以发生油析现象的溶剂体系。而辅助实验的结果表明，以上溶剂体系除了正辛烷以外均发生油析现象，证明了该方法的准确性。同理，在香兰素、乙基香兰素、苯甲酸和柠檬酸的例子中，按照上述方法，成功预测出其对应能够发生油析现象的溶剂体系见图3((a)香兰素(水)、(b)乙基香兰素(水)、(c)柠檬酸(水)和(d)苯甲酸(水 异丙醇))。
41.(6)按照图1所示的流程，接下来进行分子动力学的模拟。在获得溶剂体系后，选择其中一组溶剂体系，构建特定溶质浓度的无定形单元，每个无定形单元包含500个分子，无定形单元的三个边长均大于以消除周期性键链的影响。基于materials studio软件的forcite模块，采用compass力场，同时使用andersen恒温器被用来控制温度，berendsen恒压器被用来控制压力。首先通过对周期性单元进行100000步基于分子力学模拟的几何优化，消除不相关的相互作用。接着在不同温度下依次进行了nvt(一种由固定粒子数量(n)、固定体积(v)和固定温度(t)组成的系统)-npt(一种由固定粒子数量(n)、固定压力(p)和固定温度(t)的组成的系统)-nvt系综分子动力学模拟，以保证体系的良好松弛并达到平衡。每个分子动力学过程的模拟时间设置为500ps，时间步长设置为1fs。在截断半径为每个分子动力学过程的模拟时间设置为500ps，时间步长设置为1fs。在截断半径为下计算了范德华(vdw)相互作用，采用ewald求和方法在精度为0.001kcal/mol下计算了静电相互作用。模拟过程的能量偏差限制在50000kcal/mol。最后通过对优化后的无定形单元进行单点能计算，可以获得特定温度下油析体系的总能量。同一组成下，利用所有温度的模拟结果即可绘制总能量-温度曲线。
42.(7)按照图1所示的流程，下面以特定体系下的总能量-温度曲线为基础，结合剪应力自相关函数分析、径向分布函数分析和表面能计算的结果，通过分析特定浓度下的总能量-温度曲线即可推测出油析发生的温度。例如，对于吡唑醚菌酯在异丙醇和环己烷混合溶剂中(图4)，在ab段，体系总能量主要由分子的内能构成，其随着温度的下降而降低；在bc段，此时发生了液液相分离，体系内溶质分子二聚体构型发生转化，溶质分子之间的相互作
用强度下降，同时油滴表面能快速增加，导致体系总能量随温度的降低而逐渐升高；在cd段，由于此时溶质之间的相互作用强度和油析表面能都趋于稳定，因此体系总能量再次由分子具有的内能所主导，且随着温度的下降而降低。因此，总能量-温度曲线的局部最低点对应的温度就是当前浓度下发生油析现象的温度。
43.(8)按照图1所示的流程，将不同浓度下的油析温度点相连接即可得到吡唑醚菌酯、香兰素、乙基香兰素、柠檬酸和苯甲酸在相应溶剂体系的油析曲线。
44.(9)按照图1所示的流程，发生油析的溶剂体系和油析曲线即可组成油析热力学相图(图5吡唑醚菌酯(环己烷 异丙醇)，图6中，(a)香兰素(水)、(b)乙基香兰素(水)、(c)柠檬酸(水)和(d)苯甲酸(水 异丙醇))，可以看出其平均相对偏差(ard)值和均方根误差(rmsd)值均较小，说明预测结果和实验结果匹配程度较好，该方法预测结果较为准确。
45.以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。