清晰度函数判据在晶体结构解析中的应用

陈志红++李晖

摘 要：晶体结构解析的核心是相位恢复问题，正确的相位对应正确的结构。与通常使用残差作为评判标准不同，本文使用实空间中基于电荷密度的清晰度函数作为判据，结合密度调整迭代方法，构建出一套完整的晶体结构解析流程。通过测试，该方法简便高效，计算过程中不需要调整任何参数，对使用人员的晶体学知识依赖较少。该算法的实现为未来实现晶体结构解析自动化奠定了基础。

关键词：晶体结构解析；清晰度函数；DM迭代算法

中图分类号：O645 文献标识码：A 文章编号：1671-2064（2017）13-0237-02

在晶体学发展的过程中，晶体结构的解析问题一直处于核心地位。在X-ray单晶衍射实验中，单晶样品的衍射信号被探测器接收，并通过计算机处理，获得晶体在实空间的“图像（电荷密度图）”。晶体结构解析的实质就是通过衍射实验获得后焦面的衍射信号，并对此信号做傅里叶变换（Fourier transform，FT）得到实空间中样品图像的过程。然而现有的探测器只能记录衍射点的强度（振幅），无法记录相位信息。相比于振幅，相位对于结构解析更为重要[1]。为了获得正确的结构，必须首先恢复衍射點的相位，这就是晶体结构解析中的“相位问题”，也可以说晶体结构解析问题等价于相位恢复问题。在过去一个世纪的发展历程中，晶体学家发展了许多优秀的晶体结构解析方法，如：早期的试错法、Patterson法、直接法（包括重原子法）、最大熵最大似然函数法、分子替代法（包括硒代法）、DM法（Density Modification Method）、多波长反常散射法；利用现代优化理论，如遗传算法，蒙特卡洛，模拟退火等发展起来的现代试错算法；以及最近引起极大关注的Charge Flipping算法等等[2，3]。对于不同种类的单晶样品，实验数据的质量，需要慎重挑选不同的方法进行尝试，并调整相关的参数。

另一个重要的问题是，如何在众多相位组合中判断哪一个才是正确的，或者说如何来判别哪一个结果是正确的。最常用的判据是倒易空间中的统计，理论计算与实验数据相比较，如R值，ωR2值，GooF值等[4]。实空间的物理量是电荷密度，是否有基于电荷密度的判断依据呢？在之前的工作中[5]，我们从图像处理的角度，借鉴图像中清晰度的概念，指出正确的相位对应的图像清晰度最高，并构建出两个评价晶体结构（电荷密度）“清晰度”的函数——Tian1函数。该清晰度函数依据的是Fourier变换得到的电荷密度的自身特征，无需设定参考值，本文将结合密度调整迭代算法，对这一方法进行测试。

1 清晰度函数算法

实验收集到的单晶衍射强度与结构因子的平方呈正比关系：

（1）

结构因子与电荷密度互为Fourier关系，

（2）

同时，结构因子是一个复数，常可以写成极坐标的形式：

（3）

从衍射实验中我们只能得到其振幅，而无法记录衍射的相位信息。因此，实空间图像并不能直接通过Fourier变换获取，这就是所谓的晶体结构解析的“相位问题”。晶体结构解析问题就是在已知衍射强度条件下，如何恢复相位角，以及如何确定正确的相位的问题。首先讨论后一个问题，在一系列相位组合中，根据式每一组相位都可以计算一组理论结构因子，常用的判据是与实验观测结构因子相比较，如残差因子定义为：

（4）

其中是对应于实验观测结构因子的实验标准偏差，是计算的结构因子，M是衍射点的个数。R因子定义为：

（5）

ωR2定义为：

（6）

GooF定义为：

（7）

其中，是理论计算衍射强度，是理论结构因子的模方。结构因子和衍射强度均是倒易空间中的物理量，可以通过X-ray衍射实验获得，即公式中的。

