为了快速地计算某Ewald sum，现推导一种基于矩阵的计算策略，一种很新的算法。在时间复杂度方面本策略没有什么提升，还是。但是，比起传统方法，本策略的FLOPs大幅降低，并且支持cuda加速。
本文涉及的矩阵运算
向量的数量积
一个维的向量实际上就是大小为的矩阵，向量点乘可以定义为一个行向量与列向量的乘积。设，，则与的数量积
矩阵乘法
设是一个的矩阵，是一个的矩阵，矩阵乘法定义为 其中第行列的项的值为。
矩阵的Hadamard积
设，均为的矩阵，A与B每个相同位置元素的乘积被称为Hadamard积，记作，矩阵之间的Hadamard积的逆运算，即矩阵元素分别做除法，记作。矩阵A与B的Hadamard积
矩阵的数值平方
设矩阵是一个维的向量，定义矩阵的一种运算 ……就称为矩阵的数值平方好了。其中均为维度为N的行向量。
矩阵求和
对于大小为矩阵，规定表示矩阵所有元素的和，即。
矩阵上的函数
对于大小为的矩阵，如果其所有元素均在函数的定义域内，则
厄米矩阵与共轭转置
如果复矩阵的共轭转置等于它自己，即，则称是一个厄米矩阵。
现证明厄米矩阵的对角元素组成的矢量等于的实部的对角元素组成的矢量。
首先将厄米矩阵分解为一个实矩阵和一个只有纯虚数和0的矩阵的和 其中为的实部，为虚部。厄米矩阵的对角元素必须等于自己的复共轭，因此都是实数。所以是一个反对称矩阵，满足。厄米矩阵的对角元素组成的矢量可以拆成两个向量的和： 因为矩阵的对角元素都是0，所以。因此有 其中表示复矩阵的虚部。
引理1 现有大小相等的实对称矩阵与厄米矩阵，则。
证明
首先矩阵乘以矩阵可以分解为 根据矩阵的迹的性质有 因为反对称矩阵满足，从而有 两边取矩阵的迹得 联立得。代入式(1)得 证毕。
Ewald summation
本章中的符号与运算含义参见附录。
通常的Ewald sum由三部分组成： 其中为能量的短程部分，为能量的长程部分，为自能修正项。下面分别介绍三个能量项的矩阵求法。
能量短程部分的矩阵计算方法
能量的短程部分是一个与点的坐标和电荷相关的量，其定义为 其中表示第个点的电荷，为点的坐标，为两个点的坐标之差，表示对使用周期性边界条件化归之后求得的模的平方。
为了使用矩阵加速计算速度，需要将式(5)中的求和符号转换成矩阵乘法。构造电荷耦合矩阵表示每两个电荷的乘积： 其中表示所有点处的电荷。
矩阵表示每两个点之间的相互作用 其中为狄拉克函数，当且仅当时等于1，否则等于0。
能量的短程部分最终可表示为
能量长程部分的矩阵计算方法
能量的长程部分的计算方法如下： 其中是波矢，和均为关于波矢的函数。和的计算公式分别为： 其中是一个参数，是第个点处的电荷，是第个点处的位移。对于每一个波矢，其对应的和都是标量。
为了使用MATLAB快速高效地计算，我们试图将向量运算推广到矩阵维度。如果用矩阵表示所有的波矢，则 矩阵，其中为的数量。此时是一个列向量，每一个由式(8)代入计算得到。因此，对于矩阵，有 其中表示Hadamard积的逆运算，即对待操作的矩阵的每个元素做除法的运算。
与矩阵类似地，我们用矩阵表示所有点的坐标，的矩阵，其中为的数量。此时对于每一个波矢有 是一个数量积。也就是说，每一个计算得到的是一个数字。令是一个维的行向量，则有 将式(10)推广到矩阵上有 其中表示取矩阵对角元素构成的列向量，矩阵满足
于是能量的长程部分可以表示为 其中表示的转置，是一个行向量。为了尽可能将矩阵分离出来，构造对角矩阵，其中是一个单位矩阵。可以将式(11)变换为 其中表示对矩阵的所有对角元素求和，也就是矩阵的迹。
自能项的矩阵计算方法
如果用矢量表示所有的点的电荷，则自能可以表示为。又因为电荷，所以自能的计算公式为
后续用于系统模拟的缓存优化
当空间中的点的坐标不改变时，能量的短程部分仅与电荷有关。对式(7)中的稍作变形可以与前面的常数项合并，从而减少浮点数计算次数：
类似地，能量的长程部分也仅与电荷有关，利用矩阵的迹的性质对式(12)变形，整理得到 令，显然是一个厄米矩阵。 引理1已经证过，是一个实数，其中表示取的实部。与此同时，注意到自能项也可以用矩阵的迹表示，于是与可以合并。对式(13)变形得到 其中是单位矩阵。
于是系统的总能量可以表示为 设矩阵称为邻接矩阵，矩阵称为系数矩阵，代入上式得 其中邻接矩阵满足 系数矩阵的计算公式为
实验发现，使用MATLAB编写代码计算Ewald sum的时候，计算的时间比起计算的时间显著地多。原因是MATLAB在计算的时候，先计算矩阵的乘积然后再求对角的和，时间复杂度自然要多一个数量级（之于）。为了避免不必要的计算，现从矩阵运算的角度优化。观察发现，
于是系统的总能量可以表示为 其中矩阵表征系统的除电荷以外的所有信息，不妨称为Ewald矩阵。实际上，Ewald sum可以定义为是一个二次型，即
结论
对RPM系统进行模拟时，实际上只需要改变各个点处的电荷，而Ewald矩阵不需要改变。因此，如果预先将Ewald矩阵预先计算出来，则每次模拟时系统的总能量的求解等同于式(16)的值的求解。该式的时间复杂度为，浮点数运算次数（包括乘法和加法）为(GFlops)。
分别使用MATLAB和C++编程（调用OpenBLAS库进行矩阵计算）对的三维空间内有72对正负带电粒子的系统进行模拟。计算所使用的服务器有40个CPU，型号为Intel(R) Xeon(R) Silver 4316 CPU @ 2.30GHz。 MATLAB计算速度约为s，C++的计算速度约为s。
结果表明，将公式里的向量运算和求和改写成矩阵运算，不仅可以有效地减少重复计算，而且还能将电荷与其他矩阵分离，从而可以预先计算Ewald矩阵，进一步加快运算速度，并且使GPU加速变为可能。
Fin.
附录：符号约定
表1: Ewald summation章节里用到的符号及其含义表 | 总能量 |
| 能量的短程部分 |
| 能量的长程部分 |
| 自能修正 |
| 边界条件内某一维度的大小 |
| 反正是个参数 |
| 总粒子数 |
| 总波矢数 |
| 带电粒子个数 |
| 某特定粒子的坐标，矢量 |
| 动量空间里的一个波矢，矢量 |
| 所有粒子电荷，矢量 |
| 表示所有粒子坐标的矩阵 |
| 表示所有波矢的矩阵 |
| 描述每两个电荷之间的耦合的矩阵 |
| 向量与向量的数量积 |
| 尺寸灵活的单位矩阵 |
| 复矩阵的模 |
| 矩阵的转置 |
| 复矩阵的共轭转置 |
| 矩阵与的乘积 |
| 矩阵与的Hadamard积 |
| 矩阵与的Hadamard积的逆运算 |
| 矩阵的数值平方 |