基于超算互联网的Gromacs水结冰模拟

在材料工程领域，水的相变行为（如过冷水结冰）是影响能源储存、材料稳定性及气候模型预测的核心科学问题。

本次实操，我们详细介绍如何基于超算互联网利用分子动力学软件Gromacs进行水结冰模拟全流程：通过TIP4P/ICE水模型构建冰-水两相体系，结合能量最小化、预平衡与长时相变模拟，探究过冷条件下水分子结晶的动力学机制，并结合VMD可视化与Python数据分析，为相变研究、材料设计及气候模型提供可复现的技术实践指南。

超算互联网平台提供预编译的Gromacs软件，支持一键购买和使用。本文测试所用环境与平台提供版本一致，用户可直接通过平台进行计算。

https://www.scnet.cn/ui/mall/search/global?keyword=GROMACS+v2020.6

您可以通过“我的商品”，在“应用软件”中找到或者搜索该版本商品，点击“命令行”，即可选择区域中心、命令行快速使用软件。

一、前期准备

水模型介绍

首先我们以TIP4P四点水模型和TIP3P三点水模型为例，进行水结冰的分子动力学模拟。

TIP4P是一种四点水模型，亦即除了一个O和2个H原子外还加入了一个虚原子（在氧附近），虚原子仅具有负电荷，改善了水分子周围的静电分布。

TIP4P水模型如下所示：

相比于TIP4P水，TIP3P没有虚原子，整体电荷由3个原子分摊，这里我们不多做赘述。由于Gromacs的力场中没有自带该TIP4P-ICE水模型，需要到gromacs官网获取，这里我们提供TIP4P-ICE.itp文件。

而TIP3P.itp文件opls-aa力场中自带，我们可直接引用无需额外生成。

冰结构生成

接下来，我们利用genice2生成4×5×6结构的1h六方型冰晶结构，运行命令：

genice2 --water tip4p --rep 4 5 6 1h > ice456.gro

至此，我们得到了冰的基本结构文件，再利用gmx editconf 命令将其转化为pdb结构，从ice456.gro文件的最后一行得到冰的尺寸信息，以便我们后续构建冰水混合结构。

冰水混合结构搭建

最后，利用packmol搭建冰水混合结构，在Z轴方向上方放水分子，下方放我们生成的冰结构。再利用gmx editconf 将得到的pdb文件转化为带尺寸信息的gro文件。至此我们的结构文件准备完毕。

Packmol输入文件示例：

二、模拟流程

能量最小化

第二步，我们对体系进行能量最小化，对水采用tip3p模型，对冰采用tip4p-ice模型。

输入文件：

初始构型：box1.gro 拓扑文件：em.top（含TIP3P水模型参数） MDP参数文件：em.mdp（设置能量最小化方法、收敛阈值等）

运行命令：

gmx_mpi grompp -f em.mdp -c box1.gro -p em.top -o 01em.tpr
mpirun -np 32 gmx_mpi mdrun -ntomp 1 -v -s 01em.tpr -deffnm 01em

虚原子添加：

接着，在不改变原坐标的情况下，给所有水引入虚原子，转化为四点水模型的坐标文件。即在水分子的HW2原子坐标下添加MW虚原子，生成新的坐标文件 icebox.gro。

复制tip4p-ice.itp 为tip4p-ice-sol.itp，将resname ICE修改为SOL，保证和gro文件的残基名一致。

最后，再对体系进行能量最小化，对水采用tip4p-ice-sol模型，对冰采用tip4p-ice模型。

gmx_mpi grompp -f em.mdp -c icebox.gro -p em-ice.top -o 02em.tpr
mpirun -np 32 gmx_mpi mdrun -ntomp 1 -v -s 02em.tpr -deffnm 02em

这样我们便得到了基于TIP4P的冰/水坐标文件。

预平衡与产生相模拟

第二步，进行NPT、NVT预平衡与产生相模拟。NVT和NPT预平衡的mdp文件以及产生相的mdp文件我会上传到网站，运行以下代码：

NPT预平衡代码：gmx_mpi grompp -f npt-pre.mdp -c 02em.gro -p npt-pre.top -o 03npt.tpr
mpirun -np 32 gmx_mpi mdrun -ntomp 1 -v -s 03npt.tpr -deffnm 03npt
NVT预平衡代码：gmx_mpi grompp -f nvt-pre.mdp -c 03npt.gro -p nvt-pre.top -o 04nvt.tpr
mpirun -np 32 gmx_mpi mdrun -ntomp 1 -v -s 04nvt.tpr -deffnm 04nvt
产生相代码：gmx_mpi grompp -f npt.mdp -c 04nvt.gro -p npt.top -o 05npt.tpr
mpirun -np 32 gmx_mpi mdrun -ntomp 1 -v -s 05npt.tpr -deffnm 05npt

通过上述流程，我们完成了从能量最小化到产生相模拟的流程，接下来对结果进行可视化分析。

三、后期分析

轨迹处理

首先，打包下载相关文件，消除周期性边界条件对轨迹的影响，运行命令：

gmx trjconv -f 05npt.xtc -s 05npt.tpr -pbc mol -n new.ndx -o 05npt-notraj.xtc

命令执行时提示选择输出组，输入 0（选择所有原子）以保留完整体系信息。

可视化分析

利用vmd 载入构象和轨迹，可以很清楚的看到水转化为冰的过程。

以上，我们就完成了使用Gromacs水结冰模拟的全流程，并结合VMD可视化与Python数据分析，为研究人员提供了可复现的理论框架与技术路径。