「科研干货」一起学分子动力学模拟和结合自由能预测

一起学amber分子动力学模拟和结合自由能预测

linux操作基础

创建目录 mkdir

拷贝文件 cp

查看目录内容 ls

修改权限 chmod

检索内容 find locate

linux下的文本编辑:vi

打开文件 vi filename

插入内容i 或者 a

替换r

删除整行dd

退出:wq

计算前准备：

o 需要编辑.bashrc文件

o 添加并行环境所需的变量

source /opt/intel/composer_xe_2013/bin/compilervars.sh intel64

source /opt/intel/mkl/bin/intel64/mklvars_intel64.sh

source /opt/intel/impi/4.1.1/bin64/mpivars.sh

o 添加应用软件所需的变量

export PATH="/opt/chem/amber12/bin:$PATH"

export AMBERHOME=/opt/chem/amber12

最后运行命令source .bashrc，或者登出重新登录

做好计算准备后接下来进行

模拟二肽与环糊精主客体复合物的相互作用

以及

计算二肽与环糊精的结合自由能

1.复合物文件的预处理

o 删除文件中的H

o 删除不必要的水

o 找到文件中的CYS和HIS

– CYS可以有自由、二硫键、金属配位等三种不同形态

– HIS侧链为咪唑，可以有两个不同的质子化位点，

– 不同的状态需要修改残基名

– CYS,CYX,CYM,HID,HIE,HIP

o 利用正则表达式和grep命令可以简单地处

理复合物文件

grep -v '^..............H' ***.pdb

grep 'CYS' ***.pdb

grep 'HIS' ***.pdb

检查pdb文件，保证小分子内原子正确

生成拓扑文件和坐标文件

o 坐标文件记录模拟系统内所有原子的三维坐标

o 拓扑文件记录原子的连接关系和分子力学参数

o 需要将pdb文件转换成拓扑文件和坐标文件才能

够开展计算

o 所需工具leap

– tleap：纯命令行操作

– xleap：奇丑无比的图形界面

tleap的操作

o 登录scc主机

o 输入命令tleap

o 也可以输入tleap -f 文件名，用于读入必须

的分子力场参数文件

o source leaprc.proteins.ff03.r1

– 读入amber03力场参数，用于处理大分子

o source leparc.GLYCAM_06j-1

– 读入GLYCAM力场参数，用于处理环糊精的力场

o source leaprc.gaff

o MOL = loadmol2 ***.mol2

– 读入小分子模板

o loadamberparams ***mod

– 读入小分子模板补充的参数

o list

– 可以显示目前leap中读入的unit

o com = loadpdb ***.pdb

– 读入处理过之后的复合物文件

o saveamberparm com com.prmtop com.inpcrd

– 保存复合物的拓扑文件和坐标文件

o addions com Na+ 0

– 向复合物中添加抗衡离子，保持模拟体系电中性

o solvateoct com TIP3PBOX 12

– 为复合物添加消角八面体的周期边界溶剂盒子

o saveamberparm com cmp.prmtop cmp.inpcrd

– 保存加入抗衡离子和溶剂盒子的复合物拓扑文件和坐

标文件

o saveamberparm MOL lig.prmtop lig.inpcrd

– 保存小分子的拓扑文件和坐标文件

o 此外还需要单独读入和保存蛋白质的拓扑文件坐

标文件

可能出现的错误

o Created a new atom named: N within

residue: .R<ASP 178>

– 到pdb文件里找到对应残基号的原子，这是原子名错了

编写输入文件

Initial minimisation of our structures&cntrl

imin=1, maxcyc=2500, ncyc=1000,

cut=999, ntb=0, igb=1, saltcon=0.2,/

o Initial minimisation of our structures

–第一行是说明

o &cntrl

– 表示以下开始为输入文件的参数

o imin=1, maxcyc=2500, ncyc=1000,

– imin=1表示为能量优化计算

– maxcyc=2500表示计算2500轮

– ncyc=1000表示其中前1000轮为共阨梯度法

o cut=999, ntb=0, igb=1, saltcon=0.2,

– cut=999表示长程相互作用不设置截断值

– ntb=0表示没有周期边界

– igb=1表示使用gb模型拟合溶剂效应

– saltcon=0.2表示离子强度为0.2o /

– 表示输入文件到此为止

– 这是一个使用连续溶剂模型对结构进行优化的

输入文件其他常见键值

o imin=0, irest=1, ntx=5,

– 延续上一次模拟的输出文件继续计算

o nstlim=200000, dt=0.002,

– 总共模拟200000步，步长为2fs

o ntc=2, ntf=2,

– 应用SHAKE算法处理含H键的伸缩振动

o ntt=1, tautp=0.5,

– 应用温度耦合(NVT系综)，耦合系数为0.5ps

o tempi=310.0, temp0=310.0,

– 模拟的起止温度都是310K

o ntpr=500, ntwx=500,

– 每隔500步写入一次rst文件，每隔500步写入一次轨迹文件

下节接着分享：运行分子动力学模拟

未完待续 看完记得关注分享呀～～

有主意了