计算化学公社

标题: 用Gaussian 16做二分量赝势自旋轨道DFT（SODFT）计算 [打印本页]

---

作者
Author: beefly    时间: 2020-4-16 18:29
标题: 用Gaussian 16做二分量赝势自旋轨道DFT（SODFT）计算


能够做二分量赝势SODFT计算的程序里[1]，用的比较多的大概就是“西北大学化学系”[2]和“淘宝猫”了[3]，因为这两个程序的SODFT都有解析梯度，能够做结构优化。如果把平面波程序也算上，还有VASP。
注释

• 二分量全电子SODFT是另一个话题，这里不涉及。
• NWChem做SODFT见：http://bbs.keinsci.com/thread-5573-1-1.html
• Turbomole除了SODFT，还有二分量的TDDFT以及post-HF（MP2、CC2），不过除了SODFT之外的方法目前都没有梯度，并且用免费的Dirac也能做，而且功能更强。
  

其实从16版开始，Gaussian也支持二分量赝势SOHF、SODFT的能量和梯度计算，只是很低调，没有宣传。在Gaussian的安装包中，找到tests目录下的test1198，就是一个用二分量赝势在HF级别计算Sg原子的测试，输入文件里包含详细注释。以下是TlCl分子结构优化+振动频率计算的输入示例。
输入文件说明

• 调用二分量SCF方法需要在HF、泛函名称之前加前缀g。Gaussian的二分量SCF方法是在GHF/GKS non-collinear框架下执行的，未考虑时间反演对称性，也不能用点群对称性。如果对分子结构的对称性有严格要求，建议用Z矩阵坐标，并且opt(z-mat)。
• GHF/GKS没有二阶导数，频率计算需要用梯度做数值差分。
• GHF/GKS必须结合二分量赝势。如果用一般的标量赝势，则GHF/GKS忽略旋轨耦合效应，得到与标量HF/KS计算一样的结果。二分量赝势与标量赝势的区别在于L > 0 (s)以上的赝势函数多出一列。Gaussian基组库没有二分量赝势，必须通过genecp或gen pseudo从输入文件读入。获得二分量赝势的方法见http://bbs.keinsci.com/thread-5573-1-1.html 。需要注意的是，这些基组数据库都不支持Gaussian格式的二分量赝势输出。建议先输出Gaussian格式的标量赝势，粘贴到Gaussian输入文件里，再输出其他格式的二分量赝势，把旋轨耦合部分手工加入到Gaussian输入文件相应的部分（建议找个具有列编辑功能的文本编辑器）。
• 旋轨耦合赝势在早期文献中有两种定义，差别在于是否乘上因子2/L。如果没有经过2/L换算，需要加上输入选项pseudo=soscal（见test1198第二步计算）。不过目前的赝势基组库网站都用了2/L换算后的定义，因此一般不用加这个选项。唯一的例外是Clarkson大学的赝势基组库（https://people.clarkson.edu/~pchristi/reps.html），给了两种定义，建议用“Download for use in GAUSSIAN and COLUMBUS packages.”中给出的旋轨耦合赝势，这是换算过的。
• 在这个输入里，忽略了Cl原子的旋轨耦合效应，用了非相对论基组cc-pVDZ，当然也可以用一般的标量赝势基组，比如cep-4g、lanl2dz、sddall。不要用全电子相对论收缩基组（基组名一般带有x2c、dk、zora）。
• 在这个输入里，其他关键词的含义和一般的HF、DFT计算相同。但是GHF/GKS不支持布局分析和各种单电子性质，不要加这类关键词。
  

GHF/GKS计算开壳层体系的注意事项

• 在二分量相对论方法中，自旋多重度是没有严格物理意义的。输入文件里给出的是近似的自旋多重度（在分子坐标前一行），但是结果依赖于初猜以及SCF迭代，最终得到的自旋多重度可能并不是想要的，需要检查<S**2>和S是否近似符合要求。
• GHF/GKS计算开壳层体系一般会收敛很慢，需要加大SCF迭代次数或适当降低SCF收敛阈值。一些解决标量SCF收敛的选项可能无效。
• GHF/GKS计算开壳层过渡元素体系，初猜具有随机性，结果可能无法重复，甚至不是基态。
  

1. %mem=8GB
   
2. %nprocshared=4
   
