0 引言

目前高性能计算已成为当代科学与工程发展的关键性技术和手段.高性能计算已经能够用来模拟越来越精细和复杂的模型和问题, 并提供越来越可靠的计算结果.量子化学一直是高性能计算应用最活跃的领域之一, 量子化学的计算模拟已成为现代工业和学术研究的重要支柱之一.人们结合量子化学理论、计算机技术和编程算法等编写了大量的量子化学计算程序, 如Gaussian、NWChem、Q-Chem、MOLPRO、MPQC、ACES、GAMESS、PSI、Dalton、Turbomole、MOLCAS、PQS、ORCA和COLUMBUS等等^[1-15], 广泛用于计算分子的结构、性质以及化学物质之间的相互作用.这些计算通常应用于传统实验手段无法解决的化学和生物学问题, 例如研究和表征生物分子与材料, 研究它们之间在原子、分子、超分子和系统级的相互作用, 提供定性和定量的分析等.

早在计算机发明之前, 人们就已经开始使用计算的方式或量子力学的方法来研究化学问题.随着计算机技术和计算技术的进步, 研究人员逐渐开始依赖并行计算机作为量子化学计算模拟的主要工具.量子化学计算可以在个人电脑、工作站、部门级(T级)并行计算机和超级计算机(P级)等多种平台上进行.目前的个人电脑已经拥有了多个计算单元(多核CPU), 因而并行计算已经成为量子化学计算的主要方式, 绝大多数目前流行的量子化学计算软件都可以采用并行计算.

在计算机发展的早期, 由于计算机存储能力的匮乏, 人们想要计算模拟较大的分子体系, 主要通过改进或发展新的算法来实现, 包括直接积分法^[16]、单元分解法^[17]、从头计算法局域积分^[18]和冻结实近似^[19]等等, 将限制大分子从头计算的内存或硬盘存储能力的不足转换为相对宽松的处理器计算能力的需求.此外, 人们也充分利用特定计算平台的硬件特性来增加量子化学计算的效率.每一次硬件技术的提升, 都促使计算化学家们改进计算理论和算法以充分利用硬件的潜能.并行计算方法最早出现在量子化学领域中是1980年^[20], 在该项研究中, 微程序设计技术被用来实现内存和基本运算操作的并行使用, 并获得了2.5$\sim $2.8倍的性能提升.到了20世纪90年代后期, 单个处理器计算能力的增长变得缓慢了.提高计算机整体计算速度的最主要的方法变成将许多价格低廉的处理器以共享或分布内存的网络联结的方式来实现.人们需要在量子化学程序中采用并行算法以提升计算效率.许多传统的串行量子化学程序重新设计算法并改写了代码, 或者完全重新开发以适应大规模并行计算平台.近年来, 多核和众核处理器的出现要求计算软件进一步实现细粒度的并行, Intel Xeon Phi芯片的出现, 以及图形处理器(GPGPU)及其编程接口的完善推动其在一般科学计算中的广泛应用, 这些都使得量子化学家们再一次改写或编制新量子化学程序以充分利用硬件的性能^[21-22].

目前计算化学研究人员已经能够广泛地使用T级或P级计算系统开展化学领域的研究, 使用的计算规模从数核到数十万核, 计算模拟的体系和精度也越来越复杂.本文将介绍量子化学的计算方法和并行特性, 并对未来量子化学软件的发展和应用进行了展望.

1 计算方法与特性

量子化学计算一般是为了获取原子或分子的电子结构和能量特征, 进而获得相关的物理或化学性质.根据近似程度和能量表述的不同, 可以区分为两类方法.一类为单电子近似方法, 例如Hartree-Fock方法(HF)^[23]和密度泛函理论(DFT)^[24].此类方法通过准确地计算单电子函数, 采用近似的能量表达式来计算总能量.另外一类为后Hartree-Fock方法或相关能方法, 例如组态相互作用方法(CI)^[25]、多体扰动理论(MBPT)^[26]、耦合簇理论(CC)^[26]等.此类方法在Hartree-Fock方法的基础上, 引入激发态slater行列式计算相关能, 通过多个slater行列式的线性组合来实现双电子相互作用的计算.从计算复杂度来说, 假设需要计算的体系的电子数为$N$, 那么在单电子近似方法中, 内存的需求一般为$O(N^{2})$, 计算时间的需求一般为$O(N^{3})$; 而在后Hartree-Fock方法中, 在只考虑单、双激发的条件下, 相应的内存需求和计算时间需求一般分别为$O(N^{4})$和$O(N^{6})$. 图 1显示了主要量子化学计算方法的常用计算体系大小、计算开销和计算精度.

