一种基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法转让专利
申请号 : CN202110247229.5
文献号 : CN113012764B
文献日 : 2022-04-19
发明人 : 陈晓峰 , 何闲 , 王刚
申请人 : 华南理工大学 , 佛山今兰生物科技有限公司
摘要 :
权利要求 :
1.一种基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,包括以下步骤:
(1)选取待测生物活性玻璃并测试密度,规定该生物活性玻璃的原子总数,通过密度和原子总数得到模拟立方盒子的边长,使用Python脚本生成包含待测生物活性玻璃的各原子初始坐标的data文件,即待测生物活性玻璃初始结构的data文件;
(2)将步骤(1)所述待测生物活性玻璃初始结构的data文件和势函数导入执行Lammps软件运行命令的in文件中,在Lammps软件将待测生物活性玻璃初始结构进行能量最小化,在Lammps软件中对待测生物活性玻璃初始结构升温进行模拟加热处理,然后降温至常温,在模拟加热和降温后,收集平衡结构数据,进行XRD计算,得到XRD图;
(3)将步骤(2)所述平衡结构数据输入可视化软件VMD中,得到待测生物活性玻璃的网络结构图和原子结构图并进行结构分析,完成生物活性玻璃结构及XRD计算的模拟。
2.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(1)所述生物活性玻璃初始结构在一个模拟立方盒子中建立,每个原子初始坐标均位于模拟立方盒子内,根据实验测试得到该生物活性玻璃的密度,根据盒子边长及Lammps对输入模型格式的要求编写Python脚本,建立生物活性玻璃的初始无定形结构,并产生Lammps可识别的模型输入数据data文件。
3.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(2)所述执行Lammps软件运行命令的in文件中,还包含边界条件、模拟系综、温控方法、降温速率、求和方法及XRD计算参数的信息。
4.根据权利要求3所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,所述边界条件采用周期性边界条件,所述模拟系综选择NVT正则系综及求和方法选择Ewald求和方法,时间步设为2fs,所述温控方法采用Nose‑Hoover热浴法来调控NVT正则系综的温度。
5.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(2)所述模拟加热处理的温度为5000‑6000k,所述降温的速率为1‑5k/ps。
6.根据权利要求5所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,所述模拟加热处理的温度为6000k,所述降温的速率为5k/ps。
7.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(2)所述势函数为描述长程相互作用力的库仑势、描述对硅基玻璃的短程相互作用力的Morse势函数及高温下防止原子融合的互斥项组成。
8.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(2)所述势函数的表达式为其中,Uij是距离为rij的所有原子i和j的总势能;i和j表示体系中任意两个原子;zi和zj表示有效电荷;e表示元电荷;rij表示i和j两个原子之间的距离;Dij为原子i和j的键离解能;aij是势能阱斜率的函数;r0是键平衡距离;Cij为弹性常数。
9.根据权利要求1所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(2)所述XRD计算的过程中,设定2Theta范围,收集平衡结构数据后运行compute xrd命令以完成XRD计算。
10.根据权利要求1‑9任一项所述的基于分子动力学的生物活性玻璃结构及XRD计算的模拟方法,其特征在于,步骤(3)中,将平衡结构数据的xyz格式文件输入可视化软件VMD中,观察结构;并通过径向分布函数,近程及中程结构分析方法来判断结构与其性能的关系。
说明书 :
一种基于分子动力学的生物活性玻璃结构及XRD计算的模拟
方法
技术领域
背景技术
机理的探索成为研究者们不可避免的问题。目前,研究者已通过各种生物学方法研究了BG
在生物体中的作用机理,但是在材料结构方面,由于当前对于非晶体材料检测技术的局限
性,计算机模拟已经成为研究无定形材料结构和性能的强大工具。
与性能关系的重要信息。由于生物活性强烈依赖于离子物质在生理环境中的释放,所以直
接预测玻璃的生物活性反应的先决条件是在原子水平上正确理解构成玻璃的组分结构,以
及当改变玻璃组合物时发生在(亚)纳米尺度上的变化。利用分子动力学模拟技术,可以从
原子层面对BG结构及表面特征进行微观层面深入解释及表征,并且能合理地设计用于生物
材料的新型玻璃复合物以及实现BG与生物大分子的反应过程,以更准确的指导实验研究,
从而降低研发成本。
测试昂贵,难预约,周期长,最后只能得到结构的单一性质。但是分子动力学模拟技术在同
样时间内不仅可以得到可靠的BG结构,而且可以进行热力学、动力学等多种性能及结构的
分析。
发明内容
到生物活性玻璃结构的坐标文件与XRD数据,通过与实验对比,并对生物活性玻璃进行近程
及中程结构分析。
行结构及性能预测来直接进行实验,避免大量实验开销;另外,也可以通过此方法来深入研
究生物活性玻璃的结构及离子扩散性能与生物活性的联系等。
玻璃初始无定形结构模型并产生分子动力学模拟软件可识别的数据文件;确定分子动力学
模拟计算过程参数及条件并写入分子动力学模拟软件中;通过分子动力学模拟软件计算并
输出生物活性玻璃结构模型坐标文件及XRD数据;通过XRD数据与实验对比证明结构的可靠
性,结合可视化和数据分析软件近程、中程结构分析得到生物活性玻璃微观结构信息。
原子初始坐标的data文件,即待测生物活性玻璃初始结构的data文件;
小化,在Lammps软件中对待测生物活性玻璃初始结构升温进行模拟加热处理,然后降温至
常温,在模拟加热和降温后,收集平衡结构数据,进行XRD计算,得到XRD图;
子边长及Lammps对输入模型格式的要求编写Python脚本,建立生物活性玻璃的初始无定形
结构,并产生Lammps可识别的模型输入数据。所述将生物活性玻璃初始结构限定于一个立
方盒子中,根据实验测试的生物活性玻璃密度,计算盒子边长。
data文件为Lammps输入文件。模拟过程计算条件也需要写成Lammps输入in文件。所述data
文件为根据盒子边长及Lammps对输入模型格式的要求编写的Python脚本,建立生物活性玻
璃的初始无定形结构,并产生Lammps可识别的模型输入数据。
弛豫后,在1‑5k/ps之间选择最佳降温速率降至300k。在高温下充分弛豫100ps以平衡体系,
降温至300k以后先在NVT系综下运行,接着在微正则系综(NVE)下进行充分弛豫以消除残余
应力达到能量最小化,并收集数据用来结果分析。
和zj表示有效电荷;e表示元电荷(e=1.6×10 C);rij表示i和j两个原子之间的距离;Dij为
原子i和j的键离解能;aij是势能阱斜率的函数;r0是键平衡距离;Cij为弹性常数。
察结构,并通过径向分布函数和Q 等近程及中程结构分析方法来判断结构与其性能的关
n
系。Q是玻璃中Si或P与桥氧的连接占比,n为0‑4。
比,以验证结构的可靠性。通过VMD可视化软件打开模拟坐标数据观察结构,并通过近程及
中程结构来分析结构与其性能的关系。
用力的势函数;设定分子动力学模拟结构的计算参数,通过分子动力学模拟软件运行计算,
包括1)高温结构弛豫,2)按指定速率连续降温,3)低温结构平衡,4)计算平衡结构的XRD;输
出生物活性玻璃模型的原子坐标文件及XRD数据;通过XRD数据与实验对比证明结构的可靠
性,结合可视化及数据分析软件进行近程、中程结构分析得到生物活性玻璃微观结构信息。
拟辐射,在各倒易格点处的x射线衍射强度可由结构因子计算,由此获得模拟结构的XRD数
据;将XRD数据与实验进行对比,可进一步验证本发明的方法可行性。
XRD与实验对比验证结构可靠性。
化实验,并对生物活性玻璃批量化生产进行相关理论指导。
附图说明
具体实施方式
现或理解的。所用试剂或仪器未注明生产厂商者,视为可以通过市售购买得到的常规产品。
和zj表示有效电荷;e表示元电荷(e=1.6×10 C);rij表示i和j两个原子之间的距离;Dij为
原子i和j的键离解能;aij是势能阱斜率的函数;r0是键平衡距离;Cij为弹性常数。
的降温速率降至300k。
最小化,并收集数据用来结果分析。
峰符合即可。
性能的关系。
称58S,设计模拟采用7560原子,根据实验密度2.67g/cm计算得到模拟包含7560原子的立
方盒子边长 根据Lammps对输入的结构数据文件格式要求,使用Python脚本产生58S
结构模型data文件;
NVT系综下,在立方盒子的XYZ方向设置周期性边界条件,对初始58S结构进行加热至5000k,
时间步为2fs,并在5000k下恒温运行100ps,以保证结构充分弛豫,再以5k/ps的速率降温至
300k,并进行50ps弛豫,最后在NVE系综下运行50ps收集平衡结构数据用来数据分析,最后
使用compute xrd命令,并设置扫描角度10‑70度,进行XRD计算。
成羟基磷灰石(HAP)结晶,与样品组成无关,因此模拟产生的58S结构是可靠的,可进一步进
行结构分析。
数、键长、键角分布、环分布统计等近程及中程范围结构分析,得到实验测试无法实现的结
构剖析,进一步深入解释生物活性机理,并对生物活性玻璃的实验设计做出指导。
性玻璃,简称45S5,设计模拟采用4992原子,根据实验密度2.72g/cm 计算得到模拟包含
4992原子的立方盒子边长 根据Lammps对输入的结构数据文件格式要求,使用
Python脚本产生45S5结构模型data文件;
NVT系综下,在立方盒子的XYZ方向设置周期性边界条件,对初始45S5结构进行加热至
6000k,时间步为2fs,并在6000k下恒温运行100ps,以保证结构充分弛豫,再以5k/ps的速率
降温至300k,并进行50ps弛豫,最后在NVE系综下运行50ps收集平衡结构数据用来数据分
析,最后使用compute xrd命令,并设置扫描角度10‑70度,进行XRD计算。
析,进一步深入解释生物活性机理,并对生物活性玻璃的实验结果与预测做出指导。
设计模拟采用114原子,根据实验密度2.5g/cm 计算得到模拟包含114原子的立方盒子边长
根据Lammps对输入的结构数据文件格式要求,使用Python脚本产生其结构模型
data文件;
NVT系综下,在立方盒子的XYZ方向设置周期性边界条件,对初始结构进行加热至6000k,时
间步为2fs,并在6000k下恒温运行100ps,以保证结构充分弛豫,再以2.5k/ps的速率降温至
300k,并进行50ps弛豫,最后在NVE系综下运行50ps收集平衡结构数据用来数据分析,最后
使用compute xrd命令,并设置扫描角度10‑70度,进行XRD计算。
为实验的部分参考依据。
用2798原子,根据实验密度2.58g/cm 计算得到模拟包含原子的立方盒子边长 根
据Lammps对输入的结构数据文件格式要求,使用Python脚本产生结构模型data文件;
NVT系综下,在立方盒子的XYZ方向设置周期性边界条件,对初始结构进行加热至6000k,时
间步为2fs,并在6000k下恒温运行100ps,以保证结构充分弛豫,再以1k/ps的速率降温至
300k,并进行50ps弛豫,最后在NVE系综下运行50ps收集平衡结构数据用来数据分析,最后
使用compute xrd命令,并设置扫描角度10‑70度,进行XRD计算。
一步深入解释二元体系BG的生物活性机理,并指导实验进行组分改进。
护范围。