计算化学公社

 找回密码 Forget password
 注册 Register
Views: 6603|回复 Reply: 19

[辅助/分析程序] 计算化学键强度的利器 LModeA-nano 正式版发布

  [复制链接 Copy URL]

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

发表于 Post on 2022-3-29 10:00:29 | 显示全部楼层 Show all |阅读模式 Reading model
计算化学键强度的利器 LModeA-nano 正式版发布


文/smutao, beefly

2022年3月28日





化学键可以称得上是化学研究中最重要的概念之一,其物理本质是原子之间的作用力。在量子化学现有的理论框架下,可以用于分析化学键的手段可谓五花八门。从定态薛定谔方程的角度出发,我们可以将现有的化学键分析方法分为两大类:

一是从波函数信息中提取挖掘尽可能多的有用信息、帮助我们理解各类化学键本质。这一类的方法包含了着眼于电子密度的AIM分析、各类定域化轨道模型、各种键级、电子布居分析和电荷模型、键的极性等。这些方法在Multiwfn程序中可以计算,具体可参见卢天博士的《Multiwfn支持的分析化学键的方法一览》一文。




202203290948103918..png


二是从能量出发着重于衡量化学键的强度。这一类的方法包含了键解离能(bond dissociation energy, 即断开化学结构中单个化学键所需的能量)、键能(bond energy, 即化学结构中某类多个化学键断开时所需要的平均能量)、结合能(binding energy, 即由一个非键相互作用稳定的二体结构分开时所需要的能量)。这三种与键强有关的能量既可以在实验上测量,当然也可以通过量子化学计算获得。由这些可观测的能量衍生出来了能量分解的计算方法,其将一个分子体系划分为了A、B两个片段,它们之间的相互作用能量分解成了更易理解的物理作用能(如静电能、Pauli互斥能、轨道作用能等),如此一来我们就可以知道A、B之间的化学键是哪种物理作用占主导。目前能量分解有着不同的做法,具体可以参见 Frenking (WIREs Comput. Mol. Sci. 8, e1345, 2018)、Head-Gordon (Annu. Rev. Phys. Chem. 72, 641, 2021)、Patkowski (WIREs Comput. Mol. Sci. 10, e1452, 2020) 等人的综述文章。

然而,这些以能量来衡量化学键强度的方法一大痛点在于必须将整个体系人为地划分为A、B两个片段,而这种划分一方面可能不唯一(比如某个键可以是均裂也可以是异裂),另一方面对于比较复杂的体系这种划分没有办法进行(例如氩的唯一已知化合物HArF,不同的划分会得到完全不同的结论)。因此,比较理想的衡量化学键强度的方法应该满足这些条件:1)不需要划分片段,2)基于单个化学结构、不改变结构本身的状态,3)普遍通用。如果还有第4个,那就是可以与实验测量挂钩、有确实的物理意义。

可以满足上述4个条件的一种化学键强度度量是键的力常数,其一般定义为能量对某个键长度的二阶偏导。 力常数之所以被用来表征键强,一方面是因为它的物理含义是势能面在平衡结构处在该键伸缩方向的曲率,而这个描述势能面陡峭或平缓的量和键解离能有大致的正相关关系(即BDE大的键,力常数往往也越大);另一方面是因为它在振动光谱中被广泛应用,虽然振动光谱给出的是离域振动模式的频率信息,但是对于一些局域性较强的键(如C-H,C=O),可以很方便地给出该键对应的振动频率和力常数。

使用力常数来衡量键强的理论计算方法很多。比较早期的做法是将用来拟合振动光谱的力场参数直接作为键强度来使用。等到量子化学发展起来的时候,能量二阶导Hessian矩阵的计算成为可能。有的做法是将笛卡尔坐标系下的Hessian矩阵中涉及两个成键原子的非对角块做投影;有的是在包含了某个键的(非)冗余内坐标下构建Hessian矩阵,并提取该化学键对应的对角元;还有的是在平衡结构上,对某个键进行微小幅度的伸缩,获得能量信息后套用数值求导公式来计算简谐近似的力常数……       

