量子力学软件开发工具包：TBTK的设计与应用

|⟩|⟩SoftwareX 9（2019）205原始软件出版物量子力学软件开发工具包克里斯托弗·比约恩松丹麦哥本哈根大学尼尔斯·玻尔研究所，Juliane Maries Vej 30，DKar t i cl e i nf o文章历史记录：收到2018年8月7日2019年2月21日收到修订版，2019年保留字：量子力学SDKC++数据结构a b st ra ctTBTK是用于量子力学计算的软件开发工具包，旨在开发研究二次量子化形式问题的应用程序。它还使方法开发人员能够创建紧束缚，DFT，DMFT，量子传输等的求解器它们可以很容易地相互结合。通过开发全新的求解器，以及前端和后端到已经很好地建立的包。TBTK提供了为二次量子化量身定制的数据结构，这将鼓励可重用性并实现量子力学计算的可扩展性。©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0.3用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2018_138合法代码许可证Apache许可证2.0代码版本系统使用git软件代码语言使用C++、python编译要求，操作环境依赖性BLAS，LAPACK，CMake，（可选：ARPACK，FFTW3，OpenCV，cURL，SuperLU（v5.2.1），wxWidgets，CUDA，HDF5，OpenBLAS，OpenMP，Google Test）Linux、OS X、类Unix如果可用，链接到开发人员文档/手册www.second-quantization.com支持问题电子邮件kristofer.gmail.com1. 介绍半个多世纪以来，半导体技术的进步推动了技术进步，摩尔定律描述了指数级进步。这背后的主要驱动因素是晶体管尺寸的不断减小。国际半导体技术路线图（ITRS）的目标是在2021年实现5纳米技术（http://www.itrs2.net/），但进一步降低是困难的。在晶格常数为0.54 nm的情况下，边长为5 nm的硅立方体大约为9个晶胞宽，并且包含不超过几千个硅原子。在这个尺度上，量子力学效应开始占主导地位[1，2]。因此，需要容易获得更准确的模型，以补充过去足以用于工业目的的半经典模型[3]。电子邮件地址：kristofer. second-tech.com。https://doi.org/10.1016/j.softx.2019.02.005同时，计算能力的提高，结合算法的开发，增加了可以使用量子力学模型模拟的原子数量[4因此，学术界和工业界都可以使用的系统规模已经重叠。这一点在量子计算领域的最新进展中尤为明显，学术界和工业界正在共同取得重大进展[9为使这一努力取得成功，重要的是要开发工具和程序，使来自不同领域的专家能够有效地相互合作。特别地，需要提供量子力学量的高级抽象的数据考虑符号<$lmσ（x， y， z），它在数学上等价于使用符号在哪里h是线性希尔伯特空间指数。前代表是一种高级抽象，特别适合于特定于模型2352-7110/©2019作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx206K. Björnson/SoftwareX 9（2019）205∑H=H+H=acc+H， （1）0IijjI我h'我}IJ推理，而后者是一个低层次的表示，适合于有兴趣实现计算要求的通用算法的方法开发人员。如果在两种表示之间没有一个通用的映射，低级的设计决策必然会在代码中向上传播。它们可以一直传播到最终用户，也可以通过高级的、特定于应用程序的接口隐藏起来在前一种情况下，最终用户有责任理解低级约定，而在后一种情况下，代码的通用性可能受到限制。在这两种情况下，代码都很难与其他软件集成，因为缺乏一套普遍认可的约定。在本文中，我们提出了TBTK，一个SDK的建模和求解哈密顿的二次量化形式。在其核心，它解决了上述映射问题它还提供了一组广泛的通用数据结构，可用于实现新的应用程序，求解器以及现有软件包的前端2. 数据结构：抽象和效率TBTK的主要部分是一个C++库，其中包含的数据结构旨在简化研究特定量子力学问题的应用程序的开发，以及使开发人员能够实现通用的可重用求解器。数据结构旨在提供抽象，允许开发人员专注于物理而不是数字，并提供与高度优化的单一用途代码相同的效率。通过强调面向对象的设计，代码被划分为受强封装保护的逻辑单元，这使得开发人员能够在适合给定任务的抽象级别上工作。3. 第二量化TBTK应用的起点是第二量子化形式†我IJ其中aij是复数，i和j是离散指数，并且ct和ci分别是状态i的创建和湮灭算子。初步的支持可用于相互作用项HI，但在这个简短的介绍中，我们专注于非相互作用哈密顿量。在TBTK表示法中，复数aij被称为跳跃幅度，这是从代码最初专注于紧束缚计算中推导出来的然而，命名更普遍的动机是，当薛定谔方程ihêt|（t）|中文（简体）是用有限差分手头的物理问题，以及希尔伯特空间指数，它们是相应物理指数的线性表示。从物理索引到希尔伯特空间索引的典型映射，例如对于具有索引结构（x，y，s）的自旋的二维晶格，可以被硬编码到应用中，如h = 2*SIZE_Y*x +2*y + s。这种显式映射的问题在于，它强制应用程序的每个方面都使用这种约定，既限制了代码的适用性，又在代码的每个级别上留下了不必要的数字细节。TBTK通过灵活的索引和用于跳变幅度和索引的复杂存储结构的组合来解决这个问题，该结构自动提供物理索引和希尔伯特空间索引之间的有效映射因此，应用程序开发人员可以专门使用物理索引，而方法开发人员可以编写仅依赖于希尔伯特空间索引的通用求解器。在TBTK中，使用花括号（例如{x，y，s}）指定为了完全通用，TBTK还允许索引具有不同的同时使用索引结构，例如，系统由分别具有索引结构（x，s）和（x，y，s）的两个子系统组成。要做到这一点，唯一的要求是索引结构首先不同之处左侧的子索引中的索引不同。这很容易通过在前面添加子系统索引来解决，从而产生数字索引{0，x，s}和{1，x，y，s}。5. 创建模型跳频幅度由其复值和两个指数i和j唯一确定。从状态{x，y，s}到{x+1，y，s}的值为1的跳变幅度由下式创建：HoppingAmplitude（1，{x+1，y，s}，{x，y，s}）;此外，所有跳跃幅度都存储在称为模型的对象内，例如，可以如下设置具有最近邻居跳跃的正方形网格。模型型号;for（intX = 0的整数;X < SIZE_X; x++）{for（int y = 0的整数;y < Size_X;（y ++）{for（intS=0的整数;S<二、（s ++）{public int findDuplicate（x+1）{模型跳跃幅度（1，{ x+1，y，s }，{ x，y，s }） + HC;}public intfindDuplicate（）{模型跳跃幅度（1，{x，y+1，s}，{x，y，s}） + HC;}}|（t）|Ψ (t)⟩,(3)}可以看出，aiji. 也就是说，粒子从状态j跳到状态i。4. 物理指标与希尔伯特空间指标为了允许指定任意复杂度的Hamilton，TBTK提供了一个灵活的索引方案。一个重要的区别是在物理指标之间，如（x，y，z，亚晶格，轨道，自旋），这些指标与物理指标有直观的联系这里的通过允许用户使用具有任意结构的物理指标来指定模型，在已知aij的情况下，相对容易地指定实际上任何感兴趣的哈密顿量事实上，甚至可以通过将所谓的回调函数作为第一个参数传递给跳跃幅度来指定时间相关的或依赖于一些尚未确定的参数的哈密顿量如果例如某些参数要自洽地确定或者如果它们需要计算，K. Björnson/SoftwareX 9（2019）205207Fig. 1. 随着跳跃幅度被添加到模型中，树结构正在构建。跳变幅度沿着树的分支向下传递并存储在叶节点上，其中第二跳变幅度索引的第n个子索引确定在第n层向下走哪个分支。如果将与图右侧的几何结构一致的跳变幅度添加到模型中，则得到的树结构如图左侧所示。 调用model.construct（）时，将创建线性希尔伯特空间基通过遍历树并以从0开始递增的数字来枚举叶节点。一些未知的重叠积分有关此方面的更多信息，请参阅文档。1一旦将所有相关的跳跃幅度添加到模型中，就使用下式创建物理索引与希尔伯特空间索引之间的映射：模型return（）;仔细看看这个调用的幕后发生了什么，有助于理解为什么我们可以在不产生显著性能损失的情况下提供物理索引的额外便利。跳变幅度存储在一个树结构中，见图11。1，它可以被认为是稀疏矩阵格式，每列存储在叶节点上。与每一列相关联的物理索引被编码在树结构本身中，而线性希尔伯特空间索引则存储在叶节点上，一旦进行了model.construct（）调用。除了充当跳跃幅度的存储器之外，树结构因此还提供物理索引与希尔伯特空间索引之间的映射请注意，由于映射只包括实际包含在模型中的站点，因此这会导致最小希尔伯特空间。由于物理索引可以通过遍历树到相应的叶子节点来转换为希尔伯特空间索引，因此从物理索引到希尔伯特空间索引的转换的时间复杂度在希尔伯特空间大小上是O（1）。单个反向查找的时间复杂度稍微复杂一些，因为它取决于每层的分支数，但需要反向查找的算法可以通过迭代树结构在O（N）时间内构建查找表，其中N是希尔伯特空间大小。方法开发者特别感兴趣的是，可以在O（M）的时间内遍历所有跳频幅度并提取它们在希尔伯特空间基中的值，其中M是跳频幅度的数量。这实际上在希尔伯特空间大小中也总是线性的，因为对于具有局部算子的任何模型，跳频幅度的数量与希尔伯特空间大小成比例（情况就是这样 用于紧约束、有限差分、有限元等）。虽然我们没有详细讨论许多粒子系统的任何初步支持，但我们注意到，上述树结构描述了哈密顿量的（双线性）单粒子部分。交互术语存储在单独的结构中，但是树结构提供的映射允许1 安 装 doxygen 后 ， 可 以 使 用 make documentation 构 建 文 档 ， 然 后 可 以 在BuildFolder/doc/html/index.html中找到，其中BuildFolder是构建文件夹。最新版本的文档也可以在http://www.second-quantization.com上找到。对于相互作用项中的算子，也可互换地使用它们的物理和希尔伯特空间指数。当相关时，模型对象还可以保存诸如温度、化学势和统计数据模型setTemperature（300）; model .设置化学势（0）;模型s et S t a t i s t i c s（S t a t i s t i c s：：FermiDirac）;因此，模型对象是用于模型相关信息的通用容器，并且还可以包含此处未示出的其他信息。6. 解决者不同类型的问题需要不同的求解方法，这些方法通常在TBTK的求解器中实现。 由于求解器通常是花费大部分计算时间的地方，因此方法开发人员必须完全自由地选择内部用于实现其算法的数据结构和编程范式。此外，对于通用求解器，它们可以在对模型的最小假设下工作也很重要。从物理指标到希尔伯特空间指标的映射提供了关键来解决这两个问题。在内部，求解器可以使用希尔伯特空间基从模型中请求跳频幅度，并设置最适合该方法的数据结构。通过这种方式，方法开发人员可以创建新的求解器，而无需担心特定物理模型的细节，而应用程序开发人员可以指定模型，而无需担心特定求解器的方法特定细节。TBTK包含许多不同的求解器，例如可以执行对角化，Arnoldi迭代和Green函数的Chebyshev展开使用对角化作为示例，可以设置并执行典型的求解器如下解算器：对角化解算器;解算器。setModel（model）;solver.return（）;对于其他求解器，初始化可能需要提供更多方法特定的参数。然而，主要的想法是提供一个接口的应用程序开发人员，最小化的必要的方法特定的知识，同时提供的可能性，配置求解器的需求。由于哈密顿量可以在O（N）时间内转换为最适合算法的任何格式，并且几乎所有求解方法的执行时间都是N中的超线性，因此这种转换的代价可以忽略不计。方法开发人员在求解器内设置计算时可以应用的基本配方如下所示。/ / I t er at e超过的跳跃幅度constHoppingAmplitudeSetHoppingAmplitudeSet=模型。return（）;因为（HoppingAmplitudeSet：： C onst Iterator i t e r a t or r=跳幅设置。 cbegin（）; i t e r a t or！=跳幅设置。return（）;++ i t e r a t o r）的情况下，/ /提取的振幅和 物理/ / i indices从的跳跃幅度。复数double> 振幅208K. Björnson/SoftwareX 9（2019）205××∗∗图二. （左）设置N所花费的时间N作为Hilbert空间基函数的具有现场和最近邻跳跃项的N立方紧束缚模型尺寸TBTK大约快4倍（在基本尺寸为10 6时快3.93倍）。（右）提取相同的哈密顿量所花费的时间，作为希尔伯特的函数空间基尺寸。在TBTK中将矩阵提取为SparseMatrix对象的速度比在kwant中将其提取为稀疏numpy数组的速度快约3倍（2.86倍更快的基础大小106）。方法开发人员可以通过直接访问原始稀疏矩阵数据（快6.15倍），在TBTK中挤出额外的2倍在基本尺寸106处）。=（i t e r a t o r）.return（）;const指数 toIndex=（i t e r a t o r）.return（）;const指数 fromIndex=（i t e r a t o r）.return（）;/ /转换的物理因迪塞斯到我在耳朵里/ / i indices.行间= hoppingAmplitudeSet。获取基本索引（toIndex）的情况下;int列= hoppingAmplitudeSet。获取基本索引（fromIndex）的情况下;/ /添加的矩阵元件到/ / Hamiltonian对的格式最好 适合/ / for的给定算法/ / . . .}7. 萃取性能尝试完全不同的解决方案的方法，只需改变几行代码有关的求解器初始化。然而，我们确实注意到，不是每个求解器都可以计算每个属性，并且一些求解器可以计算一些特定的细节，其他求解器根本无法提供。因此，属性提取器只是近似一致的，有时提供简单地打印相应的求解器不能用于计算给定属性的函数的实现，而有时具有在其他属性提取器中不可用的附加函数。设置属性提取器的典型表达式如下所示。属性提取器：：对角化器属性提取器（求解器）;属性Extractor。setEnergyWindow（LOWER_BOUND，上限，分辨率）的情况下;然后，我们可以提取每个站点上的密度，对旋转求和：Property：：Density=属性提取器 . 计算密度（{{ IDX_ALL，IDX_ALL，IDX_SUM_ALL}}TBTK 定 义 了 一 些 性 质 ， 如 本 征 值 、 波 函 数 、 态 密 度（DOS）、（自旋极化）局域态密度（LDOS）等。然而，不同的求解器可以在内部使用非常不同的存储结构，并且希望将求解器的责任限制为处理使用Hilbert空间索引制定由于这个原因，TBTK提供了属性提取器，这些属性提取器在求解器的方法特定细节和呈现给应用程序开发人员的更高抽象层之间架起了桥梁。强烈建议方法开发人员与他们的求解器并行创建类似的属性提取器。属性提取器为应用程序开发人员提供了更直观的界面，允许应用程序开发人员使用物理索引从求解器中提取属性。此外，他们的目标是提供统一的接口，为求解器到外部世界。因此，使用属性提取器的代码通常可以工作，即使求解器发生更改，）的情况下;计算态密度（DOS）：Property：：DOS dos=属性Extractor。calculateDOS（）;获取特征值：属性：：EigenValues eigenValues=属性Extractor。return（）;计算所有自旋组合的延迟格林属性：：GreensFunction greensFunction=属性提取器 . 计算格林函数（{{{0、 0，IDX_ALL }， { 5，5，IDX_ALL }，Property：：GreensFunction：：Type：：Retarded）的情况下;K. Björnson/SoftwareX 9（2019）205209沿y = SIZE_Y/2提取LDOS，对旋转求和：属性：：LDOSldos=属性提取器 . 计算LDOS（{{IDX_ALL，SIZE_Y/2，IDX_SUM_ALL}}）的情况下;计算所有部位的磁化强度：产品性质：：磁化强度=属性提取器 . 计算磁化强度（{{ IDX_ALL， IDX_ALL， IDX_SPIN）的情况下;在现场获得自旋极化LDOS（2，4）和（3，5）：属性：：SpinPolarizedLDOS spinPolarizedLDOS=属性提取器.计算自旋极化LDOS（{{ 2， IDX_SPIN }，{ 3， IDX_SPIN }，}）的情况下;提取状态1、3和7的所有索引的波函数：Property：：WaveFunctionswaveFunctions=属性提取器 . 计算波函数（{{ IDX_ALL，IDX_ALL，IDX_ALL } }，{ 1， 第三、七节}）的情况下;8. 基准我们执行了一些基准测试来量化扩展的性能，并提供令人信服的证据，通用数据结构的引入带来了可以忽略不计的性能损失。特别是，考虑稀疏矩阵格式的哈密尔顿算子的模型规范和检索，因为这是将应用程序开发人员代码与方法开发人员代码分开的层。因此，这些基准提供了数据结构的开销成本的上限。基准测试是针对Kwant[15]进行的，因为它能够实现同样的事情，虽然使用了一个不太一般的形式主义，其中坐标和其他指数是在不同的基础上处理的。有关实际使用的代码，请参见补充材料中的PerformanceTest.zip为了简单起见，创建了具有最近邻跳跃的立方紧束缚模型，然后将其转换为上的稀疏矩阵。运行在单核Intel（R）Xeon（R）CPU@2.30 GHz上的Ubuntu 18.04的全新安装。在图2的左侧面板中，示出了作为希尔伯特空间大小的函数的建立模型的时间。TBTK和kwant都显示线性缩放，TBTK的性能几乎是kwant的4倍。在图的右侧面板中。图2绘制了在稀疏矩阵格式上提取哈密顿量的时间。此外，TBTK和kwant都是线性缩放的，根据所选择的特定稀疏格式，TBTK比kwant高出3到6倍。在这两种情况下，建立模型和请求数据的开销对于TBTK或kwant来说都不是特别大。事实上，对于基大小为106的问题，TBTK和kwant在稀疏矩阵格式上指定和提取模型的累积时间分别为10 s和38 s，几乎可以肯定，在任何算法中工作的系统的大小。9. 影响和总结TBTK是一个软件开发工具包，可以快速开发计算量子力学特性的应用程序。它将帮助科学界和工业界开发代码，通过鼓励可重用性的可扩展方法实现大规模合作一个更加一体化的社区对于量子技术扩展到工业水平至关重要。这也意味着错误等可以更快地被发现，并提供了一种实现科学工作应有的透明度的手段。这将使社区花费更少的时间来复制其他人已经解决的数字细节，并将更多的精力放在感兴趣的物理问题为了实现这一目标，TBTK提供了一组高效的通用数据结构和构建工具，这些工具借鉴了软件开发中的最新最佳实践其目的是开发一个求解器和工具的生态系统把不同的方法结合起来。特别是，TBTK旨在通过提供允许开发人员在更高的抽象级别上工作的数据结构来帮助这种开发数字细节。感兴趣的读者可以参考文档2了解更多信息。致谢感 谢 Anna Sinelkova 、 Igor Di Marco 、 Mahdi Mashkoori 、Oladunjoye Awoga、Andreas Theiler和Tom McClintock对代码和/或手稿进行了有益的讨论和反馈。感谢Annica M。作者声明：BrianM. Andersen的富有成效的合作，帮助塑造了代码。这项研究得到了瑞典研究委员会（Vetenskapsrådet，批准号621-2014-3721）、瑞典Göran Gustafsson基金会、瑞典战略研究基金会（SSF）、瑞典Knut和Alice Wallenberg基金会（通过Wallenberg学院研究员计划）以及丹麦独立研究基金会（批准号DFF-6108-00096）的财政支持。附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2019.02.005上找到。引用[1]达塔湾介观系统中的电子输运。北京：清华大学出版社.[2] 达塔湾量子传输：原子到晶体管。北京：清华大学出版社.[3]施SM，吴KK。半导体器件物理学。Inc：John Wiley& Sons;2007.[4]ALPS项目2.0版：强相关系统的开源软件。J Stat Mech 2011;P05001.[5]Blaha P ， et al. WIEN2k ， a augmented plane wave + local orbitalsprogramfor calculating crystal properties.奥地利：维也纳工业大学;2001年。[6]Frisch MJ，et al. GAUSSIAN 09，revision A.02.沃林福德：高斯公司; 2009年[7]杨晓萍，李晓萍. Vienne Ab-Initio Simulation Package（VASP）。2018年12月28日，中国科学院出版社。https://www.vasp.at/[8]作者：J.J. M.Phys Rev Lett 2010;105：167006。[9]Gibney E 量子计算机获得设计升级。 Nature2017;541：447.[10]卡斯泰尔韦基湾硅在量子计算竞赛中占据一席之地。Nature2018;553：136.2见脚注1。210K. Björnson/SoftwareX 9（2019）205[11]卡斯泰尔韦基湾 量子云走向商业化。 Nature2018;543：159.[12]Zeng W，et al.第一量子计算机需要智能软件。Nature2017;549：149.[13]Gibney E 数十亿欧元的量子项目正在形成。 Nature 2017;545：16.[14]波普金湾量子物理学可能会得到美国国会的大力推动。科学2018;360：1158.[15]Groth CW，Wimmer M，Akhmerov AR，Waintal X. Kwant：一个量子传输软件包。New J Phys 2014;16：063065.