谈谈做分子动力学模拟什么时候适合用NPT和NVT系综

sobereva

谈谈做分子动力学模拟什么时候适合用NPT和NVT系综

On the appropriate use of NPT and NVT ensembles in molecular dynamics simulations

文/Sobereva@北京科音

First release: 2026-Jan-11   Last update: 2026-Jan-14

之前我天天在网上答疑时看过太多初学者做分子动力学模拟时本该用NPT系综的时候却误用或乱用NVT系综，导致自取其辱，或者白浪费很多时间。甚至今日在同一天内，还看到计算化学公社论坛里两个新帖http://bbs.keinsci.com/thread-57828-1-1.html和http://bbs.keinsci.com/thread-57823-1-1.html都是乱用NVT而导致了迷惑的结果，我遂感到很有必要专门写个小文批判乱用NVT的行为、澄清什么时候才真正适合用NVT。实际上本文内容就是北京科音分子动力学与GROMACS培训班（http://www.keinsci.com/KGMX）里讲“分子动力学的基本原理与算法”中的一部分，并稍作扩展而已，但凡参加过此培训完整系统性学习过分子动力学的学员都不会在这种常识性问题上犯错误，压浴的算法在培训里也有非常全面细致的介绍。

首先明确一点：做常见体系的动力学模拟，应当默认用NPT系综，即压浴（barostat）和热浴（thermostat）同时使用，或者说同时考虑控压（压力耦合）和控温（温度耦合）。这显而易见，毕竟大多数情况模拟的体系在现实中都是在恒定温度和压力条件下的，比如烧杯中的溶液、细胞中的蛋白质。若关掉控压，即不设置压浴，就是NVT系综了。

先说几个常识。体系盒子尺寸的改变会令压力改变，盒子越大内压越小，盒子越小内压越大。控压有不同算法，共同原理都是在动力学过程中不断调节盒子尺寸并同时按比例缩放原子坐标，以使得当前体系压力的平均值与自设的参考压力相符，但允许在模拟过程中有一定程度的波动。每次控压时要求三个方向的盒子尺寸按相同比例变化称为各向同性控压，XY和Z方向独立控压称为半各向同性控压，三个方向都独立控压称为各向异性控压。压浴的效果严重依赖于人为指定的可压缩系数（注：特别常用的Berendsen和C-rescale压浴需要设，Parrinello-Rahman、MTK压浴倒是不需要设），它体现温度不变的情况下增加单位压力令体积的相对减小量。各向异性控压时可以在盒子不同方向上分别设这个参数，如果在某个方向上可压缩系数设为0，则这个方向上即便开启了控压盒子尺寸也不会改变。本文说的控压默认是指设的可压缩系数>0的情况。

有的情况用压浴是绝对必须的，例如：
• 用packmol（用法介绍见http://sobereva.com/473）、GROMACS自带的insert-molecules等程序构建的体系初始模型对应的盒子往往明显偏大，因而对应的体系初始密度往往明显偏小。如果你在模拟时不考虑控压，显然在分子之间由于相互作用聚集到一起后盒子里就会出现真空区。只有用控压时才可能通过自发调节盒子尺寸令体系密度与实际一致而消掉真空区。前面提到的http://bbs.keinsci.com/thread-57828-1-1.html和http://bbs.keinsci.com/thread-57823-1-1.html这两个帖子里的截图一看就知道都是该用NPT却误用了NVT导致盒子里出现了明显真空区。顺带一提，用VMD观看周期性体系的时候强烈建议用pbc box命令把盒子显示出来，免得被蒙在鼓里，而这俩帖子都没显示出盒子。
• 研究相变过程。相变会伴随着密度的变化，比如水变成冰体积会自发显著膨胀。显然如果你模拟水结冰，而让盒子尺寸始终不变的话，结冰过程都无法如实模拟。
• 研究变温过程。比如你对一个材料热膨胀系数为正的温度区间内研究升温过程，考虑控压的话能如实描述热膨胀效应，若不允许控压则相当于内压越来越大，跟高压锅似的。
• 有的研究专门模拟很高压力下的情况，显然也必须用控压了，除非你能根据已有的拟合关系事先知道这个压力下的盒子尺寸理应是多大，并将之作为初始盒子尺寸。