在诸多基于力常数的方法中,比较系统的理论当属由 Cremer 和 Konkoli等人于1998年开创的局域振动模式(Local Vibrational Mode)理论。其基本推导思路为:在一个分子结构中选取一个感兴趣的内坐标(比如一个化学键的键长),将和这个内坐标无关的所有其他原子质量设置为0,再求解对应的质量去耦合的欧拉-拉格朗日方程后,可以得到一种由该内坐标主导的振动模式——局域振动模式。和我们熟知的简正模式(normal vibrational mode)相比,局域振动的优势在于可以研究分子结构中任意一个部分内坐标的振动特征,而简正模式则是离域的。我们恰好可以利用这一点对分子结构中任意一对原子之间的化学键,从振动的角度去研究。

和简正模式一样,局域振动有自己的振动模式向量、振动频率、力常数和约化质量。在实际的理论计算中,局域振动可以从已知的简正模式中直接提取,也可以从简正模式所需的Hessian矩阵信息开始计算获得(以常用的量化程序高斯为例,我们只需要做一个普通的opt+freq计算,便可以往下进行局域振动分析)。在局域振动模式理论提出以后,Cremer和Kraka等人便开始将局域振动模式理论应用于烷烃分子体系的C-H化学键中,并发现计算得到的C-H键的局域伸缩振动频率与实验测量得到的第5泛音带振动频率和同位素取代烷烃的C-H振动频率线性相关。同时,他们还利用局域伸缩振动力常数(local stretching force constant)将最开始只适用于双原子分子的Badger规则拓展至多原子分子。

1963年,Decius提出将compliance 常数的倒数——柔性力常数 (relaxed force constant) 用于描述化学键的强度的方法,后来被Grunenberg等人发展为计算程序。此基于compliance矩阵的方法没有给出完整的推导路径,经常被批评为缺少物理意义。2012年,Zou等人证明了 Grunenberg的柔性力常数和局域振动模式理论中的局域模式力常数(local mode force constant)等价,并开发了 adiabatic connection scheme(ACS)方法。ACS方法可以将一个分子的所有简正模式通过一个耦合因子λ (0~1)平滑地转化为相同数目的已事先定义的非冗余内坐标的局域振动模式。

2019年,Tao等人利用固体的声子振动在Γ点可以由红外、拉曼方法测量这一特性,将局域振动模式理论从分子体系(0维)拓展到了周期性体系(1-3维——涵盖了聚合物、表面、晶体/固体),极大拓展了该理论测量化学键强度的应用领域。值得一提的是:在该理论框架下,同种化学键的强度可以在0-3维任意维度下的体系之间进行直接的计算比较。因此,这一方法具有极高的应用价值。在不久之前,这一项目获得了美国国家科学基金会(NSF)的专项资助。

