求解泊松方程

aroman

4.50/5 (2投票s)

2021 年 7 月 21 日

GPL3

14分钟阅读

7803

这篇文章也介绍了 VTK，即可视化工具包

引言

当我开始写这个博客时，我就已经预料到会有使用快速傅里叶变换的项目。实际上，我为博客写下了几个主题构想，包括求解泊松方程和本文将引出的主题。我在弛豫法一文中已经提到，可以使用离散傅里叶变换来更快地解决问题，现在它来了，作为通往至少一个关于密度泛函理论项目的中间步骤。

我不会在这里讨论快速傅里叶变换，我已经在之前的文章中讲过了，所以这里不会介绍处理 FFT 的类。我只是从上一个项目中复制了它们并直接使用。我也会在 DFT 项目中重用它们，连同泊松求解器（可能会有改动）。

本文还介绍了VTK，即可视化工具包¹，这个库我也会在未来的项目中使用。即使是对于这个项目，虽然我可以用 OpenGL 轻松实现那些曲面——但要让坐标轴像 VTK 那样表现就不那么容易了，实现起来可能要花很长时间而且很枯燥——这个库也为我节省了不少时间。VTK 能够对科学数据进行非常 впечатляющие 的可视化，从而节省大量时间（据说这个库凝聚了超过1000年的开发者时间，显然我不可能在合理的时间内超越它）。使用它能让我更专注于博客主题相关的部分，并能让我实现一些因为可视化复杂性而本会避免的项目，就像我在格子玻尔兹曼一文中避免了3D流体那样。我知道，如果能直接下载代码并轻松编译，依赖越少越好，但 VTK 依赖带来的优势太多了，不容忽视。未来某些项目我可能只用 OpenGL，但对于复杂的3D可视化甚至是图表绘制，除了像我已在本博客项目中用过的非常简单的图表外，我都会使用VTK。

在本文中，我用谱方法实现了一个泊松方程求解器。下面是程序运行的截图：

你在这里看到的是具有周期性边界条件的三维体的一半截面。在左视图中，我表示了由两个高斯函数生成的电荷密度；在右视图中是泊松方程的解。该项目的源代码在GitHub上²。

链接

我之前不想透露在博客上引入 FFT 主题的目的，因为我担心可能会因为各种原因放弃实现 DFT 项目。但现在我觉得，在我实现了一个 DFT 原型之后，分享一些我查阅过的链接会很有趣。顺便说一下，文章顶部的特色图片就是我在实现 DFT 原型时用 Octave 生成的一张图表。

起初，我只想用实空间离散化写一个简单的 DFT 项目。我以前已经用这种方式实现过一个简单的 DFT 程序——概念验证性质的。离散化实空间并使用像有限差分法这样的简单方法是一个不错的起点，这种方法易于理解和实现。但如果你想要一些性能，这可能不是最好的方式。所以这次我考虑使用高斯轨道，类似于Hartree-Fock项目那样。然后我想，也许某天我会想实现一个晶体的项目，这时周期性边界条件就派上用场了。在这种情况下，由于周期性，平面波基组非常好，而且它也为这个博客带来了一些新东西。借助 Hartree-Fock 和 DFT 的文章，应该有足够的信息来实现一个使用高斯轨道的 DFT 程序了。

于是，我决定采用平面波基组的方式，然后开始在网上寻找相关信息。理论我懂，实现方面也了解一些，但我想更深入地研究细节，毕竟细节决定成败。我开始做 Hartree-Fock 项目时，并没有把所有细节都弄清楚，结果花的时间比预期的要长。当然，关于 DFT 的讲座有很多，我会在关于 DFT 的文章中指出它们，但它们只涉及理论。在 YouTube 上稍作搜索，我找到了更合适的内容，来自康奈尔大学：

如果我没记错的话，前两节课和第三节课的一部分与本项目相关。

