CASEA的奇思妙想

Properties

发表于 更新于 分类于
本文字数： 28k 阅读时长 ≈ 25 分钟

ORCA中一些性质计算

Properties

单点能（Single point energies）

单点能是薛定谔方程的最低能量解。
有不同的方法计算单点能，同时有不同的计算精度。

Hartree-Fock(HF)

最早开发的近似计算之一，不包括所谓的动态或者静态相关性。
例：

1
2
3
4
5
6
!HF DEF2-SVP
* xyz 0 1
O         -3.56626        1.77639        0.00000
H         -2.59626        1.77639        0.00000
H         -3.88959        1.36040       -0.81444
*

DEF2-SVP为基组
HF是最老和最简单的方法，他是CCSD、CASSCF方法的基础

Unrestricted HF(UHF)

如果自旋多重度超过1，对于三元组则为3，可以执行UHF计算
在计算中，必须检查结果的自旋污染

如果 S2（写为 <S**2>）算子的期望值与理想值显着不同，则可能是您的系统只能使用多参考计算（例如 CASSCF）来处理

1
2
3
4
5
6
!UHF DEF2-SVP
* xyz 0 3
O         -3.56626        1.77639        0.00000
H         -2.59626        1.77639        0.00000
H         -3.88959        1.36040       -0.81444
*
1
2
3
4
5
6
7
8
9
10
11
## 结果

----------------------
UHF SPIN CONTAMINATION
----------------------

Expectation value of <S**2>     :     2.005700
Ideal value S*(S+1) for S=1.0   :     2.000000
Deviation                       :     0.005700

## 结果显示最小偏差

SCF加速（RIJDX和RIJCOSX）

使用`RI`选项计算库伦积分，从而加速求解SCF方程（计算能量方程）
在此情况下，还需要额外的辅助基组设置。
对于DEF2基组家族，默认设置为`DEF2/J`


1
!HF DEF2-SVP DEF2/J RIJDX


如果选择的基组没有相应的辅助基组，可以使用`AUTOAUX`标志设置自动生成辅助基组


1
!HF 6-31+G(d,p) AUTOAUX RIJDX


为了更快的加速，COSX算法可以和RI一起使用，从而有效的计算交换积分。

密度泛函理论(DFT)

通过从可用的函数中选择任何函数计算，获得DFT能量


1
2
3
4
5
6
!B3LYP DEF2-SVP
* xyz 0 1
O         -3.56626        1.77639        0.00000
H         -2.59626        1.77639        0.00000
H         -3.88959        1.36040       -0.81444
*


注：注意检查自旋污染，DFT结果可能会因为密度泛函的选择而异，在使用之前应进行检查，切勿比较不同方法的能量。
RIJCOSX仅在使用混合密度泛函时才有意义，同时还计算了HF交换

> RIJDX近似是非混合密度泛函的默认值
>
> RIJCOSX是混合密度泛函的默认值。
>
> 如不使用，可以使用`!NORI`或`!NOCOSX`

MP2微扰理论

MP2是一种后HF算法，以HF为基础，并包含丢失的动态相关性。
使用RI近似或者DLPNO算法


1
2
!RI-MP2 cc-pVTZ cc-pVTZ/C
!DLPNO-MP2 cc-pVTZ cc-pVTZ/C


因为是相关算法，对于`RI`部分需要一个特殊的基组`/C`，对于所有的基组来说不易得到，但是可以使用`AUTOAUX`获得


1
2
3
4
5
6
!DLPNO-MP2 6-311++G(2d,2p) AUTOAUX
* xyz 0 1
O         -3.56626        1.77639        0.00000
H         -2.59626        1.77639        0.00000
H         -3.88959        1.36040       -0.81444
*

Coupled Cluster (CC)耦合簇

耦合簇计算，特别是CCSD(T)变体，是单点能的黄金标准。
在常用的公式中，从计算角度看他们的成本非常高
ORCA拥有DLPNO变体，该变体效率更高，对于大型体系有着几乎线性尺度的增长。
为了使用DLPNO，只需在输入中设置DLPNO-CCSD(T)并选择合适的基组（同时需要设置/C基组或者AUTOAUX）

1
2
3
4
5
6
!DLPNO-CCSD(T) cc-pVTZ cc-pVTZ/C
* xyz 0 1
O         -3.56626        1.77639        0.00000
H         -2.59626        1.77639        0.00000
H         -3.88959        1.36040       -0.81444
*


注：检查输出值`T1 diagnostic value`
		从经验法则来说，对于`T1 diagnostic value`大于0.02，结果不可信，HF可信度变低
		在ORCA中有很多不同的耦合簇和耦合对方法
		`RIJCOSX`算法可被用于加速MP2和CC计算

几何优化(Geometry optimization)

基础用法

优化结构：通过给定的算法寻找拥有最小总能量的几何构型
第一步：使用软件给定初始构型，如：Avogadro画出结构并使用简单的优化得到可信的原子排布
第二步：使用Avogadro设置ORCA输入，保存文件为.xyz格式，或者直接在输入文件中写入坐标

如果想使用DFT优化构型，选择方法和基组并添加OPT标志。
OPT命令为几何优化命令

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
!B3LYP DEF2-SVP OPT
* xyz 0 1
N         -0.83911        0.76325       -0.31843
C          0.61442        0.72014       -0.25075
C          1.01669       -0.56167        0.49740
O          0.20095       -1.36984        0.93753
H         -1.37884        0.05803        0.17605
H          1.00414        0.66192       -1.27362
C          1.17285        1.95192        0.45211
H          0.87124        2.87150       -0.05988
H          0.81191        2.01288        1.48492
H          2.26726        1.92903        0.48485
H         -1.30551        1.57618       -0.71069
O          2.33980       -0.77979        0.66176
H          3.04559       -0.08055        0.28096
*

如果设定正确，可以在输出中看到优化参数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
                    *****************************
                    * Geometry Optimization Run *
                    *****************************

Geometry optimization settings:
Update method            Update   .... BFGS
Choice of coordinates    CoordSys .... Z-matrix Internals
Initial Hessian          InHess   .... Almoef's Model