目前在一些量子化学软件中, 也有一些时间复杂度为$O(N)$的线性标度计算方法的实现.这种方法的原理是认为能量的大小与分子的大小是呈线性相关的.而这种近似一般只对非常大的体系有效, 其要求体系中的远程相互作用随着体系大小的增长迅速下降^[27].这种思路已经大量地用于各种量子化学计算方法中, 包括Hartree-Fock、DFT、MBPT和CC等^[28-31], 并且在大量的量子化学软件中得到了实现^{[3-5, 32]}.线性标度计算方法在使用较小规模计算资源计算较大体系时具有显著优势, 但对于大规模并行计算机来说, 线性标度计算方法本身在并行计算实现和扩展性上具有一定困难.主要原因在于线性标度计算方法目前基本上只对计算本身进行线性标度扩展, 而对通讯和数据并没有进行线性标度扩展.而在量子化学并行计算中, 数据I/O和通讯占据了相当的计算时间, 因而线性标度计算方法的并行实现具有相当的数据密集性, 并导致通讯瓶颈.当然, 如果能够实现直接计算远端原子相互作用以减少数据交换和数据存储, 那么并行线性标度计算方法将具有相当高的计算效率和扩展性.就目前来说, 大规模并行线性标度计算的实现还不够成熟.

1.1 HF和DFT方法

HF或自洽场(SCF)方法是所有量子化学计算程序的核心和基础.其他高阶的计算方法包括CI、MBPT和CC等方法都是在SCF方法的基础上进行扩展的. SCF方法是一种求解全同多粒子体系的定态薛定谔方程的近似方法.其近似地用一个平均场来代替其他粒子对任意一个粒子的相互作用, 这个平均场又能用单粒子波函数表示, 从而将多粒子系的薛定谔方程简化成单粒子波函数所满足的非线性方程组来解.对于一个N电子原子来说, 其能量或者薛定谔方程为

$ \begin{align} \Big[-\dfrac{1}{2}\sum\limits^{n}_{i=1}\nabla^{2}_{i}-\sum\limits^{n}_{i=1} \dfrac{Z}{r_i}+\sum\limits^{n}_{i=1}\sum\limits_{i>j}\dfrac{1}{r_{ij}}\Big]\varphi=E \varphi. \end{align} $

(1)

在上式中, 左边第一项是电子的动能, 第二项是电子与原子核的库仑吸引能, 第三项是电子间库仑排斥能, 右边$E$是体系的总能量.在SCF方法中, 每个电子都近似为在其它电子的平均势场中运动, 于是对于每个电子来说, 其薛定谔方程近似为

$ \begin{align} \Big[-\dfrac{1}{2}\nabla^{2}_{i}-\dfrac{Z}{r_i}+ \sum\limits^{n}_{j\neq i}\int\dfrac{\varphi^{2}_{j}}{r_{ij}}{\rm d}\tau_j\Big]\varphi_i=E_{i} \varphi_{i}. \end{align} $

(2)

这样就可以列出$N$个单电子耦合方程.然后用一个波函数的假设解代入方程求出新的波函数, 通过迭代法逐次逼近, 直到前后两次计算结果满足所要求的精度为止(即达到前后自洽).体系的总波函数用一个由各个电子的轨道波函数$\varphi_i$组成的slater行列式来描述^[33].对于分子体系来说, 分子轨道波函数由原子轨道波函数的线性组合来描述, 组合的参数用变分法确定.体系的总波函数同样用一个由各个电子的分子轨道波函数组成的slater行列式来描述.

