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 complex.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 complex.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，计算个氨基酸的RMSF4. 聚类分析及找代表性结构：
gmx make_ndx -f complex.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 complex.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 <<<<<<<<<<<<<<<<<<<<<<<<<\n
>> Usage: gmx_mmpbsa [-f *.xtc ] [-s *.tpr ] [-n *.ndx ]
[-pro PROTEIN] [-lig <LIGAND|none>]
[-cou dh ] [-ts ie ] [-cas <B:E:S,#> ]\n
>> Default: gmx_mmpbsa -f traj.xtc -s topol.tpr -n index.ndx
-pro Protein -lig Ligand
-ts ie
--------------------------------------------------------------------------------
>> Option:
-f: trajectory file
-s: topology 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) with Debye-Huckel screening method
-ts: ie: calculate entropy with 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. gawk version > 3; macOS awk version > 2018
!!! 3. watch _pid.pdb to make sure PBC is correct
--------------------------------------------------------------------------------
>> Log:
TODO: parallel APBS, focus
2022-07-14: fix issues with Chinease in PBSAset
2022-07-08: fix bug, No PB SA value
2022-07-01: fix bug for awk too many open files
2022-06-18: fix bug for ndx when Lig is 1st
2022-04-19: fix bug for CAS
remove -com 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 gawk
2021-05-16: revise cmd for all Linux
see https://github.com/Jerkwin/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: awk will exit for multiple system call
2020-11-15: fix bug for resID >=1000, for awk 3.x gsub /\s/
2020-06-02: fix bug of withLig
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, whole, 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(Windows), 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/vdw parmeters
pdb=_$pid.pdb # 轨迹构型文件 to save trajectory
radType=1 # 原子半径类型 radius of atoms (0:ff; 1:mBondi)
radLJ0=1.2 # 用于LJ参数原子的默认半径(A, 主要为H) radius when 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样条离散
swin 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 comps'
# 非极性计算设置(Apolar/Non-polar)
PBAset='
temp 298.15 # 温度
srfm sacc # 构建溶剂相关表面或体积的模型
swin 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 .0072*cal2J 表面张力, 常数
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 command line options
################################################################################
useDH=1; useTS=1; isCAS=0
cas="$*";
cas=${cas##*-cas}; cas=${cas%%-*}; cas=${cas// /}
opt=($*); N=${#opt[@]}
for((i=0; i<N; i++)); do
arg=${opt[$i]}; j=$((i+1)); val=${opt[$j]}
[[ $arg =~ -f ]] && { trj=$val; }
[[ $arg =~ -s ]] && { tpr=$val; }
[[ $arg =~ -n ]] && { ndx=$val; }
[[ $arg =~ -pro ]] && { pro=$val; }
[[ $arg =~ -lig ]] && { lig=$val; }
[[ $arg =~ -cou ]] && { useDH=1; }
[[ $arg =~ -ts ]] && { useTS=1; }
[[ $arg =~ -cas ]] && { isCAS=1; }
done
if [[ -z "$lig" || $lig =~ none ]]; then
withLig=0; com=$pro; lig=$pro;
else
withLig=1; com=${pro}_${lig}
fi
if [[ $step -lt 3 ]]; then
################################################################################
# 检查所需文件, 索引组
# check needed files, index group
################################################################################
[[ ! -f "$trj" ]] && { echo -e "!!! ERROR !!! trajectory File NOT Exist !\n"; exit; }
[[ ! -f "
最低0.47元/天 解锁文章
1567
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