Convergence Tolerances:
Energy Change            TolE     ....  5.0000e-06 Eh
Max. Gradient            TolMAXG  ....  3.0000e-04 Eh/bohr
RMS Gradient             TolRMSG  ....  1.0000e-04 Eh/bohr
Max. Displacement        TolMAXD  ....  4.0000e-03 bohr
RMS Displacement         TolRMSD  ....  2.0000e-03 bohr
Strict Convergence                ....  False
[...]

然后运行SCF，梯度运算，并给出当前步骤的相关信息

1
2
3
4
5
6
7
8
9
10
11
12
                      .--------------------.
----------------------|Geometry convergence|-------------------------
Item                value                   Tolerance       Converged
---------------------------------------------------------------------
RMS gradient        0.0125516277            0.0001000000      NO
MAX gradient        0.0693728096            0.0003000000      NO
RMS step            0.0324468165            0.0020000000      NO
MAX step            0.1852949368            0.0040000000      NO
........................................................
Max(Bonds)      0.0981      Max(Angles)    4.19
Max(Dihed)        2.11      Max(Improp)    0.00
---------------------------------------------------------------------

当猜测的几何形状很差时，第一步没有达到收敛标准，会根据默认算法（拟牛顿法）修改几何构型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
---------------------------------------------------------------------------
                     Redundant Internal Coordinates
                        (Angstroem and degrees)

    Definition                    Value    dE/dq     Step     New-Value
----------------------------------------------------------------------------
 1. B(C   1,N   0)                1.4557  0.013403 -0.0155    1.4403
 2. B(C   2,C   1)                1.5377  0.004146 -0.0057    1.5320
 3. B(O   3,C   2)                1.2297  0.045759 -0.0236    1.2062
 4. B(H   4,N   0)                1.0164  0.006262 -0.0075    1.0089
 5. B(H   5,C   1)                1.0961 -0.009891  0.0141    1.1102
 6. B(C   6,C   1)                1.5242 -0.009350  0.0123    1.5365
 7. B(H   7,C   6)                1.0949 -0.003989  0.0057    1.1005
 8. B(H   8,C   6)                1.0958 -0.003520  0.0050    1.1008
 9. B(H   9,C   6)                1.0951 -0.006870  0.0097    1.1049
10. B(H  10,N   0)                1.0160  0.007909 -0.0095    1.0065

几何构型的修改会不断迭代更新，知道满足标准

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
                      .--------------------.
----------------------|Geometry convergence|-------------------------
Item                value                   Tolerance       Converged
---------------------------------------------------------------------
Energy change      -0.0000009421            0.0000050000      YES
RMS gradient        0.0000326782            0.0001000000      YES
MAX gradient        0.0001152425            0.0003000000      YES
RMS step            0.0007124988            0.0020000000      YES
MAX step            0.0017837746            0.0040000000      YES
........................................................
Max(Bonds)      0.0001      Max(Angles)    0.02
Max(Dihed)        0.10      Max(Improp)    0.00
---------------------------------------------------------------------

          ***********************HURRAY********************
          ***        THE OPTIMIZATION HAS CONVERGED     ***
          *************************************************

../_images/opt.gif

尽管最初的猜测似乎不错，但正如人们所见，B3LYP 函数预测的几何形状与给定的输入有些不同。最初的猜测是胺基团处于更平面的构型中，该构型通过优化以及键长的修正而固定

在几何优化之后，会输出一个basename.xyz文件，其中包含最终几何坐标，用于后续计算

你还可以看到分子的能量如何随着优化而进行下降，最后趋于收敛
../_images/enconv.png

如何判断达到真正的最小值

鞍点的梯度范数也为零，或者梯度可能很低，但还是不够低
可以通过对计算出的最后几何构型计算振动频率来检查是否合理。在主要输入中添加关键词 FREQ

1
2
!B3LYP DEF2-SVP OPT RIJCOSX FREQ
* xyzfile 0 1 alanine_guess.xyz

在计算出准确的Hessian矩阵后，将打印频率

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  1.000000000  (already applied!)

   0:         0.00 cm**-1
   1:         0.00 cm**-1
   2:         0.00 cm**-1
   3:         0.00 cm**-1
   4:         0.00 cm**-1
   5:         0.00 cm**-1
   6:        59.59 cm**-1
   7:       231.28 cm**-1
   8:       245.18 cm**-1
   9:       254.63 cm**-1
  10:       321.35 cm**-1
[...]

除了第一次平移和旋转外的所有频率都是正的，这意味着我们确实有一个真正的最小值。
所有正频率都保证局部拥有最小值，但并不是全局最小值

数值梯度

最终生成的几何构型质量取决于你选择的方法。如果你想为当前没有解析梯度的方法优化几何结构你可以使用NUMGRAD标志使用数值梯度计算。
使用数值梯度，耗时会较长

1
!DLPNO-CCSD(T) cc-pVTZ cc-pVTZ/C OPT NUMGRAD

如何去除负频率

在几何优化运行后，通常会出现一个或多个负频率，特别是如果你的系统具有弱分子间相互作用，可能计算卡在了鞍点上，或者优化标准不够好
如果系统拥有很多小的负频率，你需要设置更为紧密的约束!TIGHTOPT,，在更严重的情况下，使用!VERYTIGHTOPT
如果您的系统有一个大的负频率，它可能比最小值更接近鞍点。最好的解决方案是采用您的最终结构并将其移向该模式的方向。

DET和OPT讨论

DFT 在分子间或弱相互作用中表现不佳。建议使用 D3 或电荷相关的 D4修正：

1
2
!B3LYP DEF2-SVP OPT D3
!B3LYP DEF2-SVP OPT D4

设定将校正能量和梯度，输出更好的结果

振动频率（Vibrational frequencies）

在优化构型后，你将计算系统的振动频率并绘制振动动画。
使用醋酸为例，使用DFT B3LYP：

1
2
3
4
5
6
7
8
9
10
11
!B3LYP DEF2-SVP OPT FREQ
* XYZ 0 1
C         -0.81589       -0.51571       -0.02512
C          0.30690        0.49327       -0.06114
H         -0.42809       -1.56713       -0.28060
H         -1.26914       -0.51520        1.06962
H         -1.64631       -0.14518       -0.75104
O          0.16587        1.68279       -0.21470
O          1.51380       -0.07303        0.21899
H          2.16801        0.64625        0.13143
*