密度泛函理论(DFT)采用的是另外一种方式来处理多电子问题.由于它既有较为准确的计算结果, 又有较高的计算效率, 因此得到了非常广泛的应用, 尤其是在材料科学领域中处理较大规模的体系.与HF方程类似, DFT的基本方程是Kohn-Sham方程(KS方程)^[24].在DFT中, 多电子相互作用被处理成与密度相关的交换相关项, 一般来说, 这个交换相关项是非局域化的.由于KS方程和HF方程在数学上具有相似性, 因此也可以使用自洽的方式来求解. DFT中的关键问题就是处理这个交换相关项.由于交换相关项的密度表示的确切形式是未知的, 人们发展了很多近似方法, 包括纯密度形式的局域密度近似LDA和广义梯度近似GGA (如PW91、BLYP和PBE96等), 以及更复杂的杂化形式(如B3LYP、PBE0和CAM-B3LYP等).杂化形式包含了HF的双电子交换贡献, 应用更为广泛.由于HF和DFT的相似性, 这两种方法在程序实现上会共用许多核心代码. HF和DFT方程原则上都可以数值求解, 但实际上都是用基组展开的方式求解.人们已经发展了大量的基组方法, 这些方法可以分为两类, 一类是局域基组方法, 另一类是平面波方法, 两类方法各有优缺点.

在局域基组方法中, HF和KS方程可以用一个矩阵方程来表示^[34],

$ \begin{align} FC=SC\varepsilon. \end{align} $

(3)

其中$F$为FOCK矩阵, $C$为系数矩阵, $S$为重叠矩阵, $\varepsilon $为对角化轨道能量矩阵. FOCK矩阵$F$可以表示为单电子部分、双电子库仑相互作用、双电子交换贡献(HF方法和DFT方法)以及相关能部分(DFT方法)之和.给定基组后, 单电子部分很容易计算, 而双电子部分依赖于一系列的双电子积分的密度矩阵, 必须使用迭代的方式求解.因此, 在实际计算程序中, 计算耗时的部分主要包括:双电子积分计算、FOCK矩阵的双电子部分构造、DFT中交换相关贡献的计算、FOCK矩阵对角化以及使用分子轨道系数的密度矩阵的计算.双电子积分的计算是最消耗资源的, 人们发展了很多办法来在提高计算效率^[35], 包括双电子积分赋值法和单电子密度矩阵约化积分.由于双电子积分只是原子位置和选用基组的函数, 因此在求解HF方程或DFT方程中, 可以作为常数来处理.一般来说, 可以一次性计算出来并存储到硬盘上, 这种技术称之为"传统FOCK矩阵构造"法.尽管存储空间和I/O带宽需求可以通过去除精度阈值以下的积分来减少, 但在实际计算中, 这种积分文件常常非常大(约几百个GB), 会造成计算瓶颈.另外一种方法是使用时按需计算双电子积分(on-the-fly), 称之为"直接FOCK矩阵构造"法^[16].这类方法的缺点是计算部分的耗时增加了.考虑到CPU速度和I/O速度的相对快慢, 实际上这种方法适合于大型系统, 并能获得较快的计算速度.半直接方法结合了上述两种方法的特点^[36], 将一小部分比较重要的积分和密度矩阵元一次性计算出来存储到硬盘上, 其他部分在使用时按需直接计算.此外, 在上述三种方法的基础上, 还有许多其他的方法用于提高双电子积分的计算效率, 例如使用辅助基组实现$O(M^{2})$级算法($M$为辅助基组的数目)、多极近似法和大规模积分预筛选等.具体的实现包括单元分解法(RI)^[17]、Cholesky分解(CD)^[37]、电荷密度拟合法^[38]和有效FOCK矩阵构造法^[29]等等.实际上, 积分算法的选择与所计算的体系是密切相关的.例如, SCF迭代收敛的步数会直接影响传统FOCK矩阵法和直接FOCK矩阵法的实际效果.半直接FOCK矩阵法对于中等规模的体系最为有效, 而对于较大的分子, 直接FOCK矩阵法更有优势.对于使用上千个基组的大体系, 由于磁盘空间的限制, 基本上是不会采用传统FOCK矩阵法进行计算.对于这些体系, 通过利用消除四中心积分和有效积分预筛选技术, RI和CD方法比直接FOCK矩阵法更有效率.在这些方法实际计算中, 计算耗时部分的效率主要依赖于计算的体系、筛选技术对于给定体系的效果以及数据直接读取再利用和重计算之间的比率关系等.此外, 要获得好的扩展性, 计算耗时部分都需要进行并行化.

