一种基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法

本发明属于生物医用材料领域，具体涉及一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法。

背景技术：

在生物医用材料领域中，生物活性玻璃(bioactiveglasses,bg)因其具有优良的生物相容性及生物活性已经被广泛应用于组织修复与再生等临床治疗中，因此其生物活性机理的探索成为研究者们不可避免的问题。目前，研究者已通过各种生物学方法研究了bg在生物体中的作用机理，但是在材料结构方面，由于当前对于非晶体材料检测技术的局限性，计算机模拟已经成为研究无定形材料结构和性能的强大工具。

了解生物活性玻璃的结构至关重要，目的是为了更好的解释及优化实验。而且，通过计算机进行分子动力学模拟(md)可得到玻璃结构的详细内容，以及关于复杂玻璃的结构与性能关系的重要信息。由于生物活性强烈依赖于离子物质在生理环境中的释放，所以直接预测玻璃的生物活性反应的先决条件是在原子水平上正确理解构成玻璃的组分结构，以及当改变玻璃组合物时发生在(亚)纳米尺度上的变化。利用分子动力学模拟技术，可以从原子层面对bg结构及表面特征进行微观层面深入解释及表征，并且能合理地设计用于生物材料的新型玻璃复合物以及实现bg与生物大分子的反应过程，以更准确的指导实验研究，从而降低研发成本。

现有的bg原子层面的结构可使用中子衍射、魔角旋转固态核磁共振技术以及高能x射线同步辐射等检测技术进行研究，测试其原子之间键长、配位、网络连接性等，但是这些测试昂贵，难预约，周期长，最后只能得到结构的单一性质。但是分子动力学模拟技术在同样时间内不仅可以得到可靠的bg结构，而且可以进行热力学、动力学等多种性能及结构的分析。

技术实现要素：

为了克服现有技术存在的不足，本发明的目的是提供一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法。

本发明提供了一种基于分子动力学的生物活性玻璃烧结及xrd计算的模拟方法。该方法采用分子动力学模拟，选择合适的势函数，设定合理的模拟参数和方法步骤，直接得到生物活性玻璃结构的坐标文件与xrd数据，通过与实验对比，并对生物活性玻璃进行近程及中程结构分析。

首先，本发明可以为深入研究某一种生物活性玻璃微观结构及表面特征来对实验做出合理解释及表征提供一种方法；其次，可以根据此方法构建不同组分生物活性玻璃进行结构及性能预测来直接进行实验，避免大量实验开销；另外，也可以通过此方法来深入研究生物活性玻璃的结构及离子扩散性能与生物活性的联系等。

本发明的目的至少通过如下技术方案之一实现。

本发明提供的基于分子动力学的生物活性玻璃烧结及xrd计算的模拟方法，包括：选取能够反映生物活性玻璃体系所包含的各原子之间相互作用力的势函数；构建生物活性玻璃初始无定形结构模型并产生分子动力学模拟软件可识别的数据文件；确定分子动力学模拟计算过程参数及条件并写入分子动力学模拟软件中；通过分子动力学模拟软件计算并输出生物活性玻璃结构模型坐标文件及xrd数据；通过xrd数据与实验对比证明结构的可靠性，结合可视化和数据分析软件近程、中程结构分析得到生物活性玻璃微观结构信息。

本发明提供的基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法，具体包括以下步骤：

(1)选取待测生物活性玻璃并测试密度，规定该生物活性玻璃的原子总数，通过密度和原子总数得到模拟立方盒子的边长，使用python脚本生成包含待测生物活性玻璃的各原子初始坐标的data文件，即待测生物活性玻璃初始结构的data文件；

(2)将步骤(1)所述待测生物活性玻璃初始结构的data文件和势函数导入执行lammps软件运行命令的in文件中，在lammps软件将待测生物活性玻璃初始结构进行能量最小化，在lammps软件中对待测生物活性玻璃初始结构升温进行模拟加热处理，然后降温至常温，在模拟加热和降温后，收集平衡结构数据，进行xrd计算，得到xrd图；

