Applications
穆斯堡尔谱
1958 年，R. L. Mössbauer在研究γ射线共振吸收现象时，发现了 穆斯堡尔效应 。 同质异能位移（isomer shift；\(\delta^{IS}\)）是穆斯堡尔谱的重要观测参量之一，它源于具有一定尺寸的原 子核与周围电子分布之间的库仑相互作用。当原子处在不同的外界环境，导致原子核附近的库仑能发生变化， 而 \(\delta^{IS}\) 对这种变化极其敏感，因此可以用来研究原子的氧化态，自旋态，以及配位环境。 穆斯堡尔谱的另一个重要观测参量是四级分裂，来源于原子核的电四极矩与原子核周围电场梯度（electric field gradient；EFG）之间的相互作用。 此外，当原子置于外磁场中，由于原子核存在核磁矩，穆斯堡尔谱会进一步发生 Zeeman 分裂。
\(\delta^{IS}\) 可以表示为重元素在测试体系A和参照体系R中“有效接触密度”（effective contact density；ED） 或“接触密度”（contact density；CD）变化量的线性函数：
以上两个公式称为 校准方程 （calibration equation），其中 \(\delta^{IS}\) 是测试体系A相对于参照系R的同质异能位移实验值， ED或CD值 \(\rho_{A}\) 和 \(\rho_{R}\) 可以通过理论计算获得， \(\alpha\) 和 \(\beta\) 是待拟合的参数，其中 \(\alpha\) 也称为核标定常数（nuclear calibration constant）， \(C\) 是任意量，一般取ED或CD的整数部分。考虑到 \(\rho_{R}\) 理论值存在误差，一般用后一个公式进行拟合。
Note
CD假定原子核附近的电子密度是均匀分布的，因此可以用一个点的电子密度代替（一般取原子核中心位置的密度， 但也有程序取一系列样点的密度做加权平均）；ED考虑了电子密度的非均匀分布，原则上比前者更合理。 很多程序计算的是CD，而在BDF中计算ED，二者近似满足换算关系（参见文献 [73] 的Table S2、S3）， 换算因子可以吸收到 \(\alpha\) 中。
为了准确计算ED及其相对变化量，需要考虑以下两个因素：
原子核具有一定的尺寸分布，而默认的点电荷近似可能会导致几个数量级的误差！为此需要在
xuanyuan模块设定nuclear=1。需要考虑相对论效应。对于重元素这是显而易见的；即便对于某些轻元素也必须要考虑相对论效应， 这是因为非相对论情况下的 p 电子在原子核附近没有分布（参见文献 [73] 的Table S6）， 从而对 p 区元素的ED造成定性错误。在BDF中可以用sf-X2C哈密顿或其局域变体考虑标量相对论效应， 通过
xuanyuan模块中的heff=21（标准sf-X2C），22（sf-X2C-AXR），或23（sf-X2C-AU）指认。
铁（ \(\ce{^{57}Fe}\) ）化合物的有效接触密度
γ射线吸收和发射的几率正比于 \(\exp(-E_\gamma^2)\) ，当核激发能 \(E_\gamma\) 超过200 keV后，一般很难观测到穆斯堡尔谱， 由此只能对少数同位素探测穆斯堡尔谱。 \(\ce{^{57}Fe}\) 就属于适合实验测量的同位素之一，不过在理论计算中，一般不对这些同位素进行区分。
计算ED需要非常陡峭（也就是高斯指数非常大）的 s 型高斯基函数才能准确描述电子在铁原子核附近的分布；
对于存在 p 价电子的 p 区元素，还需要非常陡峭的 p 型高斯基函数，而基组库中的标准收缩基组通常不符合要求。
建议采用cc-pVnZ型或ANO型全电子相对论基组，并把其中的 s 函数（ p 区元素还有 p 函数）进行非收缩处理。
在以下的全电子相对论计算中，铁采用ANO-R2基组（具有三ζ精度），
并把 s 函数做非收缩处理，也就是删除 s 函数的收缩因子部分，并把收缩度从6改为0，然后存为其它文件名（如ANO-R2-ED）。
由于铁没有4 p 价电子， p 函数不需要修改。
把ANO-R2-ED放到执行计算的目录下，供后面的计算调用（下载链接 ano-r2-ed.zip ）。
我们对于铁的一系列模型体系化合物进行相对论密度泛函理论计算，泛函选取PBE0，相对论哈密顿用sf-X2C-AU。 除了铁以外的轻元素全部用def2-TZVPP基组，它在Kr元素之前属于全电子基组，虽然是非相对论的，但用于前18号元素的相对论计算是允许的。 自旋多重度和分子坐标来自文献 [74] 。以 \(\ce{FeF6^{4-}}\) 为例，输入如下：
$compass
title
FeF_6^4-
basis-block
def2-tzvpp
Fe = ANO-R2-ED
end basis
geometry # 分子直角坐标，单位：埃
Fe -0.000035 0.000012 0.000014
F 2.116808 -0.003546 0.032360
F -2.116824 0.001611 -0.030945
F -0.003602 2.164955 0.001902
F 0.001648 -2.165219 -0.003295
F 0.032586 0.003638 2.109790
F -0.030580 -0.001452 -2.109825
end geometry
MPEC+cosx # 使用MPEC+COSX加速
$end
$xuanyuan
heff # sf-X2C-AU；计算ED必须选21-23中的一个
23
nuclear # 高斯有限核模型；ED必须设为1
1
$end
$scf
charge
-4
spinmulti
5
uks
dft functional
pbe0
grid # DFT计算ED需要用精密格点
sg1
reled
26 # 只计算Fe的ED（对于本例，10至26的整数等价）
$end
计算完成后，在SCF布居分析信息之后可以找到ED结果：
Relativistic effective contact densities for the atoms with Za > 25
----------------------------------------------------------------
No. Iatm Za RMS (fm) Rho (a.u.)
----------------------------------------------------------------
1 1 26 3.76842 14552.68329
----------------------------------------------------------------
以此为例，完成其它铁化合物分子的ED计算（输入文件下载链接 ed-fe.zip ）。
ED结果以及 \(\delta^{IS}\) 实验值 [74] 列于下表：
分子 |
2S+1 |
\(\delta^{IS}\) (mm/s) |
ED ( \(bohr^{-3}\) ) |
|---|---|---|---|
\(\ce{FeCl4^{2-}}\)
\(\ce{Fe(CN)6^{4-}}\)
\(\ce{FeF6^{4-}}\)
\(\ce{FeCl4^-}\)
\(\ce{Fe(CN)6^{3-}}\)
\(\ce{FeF6^{3-}}\)
\(\ce{Fe(H2O)6^{3+}}\)
\(\ce{FeO4^{2-}}\)
\(\ce{Fe(CO)5}\)
|
5
1
5
6
2
6
6
3
1
|
+0.90
-0.02
+1.34
+0.19
-0.13
+0.48
+0.51
-0.87
-0.18
|
14551.76
14555.78
14552.68
14553.98
14556.08
14553.01
14554.12
14558.17
14556.37
|
用这些数据进行拟合，得到校准方程
可见拟合误差比较大，这可能是以下原因造成的：
样本太少
穆斯堡尔谱是对固态的真实体系测量的，与计算所用的气态离子模型不一致。用团簇模型、溶剂化模型 [75] 、嵌入模型 [76] 可能更合适。
铁的某些化合物存在强关联，需要测试其它泛函，或者换成适合描述强关联体系的方法
有了校准方程后，就可以对一些铁的体系预测 \(\delta^{IS}\) 。例如交错状的二环戊二烯基铁 [77] ， 通过以上密度泛函理论计算得到ED为14554.25 a.u.，代入校准方程得到 \(\delta^{IS}\) 为0.37 mm/s， 与实验值0.53 mm/s [77] 基本接近。
计算重元素化合物有效接触密度的注意事项
对于4d以上的元素，经验表明默认的高斯指数还不足以描述原子核附近的电子分布，需要额外补充一些更陡峭的高斯指数。 例如，选择cc-pVnZ型或ANO型标准基组中最陡峭的4-6个 s 型高斯指数α（ p 区重元素还要考虑 p 型高斯指数），它们近似满足以下线性关系：
通过线性拟合得到参数A、B，再通过外推（i的间隔取-0.5或-1），即可得到更陡峭的高斯指数。 一般加入2-5个更陡峭的 s 函数、1-3个更陡峭的 p 函数即可满足要求，但是要避免用10 11 以上的高斯指数， 因为这可能会造成数值不稳定。