24小时热门版块排行榜    

查看: 3539  |  回复: 6
【奖励】 本帖被评价3次,作者beefly增加金币 2.4
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

[资源] 用数值方法计算偶极矩

此文目的:
1. 介绍偶极矩计算公式
2. 对于没有密度矩阵的理论方法,通过数值差分计算偶极矩

通常情况下,偶极矩都是用解析方法得到的。以CO分子为例,键长取1.2 Angstrom,用Gaussian程序计算CCSD/3-21G级别偶极矩的输入如下:

# ccsd/3-21g density=cc

Dipole moment with CCSD density

0 1
C
O 1 1.2


计算得到的偶极矩为-0.0743 Debye,方向为-Z(即,从O指向C)。很多初学者忘了加density=cc,这种情况下得到的是Hartree-Fock级别的偶极矩。

如果用数值方法,偶极矩该怎么算呢?偶极矩定义为能量对电场强度的一阶导数:(μx, μy, μz) = (dE/dFx, dE/dFy, dE/dFz),因此可以把偶极矩计算分解为不同电场下同一理论级别的一系列单点能计算。

从很多数值分析书上可以查到3点数值差分公式,
dE/dFx = [-E(-Dx) + E(+Dx)]/(2Dx)

其中Dx是X坐标轴方向上某个自定义的电场大小;E(-Dx)、E(+Dx)分别为给定电场下计算的分子总能量。还可以采用更精确的5点数值差分公式,
dE/dFx = [E(-2Dx) -8E(-Dx)+8E(+Dx)-E(+2Dx)]/(12Dx)

类似地,可以得到dE/dFy和dE/dFz的数值差分公式。

从上面的公式可以看到,3点公式只需要算两个点,5点公式只需要4个点。这是因为在一阶导数的数值差分计算中,用不到电场为0的点。

有了公式,下面可以设计计算流程了。仍是CO分子,C-O键仍放在Z轴。因为在X和Y方向,分子是对称的,偶极矩为零,故只需要算μz = dE/dFz。Gaussian提供了几种定义电场的方法,其中一种的命令为Field=z-n或Field=z+n,z表示电场方向,-n/+n表示电场为-/+n * 0.0001 a.u.。以下是第一个点的输入,电场取-0.002 a.u.:

# ccsd/3-21g Field=z-20

CCSD energy calculation with an electric dipole field in the Z direction of -0.002 au

0 1
C
O 1 1.2


第二、三、四个点的电场分别为-0.001,+0.001,+0.002 a.u.,需要分别把Field=z-20替换为Field=z-10,Field=z+10,和Field=z+20。四个点的CCSD能量如下(单位Hartree):

-112.2970354
-112.2970453
-112.2971037
-112.2971523


代入5点公式(Dz=0.001),换算成Debey(换算因子2.542)后,为-0.0742 Debey,非常接近上面通过解析方法得到的-0.0743 Debey。如果用三点公式(取第2、3个点),仍然得到-0.0742 Debey。

应用

上面计算偶极矩的数值方法有什么用呢?我们知道,Gaussian程序目前还不支持CCSD(T),BCCD(T),QCISD(TQ),等理论级别的密度矩阵。如果想得到这些理论级别的偶极矩,要么换程序(例如免费的CFour程序),要么就用上面的数值方法。

还有些理论方法是没有密度矩阵的,例如CBS,BSSE,Gn理论等。用上面的数值方法,原则上也可以计算这些理论级别的偶极矩。

[ Last edited by beefly on 2011-6-26 at 11:56 ]

补充:
全电子相对论方法(如ZORA,DKH2,等)的哈密顿需要做图景变换,因此一般情况下不能用上面的数值方法算偶极矩,除非在电场计算程序中加入图景变换。

[ Last edited by beefly on 2012-8-2 at 21:56 ]
回复此楼

» 收录本帖的淘帖专辑推荐

@个人收集@ 我的收藏

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
引用回帖:
6楼: Originally posted by 严登鑫 at 2017-08-03 11:21:16
好的,谢谢beefly~~

我需要得到电场中的最低能量,#p B3LYP/6-311G(3df) Field=x+60 Opt=Z-Matrix NoSymm freq 可以得到zero-point energies 这样做正确吗?...

1,确定电场方向和分子方位是否一致

2,g09加电场的情况下,dft没有解析hessian。频率计算很耗时建议用g16做
7楼2017-08-03 16:53:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 7 个回答

mei1988925

金虫 (小有名气)


★★★★★ 五星级,优秀推荐

学习了!谢谢beefly!
2楼2011-06-26 13:08:43
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

lihb734

铁杆木虫 (职业作家)


★★★★★ 五星级,优秀推荐

感谢您的无私贡献
3楼2011-06-26 13:09:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

匿名

用户注销 (正式写手)


本帖仅楼主可见
4楼2017-08-02 16:47:22
已阅   申请QC强帖   回复此楼   编辑   查看我的主页
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 面上函评意见出来了,像什么等级? 20+4 Tsingking1 2026-08-27 15/750 2026-09-01 12:23 by icm639
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +13 Kittylucky 2026-08-27 14/700 2026-09-01 11:06 by feng6531
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 国社科又开始会评了,不知道这次命运如何 +7 雨打竹帘 2026-08-30 11/550 2026-08-31 23:16 by hittle2008
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 能否申诉? +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[文学芳草园] 梦想 +5 myrtle 2026-08-26 5/250 2026-08-31 14:40 by hahaboy
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 基金不中,共勉 +12 eulota 2026-08-26 12/600 2026-08-31 08:42 by ZJTJZ
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 看板上这么多中的,有点像50人群里49个人都是骗子的那种感觉…… +5 a089 2026-08-26 6/300 2026-08-27 14:05 by jonewore
[基金申请] 为什么 国际(地区)合作与交流项目 没有放榜? 10+3 majunge000 2026-08-26 11/550 2026-08-27 08:42 by 北京莱茵编辑
[基金申请] 2026年8月25日国自然放榜前突然收到列入评审专家邮件,有关系吗? +25 木水思豆 2026-08-25 28/1400 2026-08-26 14:53 by draco1987
[基金申请] 出来了 +9 trojank 2026-08-26 9/450 2026-08-26 14:25 by 宝贝虫子
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
信息提示
请填处理意见