不幸的是，它缺少了一些部分，但借助一些链接可以填补这些空白。这里有一个非常有用的链接（我在讲座中看到这个链接，才找到了它）：Arias 课程³。你应该查看 Phys 7654 Practical DFT 文件夹，最好是最新的。在那里你会找到阅读材料，以及非常重要的作业。其中一篇理论论文也可以在 arxiv 上找到：密度泛函计算的新代数表示法⁴。我刚刚用 Octave 实现了这些作业，以确保它们相对容易完成，现在我只需要把它们“翻译”成 C++，当然，如果可能的话，还要有一个更漂亮的可视化。这篇博客文章的项目对应第一个作业，但为了使项目更易于理解，计算过程被简化了。C++ 项目中没有“重叠”算符（因为基是正交的），我也去掉了各种 cI、cIdag、cJ、cJdag，对于这个项目来说，正向和反向 FFT 就足够了。

一些理论

泊松方程出现在各种情境中，但目前我们关心的是电磁场。方程如下：

$\nabla^2\varphi= -\frac{\rho}{\epsilon}$

其中 $\nabla^2=\Delta$ 是拉普拉斯算符。顺便说一下，从现在开始我将考虑库仑规范。你可以直接用数值方法处理这个方程，通过离散化空间，在二维情况下甚至可以用多重网格方法击败 FFT，详情请查看此页面。它的解是：

$\varphi(\vec{r}) = \frac{1}{4\pi\epsilon}\int\limits_V \frac{\rho(\vec{R})}{|\vec{R}-\vec{r}|} d^3R$

更多细节，请查看此页面：电磁场的数学描述。
毫不奇怪，由于线性关系，场由体积内每个“点”电荷产生的场的“和”（积分）给出。你可以尝试用数值方法求解这个积分，对于少数点电荷来说这并不困难。对于数量有限的具有球对称性的局部电荷分布，你也可以用一些技巧来简化问题，但总的来说，这个积分在实空间中很难求解。对于三维情况，如果将体积切成 $M=N^3$ 块，你必须对每个点都对所有块进行求和，也就是说，复杂度是 $O(N^6)$ 。而 FFT 的复杂度是 $O(M \log(M))$ ，速度快得多。所以我们对泊松方程使用傅里叶变换：

$\nabla_r^2 \int \varphi(\vec{k}) e^{i \vec{k} \vec{r} }d^3k = -\frac{1}{\epsilon} \int\rho(\vec{k}) e^{i \vec{k} \vec{r} }d^3k$

你可以将拉普拉斯算符与积分交换位置，而且由于该算符不作用于 k，它可以越过 $\varphi$ ，只作用于指数部分。它只是将 $i\vec{k}$ 因子带下来两次。所以在倒易空间中，拉普拉斯算符变得非常简单，即乘以 $-k^2$ 。这意味着你不需要处理耗时的对角化，算符本身就是对角的，你只需执行一个简单的乘法来应用它。要解泊松方程，你需要使用离散傅里叶变换计算倒易空间中的电荷密度 $\rho(\vec{k})$ ，然后通过将每个值除以 $\epsilon k^2$ 来求解，得到 $\varphi(\vec{k})$ ，最后再进行逆离散傅里叶变换回到实空间。

代码

FFT 代码位于 Fourier 命名空间中，我已在之前的文章中提到过。在这个项目中，我也使用了在其他项目中用过的 Vector3D<T> 类。我想它的功能应该很清楚，我也不再赘述。代码中还有其他一些东西（比如 Ewald 求和），我也不再详述，有关这些细节，请参考视频讲座和讲座文件。

泊松求解部分

