我翻译的教程是根据MS8.0的官方教程来写的。我发现新版的尤其是2018版以后的,在菜单上面可能会有一些差异,大家可以灵活变通一下。
MD计算是特别烧钱的,大部分钱花在机时费上面了。我们学习时候都是尽可能简化条件缩短时间,大家要理性面对这种演示性计算中得到违反常规认知结果的情况。另外,建议没有学明白的情况下,不要贸然用文章中的设置水平算自己的体系,容易浪费机时。
上期内容:【计算教程011期】计算气体在聚合物中的扩散率
平衡系统时,我们只对最终结构感兴趣。但是,要计算胞中甲烷分子的均方位移,我们需要许多帧frames,以便可以分析甲烷分子的移动位置。所以我们还要运行另一个分子动力学模拟,并生成一个轨迹文档,然后可以使用Forcite Analysis tool工具进行分析。为了避免过多的子文件夹,请首先将工作文档移动到文件夹树的顶部。直接把cell Forcite Dynamics文件夹中的cell.xsd拖到gas_polymer根目录下。之前,你运行了一个NVT系综,不过最好用NVE系综。因为就方法而言,NVE 动力学不会被系统的热力学过程干扰,那些干扰可能会影响稍后要计算的扩散系数。在Forcite Dynamics对话框的Dynamics选项卡上,从Ensemble下拉列表中选择NVE 。为了收集足够的数据进行分析,应该增加步数并减少帧输出间隔。就本教程而言,需要执行 5000步。把Number of steps 改为5000。把Frame output every改为250。关闭Forcite Dynamics对话框。在Forcite Calculation对话框中,单击Run按钮,然后关闭对话框。特别提醒,实际计算时候(有别于为了写软件随便鼓捣),模拟时间上最少50 ps。随着计算的进行,将更新两个图表文档。一个绘制总能量和各种成分随时间变化的曲线,另一个则是温度随时间的变化。因为这是NVE 系综,因此总能量应该是恒定的。不过温度会有涨落,将在300 K左右的平均值附近波动,直至收敛到目标温度。时间模拟太短,这部分很难看,计算实际课题时候模拟长一些。计算完成后,就可以开始分析输出文件了。cell Forcite Dynamics文件夹中包含了一个21帧的轨迹cell.xtd。双击cell.xtd然后点击Animation工具栏的播放按钮就可以观看动画。看完后可以点击 Stop按钮停止播放。要计算甲烷分子的均方位移,需要将其与聚合物分子区分开。这可以通过将它们定义为一组来实现。为了计算氧分子的均方位移,你要把它们同聚合物分子区分开来。这可以通过把它们定义成一组来达到。要选择所有甲烷分子,请按住CTRL键, 然后依次双击每个分子。将所有的甲烷分子都选上后,点击Edit | Edit Sets打开Edit Sets对话框,点击New…然后在弹出的新对话框中输入methane并点击OK。现在,已将甲烷分子定义为一组,可以分析它们的运动了。点击Modules | Forcite | Analysis打开Forcite Analysis对话框。从上图可以看出,Forcite可以做很多不同类型的分析,主要分为三类Structural, Statistics和 Dynamic。均方位移Mean square displacement在Dynamic分类中。从列表中选择Mean square displacement然后设置Sets为methane,设置Length为21。然后点击Analyze按钮并关闭对话框。Forcite Analysis tool计算均方位移并生成图表文档,cell Forcite MSD.xcd(图)和cell Forcite MSD.std(表)其中包含甲烷分子随时间的均方位移信息。在给定时间的均方位移是对所有相同长度的时间段和那个组里的所有原子作平均得到的。这句我翻译比较拗口,长难句让我头秃。大家可以自己看看原文The value of the MSD for a given time reported in the chart is the average over all time intervals of that length and over all atoms in the set.看图,看上一节那个图,可以看到均方位移具有两个阶段。气体分子在短时间仅能在一个很小的自由体积内碰撞。由于分子受到限制,因此它不会在此时间范围内扩散diffuse,MSD趋于恒定。在更长的时间尺度上,分子会跳出约束区域,到达另一个自由体积。重复跳跃的最终运动是扩散,其特征在于时间上呈线性的均方位移。实际上,统计信息会随着时间间隔而减少,通常最终会导致较大的波动。其中Nα是系统中扩散原子的数量。Forcite中的MSD已在甲烷分子中的所有原子上求平均。这个公式看着有点复杂,大部分是MSD啊。为了确定扩散系数,您需要将一条直线y = a * x + b拟合到扩散状态下的数据并提取斜率a。a的单位是Å2/ps。根据上面的定义,D于是可以简化如下:此时可以使用任何外部的表格或者图表工具进行这部分的拟合分析。本教程中使用MS自带的Plot Graph tool。删除cell Forcite MSD.std的初始非线性部分以及最后的嘈杂部分相对应的行。就是选中那几行然后右键点击行标题,然后选择Delete。然后对其余数据执行线性回归分析。选中剩余表格中的A和B列,然后Tools | Plot Graph打开Plot Graph对话框。将Graph type设为Scatter (2-D) 将 Show best fit line勾上。然后点击Plot按钮。系数a可以从图中读取。将此数除以6获得以Å2/ps为单位的扩散系数。要转换为更常用的单位cm 2 /s,须将结果值除以1e4。J. Chem. Phys., 123, 134906 (2005) 报道,在PBD中甲烷的扩散系数为 2.25× 10-6 cm2 s-1 至 7.5× 10-6 cm2 s-1之间。具体参数可以看文章里面写的,看完这个教程,大家应该能够读懂文中的相关参数说明了。您的计算值可能与报告值有很大差异,因为运行时间很短,且扩散区域的统计信息有限。