这是本文档旧的修订版！

ReaxFF-GCMC：巨正则系综蒙特卡洛模拟

应用案例：

本教程以Si表面的氧化为例。

GCMC是AMS驱动的，各个模块使用起来，除了各个方法（分子体系DFT、周期边界条件DFT、MOPAC、DFTB等）本身的参数（也就是Main窗口的设置）之外，其他参数设置是一样的。因此其他模块的GCMC模拟，也可以参考ReaxFF的GCMC模拟教程。

建模

View - Periodic - Show Unit Cell:

点击底部四边形按钮，切出表面：

Edit - Crystal - Generate Super Cell，生成超胞3*3。

周期性这里是Slab（二维周期性，上下表面外部为半无限大真空），需要改为Bulk，并在Model - Latice中设置适当的Z方向高度：

Task选择GCMC：

然后点击GCMC右侧的>按钮，详细设置GCMC参数。点击Molecules前面的＋，可以增加一种粒子源，这里只设置一种。Molecule下拉菜单选择New Molecule，左侧窗口增加了一个，用户可以将粒子创建在该窗口中。窗口底部的Mol-1、Mol-2可以切换：

用户还需要设定：

系综（Mμ-VT化学势、体积、温度恒定，Mμ-PT化学势压强温度恒定）
MC的步数，这里设定为30000
温度，这里设置1500K
Chemical Potential，这里设置为-0.8271 Hartree（后面的单位是可以更改的），注意MC化学势和一般实验测量的分子化学势定义上有所不同

原子的MC化学势：如果粒子源是原子，那么原子的MC化学势可以用统计力学从第一原理推导出来，或者从文献资料中的热化学表中查询。单个氧原子的化学势等于同温度压强下$O_2$的MC化学势的1/2。

O原子的MC化学势计算公式：

因为$O_2$是双原子分子，所以系数为1/2
μ$_{O2}^{ref}$(T,P$_{ref}$)是$O_2$在温度T、压强P$_{ref}$下，实验测得的化学势，如果查不到化学势，可以近似使用“分子的频率、红外IR、零点能、转动能量、转动惯量、熵、焓、热容与Gibbs自由能的计算“中的吉布斯自由能作为化学势
k为波尔兹曼常数，T为温度
P是$O_2$的分压。如果只有$O_2$一种气体，那么kT*ln(P/P$_{ref}$)＝kT*ln1＝0
E$_{O2}^{diss}$是分子的解离能，也就是分子解离为原子所需的能量（正值），可以使用ADF结构优化得到的Bond Energy的绝对值

分子的MC化学势：一般分子的MC化学势计算公式，可以参考上述公式，只是不用乘以系数1/2

在movie中可以看到Si的氧化中间过程：

能量曲线平稳后，系统也达到平衡状态。