经验赝势

aroman

5.00/5 (3投票s)

2021 年 7 月 21 日

GPL3

10分钟阅读

5488

这篇博文分享了一个计算具有金刚石/闪锌矿结构的晶体能带结构的工程项目，适用于各种元素。

引言

我很久没有在博客上写东西了。我非常忙，没有多少时间。不过，我还是做了一些相关的项目。其中一个项目，一个使用周期性边界条件方法（supercell method）计算分子的DFT程序，原本打算作为一篇博文的主题。我已经实现了它，它也能工作（大致如此），但还需要一些检查和清理。我使用了局部赝势，在实空间中指定，其中一个是我从参考文献中提到的Ge的一个生成出来的，还有一些是我从网上找到的。不幸的是，它的结果不如我预期的好。我可能想将其改为使用非局部赝势。

我还修改了 Hartree-Fock 项目，使其能够使用 Koopmans 定理计算和显示解离能和电离能，实际上就是最后一层能量。我这样做主要是因为我想将一些结果与DFT分子程序进行比较， Hartree-Fock 计算的解离能可能不会很好。我稍微修改的另一个项目是 Lattice Boltzmann 项目，我试图提高其速度。

我最终使用所有这些DFT项目想要达到的目标是将其用于晶体计算。所以，既然我已经无限期地推迟了那个项目，我想我至少应该有一个带有赝势的另一个项目。比DFT简单得多。但在讨论那个项目之前，这里有一个结果，这也是我在最后几篇博文中提到的讲座中呈现的。

Ge晶体 110 截面，使用DFT和赝势计算

这是使用Octave代码获得的，尽管我也可以用分子计算项目中的C++代码重现这些结果。

现在，回到这篇博文的主题：一个计算具有金刚石/闪锌矿结构的晶体能带结构的工程项目，适用于各种元素。该项目已经在GitHub上：Empirical Pseudopotential¹，在我写这篇博文之前它已经存在了一段时间。实际上，GitHub上已经有另一个项目在这里也缺少描述，但关于那个项目，稍后会讲。

一如既往，这是程序运行的YouTube视频

一些理论

由于这是一个非常容易的主题，我将主要提供一些维基百科链接给不熟悉这些概念的访问者，然后是一些关于该主题的PDF文件链接，最后我会非常简要地写几句关于理论的话。

维基百科链接

您应该熟悉布拉维晶格。由于我们处理的是周期性的无限结构，我们可以通过应用傅里叶理论来研究它，所以您也应该理解倒格矢是如何使用的。它不仅仅是一个数学技巧，倒格矢与晶格衍射产生的散射矢量有直接的对应关系。布洛赫波在研究晶体时起着非常重要的作用，所以您也应该熟悉它们。一旦您理解了一个特定的布洛赫波可以有不同的分解，您就应该理解为什么我们可以限制在只研究第一个布里渊区。

一些论文

该程序基于 Marvin L. Cohen 和 T. K. Bergstresser 的论文Band Structures and Pseudopotential Form Factors for Fourteen Semiconductors of the Diamond and Zinc-blende Structures 中给出的参数²，所以我认为我应该首先引用它。一篇很好地描述了该理论的PDF文件在nanohub上：Dragica Vasileska 的Empirical Pseudopotential Method³。一篇包含理论和C代码的硕士论文是：Adam Lee Strickland 的Empirical pseudopotential method for modern electron simulation⁴。如果您理解了理论，您可能会想将该方法应用于其他材料，使用其他结构，例如石墨烯。在这种情况下，这里有一篇关于此的硕士论文：Srinivasa Varadan Ramanujam 的Band Structure of Graphene Using Empirical Pseudopotentials⁵。如果这些链接不够，谷歌应该能为您揭示更多与该主题相关的内容。

非常简要的理论

我们将从单电子薛定谔方程开始

$\left(\frac{\hat{p}^2}{2m} + \hat{V}\right)|\Psi\rangle = E |\Psi\rangle$

