分子动力学模拟学习3-Gromacs数据处理

1. 周期性边界处理

gmx trjconv -s md.tpr -f md.xtc -o md_noPBC.xtc -pbc mol -center #选择1("Protein")作为置中的组分，选择0("System")作为输出

2. 计算RMSD

gmx make_ndx -f plex.gro -o rmsd_index.ndx #新建用于计算RMSD的ndx文件，并在此运行下输入下列命令

4 & r1-303 #新建4 & 待分析氨基酸

name 28 Enzyme_backbone #将新建组命名为Enzyme_backbone 

28 | r 270 | r275 #再新建Enzyme_backbone | 需要额外添加的氨基酸组，如非标准氨基酸（因为这样包括了非标准氨基酸的全部原子而只有非骨架原子，并不确定能否这样做，无需计算backbone而计算全部氨基酸的RMSD则没有这种情况）

q #保存并退出

gmx rms -n rmsd_index.ndx -s md.tpr -f md_noPBC.xtc -o rmsd.xvg -tu ns #期间选择Enzyme_backbone | r 270 | r 275组，以ns的时间间隔运行RMSD计算

3. 计算RMSF

gmx make_ndx -f plex.gro -o rmsf_index.ndx #和上面类似，新建ndx文件

3  r 1-303 #新建待分析氨基酸组ca组

name 28 Enzyme_ca #重命名

28 | r 270 | 275 #新建包括非标准氨基酸的组 

q #保存并退出

gmx rmsf -n rmsf_index.ndx -f md_noPBC.xtc -s md.tpr -o rmsf-per-residue.xvg -ox average.pdb -oq bfactors-residue.pdb -res #选择组29，计算个氨基酸的RMSF

4. 聚类分析及找代表性结构

gmx make_ndx -f plex.gro -o cluster_index.ndx #新建进行cluster的ndx文件，并在其运行下执行一下命令

r 1-311 #新建待分析的氨基酸组

name 28 Enzyme #重命名新建组28为Enzyme

q #保存并退出

gmx cluster -f md_noPBC.xtc -s md.tpr -cutoff 0.15 -method gromos -n cluster_index.ndx -cl cluster.pdb -g cluster.log #选择组28，选择合适的cutoff（此例中时0.15）保证获得的cluster数量合适，最终获得cluster.pdb，第一个构象最多的cluster就是代表性构象

5. gmx_mmpbsa计算结合能

gmx trjconv -f md_noPBC.xtc -o trj_50-100ns.xtc -b 50000 -e 100000 #轨迹采样，本例共模拟了100 ns，取50-100ns

gmx check -f trj_50-100ns.xtc #轨迹检查

gmx make_ndx -f plex.gro -o mmpbsa_index.ndx #新建用于计算的ndx组，包括Enzyme组和ligand组，具体操作和上述相似，就不多说了

之后计算主要依赖李继存老师的gmx_mmpbsa.bsh的脚本，可以参考以下内容结合自由能计算MM/PBSA

本例中的 gmx_mmpbsa.bsh脚本内容如下,可供参考

echo -e "
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   gmx_mmpbsa   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Jicun Li <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>>>>>>>>>>>  2022-07-14 11:17:41  <<<<<<<<<<<<<<<<<<<<<<<<>   Usage: gmx_mmpbsa [-f   .xtc  ] [-s   .tpr  ] [-n   .ndx  ]
        [-pro PROTEIN] [-lig ]
        [-cou dh  ] [-ts ie   ] [-cas  ]n
>> Default: gmx_mmpbsa -f   traj.xtc  -s   ol.tpr  -n   index.ndx
        -pro Protein   -lig Ligand
            -ts ie
--------------------------------------------------------------------------------
>> Option:
   -f: trajectory file
   -s: ology file
   -n: index file
 -pro: index group name of protein
 -lig: index group name of ligand, can be ignored using none
 -cou: dh: calculate MM(COU) ith Debye-Huckel screening method
  -ts: ie: calculate entropy ith interaction entropy method
 -cas: select res for CAS, begin:end:step or number, using global resnum
   Other: change them in the script directly
--------------------------------------------------------------------------------
>> Warning:
   !!! 1. make sure gmx is available and version is consistent
   !!! 2. gak version > 3; macOS ak version > 2018
   !!! 3. atch _pid.pdb to make sure PBC is correct
--------------------------------------------------------------------------------
>> Log:
   TODO:    parallel APBS, focus
   2022-07-14: fix issues ith Chinease in PBSAset
   2022-07-08: fix bug, No PB SA value
   2022-07-01: fix bug for ak too many open files
   2022-06-18: fix bug for ndx hen Lig is 1st
   2022-04-19: fix bug for CAS
      remove - option
   2022-02-09: CAS except GLY, PRO, CYX
      fix radius bug, use AMBER mBondi
      change to AMBER PB4
   2021-11-25: keep original PDB residue name
      output avaraged contibutions
      clean output files
   2021-05-17: add systime and fflush, only for gak
   2021-05-16: revise cmd for all Linux
      see https://github./Jerkin/gmxtool/issues/2
   2021-03-14: add -cou, -ts option, see 10.1088/0256-307X/38/1/018701
      change default PB method to npbe
   2021-03-03: fix bug: ak ill exit for multiple system call
   2020-11-15: fix bug for resID >=1000, for ak 3.x gsub /s/
   2020-06-02: fix bug of ithLig
   2020-06-01: fix bug of -skip
   2020-05-27: fix bug for sdie
   2020-05-26: fix bug for RES name
   2020-04-03: use C6, C12 directly
   2020-01-08: support protein only
   2019-12-24: fix bug for small time step
   2019-12-10: fix bug for OPLS force field
   2019-11-17: fix bug for c6, c12 of old version tpr
   2019-11-03: fix bug for time stamp
   2019-09-19: push to gmxtool
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>n"
################################################################################
# 设置运行环境, 计算参数
# setting up environmets and parameters
################################################################################
trj=trj_50-100ns.xtc		# 轨迹文件 trajectory file
tpr=md.tpr		# tpr文件  tpr file
ndx=mmpbsa_index.ndx		# 索引文件 index file
pro=Enzyme			# 蛋白索引组   index group name of protein
lig=ligand			# 配体索引组   index group name of ligand
step=0	# 从第几步开始运行 step number to run
		# 1. 预处理轨迹: 复合物完整化, 团簇化, 居中叠合, 然后生成pdb文件
		# 2. 获取每个原子的电荷, 半径, LJ参数, 然后生成qrv文件
		# 3. MM-PBSA计算: pdb->pqr, 输出apbs, 计算MM, APBS
		# 1. pre-processe trajectory, hole, cluster, center, fit, then generate pdb file
		# 2. abstract atomic parameters, charge, radius, C6/C12, then generate qrv file
		# 3. run MM-PBSA, pdb->pqr, apbs, then calculate MM, PB, SA
