新冠阳了，今天开始复建运动。先试着把小分子的力场文件构建起来吧。这是我并不熟悉的领域，所以需要学习好多教程。

获得PDB结构

这里我们主要是要制作醋酸盐(acetate)和甲铵盐(methylammonium)的力场，我们首先需要得到小分子的PDB结构，可以从PDB库中找到我们需要的小分子。醋酸盐和甲铵盐的代码分别为ACT和3P8。在网页上可以下载到.sdf文件，接下来可以使用openbabel将.sdf文件转化为.pdb坐标。

这样就获得了对应的.pdb文件。

此时由openbabel生成的acetate的.pdb文件如下:

COMPND    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
END

antechamber对此时的pdb会识别出一些问题(比如类型出错)，因此我们手动对pdb进行一些修改，删除一些无用的信息，修改残基名和原子名，修改后的pdb文件内容如下:

ATOM      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文件中的内容如下:

source leaprc.protein.ff14SB
TMP = loadpdb ACT.pdb
savepdb TMP ACT_tleap.pdb
quit

用以下的命令生成ACT_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

上面各个选项的意思分别为:

-i ACT_h.pdb: 指定输入3D结构文件名称。
-fi pdb: 指定输入文件格式为PDB格式。
-o ACT.mol2: 指定输出的文件名称。
-fo mol2: 指定输出文件类型为.mol2类型。
-c bcc: 指示antechamber使用AM1-BCC电荷模型来计算原子上的电荷。
-s 2: 指示antechamber程序提供的状态信息的冗长度。我们选择提供更多信息(2)。
-nc -1: 指定净电荷数目为-1。

ACT.mol2文件包含了acetate残基的定义，包含了所有的电荷信息以及原子类型。之后会用它来生成.prmtop和.inpcrd文件。生成的.mol2文件如下:

@<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文件如下:

Remark 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

可以看到缺失了一个反常二面角，假定antechamber建议的参数可以接受，无需修改。

现在我们已经准备好了acetate作为一个单元的所有工作，只需运行tleap保证GAFF力场可以工作即可。运行以下命令:

1	tleap -f leaprc.protein.ff14SB

之后在tleap中导入GAFF力场:

1	source leaprc.gaff

现在导入acetate单元:

1	ACT = loadmol2 ACT.mol2

此时在tleap中键入list，可以看到新增加的ACT单元:

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

检查ACT单元:

> 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

接下来就可以将生成的残基导出了:

> 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

这样我们就有了gromacs可用的.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
3

export g16root=/home/hapo
export GAUSS_SCRDIR=/home/hapo/g16/scratch
source /home/hapo/g16/bsd/g16.profile

接下来在进入/home/hapo/g16中更改文件的权限:

1	chmod 750 -R *

加下来只要重开终端即可运行g16。

生成mol2文件

这一步使用别人的脚本，使用方法记录在该Jerkwin博客网页中。

首先生成gaussian程序的输入文件:

1	python pre.py acetate.pdb -1

之后运行gaussian程序:

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文件中写下如下内容:

#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文件，盒子大小为$2.5nm \times 2.5 nm \times 25 nm$，其中填充了水分子。

同时solvate程序可以通过-maxsol来限制添加的溶液数目。因此我们可以用以下命令来生成包含一个醋酸盐和甲铵盐的水盒子:

1
2
3

gmx 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 2	gmx 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

Hapo's blog

制作小分子力场