虽然用了压浴后会增加耗时，但比NVT耗时增加程度甚微，因此若没有什么不能用控压的明确理由时都应当开着压浴。不适合用控压的情况比较有限，应当用NVT的常见情况有以下这些：
• 需要维持特定真空区的体系，典型的是气液界面。在垂直于界面方向上若用控压会导致在此方向上盒子不断收缩、真空区逐渐消失，经过足够长时间的模拟后体系就变成了纯液体盒子了。气液界面体系在平行于界面的方向上也不应该用控压，否则由于表面张力原因界面会越来越小。显然模拟气液界面要用NVT。如果要模拟液液界面，虽然平行于界面方向上不应该控压（或者设可压缩系数为0也行），但垂直于界面方向上应该控压，从而令液体密度能自发调节到合理。
• 周期边界条件下算孤立体系。目前版本的GROMACS已经不支持非周期边界条件了，因此模拟单分子或者分子团簇时，只能用一个足够大的盒子当三维周期性做NVT模拟，并且为了避免体系与其相邻镜像存在相互作用而需要令盒子足够大。此情况若用控压显然也会令盒子逐渐缩小，导致体系与其镜像最后虚假地作用上了。
• 想精确维持特定的密度。例如使用已知的实验密度进行模拟、考察特定性质随密度的变化（跑很多次模拟，每次设置不同的盒子体积），显然需要NVT。
• 模拟一些特殊体系。比如超临界气体的可压缩系数很大，用控压时盒子尺寸往往会剧烈变化，甚至导致崩溃，因而通常基于其实验已知密度来定义特定尺寸的盒子进行NVT模拟。
• 计算某些性质的时候原理上用NVT比NPT更合适，对于这种情况可以先做NPT模拟使得盒子尺寸达到平衡，然后再切换到NVT续跑。例如计算MSD曲线并得到扩散系数，用控压的话由于在模拟过程中会不断对原子坐标进行缩放，多多少少会对结果带来不利影响（不过对于液体体系模拟来说，盒子尺寸通常波动很小，因此也不能说NPT下算的结果就不能用）。再比用周期扰动法计算粘度，原理上就应当在NVT下进行，此方法的原文里也是这么做。
• 用控压时发现模拟崩溃，并怀疑是控压所致，此时可以暂时尝试不用控压看看动力学情况。这点我之前在http://sobereva.com/soft/Sobtop#FAQ8里介绍怎么排查动力学崩溃原因的时候也说过。如果改用NVT跑了一阵子都没问题，可以尝试再把控压打开续跑，看看NPT能否顺利跑下去；如果能的话，说明之前做NVT的时候令体系得到了足够的弛豫、状态被稳定化了。如果一上来就用NPT模拟虽然没崩溃，但是盒子变化异常，比如尺寸波动剧烈，也可以先跑一段时间NVT试图稳定化一下。
• 在GROMACS里利用freeze冻结原子时，同时开启压浴容易造成崩溃，因此要用NVT，或者用NPT结合位置限制势实现类似目的。

还值得一提的是做盒子以特定方式变化的模拟的情况。GROMACS里可以在mdp里用deform设置要求盒子的特定分量按照特定速度发生变化，比如你设了Z方向尺寸的变化速度，那么Z方向就不能用控压来自发调节盒子。

网上有一些教程，以及有不少文章，做NPT之前都先做NVT，我强烈不建议盲目效仿！！！这事我在网上答疑时已经说过无数次。NPT之前做NVT有>99%的可能都是多此一举的，白浪费你的时间。虽然如前所述，NPT之前做NVT有可能能令体系先稳一稳避免直接NPT模拟时异常或崩溃，但实际中必须先做NVT的情况甚少。只要你的拓扑文件合理、体系初始结构别太恶心（比如存在严重不合理接触、初始压力特别离谱）、模拟设置没有什么硬伤、跑动力学之前已做了充分的能量极小化，并且你最终就是要在NPT下模拟，那么上来直接用NPT模拟就行了。倘若万一真的发现异常，可以到时候再考虑做NVT预平衡一下。即曰此时NVT是补救时才需要用的，而不要当做默认要做的过程！

之前在我的思想家公社QQ群里谈论以上问题时，有人提到Molecules, 25, 5120 (2020)一文以试图论证NPT之前做NVT的意义。这篇文章的意思是对此文研究的糖脂这类体系，当初始结构的压力明显不合理的时候，平衡相一开始先用NVT再NPT能避免直接用NPT时由于模拟初期盒子显著震荡、相应地缩放原子位置导致水被巧合地弄进疏水区。但实际上此发现明显不足以构成预平衡先要用NVT的“普遍的必要性”，这只是属于小众体系的罕见情况。只能说对于没有人工监控的靠脚本实现的批量模拟，并且你也不在乎NVT多跑几百ps花的计算量，倒是可以加入先做NVT的过程以求绝对保险。顺带一提，也有很多动力学研究的本身就是非平衡过程，比如两相间的互溶，这种情况就不可能先加入NVT预平衡过程了。

对初学者们，本文最后留下一句话进行强调：若无明确、说得通的理由，不要用NVT！！！

Uus/pMeC6H4-/キ

可以也讨论一下（除了对正确能量交换/动力学行为要求很严格的NEMD非平衡模拟外）连热浴控温都不能开，而必须用NVE系综（哪怕实际计算中数值误差再怎么影响能量守恒性）的情形吗？

论坛上有些用M$-forcite算扩散系数的古老讨论帖会提到一种先NPT或NVT系综平衡再用NVE统计的流程，比如thread-9705-1-1.html的7楼用“NVT/NPT里的控温/控压方法可能会对轨迹产生扰动从而影响扩散行为”的观点论证其合理性，但是在thread-2750-1-1.html里的回应“从原理上说，并不要求必须NVE，用NVT就完全可以”没有细说原理。即使不看特定软件的问题，也有thread-15360-1-1.html的9楼认为只有Langevin dynamics能得到正确的ensemble，thread-9699-1-1.html的5楼解释不同时期发展的控温算法对水分子动力学的影响时又认为有些问题不大。