通过它，我们描述了非相互作用的电子在晶体原子核势场中运动。当然，这是一个多体问题，通常您无法简化，但有时即使是自由电子模型也可能足够好，或者至少是近自由电子模型，所以值得一试。

我们可以从原子的周期性排列中获益。因为哈密顿量与所有离散平移算子都对易，我们可以选择共同的特征向量来描述解，即布洛赫波。布洛赫波函数具有以下形式

$\Psi_{n,\vec{k}}\left(\vec{r}\right) = u_n(\vec{k}) e^{i \vec{k}\vec{r}}$

我们可以避免在平面波和周期函数u的分解上的歧义，通过限制在第一个布里渊区。由于靠近原子核的电子处于原子核的势阱深处，它们受其他原子核势场的影响很小，并且不容易隧穿出去（也就是说，它们不像我们试图描述的电子那样‘自由’），因此我们可以使用原子轨道线性组合来描述它们。一个问题是，布洛赫波函数与原子轨道不是正交的。这会迫使我们使用重叠矩阵并处理广义特征值问题。或者，我们可以通过通常的程序来正交化“平面波”：从向量中提取向量沿所有其他向量的投影

$|\psi\rangle = |\Psi\rangle - \sum\limits_{i} \langle \phi_i | \Psi \rangle | \phi_i \rangle$

其中 $| \phi_i \rangle$ 是原子轨道。很容易从这里得到原始波向量并将其代入薛定谔方程

$\left(\frac{\hat{p}^2}{2m} + \hat{V}\right)\left( |\psi\rangle + \sum\limits_{i} \langle \phi_i | \Psi \rangle | \phi_i \rangle \right) = E \left( |\psi\rangle + \sum\limits_{j} \langle \phi_j | \Psi \rangle | \phi_j \rangle \right)$

将哈密顿算子作用在每一项上，并将第二项移到左边，我们得到

$\left(\frac{\hat{p}^2}{2m} + \hat{V}\right) |\psi\rangle + \sum\limits_{i} \langle \phi_i | \Psi \rangle (E_i - E) | \phi_i \rangle = E |\psi\rangle$

我们可以将其强制回到薛定谔方程

$\left(\frac{\hat{p}^2}{2m} + \hat{V} + \sum\limits_{i} \langle \phi_i | \Psi \rangle (E_i - E) \frac{| \phi_i \rangle}{|\psi\rangle} \right)|\psi\rangle = E |\psi\rangle$

将排斥项称为 $\sum\limits_{i} \langle \phi_i | \Psi \rangle (E_i - E) \frac{| \phi_i \rangle}{|\psi\rangle} = \hat{V}_R$ ，并注意到我们仍然有一个薛定谔方程，具有相同的特征值但不同的特征向量，它不仅具有原子核给出的吸引势，还具有核心电子给出的排斥势。将 $\hat{V}_p = \hat{V} + \hat{V}_R$ 称为赝势，我们必须解

$\left(\frac{\hat{p}^2}{2m} + \hat{V}_p\right)|\psi\rangle = E |\psi\rangle$

除了求解它之外，唯一 remaining 的问题是确定赝势。您可以通过各种方法来做到这一点，从计算它到从实验结果中确定它。我让您从提供的链接中查找详细信息……

代码

如果您只对计算感兴趣，所有相关的类 – 除了 `Vector3D`（它已经在该博客的许多项目中用于此） – 都位于 `EmpiricalPseudopotential` 命名空间中。`SymmetryPoint` 是一个非常简单的对象类：它只有一个‘名称’和在倒空间中的位置。`SymmetryPoints` 包含一个此类对象的映射。它的构造函数只是用面心立方晶格的临界点对其进行初始化。`SymmetryPoints::GeneratePoints` 生成要绘制的点。您将一个由对称点组成的‘路径’传递给它，以及您希望在图表中拥有的点数，它就会计算 k 点坐标。它还会返回 `symmetryPointsPositions` 中对称点的索引。`Material` 类也很简单，它用于存储材料名称、基元胞原子之间的波尔距离和赝势的对象。有一个 `Materials` 类，它只包含一个材料映射。该映射在构造函数中填充，配置文件可能是一个更好的选择，但对于本项目而言，它已经足够了。顺便说一句，您会看到一些转换，这是因为参数是以埃（Angstroms）和里德堡（Rydbergs）给出的，而计算是以哈特里原子单位（Hartree atomic units）进行的。