3. #p gb3lyp/genecp opt freq=numer
   
4. 
   
5. TlCl. Tl=dhf-SVP, Cl=cc-pvdz
   
6. 
   
7. 0 1
   
8. Tl
   
9. Cl  1  r1
   
10. 
    
11. r1 2.485
    
12. 
    
13. Cl     0
    
14. cc-pvdz
    
15. ****
    
16. Tl    0
    
17. S  6  1.00
    
18.   45.356924655      0.50938157397E-02
    
19.   20.278843336     -0.14805704314
    
20.   14.420612394      0.36768453564
    
21.   6.2161909754     -0.66820394623
    
22.   1.5002511178      0.96232806743
    
23. 0.76330879305      0.33175999361
    
24. S  2  1.00
    
25.   7.8193195283     -0.18828631537E-01
    
26.   1.3176904042      0.32093431467
    
27. S  1  1.00
    
28. 0.17643545481       1.0000000000
    
29. S  1  1.00
    
30. 0.65296835238E-01   1.0000000000
    
31. P  4  1.00
    
32.   9.0323058024      0.22738839085
    
33.   7.5356419671     -0.35117975889
    
34.   1.0521421998      0.25690085604
    
35. 0.50303932656      0.84549942764E-01
    
36. P  1  1.00
    
37.   2.0067660913       1.0000000000
    
38. P  1  1.00
    
39. 0.16675605056       1.0000000000
    
40. P  1  1.00
    
41. 0.48596401464E-01   1.0000000000
    
42. D  5  1.00
    
43.   9.4421940760      0.58250776986E-01
    
44.   7.4410003381     -0.11789199301
    
45.   1.9831617772      0.34232057636
    
46. 0.91671748567      0.47557282350
    
47. 0.39647207885      0.30528029435
    
48. D  1  1.00
    
49. 0.15500000000       1.0000000000
    
50. ****
    
51. 
    
52. TL     0
    
53. TL-ECP     5     60
    
54. h POTENTIAL
    
55.   1
    
56. 2      1.00000000            0.00000000
    
57. s-h POTENTIAL
    
58.   2
    
59. 2     12.16780500          281.28466300
    
60. 2      8.29490900           62.43425100
    
61. p-h POTENTIAL
    
62.   4
    
63. 2      7.15149200            4.63340800       -9.266817
    
64. 2      5.17286500            9.34175600        9.341756
    
65. 2      9.89107200           72.29925300     -144.598506
    
66. 2      9.00339100          144.55803700      144.558037
    
67. d-h POTENTIAL
    
68.   4
    
69. 2      7.13021800           35.94303900      -35.943039
    
70. 2      6.92690600           53.90959300       35.939729
    
71. 2      5.41757000           10.38193900      -10.381939
    
72. 2      5.13868100           15.58382200       10.389215
    
73. f-h POTENTIAL
    
74.   4
    
75. 2      5.62639900           15.82548800      -10.550326
    
76. 2      5.54895200           21.10402100       10.552010
    
77. 2      2.87494600            2.91512700       -1.943418
    
78. 2      2.82145100            3.89690300        1.948451
    
79. g-h POTENTIAL
    
80.   4
    
81. 2      6.67905700           -7.49453400        3.747267
    
82. 2      6.70683500           -9.54057500       -3.816230
    
83. 2      7.20928400           -7.79799200        3.898996
    
84. 2      7.07096400           -9.25952400       -3.703809
    

复制代码

---

作者
Author: Smes    时间: 2021-6-6 11:13
老师关于把旋轨耦合部分手工加入到Gaussian输入文件相应的部分，具体怎么加呢? 是将组合系数直接加在其相同高斯函数指数的系数后面吗？C:\Users\NO.9-808\Desktop\1.jpg

---

作者
Author: snljty    时间: 2021-6-6 15:00

Smes 发表于 2021-6-6 11:13
老师关于把旋轨耦合部分手工加入到Gaussian输入文件相应的部分，具体怎么加呢? 是将组合系数直接加在其相同 ...

图没贴上。看http://bbs.keinsci.com/thread-25-1-1.html第81行了解如何正确贴图。

---

作者
Author: wzkchem5    时间: 2021-6-6 15:01

Smes 发表于 2021-6-6 04:13
老师关于把旋轨耦合部分手工加入到Gaussian输入文件相应的部分，具体怎么加呢? 是将组合系数直接加在其相同 ...

