GROMACS分子模拟基础实验教程(7-8)-nvt与npt平衡

文摘科学 2023-05-26 09:59 江苏

我们接着之前的实验文档：

GROMACS分子模拟基础实验教程(1)-运行环境准备

GROMACS分子模拟基础实验教程(2)-生成蛋白结构拓扑

GROMACS分子模拟基础实验教程(3)-检查和理解拓扑文件

GROMACS分子模拟基础实验教程(4)-定义模拟体系的容器尺寸并加水

GROMACS分子模拟基础实验教程(5)-为模拟体系添加离子

GROMACS分子模拟基础实验教程(6)-体系能量最小化

今天是GROMACS分子模拟基础实验教程(7-8)-nvt（温度）与npt（压力、密度）平衡

第7步-NVT温度平衡

7.1、NVT平衡是做什么

能量最小化（EM）的目的是促使体系初始结构在几何构型和溶剂分子取向等方面都合理，为了开始真正的动力学模拟, 我们还需要对蛋白质周围的溶剂和离子进行平衡。如果未经平衡就直接尝试进行非限制的动力学模拟，体系可能会崩溃，原因在于目前我们只是优化了溶剂分子本身, 而没有考虑溶质。我们需要将体系置于设定的模拟温度下, 以确定蛋白质溶质的合理结构取向。达到目标温度后, 再对体系施加压力，达到合适的密度。

之前我们用pdb2gmx生成了posre.itp文件，此文件的目的在于对蛋白质中的非氢原子类的重原子施加位置限制力(position restraining)。这些原子一般不会移动, 除非增加非常大的能量。位置限制的作用在于平衡蛋白质周围的溶剂分子, 而不引起蛋白质结构的变化。

平衡往往分两个阶段进行：

第一个阶段在NVT系统包括粒子数, 体积和温度都恒定的情况下进行。这个系统也被称为等温、等容系综。这个过程的需要的时间与体系的构成有关, 但在NVT系统中, 体系的温度应达到预期值并基本保持不变，如果温度仍然没有稳定, 那就需要更多的时间。一般50-100 ps就足够了, 因此在本实验中我们进行100ps的NVT平衡。

7.2、nvt mdp参数

运行之前需要准备的.mdp文件如下：

获取链接 http://www.mdtutorials.com/gmx/lysozyme/Files/nvt.mdp

title                   = OPLS Lysozyme NVT equilibration define                  = -DPOSRES  ; position restrain the protein; Run parametersintegrator              = md        ; leap-frog integratornsteps                  = 50000     ; 2 * 50000 = 100 psdt                      = 0.002     ; 2 fs; Output controlnstxout                 = 500       ; save coordinates every 1.0 psnstvout                 = 500       ; save velocities every 1.0 psnstenergy               = 500       ; save energies every 1.0 psnstlog                  = 500       ; update log file every 1.0 ps; Bond parameterscontinuation            = no        ; first dynamics runconstraint_algorithm    = lincs     ; holonomic constraints constraints             = h-bonds   ; bonds involving H are constrainedlincs_iter              = 1         ; accuracy of LINCSlincs_order             = 4         ; also related to accuracy; Nonbonded settings cutoff-scheme           = Verlet    ; Buffered neighbor searchingns_type                 = grid      ; search neighboring grid cellsnstlist                 = 10        ; 20 fs, largely irrelevant with Verletrcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)DispCorr                = EnerPres  ; account for cut-off vdW scheme; Electrostaticscoulombtype             = PME       ; Particle Mesh Ewald for long-range electrostaticspme_order               = 4         ; cubic interpolationfourierspacing          = 0.16      ; grid spacing for FFT; Temperature coupling is ontcoupl                  = V-rescale             ; modified Berendsen thermostattc-grps                 = Protein Non-Protein   ; two coupling groups - more accuratetau_t                   = 0.1     0.1           ; time constant, in psref_t                   = 300     300           ; reference temperature, one for each group, in K; Pressure coupling is offpcoupl                  = no        ; no pressure coupling in NVT; Periodic boundary conditionspbc                     = xyz       ; 3-D PBC; Velocity generationgen_vel                 = yes       ; assign velocities from Maxwell distributiongen_temp                = 300       ; temperature for Maxwell distributiongen_seed                = -1        ; generate a random seed

