Posted 2020-11-16Updated 2023-12-24Notes / Toolsan hour read (About 7226 words)

Amber MD关键字

之前了解了一遍Amber动力学输入文件中关键字的含义，但是几个月不用再一看又不熟悉了。为此做一个方便快速检索参考的中文笔记。

MD模拟

输入关键字

这是我使用的production步的输入文件

prod.in

Equilibration:constant pressure
 &cntrl
  imin  = 0, ntx = 5, irest = 1,    # 工作类型、读取设置
  ntf   = 2,  ntc = 2,              # SHAKE束缚
  nstlim= 50000000, dt= 0.002,      # 模拟时长
  cut   = 10.0,                     # 截断
  temp0 = 300.0,                    # 温度
  ntt   = 3, gamma_ln = 5.0,        # 温度控制
  ntb   = 2,  ntp = 1,              # 周期性边界的类型，系综
  ntpr  = 25000, ntwx = 5000,       # 输出：写入频率
  iwrap = 1,                        # 输出：坐标居中
  ig    = -1,                       # 随机化
 /

其中功能相关的设置被放置在相近的位置

imin
选择任务类型
1（最小化）
2（MD）
5（读取轨线进行再次计算）
ntx
坐标（-c指定）读取方式
1（仅读取笛卡尔坐标）
5（读取坐标和速度以进行重启）
irest
读取rst文件（-c指定）的行为
0 （不重启模拟，不读取速度，时间设置为0，开始新的模拟）
1 （重启模拟，从rst文件中读取坐标和速度，时间继续，不开始新的模拟）

ntpr
out和mdinfo文件的写入频率
每ntpr步写入一次，默认为50
ntwx
轨线文件的写入频率
每ntwx步写入一次，默认是0（不写入）
ntwr
rst文件的更新频率，减少任务崩溃的损失
每ntwr步更新一次，任务结束后也更新一次，默认值等于nstlim（只在任务结束写入一次）
iwrap
写入到轨线和rst文件中时分子的绝对位置和相对位置：
1 （写入每个分子在初始盒子中对应的等价分子，利用平移对称性只用初始盒子的大小描述体系）
0 （写入每个分子的绝对位置，往往用于扩散分析，在1的情况下无法分析扩散）

类似贪吃蛇
*不涉及扩散分析的任务都建议使用iwarp=1

ntf与ntc
SHAKE键振动限制，往往设置NTF=NTC
目的：
MD的基础步长受体系中最快的运动限制，限制最快的X-H振动以获取更大的步长
TIP3P也是在限制X-H振动的情况下拟合的，使用该力场需要加对应限制
1 （不使用SHAKE）
2 （所有含H的键被限制，在dt=0.002时必用）
3 （所有键长都被固定）
nstlim
MD的步数
dt
MD的步长，单位：ps，默认为0.001，加氢SHAKE往往使用0.002
cut
描述nonbonded截断的范围单位：埃
范围内静电作用使用精确的 QQ/r^2 范围外用Particle Mesh Ewald (PME)方法近似
有文献指出需要大于等于10.0，尤其是计算蛋白相互作用时

ntt
温度控制
0 （体系能量固定，绝热，NVE系综）
1 （恒定温度，使用weak-coupling算法：只保证总动能不变，将超出的动能重新分配。不能保证温度在分子各处一致，往往导致低频自由度能量过高，无法保证各自由度间的能量关系）
2 （恒定温度，使用Anderson-like算法：每隔vrand步就模拟碰撞，随机的重置体系中的速度以符合设定温度。会导致轨线没有绝对时间相关的意义，只有在vrand步之内可以进行绝对时间相关的分析，vrand值太小会导致构象空间的探索变慢(默认1000)）
3 （恒定温度，使用Langevin算法：使用gamma_ln设置的$\gamma$值，有时会崩溃？）
9 （恒定温度，使用优化的Isokinetic Nose-Hoover chain ensemble (OIN)方法：貌似是为了使用一些加速步进的方法？）
10（恒定温度，使用Stochastic Isokinetic Nose-Hoover RESPA integrator方法：貌似是为了使用一些加速步进的方法？Tuckerman组开发）
11（Berendsen控温器，之后用到了再补充）
nmropt
0 (默认)无nmr分析
1 读入NMR限制和权重变化（控制恒温器温度变化）
2 读入NMR限制、权重变化、NOESY volumes、化学位移、残基偶极限制
temp0
设定恒定温度值，大于300时需要降低步长