贴图没贴对，我们看不到

---

作者
Author: Smes    时间: 2021-6-6 21:06

snljty 发表于 2021-6-6 15:00
图没贴上。看http://bbs.keinsci.com/thread-25-1-1.html第81行了解如何正确贴图。

谢谢老师，请问老师是这样加吗？

---

作者
Author: 天月风见    时间: 2023-6-5 15:41
老师您好，请问这个sodft只能给闭壳层用吗？我发现算出来的结构没有轨道没有分alpha和beta了。电子态也没有了。

---

作者
Author: beefly    时间: 2023-6-5 19:31

天月风见 发表于 2023-6-5 15:41
老师您好，请问这个sodft只能给闭壳层用吗？我发现算出来的结构没有轨道没有分alpha和beta了。电子态也没有 ...

二分量方法没有自旋的概念

用二分量方法计算开壳层体系，初猜占据很可能对应某个激发态，导致能量偏高，或者几千次迭代不收敛。有个小技巧：先用udft做标量的单点计算，其中指定合适的自旋多重度，保存chk文件（注意做好备份）。然后在二分量dft计算中通过guess=read从标量计算产生的chk文件读取波函，大多数情况下能解决上述问题。

---

作者
Author: 天月风见    时间: 2023-6-7 10:12
好的，谢谢老师指导

---

作者
Author: shenzp    时间: 2023-6-18 22:18
老师好，gaussian的GHF/GKS如果结合int=DKHSO以及全电子基组的话是不是就可以做全电子的二分量相对论计算了？

---

作者
Author: beefly    时间: 2023-6-19 20:00

shenzp 发表于 2023-6-18 22:18
老师好，gaussian的GHF/GKS如果结合int=DKHSO以及全电子基组的话是不是就可以做全电子的二分量相对论计算了 ...

对，不过只能算单点能

---

作者
Author: shenzp    时间: 2023-6-19 20:18

beefly 发表于 2023-6-19 20:00
对，不过只能算单点能

OK谢谢老师！我就是想算一下旋轨耦合对键能的校正，看来应该是没问题

---

作者
Author: 栗悟饭与龟波功    时间: 2023-12-15 22:36
老师，我想问一下在gaussian的SODFT计算中出现以下报错是怎么回事
XCFlgR: Bad density values, RhoT/X/Y/Z=  0.000000000000D+00  0.374773466715D-04 -0.915221112058D-06 -0.862486450281D-03
Error termination via Lnk1e in /opt/gaussian/16/C02/g16/l1110.exe at Fri Dec 15 06:25:06 2023.
这是我的输入描述文件
#p opt=(maxcyc=200,recalc=3) freq gpbe1pbe/genecp scf=(vshift=400,maxcyc=500) EmpiricalDispersion=GD3BJ guess=read geom=check
请问是泛函的问题还是初猜波函数的问题，初猜波函数我是在相同泛函基组下进行udft得到的

---

作者
Author: beefly    时间: 2023-12-24 14:20

栗悟饭与龟波功 发表于 2023-12-15 22:36
老师，我想问一下在gaussian的SODFT计算中出现以下报错是怎么回事
XCFlgR: Bad density values, RhoT/X/Y ...

二分量方法不支持解析频率，改为freq=numer

---

作者
Author: 栗悟饭与龟波功    时间: 2023-12-28 18:17

beefly 发表于 2023-12-24 14:20
二分量方法不支持解析频率，改为freq=numer

谢谢老师答疑解惑，最近在算别的东西，一直没看论坛。目前问题已解决。

但是我发现光改freq=numer似乎还是会报相同的错误，即在一次几何优化循环结束后，在计算结构间的力时，还是会显示XCFlgR: Bad density values

后来我发现只有把opt里面的recalc=3去掉后才能正常进行梯度的计算，而且之前的报错都是在第一步opt的scf收敛后，计算力的时候报错的。我怀疑可能是recalc是不是会在第一步opt就开始精确计算Hessian矩阵，之后每隔n步再精确计算一次，而这个SODFT不适合用精确计算的Hessian矩阵？不太清楚是什么原因，也可能是我没有做好完善的对照。

---

作者
Author: Hilbrac    时间: 2024-1-20 03:11
老师，我想问一下，这个是否能用类似 freq(numer, raman) 这样的关键词把拉曼谱也一起算了呢？