7.3、实操15-执行NVT平衡操作

Linux下执行：

第1步：

gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr

第2步：

gmx mdrun -deffnm nvt

（根据机器的不同, 运行需要一段时间，我这里12核CPU大概10分钟运行完成）

第1步执行截图：

第1步得到 nvt.tpr文件

第2步执行截图

第2步是执行NVT平衡。可以查看输出提示需要进行100ps的运行

7.4、查看运行的详细进度

现在通过linux的top命令可以看到CPU使用率比较高

使用tail -f nvt.log可以实时查看输出结果：

我这个服务器大概10分钟完成了（12核CPU）

运行结束后，会给出相关的汇总信息：

7.5、实操16-NVT平衡结果分析

Linux执行：

gmx energy -f nvt.edr -o temperature.xvg

在提示符下键入“16 0”，选择温度代号16，选择0退出

继续使用xmgrace工具查看变化图：

xmgrace temperature.xvg

可以看出，系统的温度迅速达到目标值300K左右，并在平衡的剩余时间内保持稳定（按这个图的结果，总体波动不大，其实较短的平衡时间(大约50 ps)应该就足够了。

第8步-NPT压力平衡

8.1、NPT平衡是做什么

上一步NVT平衡建立了温度系统，这一步需要建立压力和密度系统（稳定系统的压力，以及保证合适的密度）。压力的平衡是在NPT系统下进行的，其粒子数, 压力和温度都保持不变，这个系统也被称为等温、等压系统，比较接近实验环境。

8.2、npt mdp参数

100-ps NPT equilibration 配置参数链接：http://www.mdtutorials.com/gmx/lysozyme/Files/npt.mdp

这个参数文件与NVT平衡时所用的参数文件类似。

注意添加的压力耦合部分, 其中使用了Parrinello-Rahman控压器. 其他几项改动如下:

continuation = yes: 我们将从NVT平衡阶段开始继续进行模拟

gen_vel =no: 从轨迹中读取速度（Velocities are read from the trajectory ）

$ cat npt.mdp title                   = OPLS Lysozyme NPT equilibration define                  = -DPOSRES  ; position restrain the protein; Run parametersintegrator              = md        ; leap-frog integratornsteps                  = 50000     ; 2 * 50000 = 100 psdt                      = 0.002     ; 2 fs; Output controlnstxout                 = 500       ; save coordinates every 1.0 psnstvout                 = 500       ; save velocities every 1.0 psnstenergy               = 500       ; save energies every 1.0 psnstlog                  = 500       ; update log file every 1.0 ps; Bond parameterscontinuation            = yes       ; Restarting after NVT constraint_algorithm    = lincs     ; holonomic constraints constraints             = h-bonds   ; bonds involving H are constrainedlincs_iter              = 1         ; accuracy of LINCSlincs_order             = 4         ; also related to accuracy; Nonbonded settings cutoff-scheme           = Verlet    ; Buffered neighbor searchingns_type                 = grid      ; search neighboring grid cellsnstlist                 = 10        ; 20 fs, largely irrelevant with Verlet schemercoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)DispCorr                = EnerPres  ; account for cut-off vdW scheme; Electrostaticscoulombtype             = PME       ; Particle Mesh Ewald for long-range electrostaticspme_order               = 4         ; cubic interpolationfourierspacing          = 0.16      ; grid spacing for FFT; Temperature coupling is ontcoupl                  = V-rescale             ; modified Berendsen thermostattc-grps                 = Protein Non-Protein   ; two coupling groups - more accuratetau_t                   = 0.1     0.1           ; time constant, in psref_t                   = 300     300           ; reference temperature, one for each group, in K; Pressure coupling is onpcoupl                  = Parrinello-Rahman     ; Pressure coupling on in NPTpcoupltype              = isotropic             ; uniform scaling of box vectorstau_p                   = 2.0                   ; time constant, in psref_p                   = 1.0                   ; reference pressure, in barcompressibility         = 4.5e-5                ; isothermal compressibility of water, bar^-1refcoord_scaling        = com; Periodic boundary conditionspbc                     = xyz       ; 3-D PBC; Velocity generationgen_vel                 = no        ; Velocity generation is off

