关键词：MOF材料；IRMOF-1；GCMC；巨正则系综蒙特卡洛；ReaxFF

吸附等温线表示“主体材料”在给定温度和压力下，摄入给定数量“客体分子”的性质。吸附等温线是负载与压强的函数，负载通常表示为每千克主体材料吸附多少摩尔客体分子，单位为 mol/kg）。在本教程中，我们将以通用性较好的 MOF 材料 IRMOF-1 对 CO$_2$ 的负载为例，演示如何使用巨正则系综蒙特卡洛模拟生成吸附等温线。可以使用各种力场、机器学习势，使用过程基本上没有大的区别，不过本例中以 ReaxFF/CHOZn.ff 为例进行演示。CHOZn.ff 力场在 AMS2023 中尚未加入，用户可自行下载解压使用（AMS2024 版中将加入该力场）。

在 ReaxFF 中创建并优化 CO$_2$ 如下所示：

File → Save，保存名为例如 CO2_ref_ReaxFF，并运行。运行结束后，SCM → logfile 在尾部可以看到总能量E$_{CO2}$ = -0.6434 Hartree = -1689.155kJ/mol：

点开链接，复制所有数据，在 AMSinput Ctrl v粘贴（软件会提示“Use new lattice vector?”，点击Yes），即得到 IRMOF-1 结构。

优化结构：

如果压强不为0，则优化时，应在 Details → Geometry Optimization 窗口，勾选 Optimize Lattice，且在 Pressure 栏输入压强。

保存并运行作业。例如作业名为：IRMOF-1_ref_ReaxFF。

GCMC 模拟时，需要提供在指定的这个温度、压强下，体系的化学势。化学势 μ(T,P) 的通用定义如下：

\[\begin{split}\mu(T,P) = E_{CO2}+\Delta\mu(T,P_0)+RT\ln{\frac{P}{P_0}}\\ \text{其中 }\Delta\mu(T,P_0) = H(T,P_0)-H(0,P_0)-T\times \left[S(T,P_0)-S(0,P_0)\right]\end{split}\]

CO$_2$ 在不同温度下的 H 和 S 的值可以在 NIST 网站上直接查到。

298.15 K 时，不同压力下的 H 和 S 查询结果如下：

可以看到计算 Δμ(T,P$_0$) 所需的几个数据，T = 298.15K 时：

H(T,P$_0$) - H(0,P$_0$) = H(298.15K,P$_0$) - H(0K,P$_0$) 即表中 H-H$^。$(T$_r$) 这一列的两个数据，因此二者之差为：0.0 - (-9.364) kJ·mol$^{-1}$

S(T,P$_0$) - S(0,P$_0$) = H(298.15K,P$_0$) - H(0K,P$_0$) = 0.213795 - 0 kJ/(K·mol)

注意表中的数据是 213.795 J/(K·mol) 换算单位为 kJ/(K·mol)，因此为 0.213795 kJ/(K·mol)

带入上文中 Δμ(T,P$_0$) 与 μ(T,P) 的计算公式中，压强 P = 30 bar 时：

\[\begin{split}\Delta\mu(298.15\text{ K},1.01\text{ bar}) = \left(0.0 - \left(-9.364\right)\right) - 298.15\times\left(0.213795 - 0.0\right) = -54.379\text{ kJ/mol} \\ \mu(298.15\text{ K},30\text{ bar}) = -1689.155 - 54.379 + 0.008314 \times 298.15 \times \ln{\frac{30.0}{1.01325}} = -1735.136\text{ kJ/mol}\end{split}\]

其中 E$_{CO2}$ 即前面优化 CO$_2$ 得到的能量 -1689.155kJ/mol，常数 R = 0.008314。类似地计算298.15K，压强为 5、10、15、20、25、30 bar 时，化学势为：

注：

• 上表中查阅到的 298K 以上的焓值和熵值，也可以通过 Shomate 方程计算得到，具体参考 NIST 网站（不过需要注意该网站无法算出 298K 以下温度的数据，而 0K 的数据，在我们上述计算中是需要用到的，因此仅供参考）。
• 根据焓和熵的来源，以及您使用的精度和转换因子，您独立计算的数字可能略有不同。