---

作者
Author: beefly    时间: 2024-1-22 17:30

Hilbrac 发表于 2024-1-20 03:11
老师，我想问一下，这个是否能用类似 freq(numer, raman) 这样的关键词把拉曼谱也一起算了呢？

不能直接做，因为g16目前没有二分量方法的解析极化率。可以自己写脚本，调用g16做数值极化率（polar=numer）对坐标位移的数值差分，得到极化率导数矩阵，然后存入fchk文件，剩下的交给gaussview

---

作者
Author: Hilbrac    时间: 2024-1-24 15:23

beefly 发表于 2024-1-22 17:30
不能直接做，因为g16目前没有二分量方法的解析极化率。可以自己写脚本，调用g16做数值极化率（polar=nume ...

感谢老师回复，我目前只知道g16可以用 polar=numer 算数值极化率，但是如何调用g16再进一步算极化率对核坐标的导数再存入fchk文件，对此我还完全不了解，还望老师能指条明路（这种调用脚本做法的参考资料或脚本的写法出处之类的都可以）。
另外，我还想确认一下，目前g16的sodft只有能量的一阶解析（对核坐标）导数，而没有二阶解析（对核坐标）导数对吗？那nwchem的sodft是否有能量的二阶解析（对核坐标）导数呢？

---

作者
Author: beefly    时间: 2024-1-24 17:31

Hilbrac 发表于 2024-1-24 15:23
感谢老师回复，我目前只知道g16可以用 polar=numer 算数值极化率，但是如何调用g16再进一步算极化率对核 ...

其实就是高数学的数值导数。网上搜一下一阶、二阶导数的数值求导公式。如果分子不大，也可以人工操作极化率的求导

---

作者
Author: beefly    时间: 2024-1-25 09:07

Hilbrac 发表于 2024-1-24 15:23
感谢老师回复，我目前只知道g16可以用 polar=numer 算数值极化率，但是如何调用g16再进一步算极化率对核 ...

数值差分算拉曼活性可以用freq=numer polar=numer
由于ghf/gks的收敛较慢，能量精度不够，step需要设置大一些

2024.11.17. g16做不了

---

作者
Author: Hilbrac    时间: 2024-2-1 16:56

beefly 发表于 2024-1-25 09:07
数值差分算拉曼活性可以用freq=numer polar=numer
由于ghf/gks的收敛较慢，能量精度不够，step需要设置 ...

先谢过老师的回复。
老师这里的意思是：用 polar=numer 关键词让G16用数值方法算出极化率。
那至于之后如何算极化率导数来得到拉曼活性，还是得自己用脚本调用G16对极化率做数值差分才能获得？

---

作者
Author: beefly    时间: 2024-2-1 17:10

Hilbrac 发表于 2024-2-1 16:56
先谢过老师的回复。
老师这里的意思是：用 polar=numer 关键词让G16用数值方法算出极化率。
那至于之后 ...

不用。数值频率+数值拉曼可以一起做

---

作者
Author: Hilbrac    时间: 2024-2-5 01:57

beefly 发表于 2024-2-1 17:10
不用。数值频率+数值拉曼可以一起做

我目前在g16手册里，只找到了”freq=nraman“和"freq=nnraman"这两种只做了一次数值差分的拉曼方法。
(, 下载次数 Times of downloads: 12)

上传 Uploaded

点击下载
Click to download

如果要用二分量方法的话，极化率就已经是”数值的“而非”解析的“了，也没法用这里的”nnraman“关键词了，那老师说的”一起做“要怎么写关键词呢？

---

作者
Author: beefly    时间: 5 day ago

Hilbrac 发表于 2024-2-5 01:57
我目前在g16手册里，只找到了”freq=nraman“和"freq=nnraman"这两种只做了一次数值差分的拉曼方法。

...

这些方法仅支持解析频率计算。sodft算拉曼的近似做法是先用sodft算数值频率并保存chk，然后在标量dft级别通过“polar=raman cphf=rdfreq”读取Hessian并解析计算极化率导数。示例输入见以下话题的123楼：
http://bbs.keinsci.com/thread-3545-9-1.html

---

欢迎光临 计算化学公社 (http://bbs.keinsci.com/)

Powered by Discuz! X3.3