Amber中利用nmode模块进行熵计算

最近在对一个系统进行MD计算模拟，在对生成的轨道做能量分析时，突然遇到点问题，也就是对系统熵的计算。主要原因也就是自己系统太大，没法进行熵计算。因此，就需要将系统中关键部位的残基拿出来构建一个小系统。我看目前许多发表的文章也是按照这个思路在做。如何定义关键部位呢，通常选取小分子周围8-12A范围的残基作为关键残基，并与小分子一起拿出来共同构建一个小系统。

     我具体的做法如下：

     1. 确定关键残基的编号。由于先前在利用tleap准备sander的输入文件.top和.crd后，会将系统中残基和原子的编号重新从1开始计数，因此这里的我们所确定的关键残基编号也是经tleap处理完重排后的编号。具体实现可以从tleap生成的系统PDB文件中获取，获取的时候记录编号顺序一定不能变。另外我们还需要将原系统中肽链N端和C端的残基编号加入到关键残基编号列表中，这一步非常重要，原因见后面。

     2. 生成新的rec.pdb，prot_lig.pdb和new.trj轨道文件。基于已有的编号列表，按照顺序利用脚本从tleap生成的系统PDB文件中提取坐标构成新的PDB文件，同时获得所包含的原子编号列表。基于所包含的原子编号列表，从先前生成的系统old.trj轨道文件中提取这些原子的坐标，生成新的new.trj轨道文件，格式一定要正确。

     这里强烈建议将生成的new.trj文件仔细检查下，一个简单的方法就是利用pymol将后面生成的新prot_lig.top和该new.trj文件导入，并查看残基运动坐标是否正确。

     3. 生成输入文件.top。因为amber10中mmpbsa.pl需要lig.top, rec.top 和 prot_lig.top三个文件，因此我们就得再次利用tleap对所生成新的rec.pdb和prot_lig.pdb文件进行处理。在计算的过程中需要用到lig.lib文件，这个在进行模拟阶段已经生成，导入即可，同时lig.top文件也已经生成。

     注意：这里一个经常遇到的问题就是tleap对肽链两端残基的处理。由于N-端和C-端残基的特殊性（通常会含有不同类型氢原子），因此在利用tleap处理的时候，如果那肽链内部的残基直接作为两端残基进行处理的话，tleap会报错，当然我们可以对报错残基进行人为的修改，但一个简单的办法就是直接将原来系统中的肽链端残基加入新系统同样作为两端残基就可以避免出错了，同时也可以实现程序的批量处理。

     4. 修改mmpbsa.in文件。修改相应mmpbsa.pl参数文件mmpbsa.in的相应参数，生成snapshots,并进行GB, PB和nmode计算。