计算化学公社

 找回密码 Forget password
 注册 Register
Views: 1410|回复 Reply: 2

[建模与可视化] 利用ABCluster和Multiwfn自动生成有机分子拓扑文件用于全局优化

[复制链接 Copy URL]

220

帖子

8

威望

3030

eV
积分
3410

Level 5 (御坂)

发表于 Post on 2022-5-27 11:28:25 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 coolrainbow 于 2022-5-27 11:29 编辑

ABCluster是一个广泛使用的团簇结构优化软件(ABCluster - Download (zhjun-sci.com)),对原子团簇、分子团簇、表面支撑团簇到周期性团簇都可以取得非常好的效果(案例:ABCluster - Stories (zhjun-sci.com) 图片:ABCluster - Gallery (zhjun-sci.com))。其中的 rigidmol 是一个非常好用的模块,可以对有机分子非常快速的预测团簇结构(例如ABCluster - Stories (zhjun-sci.com) 中的Atmospheric Chemistry案例中的大部分都是用rigidmol做的)。
对于rigidmol,预测结构本身非常快,但是对于大分子,拓扑文件的构造非常麻烦,在过去的版本中需要手动指认原子类型等。虽然有时可以用CHARMMGen来自动生成参数和原子类型,但是CHARMMGen是在线程序,使用不便,而且一些分子常常预测失败,即使成功了,生成的文件也需要经过一些转换才能被ABCluster使用。

从ABCluster 3.1开始,程序自带的组件topgen的功能得到了极大的增强,已经可以自动产生任意有机分子的拓扑文件。结合Multiwfn方便的电荷拟合功能,可以直接用于ABCluster的全局优化。本文就结合一个例子来展示这些功能。

1 产生有机分子的结构

本文考虑的有机分子是下面这个“钳子”:
202205270822355241..png
首先使用GaussView或者任意软件画出这个分子,然后做几何优化,保存波函数文件。如果是Gaussian,可以利用如下输入文件:

%nprocs=24
%mem=20GB
%chk=pincer.chk
#b3lyp/6-31g(d) geom=connectivity

tile

0 1
[ coordinates ]

优化完后,将chk转化为fchk文件:
formchk pincer.chk pincer.fchk

再把.out文件中优化的稳定结构结构转化为pincer.xyz (见附件)。

至此,你已经有了pincer的结构和波函数信息。

2 产生有机分子的电荷

利用Multiwfn计算RESP电荷(详见RESP拟合静电势电荷的原理以及在Multiwfn中的计算 - 思想家公社的门口:量子化学·分子模拟·二次元 (sobereva.com))。执行如下命令:
Multiwfn pincer.fchk
7
18
1
y

这样电荷就存储到了pincer.chg文件中。