gmx='gmx'								# /path/to/GMX/bin/gmx_mpi
dump="$gmx dump"						# gmx dump
trjconv="$gmx trjconv -dt 1000"			# gmx trjconv, use -b -e -dt, NOT -skip
apbs='apbs'			# APBS(Windos), USE "/", NOT ""
pid=pid				# 输出文件($$可免重复) prefix of the output files($$)
err=_$pid.err		# 屏幕输出文件 file to save the message from the screen
qrv=_$pid.qrv		# 电荷/半径/VDW参数文件 to save charge/radius/vd parmeters
pdb=_$pid.pdb		# 轨迹构型文件 to save trajectory
radType=1			# 原子半径类型 radius of atoms (0:ff; 1:mBondi)
radLJ0=1.2			# 用于LJ参数原子的默认半径(A, 主要为H) radius hen LJ=0 (H)
meshType=0			# 网格大小 mesh (0:global  1:local)
gridType=1			# 格点大小 grid (0:GMXPBSA 1:psize)
cfac=3				# 分子尺寸到粗略格点的放大系数
					# Factor to expand mol-dim to get coarse grid dim
fadd=10				# 分子尺寸到细密格点的增加值(A)
					# Amount added to mol-dim to get fine grid dim (A)
df=0.5				# 细密格点间距(A) The desired fine mesh spacing (A)
# 极性计算设置(Polar)
PBEset='
  temp  298.15   # 温度
  pdie  2     # 溶质介电常数
  sdie  78.54    # 溶剂介电常数, 真空1, 水78.54

  npbe     # PB方程求解方法, lpbe(线性), npbe(非线性), smbpe(大小修正)
  bcfl  mdh   # 粗略格点PB方程的边界条件, zero, sdh/mdh(single/multiple Debye-Huckel), focus, map
  srfm  smol  # 构建介质和离子边界的模型, mol(分子表面), smol(平滑分子表面), spl2/4(三次样条/7阶多项式)
  chgm  spl4  # 电荷映射到格点的方法, spl0/2/4, 三线性插值, 立方/四次B样条离散
  sin  0.3   # 立方样条的窗口值, 仅用于 srfm=spl2/4

  srad  1.4   # 溶剂探测半径
  sdens 10    # 表面密度, 每A^2的格点数, (srad=0)或(srfm=spl2/4)时不使用

  ion charge  1 conc 0.15 radius 0.95  # 阳离子的电荷, 浓度, 半径
  ion charge -1 conc 0.15 radius 1.81  # 阴离子

  calcforce  no
  calcenergy ps'
# 非极性计算设置(Apolar/Non-polar)
PBAset='
  temp  298.15 # 温度
  srfm  sa   # 构建溶剂相关表面或体积的模型
  sin  0.3 # 立方样条窗口(A), 用于定义样条表面

  # SASA
  srad  1.4 # 探测半径(A)
  gamma 1   # 表面张力(kJ/mol-A^2)

  #gamma const 0.0226778 3.84928
  #gamma const 0.027  0
  #gamma const 0.0301248 0   # AMBER-PB4 .0072cal2J 表面张力, 常数

  press  0  # 压力(kJ/mol-A^3)
  bconc  0  # 溶剂本体密度(A^3)
  sdens 10
  dpos  0.2
  grid  0.1 0.1 0.1

  # SAV
  #srad  1.29   # SAV探测半径(A)
  #press 0.234304  # 压力(kJ/mol-A^3)

  # WCA
  #srad   1.25     # 探测半径(A)
  #sdens  200   # 表面的格点密度(1/A)
  #dpos   0.05     # 表面积导数的计算步长
  #bconc  0.033428    # 溶剂本体密度(A^3)
  #grid   0.45 0.45 0.45 # 算体积分时的格点间距(A)

  calcforce no
  calcenergy total'
################################################################################
# 检查 gmx, apbs 是否可以运行
# check gmx, apbs availability
################################################################################
str=$($gmx --version | grep -i "GROMACS version")
[[ $step -lt 3 && -z "$str" ]] && { echo -e "!!! ERROR !!!  GROMACS NOT available !n"; exit; }

str=$($apbs --version 2>&1 | grep -i "Poisson-Boltzmann")
[[ -z "$str" ]] && { echo -e "!!! WARNING !!!  APBS NOT available !n"; }

################################################################################
# 解析命令行参数
# parse mand line options
################################################################################

useDH=1;  useTS=1; isCAS=0

cas="$";
cas=${cas##-cas}; cas=${cas%%-}; cas=${cas//