tempi
初始温度，用于heat步，初始速度使用maxwellian分布
可增加温度控制段，设置在不同阶段温度的变化，如：

heat.in >folded

&wt
 type  = 'TEMP0',
 istep1= 0, istep2=18000,
 value1= 0.0, value2=300.0,  #在前18000步升温至300K
/
&wt
 type  = 'TEMP0',
 istep1= 18001, istep2=20000,
 value1= 300.0, value2=300.0,   #在后2000步保持恒温
/
&wt
 type  = 'END',
/

ig
决定初始速度的随机种子，也决定了Langevin和Anderson-like方法中的伪随机值
-1 （时间相关的随机种子）默认
其他值 （固定种子，一般只在调试时使用）

ntp
压力控制，逐渐调整cell体积以达到设定压力
0 （无压力控制，默认）
1 （使用各向同性位置修正）
2 （使用各向异性位置修正，只能用于正交盒子，一般用于膜系统的模拟）
3 （使用半各向异性位置修正，只能用于正交盒子和恒定表面张力体系，将表面切向的两个压力耦合在一起）
pres0
恒定压力值，单位bar，默认为 1.0
ntb
周期性边界的类型, Amber20中不需要单独设置
0 （无周期性边界关闭PME 微观系综）
1 （定容）
2 （定压）
ntr
用简谐势在笛卡尔空间中限制原子位置
需要指定参考结构-ref 目标原子和力常数
0 （不限制）
1 （限制）
restraint_wt
位置限制的权重，单位：$kcal·mol^{−1}·Å^{−2}$
设置$k(\Delta x)^2$ 中的k值
restraintmask
描述目标原子，
语法：
- 描述基元
  起始符
  : 起始选择残基
  @ 起始描述原子
  序号/名称选择
  可以用 , 分隔和 - 连接的序号描述例：:1-10,12 和 @12,54-85
  可以用残基/原子名称描述例：:ARG,ALA 和 @CA,C,O
  可以用原子类型描述但是需要在 @ 后加 % 例： @%N*,N3
  可以用 :* 描述所有结构，注意！不是通配符，貌似是由于有些原子类型里包含星号如C*
  可以用 = 充当通配符，但是只能在字符之后例:@H= 描述所有H原子
- 描述基元间
  逻辑关系
  & 或 | 与 ! 非
  距离关系
  <:数字 与 >:数字 在数字范围之内/之外的所有残基
  <@数字 与 >@数字 在数字范围之内/之外的所有原子
  *可以用()描述结合的优先级
可以用ambmask测试
ambmask -p prmtop -c inpcrd -out [short| pdb| amber] -find [maskstr]
会输出选中原子到stdout

施加基于内坐标的constraint

NMR restraint的方案

Amber tutorial
http://ambermd.org/Questions/constraints.html

在.in文件中加入nmropt=1关键字，在结尾添加DISANG={path-to-rs-file}

制作.rs文件描述restraint

&rst  
iat={atom_idx},{atom_idx}, r1=抛物线势左起点, r2=抛物线势左终点, r3=抛物线势右起点, r4=抛物线势右终点, rk2=左边抛物线力常数, rk3=右边抛物线力常数,
ir6={指定一组原子时确定中心的方式}, ialtd={在势极大时的行为/是否将偏离过大的结构拖回}
&end

验证MD结果是否有效

查看out文件中的参数随轨线步数变化
使用process_mdout.perl脚本处理

perl process_mdout.pl 文件名

用cpptraj查看RMSD

cpptraj.in

parm *prmtop
trajin *nc
rms ToFirst out 名称.rmsd @CA,C,N  #restraintmask描述的目标原子
run