`Pseudopotential` 类只存储赝势的参数，并有一个函数

std::complex<double> Pseudopotential::GetValue(const Vector3D<int>& G, const Vector3D<double>& tau) const
{
    const int G2 = G * G;
    const double Gtau = 2. * M_PI * tau * G;

    double VS = 0;
    double VA = 0;

    if (3 == G2)
    {
        VS = m_V3S;
        VA = m_V3A;
    }
    else if (4 == G2)
    {
        VA = m_V4A;
    }
    else if (8 == G2)
    {
        VS = m_V8S;
    }
    else if (11 == G2)
    {
        VS = m_V11S;
        VA = m_V11A;
    }

    return std::complex<double>(cos(Gtau) * VS, sin(Gtau) * VA);
}

它相当简单，因为我们只需要在某些点的值，超过 11 个对于 $G^2$ ，该值可以很好地近似为 0，对于 $G^2=0$ ，该值可以设为零，因为它只是能量的偏移。

以上提到的所有类都应该很容易理解。我只剩下两个类要描述，它们都很重要。首先，一个用于哈密顿量的类

class Hamiltonian
{
public:
    Hamiltonian(const Material& material, const std::vector<Vector3D<int>>& basisVectors);

    void SetMatrix(const Vector3D<double>& k);
    void Diagonalize();

    const Eigen::VectorXd& eigenvalues() const { return solver.eigenvalues(); }
protected:
    const Material& m_material;
    const std::vector<Vector3D<int>>& m_basisVectors;

    Eigen::MatrixXcd matrix;

    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> solver;
};

您可能在这里认出了构造哈密顿量的材料、用于哈密顿量矩阵表示的基矢、实际的哈密顿量矩阵以及用于对角化它的求解器。您可以查看项目源码¹以获取完整的实现，这是 `SetMatrix` 的实现

void Hamiltonian::SetMatrix(const Vector3D<double>& k)
{
    const unsigned int basisSize = static_cast<unsigned int>(m_basisVectors.size());

    for (unsigned int i = 0; i < basisSize; ++i)
        for (unsigned int j = 0; j < i; ++j) // only the lower triangular of matrix is set because the diagonalization method only needs that
            // off diagonal elements
            matrix(i, j) = m_material.pseudopotential.GetValue(m_basisVectors[i] - m_basisVectors[j]);

    for (unsigned int i = 0; i < basisSize; ++i)
    {
        // diagonal elements
        // this is actually with 2 * M_PI, but I optimized it with the /2. from the kinetic energy term
        const Vector3D<double> KG = M_PI / m_material.m_a  * (k + m_basisVectors[i]);

        matrix(i, i) = std::complex<double>(2. * KG * KG); // 2* comes from the above optimization, instead of a /2
    }
}

如果您仔细看过理论部分，您会立即认出这些公式。

最后一个是 `BandStructure`，稍微复杂一些，但也不是什么大问题。`AdjustValues` 只是将哈特里转换为 eV 并将能量值移到零点，这个零点由 `FindBandgap` 找到。`Initialize` 是直接的，`GenerateBasisVectors` 只是用具有适当长度的向量填充 `basisVectors`（其长度的平方在 `G2` 中）。

这是 `Compute` 函数