也可以在优化后立即运行

1
2
!B3LYP DEF2-SVP OPT FREQ
* XYZFILE 0 1 guess_aceticacid.xyz

在几何优化成功之后，进行频率计算。当Hessian矩阵计算开始时，会输出以下信息：

1
2
3
4
5
6
7
8
9
10
11
12
-------------------------------------------------------------------------------
                               ORCA SCF HESSIAN
-------------------------------------------------------------------------------

Hessian of the Kohn-Sham DFT energy:
Kohn-Sham wavefunction type                      ... RKS
Hartree-Fock exchange scaling                    ...    0.200
Number of operators                              ...    1
Number of atoms                                  ...    8
Basis set dimensions                             ...   76
Integral neglect threshold                       ... 2.5e-11
Integral primitive cutoff                        ... 2.5e-12

然后计算所有必要的术语，完成后，可以直接看到以cm-1为单位的振动频率

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  1.000000000  (already applied!)

0:         0.00 cm**-1
1:         0.00 cm**-1
2:         0.00 cm**-1
3:         0.00 cm**-1
4:         0.00 cm**-1
5:         0.00 cm**-1
6:        82.29 cm**-1
7:       424.53 cm**-1
8:       544.40 cm**-1
9:       593.63 cm**-1
[...]

频率前几行均为0，他们对应着旋转和平移模式。线性分子应为5个0，非线性分子为6个0，剩余的对应实际振动模式。

缩放频率

如果你想在计算热力学函数之后所有频率乘以某个数字，可以使用

1
%FREQ SCALFREQ 1.035 END

所有的频率都将乘以1.035
时刻检查振动频率。如果几何构型最小，应当没有负频，如果是过渡态，有且仅有一个负频

可视化

振动模式可以使用Avogadro查看，打开.out文件，在右侧可以看到一个面板
img

你可以选择一个振动频率并单击Animate查看与该频率对应的动画。如果点击约~1800cm-1， 这是有机化学中C=O键拉伸模式，

~3150cm-1对应C-H键拉伸
img

消除负频

在计算之后，如果有一个或多个负频，意味着结构不稳定
以下时会产生负频的输入文件：

1
2
3
4
5
6
7
8
9
10
11
!B3LYP DEF2-SVP OPT FREQ
* XYZ 0 1
C         -8.23508        1.75252       -0.00446
C         -6.80437        2.18360        0.03487
H         -8.30301        0.65382       -0.15041
H         -8.75566        2.25995       -0.84350
H         -8.73273        2.02013        0.95115
O         -6.52860        3.36373        0.18904
O         -5.82021        1.26961       -0.09983
H         -4.88856        1.53058       -0.07686
*

优化后可以得到输出：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  1.000000000  (already applied!)

0:         0.00 cm**-1
1:         0.00 cm**-1
2:         0.00 cm**-1
3:         0.00 cm**-1
4:         0.00 cm**-1
5:         0.00 cm**-1
6:       -85.00 cm**-1 ***imaginary mode***
7:       430.96 cm**-1
8:       550.70 cm**-1
9:       590.13 cm**-1

使用所谓的虚模式，这意味着几何构型会收敛到鞍点，并且沿着该模式的方向移动会降低分子的能量。解决问题的方法时打开Avogadro，点击-85.0cm-1负频。
现在，如果您在分子位移时单击“暂停”，您将获得远离鞍点的几何形状。然后可以保存这些坐标并重新开始计算，最终收敛到最小值：

1
2
!B3LYP DEF2-SVP OPT FREQ
* XYZFILE 0 1 displaced_acetic.xyz

数值频率

ORCA中并不是所有方法都可以计算Hessian矩阵，但是都可以用numerical Hessian计算。使用NUMFREQ

1
!RI-B2PLYP DEF2-SVP DEF2-SVP/C NUMFREQ

在 6N 次位移后，Hessian 将通过中心差异方法计算，其中 N是您系统中的原子数。这可以完全并行化，有时实际上是比分析 FREQ 计算更好的选择，具体取决于您的内存。

重启NUMFREQ计算

如果系统较大或者计算较重，计算可能需要很长时间，如果终断，可以通过以下方法计算

1
2
!RI-B2PLYP DEF2-SVP DEF2-SVP/C NUMFREQ
%FREQ RESTART TRUE END

如果basename.hess文件与输入文件存放在同一文件夹，计算将从停止地方重新开始。

热力学（Thermodynamics）

在成功计算振动频率后，你可以进行热力学函数计算，例如焓（enthalpy H），熵（entropy S），吉布斯自由能（Gibbs free energy G）
对一些丁二烯异构体计算热力学函数

1
2
!B3LYP DEF2-TZVP OPT FREQ
* XYZFILE 0 1 gem-butene.xyz
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
!B3LYP DEF2-TZVP OPT FREQ
* XYZ 0 1 
C          0.25010        1.24571        0.30874
C          1.62155        1.73497       -0.07153
C         -0.38515        0.13754       -0.48356
C         -0.40976        1.78568        1.34351
H          0.28580       -0.19897       -1.30226
H         -1.33741        0.49224       -0.93172
H         -0.59374       -0.72860        0.17956
H          1.60282        2.11515       -1.11437
H          2.34546        0.89600       -0.00460
H          1.97134        2.55356        0.59367
H          0.02641        2.58709        1.93204
H         -1.39946        1.43084        1.61599
*

热力学功能

假设输入文件写入 !FREQ，在频率和正常模式之后会立刻看到计算热力学部分
这些属性将在给定的温度和压力下计算，假设气相分子并使用统计力学。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
!B3LYP DEF2-TZVP OPT FREQ
* XYZ 0 1 
C          0.25010        1.24571        0.30874
C          1.62155        1.73497       -0.07153
C         -0.38515        0.13754       -0.48356
C         -0.40976        1.78568        1.34351
H          0.28580       -0.19897       -1.30226
H         -1.33741        0.49224       -0.93172
H         -0.59374       -0.72860        0.17956
H          1.60282        2.11515       -1.11437
H          2.34546        0.89600       -0.00460
H          1.97134        2.55356        0.59367
H          0.02641        2.58709        1.93204
H         -1.39946        1.43084        1.61599
*