如果在prod步之前密度、温度、能量、rmsd都收敛，同时prod中这几项稳定
则认为prod产生的数据可以用于采样

关于提升GPU下的计算效率

Amber文档中给出的最大化GPU加速效率的建议：

避免使用过高的写入频率（小的NTPR, NTWX, NTWV, NTWE, NTWR）
避免$ntave \neq 0$
避免在不必要时使用NPT系综，考虑使用Monte Carlo barostat (barostat=2)
使用gbsa=3而不是gbsa=1，前者适合GPU
使用langevin温度控制会稍微降低速度

MMPBSA 与 MMGBSA

将结合自由能拆分成气相的结合自由能和三个物种各自的溶剂化自由能
$\Delta G_{solv,bind}^0 = \Delta G_{bind,vacuum}^0+\Delta G_{solv,complex}^0−(\Delta G_{solv,ligand}^0+\Delta G_{solv,receptor}^0)$
GB（Generalized Born）和PB（Poisson-Boltzman）描述不同的求解溶剂化自由能中静电部分的方法
使用LCPO方法求解疏水部分的贡献
J. Comput. Chem., 2016, 37, 2436–2446
关于参数选取

准备文件

需要准备三部分文件：

轨线（.nc）
描述配体、受体、复合物和溶剂化复合物的拓扑文件（.prmtop）
- 为PB/GB设置合适的PBRadii
  该设置(只作用用隐式溶剂模型)不影响常规(即使用显式溶剂)MD计算相关讨论
  可以通过set修改leap的全局设置之后生成新的prmtop文件
  set default PBRadii mbondi2
  也可以通过ParmEd修改现有的prmtop文件
  parmed -p <your_prmtop> -i input_file
  input_file中为
  parmed.in
  changeRadii mbondi2
  parmout new_radii.prmtop
  参数参考
  
  值 igb值
  
  bondi 7
  
  mbondi 1
  
  mbondi2 2 or 5
  
  mbondi3 8
- 分别为三部分结构生成不含溶剂的拓扑文件
- 同时需要了解体系的离子强度

值	igb值
bondi	7
mbondi	1
mbondi2	2 or 5
mbondi3	8

设置输入文件（.in）
格式类似MD的输入文件$AMBERHOME/bin/MMPBSA.py --input-file-help可以查询可设置选项

MMPBSA.in

Sample input file for GB and PB calculation
&general
  startframe=5, endframe=100, interval=5,
  verbose=2, keep_files=0,
/
&gb
  igb=5, saltcon=实际离子浓度,
/
&pb
  istrng=实际离子浓度, fillratio=4.0,
/

startframe & endframe 描述开始和结束的帧
interval 步长
keep_files 默认是1（保留所有轨迹文件和mdout）0（删除所有）2（保留所有）
(*关于巨大的reference.frc 重编译有些时候不太显示，希望之后amber能给个选项)
verbose 0（只输出差值部分）1（也输出组分的数据）2（增加bonded terms）

saltcon 盐浓度（单位：M）
istrng 离子强度（单位：M）
fillratio 矩形有限差分网格的最长维度与溶质的最长维度之比。对于与水尺度类似的溶剂一般不用改

gb和pb部分设置对应的参数

运行

$AMBERHOME/bin/MMPBSA.py -O -i MMPBSA.in -o FINAL_RESULTS_MMPBSA.dat -sp complex_solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y *.mdcrd

或者使用并行版

mpirun -np 核数 $AMBERHOME/bin/MMPBSA.py.MPI -O -i MMPBSA.in -o FINAL_RESULTS_MMPBSA.dat -sp complex_solvated.prmtop -cp complex.prmtop -rp receptor.prmtop -lp ligand.prmtop -y *.mdcrd

TIPS for COMET USER
在comet上没有直接的MMPBSA.py.MPI文件经过咨询了解到可以直接讲MMPBSA.py复制到自己目录下为MMPBSA.py.MPI
参考的提交脚本：

RunMMPBSA.sh >folded

