制作仅在xy方向有周期性的夹层结构，且不使用约束平板

转载请注明出处！！！
本篇是对这篇文章的改进，解决了xy方向的周期问题，具体原理请见上面的文章，本篇只简述过程。仿真过程用的是alpha-tubulin和beta-tubulin

制作蛋白夹水层

这里制作的收敛水层，要在z轴的负方向上给个余量，这样就可以成功进行nwall =2的模拟了。

gmx editconf -f small_water_npt.gro -c -box 5.112 4.920 5 -o small_water_npt_center.gro

然后还是用packmol做整体的结构

vi center_224.inp
packmol < center_224.inp

center_224.inp

tolerance 1.0
filetype pdb
output center_224.pdb

structure chain_A.pdb
  number 1
  inside box 0. 0. 0. 60. 60. 60.
  constrain_rotation x 0. 20.
  constrain_rotation y 0. 20.
end structure

structure small_water_npt_center.pdb
  number 1
  inside box 0. 0. 60. 60. 60. 85.
end structure

structure chain_B.pdb
  number 1
  inside box 0. 0. 85. 60. 60. 145.
  constrain_rotation x 0. 20.
  constrain_rotation y 0. 20.
end structure

然后要做成gro

gmx editconf -f small_water_npt_center.pdb -o small_water_npt_center.gro

根据packmol中的大小改gro的尺寸

  6.00000   6.00000   14.50000

再做个top

vi center_224.top

这是根据small_water.topz中有1487个sol得来的
center_224.top

#include "oplsaa.ff/forcefield.itp"
;#include "metal.itp"
;#include "mPlus.itp"
;#include "mMinus.itp"
#include "two_test_Protein_chain_A.itp"
#include "two_test_Protein_chain_B.itp"

#include "oplsaa.ff/spc.itp"
#include "oplsaa.ff/ions.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif



[ system ]
; Name
Water In Sheets and protein

[ molecules ]
; Compound       #mols

Protein_chain_A  1

SOL         1456
NA               31
Protein_chain_B  1

这里直接放离子平衡后的top了

gromacs仿真

离子平衡

gmx grompp -f ions.mdp -c center_224.gro -p center_224.top -o center_224_ions.tpr
gmx genion -s center_224_ions.tpr -o center_224_ions_finish.gro -p center_224.top -pname NA -np 31

能量最小化

gmx grompp -f em1.mdp -c center_224_ions_finish.gro -p center_224.top -o center_224_em.tpr
gmx mdrun -v -deffnm center_224_em

能量最小化之后的图

制作索引

索引在能量最小化制作，也就是nvt的前一步，这样最准确

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

最后nvt

gmx grompp -f nvt_second.mdp -c center_224_em.gro -p center_224.top -o center_224_nvt_xy.tpr
gmx mdrun -v -deffnm center_224_nvt_xy

两蛋白结合了

介电常数相关

gmx dipoles -corr total -c center_224_nvt_xy.xvg -f center_224_nvt_xy.trr -s center_224_nvt_xy.tpr
gmx dielectric -f center_224_nvt_xy.xvg -eps0 11.6619 -ffn aexp -o center_224_nvt_xy_freq.xvg -nsmooth 10

静止介电常数 11.6619
频响

相关mdp

ions.mdp

ions.mdp - used as input into grompp to generate ions.tpr 
; Parameters describing what to do, when to stop and what to save integrator  = steep        
 ; Algorithm (steep = steepest descent minimization) emtol       = 1000.0        
 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm emstep      = 0.01          
 ; Minimization step size nsteps      = 50000         
 ; 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    
 ; Buffered neighbor searching  ns_type         = grid      ; Method to determine neighbor list (simple, grid) coulombtype     = cutoff    
 ; Treatment of long range electrostatic interactions rcoulomb        = 1       
 ; Short-range electrostatic cut-off rvdw            = 1       
 ; Short-range Van der Waals cut-off pbc             = xyz       
 ; Periodic Boundary Conditions in all 3 dimensions

em1.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      = h-bonds		; No contraints except those defined in the topology

nvt_second.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		= 50000	; 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.2		; [nm] Short-range electrostatic cut-off

; Van der waals
vdwtype		= cutoff
rvdw		= 1.2		; [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
;freezedim       = Y Y Y

创作不易，转载请注明出处！！！谢谢！