std::vector<std::vector<double>> BandStructure::Compute(const Material& material, unsigned int startPoint, unsigned int endPoint, unsigned int nrLevels, std::atomic_bool& terminate)
{   
    std::vector<std::vector<double>> res;

    Hamiltonian hamiltonian(material, basisVectors);

    for (unsigned int i = startPoint; i < endPoint && !terminate; ++i)
    {
        hamiltonian.SetMatrix(kpoints[i]);
        hamiltonian.Diagonalize();

        const Eigen::VectorXd& eigenvals = hamiltonian.eigenvalues();

        res.emplace_back();
        res.back().reserve(nrLevels);

        for (unsigned int level = 0; level < nrLevels && level < eigenvals.rows(); ++level)
            res.back().push_back(eigenvals(level));
    }

    return std::move(res);
}

`kpoints` 是一个倒空间位置的向量，之前在 `Initialize` 中使用 `SymmetryPoints::GeneratePoints` 生成。您可以在上面看到 `Hamiltonian` 类是如何使用的。代码只是在 `startPoint` 和 `endPoint` 之间迭代，对于每个 k 点，设置哈密顿量矩阵，然后对其进行对角化，然后检索特征值并保存在结果中。`startPoint`、`endPoint` 和 `terminate` 用于多线程计算。

有关这是如何完成以及数据如何放入图表的信息，您需要查看项目源码¹，主要是在 `EPThread` 和 `EPseudopotentialFrame` 类中。

`Options` 用于...选项，而 `OptionsFrame` 是实现用于编辑它们的选项对话框的类。`EPseudopotentialApp` 是一个非常简单的wxWidgets⁶应用程序类，`wxVTKRenderWindowInteractor` 是一个允许从 wxWidgets 轻松使用VTK⁷的类，我从这里⁸获得了这个类，并对其进行了轻微修改，以便能够与更新版本的 wxWidgets 编译，大致就是这样。

使用的库

现在应该很明显了，我从 mfc 转向了wxWidgets⁶。我厌倦了 MFC，尽管它与 VisualStudio 一起使用更容易，但自从他们发布了一个在 Windows 上无问题编译的 wxWidgets 版本以来，我决定从现在开始将其用于本博客的项目。最大的好处是，只需付出最小的努力，就可以将项目编译为不仅在 Windows 上运行，还可以在 Linux 或 Mac 上运行。自从我切换到 wxWidgets 以来，为了保持代码的可移植性，我还不得不避免使用我在其他项目中使用的 2D 图表类。与其编写一些新的、可移植的图表代码，我宁愿使用VTK⁷来制作 2D 图表。

与许多其他项目一样，我也使用了Eigen⁹。

结论

到此，《经验赝势》项目的介绍就结束了。如果您发现任何问题，请告诉我，无论是这里还是在 GitHub 上。

由于我非常忙，我可能不会在这里再发帖了，但 GitHub 上已经有另一个项目，也用于计算能带结构，但使用紧束缚方法。该项目基于这个项目，我只需删除一些代码并在某些地方进行更改（显然，例如哈密顿量需要更改）。代码甚至比这个项目中的代码更简单。事实上，它 deceptively simple，理论如果深入细节的话并不那么简单。

Empirical Pseudopotential GitHub 上的项目 ↩ ↩ ↩
Band Structures and Pseudopotential Form Factors for Fourteen Semiconductors of the Diamond and Zinc-blende Structures，作者：Marvin L. Cohen 和 T. K. Bergstresser，Phys. Rev. 141, 789 – 1966年1月14日发布 ↩
Empirical Pseudopotential Method，作者：Dragica Vasileska ↩
Empirical pseudopotential method for modern electron simulation，作者：Adam Lee Strickland ↩
Band Structure of Graphene Using Empirical Pseudopotentials，作者：Srinivasa Varadan Ramanujam ↩
wxWidgets 跨平台 GUI 库 ↩ ↩
VTK, The Visualization Toolkit 科学（及非科学）可视化库 ↩ ↩
wxVTK ↩
Eigen 矩阵库。 ↩

这篇帖子 Empirical Pseudopotential 最初出现在 Computational Physics。