SCF方法和DFT方法的并行化主要集中在FOCK矩阵对角化、DFT交换相关贡献和双电子贡献部分.对于大分子来说, 上述3个部分的时间复杂度分别是$O(N^{3})$、$O(N^{3})$和$O(N^{2})$, $N$为基组函数的数目. FOCK矩阵对角化的并行效率会随着矩阵大小的增加而增加, 但对于大规模处理器计算, FOCK矩阵对角化的并行比较困难, 因而人们开发了一些扩展性较好的高效免对角化算法^[39].此外, 传统的HF/DFT实现都是利用并行线性代数库例如ScaLAPACK来进行矩阵操作. DFT交换相关贡献主要涉及到在原子网格格点求积分, 而这些积分是相互独立的.因此一般的做法是将这些网格格点作为独立的批处理任务在处理器中进行分发, 然后广播密度矩阵^[40].双电子贡献的并行化主要有两种方式, 数据副本方式和数据分发方式或者两者混合^[41-43].数据副本方式简单明了, 通讯量小.在这种方式下, 每个处理器都保留一份必要的数据, 也就是说每个处理器都保留一份密度矩阵和FOCK矩阵.每个处理器都对分块的四重积分进行计算, 循环主体为这些四重积分.每完成一个四重积分, 该处理器都将结果整合到它所负责的部分FOCK矩阵.待所有处理器都完成任务后, 再全局整合, 构成完整的FOCK矩阵.如果工作负载分配均衡, 全局整合算法高效, 这种方式就具有很高的并行效率.对大的计算体系来说, 由于工作负载较多, 特别是如果采用直接FOCK矩阵法, 可以获得很高的计算效率.这种方法的缺点在于内存的使用.由于每个处理器都要保留数据副本, 因此内存的使用将是$O(N^{2})$级($N$为使用的基组数目), 这对每个处理器可使用的内存有一定要求.数据分发方式是解决数据副本方式内存瓶颈的一个有效方式.在这种方式下, 每个处理器只处理一部分数据因而对本地内存的需求较小.人们已经发展了大量的数据分发方式的算法^[42-43], 这些算法基本上是通过有效再利用数据来抵消对密度或FOCK矩阵获取的消耗, 本质上是先对基组函数进行分组, 然后基于这个分组再将密度或FOCK矩阵元素分配到相应的处理器进行计算.这种方式下, 通讯的消耗是平方级的, 计算的消耗是四次方级的, 因此通过适当设置分组的大小, 可以使得计算成为主要耗时部分.另一方面, 由于数据分发方式中, 每个处理器只处理一部分数据, 因此需要从其他处理器处获得数据来评估本地数据对整个FOCK矩阵的贡献.这需要人们小心处理通讯和计算, 避免死锁或处理器空闲.此外数据分发方式中的负载均衡也需要仔细处理, 简单地将积分平均分配到处理器上会导致较差的负载均衡.这主要是因为双电子积分计算的耗时与角动量以及基组函数的收缩密切相关.人们已经发展了一些方法来处理这个问题^[42-43].数据分发方式具有很好的扩展性, 已在不少程序中得到了应用^[44].

1.2 MBPT方法

多体扰动理论(MBPT)是重要的量子化学相关能计算方法之一.所谓相关能指的是体系的实际能量与独立粒子模型(即HF计算)获得的能量之差.虽然相关能占体系总能量不到1%, 但对准确描述化学反应过程至关重要.相关能或相关效应一般是由电子的瞬间相互作用产生的. MBPT就是用来精确求解薛定谔方程以更为准确的获得体系能量的方法之一. MBPT根据扰动阶数的不同可以分为二阶微扰处理、三阶微扰处理、四阶微扰处理等等, 阶数越高, 计算结果越准确.但在实际应用中, 由于计算复杂度的增加, MBPT的应用一般只限于低阶微扰处理.例如MBPT2和MBPT3方法计算能量的时间复杂度分别为$O(N^{5})$和$O(N^{6})$, $N$为计算体系的大小.因而在日常计算中MBPT2或MP2方法^[45]使用地最为广泛.

