lingyuzeng/gromacs_docker
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

big update

Browse Source
This commit is contained in:
root
2024-10-05 14:41:52 +08:00
parent dbd3e452be
commit 06e5a0f94b
1 changed files with 202 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files

305
script/S112D/S112D.sh
View File

@@ -10,13 +10,14 @@ if [ -n "$GMXRC_PATH" ]; then
fi fi
GPU_ID=${GPU_ID:-"0"} GPU_ID=${GPU_ID:-"0"}
LIGAND_GRO="ligand.gro" # 配体的 .gro 文件 NAME="protein"
LIGAND_ITP="ligand.itp" # 配体的 .itp 文件 LIGAND_GRO="mole.gro" # 配体的 .gro 文件
LIGAND_ITP="mole.itp" # 配体的 .itp 文件
# FORCEFIELD=${FORCEFIELD:-"amber99sb-ildn"} # FORCEFIELD=${FORCEFIELD:-"amber99sb-ildn"}
FORCEFIELD=${FORCEFIELD:-"amber99sb"} FORCEFIELD=${FORCEFIELD:-"amber99sb"}
WATERMODEL=${WATERMODEL:-"tip3p"} WATERMODEL=${WATERMODEL:-"tip3p"}
WATERTOPFILE=${WATERTOPFILE:-"spc216.gro"} WATERTOPFILE=${WATERTOPFILE:-"spc216.gro"}
BOXTYPE=${BOXTYPE:-"tric"} BOXTYPE=${BOXTYPE:-"tric"} # or dodecahedron
BOXORIENTATION=${BOXORIENTATION:-"1.0"} BOXORIENTATION=${BOXORIENTATION:-"1.0"}
# NSTEPS=${NSTEPS:-500000} # 50,000 steps for 1 ns # NSTEPS=${NSTEPS:-500000} # 50,000 steps for 1 ns
NSTEPS=${NSTEPS:-50000000} # 100 ns 模拟 NSTEPS=${NSTEPS:-50000000} # 100 ns 模拟
@@ -62,11 +63,11 @@ gmx_mpi pdb2gmx -f $NAME.pdb -o $NAME.gro -ff $FORCEFIELD -water $WATERMODEL -p
# 使用 Sobtop（http://sobereva.com/soft/Sobtop）提取 # 使用 Sobtop（http://sobereva.com/soft/Sobtop）提取
# Step 2: Merge protein and ligand .gro files into a complex # Step 2: Merge protein and ligand .gro files into a complex
cat $NAME.gro mol.gro > complex.gro cat $NAME.gro mole.gro > complex.gro
# 获取蛋白质和小分子的原子数 # 获取蛋白质和小分子的原子数
PROTEIN_ATOMS=$(head -n 2 $NAME.gro | tail -n 1 | tr -d ' ') PROTEIN_ATOMS=$(head -n 2 $NAME.gro | tail -n 1 | tr -d ' ')
LIGAND_ATOMS=$(head -n 2 mol.gro | tail -n 1 | tr -d ' ') LIGAND_ATOMS=$(head -n 2 mole.gro | tail -n 1 | tr -d ' ')
# 打印原子数，确保正确 # 打印原子数，确保正确
echo "Protein atoms: $PROTEIN_ATOMS" echo "Protein atoms: $PROTEIN_ATOMS"
@@ -75,41 +76,23 @@ echo "Ligand atoms: $LIGAND_ATOMS"
# 计算总原子数 # 计算总原子数
TOTAL_ATOMS=$(($PROTEIN_ATOMS + $LIGAND_ATOMS)) TOTAL_ATOMS=$(($PROTEIN_ATOMS + $LIGAND_ATOMS))
# 打印总原子数 # 打印总原子数 7482
echo "Total atoms: $TOTAL_ATOMS" echo "Total atoms: $TOTAL_ATOMS"
# Update the atom count in the complex.gro file # 更新原子数目 # ! 手动实现
sed -i "2s/\(^[[:space:]]*\)[0-9]\+/\1$TOTAL_ATOMS/" complex.gro sed -i "2s/\(^[[:space:]]*\)[0-9]\+/\1$TOTAL_ATOMS/" complex.gro
echo -e "0\nq" | gmx_mpi genrestr -f mol.gro -o posre_mol.itp # 定义盒子、添加溶剂和离子、能量最小化
cat << EOF >> mol.itp
#ifdef POSRES
#include "posre_mol.itp"
#endif
EOF
# Modify the topology file to include the ligand
# echo "#include \"$LIGAND_ITP\"" >> topol.top
# Step 6: Modify the topology file to include the ligand and its parameters
# 6.1 在topol.top顶部添加小分子的参数文件
sed -i '/^; Include forcefield parameters/a #include "jz4.prm"' topol.top
# 6.2 在position restraint后添加小分子的拓扑文件
sed -i '/^; Include Position restraint file/!b;n;/#endif/a ; Include ligand topology\n#include "jz4.itp"' topol.top
# 6.3 在[molecules]部分添加小分子
sed -i '/^Protein_chain_A/a JZ4 1' topol.top
# 定义盒子并向盒子充满水
# Step 3: Define the simulation box # Step 3: Define the simulation box
gmx_mpi editconf -f complex.gro -o complex-box.gro -bt $BOXTYPE -c -d $BOXORIENTATION gmx_mpi editconf -f complex.gro -o complex-box.gro -bt $BOXTYPE -c -d $BOXORIENTATION
# gmx_mpi editconf -f complex.gro -o complex-box.gro -bt dodecahedron -c -d 1.0
# Step 4: Add solvate # Step 4: Add solvate
gmx_mpi solvate -cp complex-box.gro -cs $WATERTOPFILE -o complex-solv.gro -p topol.top gmx_mpi solvate -cp complex-box.gro -cs $WATERTOPFILE -o complex-solv.gro -p topol.top
# gmx_mpi solvate -cp complex-box.gro -cs spc216.gro -o complex-solv.gro -p topol.top
# 添加离子
# add icons # ! ions.mdp add by manual # add icons # ! ions.mdp add by manual
# --- ions.mdp file content --- # # --- ions.mdp file content --- #
cat << EOF > ions.mdp cat << EOF > ions.mdp
@@ -129,7 +112,14 @@ rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
EOF EOF
gmx_mpi grompp -f ions.mdp -c complex-solv.gro -p topol.top -o ions.tpr -maxwarn 1 gmx_mpi grompp -f ions.mdp -c complex-solv.gro -p topol.top -o ions.tpr -maxwarn 1
# 添加NaCl
echo SOL | gmx_mpi genion -s ions.tpr -o complex-solv-ions.gro -p topol.top -pname NA -nname CL -conc 0.125 -neutral echo SOL | gmx_mpi genion -s ions.tpr -o complex-solv-ions.gro -p topol.top -pname NA -nname CL -conc 0.125 -neutral
# -pname 阳离子的名称
# -nname 阴离子的名称
# -np 阳离子个数
# -nn 阴离子个数
# 能量最小化
# energy minimization of the structure in solvate # ! minim.mdp add by manual # energy minimization of the structure in solvate # ! minim.mdp add by manual
# --- minim.mdp file content --- # # --- minim.mdp file content --- #
cat << EOF > minim.mdp cat << EOF > minim.mdp
@@ -150,46 +140,117 @@ pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
EOF EOF
gmx_mpi grompp -f minim.mdp -c complex-solv-ions.gro -p topol.top -o em.tpr gmx_mpi grompp -f minim.mdp -c complex-solv-ions.gro -p topol.top -o em.tpr
gmx_mpi mdrun -v -deffnm em gmx_mpi mdrun -v -deffnm em
# optional em, you will need the Xmgrace plotting too # optional em, you will need the Xmgrace plotting too
#gmx_mpi energy -f em.edr -o potential.xvg #gmx_mpi energy -f em.edr -o potential.xvg
#position restrain # 查看系统的系统的总势能变化情况
echo "12\n" | gmx_mpi energy -f em.edr -o potential.xvg
# 系统的总势能变化绘图
./plot_potential.sh potential.xvg potential.svg
# 位置约束（Position Restraint）
# 位置约束是指在进行模拟时，约束某些原子的位置，使它们不会偏离初始坐标。在能量最小化之后的位置约束步骤用于对溶剂和其他环境进行平衡，同时保持蛋白质或其他分子的位置不变。通常用于溶液系统的准备过程，以避免对生物分子结构造成过大影响。
# 位置约束步骤解释
# posre.mdp 文件是位置约束的输入文件，它会约束特定原子的运动。
# grompp -r em.gro：此命令中的 -r em.gro 参数指明位置约束应基于已经能量最小化的结构 em.gro 进行。
# mdrun -deffnm posre：执行模拟时，会生成位置约束下的模拟结果。
# 如何使用
# 首先准备好 posre.mdp 文件，然后执行以下命令：
# gmx_mpi grompp -f posre.mdp -c em.gro -p topol.top -o posre.tpr -r em.gro # gmx_mpi grompp -f posre.mdp -c em.gro -p topol.top -o posre.tpr -r em.gro
# gmx_mpi mdrun -v -deffnm posre # gmx_mpi mdrun -v -deffnm posre
# 这将生成一个应用位置约束后的系统，用于后续的分子动力学模拟步骤。
# 位置限制的目的: 加速体系平衡, 减少计算资源浪费
# 对配体施加位置限制（position restraint），通常是为了在分子动力学模拟的早期阶段保持配体在特定位置，以帮助系统更快达到平衡，尤其是在 NVT（恒体积恒温）和 NPT（恒压恒温）平衡的过程中。这种约束有助于避免在系统能量未平衡时，配体位置大幅偏离而影响蛋白-配体复合物的稳定性。
# 位置限制通常在以下情况中使用：
# 能量最小化后，NVT 和 NPT 平衡之前：
# 在模拟的平衡阶段（NVT 和 NPT），施加配体限制可以保证配体在溶剂环境中不产生大幅度的偏移，特别是在系统温度和压力稳定之前。这时配体的限制作用是关键的。
# 平衡后，生产阶段去除限制：
# 在生产阶段（即实际的分子动力学生产模拟），通常应去除位置限制，以观察蛋白-配体复合物的真实动态行为。如果限制保留在生产阶段，配体的运动将受到不必要的约束，无法得到真实的配体运动轨迹。
# 限制可做可不做的情况：
# 如果体系比较简单，且配体与蛋白质的相互作用较强，位置限制可以选择不做，直接让配体在系统中自由运动。
# 如果配体是比较松散或弱结合到蛋白质上的，或者你关心配体的初期运动轨迹，初期的短时间平衡阶段可以不施加限制。
# 限制配体
# 平衡蛋白配体复合物与平衡水中只包含蛋白的体系非常类似，值得注意的是：
# 1、对配体施加限定。
# 2、对于温度耦合群的处理。
# 为了限定配体，需要为他生成一个位置限定拓扑。首先，为JZ4创建一个只包含非氢原子的index group。
gmx_mpi make_ndx -f mole.gro -o index_mole.ndx << EOF
0 & ! a H*
q
EOF
# 执行genrestr模块，选择新创建的index group（在index_mole.ndx文件中是group 3）
# echo -e "0\nq" | gmx_mpi genrestr -f mol.gro -o posre_mol.itp
gmx_mpi genrestr -f mole.gro -n index_mole.ndx -o posre_mole.itp -fc 1000 1000 1000 << EOF
3
EOF
# 下一步需要把这个信息包括在拓扑文件中。有几种方式，取决于我们希望使用的条件。这里介绍第一种：如果我们想在无论蛋白是否被限定的条件下对配体施加限定，在拓扑文件的Include ligand topology之后中加入如下内容：
# 在拓扑文件的Include ligand topology之后中加入如下内容
: << EOF
; Ligand position restraints
#ifdef POSRES
#include "posre_mole.itp"
#endif
EOF
# 必须要把限定这段话加在Include ligand topology之后，Include water topology之前，加在其他任何地方都会报错。
# 由于对只有几个原子的组在控制其动能的涨落时, 温度耦合算法不够稳定，不要对体系中的每一单个物种使用独立的耦合。创建一个特殊的组, 其中包含蛋白质和配体。
gmx_mpi make_ndx -f em.gro -o index.ndx << EOF
1 | 13
q
EOF
# NVT、NPT平衡
# nvt # nvt
# gmx_mpi grompp -f nvt.mdp -c posre.gro -p topol.top -o nvt.tpr # gmx_mpi grompp -f nvt.mdp -c posre.gro -p topol.top -o nvt.tpr
# --- nvt.mdp file content --- # # --- nvt.mdp file content --- #
cat << EOF > nvt.mdp cat << EOF > nvt.mdp
title = OPLS Lysozyme NVT equilibration title = Protein-ligand complex NVT equilibration
define = -DPOSRES ; position restrain the protein define = -DPOSRES ; position restrain the protein and ligand
; Run parameters ; Run parameters
integrator = md ; leap-frog integrator integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs dt = 0.002 ; 2 fs
; Output control ; 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 nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters ; Bond parameters
continuation = no ; first dynamics run continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy lincs_order = 4 ; also related to accuracy
; Nonbonded settings ; Neighbor searching and vdW
cutoff-scheme = Verlet ; Buffered neighbor searching cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet nstlist = 20 ; largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff （in nm） rlist = 1.2
rvdw = 1.0 ; short-range van der Waals cutoff （in nm） vdwtype = cutoff
DispCorr = EnerPres ; account for cut-off vdW scheme vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics ; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
pme_order = 4 ; cubic interpolation pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on ; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate tc-grps = Protein_MOL Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310 310 ; reference temperature, one for each group, in K ref_t = 310 310 ; reference temperature, one for each group, in K
; Pressure coupling is off ; Pressure coupling is off
@@ -201,118 +262,156 @@ gen_vel = yes ; assign velocities from Maxwell distributio
gen_temp = 300 ; temperature for Maxwell distribution gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed gen_seed = -1 ; generate a random seed
EOF EOF
gmx_mpi grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr gmx_mpi grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -n index.ndx
gmx_mpi mdrun -v -deffnm nvt gmx_mpi mdrun -v -deffnm nvt
# optional : Let's analyze the temperature progression, again using energy: # optional : Let's analyze the temperature progression, again using energy:
# gmx_mpi energy -f nvt.edr -o temperature.xvg # 提取并保存温度数据
# npt gmx_mpi energy -f nvt.edr -o temperature.xvg << EOF
17 39 40
EOF
# 拆分 .xvg 文件
# 提取整体温度
grep -v '^@' temperature.xvg | grep -v '^#' | awk '{print $1, $2}' > system_temperature.xvg
# 提取 Protein_MOL 组的温度
grep -v '^@' temperature.xvg | grep -v '^#' | awk '{print $1, $3}' > protein_mol_temperature.xvg
# 提取 Water_and_ions 组的温度
grep -v '^@' temperature.xvg | grep -v '^#' | awk '{print $1, $4}' > water_and_ions_temperature.xvg
# 绘制温度变化曲线
./plot_potential.sh system_temperature.xvg system_temperature.svg
./plot_potential.sh protein_mol_temperature.xvg protein_mol_temperature.svg
./plot_potential.sh water_and_ions_temperature.xvg water_and_ions_temperature.svg
# NPT平衡
# gmx_mpi grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr # gmx_mpi grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
# --- npt.mdp file content --- # # --- npt.mdp file content --- #
cat << EOF > npt.mdp cat << EOF > npt.mdp
title = OPLS Lysozyme NPT equilibration title = Protein-ligand complex NPT equilibration
define = -DPOSRES ; position restrain the protein define = -DPOSRES ; position restrain the protein and ligand
; Run parameters ; Run parameters
integrator = md ; leap-frog integrator integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs dt = 0.002 ; 2 fs
; Output control ; 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 nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps nstlog = 500 ; update log file every 1.0 ps
nstxout-compressed = 500 ; save coordinates every 1.0 ps
; Bond parameters ; Bond parameters
continuation = yes ; Restarting after NVT continuation = yes ; continuing from NVT
constraint_algorithm = lincs ; holonomic constraints constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy lincs_order = 4 ; also related to accuracy
; Nonbonded settings ; Neighbor searching and vdW
cutoff-scheme = Verlet ; Buffered neighbor searching cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme nstlist = 20 ; largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm) rlist = 1.2
rvdw = 1.0 ; short-range van der Waals cutoff (in nm) vdwtype = cutoff
DispCorr = EnerPres ; account for cut-off vdW scheme vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics ; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on ; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate tc-grps = Protein_MOL Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on ; Pressure coupling
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupl = Parrinello-Rahman ; Parrinello-Rahman or C-rescale for NPT production
pcoupltype = isotropic ; uniform scaling of box vectors pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
refcoord_scaling = com refcoord_scaling = com
; Periodic boundary conditions ; Periodic boundary conditions
pbc = xyz ; 3-D PBC pbc = xyz ; 3-D PBC
; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = no
; Velocity generation ; Velocity generation
gen_vel = no ; Velocity generation is off gen_vel = no ; velocity generation off after NVT
EOF EOF
gmx_mpi grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr gmx_mpi grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -n index.ndx
gmx_mpi mdrun -v -deffnm npt gmx_mpi mdrun -v -deffnm npt
# Optional: Let's analyze the pressure progression, again using energy: type 18 0 # Optional: Let's analyze the pressure progression, again using energy: type 18 0
# gmx energy -f npt.edr -o pressure.xvg gmx_mpi energy -f npt.edr -o pressure.xvg << EOF
18
0
EOF
# Optional: Let's take a look at density as well, this time using energy and entering "24 0" at the prompt. # Optional: Let's take a look at density as well, this time using energy and entering "24 0" at the prompt.
# gmx energy -f npt.edr -o density.xvg gmx_mpi energy -f npt.edr -o density.xvg << EOF
24
0
EOF
# 绘制 SVG 图像
# 绘制压力图像
./plot_potential.sh pressure.xvg pressure.svg
# 绘制密度图像
./plot_potential.sh density.xvg density.svg
# md # md
# --- md.mdp file content --- # # --- md.mdp file content --- #
cat << EOF > ${MDRUN_NAME}.mdp cat << EOF > ${MDRUN_NAME}.mdp
title = OPLS Lysozyme NPT equilibration title = Protein-ligand complex MD simulation
; Run parameters ; Run parameters
integrator = md ; leap-frog integrator integrator = md ; leap-frog integrator
nsteps = ${NSTEPS} ; steps for simulation nsteps = 5000000 ; 2 * 5000000 = 10000 ps (10 ns)
dt = ${DT} ; time step in fs dt = 0.002 ; 2 fs
; Output control ; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 5000 ; save energies every 10.0 ps nstenergy = 5000 ; save energies every 10.0 ps
nstlog = 5000 ; update log file every 10.0 ps nstlog = 5000 ; update log file every 10.0 ps
nstxout-compressed = 5000 ; save compressed coordinates every 10.0 ps nstxout-compressed = 5000 ; save coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters ; Bond parameters
continuation = yes ; Restarting after NPT continuation = yes ; continuing from NPT
constraint_algorithm = lincs ; holonomic constraints constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained constraints = h-bonds ; bonds to H are constrained
lincs_iter = 1 ; accuracy of LINCS lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy lincs_order = 4 ; also related to accuracy
; Neighborsearching ; Neighbor searching and vdW
cutoff-scheme = Verlet ; Buffered neighbor searching cutoff-scheme = Verlet
ns_type = grid ; search neighboring grid cells ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme nstlist = 20 ; largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff （in nm） rlist = 1.2
rvdw = 1.0 ; short-range van der Waals cutoff （in nm） vdwtype = cutoff
vdw-modifier = force-switch
rvdw-switch = 1.0
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; Electrostatics ; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2
pme_order = 4 ; cubic interpolation pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on ; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate tc-grps = Protein_MOL Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310 310 ; reference temperature, one for each group, in K ref_t = 310 310 ; reference temperature, one for each group, in K
; Pressure coupling is on ; Pressure coupling
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupl = Parrinello-Rahman ; pressure coupling is on for NPT
pcoupltype = isotropic ; uniform scaling of box vectors pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1 compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions ; Periodic boundary conditions
pbc = xyz ; 3-D PBC pbc = xyz ; 3-D PBC
; Dispersion correction ; Dispersion correction is not used for proteins with the C36 additive FF
DispCorr = EnerPres ; account for cut-off vdW scheme DispCorr = no
; Velocity generation ; Velocity generation
gen_vel = no ; Velocity generation is off gen_vel = no ; continuing from NPT equilibration
EOF EOF
# Generate GROMACS .tpr file for the simulation # Generate GROMACS .tpr file for the simulation
gmx_mpi grompp -f ${MDRUN_NAME}.mdp -c npt.gro -t npt.cpt -p topol.top -o ${TPR_FILE} gmx_mpi grompp -f ${MDRUN_NAME}.mdp -c npt.gro -t npt.cpt -p topol.top -o ${TPR_FILE} -n index.ndx
# Run the simulation # Run the simulation
gmx_mpi mdrun -deffnm ${MDRUN_NAME} -update gpu -ntomp 1 -gpu_id ${GPU_ID} gmx_mpi mdrun -deffnm ${MDRUN_NAME} -update gpu -ntomp 1 -gpu_id ${GPU_ID}