烷基磺酸钠分子层以及融入蛋白后的所有流程"/>
重复双层十二烷基磺酸钠分子层以及融入蛋白后的所有流程
制作水的夹层
制作的xoy面积和文献一样5.217 nm平方，一共有4877个水分子
gmx grompp -f em_water.mdp -c water_330.gro -p water_330.top -o water_330_em.tpr
gmx mdrun -v -deffnm water_330_em
依次能量最小化、nvt、npt
设计到的文件都是water_330标注的，这里粘几个mdp和重要的top，水是用solvate做的spc216.gro，然后top里面写spce
#include "oplsaa.ff/forcefield.itp"#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz1 1 1000 1000 1000
#endif[ system ]
; Name
water 330[ molecules ]
; Compound #mols
SOL 4877em_water.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
;define = -DPOSRES
; 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
energygrps = system
; 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 = 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; 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 = all-bonds ; No contraints except those defined in the topology
;constraint-algorithm = lincs ; holonomic constraints ;freezegrps = a_15238-15804
;freezedim = Y Y Ynvt_water.mdp
title = Protein-ligand complex NVT equilibration
;define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 100000 ; 1 * 100 = 100 ps
dt = 0.001 ; 1 fs
; Output control
nstxout = 1000 ; save coordinates every 1.0 ps
nstvout = 1000 ; save velocities every 1.0 ps
nstenergy = 1000 ; save energies every 1.0 ps
nstlog = 1000 ; update log file every 1.0 ps
energygrps = system
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1 ; short-range electrostatic cutoff (in nm)
rvdw = 1 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 298 ; 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;ewald-geometry = 3dc
;nwall = 2
;wall-atomtype = MW MW
;wall-type = 12-6; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 298 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seednpt_water.mdp
title = Protein-ligand complex NPT equilibration
;define = -DPOSRES -DPOSRES_LIG ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 1000000 ;20ns
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
energygrps =
; Bond parameters
continuation = yes ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.4 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling
pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; velocity generation off after NVT 然后由生成的water_330.npt的相关文件计算静态介电常数为69.65
gmx dipoles -corr total -c water_npt_330.xvg -f water_npt_330.trr -s water_npt_330.tpr
gmx dielectric -f water_npt_330.xvg -eps0 69.65 -ffn aexp -o water_npt_330_freq.xvg
xmgrace water_npt_330_freq
C12H25SO3Na制作分子层
用十二烷基磺酸钠做的分子层，分为上层和下层，上层和下层都是亲水端对着水，而整体的结构结构尽量垂直于水的表面，实现的方法是用gaussview画的时候在z轴方向注意连接H的方式，然后用gromacs复制156个，插入的时候虽然gro文件不同，但是itp用一个，itp用LigParGen自动生成后修正0.0001的电荷误差，全部加到S原子身上。
上下层的itp分别是
sds_3.itp和sds_4.itp，但是我们只有sds_3.itp
红色是亲水端 C12H25SO3-，然后将NA加入水中，参考文献中的方法
做好了itp和gro文件，就是插入分子
这里我最后让仿真盒子的大小大概是10.4
但是后面蛋白好像要到20才能有正确结构，所以现在的工作正在尝试h-bonds和20的盒子如何能重复
我做的时候是nvt过后才editconf -c -box，但是也可以在这里就这样做，我就直接写在这里了
计算gro中盒子尺寸的方法是把水的盒子x大小/13，y轴/6，因为以后跑起来就紧密了，所以无所谓大小，主要是数量。
制作分子层
gmx genconf -f sds_3.gro -nbox 13 6 1 -o sds_3_level.gro
gmx genconf -f sds_4.gro -nbox 13 6 1 -o sds_4_level.gro
居中并扩大盒子，能放下分子层，同时也为了后面考虑z轴的镜像作用
gmx_mpi editconf -f water_330_npt.gro -c -box 5.16905 5.16905 10.4 -o water_330_npt_center.gro
然后插入
gmx insert-molecules -f water_330_npt_center.gro -ip position1.dat -ci sds_3_level.gro -o one_sds.gro -rot nonegmx insert-molecules -f one_sds.gro -ip position2.dat -ci sds_4_level.gro -o two_sds.gro -rot none
再做top two_sds.top
#include "oplsaa.ff/forcefield.itp"
#include "sds_3.itp"#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz1 1 1000 1000 1000
#endif[ system ]
; Name
water 324 ch4[ molecules ]
; Compound #mols
SOL 4721
NA 156H 156离子平衡，加入na+
gmx grompp -f ions.mdp -c two_sds.gro -p two_sds.top -o two_sds_ions.tprgmx genion -s two_sds_ions.tpr -o two_sds_ions.gro -p two_sds.top -pname NA -np 156到这里分子层就做好了，接下来就是依次能量最小化、nvt、npt过程
#!/bin/bash
#SBATCH -p gsc # Queue
#SBATCh -N 2 # Node count required for the job
#SBATCH -n 48 # Number of tasks to be launched
#SBATCH -J jobname # Job name
#SBATCH -o %J.out # Standard output
#SBATCH -e %J.err # File in which to store job error messagesexport PATH=/THFS/opt/gcc/7.3.0/bin:/THFS/opt/anaconda/3/2019.3/bin:/THFS/home/nuc_whf/software/2020.6/bin:$PATH
export LD_LIBRARY_PATH=/THFS/opt/gcc/7.3.0/lib64/:/THFS/opt/gcc/7.3.0/lib:$LD_LIBRARY_PATH
source /THFS/opt/intel/parallel_studio_xe_2018_update2/bin/compilervars.sh intel64
source /THFS/home/nuc_whf/software/2020.6/bin/GMXRC
export LD_LIBRARY_PATH=/THFS/home/nuc_whf/software/2020.6/lib64:$LD_LIBRARY_PATHyhrun -N 2 -n 48 /THFS/home/nuc_whf/software/2020.6/bin/gmx_mpi grompp -f /THFS/home/nuc_whf/sds/npt_sds_10ns_1.2r.mdp -c /THFS/home/nuc_whf/sds/two_sds_nvt_xy_10.gro -p /THFS/home/nuc_whf/sds/two_sds.top -o /THFS/home/nuc_whf/sds/two_sds_npt_be_10ns_1.2r.tpr
yhrun -N 2 -n 48 /THFS/home/nuc_whf/software/2020.6/bin/gmx_mpi mdrun -v -deffnm /THFS/home/nuc_whf/sds/two_sds_npt_be_10ns_1.2r
相关mdp
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
;define = -DPOSRES
; 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
energygrps = system
; 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 = 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; 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 = all-bonds ; No contraints except those defined in the topology
;constraint-algorithm = lincs ; holonomic constraints ;freezegrps = a_15238-15804
;freezedim = Y Y Ynvt_define2.mdp
title = Protein-ligand complex NVT equilibration
define = -DPOSRES ; position restrain the protein and ligand
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 100000 ; 1 * 100 = 100 ps
dt = 0.001 ; 1 fs
; Output control
nstxout = 1000 ; save coordinates every 1.0 ps
nstvout = 1000 ; save velocities every 1.0 ps
nstenergy = 1000 ; save energies every 1.0 ps
nstlog = 1000 ; update log file every 1.0 ps
energygrps = system
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1 ; short-range electrostatic cutoff (in nm)
rvdw = 1 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 298 ; 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;ewald-geometry = 3dc
;nwall = 2
;wall-atomtype = MW MW
;wall-type = 12-6; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 298 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
npt_20ns.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 = 10000000 ; 20 ns; Output Control
nstxout = 5000 ; save coordinates every 10 ps
nstenergy = 500 ; save energies every 1 ps
nstlog = 50000 ; update log file every 100 ps; Bond parameters
constraints = h-bonds ; convert all bonds to constraints
constraint-algorithm = lincs ; holonomic constraints
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; 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 = xyz ; Use periodic boundary conditions in all directions
rlist = 1.2 ; [nm] Cut-off for making neighbor list (short range forces); Electrostatics
coulombtype = PME ; Treatment of long range electrostatic interactions
pme_order = 4 ;
fourierspacing = 0.12
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; 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 = 298 ; [K] Reference temperature for coupling; Pressure coupling is off for NVT
Pcoupl = berendsen ; Pressure coupling on in NPT
pcoupltype = semiisotropic
tau_p = 0.5 ; [ps] time constant for coupling
ref_p = 1.0 1.0 ; [bar] reference pressure for coupling
compressibility = 0.0 4.5e-05 ; [bar^-1] isothermal compressibility of water注意我当时做的是em和nvt都是xy周期的，而没有居中盒子，然后npt是一样的，但是所有的约束都all-bonds！
为了和溶剂蛋白后的mdp相同才改的h-bonds，目前的问题是造成很慢的平衡速度，可以就不做水的h-bonds了，直接做蛋白的h-bonds，如果证明只是平衡时间的问题话。
最终的结果
计算结果42.1287
298K mdp设置与他一样
跑20ns 需要调整盒子的尺寸
经过测试20ns前变化比较大，大概18ns后开始趋于平稳，20ns是必须保证的
左边是文献，右边是自己的，平衡时间越长压的越紧。
讨论蛋白的仿真流程
拟进行的过程是
em_water nvt_define npt_define npt_water
加入分子层后 em1 nvt_define2 npt_define2 npt_20ns
目前准备去掉npt_define2，如果直接可以跑的话，deifne都是用之前nvt_38和npt_xyz改的，只是加了define和h-bonds，其他的没有变，现在就等水算出来然后如果证明是仿真时间的问题就直接跑那4个就可以了。
alpha beta的特殊过程
要用top做itp，具体过程就是去掉下面位置限定后的部分和删除力场文件，然后在双层的top中加itp，这里给一个1zvb的示例
two_1zvb.top
#include "oplsaa.ff/forcefield.itp"
#include "sds_3.itp"
#include "1zvb_Protein_chain_A.itp"
#include "1zvb_Protein_chain_B.itp"
#include "1zvb_Protein_chain_C.itp"
#include "oplsaa.ff/spce.itp"
#include "oplsaa.ff/ions.itp"#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz1 1 1000 1000 1000
#endif[ system ]
; Name
water 324 ch4[ molecules ]
; Compound #mols
Protein_chain_A 1
Protein_chain_B 1
Protein_chain_C 1
SOL 28
SOL 36
SOL 35
SOL 4128
NA 156
NA 3
H 156还有就是chain_A的过程想扩大一下水的体积，否则genion加不够156个na。
更多推荐
重复双层十二烷基磺酸钠分子层以及融入蛋白后的所有流程
发布评论