穆斯堡尔谱
1958 年，R. L. Mössbauer在研究γ射线共振吸收现象时，发现了 穆斯堡尔效应 。 同质异能位移（isomer shift；\(\delta^{IS}\)）是穆斯堡尔谱的重要观测参量之一，它源于具有一定尺寸的原 子核与周围电子分布之间的库仑相互作用。当原子处在不同的外界环境，导致原子核附近的库仑能发生变化， 而 \(\delta^{IS}\) 对这种变化极其敏感，因此可以用来研究原子的氧化态，自旋态，以及配位环境。 穆斯堡尔谱的另一个重要观测参量是核四级分裂（nuclear quadrupole splitting；\(\Delta E_{Q}\)）， 它来源于核自旋量子数I > 1/2的原子核产生的电四极矩（eQ；这里e是一个质子电荷，Q是核四极矩NQM；nuclear quadrupole moment） 与原子核周围电场梯度（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}\) 理论值存在误差，一般用后一个公式进行拟合。
备注
CD假定原子核附近的电子密度是均匀分布的，因此可以用一个点的电子密度代替（一般取原子核中心位置的密度， 但也有程序取一系列样点的密度做加权平均）；ED考虑了电子密度的非均匀分布，原则上比前者更合理。 很多程序计算的是CD，而在BDF中计算ED，二者近似满足换算关系（参见文献 [95] 的Table S2、S3）， 换算因子可以吸收到 \(\alpha\) 中。
为了准确计算ED及其相对变化量，需要考虑以下两个因素：
原子核具有一定的尺寸分布，而默认的点电荷近似可能会导致几个数量级的误差！为此需要在
xuanyuan模块设定nuclear=1。需要考虑相对论效应。对于重元素这是显而易见的；即便对于某些轻元素也必须要考虑相对论效应， 这是因为非相对论情况下的 p 电子在原子核附近没有分布（参见文献 [95] 的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复制到$BDF_WORKDIR目录下，供后面的计算调用（下载链接 ano-r2-ed.zip ）。
注意：对于非标准基组，基组名的字母必须全部大写。
此外，也可以用我们专门为Fe元素优化的x2c-TZVPPall-CD基组，可以从基组库调用。
我们对于铁的一系列模型体系化合物进行相对论密度泛函理论计算，泛函选取PBE0，相对论哈密顿用sf-X2C-AU。 除了铁以外的轻元素全部用def2-TZVPP基组，它在Kr元素之前属于全电子基组，虽然是非相对论的，但用于前18号元素的相对论计算是允许的。 自旋多重度和分子坐标来自文献 [96] 。以 \(\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需要用精密格点
ultra fine
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.65555
----------------------------------------------------------------
以此为例，完成其它铁化合物分子的ED计算（输入文件下载链接 ed-fe.zip ）。
ED结果以及 \(\delta^{IS}\) 实验值 [96] 列于下表：
分子 |
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
|
用这些数据进行拟合，得到校准方程
可见拟合误差比较大，这可能是以下原因造成的：
样本太少
穆斯堡尔谱是对固态的真实体系测量的，与计算所用的气态离子模型不一致。用团簇模型、溶剂化模型 [97] 、嵌入模型 [98] 可能更合适。
铁的某些化合物存在强关联，需要测试其它泛函，或者换成适合描述强关联体系的方法
有了校准方程后，就可以对一些铁的体系预测 \(\delta^{IS}\) 。例如交错状的二环戊二烯基铁 [99] ， 通过以上密度泛函理论计算得到ED为14554.25 a.u.，代入校准方程得到 \(\delta^{IS}\) 为0.37 mm/s， 与实验值0.53 mm/s [99] 基本接近。
计算重元素化合物有效接触密度的注意事项
对于4d以上的元素，经验表明默认的高斯指数还不足以描述原子核附近的电子分布，需要额外补充一些更陡峭的高斯指数。 例如，选择cc-pVnZ型或ANO型标准基组中最陡峭的4-6个 s 型高斯指数α（ p 区重元素还要考虑 p 型高斯指数），它们近似满足以下线性关系：
通过线性拟合得到参数A、B，再通过外推（i的间隔取-0.5或-1），即可得到更陡峭的高斯指数。 一般加入2-5个更陡峭的 s 函数、1-3个更陡峭的 p 函数即可满足要求，但是要避免用10 11 以上的高斯指数， 因为这可能会造成数值不稳定。
铁（ \(\ce{^{57}Fe}\) ）化合物的EFG计算
EFG计算对相对论哈密顿的要求与有效接触密度的计算类似，但对基组的要求不同。
只有 s 电子以及少量 p 电子在原子核附近存在非零的分布，因此有效接触密度计算中只需要修改 s 、 p 基函数
原子核形变产生的电四极矩只能与轨道角动量 L > 0 电子的EFG发生相互作用，因此不必修改 s 基函数。 首先把 p 函数进行非收缩处理（删除 p 函数的收缩因子部分，并把收缩度改为0），并酌情添加1-2个陡峭的 p 型高斯函数。 对于存在 d 价轨道的过渡元素（镧系、锕系还有 f ），需要把 d （ f ）函数进行非收缩处理，并适当添加更陡峭的 d 、 f 函数。 我们已经为Fe元素的EFG计算专门优化了 x2c-TZVPPall-EFG 基组，可以从基组库调用。
若同时计算有效接触密度与EFG，对基函数的修改要满足以上两条要求。对于铁元素，可以用基组库中的 x2c-TZVPPall-CDEFG 基组。
EFG计算的关键词为 relefg 。
例如，同时计算有效接触密度与EFG，以上算例的 SCF 模块输入需要改为：
$scf
charge
-4
spinmulti
5
uks
dft functional
pbe0
grid # DFT计算EFG需要用精密格点
ultra fine
relefg
26 # 只计算Fe的EFG张量
reled
26 # 只计算Fe的ED
$end
计算完成后，在SCF布居分析信息以及ED结果之后，可以找到EFG张量的结果：
Relativistic electric field gradients for the atoms with Za > 25
-----------------------------------------------------------------------------
No. Iatm Za RMS (fm) EFG tensor (a.u.)
-----------------------------------------------------------------------------
1 1 26 3.76842 -0.1061 -0.0023 0.1850
-0.0023 0.0395 -0.0018
0.1850 -0.0018 0.0666
eta Vaa Vbb Vcc
0.64736 0.0395 0.1844 -0.2239
NQCC = -8.4172 MHz with Q(ISO-057) = 160.00 mbarn
-----------------------------------------------------------------------------
在EFG 张量的9 个分量中，6 个非对角元是行列对称的；3 个对角线元之和为零。如果选择一个特殊的坐标系 \(\{\vec{a},\vec{b},\vec{c}\}\) （即EFG 张量的主轴或特征矢量）， 使得非对角元为零，而对角元（即特征值）满足 \(|V_{aa}| \le |V_{bb}| \le |V_{cc}|\) ，此时EFG 张量只需要两个非独立参数来表示就可以了， 即主值 \(V_{cc}\) 和非对称参数 \(\eta = |(V_{aa} − V_{bb})/V_{cc}|\) （\(0 \le \eta \le 1\)）。当η = 0 时，EFG 张量为轴对称。 在本例中，η = 0.64736， \(V_{cc}\) = -0.2239 a.u. 。
注意
非阿贝尔群分子的简并态在计算EFG时，单个分支的 \(V_{cc}\) 和η 一般没有意义。必须对简并态的所有分支（通过在SCF中指定占据数）分别计算EFG张量， 对它们做平均后再计算 \(V_{cc}\) 和η。
对于孤立原子， \(V_{aa} = V_{bb} = V_{cc} = 0\) ；对于线形分子（包括双原子分子）， \(V_{cc} = V_{zz}\) （分子轴为z）。 利用这一特点，BDF可以对开壳层原子、线形分子简并态的EFG结果进行校正。
\(V_{cc}\) 的正负号是有实际意义的，但有两种例外情况。(1) 当 η 接近1时（例如 η > 3/4）， \(V_{bb}\) 、 \(V_{cc}\) 受计算误差的影响极易发生互换，导致与实验值相比差一负号。特别是当 η = 1 时， \(V_{cc}\) 的正负号没有实际意义。 (2) \(V_{cc}\) 很小时（例如在 0.25 a.u.以下），计算误差很可能导致三个特征值发生互换， \(V_{cc}\) 的正负号难以确定。
核四极矩与EFG 之间的相互作用通常用核四极耦合常数（NQCC；nuclear quadrupole coupling constant） \(eQq\) 来衡量（在一些文献中也写作 \(eqQ\) ），定义为
其中 \(V_{cc}\) 仍取原子单位，核四极矩Q的单位是Barn（1 Barn = 1.0e-28 平方米）， \(eQq\) 的单位是MHz。 当同位素的Q实验值或理论预测值已知时（对于 \(\ce{^{57}Fe}\) ，I=1/2 核基态的 Q 值为零，只有 I=3/2第一激发态的核四极矩可以与EFG发生相互作用， 该态的Q值为 0.160 barn），程序会打印 \(eQq\) ，在本例中是-8.4172 MHz。
穆斯堡尔谱测量的核四极分裂 \(\Delta E_{Q}\) 与NQCC满足一定的关系。例如，对于 \(\ce{^{57}Fe}\) 的I=1/2 → I=3/2核激发跃迁，有
其中， \(\Delta E_{Q}\) 的单位取 mm/s， \(E_\gamma = 14.412497\) KeV 为 γ 线的激发能， \(c / E_\gamma = 11.6248\) 为 MHz 到 mm/s 的换算因子。
对于其它核素，习惯上仍然选取 I=3/2 的核基态或第一激发态作为参考（如果存在的话），这样上面的 \(\Delta E_{Q}\) 公式依旧成立， 只不过 \(E_\gamma\) 要改为该核素第一激发能的实验值。
\(\Delta E_{Q}\) 以及之前ED的理论结果可以用于验证实验上确定的原子价态， 而高精度的 EFG 理论结果结合 \(eQq\) 的实验值可以用来确定核素的Q值。