(3)将步骤(2)所述平衡结构数据输入可视化软件vmd中，得到待测生物活性玻璃的网络结构图和原子结构图并进行结构分析，完成生物活性玻璃结构及xrd计算的模拟。

进一步地，步骤(1)所述生物活性玻璃初始结构在一个模拟立方盒子中建立，每个原子初始坐标均位于模拟立方盒子内，根据实验测试得到该生物活性玻璃的密度，根据盒子边长及lammps对输入模型格式的要求编写python脚本，建立生物活性玻璃的初始无定形结构，并产生lammps可识别的模型输入数据。所述将生物活性玻璃初始结构限定于一个立方盒子中，根据实验测试的生物活性玻璃密度，计算盒子边长。

进一步地，步骤(1)所述包含待测生物活性玻璃的各原子初始坐标的data文件中，还包含边界条件、模拟系综、温控方法、降温速率、求和方法及xrd计算参数的信息。所述data文件为lammps输入文件。模拟过程计算条件也需要写成lammps输入in文件。所述data文件为根据盒子边长及lammps对输入模型格式的要求编写的python脚本，建立生物活性玻璃的初始无定形结构，并产生lammps可识别的模型输入数据。

进一步地，所述边界条件采用周期性边界条件，选择正则系综(nvt)及ewald求和方法，时间步设为2fs，采用nose-hoover热浴法来调控nvt正则系综的温度，在5000-6000k弛豫后，在1-5k/ps之间选择最佳降温速率降至300k。在高温下充分弛豫100ps以平衡体系，降温至300k以后先在nvt系综下运行，接着在微正则系综(nve)下进行充分弛豫以消除残余应力达到能量最小化，并收集数据用来结果分析。

进一步地，步骤(2)所述模拟加热处理的温度为5000-6000k，所述降温的速率为1-5k/ps。

进一步地，步骤(2)所述势函数为morse势函数，所述势函数的表达式为

其中，uij是距离为rij的所有原子i和j的总势能；i和j表示体系中任意两个原子；zi和zj表示有效电荷；e表示元电荷(e＝1.6×10^-19c)；rij表示i和j两个原子之间的距离；dij为原子i和j的键离解能；aij是势能阱斜率的函数；r0是键平衡距离；cij为弹性常数。

进一步地，步骤(2)所述势函数为描述长程相互作用力的库仑势、描述对硅基玻璃的短程相互作用力的morse势函数(morse项)及高温下防止原子融合的互斥项组成

进一步地，步骤(3)所述xrd计算的过程中，根据需求设定2theta范围以计算xrd，平衡结构后运行computexrd命令以完成xrd计算。

进一步地，步骤(3)中，将平衡结构数据的xyz格式文件输入可视化软件vmd中，观察结构，并通过径向分布函数和qⁿ等近程及中程结构分析方法来判断结构与其性能的关系。qⁿ是玻璃中si或p与桥氧的连接占比，n为0-4。

本发明通过lammps运行分子动力学模拟计算，并得到生物活性玻璃模拟结构后的平衡结构模型的坐标数据文件及xrd数据。通过数据处理软件将xrd数据作图与实验进行对比，以验证结构的可靠性。通过vmd可视化软件打开模拟坐标数据观察结构，并通过近程及中程结构来分析结构与其性能的关系。

本发明的方法通过编写脚本构建分子动力学模拟软件可识别的数据文件，以建立生物活性玻璃无定形结构的初始模型；选取计算生物活性玻璃所包含的各原子之间相互作用力的势函数；设定分子动力学模拟结构的计算参数，通过分子动力学模拟软件运行计算，包括1)高温结构弛豫，2)按指定速率连续降温，3)低温结构平衡，4)计算平衡结构的xrd；输出生物活性玻璃模型的原子坐标文件及xrd数据；通过xrd数据与实验对比证明结构的可靠性，结合可视化及数据分析软件进行近程、中程结构分析得到生物活性玻璃微观结构信息。

