用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

复制代码

冬天一起堆雪人 · 发表于 Post on 2025-4-1 11:38:28

冬天一起堆雪人发表于 2025-3-31 17:24
为什么没有红色框这部分呀？如果加这部分格式怎么改呀？

已解决

冬天一起堆雪人 · 发表于 Post on 2025-3-31 17:24:32

为什么没有红色框这部分呀？如果加这部分格式怎么改呀？

冬天一起堆雪人 · 发表于 Post on 2025-3-31 16:59:53

我搜出来的赝势有红色框部分，为什么您这边没有加这部分呀

beefly · 发表于 Post on 2024-12-8 17:43:27

本帖最后由 beefly 于 2024-12-8 17:44 编辑

osric 发表于 2024-12-8 14:56
各位老师好：
想请问我按照教程设置好，进行了SODFT的单点计算。计算是正常结束的，但是用gaussview打开fc ...

包括GaussView在内的几乎所有轨道显示程序均不支持二分量或四分量计算的复轨道。唯一的例外是ADF的GUI，支持二分量轨道，而且必须是自己的数据格式，不过十年以前的版本显示的轨道是错的，并且很难看，不知道后来的版本改对了没有。

文献中有人用PyMol（https://doi.org/10.1002/qua.25683），修改版的IboView（http://dx.doi.org/10.1021/acs.jpca.9b01134）实现过复轨道显示。此外，这两个程序似乎也可以做
https://github.com/rbeck4/vmd_plot_complex
https://gitlab.com/quantum_rose/rose

osric · 发表于 Post on 2024-12-8 14:56:48

各位老师好：
想请问我按照教程设置好，进行了SODFT的单点计算。计算是正常结束的，但是用gaussview打开fchk文件后,发现看不了MOs轨道。Multiwfn也打开不了fchk文件。所以想问是哪里需要特殊设置吗？

beefly · 发表于 Post on 2024-11-18 00:27:47

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

Hilbrac · 发表于 Post on 2024-2-5 01:57:30

本帖最后由 Hilbrac 于 2024-2-23 13:44 编辑

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

我目前在g16手册里，只找到了”freq=nraman“和"freq=nnraman"这两种只做了一次数值差分的拉曼方法。

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

beefly · 发表于 Post on 2024-2-1 17:10:42

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

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

Hilbrac · 发表于 Post on 2024-2-1 16:56:00

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

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

beefly · 发表于 Post on 2024-1-25 09:07:56

本帖最后由 beefly 于 2024-11-17 23:39 编辑

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

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

2024.11.17. g16做不了

beefly · 发表于 Post on 2024-1-24 17:31:14

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

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

Hilbrac · 发表于 Post on 2024-1-24 15:23:56

本帖最后由 Hilbrac 于 2024-1-24 16:24 编辑

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

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

beefly · 发表于 Post on 2024-1-22 17:30:07

本帖最后由 beefly 于 2024-1-22 17:32 编辑

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

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

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

评分 Rate

评分 Rate

评分 Rate

浏览过的版块