Gromacs安装、蛋白配体复合物模拟操作流程记录

彼此7788

于 2025-02-21 22:09:24 发布

阅读量1.3k

点赞数 29

文章标签：数据库

本文链接：https://blog.youkuaiyun.com/Lxqjy/article/details/145183893

版权

分子动力学模拟概述：分子动力学模拟大致包括准备结构文件、预处理、能量最小化、平衡阶段和生产模拟等步骤，步骤要点如下：

数据准备：下载蛋白的pdb结构文件，分子的结构文件，准备相应的拓扑文件等进行后续的模拟。
预处理阶段：将蛋白和配体组合在一起，生成复合物的结构文件，然后需要确定模拟的体系，比如水模型（SPC、TIP3P、TIP4P等）和添加离子以中和电荷
能量最小化：为了消除结构中的不合理接触，比如原子间距离过近导致的高能量。通常使用最速下降法或共轭梯度法来优化结构，使得体系的总能量达到一个局部最小值，避免模拟过程中因为结构不合理而崩溃。
平衡阶段：分为NVT和NPT平衡。NVT平衡是为了让体系在恒定粒子数、体积和温度下达到热平衡，通常使用Berendsen或V-rescale温度耦合方法。而NPT平衡则是在恒定粒子数、压力和温度下进行，让体系的密度趋于稳定，常用的压力耦合方法有Berendsen和Parrinello-Rahman。这一步可能需要运行一段时间，确保温度和压力稳定。
生产模拟：长时间的数据采集阶段，生成轨迹文件用于后续分析。分析部分包括RMSD、RMSF、氢键分析、相互作用能计算等，可能需要使用GROMACS自带的工具或者第三方软件如VMD、PyMOL、MDTraj等。

基本信息

服务器内核

uname -a

Linux amax 5.4.0-59-generic #65-Ubuntu SMP Thu Dec 10 12:01:51 UTC 2020 x86_64 x86_64 x86_64 GNU/Linux

cat /etc/issue

Ubuntu 20.04.1 LTS \n \l

安装

根据网上的传统教程先安装cmake（成功），FFTW（失败，原因 ./configure 一直报错），

configure: error: working directory cannot be determined

而且没有其他的报错原因，所有能排查的错误原因我都一一确认无误，所以这种安装方法对我来说不可行

最后成功的过程记录如下：

1. 创建一个gmx环境 conda create -n gmx

2. 激活环境 conda activate gmx

3. 安装依赖库conda-forge conda config --add channels conda-forge

4. 设置安装通道 conda config --set channel_priority strict

5. 安装相关的包 conda install pocl oclgrind intel-compute-runtime beignet ocl-icd-system

6. 安装network conda install networkx=2.3

7. 安装gromacs（非gpu版，过程很慢，但能成功）conda install gromacs=2022.2

对蛋白配体复合物进行分子动力学模拟的操作流程

一、蛋白拓扑数据准备

安装SPDBV优化蛋白数据，将pdb文件存放在该文件夹下（应该是非必要，个人习惯存放在一起），在下图三的file按钮中进行导入，修复蛋白质结构，并导出保存为pdb格式

首先，准备蛋白的pdb文件，这里我使用的是FGSG678.pdb

gmx pdb2gmx -f FGSG678.pdb -o FGSG678_processed.gro -water spc

提示选择立场，这里选择Amber99sb-ildn力场，输入6

生成三部分数据 gro itp top

二、配体拓扑数据准备

Gromacs无法对小分子配体进行处理生成拓扑文件，这里使用Sob老师开发的处理小分子软件Sobtop，这个软件仅需本地安装，安装网址Sobtop，找到Download下载，解压安装即可，运行exe开始使用。

需要准备mol2文件，打开pymol从file打开

将pdbqt格式的分子导出成mol2格式，File —> Export Molecule，跳出save窗口，选项默认，保存的时候选择mol2格式。

对.mol2文件进行预处理：使用文本编辑器打开.mol2文件，将图中位置（只需要修改第三行阴影处以及第八列）改成配体名称，并保存mol2文件

运行sobtop，添加分子mol2的路径，添加后回车

接下来有几个选项，第一步选择1，生成top文件；第二步选择3，尽可能使用GAFF力场；第三步选择0，进入下一步；第四步选择4，如果可能，预先构建成键参数，任意猜测缺少的选项；第五步回车，生成的top文件生成在sobtop软件根目录下；第六步回车，生成的itp位置限制文件在sobtop软件根目录下；第七步，选择2回车，生成gro文件；第八步，回车，生成的gro文件在sobtop软件根目录下；回车，退出sobtop软件（一步回车不行再多按一次回车）

生成gro、itp、top文件，再把这三个文件上传至服务器

三、数据合并（！！！）

打开配体的gro文件，复制分子的坐标（第三行至倒数第二行），也就是根据分子数18，复制第三列1~18的数据。

复制蛋白的gro文件，重命名为complex.gro，在倒数第二行后插入分子拓扑坐标信息，并修改第二行的原子总数（7195+18），保存

构建拓扑，将分子的拓扑文件引入到topol.top中，在include forcefield和molecule type中间添加如下内容

在文档末尾Molecules部分加上配体名称，保存文件。

四、体系搭建

添加盒子，

gmx editconf -f complex.gro -o newbox.gro -bt dodecahedron -d 1.0

添加水模型，

gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro

运行完命令后，打开.top文件，将最后一行SOL和后面的数字成为单独一行，保存。

添加离子

 gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

选择溶剂替换分子

gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

提示选择替换溶剂分子的group时选择“SOL”(15)，得到得到solv_ions.gro文件

