能够做二分量赝势SODFT计算的程序里[1],用的比较多的大概就是“西北大学化学系”[2]和“淘宝猫”[3]了,因为这两个程序的SODFT都有解析梯度,能够做结构优化。如果把平面波程序也算上,还有VASP。
注释
1. 二分量全电子SODFT是另一个话题,这里不涉及。
2. NWChem做SODFT见
http://bbs.keinsci.com/thread-5573-1-1.html
3. Turbomole除了SODFT,还有二分量的TDDFT以及post-HF(MP2、CC2),不过除了SODFT之外的方法目前都没有梯度,并且用免费的Dirac也能做,而且功能更强。
其实从16版开始,Gaussian也支持二分量赝势SOHF、SODFT的能量和梯度计算,只是很低调,没有宣传。在Gaussian的安装包中,找到tests目录下的test1198,就是一个用二分量赝势在HF级别计算Sg原子的测试,输入文件里包含详细注释。以下是TlCl分子结构优化 振动频率计算的输入示例。
输入文件说明
1. 调用二分量SCF方法需要在HF、泛函名称之前加前缀g。Gaussian的二分量SCF方法是在GHF/GKS non-collinear框架下执行的,未考虑时间反演对称性,也不能用点群对称性。如果对分子结构的对称性有严格要求,建议用Z矩阵坐标,并且opt(z-mat)。
2. GHF/GKS没有二阶导数,频率计算需要用梯度做数值差分。
3. 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输入文件相应的部分(建议找个具有列编辑功能的文本编辑器)。
4. 旋轨耦合赝势在早期文献中有两种定义,差别在于是否乘上因子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.”中给出的旋轨耦合赝势,这是换算过的。
5. 在这个输入里,忽略了Cl原子的旋轨耦合效应,用了非相对论基组cc-pVDZ,当然也可以用一般的标量赝势基组,比如cep-4g、lanl2dz、sddall。不要用全电子相对论收缩基组(基组名一般带有x2c、dk、zora)。
6. 在这个输入里,其他关键词的含义和一般的HF、DFT计算相同。但是GHF/GKS不支持布局分析和各种单电子性质,不要加这类关键词。
GHF/GKS计算开壳层体系的注意事项
1. 在二分量相对论方法中,自旋多重度是没有严格物理意义的。输入文件里给出的是近似的自旋多重度(在分子坐标前一行),但是结果依赖于初猜以及SCF迭代,最终得到的自旋多重度可能并不是想要的,需要检查<S**2>和S是否近似符合要求。
2. GHF/GKS计算开壳层体系一般会收敛很慢,需要加大SCF迭代次数或适当降低SCF收敛阈值。一些解决标量SCF收敛的选项可能无效。
3. GHF/GKS计算开壳层过渡元素体系,初猜具有随机性,结果可能无法重复,甚至不是基态。
代码语言:javascript复制%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