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//