乙腈分子力场参数的获取

乙腈分子力场参数的获取

最近计算一个体系要用到乙腈的模型，这里就想顺手在生成乙腈分子拓扑文件的时候记录一下，毕竟教程肯定不会嫌多的。

建模并优化—Gaussian

通过Gaussian先对建模的乙腈（AN）分子进行结构优化，并获取它的fchk文件，一下是建立的模型：

导出gjf格式的文件，然后修改一下gjf的文件内容，一定要记得加上freq关键词：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
%chk=./AN.chk
%mem=10000MB
%nprocshared=32
# 6-311+G(d,p) em=gd3 b3lyp opt freq

AN

0 1
 C                 -1.66216219    0.04054044    0.00000002
 H                 -1.30550794   -0.96826963   -0.00000005
 H                 -1.30548929    0.54493894   -0.87365128
 H                 -2.73216219    0.04055361   -0.00000001
 C                 -1.14881994    0.76649699    1.25740481
 N                 -0.76661338    1.30700424    2.19360010



这个分子还是相对大一点，所以我是上传到服务器集群上计算的，优化完成之后用formchk转化一下格式，将chk格式的文件转化为fchk文件，在工作目录下输入：

1
formchk AN.chk

注意，用formchk时需要注意你的环境变量里面是否添加了这个脚本，没有的话记得联系你们的管理员。

得到了fchk文件之后我们就可以开始下一步操作了。

Multiwfn生成电荷信息（RESP）

对于乙腈这样的分子结构来说，极性比较强，而且在我计算的体系中它是作为溶剂，所以获得每个原子的电荷再来编写它的拓扑文件是非常有必要的，Multiwfn提供了非常便捷的操作来生成分子的RESP电荷。

安装好Multiwfn之后，直接再命令行敲Multiwfn，然后输入./AN.fchk，将刚刚得到的fchk文件输入，再进行后续操作，一下是输入的内容：

1
2
3
4
5
6
Multiwfn
./AN.fchk #输入文件
7
18
1
y

后续的操作在Gromacs的学习以及应用（五） - PhoenixJason这篇博客中有解释，可以去看看。

这样我们会获得一个AN.chg文件，cat一下这个文件看一下：

1
2
3
4
5
6
7
[jason@cabc1 01Gaussian]$ cat AN.chg
C    -0.000000    0.000000   -1.176511  -0.4202575143
H     0.000000    1.025236   -1.552823   0.1595174348
H    -0.887881   -0.512618   -1.552823   0.1595174348
H     0.887881   -0.512618   -1.552823   0.1595174348
C     0.000000   -0.000000    0.280675   0.4424515331
N    -0.000000   -0.000000    1.433355  -0.5007463231

这里面保存了AN每个原子的元素、坐标信息和所带电荷。

这样我们目录下就会有.fchk，.chg，.log文件，为了方便操作，接下来我们可以将这三个文件全部输出到自己的windows系统上进行下一步操作。

Sobtop生成带电荷信息的top、itp文件

拿到.log文件之后首先在GaussView中看一下计算的收敛情况，看有没有虚频之类的，一般这种小分子都不会出现什么大问题。然后就可以输出一份mol2文件作为Sobtop的输入文件。

打开Win版的Sobtop，以刚刚输出的mol2文件作为结构文件输入(最简单的办法就是直接将那个文件拖到sobtop的命令行页面中)，然后先选择7，为每个原子分配电荷，再选择10，用chg文件来分配，然后输入刚刚我们得到的chg文件的路径（也可以直接拖入），

成功分配之后如下图：

接下来就可以进行下一步，生成top文件了，再输入0，return到开始的界面，然后输入1生成top文件，这里我们用GAFF力场，输入3，所有原子都被指定了类型，那么就输入0，进入下一步，设定成键信息，因为我们有AN分子的fchk文件，所以这里我们可以直接选2，这是Sobtop的一种根据Gaussian计算结果获取成键参数的方法，然后拖入刚刚得到的fchk文件，最后就是连续两下Enter键，将得到的top文件和itp文件输出。