> #!/usr/bin/env bash
> #SBATCH --job-name=name
> #SBATCH --partition=compute
> #SBATCH -N 1
> #SBATCH --ntasks-per-node=24
> #SBATCH --mem=100GB
> #SBATCH -t 48:00:00
> #SBATCH --signal=B:SIGUSR1@60
> #SBATCH --verbose
> #SBATCH --hint=nomultithread
> #SBATCH --mail-user=mail
> #SBATCH --mail-type=all
> 
> #Run MMPBSA
> module purge
> module load amber
> module load intel
> module load mpi4py
> 
> ibrun python3 /path/to/MMPBSA.py.MPI .... options ...
>

计算结果

计算会在当前目录下输出


*.mdcrd	GB和PB计算过程中分析过的坐标
*.mdout	所有指定帧的能量 –> 并行按后缀依次合并即可
*.pdb	平均结构的PDB文件,可进行准简谐熵计算
FINAL_RESULTS_MMPBSA.dat	结果输出文件

其中结果输出文件中会输出组分能量项：


`VDWAALS`	$H_{MM,vdw}$分子力场能的vdw部分
`EEL`	$H_{MM,ele}$分子力场能的静电部分 (不包含周期性边界条件cut=999.0)
`EPB/EGB`	$E_{solv,ele}$溶剂化自由能的静电部分，由PB或GB计算
`ESURF`	$E_{solv,nonpolar}$ ？？有什么区别 (GB only)
`ENPOLAR`	$E_{solv,nonpolar}$ 非极化作用的空穴部分 // 在mdout文件里也被写做ECAVITY (PB only)
`EDISPER`	$E_{solv,nonpolar,disper}$ 非极化作用的色散部分 (PB only)

其中并行相关的帧分配方式：将所有帧连续的分成MPI.COMM_WORLD.get_size()块之后，块的大小通过平均后渐进的方式决定

MMPBSA_mods.make_trajs.py >folded

# Sanity check
if traj.processed_frames < size: 
  raise MMPBSA_Error('Must have at least as many frames as processors!')

# We now know how many frames we have in total, so make a list that lists the
# number of frames found for each rank, and assign extra frames incrementally
frames_per_rank = traj.processed_frames // size
extras = traj.processed_frames - frames_per_rank * size
frame_count = [frames_per_rank for i in range(size)]
for i in range(size):
  if i < extras: frame_count[i] += 1

# Dump our complex trajectories
if INPUT['full_traj'] or INPUT['entropy']:
  traj.Outtraj(pre + 'complex.%s' % trj_suffix, filetype=INPUT['netcdf'])
traj.Outtraj(pre + 'complex.pdb', frames='1', filetype='pdb')
traj.Outtraj(pre + 'dummycomplex.inpcrd',frames='1',filetype='restart')

# Now dump thread-specific trajectories
last_frame = 1
for i in range(size):
  frame_string = '%d-%d' % (last_frame, last_frame+frame_count[i]-1)
  traj.Outtraj(pre + 'complex.%s.%d' % (trj_suffix, i), 
                frames=frame_string, filetype=INPUT['netcdf'])
  last_frame += frame_count[i]

根据mdout文件的逐帧分析/作图脚本
https://github.com/shaoqx/CalcKit/tree/master/Misc_scripts/MMPBSA
提取每一帧各能量项的数值，可以容易的分析和作图

MMPBSA中包含的熵

分子力场项 - 不包含熵
溶剂静电PB项 - 不包含熵
溶剂非极性部分(inp=2) - 1. 形成空穴（溶剂内聚能破坏，溶剂熵减）2. 溶质溶剂吸引（VDW作用）

MMPBSA中的周期性边界条件

查看具体的模块发现计算MM部分能量时使用的是无周期性边界条件(igb>0时默认关闭周期性边界条件ntb=0，cut=999.0)
这就意味着盒子的大小和形状会对结果产生影响（离子的位置？）
！！如果用长方形盒子必然会导致GPU在equi步报错！！

关于各种选项

根据Expert Opin Drug Discov. 2015 May 4; 10(5): 449–461.一文，总结一些在用MMPBSA研究新体系的思考：