热力学开始

1
2
3
4
5
6
7
--------------------------
THERMOCHEMISTRY AT 298.15K
--------------------------

Temperature         ... 298.15 K
Pressure            ... 1.00 atm
Total Mass          ... 56.11 AMU

内能（inner energy U）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
------------
INNER ENERGY
------------

The inner energy is: U= E(el) + E(ZPE) + E(vib) + E(rot) + E(trans)
    E(el)   - is the total energy from the electronic structure calculation
              = E(kin-el) + E(nuc-el) + E(el-el) + E(nuc-nuc)
    E(ZPE)  - the the zero temperature vibrational energy from the frequency calculation
    E(vib)  - the the finite temperature correction to E(ZPE) due to population
              of excited vibrational states
    E(rot)  - is the rotational thermal energy
    E(trans)- is the translational thermal energy

Summary of contributions to the inner energy U:
Electronic energy                ...   -157.17473665 Eh
Zero point energy                ...      0.10741046 Eh      67.40 kcal/mol
Thermal vibrational correction   ...      0.00247510 Eh       1.55 kcal/mol
Thermal rotational correction    ...      0.00141627 Eh       0.89 kcal/mol
Thermal translational correction ...      0.00141627 Eh       0.89 kcal/mol
-----------------------------------------------------------------------
Total thermal energy                   -157.06201854 Eh

Electronic energy 是从最终单点能量得到的
Zero point Energy 零点能量

thermal correction 热力学校正

rotational 旋转
vibrational 振动
translational degrees of freedom 平移自由度

总校正能量是用于将Eel转换为U时添加进最终单点能的能量

焓（Enthalpy H）= U + kBT

热焓校正

1
2
3
4
5
6
7
8
9
10
--------
ENTHALPY
--------

The enthalpy is H = U + kB*T
                kB is Boltzmann's constant
Total free energy                 ...   -157.06201854 Eh
Thermal Enthalpy correction       ...      0.00094421 Eh       0.59 kcal/mol
-----------------------------------------------------------------------
Total Enthalpy                    ...   -157.06107433 Eh

熵（Entropy S）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
-------
ENTROPY
-------

The entropy contributions are T*S = T*(S(el)+S(vib)+S(rot)+S(trans))
     S(el)   - electronic entropy
     S(vib)  - vibrational entropy
     S(rot)  - rotational entropy
     S(trans)- translational entropy
The entropies will be listed as mutliplied by the temperature to get
units of energy

Electronic entropy                ...      0.00000000 Eh      0.00 kcal/mol
Vibrational entropy               ...      0.00394490 Eh      2.48 kcal/mol
Rotational entropy                ...      0.01155172 Eh      7.25 kcal/mol
Translational entropy             ...      0.01805279 Eh     11.33 kcal/mol
-----------------------------------------------------------------------
Final entropy term                ...      0.03354941 Eh     21.05 kcal/mol

吉布斯自由能（Gibbs free energy G）

1
2
3
4
5
6
7
8
9
10
11
12
13
-------------------
GIBBS FREE ENERGY
-------------------

The Gibbs free energy is G = H - T*S

Total enthalpy                    ...   -157.06107433 Eh
Total entropy correction          ...     -0.03354941 Eh    -21.05 kcal/mol
-----------------------------------------------------------------------
Final Gibbs free energy         ...   -157.09462374 Eh

For completeness - the Gibbs free energy minus the electronic energy
G-E(el)                           ...      0.08011291 Eh     50.27 kcal/mol

G-E(el)是在特定温度压力下，添加进电子能量项用于将Eel转换为G0
G-E(el)仅取决于几何结构和频率，能用低等级运算例如DFT，然后结合更高等级的方法计算Eel，例如DLPNO-CCSD(T)获得更好的预测结果

实验比较

DFT直接计算吉布斯自由能

1
2
!B3LYP DEF2-TZVP OPT FREQ
* XYZFILE 0 1 butene.xyz

DLPNO-CCSD(T)计算Eel

1
2
!DLPNO-CCSD(T) aug-cc-pVTZ AUTOAUX RIJCOSX
* XYZFILE 0 1 butene_optimized_DFT.xyz

将DFT算出的G-E(el)添加到更高级的算法结果中，形成更好的吉布斯自由能

不同温度的热力学

计算不同温度下的热力学性质

1
%FREQ TEMP 77, 298, 330, 450 END

热力学函数被不同的温度标题区分

1
2
3
4
5
6
7
8
--------------------------
THERMOCHEMISTRY AT 77K
--------------------------
[...]
--------------------------
THERMOCHEMISTRY AT 298K
--------------------------
[...]

在计算结束之后改变温度

如果在计算结束之后，还需进一步计算热力学函数，可以设置如下输入：

1
2
3
4
5
6
7
!PRINTTHERMOCHEM
%GEOM
         INHESSNAME "basename.hess"
END
%FREQ TEMP 290, 295, 300
END
* XYZFILE 0 1 geometry.xyz

可以从重新计算中节省大量时间

结构

gem-butadiene

1
2
3
4
5
6
7
8
9
10
11
12
C         -1.43336        0.99760       -0.11833
C         -0.17392        0.62436       -0.39435
C          0.14450       -0.70205       -1.02166
C          1.00254        1.50776       -0.09468
H          0.65348       -0.55631       -1.97995
H         -0.75651       -1.29570       -1.20886
H          0.79819       -1.28546       -0.36521
H         -2.27870        0.35250       -0.33839
H         -1.65494        1.95886        0.33562
H          1.55177        1.73378       -1.01439
H          0.69898        2.45896        0.35516
H          1.68283        1.00994        0.60389

trans-butadiene

1
2
3
4
5
6
7
8
9
10
11
12
C          1.30292        0.51517        0.05694
C         -0.04707       -0.10202        0.22754
H          1.76922        0.66709        1.03514
H          1.25468        1.48282       -0.45237
H          1.94420       -0.14853       -0.53097
C         -1.18840        0.46276       -0.19477
C         -2.53838       -0.15441       -0.02414
H         -3.00427       -0.30734       -1.00238
H         -3.17998        0.50979        0.56286
H         -2.49026       -1.12154        0.48617
H         -0.07366       -1.06762        0.72803
H         -1.16180        1.42834       -0.69529