这里遇到一个蜜汁bug，之前还可以正常输出top和itp文件的，我这想截个图就会在输出文件的那一步直接关闭窗口，而且没有top文件的输出，很奇怪，试了好几遍都不行。

后面测试了一下，初步估计是sobtop目录下已经存在了AN.top和AN.itp文件导致的，而且不知道为啥，我尝试几遍之后fchk文件会变得不能用，原来300多kb的文件会变得只剩下1kb的大小，需要重新从服务器上下载，很奇怪，但是最后还是弄好了，一下是完成之后的截图：

这样我们即可以成功的获得AN分子的GAFF力场的拓扑文件了，接下来我把这俩文件的内容贴出来吧，需要的可以自取：

AN.top:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
; Created by Sobtop (http://sobereva.com/soft/sobtop) Version 1.0(dev3.1) on 2023-04-18

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ    fudgeQQ
     1              2              yes            0.5       0.8333

#include "AN.itp"

[ system ]
AN

[ molecules ]
; Molecule      nmols
AN     1

AN.itp:

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
32
33
34
35
36
37
38
39
40
41
42
43
; Created by Sobtop (http://sobereva.com/soft/sobtop) Version 1.0(dev3.1) on 2023-04-18

[ atomtypes ]
; name   at.num      mass       charge   ptype     sigma (nm)    epsilon (kJ/mol)
c3           6    12.010736    0.000000    A      3.399670E-01    4.577296E-01
hc           1     1.007941    0.000000    A      2.649533E-01    6.568880E-02
c1           6    12.010736    0.000000    A      3.399670E-01    8.786400E-01
n1           7    14.006703    0.000000    A      3.249999E-01    7.112800E-01

[ moleculetype ]
; name          nrexcl
AN       3

[ atoms ]
;  Index   type   residue  resname   atom        cgnr     charge       mass
     1     c3         1      MOL     C1            1   -0.42025751   12.010736
     2     hc         1      MOL     H2            2    0.15951743    1.007941
     3     hc         1      MOL     H3            3    0.15951743    1.007941
     4     hc         1      MOL     H4            4    0.15951743    1.007941
     5     c1         1      MOL     C5            5    0.44245153   12.010736
     6     n1         1      MOL     N6            6   -0.50074632   14.006703
 
[ bonds ]
; atom_i  atom_j  functype      r0 (nm)    k (kJ/mol/nm^2)
    1       2         1        0.109212     2.933368E+05     ; C1-H2, mSeminario method
    1       3         1        0.109212     2.933369E+05     ; C1-H3, mSeminario method
    1       4         1        0.109212     2.933369E+05     ; C1-H4, mSeminario method
    1       5         1        0.145719     2.582536E+05     ; C1-C5, mSeminario method
    5       6         1        0.115268     1.169430E+06     ; C5-N6, mSeminario method
 
[ angles ]
; atom_i  atom_j  atom_k  functype    a0 (Deg.)  k (kJ/mol/rad^2)
    2       1       3         1       108.778      2.808113E+02     ; H2-C1-H3, mSeminario method
    2       1       4         1       108.778      2.808113E+02     ; H2-C1-H4, mSeminario method
    2       1       5         1       110.156      3.116111E+02     ; H2-C1-C5, mSeminario method
    3       1       4         1       108.778      2.808140E+02     ; H3-C1-H4, mSeminario method
    3       1       5         1       110.156      3.362116E+02     ; H3-C1-C5, mSeminario method
    4       1       5         1       110.156      3.362118E+02     ; H4-C1-C5, mSeminario method
    1       5       6         1       180.000      1.662352E+03     ; C1-C5-N6, mSeminario method

; No pair needs to generate


总结

总的来收，整个力场文件的获取就是Gaussian —> Multiwfn —> Sobtop，后面两个软件Sobereva都有非常详细的教程，具体可以参考他的教程，这样应该是能够解决大部分小分子气体的力场参数的。

#Study #Computer #Gromacs
乙腈分子力场参数的获取
http://phoenixjason.cn/2023/04/18/20230418乙腈分子力场参数的获取/
作者
Jason
发布于
2023年4月18日
许可协议