Gromacs用alpha-A链和beta-B链制作电荷平衡的测试蛋白

如果直接将蛋白质不经过任何处理就作为输入使用insert插入，在进行到最后nvt平衡的grompp这一步就会报错初始能量过高，所以这里尝试使用离子平衡和能量最小化的方法先处理蛋白，然后再插入整体结构。

这里用chain-A和chain-B测试，具体结构是先将chain-A和chain-B溶于水，然后做离子平衡，将部分水替换成NA+后，每个单层蛋白质盖成电中性，然后按照这里的步骤按步骤进行。

制作alpha-A链结构

gmx pdb2gmx -f chain_A.pdb -o chain_A.gro -p chain_A.top -i chain_A.itp -water spc -ignh

生成相关文件后，用制作一个仿真盒子

gmx editconf -f chain_A.gro -o chain_A_box.gro -c -d -1.0 -bt cubic

用水溶解

gmx solvate -cp chain_A_box.gro -cs spc216.gro -p chain_A.top -o chain_A_water.gro

然后做离子平衡

 gmx grompp -f ions.mdp -c chain_A_water -p chain_A.top -o chain_A_ions.tpr

这里我观察到系统成-16电荷，所以加16个NA+离子

gmx genion -s chain_A_ions.tpr -o chain_A_ions.gro -p chain_A.top -pname NA -np 16

ions.mdp文件如下

ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to saveintegrator  = steep         
; Algorithm (steep = steepest descent minimization)emtol       = 1000.0        
; Stop minimization when the maximum force < 1000.0 kJ/mol/nmemstep      = 0.01          
; Minimization step sizensteps      = 50000         
; Maximum number of (minimization) steps to perform; Parameters describing how to find the neighbors of each atom and how to calculate the interactionsnstlist         = 1        
 ; Frequency to update the neighbor list and long range forcescutoff-scheme = Verlet   
  ; Buffered neighbor searching ns_type         = grid     
   ; Method to determine neighbor list (simple, grid)coulombtype     = cutoff    
   ; Treatment of long range electrostatic interactionsrcoulomb        = 1.0       
   ; Short-range electrostatic cut-offrvdw            = 1.0      
    ; Short-range Van der Waals cut-offpbc             = xyz       
    ; Periodic Boundary Conditions in all 3 dimensions

然后打开vmd看一下图形

效果是正确的，可以继续进行。

然后就进行能量最小化
先放上em.mdp

; Run Control
integrator	= steep		; A steepest descent algorithm for energy minimization
emtol 		= 100		; Stop minimization when the maximum force < 100.0 kJ/mol/nm
emstep 		= 0.01		; [nm] initial step-size
nsteps 		= 50000		; Maximum number of (minimization) steps to perform

; Output Control
nstlog 		= 1		; # of steps that else between writing energies to the log file
nstenergy 	= 1		; # of steps that else between writing energies to energy file

