【摘要】 本教程介绍如何使用GROMACS对分子动力学轨迹进行一维数密度分布分析，以NaCl水溶液气液界面为例，讲解gmx density命令、参数设置、Na⁺和Cl⁻分析组选择及结果绘图方法。

科研人必备干货｜全套科研知识库

常用科研资料汇总，点击领取即可获取。

 

本教程默认用户已经完成能量最小化和NVT模拟，并获得prod.xtc与prod.tpr文件；模拟体系的气液界面法向沿Z轴，且拓扑或索引中已包含Na⁺和Cl⁻分析组。

 

今天我们将介绍分子动力学模拟中关于数密度的分析方法。

使用软件：GROMACS

 

一、原理介绍

GROMACS 中一维数密度剖面的核心思路是：沿选定方向（X、Y 或 Z）把模拟盒子划分成若干薄层，统计每一帧中落入各薄层的目标粒子数，再除以薄层体积并对轨迹帧求平均。

 

二．案例展示

下面以氯化钠水溶液气液界面为例，统计该体系中钠离子和氯离子的数密度分布情况。初始模型如图 1 所示，体系包含 20 个钠离子、20 个氯离子和 844 个水分子。氯化钠水溶液位于模拟盒子中部，界面法向两端留有真空区，用于形成两个气–液界面。

图1. 体系初始模型。

做完能量最小化后进行 5 ns 的 NVT 模拟。模拟结果如图 2 所示：

图2. 模拟 5 ns 后体系结构。

模拟结束后输入以下分析命令，统计 Na⁺ 和 Cl⁻ 在 5 ns 模拟期间沿着 Z 轴方向的数密度分布：

gmx density -f prod.xtc -s prod.tpr -dens number -ng 2 -sl 300 -symm

依次选择 system, NA, CL，得到数密度分布图 density.xvg。

图3. 命令行界面。

输入各命令所代表的含义分别如下：

  • gmx density：密度统计分析

  • prod.xtc：压缩的轨迹文件

  • prod.tpr：mdrun 的输入文件

  • -dens：设定统计什么密度，默认是质量密度，number 则是数密度，但需注意统计的是原子数密度，而不是分子数密度

  • -ng：设定考察几个组，默认是 1

  • -sl：设定把盒子划分为多少格，默认是 50

  • -symm：参照某个组的中心来把曲线对称化

之后绘制图像：

xmgrace -nxy density.xvg

得到如下图 4：

图4. 离子数密度分布图。

可以看到，随着模拟的进行，钠离子和氯离子会逐渐往气液界面富集，界面处的离子数密度明显高于体相中离子的数密度。同理，也可以得到水分子的数密度，由于水分子数目更多，得到的数密度分布曲线自然也更光滑。

图5. 水分子数密度分布图。

 

常见问题

1.GROMACS如何计算一维数密度分布？​

GROMACS会沿指定方向将模拟盒子划分为多个薄层，统计每一帧中目标粒子落入各薄层的数量，并根据薄层体积和轨迹帧进行平均，得到数密度剖面。

2.如何用GROMACS分析Na⁺和Cl⁻的数密度？​

在完成模拟并获得轨迹文件后，可以使用：

gmx density -f prod.xtc -s prod.tpr -dens number -ng 2 -sl 300 -symm

随后按提示选择Na⁺和Cl⁻对应的分析组，并使用Xmgrace绘制输出的.xvg文件。

3.-dens number统计的是分子数密度吗？

不一定。-dens number统计的是所选分析组中粒子的数量密度。对于Na⁺和Cl⁻这类单原子离子，可近似理解为离子数密度；对于水分子，如果选择的是氧原子组，则得到的是氧原子数密度，而不是完整水分子数密度。

4.为什么要使用-sl 300？

-sl 300表示将分析方向划分为300个薄层。分层越细，空间分辨率通常越高，但每层包含的粒子数可能减少，曲线噪声也可能增大。应根据盒子尺寸、粒子数量和轨迹长度选择合适的分层数。

5.什么时候适合使用-symm？

当体系沿界面法向具有合理对称性，且希望将对称位置的统计结果进行整合时，可以考虑使用-symm。对于不对称界面或存在外场、浓度梯度的体系，使用前应谨慎检查其是否会掩盖真实差异。