在过去的几十年里, 人们做了许多工作使得MP2能够应用到较大的分子体系, 包括拉普拉斯变换法消除能量分母^[46]、双电子积分分解的直接实现(RI-MP2等)^[47]和局域MBPT2方法^[18]等. MBPT2方法的并行实现需要仔细处理通讯和内存的消耗.最早的数据分发并行模式的MP2算法出现在1995年^[48], 但这个方法具有较高的通讯和内存消耗. NWChem基于数据分发模型和Global Array环境(GA)也实现了MP2和RI-MP2的并行^[49-50], 该实现属于计算密集型, 可以扩展到较多的处理器.另外一种MBPT2的并行实现采用了直接积分数据分发模型, 并获得了超线性加速性能^[51].其他通讯密集型^[52]和I/O密集型MP2并行算法^[53]的实现也可以满足涉及几千个基组函数的化学体系的计算.近期的一些MBPT2并行实现的进展包括可以扩展到32核, 计算4 000个基组函数的RI-MP2方法^[54]; 使用MPI通讯, 并可以扩展到数百核的局域MP2并行实现(LMP2)^[30]; 可以在一个处理器上计算由一万个基组函数描述的含一千个原子体系的线性标度原子轨道MP2方法及其并行实现^[55]等等.

1.3 CC方法

耦合簇(CC)方法是另外一种广泛使用的相关能计算方法, 其与MBPT方法具有一定的关联性. CC方法可以准确地描述多电子体系的相关效应, 包括弱相关体系和强相关体系.根据CC簇算符激发阶数近似的不同, 可以分为CCSD方法(单、双激发CC方法)、CCSDT (单、双、三激发CC方法)和CCSDTQ (单、双、三、四激发CC方法)等.高阶的CC方法意味着算法复杂性和计算消耗的快速增加.例如CCSD的计算复杂性为$O(N^{6})$, 而CCSDT的计算复杂性则为$O(N^{8})$, $N$为计算体系的大小.因此, 为了减少计算的复杂性, 人们发展了一些非迭代性的高阶CC方法, 以平衡计算精度与计算代价.例如人们在日常计算中经常使用计算复杂性为$O(N^{7})$的CCSD(T)方法^[56].对于激发态问题, CC方法分为线性响应CC方法(LR-CC)^[57]和运动方程CC方法(EOMCC)^[58]两类处理方式.与基态类似, 人们也发展了各种迭代性(EOMCCSD^[58]、EOMCCSDT^[59]、EOM-SF-CCSD^[60]等等)和非迭代性的CC方法(EOMCCSD(T)^[61]、CCSDR(T)^[62]、CR-EOMCCSD(T)^[63]等等)处理激发态.

与其他方法相比, CC代码的并行化具有相当的困难.早期人们也对CC代码的并行化做了一些工作^[64], 目前能够提供CC并行计算的量子化学计算程序主要包括NWChem、MOLPRO、GAMESS、ACES Ⅲ和PQS等. NWChem中的CC代码的并行实现主要分为两类, 适用于闭壳层体系的非自旋代码和使用自旋轨道的TCE实现.闭壳层CC代码是基于GA库实现的, 可以计算CCSD和CCSD(T)方法, 对处理器之间互联通讯的性能要求较高. TCE产生的CC代码可广泛地使用与限制性HF(RHF)、非限制性HF(UHF)、限制性开壳层(ROHF)、DFT模型等有关的参考态函数, 其扩展性可达到数万核, 计算的体系可使用超过一千个基组函数. MOLPRO中的CC代码也是基于GA库实现了基于RHF和ROHF的CCSD和CCSD(T)方法的并行计算. GAMESS-US^[7]的CC代码利用的是分布式数据界面(DDI)实现了CCSD、CCSD(T)和CR-CCSD(T)的并行计算, 其扩展性可达到数十个核, 计算的体系可使用数百个基组. ACES Ⅲ使用超级指令架构语言(Super Instruction Architecture Language, SIAL)来编写CC并行代码.目前已实现了基于RHF和UHF参考态函数的CCSD、CCSD(T)和EOMCCSD方法, 其并行规模可达数千核. ASES Ⅲ的设计分为两个领域, 一个领域是使用SIAL控制电子结构方法的理论部分, 而不需要控制计算细节, 另一个领域是使用超级指令处理器(Super Instruction Processor, SIP)执行、解释SIAL代码, 并为给定的计算机架构提供必要优化和调试. PQS实现的闭壳层CCSD和CCSD(T)并行计算主要是针对中等规模的计算机群, 计算的并行规模在数十核左右.