cis-butadiene

1
2
3
4
5
6
7
8
9
10
11
12
C          1.48629       -0.48153        0.16062
C          0.03276       -0.80512        0.28503
H          1.68200        0.56889       -0.06324
H          1.92979       -1.08298       -0.63914
H          1.99954       -0.72748        1.09552
C         -1.01110        0.03009        0.16420
H         -0.18277       -1.85160        0.49931
C         -0.98494        1.49624       -0.12334
H         -2.00868       -0.39072        0.28765
H         -1.47672        2.03930        0.68986
H          0.02271        1.90131       -0.23449
H         -1.53322        1.69981       -1.04849

隐式解模型(inplicit solvation models)

隐式溶剂模型：在计算中包含溶剂效应的一种方法
主要讨论ORCA中两个主要的选项

模型解释：溶质被放置在一个大致分子形状的空腔中。溶剂由与空腔表面上电荷相互作用的连续体描述，而连续体又由溶质决定。
img

probe是理想的溶剂分子
分子空腔被定义为

  1. Van der Waals surface(VDW)，从vdw半径构建
  2. Solvent-Accessible Surface(SAS)，溶剂可及表面，由probe中心定义
  3. Solvent-Excluded Surface(SES)，无溶剂表面，沿着在 vdW 表面滚动的溶剂探针球体的内部部分获得

溶剂能由两部分组成，静电势electrostatic(ΔGENP)和腔色散cavitydispersion(ΔGCDS)
ΔG0solv=ΔGENP+ΔGCDS

Conductor-like Polarizable Continuum Model (CPCM)

在 CPCM 方法 中，本体溶剂被视为类似导体的可极化连续体，定义该方法的主要参数是介质的__折射率__和__介电常数__
由介质和分子表面电荷相互作用产生的静电贡献 (ΔGENP) 包含在 SCF 计算中 - 这样您甚至可以获得“溶剂化”轨道 - 并且腔项 (ΔGCDS) 可以从更复杂的方案中获得。
ORCA 有一个预定义溶剂列表，可以使用以下命令调用这些溶剂：!CPCM(solvent)

Solvent Name Dielec. Const.介电常数 Refrac. Index折射率
Water 80.4 1.33
Acetone 20.7 1.359
Acetonitrile 36.6 1.344
Ammonia 22.4 1.33
Benzene 2.28 1.501
CCl4 2.24 1.466
CH2Cl2 9.08 1.424
Chloroform 4.9 1.45
Cyclohexane 2.02 1.425
DMF 38.3 1.430
DMSO 47.2 1.479
Ethanol 24.3 1.361
Hexane 1.89 1.375
Methanol 32.63 1.329
Octanol 10.3 1.421
Pyridine 12.5 1.510
THF 7.25 1.407
Toluene 2.4 1.497

如果想包括其他溶剂，可以输入下列参数：

1
2
3
%CPCM       EPSILON      80.4
            REFRAC       1.33
END

例：使用单点能计算水中阿司匹林分子

1
2
!BP86 DEF2-SVP CPCM(WATER)
* XYZFILE 0 1 aspirin.xyz
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 阿司匹林
C         -2.64076        2.23326        0.00014
C         -3.28499        1.00146       -0.00001
C         -2.53276       -0.17323       -0.00006
C         -1.24586        2.29692        0.00023
C         -0.48687        1.12113        0.00013
C         -1.12591       -0.13504        0.00002
C         -0.44272       -1.46662       -0.00003
O         -1.06358       -2.51853       -0.00017
O          0.90426       -1.44815        0.00010
O          0.94078        1.15087        0.00026
C          1.78023        2.27870       -0.00045
O          1.43411        3.44948       -0.00130
C          3.21351        1.83099        0.00005
H          3.28673        0.73990       -0.00017
H          3.70751        2.20906        0.89847
H          3.70821        2.20936       -0.89786
H          1.28154       -0.55226        0.00027
H         -3.22416        3.15173        0.00019
H         -4.37149        0.95196       -0.00007
H         -0.79315        3.28276        0.00040
H         -3.05624       -1.12937       -0.00015

CPCM可以使用解析梯度以及Hessian矩阵，所以CPCM可以和OPT以及FREQ一起使用
在SCF之前，打印一个包含所有CPCM相关信息的标题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
--------------------
CPCM SOLVATION MODEL
--------------------
CPCM parameters:
  Epsilon                                         ...      80.4000
  Refrac                                          ...       1.3300
  Rsolv                                           ...       1.3000
  Surface type                                    ... GAUSSIAN VDW
  Epsilon function type                           ...         CPCM
Radii:
 Radius for C  used is    3.8550 Bohr (=   2.0400 Ang.)
 Radius for O  used is    3.4469 Bohr (=   1.8240 Ang.)
 Radius for H  used is    2.4944 Bohr (=   1.3200 Ang.)
Calculating surface                               ...        done! (  0.0s)
GEPOL surface points                              ...         1160
GEPOL Volume                                      ...    1403.2346
GEPOL Surface-area                                ...     757.5566
Calculating surface distance matrix               ...        done! (  0.0s)
Performing Cholesky decomposition & store         ...        done! (  0.1s)
Overall time for CPCM initialization              ...                 0.2s

在SCF计算结束之后，会有CPCM相关信息的汇总

1
2
3
4
5
6
7
8
9
10
11
Components:
Nuclear Repulsion  :          765.01084602 Eh           20817.00344 eV
Electronic Energy  :        -1413.26452682 Eh          -38456.88288 eV
One Electron Energy:        -2400.58284083 Eh          -65323.18007 eV
Two Electron Energy:          987.31831400 Eh           26866.29718 eV
CPCM Dielectric    :           -0.01874113 Eh              -0.50997 eV
[...]
CPCM Solvation Model Properties:
Surface-charge          :           -0.01502552
Charge-correction       :           -0.00002872 Eh              -0.00078 eV
Free-energy (cav+disp)  :  This term is not implemented in the current solvation scheme

ORCA默认不会为CPCM计算或添加空腔

将DFT和后HF结合

