gmx蛋白每个碱基与材料相互作用能计算

CruiseBend

老师们好，在gmx全原子计算蛋白吸附在二维材料时，如何计算每一个碱基和蛋白的vdW和Coul力？一般计算interaction energy，都是设置tc_grp，然后rerun轨迹，gmx energy提取LJ和Coul进行计算。但是一个蛋白碱基比较多，肯定不能一个一个设置group跑，如何操作比较方便呢？好像gmx不支持设置多个group然后rerun

sobereva

gmx当然支持设置多个能量组来rerun
当前问题和tc_grp毫无关系，要设的显然是energygrps

我的gmx培训班（http://www.keinsci.com/workshop/KGMX_content.html）的两页ppt






不要随便用力字，仔细看
计算化学中的一些常见不良写法和用词
http://sobereva.com/298（http://bbs.keinsci.com/thread-1358-1-1.html）

CruiseBend

是的老师，是energygrp。但如果要统计每个碱基和材料的相互作用，就只能对每个碱基重复energygrp=residue nano，grompp，rerun，energy的流程是吗？这样会很麻烦，是否有什么scripts能方便这一流程呢？
PS：我想得到类似这样的图

sobereva

CruiseBend 发表于 2022-2-13 20:39
是的老师，是energygrp。但如果要统计每个碱基和材料的相互作用，就只能对每个碱基重复energygrp=residue n ...

自己写shell脚本循环就完了，没什么难的

shell脚本编写的一些基础知识看
详谈Multiwfn的命令行方式运行和批量运行的方法
http://sobereva.com/612（http://bbs.keinsci.com/thread-24929-1-1.html）

CruiseBend

#!/bin/bash

# make sure ^Pro900.gro #Ene.mdp topol.top——ps ^Pro_900Bp.xtc in the last file
# edit line

conf=Pro900.gro
mn=b
mkdir rerun$mn
cd rerun$mn

ma=GRA
cp /home/junfan/zyliao/eneref/np.inp np.inp
gmx_mpi make_ndx -f ../$conf -o index.ndx < np.inp
if [ $mn == 'a' ]; then
        cp /home/junfan/zyliao/eneref/Enea.mdp Enea.mdp
  sed -i 's/MAT/GRA/' Enea.mdp
else
        cp /home/junfan/zyliao/eneref/Eneb.mdp Eneb.mdp
        sed -i 's/MAT/GRA/' Eneb.mdp
fi
gmx_mpi grompp -f Ene$mn.mdp -o Ene$mn.tpr -c ../$conf -r ../$conf -p ../../topol.top -n index.ndx
gmx_mpi mdrun -v -deffnm Ene$mn -rerun ../Pro_900$ma.xtc
echo "1000\n" | gmx_mpi energy -f Ene$mn.edr -b 200000 > output.txt 2>output2.txt


for ((i=1; i<=85;i=i+1))
do
nu=$(grep "LJ-SR:r_$i-$ma" output2.txt | awk -F LJ-SR:r_$i-$ma '{print $1}' | awk '{print $NF}')
echo "$nu\n" | gmx_mpi energy -f Ene$mn.edr -b 200000 > output.txt 2> /dev/null
t1=$(grep "LJ-SR" output.txt | awk '{print $2}')
t2=$(grep "LJ-SR" output.txt | awk '{print $4}')
echo "$i $t1 $t2" >> result.txt
rm -r output.txt
done

老师建议很好，很快的就写了出来，供大家参考

CruiseBend

老师好，我使用rerun计算了氮化硼片和蛋白的相互作用，发现Coul interaction很大，有+100KJ/mol。这是怎末回事呢？一般BN即使有电荷，他和蛋白的Coul力应该很小接近0，主要是pi-pi作用，体现在vdW interaction。
请老师帮忙解答，谢谢。

Kelly00

sobereva 发表于 2022-2-13 03:38
gmx当然支持设置多个能量组来rerun
当前问题和tc_grp毫无关系，要设的显然是energygrps

老师，我按照这个方法计算的组内的非键相互作用项，其中组内的静电相互作用总是表现为很大的正值（静电排斥），我算了好几组物质都是，感觉没有道理。想问一下老师计算像r_2-r_2这样组内的静电作用，是仅仅包含组内各个分子间的非键相互作用吗，是不是就是每两个分子间的静电相互作用相加的累计和。

sobereva

Kelly00 发表于 2023-2-20 11:57
老师，我按照这个方法计算的组内的非键相互作用项，其中组内的静电相互作用总是表现为很大的正值（静电排 ...

假设你用的是cutoff方式算静电作用，你这么定义得到的是r_2内每一对原子间静电相互作用的总和（1-4静电作用另计）
对于用Ewald、PME等情况，没法这么分解

Kelly00

sobereva 发表于 2023-2-21 02:28
假设你用的是cutoff方式算静电作用，你这么定义得到的是r_2内每一对原子间静电相互作用的总和（1-4静电作 ...

谢谢老师，也就是用cutoff这样算组内静电，实际上也包括了同一个分子的两两原子间的静电作用是吗？那同样用cutoff算范德华，也会包括同个分子的原子间的范德华作用？

sobereva

Kelly00 发表于 2023-2-23 16:23
谢谢老师，也就是用cutoff这样算组内静电，实际上也包括了同一个分子的两两原子间的静电作用是吗？那同样 ...

包含
是