请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 7844|回复: 21

[量化理论] 谈谈约束性DFT (CDFT)

[复制链接]

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

发表于 2014-12-29 10:17:09 | 显示全部楼层 |阅读模式
谈谈约束性DFT (CDFT)

文/Sobereva @北京科音   2014-Dec-29


1 原理

经常有人,特别是初学者会问量化计算中如何设定原子的电荷,他们想让某些原子或片段的电荷成为理想的整数。会问这个问题,大多都是对量化计算原理没搞透彻所致,笔者在《谈谈片段组合波函数与自旋极化单重态》(http://sobereva.com/82)一文的最后一节对此有过专门的讨论。一些初学者误以为通过设定初猜、利用gview里的片段设定就能达到目的,显然是错误的,不管初猜怎么定义,如果在SCF过程中不加约束的话,迭代至收敛后的电荷分布会和预期的初猜有明显不同。若想令最终的电荷分布与自己设定的一致,唯一办法就是在SCF迭代中进行约束,通常用的是拉格朗日乘子方法。对于DFT而言,这叫做约束性DFT(Constrainted DFT, CDFT),已经被广泛使用。

如今CDFT的标准实现方法是在PHYSICAL REVIEW A 72, 024502 (2005)中提出的。它在对能量变分时引入拉格朗日乘子Vc使得电子分布与期望的状态一致,这等价于在KS算符中额外引入了一个外势项Vc*w(r),其中w(r)是用来定义约束的权重函数。比如令w(r)为A原子的权重函数用来定义A原子所属的空间范围,那么就可以设定约束令w(r)*ρ(r)的全空间积分等于人为指定的A原子的电子数。在SCF迭代求解过程中,轨道和Vc交替优化。即每一步,先在当前的轨道下优化Vc直到满足预设的约束条件,然后用Vc再构建KS算符产生新的轨道,如此往复。CDFT的过程,既可以看成是在约束条件下寻找能量最低的解,也可以视为寻找一个外势,最终得到的是在这个外势下的基态的解。

CDFT不仅可以约束某些原子的布居数等于特定值,也可以让自旋布居数,即alpha和beta电子数的差值等于特定值(对alpha和beta布居数分别做约束),或者让两个片段间的电子布居数或自旋布居数的差值恰为特定值。权重函数也可以有多种选择,而且原理上可以同时设定多个约束。

用CDFT时切不可盲目,一定要搞清楚原理,弄清楚自己使用CDFT的目的何在、所设定的约束有什么物理意义。CDFT用对了可以讨论一些常规方式DFT计算研究不了的问题,比如讨论某些电荷转移、电子离域问题,但用得不合理得话只会得到毫无意义的结果。


2 NWChem中CDFT的使用

CDFT在NWChem和Q-Chem中都已经实现了。由于Q-Chem是收费的,下面就用免费的NWChem来演示一下CDFT。本文使用的是NWChem 6.5,编译方法可参考《NWChem 6.6编译方法》(http://sobereva.com/270

NWChem里实现CDFT很简单,在dft...end段落中加上比如cdft 2 5 charge 0.5 pop becke即设定2~5号原子的Becke原子电荷之和为0.5。也可以写诸如cdft 2 5 10 13 charge 0.5 pop lowdin来让2~5与10~13的Lowdin电荷差值为0.5。如果设定自旋布居数约束,把charge改为spin即可。NWChem的CDFT效率比较好,不比普通DFT计算多耗时多少。

注意计算原子电荷(或布居数)的方法很多,比如《原子电荷计算方法的对比》(http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml)一文就讨论了很多种。因此,对于同一个约束条件,使用不同的布居分析方法最终得到的波函数并不相同。NWChem的CDFT支持Mulliken、Lowdin和Becke三种布居方法。默认是Lowdin,虽然从原理上并不比Mulliken可靠,但用在CDFT上比Mulliken要好不少,CDFT原文中也说了Mulliken布居下结果有问题。至于Becke布居,并不是严格定义的方法,它的合理性和锐度参数k以及选用的原子半径密切相关,我不清楚NWChem内部是怎么具体设的,因此怀疑其可靠性。注意,Mulliken和Lowdin方法都不适合含有弥散函数的情况,此时要么去掉弥散函数,要么就只能用不怕弥散函数的Becke布居了。

碰见CDFT不收敛时,建议在dft...end里加入一行convergence nolevelshifting,如果还不行,可以尝试其它布居方法,或者改改约束方式,或者改改基组和泛函。

NWChem中对电荷约束可以做几何优化(且是解析梯度),对自旋约束不行。NWChem可以施加多个约束,不过正常收敛的几率很低。使用Becke布居时,对于有对称性的体系,如果计算时出现vectors were symmetry contaminated,需要用noautosym关键词避免使用对称性,否则结果可能有问题或无法正常结束,使用Becke布居时几乎总要加noautosym。用Becke布居时有时候并行运算没法正常结束,此时应改用串行计算。


3 实例1:乙酸钠

这是个最简单的例子。我们先用以下输入文件用NWChem算一下乙酸钠在B3LYP/6-31G*下的能量以及原子电荷。结构已在此级别下优化过。
  1. start
  2. geometry
  3. C                  0.00000000    0.52007900    0.00000000
  4. O                 -1.13032800   -0.06349600    0.00000000
  5. O                  1.11739200   -0.09103900    0.00000000
  6. C                  0.03308100    2.04494400    0.00000000
  7. H                  0.58076900    2.39854700    0.88048700
  8. H                  0.58076900    2.39854700   -0.88048700
  9. H                 -0.97577100    2.46220700    0.00000000
  10. Na                -0.02552400   -1.94665100    0.00000000
  11. end
  12. basis
  13. * library 6-31G*
  14. end
  15. dft
  16. xc b3lyp
  17. mulliken
  18. end
  19. TASK dft
复制代码

得到的能量为-390.839634985795,Na的Lowdin布居数为10.63,对应于原子电荷为11-10.63=0.37。其数值与传统化学概念上Na为1.0电荷相差不少,说明Na并没有把一个电子完完全全给了CH3COO-,或者说CH3COO-的电子有一部分转移向了Na+离子。

然后在dft...end当中插入
convergence nolevelshifting
cdft 8 8 charge 1.0
再次运算,从输出文件中也可以看到此时Na的Lowdin布居数为10.0,即原子电荷确实为设定的1.0。并且能量为-390.589504441073,比之前没加约束时高了656.7KJ/mol。

如果想优化的话,把TASK dft改为TASK dft optimize即可。

可见,我们可以通过CDFT和常规DFT计算的能量差,来讨论离子之间或分子之间,电荷发生一定幅度转移前后能量改变了多少。


4 实例2:叔丁基碳正离子

叔丁基碳正离子,即C(CH3)3+这个体系,一般画结构式的时候正电荷都画在中间的碳上,并认为这个碳用sp2杂化,垂直于碳平面的p轨道(以下简称pz轨道)是空着的。实际上这种说法完全忽略了三个甲基的电子向那个空pz轨道的离域。我们这里也算算玩玩。

运行以下输入文件做常规B3LYP/6-31G*计算。体系结构已在此级别下优化。
  1. start
  2. charge 1
  3. geometry
  4. C                  0.00000000    0.00000000    0.00000000
  5. C                  0.00000000    1.46600800    0.00000000
  6. H                  0.58581400    1.81844700    0.86660600
  7. H                  0.58581400    1.81844700   -0.86660600
  8. H                 -0.99127200    1.92058400    0.00000000
  9. C                 -1.26960000   -0.73300400    0.00000000
  10. H                 -1.16763900   -1.81875900    0.00000000
  11. H                 -1.86772800   -0.40189400    0.86660600
  12. H                 -1.86772800   -0.40189400   -0.86660600
  13. C                  1.26960000   -0.73300400    0.00000000
  14. H                  1.28191500   -1.41655300    0.86660600
  15. H                  1.28191500   -1.41655300   -0.86660600
  16. H                  2.15891100   -0.10182600    0.00000000
  17. end
  18. basis
  19. * library 6-31G*
  20. end
  21. dft
  22. xc b3lyp
  23. mulliken
  24. end
  25. TASK dft
复制代码

能量为-157.554198182446,中心碳的Lowdin电荷为+0.30,可见,把1个正电荷画在中心碳上的做法和实际的差距很大,那个正电荷明显没有完全定域在中心碳上。

然后在dft...end当中插入
cdft 1 1 charge 1.0
再次运算,看看人为让一个正电荷定域在中心碳上的结果怎样。结果为-157.460833435721,能量比不加约束时高了245.1KJ/mol。可以说,三个甲基的电子向中心碳离域,使得能量降低了245.1KJ/mol。此例表明恰当地利用CDFT可以用来计算离域化能。

需要特别注意的是,不同布居方法下得到的结果差异可能很大。比如此例,三种布居方法结果分别为
Mulliken -157.551751877858
Lowdin -157.460833435721
Becke -157.460056500370
Lowdin和Becke的差异不大,但Mulliken的数值却不太合理,能量仅比不加约束时高了6.4KJ/mol,明显偏小了。

CDFT计算时虽然使得中心碳的电荷为理想化的+1,但是,此时中心碳的pz轨道并非是空着的。如果在dft...end中使用
mulliken
print "mulliken ao"
来查看每个基函数的布居数的话,就会看到CDFT计算后中心碳的两个pz基函数的Lowdin布居数之和为0.237,远非为0,尽管比普通计算时的0.437有所减小。原理上讲,CDFT也能使得某些基函数的布居数为指定的值,但是NWChem的CDFT尚不支持这种约束方式。

顺带一提,如果就是想考察周围原子向中心碳的空pz轨道的离域的话,可以用NBO分析,计算当前体系会看到中心碳有个占据数为0.409的LP*轨道,这就对应于那个空pz(由于周围电子已经明显向此轨道离域了,所以其占据数远比0大)。从E2分析中可以看到这三个甲基中都有两个C-H键的BD轨道向此LP*有很强的超共轭作用,每对儿的E2作用能皆为11.5kcal/mol。虽然E2是不能简单相加的,但是如果姑且这么算一下,3*2*11.5*4.186=288.8KJ/mol,倒也和CDFT给出的离域化能在数量级上相仿佛,但我认为有这种程度的相符多半是巧合。


5 实例3:LiF

这次我们计算LiF,令Li的Lowdin电荷逐渐从0变为1,步长0.05,看看能量是如何变化的。手写那么多输入文件显然麻烦,我们用shell脚本来实现。把以下内容拷贝到一个文本文件里,然后执行之。
  1. for ((i=0;i<=20;i=i+1))
  2. do
  3. chg=`echo "$i*0.05"|bc`
  4. file=`printf "%4.2f\n" $chg`
  5. echo processing $chg...
  6. cat << EOF > LiF.nw
  7. start
  8. geometry
  9. Li                  0.00000000    0.00000000    1.52885
  10. F                   0.00000000    0.00000000    0.00000
  11. end
  12. basis
  13. * library def2-svp
  14. end
  15. dft
  16. XC b3lyp
  17. convergence nolevelshifting
  18. cdft 1 1 charge $chg
  19. end
  20. TASK dft
  21. EOF
  22. nwchem LiF.nw > $file.out
  23. done
复制代码

执行过后,就出现了0.00.out、0.05.out ... 1.00.out一系列文件,文件名就是约束的Li的电荷。然后运行
grep "Total DFT energy" *.out |tee t.txt
把t.txt用ultraedit等程序的列模式进行编辑,去掉多余的列,然后导入到origin里作图,得到下图(没有0.9的值是因为此时计算没有收敛)
1.gif


这个LiF体系正常DFT计算下Li的Lowdin电荷为0.46,这个数值也正是图中的最低点,约束设定的电荷偏离这个值越大,显然能量也越高,图中很好地表现了这点。

如果将Li的电荷约束为1.0并进行优化,键长会由1.52885埃增加到1.84477埃。虽说经过优化,CDFT能量有所降低,为-107.12133927a.u.,但也比不加约束时在平衡结构时的-107.33859717a.u.要高得多。


6 实例4:C4H8双自由基

双自由基的对称破缺DFT计算,在《谈谈片段组合波函数与自旋极化单重态》(http://sobereva.com/82)谈了很多了,CASSCF计算在《CASSCF计算双自由基以及双自由基特征的计算》(http://sobereva.com/264)中也谈了很多了,都用到了C4H8双自由基这个例子,不熟悉的话建议先看看。

C4H8双自由基直接用对称破缺B3LYP/6-31G*计算的话,末端两个碳的Lowdin自旋布居分别为0.919和-0.919。这里我们通过CDFT计算,让这两个碳原子自旋布居数之差恰为2.0。运行以下输入文件。odft必须写,否则会被当成闭壳层单重态计算。由于是对称破缺计算而当前体系又有对称性,故应当写noautosym避免使用对称性。
  1. start
  2. geometry noautosym
  3. C                 -0.74400100    1.78566400    0.00000000
  4. H                 -0.60282700    2.33865300    0.92499500
  5. H                 -0.60282700    2.33865300   -0.92499500
  6. C                 -0.74400100    0.30988100    0.00000000
  7. H                 -1.25452600   -0.08746700    0.88463900
  8. H                 -1.25452600   -0.08746700   -0.88463900
  9. C                  0.74400100   -0.30988100    0.00000000
  10. H                  1.25452600    0.08746700   -0.88463900
  11. H                  1.25452600    0.08746700    0.88463900
  12. C                  0.74400100   -1.78566400    0.00000000
  13. H                  0.60282700   -2.33865300   -0.92499500
  14. H                  0.60282700   -2.33865300    0.92499500
  15. end
  16. basis
  17. * library 6-31G*
  18. end
  19. dft
  20. xc b3lyp
  21. odft
  22. mulliken
  23. cdft 1 1 10 10 spin 2.0
  24. end
  25. TASK dft
复制代码

此例约束了两个末端的碳的自旋布居之差为2.0,由于体系本身就是对称的,这使得其中一个碳的自旋布居数恰为1.0,另一个为-1.0。从输出的Lowdin自旋布居上可以确认这一点。

我们也可以通过设定多个约束让体系末端两个CH2的自旋布居数分别为1和-1,即约束设定改为
cdft 1 3 spin 1.0 pop becke
cdft 10 12 spin -1.0 pop becke
(Mulliken和Lowdin布居时这种约束设定下无法正常计算,故改用Becke布居)


7 实例5:H+ H-体系

最后,我们用CDFT和背景电荷两种办法计算一下H+ H-这个体系,相距5埃。NWChem输入文件如下
  1. start
  2. geometry noautosym
  3. H                  0.00000000    0.00000000    0.00000000
  4. H                  0.00000000    0.00000800    5.00000000
  5. end
  6. basis
  7. * library 6-31G*
  8. end
  9. dft
  10. xc b3lyp
  11. cdft 1 1 charge -1
  12. end
  13. TASK dft
复制代码
结果为-0.567651263156。

然后在Gaussian里,计算H-体系,在与之距离5埃处放一个单位正电荷,输入文件如下
  1. #p b3lyp/6-31G* charge

  2. test

  3. -1 1
  4. H 0. 0. 0.

  5. 0. 0. 5. 1.0


复制代码
结果为-0.5676521。可见,用CDFT计算H+ H-和计算H-但在附近放一个质子,两种方式计算结果是等同的,也体现了CDFT的合理性。

评分

参与人数 10eV +47 收起 理由
ChunLinX + 2 好物!
ter20 + 10 好物!
liyuanhe211 + 5 好物!
stm8150 + 1
小苹果 + 4 赞!
aqhuangry + 5
katrina + 5 谢谢,Sob
Shannon + 5 谢谢,学习了
卡开发发 + 5 好物!
zhanfei + 5 学习了

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

20

帖子

0

威望

631

eV
积分
651

Level 4 (黑子)

发表于 2014-12-29 10:44:20 | 显示全部楼层
学习了。非常感谢~

41

帖子

0

威望

1208

eV
积分
1249

Level 4 (黑子)

发表于 2014-12-31 14:35:31 | 显示全部楼层

学习了。非常感谢~

1

帖子

0

威望

4

eV
积分
5

Level 1 能力者

发表于 2015-1-4 21:26:25 | 显示全部楼层
学习了,这个可以计算周期性吗?

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2015-1-4 21:38:16 | 显示全部楼层
sinap 发表于 2015-1-4 21:26
学习了,这个可以计算周期性吗?

CDFT从原理上讲也能用于周期性,但是程序可能不支持
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

66

帖子

0

威望

764

eV
积分
830

Level 4 (黑子)

发表于 2015-1-6 10:54:24 | 显示全部楼层

作为一头蒙头乱撞的菜鸟,贡献一个错误案例, shame

      做一个双核金属配合物的结构优化,由于这个配合物在化学结构上是对称的(L-M-O2-M-L),而晶体解析发现是不对称的,于是将分子划成4个片段,(L-M)(O)(O)(M-L),分别给两边的 M-L指定不同的电荷(也就是电子数) 和单电子的自旋状态进行计算,期望了解为何两边的金属会出现不同的配位模式。
   结果当然是可耻的失败了,计算结果表明,只要分子总的 自旋多重度一样,两边片段的划分不影响优化结果。不过有意思的倒是,某些指定的初始猜测会加速优化速度,甚至还有一个例子,直接指定S=3,不进行划分,结果一上来就SCF不收敛报错了,划成(1,5  -2,0  -2,0 3,-3)之后,总的S=3,就顺利优化下去了,可能是相当于给了一个好的初始猜测?

   顺便请教一下sob老师,这样的双核结构,如果要了解两边M的电子组态方面的信息,应该进行什么样的计算呢?

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2015-1-6 13:28:43 | 显示全部楼层
nkallwar 发表于 2015-1-6 10:54
作为一头蒙头乱撞的菜鸟,贡献一个错误案例, shame

      做一个双核金属配合物的结构优化,由于这个 ...

如果是那种反铁磁性耦合的体系,不管是否合理地设定了片段初猜,起码初猜是得对称破缺的才行,才有可能得到正确的对称破缺的电子结构。

如果体系不是这种情况,那么初猜通常无所谓,只要不是很离谱的就行了。不管怎么初猜只要最后收敛了都一样,片段初猜没太大必要。实际上,片段初猜可能还不如默认的初猜好,毕竟片段初猜是拼接起来的,离最终波函数还比较远。至于片段初猜收敛了而默认情况下没收敛,这估计只是巧合,尝试使用帮助收敛的关键词通常能解决。

至于电子组态,等算完了再做波函数分析就行了。比如获得原子轨道的电子/自旋布居数可以做NPA分析。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

66

帖子

0

威望

764

eV
积分
830

Level 4 (黑子)

发表于 2015-1-16 12:13:23 | 显示全部楼层
sobereva 发表于 2015-1-6 13:28
如果是那种反铁磁性耦合的体系,不管是否合理地设定了片段初猜,起码初猜是得对称破缺的才行,才有可能得 ...

多谢sob指教
最近在看Karl Wieghardt组做的一些关于ligand-Metal 之间自旋耦合的例子,他们为了解释某些金属配合物展现出的低自旋状态,经常采用broken symmetry approach,比如这篇文献(dx.doi.org/10.1021/ic302743s | Inorg. Chem. 2013, 52, 4472−4487),他们的计算结果似乎也表明了只有在S=0(论文里的S指的是总自旋量子数),也就反铁磁性的体系中,BS策略才是奏效的,计算出来的能量要比RKS明显的要低,而其它的情况,划分片段和直接用UKS计算结果几乎是一样的。 下面摘抄supporting information文件里的几段数据:

BS(m,n)他们指的是两个片段中含有m个a电子和n个β电子,
Table S5  S=0  BS(3,3) 的能量比 RKS 低 40kcal/mol,类似的还有Table S6

Table S2-S3及S7中,当S!=0时,UKS与BS(m,n)计算得到的能量几乎一致,差距都在1kcal/mol 以下,不过似乎也稍有区别,后面的几何结构(Table S8) 也非常接近,但也有细微的区别

关于计算部分的说明,似乎他们的BS片段划定,就是用来做了一个初始猜测,后面应该是没有做约束性DFT计算
Throughout this study, our computational results are described using the broken symmetry (BS) approach.23 The following notation is used to describe the broken symmetry solutions, where the given system is divided into two fragments. The notation BS(m,n) refers to a broken symmetry state with m unpaired α-spin electrons localized on fragment 1 and n unpaired β-spin electrons localized on fragment 2. In most cases, fragments 1 and 2 correspond to the metal and the ligands, respectively. In this notation, the standard high-spin, open-shell solution is written as BS(m + n,0). The BS(m,n) notation refers to the initial guess for the wave function.

在论文关于Co(AP*)3 电子结构的讨论部分,关于S=1/2 的情况,文章中采用了BS(2,1)的计算数据来说明三个配体自由基之间存在的耦合关系,见论文Figure13及其下讨论 :The UKS calculation for the S = 1/2 excited state converged to a BS(2,1) solution in which a through-space ligand radical−radical coupling was observed (Figure S8, Supporting Information).  而在前面那个supporting information 中Table S2 中,UKS的计算能量似乎还要比BS得到的结果略低一点,所以为什么要采用BS(2,1)的结果呢?


tableS5.JPG
tableS7.JPG
tableS2-S3.JPG

ic302743s.pdf

796.69 KB, 下载次数: 11

ic302743s_si_001.pdf

2.63 MB, 下载次数: 1

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2015-1-16 19:53:42 | 显示全部楼层
nkallwar 发表于 2015-1-16 12:13
多谢sob指教
最近在看Karl Wieghardt组做的一些关于ligand-Metal 之间自旋耦合的例子,他们为了解释某些 ...


表S2里S=0.5在UKS和BS(2,1)下的结果的差异,感觉很可能只是收敛精度问题造成的。应当做一下布居分析,看看电子结构能差出多少。
表S3的差异更像是收敛精度造成的。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

66

帖子

0

威望

764

eV
积分
830

Level 4 (黑子)

发表于 2015-1-17 12:43:08 | 显示全部楼层
sobereva 发表于 2015-1-16 19:53
表S2里S=0.5在UKS和BS(2,1)下的结果的差异,感觉很可能只是收敛精度问题造成的。应当做一下布居分析, ...


嗯,从能量和几何结构参数上看,几乎就是一样的优化结果。

那么这带来两个问题: 1,文章中的BS策略似乎也就是做了一个分片段初始猜测 “The BS(m,n) notation refers to the initial guess for the wave function.“  几何优化的结果几乎已经表明了做不做这个初始猜测没有什么意义,那么,对于布居分析来说,指定片段的电荷和自旋多重度拼起来的初始猜测对于布居分析会有很大的影响么?
2,即便是UKS和BS得到的布居结果有较大的差异,如何确定哪一个的电子组态才是真正符合实际情况的呢? 毕竟从几何优化的结果来看,这两者几乎没有差别。文章中似乎也没有提到使用BS数据的理由是什么

对我前面提到的那个双核金属的问题,是不是有必要专门做一个片段初猜的布居分析呢?



1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2015-1-17 15:16:57 | 显示全部楼层
nkallwar 发表于 2015-1-17 12:43
嗯,从能量和几何结构参数上看,几乎就是一样的优化结果。

那么这带来两个问题: 1,文章中的BS策略 ...

1 也没什么影响,毕竟收敛后的波函数是相同的(或者说,即便确实有差异,也小到察觉不到)
2 如果确实差异大,那么谁的能量最低就用哪个。如果UKS明摆着能量比BS的高很多,作者却反倒用UKS,那只能问作者是否有什么特殊的考虑。

可以试试,但我估计和UKS的是一样的。如果能量极为相近,但电子结构差异却很大,那就太巧了,除非是两个状态由于什么原因能量是简并的。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

66

帖子

0

威望

764

eV
积分
830

Level 4 (黑子)

发表于 2015-1-26 12:48:24 | 显示全部楼层
sobereva 发表于 2015-1-17 15:16
1 也没什么影响,毕竟收敛后的波函数是相同的(或者说,即便确实有差异,也小到察觉不到)
2 如果确实差 ...


选了一个小的结构尝试了一下,BS和UKS几乎是一样的

不知道作者是出于什么样的考虑

paper里后面还有一段是这么写的:

UKS calculations for the S = 1 dication 6 underwent spontaneous symmetry breaking to yield a BS(3,1) solution.

Similarly, UKS calculations for S = 1/2 [Cr(PDI)2]1+ undergo spontaneous symmetry breaking to yield BS(3,2) solutions.

UKS calculations for neutral complex 2 (S = 1) converged exclusively to BS(3,1) solutions.

总之就是用UKS算出来和BS是一样的,于是认为他的BS片段划分是合理的, 可是这样的话,那他费这么大劲弄出来的片段划分到底有什么意义呢?

我弄的那个双核, 划成S=1,BS(4,2)(aaaa, bb) 与 S=1 BS(0,2) (abab,aa),以及直接UKS,S=1,结果完全是一样的,应该还是最后从布居分析中看实际更符合哪种“BS” 情况吧?

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2015-1-26 13:35:52 | 显示全部楼层
nkallwar 发表于 2015-1-26 12:48
选了一个小的结构尝试了一下,BS和UKS几乎是一样的

不知道作者是出于什么样的考虑


他也就是为了严谨些吧,毕竟BS设定下心里有谱,而直接UKS说不准收敛到对应于哪种情况。一般来说就直接UKS就行了,没必要这么复杂。

你说得对。

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

105

帖子

0

威望

273

eV
积分
378

Level 3 能力者

发表于 2016-10-12 17:04:30 | 显示全部楼层
叔丁基碳正离子 的例子 在NWCHEM-6.6下计算CDFT时,要在geometry 后面加 noautosym 变成 geometry noautosym
否则会出现vectors were symmetry contaminated的问题。想请教sob大大,这个vectors are symmetry contaminated是什么意思?是什么vectors啊?

1万

帖子

25

威望

1万

eV
积分
34090

管理员

公社社长

 楼主| 发表于 2016-10-12 23:00:00 | 显示全部楼层
wswrpd 发表于 2016-10-12 17:04
叔丁基碳正离子 的例子 在NWCHEM-6.6下计算CDFT时,要在geometry 后面加 noautosym 变成 geometry noautosy ...

这只是个抽象的说法,可能指波函数什么的,具体问作者
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)
计算化学公社论坛:http://bbs.keinsci.com(高水平、高人气、综合性计算化学交流论坛)
思想家公社QQ群1号:18616395,2号:466017436。用于讨论理论、计算化学,两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。研究方向和理论、计算化学无关者勿加,以免浪费宝贵的空位

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2018-9-20 01:41 , Processed in 0.189380 second(s), 27 queries .

快速回复 返回顶部 返回列表