HF密度质量低于从DFT获得的质量，因此使用HF获得的CPCM校正也是如此
通常做法：从某些计算(如 DFT)中获取 ΔGENP（CPCM 电介质）和 ΔGCDS（自由能（cav+disp）），并将其添加到使用该方法获得的总能量中，例如如果您想将 CPCM 与 DLPNO-CSSD 结合使用。

Gaussian point charges

通常点电荷会导致能量不稳定，例如两点太近而崩溃。在ORCA，使用gaussian电荷来获得一个平滑的势能面。
在ORCA5.0中，默认设置为

1
%CPCM SURFACETYPE VDW_GAUSSIAN END

但是SES表面也可以被设置为：

1
%CPCM SURFACETYPE GEPOL_SES_GAUSSIAN END

Universal Solvation Model (SMD)全溶剂模型

SMD方法可被认为是CPCM的改进版
SMD使用完整的溶质电子密度来计算腔色散贡献，而不仅仅是区域
SMD需要更多的参数，这使得SMD对未知溶剂的灵活性兼容性较差，但是目前ORCA拥有179种溶剂可用。
使用方法：

1
2
3
%CPCM SMD TRUE
      SMDSOLVENT "SOLVENT"
END

例：

1
2
3
4
5
!BP86 DEF2-SVP
%CPCM SMD TRUE
      SMDSOLVENT "WATER"
END
* XYZFILE 0 1 aspirin.xyz

SMD最开始输出与CPCM相似

1
2
3
4
5
6
7
8
SMD-CDS solvent descriptors:
  Soln                                       ...    1.3328
  Soln25                                     ...    1.3323
  Sola                                       ...    0.0000
  Solb                                       ...    0.0000
  Solg                                       ...    0.0000
  Solc                                       ...    0.0000
  Solh                                       ...    0.0000

可以使用ΔGCDS (Gcds) 项和ΔGENP (CPCM Dielectric) 校正进一步计算的能量。

Example-辛醇/水分配系数

使用SMD模型，预测药物Diazepam在octanol/water中的分配系数
$$
P_{o/w}=\frac {[solute]{octanol}} { [solute]{water}}
$$
P可以从平衡常数以及溶剂在辛醇和水中的吉布斯自由能差异关系中得到 ΔGoo/w=ΔGoo−ΔGow
$$
logP_{o/w}=\frac{-\Delta G^o_{o/w}}{2.303RT}
$$
首先，优化在不同溶剂中的几何构型

1
2
3
4
5
6
7
8
9
10
11
12
!B3LYP DEF2-SVP OPT NUMFREQ D4
%CPCM SMD TRUE
      SMDSOLVENT "1-OCTANOL"
END
* XYZFILE 0 1 diazepam.xyz

or
!B3LYP DEF2-SVP OPT NUMFREQ D4
%CPCM SMD TRUE
      SMDSOLVENT "WATER"
END
* XYZFILE 0 1 diazepam.xyz

得到每个溶剂化化合物的吉布斯自由能及其几何形状
经过计算，得到3.41kcalmol-1的能量差，对应logP=2.5，实验值为2.99.

气液相转变

$$
\Delta G^o_{solv}=\Delta G {ENP}+\Delta G{CDS}+\Delta G^o_{conc}
$$

在1atm 298k时，从1molL-1气相转变为液相(需要计算G0，默认使用FREQ)
$$
\Delta G^o_{conc}=RTln(24.5)=1.89kcalmol^{-1}
$$
在预测溶液热力学时，不要忘记加G0.

CPCM 和COSMO

ORCA没有COSMO接口，但是可以使用CPCMC关键词来调用COSMO epsilon

1
!B3LYP DEF2-SVP CPCMC(WATER)

Correlation Energy

Correlation Energy被定义为高等级理论方法（MP2，CCSD）与参考HF方法之间的能量差
他是精确总能量的重要指标，对于预测特性或能量差有重要意义
在ORCA中，有几种方法来计算或者考虑Correlation Energy。
在此讨论三种

  1. Møller–Plesset perturbation theory (MP2);
  2. double-hybrid DFT (DHDF)
  3. coupled cluster (CC), using the singlet-triplet gap of methylene as a benchmark

讨论的范例，三重态比单重态更加稳定，HF方法误差极大，并且实验特性测量的非常准确

Møller–Plesset perturbation theory (MP2)

使用DHDF B2PLYP计算单重态和三重态的参考几何，能够准确的体现实验得到的几何坐标

1
2
3
4
5
6
!RI-B2PLYP DEF2-TZVPP DEF2-TZVPP/C OPT
* xyz 0 1
  C      0.000000    0.000000  -0.092116
  H     -0.905815    0.000000   0.548403
  H      0.905815    0.000000   0.548403
*

计算得出的三重态键角为134°，与实验相符
接着，使用MP2和先前针对单重态和三重态优化的结构计算绝热单重态-三重态之间的能量差

1
2
3
4
5
6
7
# 单重态
!MP2 DEF2-QZVPP
* xyz 0 1
  C      0.000000    0.000000   -0.129289
  H     -0.859359    0.000000    0.566990
  H      0.859359    0.000000    0.566990
*
1
2
3
4
5
6
7
# 三重态
!MP2 DEF2-QZVPP
* xyz 0 3
  C     -0.000000    0.000000    0.056545
  H     -0.990707    0.000000    0.474072
  H      0.990707    0.000000    0.474072
*

相关方法的准确性很大程度上取决于基组的大小，原则上选取更大的基组作为目标。

SCF之后将显示MP2header

1
2
3
4
5
6
7
8
9
10
11
12
13
14
------------------------------------------------------------------------------
                                ORCA  MP2
------------------------------------------------------------------------------

Freezing NCore=2 chemical core electrons

----------
MP2 ENERGY (disk based algorithm)
----------

Dimension of the basis                    ...  117
Memory devoted to MP2                     ...  256 MB
Data format for buffers                   ... DOUBLE
Compression type for matrix containers    ... UNCOMPRESSED

积分详细信息

1
2
3
4
5
6
7
8
9
10
11
-------------------------------
PARTIAL EXCHANGE TRANSFORMATION
-------------------------------