本发明提供的方法中，分子动力学模拟结果可以计算得到结构因子；计算中所述的x射线衍射(xrd)强度的计算由整个模拟域定义的倒易格点网格上使用特定波长射线模拟辐射，在各倒易格点处的x射线衍射强度可由结构因子计算，由此获得模拟结构的xrd数据；将xrd数据与实验进行对比，可进一步验证本发明的方法可行性。

与现有技术相比，本发明具有如下优点和有益效果：

(1)本发明提供了一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法，该方法能够采用lammps进行分子动力学模拟出生物活性玻璃结构及xrd数据，并通过xrd与实验对比验证结构可靠性。

(2)本发明采用的生物活性玻璃结构表征方法与传统实验表征相比，可从原子尺度详细剖析生物活性玻璃近程及中程结构特点与生物活性之间的关系，以进一步解释并优化实验，并对生物活性玻璃批量化生产进行相关理论指导。

(3)可以根据本发明提供的方法构建不同组分的生物活性玻璃进行结构性能预测来进行针对性实验，避免昂贵的实验开销。

(4)通过本发明提供的方法构建的结构模型，可以进行生物活性玻璃的离子扩散性能及表面结构研究，对实验结果做出合理解释。

附图说明

图1是本发明一典型实施案例中一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法的流程图；

图2是本发明一典型实施案例中模拟生物活性玻璃结构的xrd(58s-md)与实验xrd(58s-exp)的对比图；

图3是本发明一典型实施案例中的经过弛豫平衡后的58s生物活性玻璃无定形网络结构示意图；

图4是本发明一典型实施案例中的经过弛豫平衡后的58s生物活性玻璃原子结构示意图。

具体实施方式

以下结合实例对本发明的具体实施作进一步说明，但本发明的实施和保护不限于此。需指出的是，以下若有未特别详细说明之过程，均是本领域技术人员可参照现有技术实现或理解的。所用试剂或仪器未注明生产厂商者，视为可以通过市售购买得到的常规产品。

本发明实施例的一个方面提供了基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法，该方法流程图如图1所示，其包括：

步骤一：选取能够反映生物活性玻璃体系所包含的各原子之间相互作用力的势函数；

步骤二：构建生物活性玻璃初始无定形结构模型并产生分子动力学模拟软件可识别的数据文件；

步骤三：确定分子动力学模拟计算过程参数及条件并写入分子动力学模拟软件中；

步骤四：通过分子动力学模拟软件计算并输出生物活性玻璃结构模型坐标文件及xrd数据。

步骤五：通过xrd数据与实验对比证明结构的可靠性，结合可视化和数据分析软件近程、中程结构分析得到生物活性玻璃微观结构信息。

步骤一中，选取的势函数由描述长程相互作用力的库仑势与对硅基玻璃的短程相互作用力有良好描述的morse势函数，以及高温下防止原子融合的互斥项组成。

所以，描述原子间作用的势函数为：

其中，uij是距离为rij的所有原子i和j的总势能；i和j表示体系中任意两个原子；zi和zj表示有效电荷；e表示元电荷(e＝1.6×10^-19c)；rij表示i和j两个原子之间的距离；dij为原子i和j的键离解能；aij是势能阱斜率的函数；r0是键平衡距离；cij为弹性常数。

步骤二中，生物活性玻璃的初始结构模型置于一个立方盒子中，根据实验测试的生物活性玻璃密度,计算盒子边长。

编写python脚本，根据盒子边长建立生物活性玻璃的初始无定形结构，使每个原子初始坐标位于盒子内，并产生lammps可识别的data格式数据。

步骤三中，模拟过程的条件及参数通过lammps模拟控制文件in文件的编写输入，主要包括：边界条件，模拟系综，温控方法，降温速率，求和方法，xrd计算参数。