所有与求解泊松方程相关的代码都在 Poisson 命名空间中。RealSpaceCell 和 ReciprocalSpaceCell 即使粗略地看一下代码也应该很明显，我认为甚至命名本身就足够说明问题了。如果还不清楚，这两个维基百科页面应该有相关信息：布拉维晶格、倒易晶格。此外，跟随我指出的讲座并查看作业也会非常有帮助。顺便说一下，如果你对这个主题感兴趣，你应该自己尝试用 Octave 完成这些作业，这里的代码可能与作业要求的并不完全一致。例如，你会在那里看到一个矩阵 S，我已经去掉了它（比如，晶胞体积是用 det(S) 计算的）。通常，你会把晶胞基向量放在该矩阵的列中，它的行列式就是晶胞体积，但由于在计算分子之前你只需要正交向量，所以我暂时去掉了它。不过，对于晶体来说，这通常是不成立的。虽然会增加复杂性，但从更普遍的角度看待事物还是有用的。我也去掉了讲座和作业中出现的重叠矩阵，因为平面波基是正交的，所以它是对角的，但如果你使用高斯基组，矩阵 O 就会非常重要。有一件事可能看起来很奇怪，那就是倒易晶胞中点的索引方式，这是为了缓解混叠⁵。如果你看看索引的样子以及晶胞是周期性的事实，就会发现那里发生的事情并不是什么大问题。

目前，代码只使用高斯电荷分布，GaussianChargeDistribution 类代表这样一个电荷。它内部有一个 position 向量和电荷 Z。一些这样的电荷被打包在 Charges 类的向量中，同时还有一些对计算有用的附加信息。现在，如果你看 ComputeStructureFactorAndChargeDensity 内部，你可能会注意到从傅里叶空间返回到实空间计算电荷密度的过程，这应该会提示一种优化方法，因为用这种方法求解泊松方程需要将电荷密度转换到傅里叶空间。我没有使用这种优化的原因是为了以一种通用的方式求解泊松方程，通常你是从实空间中指定的电荷密度开始的。但我添加了一段应该足够清晰的注释：你只需注释掉两行代码并取消另一行的注释，就可以切换到更快但通用性较差的解决方案。你可以利用电荷分布的生成方式：先生成一个中心在晶胞中心的高斯函数，然后将其转换到倒易空间，通过简单地“平移”（这就是 StructureFactor 存在的原因）该分布来生成所有电荷分布。

现在，这里是 PoissonSolver 类中关于求解泊松方程的相关代码部分，首先是将解带回实空间的方法：

static inline Eigen::VectorXcd SolveToRealSpace(Fourier::FFT& fftSolver, Poisson::RealSpaceCell& realSpaceCell, Poisson::ReciprocalSpaceCell& reciprocalCell, Charges &charges)
{
    Eigen::VectorXcd fieldReciprocal = SolveToReciprocalSpace(fftSolver, realSpaceCell, reciprocalCell, charges);

    Eigen::VectorXcd field(realSpaceCell.Samples());

    fftSolver.inv(fieldReciprocal.data(),field.data(),realSpaceCell.GetSamples().X, realSpaceCell.GetSamples().Y, realSpaceCell.GetSamples().Z);

    return field;
}

它只是调用了 SolveToReciprocalSpace，该方法返回倒易空间中的解，然后对其执行逆傅里叶变换，得到实空间中的解。下面是获取傅里叶空间中解的代码，连同前面提到的注释：

static inline Eigen::VectorXcd SolveToReciprocalSpace(Fourier::FFT& fftSolver, Poisson::RealSpaceCell& realSpaceCell, Poisson::ReciprocalSpaceCell& reciprocalCell, Charges &charges)       
{
    Eigen::VectorXcd fieldReciprocal(realSpaceCell.Samples());

    fftSolver.fwd(charges.ChargeDensity.data(), fieldReciprocal.data(), realSpaceCell.GetSamples().X, realSpaceCell.GetSamples().Y, realSpaceCell.GetSamples().Z);

    // uncomment this line and comment the two above if you want a faster solution          
    //Eigen::VectorXcd fieldReciprocal = realSpaceCell.Samples() * charges.rg;

    fieldReciprocal(0) = 0;
    for (int i = 1; i < realSpaceCell.Samples(); ++i)
    {
        // inverse Laplace operator
        fieldReciprocal(i) *= 4. * M_PI / realSpaceCell.Samples() / reciprocalCell.LatticeVectorsSquaredMagnitude(i);
    }

    return fieldReciprocal;
}