Transformation type                        ... one-operator
Generation of integrals (i,mue|j,nue)      ... ON
Generation of integrals (mue,kappa|nue,tau)... OFF
Generation of integrals (i,mue|a,nue)      ... OFF
Dimension of the basis                     ...  117
Number of internal MOs                     ...    3 (   1-   3)
Pair cutoff                                ... 1.000e-11 Eh

最后输出单重态相关能量

1
2
3
-----------------------------------------------
MP2 CORRELATION ENERGY   :     -0.147755815 Eh
-----------------------------------------------

生成最终能量

1
2
3
-------------------------   --------------------
FINAL SINGLE POINT ENERGY       -39.043423495766
-------------------------   --------------------

单重态减去三重态能量得到 ΔEMP2=14.24kcal/mol，比实验值高5
ΔEHF=28.25kcal/mol

Spin-component-scaled MP2 (SCS-MP2)

自旋分量级MP2
另一个版本的MP2理论为SCS-MP2，总体上由于MP2，可以如下调用

1
2
3
4
5
!SCS-MP2 DEF2-QZVPP
* xyzFILE 0 1 CH2_optimized.xyz

!SCS-MP2 DEF2-QZVPP
* xyzFILE 0 3 CH2_optimized.xyz

ΔESCS−MP2=7.92kcal/mol，比MP2更好

使用 resolution of identity(RI)

MP2或SCS-Mp2计算可以使用库伦积分的恒等近似分辨率来加速

1
2
3
4
5
!RI-MP2 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

!RI-SCS-MP2 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

此例中，相关方法需要使用特殊的辅助基/C，可以随时使用AUTOAUX自动生成
任何近似值都存在固有误差，但误差相对较小。单重态RI-MP2最终能量为：

1
2
3
-------------------------   --------------------
FINAL SINGLE POINT ENERGY       -39.054461324414
-------------------------   --------------------

RIJDX 和 RIJCOSX 近似也可用于 SCF 部分，通过在主输入中使用 !RIJDX 或 !RIJCOSX 进一步加速计算。

DLPNO scheme

Domain-based Local Pair Natural Orbital (DLPNO): ORCA团队开发了相关方法的特殊近似
相关能量可以用电子对能量计算，并且这些能量只能在每对的某个域内有效计算，使计算更加高效，并恢复了大约99%的相关能量。

1
2
3
4
5
!DLPNO-MP2 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

!DLPNO-SCS-MP2 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

输出有额外的输出信息，与轨道定位有关

1
2
3
4
5
6
7
8
9
10
11
12
------------------------------------------------------------------------------
                           ORCA ORBITAL LOCALIZATION
------------------------------------------------------------------------------

Input orbitals are from                  ... carb.gbw
Output orbitals are to                   ... carb.loc
Max. number of iterations                ... 128
Localizations seeded randomly            ... off
Convergence tolerance                    ... 1.000e-06
Treshold for strong local MOs            ... 9.500e-01
Treshold for bond MOs                    ... 8.500e-01
[...]

MP2运算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
----------------------------
Local MP2 Energy Calculation
----------------------------
Truncation parameters:
    PAOOverlapThresh    =        1.000e-08
    TCutPNO             =        1.000e-08


    Spin component scaling                      .... not used
    MP2 energy scaling factor                   .... 1.000

    Formation of semicanonical amplitudes.
     10% done.
     20% done.
     30% done.
     40% done.
     [...]

相关能量输出，与常规MP2 99.99%相似

1
2
3
------------------------------------------------------
 DLPNO-MP2 CORRELATION ENERGY:      -0.158789484960 Eh
 ------------------------------------------------------

Calculations using DLPNO also always require a /C basis.

Double-hybrid density functionals (DHDF) 双混合密度泛涵

DFT already has a component of correlation energy to them, added by the exchange-correlation functional, but it was shown that by adding a certain amount of the MP2 correlation energy the results can be improved even further
这些具有MP2组件的泛函成为双杂化泛函，他们同时包括HF交换和MP2相互关系，他们可以被调用，同时使用RI- DLPNO-近似

1
2
3
4
5
6
7
8
!B2PLYP DEF2-QZVPP
* xyzFILE 0 1 CH2_optimized.xyz

!RI-B2PLYP DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

!DLPNO-B2PLYP DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

The output follows a regular DFT calculation, with an MP2 calculation being done after the SCF. Finally all energies are added according the functional definition and the final energy is printed. For RI-B2PLYP, the singlet-triplet gap for methylene is ΔEB2PLYP=11.73kcal/mol, better than the simple MP2.

Double-hybrid spin-component-scaled density functionals (DSD-DF)

Following the same logic used for SCS-MP2, the different spin components of the MP2 energy can be scaled for the DHDF, which usually results in even better predictions. These functionals have DSD- in front of their name, for instance:

1
2
3
4
5
6
7
8
!DSD-PBEB95 DEF2-QZVPP
* xyzFILE 0 1 CH2_optimized.xyz

!RI-DSD-PBEB95 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

!DLPNO-DSD-PBEB95 DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

所有的方法均采用解析梯度，对于几何优化十分有用

Coupled cluster methods (CC)

耦合聚类方法
量子力学中获得相关能量的最准确方法之一为Coupled cluster methods（CC），特别是单参考问题中能准确计算能量的CCSD(T)方法
然而，尺度会随着电子数量的增加而增加，并且很快会变得不现实，其应用范围限制在20个原子左右
在ORCA中，使用DLPNO方案降低尺度问题，现在可以计算具有数百个原子的真实大小的分子

1
2
!DLPNO-CCSD DEF2-QZVPP DEF2-QZVPP/C
* xyzFILE 0 1 CH2_optimized.xyz

after the regular HF output, the Matrix Driven CI module starts (MDCI)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
-------------------------------------------------------------------------------
                              ORCA-MATRIX DRIVEN CI
-------------------------------------------------------------------------------


