24小时热门版块排行榜    

查看: 3540  |  回复: 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的回帖
引用回帖:
4楼: Originally posted by 严登鑫 at 2017-08-02 16:47:22
您好请教一个问题,我想用gaussian计算外加电场对分子能量的影响。
具体步骤如下:
                              1.提交Gaussian 程序进行分子优化及频率计算,得到该分子的最稳定结构
         
             ...

1,步骤1没问题,步骤2未必需要opt,看你要的是电场中的单点能,还是电场中的最低能量

2,有用的能量值有好几个,真空中绝对零度下的电子总能量,前者的零点能校正,以及考虑环境影响的吉布斯自由能

3,频率计算才有zpe

» 本帖已获得的红花(最新10朵)

5楼2017-08-03 09:13:56
已阅   回复此楼   关注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强帖   回复此楼   编辑   查看我的主页
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 为什么资助数各大高校都创新高,自己申请怎么就这么难 +13 Kittylucky 2026-08-27 14/700 2026-09-01 11:06 by feng6531
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 怎么看青基中了没有啊 +6 叶九微 2026-08-26 6/300 2026-08-31 23:54 by yudaoqian88
[基金申请] 基金未中,这种答复是模板吗? +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[文学芳草园] 梦想 +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 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 系统查不到 +11 董八千 2026-08-26 11/550 2026-08-28 18:06 by Leogzhya
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +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
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
信息提示
请填处理意见