1.4 MCSCF方法

多构型SCF方法(MCSCF)的波函数是自旋匹配构型态函数(CSF)的线性展开, CSF的展开系数与分子轨道系数同时优化. CSF空间常常通过将初始的$N$个轨道划分为具有占据限制的轨道子空间来确定.例如一个最简单的CSF空间包括非活跃轨道、活跃轨道和虚轨道, 其中非活跃轨道是双电子占据, 虚轨道没有电子占据, 活跃电子在遵守对称和自旋限制的前提下以所有可能的方式占据活跃轨道并组成了全活跃空间(CAS).实际上CAS等价于在活跃轨道上的全CI展开.因此, 一个给定自旋多重态的CSF数随着活跃电子和活跃轨道的增加呈指数性增长.如果用行列式的形式来表示, 则数目更为庞大.人们通过将活跃空间划分为若干个具有不同占据和自旋耦合限制的子空间, 例如限制性活跃空间SCF方法(RASSCF)^[65], 来减小CSF空间.此外, 直积空间法^[66]也可以简化波函数展开.

在MCSCF计算中, 最耗时的部分包括: ①原子轨道---分子轨道(AO-MO)积分变换; ②组态相互作用(CI)本征值问题包括密度矩阵; ③求解Newton-Raphson方程(NR方程).对于第一部分, 数据副本模式实现简单, 但会导致内存和I/O负载过高.但如果有高效的并行积分代码, AO积分可以在使用时计算, 通过增加CPU时间来减少I/O负载, 同时产生的MO和中间量可以保存在分布式内存中.此外, 通过利用辅助基组扩展的方式^[67], 可以显著减少AO-MO变换计算消耗.对于第二部分, 将波函数限制在CASSCF比使用一般性的MCSCF波函数更有效率.此外, 将波函数在CAS和RAS上用行列式展开^[68]更有利于减少计算消耗.对于第三部分, 一般以亚空间展开的形式迭代求解NR方程.在亚空间展开只需要计算部分矩阵矢量与猜测矢量的积, 这样就不需要存储所有的大型矩阵, 以免在NR迭代循环中增加额外的计算消耗.

目前, 已有不少MCSCF方法的并行实现^[69], 例如GMAESS使用并行MCSCF方法在128CPU上并行, 实现了58%的效率^[8].但是, 由于CSF空间数目的指数性增长, CASSCF方法的扩展性能和计算效率已接近极限.减少昂贵的四坐标变换的数目, 同时保留二阶收敛的并行实现, 以及限制性组态相互作用展开将是可行的解决方案.

1.5 CI方法

组态相互作用CI方法在概念上有不同的实现方式, 从而影响其并行实现方式和应用方式.在CI方法中, 波函数的主要形式是CSF的线性展开或行列式的线性展开.由于一般情况下, 人们只需要最低的几个本征值, CI方法可以在一些较小的子空间上进行迭代, 从而避免存储较大的CI稀疏矩阵, 计算消耗可以减少一个级别. CI方法中每一次迭代循环的计算处理过程一般分为三步: ①组成波函数矢量; ②计算哈密顿量的子空间表示和重叠矩阵; ③求解非正交子空间本征值问题.组成波函数矢量是最耗时的部分, 其耦合系数是一个酉阵的矩阵元素, 这个矩阵非常稀疏, 一般在使用时直接计算出来.波函数矢量的计算需要遍历所有的积分, 主要的计算包括矩阵和矩阵乘、矩阵和矢量乘以及矢量操作.第二和第三步对数据的吞吐要求比较大.如果将数据保留在内存中, 将极大地减少第二和第三步的通讯和I/O的压力, 但内存的消耗会增加很多.