一般可选择的有：
- PB/GB
- 3A/1A
- 关键水的处理
- 计算熵的方法
  根据文中的建议，这些选择和对研究的具体体系有一定的依赖性，最好在研究之前benchmark
采样方法：
文中建议使用多次短模拟合并采样的方法，以降低SD为标准，正在测试其可靠性
使用循环多次跑prod步
for i in {1..20}; do pmemd.cuda -O -i ./MD/prod.in -o ./MD/prod${i}.out -p *.prmtop -c ./MD/equi.rst -ref ./MD/equi.rst -r ./MD/prod_${i}.rst -x ./MD/prod_${i}.nc; done
之后用cpptraj合并
parm ../*.prmtop
trajin prod*.nc 200 last
autoimage
trajout prod.nc
go
对比；
1. 跑100ns取后90ns，每10ps一帧取9000帧分析（结果比较接近实验值）
2. 跑100ns取后90ns，每30ps一帧取3000帧分析（结果比较接近实验值）
3. 跑20条5ns 每个轨迹取后4.5ns，每10ps一帧合起来共9000帧（结果很差）
4. 跑20条5ns 每个轨迹取后4.5ns，每30ps一帧合起来共3000帧（结果很差）

非标准残基的参数化

相关参考教程
例子 http://ambermd.org/antechamber/pro4.html

从PDB中复制非天然氨基酸的几行保存为新文件(每类一个即可) 注意需要把空缺的化合价补全，一般最后是一个独立氨基酸的结构，否则antechamber难以正确判断原子类型
使用antechamber生成ac文件
antechamber -fi pdb -i {输入路径.pdb} -fo ac -o {输出路径.ac} -c bcc -at amber -nc {净电荷数}
“Note that antechamber is mostly set up to use gaff atom types, so it occasionally makes a mistake with the Amber atom types.”
这要怎么自动化啊，问题看起来是随机的，要怎么确定什么原子类型是正确的？貌似每次都会把N给错
*注意生成的ac文件会修改原子名称，会将所有数字开头的原子名称修改

制作描述和蛋白其他部分连接方式的mc文件(mianchain)

XXX.mc >folded

HEAD_NAME {连接前一个氨基酸的原子名，如N1}
TAIL_NAME {连接后一个氨基酸的原子名，如C3}
MAIN_CHAIN {从head到tail依次的原子名}
MAIN_CHAIN {...}
MAIN_CHAIN {...}
MAIN_CHAIN {...}
OMIT_NAME {需要从最终结构中删除的原子，如某些构建方式中会多出来的H}
OMIT_NAME {...}
OMIT_NAME {...}
OMIT_NAME {...}
PRE_HEAD_TYPE {前一个氨基酸用来连接改氨基酸的原子类型，如C}
POST_TAIL_TYPE {后一个氨基酸用来连接改氨基酸的原子类型，如N}
CHARGE {净电荷数}

使用prepgen生成prepin文件

prepgen -i {输入路径.ac} -o {输出路径.prepin} -m {mc文件路径.mc} -rn {残基名}

使用frcmod第一次使用天然氨基酸的参数(ff14SB)参数化

parmchk2 -i {输入路径.prepin} -f prepi -o {输出路径.frcmod} -a Y -p $AMBERHOME/dat/leap/parm/parm10.dat

将没有参数化的部分用gaff生成的参数补全

grep -v "ATTN" {路径.frcmod} > {路径.frcmod1} # Strip out ATTN lines
parmchk2 -i {路径.prepin} -f prepi -o {路径.frcmod2}

最后得到两个参数文件

使用prepin和两个frcmod文件通过tleap完成最后的参数整合

MCPB.py的使用

参考官网教程
http://ambermd.org/tutorials/advanced/tutorial20/mcpbpy.htm
相关教程
http://bbs.keinsci.com/thread-26531-1-1.html

MCPB.py [flags] [options]
-i 输入文件
-s 步骤数

目标

参数化金属中心(MetalUnit)以及附近的配体残基（非标准质子态）（非标准连接方式）

流程

大致流程是：

制作目标蛋白模型
从目标模型中提取出各种参数化相关的组件(mol2格式)
使用MCPB.py根据提供的组件自动生成Gaussian输入文件
Gaussian计算，利用MCPB.py自动读取输出文件中的信息完成参数化

I. 制作目标蛋白模型(preparation)

根据研究需求制作目标蛋白的MD-ready模型
要求：

无原子缺失(蛋白，cofactor，底物，水的轻重原子)
正确的质子态(确定配体、水的质子朝向合理)
cofactor处于合理的binding位置
底物处于合理的binding位置
包含必要的active site water

转换为amber格式的pdb并重置index

pdb4amber -i {输入路径} -o {输出路径}

II. 提取目标蛋白模型组件(mol2文件)

以及非天然氨基酸和配体的frcmod参数文件
“MCPB.py uses mol2 files based on the residue names” - 只要保证各组件残基名称和PDB中保持一致即可
“You need to create different mol2 files for different TYPES of metals, water, and NCAA in the PDB file” - 每类组件只需要一个mol2文件

金属离子
1. 从蛋白的PDB中，复制金属离子的一行保存为新文件(一个原子即可)
2. 运行metalpdb2mol2.py脚本，指定电荷，并生成mol2文件
  python metalpdb2mol2.py -i {输入路径} -o {输出路径} -c {电荷数}
Active site water 水
1. 从PDB中复制水的几行保存为新文件(一个水即可)(*需要氢和朝向正确)
2. 使用antechamber生成mol2文件
  antechamber -fi pdb -fo mol2 -i WAT_H.pdb -o WAT.mol2 -at amber -c bcc -pf y
3. 将mol2中的原子类型从oh和ho修改为OW和HW (例)
配体（cofactor和底物）
1. 从PDB中复制配体的几行保存为新文件(每类一个配体即可)
2. 使用antechamber生成mol2文件
  antechamber -fi pdb -fo mol2 -i {输入路径} -o {输出路径} -c bcc -pf y -nc {配体净电荷数}
  *注意最后保存的mol2文件的名称要和PDB中配体的残基名一致
3. 生成frcmod文件
  parmchk2 -i {输入路径} -o {输出路径} -f mol2
非天然氨基酸(modified AA)
1. 参考”非标准残基的参数化”一节
  生成prepin,frcmod1,frcmod2文件
2. 将prepin文件转换为mol2文件
  antechamber -fi prepi -i {输入路径.prepin} -fo mol2 -o {输出路径.mol2} -nc {净电荷数}
  *注意最后保存的mol2文件的名称要和PDB中配体的残基名一致
  *检查mol2中的原子类型名称和prepin中的是否一致

III. 使用MCPB.py生成Gaussian输入文件

配置MCPB计算
制作.in配置文件

mcpb.in >folded

original_pdb {完整模型路径.pdb}
group_name {自定义一个名字}
cut_off {参数化区域选择范围，默认填2.8}
ion_ids {金属的原子序号  (空格分隔可以多个,但是只能是一个金属complex里的多个)}
ion_mol2files {金属的mol2文件  (空格分隔可以多个,但是只能是一个金属complex里的多个)}
naa_mol2files {非天然氨基酸的mol2文件 *配体、水、cofactor、modified AA  (空格分隔可以多个)}
frcmod_files {补充的参数文件 *配体、cofactor、modified AA  (空格分隔可以多个)}
large_opt {是否优化大模型，默认填1 (1：只优化H、0：不优化)}
software_version {g16或g09}

不在金属中心附近的mol2和frcmod文件影响什么我不小心没加上但是也没有报错？

不在目标金属中心附近的NCAA的mol2和frcmod文件貌似不影响结果，测试了在mcpb.in中加和不加两种情况得到的最终的mol2文件和frcmod文件没有区别
貌似只影响最后生成的tleap.in文件

运行MCPB生成com文件
MCPB.py -i {路径.in} -s 1
此步会生成小模型、标准模型和大模型的PDB和fingerprint文件。还会生成小模型和大模型的高斯输入文件。

大模型的高斯输入文件将首先对氢原子进行优化，以纠正任何放置不当的氢原子。
标准模型fingerprint文件（1mlw_standard.fingerprint）包含标准模型中原子的原子类型信息（第3列是原始原子类型，最后一列是最终分配的原子类型）。金属离子的原子类型（名称以 M 开头）和连接原子（名称以 Y/Z 开头）会自动重命名，以将它们与 AMBER 力场中的其他原子类型区分开来。(为了参数化后的识别)
- 可以将 step_number 设为 1n 来强制不重命名任何原子类型并生成fingerprint文件，然后手动修改（在此之前，请检查 AMBER force field parm*.dat 文件以确保金属离子和连接原子类型不与amber已有的原子类型冲突）。
- fingerprint文件在文件末尾还包含金属离子与周围原子之间的链接信息，开头字母为“LINK”。用户还可以根据自己的需要手动修改它们（例如，如果 QM 几何优化后有任何连接更改）

运行Gaussian计算

MCPB.py生成的输入文件默认用的是B3LYP/6-31G*、2核、3000MB。建议至少改为B3LYP-D3并根据实际资源修改核数和内存设置

g16 < 1OKL_small_opt.com > 1OKL_small_opt.log

g16 < 1OKL_small_fc.com > 1OKL_small_fc.log

formchk 1OKL_small_opt.chk 1OKL_small_opt.fchk

g09 < 1OKL_large_mk.com > 1OKL_large_mk.log

这步得到的3个log文件和fchk文件需要被放置在与第一次运行MCPB.py相同的路径

*需要检查优化得到的结构中配位键是否断裂，如果断裂需尝试：

更换理论层级
先固定优化其他坐标之后放开优化配位键

IV. 使用MCPB.py读取Gaussian输出文件并生成参数

生成frcmod参数文件
在运行第一步的目录运行
MCPB.py -i 1OKL.in -s 2
*注意运行该步的时候目录下不能存在temp文件夹，该脚本貌似会生成一个temp文件
拟合电荷并生成mol2文件
MCPB.py -i 1OKL.in -s 3
默认使用charge model B
生成tleap输入文件和修改残基名后的PDB
MCPB.py -i 1OKL.in -s 4
获得
- {名称}_mcpbpy.pdb
- {名称}_tleap.in

V. 运行tleap.in完成参数化

tleap -s -f {名称}_tleap.in

需要在tleap.in文件中加入所有非天然氨基酸、配体的参数对应命令

可以改为用protein.ff14SB和water.TIP3P吗？
为什么会报错Residues lacking connect0/connect1？貌似是使用mol2和bond组合的正常表现

V-1. （对于多聚体）修改pdb和tleap以适配其他金属中心

对于多聚体，其他对称的金属中心可以使用其中一个金属中心MCPB得到的参数

但是需要修改tleap.in和PDB使tleap可以讲其识别到对应的参数

修改_mcpbpy.pdb中其他金属中心的残基名
参考MCPB.py第四步中输出的改名位点，例：

XXX.mc >folded

Renamed residues includes: 519-HID
Renamed residues includes: 521-HID
Renamed residues includes: 684-HIE
Renamed residues includes: 1441-CU
Renamed residues includes: 1445-WAT

修改其他链中对应的名称

修改_tleap.in中的bond语句
为其他链的金属中心和配体添加对应的键，例：

XXX.mc >folded

bond mol.519.NE2 mol.1441.CU
bond mol.521.NE2 mol.1441.CU
bond mol.684.ND1 mol.1441.CU
bond mol.1441.CU mol.1445.O
bond mol.518.C mol.519.N
bond mol.519.C mol.520.N
bond mol.520.C mol.521.N
bond mol.521.C mol.522.N
bond mol.683.C mol.684.N
bond mol.684.C mol.685.N
bond mol.1239.NE2 mol.1446.CU
bond mol.1241.NE2 mol.1446.CU
bond mol.1404.ND1 mol.1446.CU
bond mol.1446.CU mol.1450.O
bond mol.1238.C mol.1239.N
bond mol.1239.C mol.1240.N
bond mol.1240.C mol.1241.N
bond mol.1241.C mol.1242.N
bond mol.1403.C mol.1404.N
bond mol.1404.C mol.1405.N

V-2. （对于多个不同金属中心）修改pdb和tleap.in

MCPB.py新加的原子类型名称都使用默认名称导致重复，需要更改重复的原子类型名称以适配多次MCBP生成的mol2和frcmod文件
- 将重复的M1,Y1等原子类型重命名之后合并到tleap.in中
- 在frcmod, mol2文件中查找替换成修改后的原子类型名称
将mol2文件和frcmod文件合并至tleap.in中
在此步注意如果有重复的残基名（例如有两个不同的Ca2+中心时都会被命名为CA1）需要在mol2文件和tleap.in都更改为彼此独特的残基名（如CA2）
**注意 CA2 = loadmol2 xxxx.mol中等号左边的名字符合和pdb配对一定要修改
将bond指令合并至tleap.in中
将更新后的残基名合并到同一个pdb文件中
如果两个不同的金属中心彼此比较接近/有共享的配体
则难以以此方式合并参数，会引起一系列的参数确实错误
有些时候对于辅助折叠的金属中心可以使用nonbonded model来描述

Misc

注意antechameber生成的mol2文件中的原子类型可能不规范 pymol可能会不识别从而转存的时候将原子名称搞错

教程中提到需要修复原子序号的问题，但是没有修复也不怎么影响模拟，跑了100ns都没问题

轨线输出文件格式

mdcrd： ASCII编码的纯坐标每帧由中心的三个坐标组成的一行分割帧内由连续的每个原子的xyz坐标组成（每行10个数值）

amber.mdcrd >folded

Cpptraj Generated trajectory                                                    
  17.357  12.865  14.673  17.869  12.805  13.804  16.361  12.909  14.514  17.650
  12.025  15.151  17.683  14.091  15.410  17.493  14.988  14.820  17.057  14.204
  ...
  33.227  24.489   1.071  33.789  24.515   1.846  32.350  24.321   1.417
  58.637  58.637  58.637
  ...
  [next frame]

NetCrd(.nc/.crd): 二进制的坐标文件
可以通过cpptraj将这些格式以及pdb相互转换

Misc

tleap

tleap自动加氢甚至不会考虑周围的金属离子(如Zn2+等配位中心)，会把旁边的残基质子化而不考虑配位

prepin格式

https://ambermd.org/doc/prep.html
利用树和LOOP的方式描述连接

GPU在equi步报错

用solvatebox的话就会用solvateOct则不会

迷惑的内存错误

Error: an illegal memory access was encountered launching kernel kClearForces
cudaMemcpy GpuBuffer::Download failed an illegal memory access was encountered

参考

https://blog.csdn.net/Litedg/article/details/124230876
http://archive.ambermd.org/202203/0149.html 说了部分原理

主要的原因：

GPU相关的MD精度有限，在涉及巨大的数值的时候会出现该问题
往往意味着结构中有bad contact从而存在巨大力

解决方案：

用CPU跑min步可以解决大部分问题
有一次遇到CPU跑20000步min还遇到这个问题改成100000步就好了

Amber MD关键字

MD模拟

输入关键字

施加基于内坐标的constraint

NMR restraint的方案

验证MD结果是否有效

关于提升GPU下的计算效率

MMPBSA 与 MMGBSA

准备文件

运行

计算结果

MMPBSA中包含的熵

MMPBSA中的周期性边界条件

关于各种选项

非标准残基的参数化

MCPB.py的使用

目标

流程

I. 制作目标蛋白模型(preparation)

II. 提取目标蛋白模型组件(mol2文件)

III. 使用MCPB.py生成Gaussian输入文件

IV. 使用MCPB.py读取Gaussian输出文件并生成参数

V. 运行tleap.in完成参数化

V-1. （对于多聚体）修改pdb和tleap以适配其他金属中心

V-2. （对于多个不同金属中心）修改pdb和tleap.in

Misc

轨线输出文件格式

Misc

tleap

prepin格式

GPU在equi步报错

迷惑的内存错误

Catalogue

Categories