在一些实例中，边界条件采用周期性边界条件，选择正则系综(nvt)及ewald求和方法，时间步设为2fs，采用nose-hoover热浴法来调控nvt的温度，在6000k弛豫后，以5k/ps的降温速率降至300k。

在一些实例中，在6000k充分弛豫100ps以平衡体系，降温至300k以后先在nvt系综下运行50ps，接着在微正则系综(nve)下运行50ps进行充分弛豫以消除残余应力达到能量最小化，并收集数据用来结果分析。

在一些实例中，xrd计算参数设定2theta范围为10-70度，使用lammps语句computexrd命令，平衡结构后运行1000步以完成xrd计算。

步骤四中，采用lammps运行in文件进行分子动力学模拟计算，得到生物活性玻璃模拟的平衡结构模型的xyz格式坐标文件及包含xrd数据的xrd文件。

步骤五中，通过数据处理软件分析处理xrd数据与实验进行对比，进行基本的结构可靠性验证，需要考虑模拟是在真空条件下进行，xrd不会有任何杂峰，与结构成分相关主峰符合即可。

步骤五中，通过vmd可视化软件打开模拟结构坐标数据观察结构，并通过数据处理软件isaacs对生物活性玻璃的近程及中程结构进行处理分析，由此结合实验解析其结构与性能的关系。

实施例1

本实施实例使用了一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法进行实验，主要包括以下步骤：

(1)首先确定一种成分为60％sio2-36％cao-4％p2o5(mol％)的生物活性玻璃，简称58s，设计模拟采用7560原子，根据实验密度2.67g/cm³计算得到模拟包含7560原子的立方盒子边长根据lammps对输入的结构数据文件格式要求，使用python脚本产生58s结构模型data文件；

(2)将已确定的morse势函数参数及相关计算条件写入执行lammps运行命令的in文件，对于高温互斥项采用lammps中的线性表引入，然后对初始结构进行能量最小化，在nvt系综下，在立方盒子的xyz方向设置周期性边界条件，对初始58s结构进行加热至5000k，时间步为2fs，并在5000k下恒温运行100ps，以保证结构充分弛豫，再以5k/ps的速率降温至300k，并进行50ps弛豫，最后在nve系综下运行50ps收集平衡结构数据用来数据分析，最后使用computexrd命令，并设置扫描角度10-70度，进行xrd计算。

(3)如图2所示，对xrd处理后作图，与实验测试对比，其玻璃无定形结构主峰21-30度与实验相符，实验测试32度左右结晶峰是由于样品反应活性强，与空气中co2和水反应形成羟基磷灰石(hap)结晶，与样品组成无关，因此模拟产生的58s结构是可靠的，可进一步进行结构分析。

(4)通过可视化软件vmd打开xyz格式的58s结构数据文件，可看到如图3所示的58s网络结构图及图4所示的原子结构图，将xyz格式文件导入isaacs软件即可进行径向分布函数、键长、键角分布、环分布统计等近程及中程范围结构分析，得到实验测试无法实现的结构剖析，进一步深入解释生物活性机理，并对生物活性玻璃的实验设计做出指导。

实施例2

本实施实例使用了一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法进行实验，主要包括以下步骤(可参照图1所示)：

(1)首先确定一种成分为46.1sio2,-24.4na2o-26.9cao-2.6p2o5(mol％)的生物活性玻璃，简称45s5，设计模拟采用4992原子，根据实验密度2.72g/cm³计算得到模拟包含4992原子的立方盒子边长根据lammps对输入的结构数据文件格式要求，使用python脚本产生45s5结构模型data文件；

(2)将已确定的morse势函数参数及相关计算条件写入执行lammps运行命令的in文件，对于高温互斥项采用lammps中的线性表引入，然后对初始结构进行能量最小化，在nvt系综下，在立方盒子的xyz方向设置周期性边界条件，对初始45s5结构进行加热至6000k，时间步为2fs，并在6000k下恒温运行100ps，以保证结构充分弛豫，再以5k/ps的速率降温至300k，并进行50ps弛豫，最后在nve系综下运行50ps收集平衡结构数据用来数据分析，最后使用computexrd命令，并设置扫描角度10-70度，进行xrd计算。

