提取GMX轨迹中两个group之间的静电和范德华相互作用

使用Gromacs进行分子动力学模拟之后,有时我们需要从轨迹文件中提取两个group之间的静电和范德华相互作用,可以使用如下方法:
下面以计算128号残基和228号残基之间的相互作用为例,介绍本人曾尝试过的三种方法及遇到的问题。

  1. 使用 g_energy

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    a) 先计算单独的group 

    i=128 or i=228

    # 预先创建包含残基128、残基228以及两个残基为一个整体的index文件
    echo r_$i | trjconv -s md.tpr -f md.xtc -n inter.ndx -o inter-${i}.xtc -b 60000 -e 100000
    echo r_$i | tpbconv -s md.tpr -n inter.ndx -o del-${i}.tpr
    mdrun -rerun inter-${i}.xtc -s inter-${i}.tpr -o inter-${i}.trr -e inter-${i}.edr -g inter-${i}.log
    echo -e "LJ-(SR) \n Coulomb-(SR)" | g_energy -f inter-${i}.edr -o inter-${i}.xvg |tail -2 > inter-r${i}.log
    awk 'NR==1 {print $3}' inter-r${i}.log > LG-SR-${i}.dat
    awk 'NR==2 {print $3}' inter-r${i}.log > Coulomb-SR-${i}.dat

    b) 然后计算两个group之和

    echo 13 r_128_228 | trjconv -s md.tpr -f md.xtc -n inter.ndx -o inter-2.xtc -b 60000 -e 100000
    echo r_128_228 | tpbconv -s md.tpr -n inter.ndx -o inter-2.tpr
    mdrun -rerun inter-2.xtc -s inter-2.tpr -o inter-2.trr -e inter-2.edr -g inter-2.log
    echo -e "LJ-(SR) \n Coulomb-(SR)" | g_energy -f inter-2.edr -o inter-2.xvg |tail -2 > inter-2.log
    awk 'NR==1 {print $3}' inter-2.log > LG-SR-2.dat
    awk 'NR==2 {print $3}' inter-2.log > Coulomb-SR-2.dat

    c) 二者相减

    实际使用时,如果用来计算蛋白与小分子之间的相互作用,则计算值与MMPBSA(介电常数选1)值相近,但用来计算两个残基之间的相互作用时,计算的结果数值特别大!不知为何

  2. 使用 g_mmpbsa

    1
    2
    3
    4
    5
    i=128

    echo r_128 r_228 | g_mmpbsa -f md_pbc.xtc -s mmpbsa.tpr -n del.ndx -pdie 2 -decomp -b 60000 -e 100000 -mm res-${i}.xvg -mmcon contrib_res-${i}.dat
    sed -i '1,25d' res-${i}.xvg
    awk '{sum+=$8} END {print "Average = ", sum/NR}' res-${i}.xvg > total-${i}.dat

    实际使用时注意介电常数的选择

  3. 使用 Energy group
    在 md.mdp 文件中的 energygrps 选项中添加感兴趣的group,然后进行rerun,即可在生成的edr文件中写入对应的group之间的相互作用信息,进而通过g_energy命令提取出来。

    实际使用时,如果用来计算蛋白与小分子之间的相互作用,则计算值与MMPBSA(介电常数选1)值相近。这个应该是相对准确的一种吧。