新冠阳了,今天开始复建运动。先试着把小分子的力场文件构建起来吧。这是我并不熟悉的领域,所以需要学习好多教程。
获得PDB结构
这里我们主要是要制作醋酸盐(acetate)和甲铵盐(methylammonium)的力场,我们首先需要得到小分子的PDB结构,可以从PDB库中找到我们需要的小分子。醋酸盐和甲铵盐的代码分别为ACT和3P8。在网页上可以下载到.sdf
文件,接下来可以使用openbabel将.sdf
文件转化为.pdb
坐标。 这样就获得了对应的.pdb
文件。
此时由openbabel生成的acetate的.pdb
文件如下: 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18COMPND ACT
AUTHOR GENERATED BY OPEN BABEL 3.1.1
HETATM 1 C UNL 1 -0.072 0.000 0.000 1.00 0.00 C
HETATM 2 O UNL 1 -0.682 1.056 0.000 1.00 0.00 O
HETATM 3 O UNL 1 -0.682 -1.056 0.000 1.00 0.00 O1-
HETATM 4 C UNL 1 1.435 0.000 0.000 1.00 0.00 C
HETATM 5 H UNL 1 1.799 0.000 1.028 1.00 0.00 H
HETATM 6 H UNL 1 1.799 -0.890 -0.514 1.00 0.00 H
HETATM 7 H UNL 1 1.799 0.890 -0.514 1.00 0.00 H
CONECT 1 2 2 3 4
CONECT 2 1 1
CONECT 3 1
CONECT 4 1 5 6 7
CONECT 5 4
CONECT 6 4
CONECT 7 4
MASTER 0 0 0 0 0 0 0 0 7 0 7 0
ENDantechamber
对此时的pdb会识别出一些问题(比如类型出错),因此我们手动对pdb进行一些修改,删除一些无用的信息,修改残基名和原子名,修改后的pdb文件内容如下: 1
2
3
4
5
6
7ATOM 1 CA ACT 1 -0.072 0.000 0.000 1.00 0.00 C
ATOM 2 OA1 ACT 1 -0.682 1.056 0.000 1.00 0.00 O
ATOM 3 OA2 ACT 1 -0.682 -1.056 0.000 1.00 0.00 O
ATOM 4 CB ACT 1 1.435 0.000 0.000 1.00 0.00 C
ATOM 5 HB1 ACT 1 1.799 0.000 1.028 1.00 0.00 H
ATOM 6 HB2 ACT 1 1.799 -0.890 -0.514 1.00 0.00 H
ATOM 7 HB3 ACT 1 1.799 0.890 -0.514 1.00 0.00 H
AMBER力场(不使用Gaussian)
对于amber力场,可以使用AmberTools中的antechamber
软件生成小分子的gaff力场。在生成.mol2
文件之前,我们先对.pdb
文件进行预处理。
首先先用tleap转化下原子的名字, tleap.in
文件中的内容如下: 1
2
3
4source leaprc.protein.ff14SB
TMP = loadpdb ACT.pdb
savepdb TMP ACT_tleap.pdb
quitACT_tleap.pdb
: 1
tleap -f tleap.in
reduce
补充H原子(实际上对于小分子即使手动删除貌似也不会补充): 1
reduce ACT_tleap.pdb >ACT_h.pdb
antechamber
软件生成.mol2
文件,命令如下: 1
antechamber -i ACT_h.pdb -fi pdb -o ACT.mol2 -fo mol2 -c bcc -s 2 -nc -1
1. -i ACT_h.pdb
: 指定输入3D结构文件名称。
2. -fi pdb
: 指定输入文件格式为PDB格式。
3. -o ACT.mol2
: 指定输出的文件名称。
4. -fo mol2
: 指定输出文件类型为.mol2
类型。
5. -c bcc
: 指示antechamber使用AM1-BCC电荷模型来计算原子上的电荷。
6. -s 2
: 指示antechamber程序提供的状态信息的冗长度。我们选择提供更多信息(2
)。
7. -nc -1
: 指定净电荷数目为-1
。
ACT.mol2
文件包含了acetate残基的定义,包含了所有的电荷信息以及原子类型。之后会用它来生成.prmtop
和.inpcrd
文件。生成的.mol2
文件如下: 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24@<TRIPOS>MOLECULE
ACT
7 6 1 0 0
SMALL
bcc
@<TRIPOS>ATOM
1 CA -0.0720 0.0000 0.0000 c 1 ACT 0.901600
2 OA1 -0.6820 1.0560 0.0000 o 1 ACT -0.861300
3 OA2 -0.6820 -1.0560 0.0000 o 1 ACT -0.861300
4 CB 1.4350 0.0000 0.0000 c3 1 ACT -0.200100
5 HB1 1.7990 0.0000 1.0280 hc 1 ACT 0.007033
6 HB2 1.7990 -0.8900 -0.5140 hc 1 ACT 0.007033
7 HB3 1.7990 0.8900 -0.5140 hc 1 ACT 0.007033
@<TRIPOS>BOND
1 1 2 1
2 1 3 1
3 1 4 1
4 4 5 1
5 4 6 1
6 4 7 1
@<TRIPOS>SUBSTRUCTURE
1 ACT 1 TEMP 0 **** **** 0 ROOT.mol
文件的第一列是原子序号,第二列为原子名称,第三到五列为三维坐标,第六列为原子类型,最后一列为原子的电荷。并且.mol2
文件还记录了成键信息,但是并不包含对应的参数,GAFF的参数都在$AMBERHOME/dat/leap/parm/gaff.dat
中定义。
参数文件可能会并未包含所需要的参数,因此需要参数进行检查,可以使用parmchk2
检查缺失的参数,命令如下: 1
parmchk2 -i ACT.mol2 -f mol2 -o ACT.frcmod
ACT.frcmod
文件。这是一个参数文件, 能够载入LEaP中用于添加缺失的参数, 这样就能包含所有缺失的参数。在模拟之前,需要仔细检查ACT.frcmod
文件。生成的ACT.frcmod
文件如下: 1
2
3
4
5
6
7
8
9
10
11
12
13Remark line goes here
MASS
BOND
ANGLE
DIHE
IMPROPER
c3-o -c -o 1.1 180.0 2.0 Using general improper torsional angle X- o- c- o, penalty score= 3.0)
NONBON
现在我们已经准备好了acetate作为一个单元的所有工作,只需运行tleap保证GAFF力场可以工作即可。运行以下命令: 1
tleap -f leaprc.protein.ff14SB
之后在tleap中导入GAFF力场: 1
source leaprc.gaff
1
ACT = loadmol2 ACT.mol2
tleap
中键入list
,可以看到新增加的ACT
单元: 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15> ACT = loadmol2 ACT.mol2
Loading Mol2 file: ./ACT.mol2
Reading MOLECULE named ACT
> list
ACE ACT ALA ARG ASH ASN ASP CALA
CARG CASN CASP CCYS CCYX CGLN CGLU CGLY
CHID CHIE CHIP CHIS CHYP CILE CLEU CLYS
CMET CPHE CPRO CSER CTHR CTRP CTYR CVAL
CYM CYS CYX GLH GLN GLU GLY HID
HIE HIP HIS HYP ILE LEU LYN LYS
MET NALA NARG NASN NASP NCYS NCYX NGLN
NGLU NGLY NHE NHID NHIE NHIP NHIS NILE
NLEU NLYS NME NMET NPHE NPRO NSER NTHR
NTRP NTYR NVAL PHE PRO SER THR TRP
TYR VAL frcmod14SBgaff parm10ACT
单元: 1
2
3
4
5
6
7
8
9> check ACT
Checking 'ACT'....
Warning: The unperturbed charge of the unit (-1.000001) is not zero.
Checking parameters for unit 'ACT'.
Checking for bond parameters.
Checking for angle parameters.
check: Warnings: 1
Unit is OK.ACT.frcmod
文件: 1
loadamberparams ACT.frcmod
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31> saveoff ACT ACT.lib
Creating ACT.lib
Building topology.
Building atom parameters.
> saveamberparm ACT ACT.prmtop ACT.inpcrd
Checking Unit.
Warning: The unperturbed charge of the unit (-1.000001) is not zero.
Note: Ignoring the warning from Unit Checking.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
total 1 improper torsion applied
Building H-Bond parameters.
Incorporating Non-Bonded adjustments.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
(Residues lacking connect0/connect1 -
these don't have chain types marked:
res total affected
ACT 1
)
(no restraints)
>ACT.prmtop
和ACT.inpcrd
, 可以用这两个文件生成gromacs需要的ACT.top
和ACT.gro
。这可以用amb2gro_top_gro.py
完成: 1
amb2gro_top_gro.py -p ACT.prmtop -c ACT.inpcrd -t ACT.top -g ACT.gro -b ACT_out.pdb
.top
文件和.gro
文件。
TIPS: 在生成甲铵盐(methylammonium)的力场的时候除了一些识别出错的问题,那是因为我把H原子的名字写成了NH1,在没有最后一列表明原子类型时,它会识别称N原子,所以警告。
1
2
3 Warning: Small distance for BOND 5 N NH1 0 1.01 [1.01 - 1.69]
Warning: Small distance for BOND 6 N NH2 0 1.01 [1.01 - 1.69]
Warning: Small distance for BOND 7 N NH3 0 1.01 [1.01 - 1.69]
AMBER力场(使用Gaussian)
使用gaussian来生成精确的静电分布。这里需要提及, g09存在一些bug导致无法用于RESP, 所以我们使用g16来进行模拟。
gaussian的安装
首先进行gaussian的安装。gaussian的安装方法在该网页中。Linux的gaussian安装包是已经编译好的,因此只要设置好环境变量即可。首先先进行解压,假设压缩包为g16.tar.gz
,则使用以下命令进行解压: 1
tar -xvf g16.tar.gz
g16
文件夹放置在/home/hapo/
文件夹下。先在/home/hapo/g16
文件夹下新建文件夹scratch
。接下来在.bashrc
文件中添加如下内容: 1
2
3export g16root=/home/hapo
export GAUSS_SCRDIR=/home/hapo/g16/scratch
source /home/hapo/g16/bsd/g16.profile/home/hapo/g16
中更改文件的权限: 1
chmod 750 -R *
生成mol2文件
这一步使用别人的脚本,使用方法记录在该Jerkwin博客网页中。
首先生成gaussian程序的输入文件: 1
python pre.py acetate.pdb -1
1
g16 <ACT.gjf >ACT.out
post.py
程序: 1
python post.py
ACT.mol2
文件,其中的静电荷是用gaussian计算后拟合得到的。
最后我们再用以上提及的tleap
和amb2gro_top_gro.py
即可生成gromacs的.top
和.gro
文件。
OPLS-AA力场(使用LigParGen)
LigParGen可以用于生成OPLS-AA力场,只要上传相应的PDB文件进行相应的设置即可。需要注意的是,该网站需要用Edges打开才能正常运行。 点击Submit Molecule
之后,只要下载gromacs的.top
和.gro
文件即可。关于使用LigParGen生成lammaps力场文件的方法在此。
OPLS-AA力场(使用TPPMKTOP)
TPPMKTOP也可以用于生成OPLS-AA力场,并且和文献中的静电值最为接近,该文献可能是使用该程序生成的。 同样,只要上传相应的PDB文件即可。
构建gromacs拓扑文件
为了能够用gromacs进行模拟,先构建gromacs可用的.top
文件。小分子的力场放在molecule
文件夹中。除了以上生成的力场文件外,水分子使用重水进行模拟,水分子的力场放在heavywater
文件夹中。首先先在topol.top
文件中写下如下内容: 1
2
3
4
5
6
7
8
9
10
11#include "amber14sb_parmbsc1.ff/forcefield.itp"
#include "heavywater/tip3p-hw.itp"
#include "molecule/ACT.itp"
#include "molecule/MLM.itp"
#include "amber14sb_parmbsc1.ff/ions.itp"
[ System ]
small molecule in water
[ Molecules ].top
和.gro
的方法。solvate
可以用来给空盒子填充水,可以这样写: 1
gmx solvate -cs spc216.gro -o conf.gro -box 2.5 2.5 2.5 -p topol.top
conf.gro
文件,盒子大小为
同时solvate
程序可以通过-maxsol
来限制添加的溶液数目。因此我们可以用以下命令来生成包含一个醋酸盐和甲铵盐的水盒子: 1
2
3gmx solvate -cs ./molecule/ACT.gro -o cp.gro -maxsol 1 -box 2.5 2.5 2.5 -p topol.top
gmx solvate -cs ./molecule/MLM.gro -cp cp.gro -o cp2.gro -maxsol 1 -p topol.top
gmx solvate -cs spc216.gro -cp cp2.gro -o conf.gro -p topol.top
比较正确的做法应该用insert-molecules
来插入小分子,用solvate
来填充水分子: 1
2gmx insert-molecules -ci molecule/ACT.gro -o box.gro -nmol 1 -box 3.2 3.2 3.2
gmx solvate -cs spc216.gro -cp box.gro -o conf.gro -p topol.top