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

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

%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