用Gaussian 16做二分量赝势自旋轨道DFT（SODFT）计算

beefly · 发表于 Post on 2020-4-16 18:29:30

本帖最后由 beefly 于 2020-4-16 19:39 编辑

能够做二分量赝势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计算开壳层过渡元素体系，初猜具有随机性，结果可能无法重复，甚至不是基态。

%mem=8GB
%nprocshared=4
#p gb3lyp/genecp opt freq=numer
TlCl. Tl=dhf-SVP, Cl=cc-pvdz
0 1
Tl
Cl 1 r1
r1 2.485
Cl 0
cc-pvdz
****
Tl 0
S 6 1.00
45.356924655 0.50938157397E-02
20.278843336 -0.14805704314
14.420612394 0.36768453564
6.2161909754 -0.66820394623
1.5002511178 0.96232806743
0.76330879305 0.33175999361
S 2 1.00
7.8193195283 -0.18828631537E-01
1.3176904042 0.32093431467
S 1 1.00
0.17643545481 1.0000000000
S 1 1.00
0.65296835238E-01 1.0000000000
P 4 1.00
9.0323058024 0.22738839085
7.5356419671 -0.35117975889
1.0521421998 0.25690085604
0.50303932656 0.84549942764E-01
P 1 1.00
2.0067660913 1.0000000000
P 1 1.00
0.16675605056 1.0000000000
P 1 1.00
0.48596401464E-01 1.0000000000
D 5 1.00
9.4421940760 0.58250776986E-01
7.4410003381 -0.11789199301
1.9831617772 0.34232057636
0.91671748567 0.47557282350
0.39647207885 0.30528029435
D 1 1.00
0.15500000000 1.0000000000
****
TL 0
TL-ECP 5 60
h POTENTIAL
1
2 1.00000000 0.00000000
s-h POTENTIAL
2
2 12.16780500 281.28466300
2 8.29490900 62.43425100
p-h POTENTIAL
4
2 7.15149200 4.63340800 -9.266817
2 5.17286500 9.34175600 9.341756
2 9.89107200 72.29925300 -144.598506
2 9.00339100 144.55803700 144.558037
d-h POTENTIAL
4
2 7.13021800 35.94303900 -35.943039
2 6.92690600 53.90959300 35.939729
2 5.41757000 10.38193900 -10.381939
2 5.13868100 15.58382200 10.389215
f-h POTENTIAL
4
2 5.62639900 15.82548800 -10.550326
2 5.54895200 21.10402100 10.552010
2 2.87494600 2.91512700 -1.943418
2 2.82145100 3.89690300 1.948451
g-h POTENTIAL
4
2 6.67905700 -7.49453400 3.747267
2 6.70683500 -9.54057500 -3.816230
2 7.20928400 -7.79799200 3.898996
2 7.07096400 -9.25952400 -3.703809

复制代码

Smes · 发表于 Post on 2021-6-6 11:13:19

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

snljty · 发表于 Post on 2021-6-6 15:00:51

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

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

wzkchem5 · 发表于 Post on 2021-6-6 15:01:23

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

贴图没贴对，我们看不到

Smes · 发表于 Post on 2021-6-6 21:06:18

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

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

天月风见 · 发表于 Post on 2023-6-5 15:41:27

老师您好，请问这个sodft只能给闭壳层用吗？我发现算出来的结构没有轨道没有分alpha和beta了。电子态也没有了。

beefly · 发表于 Post on 2023-6-5 19:31:53

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

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

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

天月风见 · 发表于 Post on 2023-6-7 10:12:19

好的，谢谢老师指导

shenzp · 发表于 Post on 2023-6-18 22:18:12

老师好，gaussian的GHF/GKS如果结合int=DKHSO以及全电子基组的话是不是就可以做全电子的二分量相对论计算了？

beefly · 发表于 Post on 2023-6-19 20:00:06

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

对，不过只能算单点能

shenzp · 发表于 Post on 2023-6-19 20:18:03

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

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

栗悟饭与龟波功 · 发表于 Post on 2023-12-15 22:36:00

老师，我想问一下在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得到的

beefly · 发表于 Post on 2023-12-24 14:20:32

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

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

栗悟饭与龟波功 · 发表于 Post on 2023-12-28 18:17:40

本帖最后由栗悟饭与龟波功于 2023-12-28 18:23 编辑

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矩阵？不太清楚是什么原因，也可能是我没有做好完善的对照。

Hilbrac · 发表于 Post on 2024-1-20 03:11:49

老师，我想问一下，这个是否能用类似 freq(numer, raman) 这样的关键词把拉曼谱也一起算了呢？

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[Gaussian/gview] 用Gaussian 16做二分量赝势自旋轨道DFT（SODFT）计算

评分 Rate