五、能量最小化

准备em.mdp文件，其中第二行integrator=steep代表最速下降法，把 steep改为cg就代表使用共轭梯度法。可以准备两种下降法的文件，比如最速下降法保存为em1.mdp，共轭梯度法保存为em2.mdp。

; minim.mdp - used as input into grompp to generate em.tpr
integrator  = steep     ; Algorithm (steep = steepest descent minimization; cg)
emtol       = 1000.0    ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01      ; Energy step size
nsteps      = 1000     ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
rlist           = 1.0       ; cut-off for making neighbor list (short range forces)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions (yes/no)

参考B站up方法，先使用steep最速下降法进行能量最小化 nsteps=10000，再进行cg共轭梯度法进行能量最小化nsteps=10000。这里使用同一个em.tpr就相当于在最速下降法处理完之后继续进行共轭梯度法处理了。

##（1）最速下降法 下降速度快但不收敛

gmx grompp -f em1.mdp -c solv_ions.gro -p topol.top -o em.tpr –maxwarn 7

gmx mdrun -v -deffnm em

##（2）共轭梯度法 下降速度慢但收敛

gmx grompp -f em2.mdp -c em.gro -p topol.top -o em.tpr –maxwarn 7

gmx mdrun -v -deffnm em

在使用steep进行1000次模拟之后直接进行cg模拟出现报错，错误原因有可能是初始结构中存在原子冲突，或者steep最小化次数太少，没有达到效果。

再次设置5000次的能量最小化，进行到2000多次的时候最大力就已经低于设置的阈值1000了，所以运行结束，

最小化与最大力

在这里明确两个概念，也就是衡量能量最小化的两个指标，最小化（Energy Minimization）和最大力（Maximum Force, Fmax），这两个指标关系到模拟的稳定性和结果的可靠性，

1.最小化（Energy Minimization）

定义：

最小化是通过调整原子位置，使系统势能（Potential Energy）达到局部最小值的过程。
在模拟开始时，初始结构中可能存在原子重叠（steric clashes）或不合理的几何构型（如键长/键角偏差），导致系统势能极高。直接运行动力学模拟会因能量爆炸而失败，因此需要先进行能量最小化。

2. 最大力（Maximum Force, Fmax）

定义：

最大力是系统中所有原子所受力的最大值（单位为 kJ/(mol·nm)）。
在最小化过程中，算法会迭代调整原子位置，直到最大力低于设定阈值（如 emtol = 1000 或 100），此时认为系统已收敛。

在上图我的输出结果中，输出信息显示 Steepest Descents converged to Fmax < 1000，说明最小化已收敛，但最大力仍较高（941.02）。

此时：

如果目标是初步稳定系统（例如后续进行 NVT/NPT 平衡），可以接受此结果并继续。
如果需进一步降低能量（例如发表级模拟或处理敏感体系），需继续优化。

在此基础上我又尝试了一下cg的能量最小化过程，发现还是出现报错，说明在cg前还是需要进行steep的能量最小化，因为只是想走一遍动力学模拟的流程，这个结果应该对后续影响不大，我就继续进行平衡阶段的模拟。

六、平衡阶段

1.限制配体

平衡蛋白配体复合物与平衡水中只包含蛋白的体系非常类似，值得注意的是：
① 对配体施加限定。
② 对于温度耦合群的处理。

为了限定配体，需要为他生成一个位置限定拓扑。首先，为zinctest创建一个只包含非氢原子的index group，逐步在命令窗口输入以下三行

gmx make_ndx -f jz4.gro -o index_jz4.ndx
 > 0 & ! a H*
> q

执行genrestr模块，选择新创建的index group（在index_zinctest.ndx文件中是group 3）

gmx genrestr -f zinctest.gro -n index_zinctest.ndx -o posre_zinctest.itp -fc 1000 1000 1000

选择group3

2.添加限制配体到拓扑文件中

如果我们想在无论蛋白是否被限定的条件下对配体施加限定，在拓扑文件的Include ligand topology之后中加入如下内容：

; Ligand position restraints
#ifdef POSRES
#include "posre_zinctest.itp"
#endif

把这段内容添加到这两部分之间

必须要把限定这段话加在Include ligand topology之后，Include water topology之前，加在其他任何地方都会报错。

由于对只有几个原子的组在控制其动能的涨落时, 温度耦合算法不够稳定，不要对体系中的每一单个物种使用独立的耦合。创建一个特殊的组, 其中包含蛋白质和配体。

gmx make_ndx -f em.gro -o index.ndx
> 1 | 13
> q

3.NVT平衡

创建nvt.mdp文件，将以下部分复制入文件中：

title                   = Protein-ligand complex NVT equilibration 
define                  = -DPOSRES  ; position restrain the protein and ligand
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstenergy               = 500   ; save energies every 1.0 ps
nstlog                  = 500   ; update log file every 1.0 ps
nstxout-compressed      = 500   ; save coordinates every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds to H are constrained 
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Neighbor searching and vdW
cutoff-scheme           = Verlet
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 20        ; largely irrelevant with Verlet
rlist                   = 1.2
vdwtype                 = cutoff
vdw-modifier            = force-switch
rvdw-switch             = 1.0
rvdw                    = 1.2       ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
rcoulomb                = 1.2       ; short-range electrostatic cutoff (in nm)
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling
tcoupl                  = V-rescale                     ; modified Berendsen thermostat
tc-grps                 = Protein_MOL Water_and_ions    ; two coupling groups - more accurate
tau_t                   = 0.1   0.1                     ; time constant, in ps
ref_t                   = 300   300                     ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr                = no 
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed