计算分子动力学直径

高斯计算分子动力学直径

前期准备

主要参考的是这篇博文

主要参考文献为 J. Phys. Chem. A, 118, 1150 (2014)

基于分子电子密度等值面的估计分子动力学直径的方法，对于文中测试的一系列小分子，计算结果和被普遍采用的Breck书里的那些值有不错的线性相关性。

在Multiwfn中输入以下命令

1
2
3
4
5
6
7
12  //定量分子表面分析
1   //选择表面的定义
1   //用电子密度等值面作为表面
0.0015  //用0.0015 a.u.作为等值面数值
6  //开始分析（但不考虑映射的函数，因为这不是我们目前感兴趣的）
6  //将表面顶点坐标导出为当前目录下的vtx.pdb

通过vtx.pdb来考察分子动力学直径，主要有两个方法

  1. VMD，肉眼判断，再测距离

  2. 通过Multiwfn直接统计表面顶点的最大、最小坐标下

软件安装

以下两个软件我都是安装的win版的

Multiwfn的安装

在官网上下载软件，官网链接

下载的是win版的，下载下来解压即可以直接使用。有个小问题是路径中不能有中文，不然无法正确读取文件。

VMD的安装

下载了VMD，准备安装了才发现我以前安装了这个软件，安装这个软件需要登记一下邮箱等信息，如实填写就行。下载链接

复现博文中的CO2分子动力学直径

高斯计算

建模，gjf的编写

按照博文中的数据写的gjf文件，对开头的chk和内存等设置做了修改。

1
2
3
4
5
6
7
8
9
10
11
%chk=./CO2.chk
%mem=50000MB
%nprocshared=32
# PBE1PBE/def2TZVPopt

Title Card Required

0 1
C 0.0 0.0 0.0
O 0.0 0.0 1.2584
O 0.0 0.0 -1.2584

计算与后处理

计算结果出得很快，几秒就计算完毕了, 以下是查看log文件的信息：

1
2
3
4
5
6
7
8
9
10
11
12
[jason@cabc1 new]$ tail CO2.log


 I WANT TO KNOW HOW GOD CREATED THE WORLD.  I AM NOT INTERESTED IN
 THIS OR THAT PHENOMENON, IN THE SPECTRUM OF THIS OR THAT ELEMENT.
 I WANT TO KNOW HIS THOUGHTS, THE REST ARE DETAILS.
     -- ALBERT EINSTEIN
 Job cpu time:       0 days  0 hours  3 minutes 46.0 seconds.
 Elapsed time:       0 days  0 hours  0 minutes  9.3 seconds.
 File lengths (MBytes):  RWF=     10 Int=      0 D2E=      0 Chk=      2 Scr=      1
 Normal termination of Gaussian 16 at Fri Nov 11 17:17:20 2022.

对chk文件用formchk进行转化，将其转化为fchk文件：

1
2
3
4
5
6
7
8
9
[jason@cabc1 new]$ formchk CO2.chk
 Read checkpoint file CO2.chk
 Write formatted file CO2.fchk
 Initial coordinates match /B/.
 Warning:  fchk array  MicOpt element           1 has value 4622945017495814144 set to -1.
 Warning:  fchk array  MicOpt element           2 has value 4625193954506324261 set to -1.
 Warning:  fchk array  MicOpt element           3 has value 4625193954506324261 set to -1.
 Rotating derivatives, DoTrsp=T IDiff= 1 LEDeriv=    485 LFDPrp=       0 LDFDPr=       0.

出现了warning，但是我这个时候没有注意到这个warning，后面了解到这个是因为版本冲突引起的问题，在后文中会给出解决办法，这样依然可以生成fchk文件，并且后续处理我也没有遇到什么问题。

Multiwfn处理计算数据

这里的处理就相当简单了，直接按照前面博文中所写的步骤运行就行，Multiwfn对新手上手来说还是很友好的，尤其是win版，开始直接按enter还可以图形化窗口选择输入文件，唯一需要注意的是文件的路径中不能有中文，不然会读取失败

感觉好多计算的软件都不支持中文路径，MS里面有时候报错也是因为路径中有中文的问题，所以还是在盘符里面单独创建一个完全没有中文的文件来存放计算相关的文件比较合适，计算文件也不要用中文命名，避免出现各种奇葩报错。

VMD查看动力学直径

根据前一步得到的fchk文件，直接拖进VMD里面，然后按照文章中的设置，调整一下，很容易就读出了动力学直径。

原始处理

转化为点图

经过选点处理，得到结果：

1
Added new Bonds label MOL1:C/MOL1:C = 3.346498

与文章中的结果符合得很好，反而跟博文中的结果有点差异，看来选点对结果的影响还是很大的，这一步或许需要用程序来处理了（挖个坑，以后科研有需要了可以想想怎么写），这里跟博文中的结果差异，不知道是选点的问题还是formchk文件的问题，不过总体来说无伤大雅。

对马来酸酐的分子动力学直径的计算

高斯计算

建模，生成gjf文件

马来酸酐的化学式如图：

这个我是用Chem Draw 20.0画的，用这个软件画完就可以直接生成mol文件，然后就可以导入到GaussianView里面生成gjf文件。对开头的chk以及mem和方法等自己输入设置之后，最后得到的gjf文件如下：

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

MA.mol

