【计算教程011期】计算气体在聚合物中的扩散率

本期尝试翻译Forcite的一个教程，个人不建议用该模块做动力学模拟，难道Gromacs不香吗？当然啦，优点还是有的，点点鼠标就能搞计算。此次翻译，还是口语化翻译，相关的专业词汇我尽量给英文，如果还是有让大家费解的说法存在可以参考英文原版教程。

这次先放上半部分，下期推怎么分析动力学模拟的结果，以及如何根据计算结果计算扩散系数。

背景：气体在有机溶剂，聚合物或沸石中的扩散率可通过运行分子动力学模拟并确定材料中气体的均方位移来计算。这可以让你计算气体的自扩散系数，并进而可以研究全扩散系数。当你进行分子动力学计算的时候，你可以分析温度，压力，密度，渗透剂尺寸和结构对扩散的影响。

在本教程中，你将通过构造一个包含甲烷和聚（顺-1,4-丁二烯）poly(cis-1,4-butadiene)PBD的无定形晶胞来计算甲烷在PBD中的扩散率。构建完晶胞后，您将执行分子动力学模拟并计算甲烷分子的均方位移。

尽管本教程仅演示了一个简短的计算过程，但是您将熟悉所涉及的方法。该教程基于Meunier发表的论文(https://aip.scitation.org/doi/10.1063/1.2049274 )，该论文研究了气体在二烯聚合物中的扩散。

1 建立一个叫做gas_polymer的project

这步不说了。

2 搭建初始结构

第一步是构建并优化甲烷分子和PBD聚合物来构建无定形原胞。为了便于以后选择甲烷分子，请使用其他显示样式。

使用均聚物Homopolymer构建工具（Build | Build Polymers | Homopolymer）从dienes library中生成c_butadiene的20个重复单元聚合物。

新建一个3D Atomistic Document并绘制一个甲烷分子。将原子的显示样式更改为CPK。重命名为methane。这个也不细说了，之前的教程讲过怎么画的。

一个经验力场forcefield计算（能量最小化或分子动力学）中花费最大的部分是非键参数的确定（库仑相互作用和范德华力）。涉及力场的计算会用各种方法来计算非键参数。说句不负责任的话，力场开发者大部分头发都在非键参数计算上脱落的。对范德华力默认的方法是原子间模拟，对库仑相互作用则是Ewald方法（可以精确计算长程静电作用，但是贵啊）模拟。

对某些聚合物，可以用一组原子而不是单个原子来逼近非键参数。这种方法叫做电荷组charge groups。每个电荷组包含数个原子，里面所有原子电荷加和为整数。利用了charge groups概念可以来降低cut-off方式计算静电作用的误差。这种方法可以在不损害精度的情况下加速计算

在本教程中，我们将在静电和范德华相互作用的计算中使用charge groups，通过基于距离的直接截断来评估两组之间的静电相互作用，从而避免了计算上更昂贵的Ewald方法。

由于电荷基团是在重复单元上定义的，因此它们可自动用于聚合物。您可以通过按charge group着色来验证这一点。

使Polyc_butadiene.xsd成为活动文档。右键单击3D Viewer，然后从快捷菜单中选择“ Display Style ”以打开“Display Style”对话框。将颜色按charge group着色。验证后，再改回来。

稍后将使用力场分配的电荷，假定这些电荷对每个组总计为零。否则，计算将失败，必须重新计算费用组电荷组。

画好的分子上不存在电荷基团，因此必须首先计算甲烷分子的电荷基团，并将Forcite设置为使用电荷组而不是默认的求和方法进行非键作用计算。

先将methane.xsd设置为活动文档，然后打开Forcite计算设置(Modules | Forcite | Calculation)

在“ Energy”选项卡上，从“ Forcefield” 下拉列表中选择“ COMPASS ” 。将“ Electrostatic”和“ van der Waals”求和方法都设置为“ Group based”。

单击“ Forcefield”的“ More... ”按钮以打开“ Forcite Preparation Options”对话框。取消选中“  Forcefield types”的“ Calculate automatically ”复选框，然后单击“ Calculate” 按钮。将“ Charges ”设置为“ Forcefield assigned ”，然后单击“Calculate”按钮。

单击“ Charge groups”的“ More...”按钮以打开“Forcite Charge Groups”对话框，然后单击“ Calculate”按钮。在“ Forcite Preparation Options”对话框中，选中“Calculate automatically”复选框，然后关闭对话框。

在构建非晶晶胞（无定形晶胞）amorphous cell之前，先要优化两个分子的几何形状。在“ Forcite Calculation”对话框的“ Setup”选项卡上，从“ Task”下拉列表中选择“ Geometry Optimization ” 。更改Max. iterations到1000。

先优化甲烷分子

点击methane.xsd成为活动文档，然后单击Forcite calculation对话框Run按钮。

计算很快完成，项目资源管理器中生成一个methane Forcite GeomOpt文件夹，计算完成后，最小化的结构将存储在新文件夹中。

然后优化聚合物

使Polyc_butadiene.xsd成为活动文档，然后单击“ Run” 按钮。关闭“Forcite calculation”对话框。

重复相同的过程，最小化的结构将被存储到Polyc_butadiene Forcite GeomOptPolyc_butadiene.xsd。

聚合物中的扭转自由度将通过Amorphous Cell进行修改。这些将在构建晶胞之后进行优化。

在继续进行非晶晶胞的构造之前，先清理工作空间。File | Save Project然后Window | Close All 。

3 建立非晶晶胞

一旦准备好两个结构，就可以使用Amorphous Cell模块在一个晶胞中构建它们的多个副本。单击“ Modules” 工具栏上的“ Amorphous Cell”按钮，然后从下拉列表中选择“ Calculation ”。这将打开“Amorphous Cell Calculation”对话框。

第一步是根据每种组分的分子数定义组成。

此次需要构建的晶胞是四个甲烷分子和十个PBD组成的，密度为0.95 g/cm3。

将Density设置为0.95  g/cm3。

在“ Composition ”格子的“ Molecule”列中，选择包含优化的甲烷结构的甲烷文件 Forcite GeomOptmethane.xsd。在“Loading”列中，输入4。

在下一行中，选择Polyc_butadiene Forcite GeomOptPolyc_butadiene.xsd并加载10个。

根据载荷和密度预估的晶胞尺寸将出现在对话框底部。

在这种情况下，将构建一个晶胞长度约为27 Å的Cubic。Tetragonal和Orthorhombic类型也可用，但在本教程中未使用。

Amorphous Cell可以优化结构，作为构造的一部分。我们将使用Forcite分别进行优化和平衡，而不使用Amorphous Cell。

单击Options...按钮以打开“Amorphous Cell Options”对话框。取消选中Optimize geometry。关闭对话框。

选择和Forcite一样的力场。在“ Amorphous Cell Calculation”对话框的“ Energy”选项卡上，从“ Forcefield”下拉列表中选择“ COMPASS ” 。

Amorphous Cell中的默认作业描述与第一个组份的名称相对应，在本例中为methane，它在所有输出文档中均用作种子名称。在本教程中，您将默认设置更改为cell。设置方法如下图。

在Project Explorer中点击gas_polymer，然后点击“Run”按钮。很快生成完毕，如下图。在Project Explorer中生成一个新文件夹cell AC Construct，里面的cell.xtd就是生成包含非晶晶胞的轨迹文档。

进行接下来的操作前，需要把xtd转化为xsd。具体的操作方法是：首先在project根目录下新建一个cell.xsd，然后双击上面做好的cell.xtd，在3D viewer界面里面右击，点击copy，然后双击cell.xsd，在右边的3D viewer窗口中右击点击paste。

保存project，然后关闭所有窗口。

4 晶胞的弛豫

当一个无定形晶胞生成时，分子可能不会在整个晶胞中均匀分布，这样就造成了低密度真空区。为了矫正这个问题，要进行能量最小化来优化晶胞。最小化过后，要进行分子动力学模拟来平衡晶胞。当你构建无定形晶胞时，都要用能量最小化和分子动力学来进行结构弛豫。PS 在能量最小化之前，清空工作区。

要执行几何优化，必须首先将Forcite设置为对3D周期结构使用电荷组。打开cell.xsd。打开“ Forcite Calculation”对话框，任务类型选择“Geometry Optimization”然后选择“Energy”选项卡。将Electrostatic和van der Waals改为Group Based。然后单击“Run”按钮。

计算完成后，最终结构存储在cell Forcite GeomOpt文件夹中。

然后要进行弛豫，就是在周期性变化的温度下运行分子动力学来继续弛豫结构，即对系统退火。本教程而言，我们将只运行一个退火周期。具体操作如下：

在“ Forcite Calculation”对话框的“ Setup”选项卡上，从“ Task”下拉列表中选择“ Anneal ” ，然后单击“ More...”按钮以打开“Forcite Anneal Dynamics”对话框。 将Annealing cycles数设置为1，将Initial temperature设置为300  K，将Mid-cycle temperature设置为500K  。关闭对话框。

现在，点击cell.xsd使其为激活状态，点击“ Forcite Calculation”对话框的“Run”对优化的结构执行退火计算。

计算完成后，生成cell Forcite Anneal文件夹，退火任务的最终结构是cell.xsd。感兴趣的话，可以双击cell.xtd然后点击play键观看动画。如果找不到这个键，可以View-Toolbars-Animation。

然后再在此结构上做一个恒温动力学模拟。

先把cell Forcite Anneal文件夹中的cell.xsd点一下变成激活状态。在“ Forcite Calculation”对话框的“ Setup”选项卡上，从“ Task”下拉列表中选择“ Dynamics ” ，然后单击“ More...”按钮以打开“Forcite Dynamics”对话框。

首先是要设定系综(ensemble)，强烈建议找找相应百科学习一下，这个是动力学上面很重要的概念。有各种不同的分子动力学模拟，以系综分类，分别为 NVE, NVT, NPT, 和 NPH。这些字母是指：

N=固定粒子数

V=固定体积

E=固定能量

T=固定温度

P=固定压强

H=固定焓

如果在构造时选择的密度（体积）需要调整到外部压力（通常是大气压），则应使用NPT系综； 如果系统以合理的密度构建，则可以使用NVT，并且压力平均应为1 atm或0.0001 GPa。要平衡一个准备进行扩散计算的晶胞，NPT 系综是最好的选择。不过，本教程中采用最快的 NVT 系综。

此前以0.95 g / cm 3的密度构建了该单元晶胞。由于这也是该系统在300 K和1 atm的平均密度，并且具有选定的力场，因此无需使用NPT进一步弛豫密度。可以使用NVT进行动力学模拟。

从“ Ensemble”下拉列表中选择NVT，然后将“ Temperature”更改为300。

本教程而言，我们将步骤数减少到2000。对于这样的短时间模拟，Velocity Scale温控器比默认温控器更合适。然后单击Forcite Calculation对话框上的“Run”开始运行。

这里只跑了2 ps真的只是演示，没人实际计算时候跑这么短，一般都是几ns了。至少运行50 ps才能正确平衡晶胞。您可以通过在实时更新图表中查找能量来监视平衡进度，该能量应保持恒定（除了小波动）。

模拟完成后，将返回许多文档。cell Forcite Dynamics文件夹中的cell.xsd是最终结构。

然后保存project并关闭窗口。

---

动动小手加星标，浏览文章不迷路！不用每天花费时间刷信息流也可以随时看到自己喜欢的内容啦！

往期推荐
万事屋告示牌
招聘话题小编

⭐我们希望你

擅长写各类硕博热点话题，内容有趣，有料，有思想

⭐简历投递

1. 投递邮箱：[e-mail]2.扫描左侧海报中的二维码，添加负责人微信
关注我们

点了“在看”的小哥哥小姐姐

今年发IF>10一作