注意,对于分子的定量/半定量模拟,RESP电荷或者某类其它电荷都可以尝试,但是类似Mulliken电荷或者EEQ这种绝对不能使用!
对于大分子的电荷拟合,建议读者引用静电势计算的算法文献:Phys. Chem. Chem. Phys. 2021, 23, 20323-20328 (详见Multiwfn使用的高效的静电势算法的介绍文章已于PCCP期刊发表! - 思想家公社的门口:量子化学·分子模拟·二次元 (sobereva.com)

3 产生有机分子的拓扑

利用ABCluster中的topgen可以一键生成拓扑:
topgen pincers.xyz

输出如下:
Analyze topology for the molecule from "pincer.xyz" ... done.

Exporting results to "pincer" +
.gjf:         GJF file with bonding information.
-cycles.txt:  Containing cycle texts that can be used for ABCluster:geom calculation.
-bonding.xyz: XYZ file with bonding information for ABCluster:geom calculation.
-rigid.xyz:   XYZ file with bonding information for ABCluster:rigidmol calculation.
.psf:         PSF file with topology informaition for ABCluster:geom/NAMD calculations.

In "pincer-rigid.xyz" and "pincer.psf":
"X" means unknown atomic type.
Total charge: 1.8360
-----------------------------------------------------------
| You can adjust charges to meet the target total charge, |
| or re-fit RESP charges using, e.g., Multiwfn.           |
-----------------------------------------------------------

Atoms are typed using graph representation learning. Please cite:
Zhang, J. Atom Typing Using Graph Representation Learning: How Do Models Learn Chemistry?
J. Chem. Phys. 2022, 156, 204108

输出包含5个文件:
  • pincer.gjf: Gaussian输入文件,可以用于查看生成的键连情况。
  • pincer-cycles.txt:分子成环情况。比如“1 21 26 27 28 29 30 31 33 ”就是由这几个原子生成的九元环。
  • pincer-bonding.xyz:带有成键信息的xyz,用于ABCluster中geom的全局优化、构象搜索。详见7.1. Flexibility — ABCluster 3.1 documentation (zhjun-sci.com)
  • pincer-rigid.xyz:带有参数信息的XYZ,用于ABCluster中rigidmol的全局优化。见下。
  • pincer.psf:拓扑信息。可以用于NAMD或者CHARMM的分子动力学模拟。

topgen在产生-rigid.xyz和.psf文件时,需要猜测原子的类型。这个类型推断用图表示学习的算法完成,请引用J. Chem. Phys. 2022, 156, 204108

对于科研级别的有机分子的自动类型推断,通常都需要人工再检查一遍,甚至重新参数化。无论是用topgen还是其它自动类型推导方法(如CGenFF或者GAFF)。在这里就直接用了。

现在,把Multiwfn计算的电荷放到pincer-rigid.xyz和pincer.psf中,也即,把pincer.chg的最后一列
O    -1.498911   -3.700847    1.522344  -0.1977931830

替换到pincer-rigid.xyz的下面第一列
par_all36_prot   q   epsilon (kJ/mol) sigma (AA)
-0.1977931830  0.4184   2.9400 # OG3C51

同时替换到pincer.psf的下面第七列
         1 U         1        PINC     O        OG3C51  -0.1977931830       15.9990           0

这样就完成了相关的拓扑文件。

4 全局优化实例

利用rigidmol对分子做全局优化的细节可以参考4. rigidmol: Rigid Molecular Clusters Using Force Fields — ABCluster 3.1 documentation (zhjun-sci.com).这里只罗列一下步骤。这里想考虑一个pincer和13个水的全局极小点。水就用tip4p.xyz,pincer就用pincers-rigid.xyz。

首先构造pincer-wat.cluster,团簇成分文件:
2
pincer-rigid.xyz 1
tip4p.xyz        13
* 20.0000

然后是pincer-wat.inp,优化文件:
pincer-wat.cluster  # cluster file name
100                 # population size
100                 # maximal generations
3                   # scout limit
10.0                # amplitude
pincer-wat          # save optimized configuration to .xyz and .gjf
30                  # number of LMs to be saved

这里,用一个100的种群迭代100步,最终存储30个isomer。现在可以运行:
rigimol pincer.inp > pincer.out

这样,在pincer-wat-LM下的0.xyz就是全局极小点,如下所示:

202205270946184367..png

可见,水分子分布在钳子之间,并且形成了一系列氢键网络。

顺便,你还可以预测“钳子”二聚体的结构,如下所示:

202205270949332615..png
所有相关的文件可以在附件pincers.zip中下载。



pincer.zip

21.99 KB, 下载次数 Times of downloads: 8

评分 Rate

参与人数
Participants 6
威望 +1 eV +22 收起 理由
Reason
zsu007 + 5 好物!
Acee + 2 谢谢分享
sobereva + 1
joeson + 5 好物!
chands + 5 好物!
妙角不脆 + 5 好物!

查看全部评分 View all ratings

1

帖子

0

威望

841

eV
积分
842

Level 4 (黑子)

发表于 Post on 2022-5-27 13:14:44 | 显示全部楼层 Show all
好物,点赞!

14

帖子

1

威望

847

eV
积分
881

Level 4 (黑子)

发表于 Post on 2022-5-29 17:45:51 | 显示全部楼层 Show all
非常好。 一个问题是geom可以预测单个柔性分子的构象,如果两个或者多个这样的分子放在一起,可以预测其结合模式吗?就是结合时让每一个分子的构象都自由变化

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-7 03:24 , Processed in 0.646551 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list