晶体结构在实空间中的图像是电荷密度分布，类比图像处理中，衡量一张图像的质量的标准是其清晰度，图像越清晰，说明越符合真实场景。相类似，我们提出了晶体清晰度的概念，或者说电荷密度的清晰度概念。正确的电荷密度清晰度最高。依据Fourier变换后的电荷密度的性质，我们构建出描述电荷密度的清晰度的函数：

（8）

公式中为第i个电荷密度大于零的格点。为第j个电荷密度小于零的格点。其中m和n分别是正电荷密度格点总数与负电荷密度格点总数。它们与电荷密度为零的网格数加和等于实空间网格总数。

下面是基于T1函数的迭代算法，基本流程如下：1）随机产生一套初始相位；2）与实验观测结构因子结合，Fourier变换计算电荷密度，以及T1函数，判断是否达到收敛；3）如果是，结束计算；如果否，保持电荷密度值为负的格点不变，将为正值的扩大倍，n为一个足够大的正数，例如100；第五步，将上述新的电荷密度做逆Fourier变换，得到新结构因子及新的相位；第六步，将新得到的相位与观测的结构因子的模相结合，即返回第二步，进入下一轮迭代结算过程。判断是否收敛这里使用是是比较相邻两轮T1值的差小于一个设定的阈值（如0.0001）。整个流程图如图1所示。

2 计算结果

为了验证上面所提出的解决方案有效性，我们以NaB3FO5予以证明，单晶衍射数据分辨率为0.8。其晶胞参数为a=6.6888、b=4.60913、c=4.00587、α=90.0、β=113.65、γ=90.0°，空间群设定为P1。基于上述算符我们编写了一套程序，计算时只需要输入晶胞参数和衍射数据文件（.hkl），计算过程中无需人为参与，结果输出电荷密度文件（.rho）用以显示晶体结构。

图2是T1函数的迭代计算图，经过200轮迭代后，T1函数值达到收敛，最终T1值为9.55057。每一轮迭代只需要做两次Fourier变换，使用快速Fourier变换（FFT），因此计算效率非常高。本实例计算过程在普通台式计算机上用时小于5分钟。

为了证明该清晰度函数方法给出的结构是正确的，我们也用SHELX[6]解析该结构，残差因子R=0.0549，wR2=0.1626，GooF=1.201。结果如图3所示，图3（a）为SHELX给出的结构模型，图中红色的是O原子，绿色的是Be原子，黄色的是Na原子，灰色的是F原子；图3（b）为T1函数得到的电荷密度图。从结构上看，两者符合的很好。

3 结语

本文讨论了清晰度函数法在晶体解析中的应用，以NaB3FO5为例进行结构解析，与SHEXL程序给出的结构符合得很好。該方法简单高效，无需使用人员拥有较多的晶体学知识，需调节的参数很少。

我们相信，通过数字图像处理角度来研究晶体结构解析方法将为晶体结构解析方法研究开拓了新的研究思路。

参考文献

[1]Aubert E， Lecomte C. Illustrated Fourier transforms for crystallography[J]. Journal of Applied Crystallography， 2007， 40（6）： 1153-1165.

[2]Hauptman H A. The phase problem of x-ray crystallography[J]. Reports on Progress in Physics， 1991， 54（11）： 1427.

[3]Oszlányi G， Süt A. Ab initio structure solution by charge flipping[J]. Acta Crystallographica Section A： Foundations of Crystallography， 2004， 60（2）： 134-141.

[4]Wall M E. Quantum crystallographic charge density of urea[J]. IUCrJ， 2016， 3（4）.

[5]Li H， He M， Zhang Z. Image definition evaluation functions for X-ray crystallography： a new perspective on the phase problem[J]. Acta Crystallographica Section A： Foundations and Advances， 2015， 71（5）： 526-533.

[6]Sheldrick G M. A short history of SHELX[J]. Acta Crystallographica Section A： Foundations of Crystallography，2008， 64（1）：112-122.