; Neighbor searching and nonbonded interaction
cutoff-scheme	= verlet	; Cutoff-related parameters
nstlist		= 1		; Frequency to update the neighbor list and the long-range forces
ns_type		= grid		; Method to determine neighbor list (simple, grid)
pbc		= xyz		; Use periodic boundary conditions in all directions
rlist		= 0.5		; [nm] Cut-off for making neighbor list (short range forces

; Electrostatics
coulombtype	= PME		; Treatment of long range electrostatic interactions
rcoulomb	= 0.5		; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype		= cutoff
rvdw		= 0.5		; [nm] Short-range Van der Waals cut-off

; Temperature and pressure coupling are off during EM
tcoupl          = no		; No temperature coupling
pcoupl          = no		; No pressure coupling (fixed box size)

; No velocities during EM 
gen_vel          = no		; Do not generate velocities

; options for bonds
constraints      = none		; No contraints except those defined in the topology

这其中rlist,rcoulomb, rvdw,这几个值需要根据蛋白链的体积调整，gromacs会提示增加增加box size or 减少rlist，对了如果只减少这些值不行的话需要调整gro文件中box size，我都调整成2了。

gmx grompp -f em.mdp -c chain_A_ions.gro -p chain_A.top -o chain_A_ions_em.tpr
gmx mdrun -v -deffnm chain_A_ions_em

然后再打开看一下

看来也是正确的，没有出现类似聚集成一团或者空洞或者管状的结构。

制作beta-B链结构

制作通过同alpha-A链，只是参数不同，这里不再重复了。

制作最终结构

首先要做个水层，这里只有参数不同，不再赘述，详细过程同这篇
我们最终得到的文件是

 gmx editconf -f small_water_npt.gro -o small_water_npt_center

然后就要添加蛋白层了

首先加下层

gmx insert-molecules -f small_water_npt_center.gro -ip position.dat -ci chain_A_ions_em.gro -o down2.gro -rot none

这个-rot none一定要加，否则结构不正确

然后先居中一下

 gmx editconf -f down2.gro -o down2.gro

然后加上层

gmx insert-molecules -f down2.gro -ip position2.dat -ci chain_B_ions_em.gro -o down3.gro -rot none

两个position是有特定要求的，太远了达不到目的，太近了无法添加进去

position.dat

0 0 -4.8

position2.dat

0 0 2.1

这样就得到了最终的gro文件


下面我们还要制作top和ndx文件。

制作top

首先要先制作两个itp文件
由于内容过多，这里只叙述原理，就是把chain-A.top和chain-B.top中的[moleculetype] [atoms]等保留，而删除top文件中描述系统分子数量的语句，具体资源上传到这里了。
然后自己制作down3.top

#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"
;#include "oplsaa.ff/spc.itp"
#include "chain_A_ions_em.itp"
#include "chain_B_ions_em.itp"
[ system ]
; Name
system

[ molecules ]
; Compound             #mols
SOL                1487
Protein_chain_A     1
SOL                 1401   
NA                  16
Protein_chain_B     1
SOL                 1862
NA                  15

这其中[ molecules ]的内容是分别从上面3个top中复制来的。

制作索引文件

使用gromacs的make_ndx制作索引文件

gmx make_ndx -f down3.gro -o down3.ndx

随后根据gromacs的提示依次选择

system

name

q

这样就完成索引文件的制作了。

尝试nvt预平衡

gmx grompp -f nvt1.mdp -c down3.gro -p down3.top -n down3.ndx -o fd_nvt.tpr

这条代码是可以成功的，说明理论上gromacs是可以做这个工作的
nvt1.mdp文件如下

; Run Control
integrator = md  ; leap-frog integrator of Newton's equations of motion
tinit  = 0  ; [ps] starting time for the run
dt  = 0.002  ; 2 fs - [ps] time step for integration
nsteps  = 5000 ; 0.1 ns

; Output Control
nstxout  = 1000  ; save coordinates every 2 ps
nstenergy = 500  ; save energies every 1 ps
nstlog  = 50000  ; update log file every 100 ps

; Bond parameters
continuation    = no  ; apply constraints to the start configuration and reset shells
constraints     = h-bonds ; convert all bonds to constraints
constraint-algorithm = lincs ; holonomic constraints 

; Neighbor searching and nonbonded interactions
cutoff-scheme = verlet ; Cutoff-related parameters
nstlist  = 20  ; Frequency to update the neighbor list and the long-range forces
ns_type  = grid  ; Method to determine neighbor list (simple, grid)
pbc  = xy  ; Use periodic boundary conditions in all directions
rlist  = 1.2  ; [nm] Cut-off for making neighbor list (short range forces)
ewald-geometry  = 3dc
nwall           = 2
wall-atomtype   = MW MW
wall-type       = 12-6

; Electrostatics
coulombtype = PME  ; Treatment of long range electrostatic interactions
rcoulomb = 1  ; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype  = cutoff
rvdw  = 1  ; [nm] Short-range Van der Waals cut-off

;                    Do I need pme-order and fourier spacing?

;pme_order = 4  ; cubic interpolation
;fourierspacing = 0.16  ; grid spacing for FFT

; Temperature coupling
tcoupl  = v-rescale ; Temperature coupling, modified Berendsen thermostat‬
tc_grps         = system ; Group to couple to separate temperature baths
tau_t           = 0.1  ; [ps] Time constant for coupling
ref_t           = 300  ; [K] Reference temperature for coupling

; Pressure coupling is off for NVT
Pcoupl          = No  ; No pressure coupling (fixed box size)

; Generate velocities to start
gen_vel  = yes  ; assign velocities from Maxwell distribution
gen_temp = 300  ; temperature for Maxwell distribution
gen_seed = -1  ; generate a random seed

; Freeze Sheets
freezegrps = Backbone ;freezing the protein, so limit the water ;
freezedim = Y Y Y

但是这个还有个问题就是目前会报错

Step 0, time 0 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
rms 0.140735, max 4.433387 (between atoms 16877 and 16880)
bonds that rotated more than 30 degrees:
 atom 1 atom 2  angle  previous, current, constraint length
 段错误 (核心已转储)

这个错误可以通过调整dt来解决。