8.3、实操17-执行NPT平衡操作

NPT平衡也是使用gmx的grompp和mdrun模块

1、-t选项，指定NPT平衡过程中的产生的检查点文件名（这个文件包含了继续模拟所需要的所有状态变量）

2、-c选项，指定NPT模拟最终输出的坐标文件

第1步：

gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr

第2步：

gmx mdrun -deffnm npt

提示时输入【16 0】来选择体系压力并退出

执行gmx mdrun后给出的提示：

Reading file npt.tpr, VERSION 2020.7 (single precision)Changing nstlist from 10 to 50, rlist from 1 to 1.115
Using 1 MPI thread
Non-default thread affinity set, disabling internal thread affinity
Using 12 OpenMP threads 
starting mdrun 'LYSOZYME in water'50000 steps,    100.0 ps.

使用了1个并行线程，使用12核CPU线程进行计算，需要执行5万步，完成100皮秒的模拟。

8.4、查看运行的详细进度

Linux下可以使用top命令，再按1，查看12核CPU使用率，基本跑满

可以通过日志tail -100f npt.log 实时查看进行的进度：

我的测试服务器Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz 12核，大概10分钟跑完

tail -100f npt.log

最后显示的日志：

生成的文件：

8.5实操18-压力分析

Linux执行：

gmx energy -f npt.edr -o pressure.xvg

根据提示输入18 0 获得压力的数据结果

查看压力变化曲线图：

xmgrace pressure.xvg

在100 ps的平衡过程中压力值涨落范围较大, 这其实也是意料中的. 图中的红线为数据的移动平均值. 在整个平衡过程中, 压力的平均值为7.5 ± 160.5 bar（参考压力值设置的是1 bar）这个结果是否可靠？从这个160.5 bar的大均方根波动可以看出，模拟过程中，压力指标变化范围是比较大的, 从统计学上讲，并不能区分 (7.5 ± 160.5 bar) 与目标参考值 (1 bar)之间的差异。

8.6、实操19-密度分析

执行：

gmx energy -f npt.edr -o density.xvg

根据提示输入 24 0 获得结果

xmgrace density.xvg

跟压力一样, 红线是密度的移动平均值. 100ps过程中密度的平均值为1019±3kg每立方米, 比较接近实验值1000 kg每立方米，SPC/E水模型的密度预期值1008 kg每立方米。SPC/E水模型的参数给出的密度值接近水的实验值。在整个过程中密度值都很稳定，随着时间推移，意味着体系的压力和密度都达到了较好的平衡

8.7 xmgrace显示补充

中间这条红色的线条如何画出：

在xmgrace中, 点击菜单Data -> Transformations -> Running averages, 在弹出的对话框中设定Length of average, Accept。

更多模拟前的基础准备文章可以参考如下链接：

基础准备1-分子动力学模拟理论概述

基础准备2-分子动力学模拟的主要步骤

基础准备3-分子模拟基础-获取Gromacs软件

基础准备4-Ubuntu16版Linux下编译安装Gromacs2020

基础准备5-Linux7.6编译安装Gromacs2018版并行版

今天的分享就到这里了

欢迎点赞、在看、分享转发朋友圈 ^.^ 公众号还可以“设为星标”喔~

联系方式：微信号 coding_gene QQ 391399793 QQ群 338633443

‍

http://mp.weixin.qq.com/s?__biz=MzI4NzQ0NDg0NA==&mid=2247485212&idx=1&sn=8c3e5e4e9d58c1f118ec34f65b222f49

生信实战

十年前你还不懂生信，看那一行行代码和报错，只是朦胧的心痛...