这里有两件事值得一提，首先是 4. * M_PI，这是因为使用了原子单位，更具体地说，是因为 $1/(4\pi\epsilon_0) = 1$ 。
另一件事是设置 fieldReciprocal(0) = 0;。如果你尝试用傅里叶空间中的拉普拉斯算符来计算它，会得到一个除以零的错误。使用这个方法的原因是零频率对应于实空间中的一个常数，我们方便地将其设置为零。我已经提到了库仑规范……

无论如何，关于泊松求解代码，我现在有耐心评论的就这么多了，其中还有很多内容，但如果你想了解更多，讲座和它们附带的文档应该足够了。

可视化部分

这个项目的第二个目标是介绍VTK，即可视化工具包¹，从现在起，我将在本博客的项目中使用它。他们有一本很好的教科书和一个很好的用户指南，请查阅它们。
该库可以在各种平台上使用，并支持并行处理。它不仅可以从 C++ 调用，还可以从其他语言，如 Python 或 Java 调用，尽管这对于本博客的项目来说不太重要。它使用管道式架构，一端是数据源（如文件、几何对象或像我在这个项目中使用的简单数据结构），另一端是在渲染窗口中进行渲染的渲染器。你可以使用多个渲染器在同一个窗口中拥有多个视图，我在这个项目中就利用了这一特性。在管道中，你可以有各种“过滤器”来处理你的数据并将其从一个过滤器传递到另一个，你已经有很多现成的算法（例如德劳内三角剖分）。这些对象可以有多个输入和输出。这个库能做很多事情，值得一看！

虽然该库具有强大的图表功能，例如 vtkChartXY 类，但我未来的博客项目可能会使用比图表更多的功能，尽管我也可能使用相对简单的图表。因此，在这个项目中，我生成并显示了3D曲面，而不是直接使用图表功能。处理 VTK 的代码位于 CPoissonDoc 中（用于数据源），其余部分在 CPoissonView 中。尽管 VTK 提供了智能指针，但在对象只创建一次（通常在构造函数中）并在最后于析构函数中删除的简单情况下，我没有使用它们。为了将数据暴露给 VTK，我将其放入“图像”对象中，这些对象表示均匀分布的数据。在这种情况下是2D，但也可以是3D（或1D）。例如，在文档对象中，你有 vtkImageData* fieldImage;，它在构造函数中像这样初始化（为清晰起见，类似的电荷密度代码已被移除）：

fieldImage = vtkImageData::New();

fieldImage->SetSpacing(realSpaceCell.GetSize().Y/realSpaceCell.GetSamples().Y, realSpaceCell.GetSize().Z/realSpaceCell.GetSamples().Z, 0);
fieldImage->SetDimensions(realSpaceCell.GetSamples().Y, realSpaceCell.GetSamples().Z, 1); //number of points in each direction

fieldImage->AllocateScalars(VTK_FLOAT, 1);

VTK 对象在文档析构函数中通过 fieldImage->Delete(); 删除。需要 ::New() 和 ->Delete() 是因为在不同平台上，实际创建的对象可能不同（例如，渲染窗口就是这种情况），New 是一个工厂方法。此外，这些对象是引用计数的，Delete 方法不仅仅是删除对象。各种对象可能会保留指向该对象的指针，并在传递给它们时（例如作为输入连接）增加其引用计数。

数据在求解泊松方程后在 Calculate 中设置：

// slice the result - put the values in the 'image' data for VTK

for (unsigned int i = 0; i < realSpaceCell.GetSamples().Y; ++i)
    for (unsigned int j = 0; j < realSpaceCell.GetSamples().Z; ++j)
    {
        unsigned int pos = start + realSpaceCell.GetSamples().Z * i + j;

        fieldImage->SetScalarComponentFromDouble(i, j, 0, 0, field(pos).real());
    }

