在双金属中心化合物中，若近似认为单电子局域在A、B两个中心，则该体系的磁性质可以用Heisenberg-Dirac-van Vleck自旋哈密顿来描述：
其中，J为交换耦合常数，文献中有多种计算方式，其中一种为
上式中，HS表示高自旋态（high-spin），BS表示对称破缺态（broken-symmetry），对应UDFT计算中的低自旋状态。由于<S**2>_HS大于<S**2>_BS，因此，若J为正值，则高自旋态的能量低，称为铁磁耦合；反之，低自旋态的能量更低，称为反铁磁耦合。由上式可知，要计算磁交换耦合常数，需要计算体系的高自旋态和对称破缺态。
本文介绍如何用ORCA来计算耦合常数。在ORCA中，可以在一次任务中完成两个计算，并且程序会自动给出耦合常数的数值；而使用Gaussian时，则需要手动计算两个态。关于Gaussian中的对称破缺计算，可参考以下文章：
- 谈谈Gaussian软件中的guess=mix
- 片段组合波函数实例1.双原子分子
下面以bisphenolate-bridged dicopper(II)配合物为例，展示如何用ORCA计算磁交换耦合常数，所用版本为ORCA 5.0.3。该化合物来自J. Phys. Chem. A, 2009, 113, 6751一文，该文中作者用Jaguar程序计算了20多个化合物的交换耦合常数。【注：该文中的公式(2)（即上文计算J的公式）有误，分母上应该是S2的期望值，而非S期望值的平方。】该分子的结构如下图所示：
首先在高自旋（自旋多重度为3）下进行结构优化，所用级别为B3LYP/def2-SVP，用Gaussian 16完成，此处略去输入文件。然后在ORCA中做对称破缺计算，输入文件如下：
代码语言：javascript复制%maxcore 2000
%pal nprocs 52 end
! B3LYP def2-TZVP
%scf
brokensym 1,1
STABPerform true
STABRestartUHFifUnstable true
end
*xyz 0 3
C -3.419045 3.890224 -0.896573
C -2.104897 4.299474 -0.653920
C -1.146701 3.373287 -0.238256
C -1.505122 2.032186 -0.047223
C -2.856474 1.608591 -0.235467
C -3.789398 2.557507 -0.695709
H -4.163474 4.604448 -1.256267
H -1.816838 5.342584 -0.808902
H -0.104161 3.651312 -0.076621
H -4.806705 2.260928 -0.936038
C -3.273567 -2.752313 0.231133
C -4.306284 -1.817816 0.337076
C -4.215422 -0.399466 0.235589
H -5.305858 -2.219693 0.500842
O -2.025360 -2.479533 0.107186
O -0.613580 1.096361 0.302626
Cu -1.339498 -0.705432 0.154426
N -3.072796 0.227475 -0.010568
C -3.591605 -4.229498 0.274193
H -4.663995 -4.430277 0.397392
H -3.240335 -4.702811 -0.657326
H -3.034844 -4.700651 1.100503
C -5.519322 0.343278 0.448761
H -5.361933 1.283257 0.995617
H -6.003977 0.590388 -0.509917
H -6.217149 -0.286483 1.015989
C 3.419103 -3.890260 -0.896490
C 2.104932 -4.299494 -0.653938
C 1.146698 -3.373283 -0.238416
C 1.505099 -2.032169 -0.047433
C 2.856478 -1.608593 -0.235545
C 3.789444 -2.557537 -0.695648
H 4.163566 -4.604507 -1.256068
H 1.816889 -5.342615 -0.808876
H 0.104153 -3.651315 -0.076821
H 4.806781 -2.260979 -0.935877
C 3.273578 2.752312 0.231112
C 4.306270 1.817799 0.337156
C 4.215401 0.399450 0.235651
H 5.305829 2.219660 0.501053
O 2.025383 2.479550 0.107007
O 0.613522 -1.096321 0.302263
Cu 1.339493 0.705455 0.154128
N 3.072797 -0.227479 -0.010642
C 3.591630 4.229492 0.274241
H 4.664020 4.430248 0.397477
H 3.240396 4.702847 -0.657270
H 3.034865 4.700631 1.100560
C 5.519264 -0.343309 0.449000
H 5.361792 -1.283285 0.995836
H 6.004048 -0.590423 -0.509611
H 6.217020 0.286447 1.016322
*在%scf部分，用brokensym NA,NB来开启对称破缺计算，其中NA和NB表示“A”、“B”两个位点的单电子数目。本例中两个金属为二价铜，显然NA和NB均为1。STABPerform true表示检查波函数的稳定性。如果在输出文件中看到波函数不稳定，则该计算结果不可信。若碰到这种情况，我们可以加上STABRestartUHFifUnstable true关键词，表示若遇到不稳定波函数，则重新做一次SCF计算。需要说明的是，ORCA中的这个功能只对波函数优化一次，并不保证能得到稳定的波函数，而Gaussian中的stable=opt则会反复优化波函数直至稳定。
在输出文件中，会显示如下信息：
代码语言：javascript复制The Broken Symmetry feature is operative
Number of electrons in the 'A' domain ... 1
Number of electrons in the 'B' domain ... 1
We will first do a calculation with S= 1.0
and then try to converge to the broken symmetry state with Ms= 0.0也就是说程序会先做S=1的非限制性计算，再做对称破缺计算。因此在输出文件中我们也可以看到两处SCF迭代过程。由于本例中我们加入了检查波函数的稳定性，程序会在两步SCF后面均进行波函数稳定性分析。
为确保结果的可靠，还应检查自旋布居，即Spin population。ORCA中会给出用Mulliken和Loewdin方法得到的自旋布居。本例中，两个铜原子的Mulliken自旋布居分别为0.5946和−0.5946，接近1，可认为结果合理。
在程序的最后，会给出耦合常数的结果：
代码语言：javascript复制------------------------------------------
BROKEN SYMMETRY MAGNETIC COUPLING ANALYSIS
------------------------------------------
S(High-Spin) = 1.0
<S**2>(High-Spin) = 2.0058
<S**2>(BrokenSym) = 0.9816
E(High-Spin) = -4542.777216 Eh
E(BrokenSym) = -4542.778211 Eh
E(High-Spin)-E(BrokenSym)= 0.0271 eV 218.386 cm**-1 (ANTIFERROMAGNETIC coupling)
---------------------------------------------------------
| Spin-Hamiltonian Analysis based on H(HDvV)= -2J*SA*SB |
------- -----------
| J(1) = -218.39 cm**-1 (from -(E[HS]-E[BS])/Smax**2) |
| J(2) = -109.19 cm**-1 (from -(E[HS]-E[BS])/(Smax*(Smax 1)) |
| J(3) = -213.21 cm**-1 (from -(E[HS]-E[BS])/(<S**2>HS-<S**2>BS)) |
-------------------------------------------------------------------------程序给出了三种形式的J值，计算公式也写在数值之后，其中J(3)为本文开头的公式，此处，我们取该值。该体系的实验值为−298 cm−1，计算结果与之基本吻合。
磁交换耦合常数的值与几何结构和所选用的泛函有较大关系，因此在计算中要适当选取计算级别。一般来说，纯泛函的结果不够好，需要使用杂化泛函。感兴趣的读者可以用不同的泛函来尝试本文中的例子，对ORCA不熟的读者也可以用Gaussian来尝试。需要注意的是，由于ORCA 5.0版本开始，在SCF计算中默认开启RIJCOSX，且B3LYP泛函在Gaussian和ORCA中的实现不完全相同，因此结果会有微小差别。
致谢：感谢jxzou在成文中的建议。