(3)将xrd与实验测试对比，其玻璃结构主峰与实验相符，因此模拟产生的45s5结构是可靠的，可进一步进行结构分析。

(4)通过可视化软件vmd打开xyz格式的45s5结构数据文件，将xyz格式文件导入isaacs软件即可进行径向分布函数、键长、键角分布、环分布统计等近程及中程范围结构分析，进一步深入解释生物活性机理，并对生物活性玻璃的实验结果与预测做出指导。

实施例3

本实施实例使用了一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法进行实验，主要包括以下步骤(可参照图1所示)：

(1)首先确定一种成分为75％sio2-15％cao-10％na2o(mol％)的无活性生物玻璃，设计模拟采用114原子，根据实验密度2.5g/cm³计算得到模拟包含114原子的立方盒子边长根据lammps对输入的结构数据文件格式要求，使用python脚本产生其结构模型data文件；

(2)将已确定的morse势函数参数及相关计算条件写入执行lammps运行命令的in文件，对于高温互斥项采用lammps中的线性表引入，然后对初始结构进行能量最小化，在nvt系综下，在立方盒子的xyz方向设置周期性边界条件，对初始结构进行加热至6000k，时间步为2fs，并在6000k下恒温运行100ps，以保证结构充分弛豫，再以2.5k/ps的速率降温至300k，并进行50ps弛豫，最后在nve系综下运行50ps收集平衡结构数据用来数据分析，最后使用computexrd命令，并设置扫描角度10-70度，进行xrd计算。

(3)将xrd与实验测试对比，其玻璃结构主峰与实验一致，因此模拟产生的结构可靠，可进一步进行结构分析，判断无活性与活性玻璃的结构区别。

(4)通过可视化软件vmd打开xyz格式的结构数据文件，将xyz格式文件导入isaacs软件即可进行径向分布函数、键长、键角分布、环分布统计等近程及中程范围结构分析，作为实验的部分参考依据。

实施例4

本实施实例使用了一种基于分子动力学的生物活性玻璃结构及xrd计算的模拟方法进行实验，主要包括以下步骤(可参照图1所示)：

(1)首先确定一种成分为70％sio2-30％cao(mol％)的生物活性玻璃，设计模拟采用2798原子，根据实验密度2.58g/cm³计算得到模拟包含原子的立方盒子边长根据lammps对输入的结构数据文件格式要求，使用python脚本产生结构模型data文件；

(2)将已确定的morse势函数参数及相关计算条件写入执行lammps运行命令的in文件，对于高温互斥项采用lammps中的线性表引入，然后对初始结构进行能量最小化，在nvt系综下，在立方盒子的xyz方向设置周期性边界条件，对初始结构进行加热至6000k，时间步为2fs，并在6000k下恒温运行100ps，以保证结构充分弛豫，再以1k/ps的速率降温至300k，并进行50ps弛豫，最后在nve系综下运行50ps收集平衡结构数据用来数据分析，最后使用computexrd命令，并设置扫描角度10-70度，进行xrd计算。

(3)将xrd与实验测试对比，忽略真实实验环境带来的误差，其玻璃结构主峰与实验结果基本相符，因此模拟产生的玻璃结构是可靠的，可进一步进行结构分析。

(4)通过可视化软件vmd打开xyz格式的结构数据文件，将xyz格式文件导入isaacs软件即可进行径向分布函数、键长、键角分布、环分布统计等近程及中程范围结构分析，进一步深入解释二元体系bg的生物活性机理，并指导实验进行组分改进。

以上实施例仅为本发明较优的实施方式，仅用于解释本发明，而非限制本发明，本领域技术人员在未脱离本发明精神实质下所作的改变、替换、修饰等均应属于本发明的保护范围。