它只是穿过“晶胞”的一个切片。

这就是文档类中的全部内容，大部分代码都在 CPoisonView 视图类中。我应该提到，他们实际上在“GUI支持”中有一个 vtkMFCWindow，但我没有使用它，我觉得我用我实现视图的方式有更多的控制权。

与文档类的情况一样，对象使用 ::New() 创建，一些在视图构造函数中，一些在 CPoissonView::OnInitialUpdate 中，并在视图析构函数中销毁。这些对象的一些属性也在 OnInitialUpdate 中设置。绘图在 CPoissonView::OnDraw 中完成，在窗口绘图的情况下，这非常简单，只需请求渲染窗口进行渲染（由于打印和打印预览，情况会更复杂）。顺便说一下，渲染窗口是由 VTK 嵌入在 MFC 窗口中的。你可以直接使用 MFC 窗口——也就是说，直接在其中绘图，而不是在子窗口中——但我无法用这种方式让交互器工作。

在同一个窗口中拥有两个视图很容易，下面是它在 Pipeline 实现中的工作方式：

void CPoissonView::Pipeline()
{
    PipelineView(ren1, geometryFilter1, warp1, mapper1, chartActor1, axes1);
    PipelineView(ren2, geometryFilter2, warp2, mapper2, chartActor2, axes2);
}

对 PipelineView 的第一次调用是针对左侧视图，第二次调用是针对右侧视图。下面是 PipelineView 的实现方式：

void CPoissonView::PipelineView(vtkRenderer *ren, vtkImageDataGeometryFilter* geometryFilter, vtkWarpScalar* warp, vtkDataSetMapper* mapper, vtkActor* chartActor, vtkCubeAxesActor2D* axes)
{
    warp->SetInputConnection(geometryFilter->GetOutputPort());    
  
    //mapper->SetInputConnection(warp->GetOutputPort());
    
    //***************************************************************************************
    // Gouraud shading needs normals. Just uncomment the above and comment what follows up to //****
    // to remove shading
    
    vtkSmartPointer<vtkPolyDataNormals> normals = vtkSmartPointer<vtkPolyDataNormals>::New();
    normals->SetInputConnection(warp->GetOutputPort());
    normals->SplittingOff();
    normals->ComputePointNormalsOn();
    normals->ComputeCellNormalsOff();
    normals->ConsistencyOn();
    mapper->SetInputConnection(normals->GetOutputPort());
    
    //****************************************************************************************

    chartActor->SetMapper(mapper);
    chartActor->GetProperty()->SetInterpolationToGouraud();
                    
    ren->AddActor(chartActor);
    
    // add & render CubeAxes
    
    axes->SetInputData(warp->GetOutput());
    axes->SetCamera(ren->GetActiveCamera());  
    ren->AddViewProp(axes);
}

你还应该查看 CPoissonView::OnInitialUpdate()。例如，数据源就在那里设置，但我会让你在 GitHub² 上查看其余的代码。

使用的库

除了 VTK¹，该项目还使用了 FFTW⁶ 和 Eigen⁷，当然还有 MFC。

结论

这是实现一个密度泛函理论程序的第二步。与前几步相比，下一步将会相当复杂，所以在看下一个项目之前，你应该对这个项目有一个很好的理解。我仍然不知道下一个什么时候能准备好，但最终我会完成它，目前它还只是一个 Octave 原型。

VTK，即可视化工具包卓越的科学（及其他领域）可视化库 ↩ ↩ ↩
泊松方程项目 GitHub 上的项目 ↩ ↩
Arias 课程包含康奈尔大学 Phys 7654 文件夹的谷歌云盘 ↩
密度泛函计算的新代数表示法 Sohrab Ismail-Beigi, T. A. Arias ↩
混叠与离散时间傅里叶变换作者 Steve Eddins ↩
FFTW 快速傅里叶变换库。 ↩
Eigen 矩阵库。 ↩

文章求解泊松方程首发于计算物理。