双层膜体系模型的构建

对于双层膜这种较复杂的体系，使用GROMACS中自带的简单工具构建较为繁琐。因此选取Packmol用来构建双层膜部分，同时利用GROMACS中的gmx solvate命令来向体系中加入溶剂，这样组合使用能大大节省构建时间。我的体系的双层膜部分的Packmol输入文件如下，没有很复杂，基本与Packmol官方教程一致。但是应注意resnumbers选项的添加，否则生成的结构文件中分子序号将会不连续。

<pre>

construction of SDS-C16-CHOL-SOL system

tolerance 2.0
filetype pdb
output C16DS-CHOL128.pdb
structure SDS-C16min.pdb
number 64
resnumbers 2
inside box 0. 0. 0. 60. 60. 20.
atoms 1 45
over plane 0. 0. 1. 16.
end atoms
atoms 17
below plane 0. 0. 1. 3
end atoms
atoms 64
below plane 0. 0. 1. 1.
end atoms
end structure

structure SDS-C16min.pdb
number 64
resnumbers 2
inside box 0. 0. -20. 60. 60. 0.
atoms 1 45
below plane 0. 0. 1. -16.
end atoms
atoms 17
over plane 0. 0. 1. -3
end atoms
atoms 64
over plane 0. 0. 1. -1.
end atoms
end structure
</pre>

快速之余，这种方法也会带来一点小问题，就是在溶剂化的过程中会有一些溶剂分子进入到双层膜部分，通常采用手动删除或用脚本删除，GROMACS官网里有分享一个自动删除不合理位置溶剂（水）分子的脚本，现将具体操作教程翻译一下，并附上我个人使用时遇到的一些问题及解决方法：

膜模拟

运行膜模拟

GROMACS用户往往会在运行双脂质层模拟时遇到各种问题，尤其是体系中还有蛋白质的时候。推荐模拟双脂质层的用户参考教程KALP-15 in DPPC

完成一次对膜蛋白体系的模拟需要执行一下步骤：

1. 为你选用的蛋白质和膜分子选取合适的力场参数;
2. 将蛋白质插入到膜结构中;（例如，在预备好的双分子层上使用g_membed命令或做一个粗粒化自组装模拟然后再转回全原子）
3. 向体系中加入溶剂并加入离子确保整个体系最终呈电中性;
4. 运行能量最小化;
5. 让膜与蛋白质相适应，在限制蛋白质所有重原子的情况下运行约5～10ns的MD;
6. 去掉限制进行平衡;
7. 运行MD;

使用genbox命令加入水

当使用genbox（5.0以后版本中为gmx solvate）向双层膜体系中添加水时，水分子往往会进入到双层膜中，一般来说解决办法有一下几种：

• 将vdwradii.dat文件从$GMXLIB拷贝到当前目录下，调大脂质层分子的范德华半径（建议调节碳原子的半径于0.35～0.5nm）以阻止genbox命令将水分子加入这些空隙;
• 运行一段短暂的MD，让脂质层分子的憎水作用把水分子排开; * 手动的去删除这些位置上的水分子;
• 借助本教程中的脚本来修改;

增大水相的z轴坐标

假设整个双层膜体系沿z轴展开，这个由Chirs Neale提供的脚本可以排除不需要的水分子，具体操作如下（如果需要，可以在第41行修改以适配脚本用于x、y方向展开的体系）：

1. 对初始双层膜构型initial.gro运行genbox命令，得到solvated.gro;
2. cp solvate.gro new_waters.gro;
3. 删除new_waters.gro中所有不为新添加进来的水的内容;（若初始构型中有水分子，则这些水分子也要删去）
4. 修改脚本keepbyz.sh中的upperz和lowerz参数为需要的值，并将sol参数改为溶剂名称;（upperz和lowerz分别为双层膜的上下边界的z轴坐标，sol对应溶剂的名称，默认为SOL）
5. ./keepbyz.sh new_waters.gro > keep_these_waters.gro;（先执行chmod +x keepbyz.sh确保脚本可以执行）
6. tail -1 initial.gro > last_line.gro;
7. head -$(expr $(cat initial.gro | wc -l | awk '{print $1}') - 1 ) initial.gro > not_last_line.gro;
8. cat not_last_line.gro new_waters.gro last_line.gro > new_system.gro;
9. editconf -f new_system.gro -o new_system_sequential_numbers.gro;（5.0以后版本中为gmx editconf)

脚本内容如下：
<pre>

!/bin/bash

give x.gro as first command line arguement

upperz=6.417
lowerz=0.820
sol=SOL
count=0 cat $1 | grep "$sol" | while read line; do
for first in $line; do
break
done
if [ "$count" = 3 ]; then
count=0
fi
count=$(expr $count + 1)
if [ "$count" != 1 ]; then
continue
fi
l=${#line}
m=$(expr $l - 24) // would use -48 if velocities are also in .gro and -24 otherwise
i=1
for word in ${line:$m}; do
if [ "$i" = 1 ]; then
popex=$word
else
if [ "$i" = 2 ]; then
popey=$word
else
if [ "$i" = 3 ]; then
popez=$word
break
fi
fi
fi
i=$(expr $i + 1)
done
nolx=echo "$popez > $upperz" | bc
nohx=echo "$popez < $lowerz" | bc
myno=$(expr $nolx + $nohx)
if [ "$myno" != 0 ]; then
z=${#first}
if [ "$z" != 8 ]; then
sfirst="[[:space:]]$first"
else
sfirst=$first
fi
echo grep $sfirst $1
fi
done
</pre>

脚本中默认使用3原子的水分子作为溶剂，若使用TIP4P,则应将：
<pre>
if [ "$count" = 3 ];
then count=0
fi
</pre>

修改为：
<pre>
if [ "$count" = 4 ];
then count=0
fi
</pre>

同理，TIP5P应修改为5。为提高效率,该脚本的作者还提供了一个同样功能的C程序，详见原文链接

这里列出我遇到的问题来供大家参考：

1. 按照我的理解上述教程的第八步应该为:cat not_last_line.gro keep_these_waters.gro last_line.gro > new_system.gro，因为keep_these_waters.gro才是删除了错误位置的水分子后的构型;
2. 在执行第九步时会出现"bad box error"，这是因为第七步产生的gro文件中第二行的原子数与总体系的总原子数不相等，这里可以手动更改一下;
3. 提醒一下新手，脚本中第十九行//后的内容为注释，告诉你如果gro文件中包含速度信息，这一数值应改为48。bash中注释不以//打头，在执行前可以删去;
4. 在vim中删除多行可以用指令ndd，n为行数;
5. 修改完分子数后记得修改拓扑文件，确保两者中分子数一致;

**上文中的脚本已改进！更高效的教程请参考李老师博客GROMACS之如何做膜模拟! **