0 1
 C                 -0.71360000    0.29750000    0.00000000
 C                 -0.71360000   -0.52750000    0.00000000
 C                  0.07110000   -0.78240000    0.00000000
 C                  0.55600000   -0.11500000    0.00000000
 O                  0.07110000    0.55240000    0.00000000
 O                 -1.38100000    0.78240000    0.00000000
 O                  1.38100000   -0.11500000    0.00000000
 H                 -1.57922481   -1.15646239    0.00000000
 H                  0.40178325   -1.80001908    0.00000000

这里用的是另外的基组进行计算，这个是以后我打算在自己的研究上用的泛函与基组，目前来说不清楚这样的设置是否能够得到准确的信息，但练手嘛，暂且先这样用着。

计算与后处理

这个虽然多了几个原子，但是计算结果依然很快，不到一分钟就算完了，以下是tail MA.log得到的信息：

1
2
3
4
5
6
7
8
9
10
11
12
[jason@cabc1 MDD]$ tail MA.log

 ... A MOLECULAR SYSTEM ... (PASSES) ... FROM ONE STATE OF EQUILIBRIUM
 TO ANOTHER ... BY MEANS OF ALL POSSIBLE INTERMEDIATE PATHS,
 BUT THE PATH MOST ECONOMICAL OF ENERGY WILL BE THE MORE OFTEN TRAVELED.

                      -- HENRY EYRING, 1945
 Job cpu time:       0 days  0 hours  8 minutes 24.3 seconds.
 Elapsed time:       0 days  0 hours  0 minutes 16.7 seconds.
 File lengths (MBytes):  RWF=     32 Int=      0 D2E=      0 Chk=      4 Scr=      1
 Normal termination of Gaussian 16 at Fri Nov 11 20:18:16 2022.

随后用formchk转化chk文件的格式：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[jason@cabc1 MDD]$ formchk MA.chk
 Read checkpoint file MA.chk
 Write formatted file MA.fchk
 Coordinates translated and rotated.
 Coordinates match /B/ after translation and rotation.
 Warning:  fchk array  MicOpt element           1 has value 4622945017495814144 set to -1.
 Warning:  fchk array  MicOpt element           2 has value 4622945017495814144 set to -1.
 Warning:  fchk array  MicOpt element           3 has value 4622945017495814144 set to -1.
 Warning:  fchk array  MicOpt element           4 has value 4622945017495814144 set to -1.
 Warning:  fchk array  MicOpt element           5 has value 4625193954506324261 set to -1.
 Warning:  fchk array  MicOpt element           6 has value 4625193954506324261 set to -1.
 Warning:  fchk array  MicOpt element           7 has value 4625193954506324261 set to -1.
 Warning:  fchk array  MicOpt element           8 has value 4607217659633734768 set to -1.
 Warning:  fchk array  MicOpt element           9 has value 4607217659633734768 set to -1.
 Rotating derivatives, DoTrsp=T IDiff= 2 LEDeriv=   1349 LFDPrp=       0 LDFDPr=       0.

这里的warning就很多，我这次就注意到了这个问题，于是我上网查了一下原因，可以参考的文章有：链接1 链接2，大致上就是g06和g16的formchk文件是不一致的，解决办法也很简单，拿对应版本的formchk文件就行，我老师之前给我的应该是g09版本的formchk文件，后面我又重新找他要了g16版本的formchk文件，然后重新配置一下环境变量，运行起来就正常了，以下是正常输出结果的信息：

1
2
3
4
5
6
[jason@cabc1 MDD]$ formchk MA.chk
 Read checkpoint file MA.chk type G16
 Write formatted file MA.fchk
 FChkPn:  Coordinates translated and rotated.
 FChkPn:  Coordinates match /B/ after translation and rotation.

Multiwfn处理

这一步的处理与之前类似，很容易就得到了处理好的vtx.pdb文件

VMD查看分子动力学直径

导入之后的图经过点处理是这样的，还是很难去判断动力学直径的。

大致调一下视角，选点并测距离。

1
Info) Added new Bonds label MOL1:C/MOL1:C = 3.841880

最后得出MA的分子动力学直径为3.8419埃

总结

一些思考

关于计算过程

这个获得的分子动力学直径正如博文中所说，还是比较粗糙的，并没有考虑极化等对电子密度的影响，感觉这个对分子动力学的影响还挺大的，还有溶剂化结构什么的，这一方面的计算还有提升的空间，可能需要更加科学合理的定义以及开创新的计算方法才行。

关于取点

在取点的过程中真的是比较折磨的，对其，取点，真的是一个只可意会不可言传的东西，主观性还是有的，虽然最后取出的值不会有太大的差异，但是终究来说还是会让人比较抓狂，自我怀疑，这一部分的工作我觉得可以交给程序去实现，但是具体的实现步骤就还得好好学习一下编程了。任重而道远呀。

关于报错

formchk无法运行，没有这个command

1
formchk: command not found

参考链接1

问题点在于软件的安装配置问题，这个需要找管理员解决，或者知道高斯的安装目录，自行设置环境变量。

  • 解决办法

    找老师拿了一个formchk的文件，自己配置了一下环境变量

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # 创建了一个目录 software/formchk_bin,将老师给的formchk文件置于其中
    # 给formchk添加运行权限
    chmod +x formchk

    # 配置环境变量
    vim ~/.bashrc
    # 在其中添加
    export PATH="/home/jason/software/formchk_bin:$PATH"

    #保存退出，刷新shell的环境
    source ~/.barshrc
#Study #Gaussian #MD #Chemistry
计算分子动力学直径
http://phoenixjason.cn/2022/11/12/20221112高斯计算分子动力学直径/
作者
Jason
发布于
2022年11月12日
许可协议