最近,Tao等人开发了一套用来进行局域振动模式分析的PyMOL插件代码——LModeA-nano,该程序既可以用来分析0维的分子体系,又可以分析1-3维的周期性体系。目前该程序不仅支持 VASP、CP2K、CRYSTAL、CASTEP和Quantum ESPRESSO这些第一性原理计算程序,还原生支持Gaussian和Q-Chem量子化学计算程序,并通过UniMoVib (https://github.com/zorkzou/UniMoVib)支持ORCA、XTB、CFOUR、BDF、PSI4等三十余款量子化学计算程序。介绍该程序的文章已发表(Y. Tao, W. Zou, S. Nanayakkara, et al., LModeA-nano: A PyMOL Plugin for Calculating Bond Strength in Solids, Surfaces, and Molecules via Local Vibrational Mode Analysis, J. Chem. Theory Comput. 18, 1821 (2022) doi:10.1021/acs.jctc.1c01269).  LModeA-nano 的代码可以在GitHub (https://github.com/smutao/LModeA-nano)以及 GitLab (国内镜像 https://gitlab.com/smutao/LModeA-nano)下载获得。有关该代码的安装、使用方法,可以参见专门的手册站点 (https://lmodea-nano.readthedocs.io/).

使用LModeA-nano 进行局域振动模式分析的大致流程和常见的振动分析计算一样,都需要已经优化到势能面极小点的分子或固体原胞结构以及对应的Hessian矩阵(见下图)。与普通的振动分析计算不同之处在于,使用LModeA-nano进行局域振动模式分析需要用户通过选取结构中的原子以指定需要分析的内坐标,这个选取过程可以很方便地通过PyMOL的界面完成。经过计算后,LModeA-nano 会给出键长和键角对应的局域振动模式的力常数和频率信息,键长对应力常数的单位为 mdyn/Angstrom而键角力常数的单位为 mdyn*Angstrom/(rad^2), 频率的单位为波数 (cm^-1).

202203290950375897..png


关于本文中提及的局域振动模式理论的详细信息和应用案例,可以参见Kraka等人的综述文章(WIREs Comput. Mol. Sci. 10, e1480, 2020 doi:10.1002/wcms.1480).



关于LModeA-nano使用以及局域振动模式理论的常见问题Q&A

1. 我以前没有接触过该方法来计算化学键的强度,这个方法靠谱吗?在文章中使用这个方法会不会被审稿人批评?

答:局域振动模式理论的力常数是定量描述键强(尤其是相同类型的化学键)的一个可靠方法。经过多年的测试和探索,该方法已经被用来研究众多类型的共价键和非键相互作用(氢键、卤键等),并为化学键研究领域作出了诸多贡献。下表是综述文章中该方法的应用一览(截至到2020年4月)。



202203290953352013..png


值得一提的是,在能量分解计算方法领域颇有建树的Gernot Frenking教授等人近期在 Int. J. Quantum Chem. 期刊撰文发表了题为 “The Strength of A Chemical Bond” 的文章 (IJQC 122, 26773, 2022 doi:10.1002/qua.26773) ,其指出基于局域振动模式理论的力常数是目前现有方法中可以应用范围最广、最为通用的化学键强度计算方法。
英文原话 (1):   … the force constant is the most general measure for determining the strength of a chemical bond in molecules; (2):  … the best broadly applicable method for estimating the strength of chemical bonds that is presently available.


2. 程序手册和很多使用了该方法的文章里都要求计算的分子或固体结构应处于势能面的极小点才能被分析,这个要求是硬性规定吗?非稳定点(梯度不为0)或者过渡态结构行不行?

答:这个是硬性规定,同时也是应用局域振动模式理论的一个局限性。从振动光谱的角度也不难理解,仪器测得的是一个平均的、稳定结构的振动。 非稳定点和过渡态结构的寿命都很短,一般的振动光谱测不到。
但是也有例外,假如有一系列柔性扫描路径下的结构,可以先用广义子系统振动分析(GSVA)方法提取柔性扫描中发生驰豫的部分的Hessian矩阵,然后用这个Hessian对发生驰豫的片段做局域振动模式分析。

至于对过渡态结构的局域模式振动分析,目前还在开发测试中,文章尚未发表。

3. 我随意选取了结构中的两个原子(如水分子的两个氢原子),在LModeA-nano程序中计算得到了一个力常数,这个结果可以作为这两个原子之间化学键的键强吗?

答:不可以。两个原子之间有一个力常数的值不能推出两个原子之间存在化学键或非共价相互作用,因此使用力常数来计算键强弱之前,必须了解这个是不是已经确定是众所周知的化学键。对于还没有文献报道的化学键,最好先通过波函数分析的诸多方法确认。

一个反面的例子是有人声称用柔性力常数(等价于局域振动模式的力常数)错误地得到了HCN分子中本不存在的强H-N键 (WIREs Comput. Mol. Sci. 4, 111, 2014 doi:10.1002/wcms.1155),这是不了解力常数计算流程所致。

4. 已经有了Grunenberg的Compliance 程序,为什么还要用这个LModeA-nano程序呢?

答:Compliance程序给出的柔性力常数虽然在理论上是等价于局域振动模式理论的力常数的,但LModeA-nano有以下优势:

- 支持 分子和周期性体系的分析 (支持5款主流的第一性原理程序和30余款量子化学计算程序)
- 除了力常数之外,还可计算相应的局域模式振动频率
- 依托PyMOL程序支持三大操作系统
- 代码、手册和算例持续更新

而 Compliance程序 只能分析高斯程序计算的分子体系、没法计算振动频率、只支持Ubuntu Linux系统。

5. LModeA-nano 可以同时给出力常数和频率,我该用哪个作为力常数的度量呢?

答:多数情况下应使用力常数,因为力常数是不依赖原子质量的物理量。振动频率在比较同一种化学键的情况下是有用的,因为同种化学键的成键原子质量相同。

6. 除了键长、键角,LModeA-nano可以分析其他内坐标的局域振动模式吗?

答:目前的LModeA-nano代码只支持这两种,但理论上可以支持多种内坐标包括:锥化角 (pyramidalization angle)、面外角 (out-of-plane angle)、Cremer-Pople ring puckering、ring deformation等。在以后的LModeA-nano版本中可能会支持。

7. 除了量子化学计算方法描述的力场之外,能否对神经网络学习势(neural network potential)力场描述下的化学结构使用LModeA-nano进行局域振动模式分析?

答:这是完全可以的,只要相应的NNP可以提供体系的Hessian矩阵信息(目前的NNP多可以通过神经网络引擎的自动求导功能计算Hessian)。

一个比较方便LModeA-nano分析的做法是:将NNP的计算引擎(如TorchANI)通过External关键词接入高斯程序,然后进行opt+freq计算得到fchk文件,给到LModeA-nano程序后在程序列表中选择Gaussian/Q-Chem即可进行分析。

8. 目前局域振动模式理论方法已经可以分析气态的分子体系和固体,以后可以用这个方法来分析液体吗?

答:对于有溶剂分子包裹的团簇分子,按一般的分子体系分析即可。

涉及显式溶剂以及动力学的方向,相关方法开发还在进行中。

9. 局域振动模式的频率可以直接由实验测得吗?

答:除了双原子分子、同位素取代、泛频等特殊情况外,一般无法在实验上测量,因为振动光谱测量的是简正模式的频率。但是可以把简正模式理解为混合的局域振动模式,因此局域振动模式的频率还是有价值的。

10. LModeA-nano和LModeA程序是什么关系?两者有何不同?

答:Kraka等人的综述文章(WIREs Comput. Mol. Sci. 10, e1480, 2020)提及的LModeA程序是由Zou等人开发的主要适用于分子体系的局域振动模式分析代码。相比LModeA-nano,LModeA支持额外的功能,包括(1)更多种的内坐标定义、(2)文中提到的ACS方法、(3)简正模式分解为局域模式等。目前LModeA代码需要申请获取,将在以后开源。



如有更多关于LModeA-nano的问题,欢迎跟帖交流。也可扫码加入微信群,二维码在下方置顶贴中。

最后希望LModeA-nano对各位的科研工作有所助益!

评分 Rate

参与人数
Participants 15
威望 +1 eV +70 收起 理由
Reason
qchem + 5
Marc0 + 5 精品内容
ghifi37 + 5 精品内容
邓苏微 + 5 牛!
hebrewsnabla + 5 GJ!
cokoy + 5 好物!
zsu007 + 5 赞!
ggdh + 5 猛!
丁越 + 5 赞!
卡开发发 + 5 谢谢分享
hdhxx123 + 5 GJ!
biogon + 5
yygong + 5 牛!
Novice + 5 赞!
sobereva + 1

查看全部评分 View all ratings

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-29 11:03:35 | 显示全部楼层 Show all
本帖最后由 smutao 于 2022-3-30 23:26 编辑

LModeA-nano使用微信群二维码如下

261

帖子

0

威望

2073

eV
积分
2334

Level 5 (御坂)

计算化学路人甲

发表于 Post on 2022-3-29 14:01:33 | 显示全部楼层 Show all
本帖最后由 Novice 于 2022-3-29 14:05 编辑

Q&A第5个相关的问题,使用力常数,可以直接对比不同化学键,如C-O和C-N键的键强度吗?

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-30 06:00:13 | 显示全部楼层 Show all
Novice 发表于 2022-3-29 14:01
Q&A第5个相关的问题,使用力常数,可以直接对比不同化学键,如C-O和C-N键的键强度吗?

对于两个键  A-B和A-C,假如B和C原子在元素周期表的同一行或列,直接将力常数作为键强度来比较是可以的。
O和N在元素周期表的同一行,所以符合上面的情况。这一点在LModeA-nano的JCTC原文中有提到。

根据JCP 137,084114,2012 这篇文章的计算,B3LYP/6-31G(d,p)级别下计算得到的甲醇和甲胺分子中的C-O、C-N单键的力常数分别为 4.908、4.426 mdyn/A.

评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
Novice + 5 谢谢。赞!

查看全部评分 View all ratings

1020

帖子

0

威望

3508

eV
积分
4528

Level 6 (一方通行)

发表于 Post on 2022-3-30 08:34:37 | 显示全部楼层 Show all
强烈支持! !!

说明中的compliance啥意思能否简单解释一下?  能否处理非谐性?

145

帖子

0

威望

1265

eV
积分
1410

Level 4 (黑子)

发表于 Post on 2022-3-31 17:01:16 | 显示全部楼层 Show all
Novice 发表于 2022-3-29 14:01
Q&A第5个相关的问题,使用力常数,可以直接对比不同化学键,如C-O和C-N键的键强度吗?

+1

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-31 22:57:31 | 显示全部楼层 Show all
granvia 发表于 2022-3-30 08:34
强烈支持! !!

说明中的compliance啥意思能否简单解释一下?  能否处理非谐性?

compliance 矩阵 简单来说是力常数矩阵的逆矩阵,具体的数学定义可以参见 [1] Chem. Soc. Rev. 37,1558,2008, [2] JCP 137,084114,2012 等文章

非谐性包含多个方面,有势能面的非谐性,也包含了泛频、合频、费米共振等方面。由于局域振动模式和简正模式一样都是基于势能面的二次估计(对势能面进行泰勒展开后截断到二次项)的,因此理论本身不考虑非谐性。如果需要考虑势能面的非谐性,有两种做法:
[1] 使用经验性的频率校正因子,将该因子取平方后乘以力常数;
[2] 使用实验测量的频率代替计算得到的简正模式频率,重新计算校正后的Hessian矩阵,然后做局域模式分析。这个方法只适用于非常简单的体系。
这两点在 J. Chem. Theory Comput. 18,1821,2022 一文中有讨论到。

1020

帖子

0

威望

3508

eV
积分
4528

Level 6 (一方通行)

发表于 Post on 2022-4-1 13:58:37 | 显示全部楼层 Show all
smutao 发表于 2022-3-31 22:57
compliance 矩阵 简单来说是力常数矩阵的逆矩阵,具体的数学定义可以参见 [1] Chem. Soc. Rev. 37,1558,2 ...

多谢解惑,学习了!

15

帖子

0

威望

15

eV
积分
30

Level 2 能力者

发表于 Post on 2022-4-16 20:55:49 | 显示全部楼层 Show all
谢谢分享

176

帖子

1

威望

1005

eV
积分
1201

Level 4 (黑子)

发表于 Post on 2022-4-22 20:28:01 | 显示全部楼层 Show all
smutao 发表于 2022-3-29 11:03
LModeA-nano使用微信群二维码如下

老师您好 二维码过期了 能否再分享一个呢?谢谢您
计算小白 水平很菜
随时可喷 万望指点

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-4-23 05:26:14 | 显示全部楼层 Show all
WaterMirror 发表于 2022-4-22 20:28
老师您好 二维码过期了 能否再分享一个呢?谢谢您

已更新

176

帖子

1

威望

1005

eV
积分
1201

Level 4 (黑子)

发表于 Post on 2022-4-23 06:32:08 | 显示全部楼层 Show all

谢谢老师
计算小白 水平很菜
随时可喷 万望指点

16

帖子

0

威望

101

eV
积分
117

Level 2 能力者

发表于 Post on 2022-6-17 16:02:33 | 显示全部楼层 Show all
smutao 发表于 2022-3-29 11:03
LModeA-nano使用微信群二维码如下

老师好!能否麻烦您更新一下微信群二维码?

287

帖子

8

威望

1664

eV
积分
2111

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-6-21 00:19:42 | 显示全部楼层 Show all
TinnKuka 发表于 2022-6-17 16:02
老师好!能否麻烦您更新一下微信群二维码?

已更新

2

帖子

0

威望

65

eV
积分
67

Level 2 能力者

发表于 Post on 2022-7-14 22:18:47 | 显示全部楼层 Show all
老师您好,二维码过期了,能否麻烦您再更新一下,万分感谢。

本版积分规则 Credits rule

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

GMT+8, 2023-2-2 23:46 , Processed in 0.261573 second(s), 26 queries .

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