组成波函数矢量的并行实现具有相当的困难, 需要同时考虑负载均衡、带宽和内存限制.在得到耦合系数的前提下, 波函数矢量通过遍历所有积分来计算, 并可随机读取.这需要将积分和波函数矢量保存在本地内存中, 并可通过全局操作求和.这种方法严重受限于CI方法的维度和本地内存的容量, 但很容易实现动态的负载均衡.另外一种方式是将CI矩阵分区为小块, 同时波函数矢量也相应的分成小段.矩阵矢量乘可以分解成独立的部分矩阵矢量乘(独立的任务), 并可随机读取相应的小段波函数矢量.为了减少数据传输量随着处理器的增加而增加, 波函数矢量可以保存在分布式内存中, 而积分可在本地内存保留副本.每一个独立的部分矩阵矢量乘需要拷贝波函数矢量到本地内存, 遍历积分子集来计算部分波函数矢量, 同时更新整个波函数矢量. CI矩阵块计算的消耗变化范围非常大.工作负载的空间分布很大程度上取决于CSF空间的选择, 并且反映了CI矩阵的稀疏性.因而, 如果不使用性能模型, 平衡内存消耗、通讯负载以及工作负载将会非常困难.足够灵活的任务分配机制和合理准确的性能模型对在当前并行计算机上进行CI计算非常重要.例如, 在CULUMBUS软件中, 其任务分配机制使得CI方法计算在1~024核的时候都具有非常理想的加速比^[70].在该软件中, 每个CPU处理的任务一般不超过4$\sim $8个, 当使用CPU数少的时候, 内存的限制导致通讯负载的增加, 当使用大规模CPU的时候, 工作负载均衡限制导致通讯量增加.

此外, I/O是制约CI大规模计算的并行效率的瓶颈之一, 尽管在获取非本地数据的时候采用压缩技术可以提高有效带宽.例如, 根据文献[71]报道, 一个基于行列式的MRCI并行代码, 在带本地盘的四核Linux机群系统上实现了16亿个行列式的计算.该实现在16个节点内, 每个节点使用一个核时, 都具有良好的加速比, 加速性能达到0.88以上.而如果增加核数, 则加速性能迅速下降到0.3.尽管使用了本地盘, I/O仍然是一个主要瓶颈, 如果使用并行文件系统, 则性能会更差.

目前, CULUMBUS软件是唯一有能力使用上千核实现CI方法的高效并行计算. I/O、负载均衡等困难导致目前各种CI方法的实现缺乏可扩展性.人们还需进一步修改计算模型, 以获得更高的可扩展性.

2 结语

随着高性能计算技术的进步, 尤其是协处理器技术(如GPGPU、Intel Xeon Phi等)的进步, T级的计算机已从十多年前的稀缺性发展到目前的桌面平台, P级系统的计算资源也已逐渐丰富, 未来几年E级系统的出现将可能导致P级系统更加大众化.这些都为大幅度提高量子化学计算模拟的广度和精度创造了机会.计算化学家们已逐渐可以使用更准确同时也更消耗计算资源的计算方法来进行更大体系的量子力学计算模拟.当前为高端设备研制的量子化学计算程序将在未来成为日常设备使用的计算工具.同时, 当前的高性能计算平台已不再是单一的CPU计算平台, 而是由CPU和协处理器共同组成的混合平台, 未来的E级计算机将可能包含百万以上的混合计算单元, 互联设备也将可能采用混合结构.人们需要开发新的量子化学计算模型和算法, 调整量子化学软件和设计思路, 充分利用硬件的特性, 进一步提高量子化学软件的计算效率和并行规模.一些可行的思路包括: ①修改计算模型和算法.借鉴结构力学或流体力学模拟中的有限元分析的方式, 将大的研究体系网格化或碎片化, 进行分布式量子力学计算, 同时处理好边界关系. ②采用混合编程模型.考虑到未来的CPU和协处理器均含有大量的核处理单元, 采用多线程或openmp与mpi混合编程模型将能更好的利用这些设备, 提高计算效率.目前广泛使用的量子化学计算软件如Gaussian (openmp并行)、NWChem (mpi并行)、Q-Chem (mpi并行)等尚未实现这一点. ③增强软件的容错性和健壮性.未来E级平台计算时涉及的设备单元非常多, 平均无故障时间将会减少, 因而软件的容错性和健壮性将会被进一步考验.量子化学软件将需要更细致的考虑容错机制, 如线程级的检查点、回滚、进程迁移等, 减少由于硬件故障导致的计算失败, 或者尽快从硬件故障中恢复之前的计算.这也是目前大多数量子化学计算软件所欠缺的.