Wavefunction type
-----------------
Correlation treatment                      ... CCSD
Single excitations                         ... ON
Orbital optimization                       ... OFF
Calculation of Z vector                    ... OFF
Calculation of Brueckner orbitals          ... OFF
Perturbative triple excitations            ... OFF
Calculation of F12 correction              ... OFF
Frozen core treatment                      ... chemical core (2 el)
Reference Wavefunction                     ... RHF
  Internal Orbitals:     1 ...    3 (  3 MO's/  6 electrons)
  Virtual  Orbitals:     4 ...  116 (113 MO's              )
Number of AO's                             ...    117
Number of electrons                        ...      8
Number of correlated electrons             ...      6

接下来输出位置、筛选、积分变换，CCSD计算

1
2
3
4
5
6
7
8
9
10
11
12
------------------------------------------------
                  RHF COUPLED CLUSTER ITERATIONS
------------------------------------------------

Number of PNO amplitudes to be optimized   ...        22277
Number of non-PNO amplitudes               ...        76614
Untruncated number of regular amplitudes   ...        76614

Iter       E(tot)           E(Corr)          Delta-E          Residual     Time
  0    -39.042828320     -0.146989446      0.000167865      0.017754515    4.35
                           *** Turning on DIIS ***
  1    -39.056621207     -0.160782333     -0.013792887      0.014559282    4.14

当收敛之后，汇总输出不同的能量组成以及诊断信息

1
2
3
4
5
6
7
8
9
10
11
----------------------
COUPLED CLUSTER ENERGY
----------------------

E(0)                                       ...    -38.895667683
E(CORR)(strong-pairs)                      ...     -0.170828218
E(CORR)(weak-pairs)                        ...     -0.000171190
E(CORR)(corrected)                         ...     -0.170999408
E(TOT)                                     ...    -39.066667092
Singles Norm <S|S>**1/2                    ...      0.022256446
T1 diagnostic                              ...      0.009086156

注意检查T1 diagnostic信息，T1诊断信息值大于0.02，结果不可信，并且HF参考会变差。

The predicted energy difference here is ΔECCSD=10.54kcal/mol，better than the MP2, but still with some error. If one runs the same calculation using the even more accurate DLPNO-CCSD(T) instead, that difference goes down to ΔECCSD(T)=9.70kcal/mol, which is almost exactly like the experimental value!

Extrapolation to the CBS limit 外推到CBS限制

As mentioned before in the MP2 section, these methods require large basis set to be really accurate, and their true power show up only when close to the basis set (CBS) limit. Even so, one does not need to use an immense basis to get there, it is also possible to extrapolate the CBS from lower basis set calculations
之前在 MP2 部分中提到的，这些方法需要大的基组才能真正准确，并且只有在接近基组 (CBS) 限制时才会显示它们的真实功效。即便如此，也不需要使用巨大的基来达到目标，也可以从较低的基组计算中推断出 CBS
计算接近CBS的DLPNO-CCSD(T)能量

1
2
!DLPNO-CCSD(T) EXTRAPOLATE(3) cc-pVQZ/C
* xyzFILE 0 1 CH2_optimized.xyz

默认不自动/C基组，选择大的基组使合适的。
Going back to methylene triplet-singlet gap, the extrapolated energy difference is ΔECCSD(T)−CBS=9.26kcal/mol, which is actually within the experimental error.

总结计算结果

Method ΔETS Error
HF 28.25 19.13
B3LYP 11.74 2.62
MP2 14.24 5.11
SCS-MP2 7.92 -1.22
RI-B2PLYP 11.73 2.60
RI-DSD-PBEB95 11.63 2.50
DLPNO-CCSD 10.54 1.42
DLPNO-CCSD(T) 9.70 0.58
DLPNO-CCSD(T)-CBS 9.26 0.14
Experiment 9.12 ± 0.2

Relativistic corrections 相对论校正

当计算包含重原子系统的性质或者几何构型时，相对论效应有很大影响，并且不能忽略。
两种方法相对论校正，零阶正则近似(Zero-Order Regular Approximation ZORA)和Douglas-Kroll-Hess (DKH) Hamiltonian。
他们有些许不同，并且具有特定的优势，必须仔细考虑结果，使用最广泛的是ZORA和其改进型

Including relativistic corrections

使用输入：

1
2
3
4
5
!B3LYP ZORA ZORA-DEF2-TZVP

or

!B3LYP DKH DKH-DEF2-TZVP

Note that here we are using specific basis for each method, named ZORA- or DKH-DEF2-TZVP. This is necessary because these basis have been specifically designed for these all-electrons calculations, and the relativistic correction should NOT be used together with the regular basis or pseudopotentials.

RI and ZORA/DKH

使用RI方法加速SCF，另一个特殊基组/J被使用

1
2
3
4
5
!B3LYP ZORA ZORA-DEF2-TZVP RIJDX SARC/J

or

!B3LYP DKH DKH-DEF2-TZVP RIJDX SARC/J

For instance, the SARC/J auxiliary basis can be used for all the ZORA or DKH-DEF2 basis. If no specific basis is available, then one can always use AUTOAUX to automatically generate one.

The SARC/J basis were optimized for RI on the SCF part, not the MP2 or higher-level correlated methods! For correlation specific /C basis consult the manual and in case of abscence use !AUTOAUX.

Example : Hg dimer

Let’s test the impact of these effects on the Hg dimer, that has an experimental bond length of 3.69Å. This a rather extreme case of a very heavy element in a homodiatomic molecule, however it highlights the importance and magnitude of these effects.
The geometry can be optimized at a regular Double-hybrid density functionals (DHDF) level using the DEF2-TZVP basis, that makes use of pseudopotentials for the core orbitals that try to simulate relativistic effects or using the all-electron ZORA-DEF2-TZVP basis:

1
2
3
4
5
6
7
8
9
10
11
12
13
!B2PLYP DEF2-TZVP DEF2-TZVP/C OPT
* XYZ 0 1
Hg 0 0 0
Hg 0 0 3
*

or

!B2PLYP ZORA SARC-ZORA-TZVP SARC/J AUTOAUX OPT
* XYZ 0 1
Hg 0 0 0
Hg 0 0 3
*

The first results in a molecule with a 3.84Å bond, and the relativistic ZORA results in 3.78Å , which is in better agreement with the experimental results!
The auxiliary basis for the RIJ approximation used during the relativistic case here was chosen as the appropriate SARC/J.
Geometry optimizations using relativistic corrections turn on by default a one-center approximation by default, that changes the energy values. Do not compare single point energies from those you obtain from an !OPT run, these numbers are be incompatible.