一些讨论NVT系综的文献要具体区分热浴控温算法是velocity-randomized还是velocity-rescaling的，例如
J. Chem. Theory Comput. 2013, 9, 7, 2887–2899, DOI: 10.1021/ct400109a
J. Mol. Liquid, 2022, 365, 120116, DOI: 10.1016/j.molliq.2022.120116
J. Phys. Chem. B 2022, 126, 48, 10172–10184, DOI: 10.1021/acs.jpcb.2c06035

另外，Phys. Chem. Chem. Phys., 2020,22, 10462-10479, DOI: 10.1039/C9CP05610F在模拟不同密度的亚临界/超临界水时用的了比较复杂的采样方式，在用NVT系综跑SPC水的经典力场MD后取样用NVT系综跑RPBE-D3水的AIMD，再取样用NVE系综跑AIMD并统计平均……

We conclude that all properties reported in this paper were obtained by averaging over 40 statistically independent NVE AIMD simulations at each state point, which were drawn from corresponding NVT AIMD simulations at each state point, thus providing proper and rigorous access to time-correlation functions in the canonical ensemble (rather than using thermostatted AIMD trajectories to compute such time-correlation functions).

sobereva

Uus/pMeC6H4-/キ 发表于 2026-1-11 11:45
可以也讨论一下（除了对正确能量交换/动力学行为要求很严格的NEMD非平衡模拟外）连热浴控温都不能开，而必 ...

那些严格要求能量守恒的情况不在当前文章涉及范围内。

算扩散系数结合主流的热浴用NVT没实际问题（如果对于鸡毛蒜皮程度的原理正确性都过于较真，大量研究都没法做了，或者做得过于辛苦）。切换成NVE后，由于数值计算误差，温度在MD过程中都往往没法很好守恒，这反倒可能是更大问题。

Uus/pMeC6H4-/キ

sobereva 发表于 2026-1-12 00:49
那些严格要求能量守恒的情况不在当前文章涉及范围内。

算扩散系数结合主流的热浴用NVT没实际问题（如 ...

谢谢。我还有几个关于压浴的问题：

1. 文中说“压浴的效果严重依赖于人为指定的可压缩系数”，这与模拟使用的理论方法是经典分子力场还是量子化学/第一性原理方法有关吗？我看像GROMACS有compressibility的设置而像CP2K就没找到，是不是意味着在NPT系综跑AIMD到平衡后要自行从轨迹的盒子尺寸或应力张量涨落统计等温可压缩系数呢，如果是的话有什么算法统计？
2. 控压时调节盒子尺寸、缩放原子坐标的操作是不是意味着无法定义绝对空间坐标固定的参考点？可否采用一种类似2楼图的策略，即从NPT系综跑平衡的轨迹中取样（间隔大于self-correlation时间使得样帧是statistically independent的），从各样帧的瞬时原子坐标、速度和盒子尺寸出发用NVT系综平行跑多个任务，且只在后面NVT系综的模拟中定义绝对空间坐标固定的参考点？CP2K的&MOTION/&MD/ENSEMBLE有个NPT_IA选项可以定义"NPT_I ensemble with frozen atoms in absolute coordinate"，这在GROMACS里有对应设置吗？
3. CP2K的&MOTION/&MD/ENSEMBLE关键词还有NPE_F和NPE_I的选项，这种只开压浴控压不开热浴控温的NPE模拟有什么实际意义呢？

sobereva

Uus/pMeC6H4-/キ 发表于 2026-1-12 13:10
谢谢。我还有几个关于压浴的问题：

1. 文中说“压浴的效果严重依赖于人为指定的可压缩系数”，这与模 ...

1 我这里说的是常用的Berendsen压浴情况。PR压浴不依赖于可压缩系数，CP2K只能用MTK压浴，和PR很相近，故不需要定义可压缩系数。

2 GROMACS中，在控压时被冻结的组的坐标保持不变（使用PR压浴时被冻结的组的坐标会随着盒子尺寸的缩放按比例调节，此bug从2021.5版开始被修正）

3 应该极少能派上用场，开发者提供它们的目的可能在于功能的完整性

喵星大佬

Uus/pMeC6H4-/キ 发表于 2026-1-12 13:10
谢谢。我还有几个关于压浴的问题：

1. 文中说“压浴的效果严重依赖于人为指定的可压缩系数”，这与模 ...

第三个，燃烧爆炸的情况

Uus/pMeC6H4-/キ

喵星大佬 发表于 2026-1-13 04:45
第三个，燃烧爆炸的情况

这应该也算是非平衡模拟了吧，毕竟反应放热是关键。我感觉论坛很多人研究的剪切/打磨/抛光/润滑过程也是类似，既需要控压保证接触紧实，又因为摩擦生热而不该维持特定温度，更何况恒速/恒外力/恒加速度/恒压等载荷驱动控制的原子运动也不太兼容热浴调节原子速度。

编辑：我想问的可能更接近于，NPE（或者再加上等焓条件的NPH）模拟平衡以后，可以获得什么样的统计平均量，是NPT或NVT模拟得不到的？