@jmbr thanks for implementing SPME, I discovered this repository from http://superadditive.com/notes/particle-mesh-ewald/.

I am trying to test it on a known problem. I took a unit cell as a box with L = 2/sqrt(3) = 1.1547005383792517 a.u. (Hartree atomic units). I have two unit charges (i.e. q = 1) at a corner and a center of the unit cell, i.e. positions [0, 0, 0] and [1/sqrt(3), 1/sqrt(3), 1/sqrt(3)]. The reference Ewald energy for this problem is E =-1.575834308 a.u..

I am trying to reproduce this using your code. So I created the following *.mat files:

$ cat cell.mat 
1.1547005383792517,0,0
0,1.1547005383792517,0
0,0,1.1547005383792517
$ cat positions.mat 
0,0.57735026918962584
0,0.57735026918962584
0,0.57735026918962584
$ cat charges.mat 
1
1

And run it like this:

$ ./run-pme 64 cell.mat positions.mat charges.mat 10 1e-10
L = 
   1.1547        0        0
        0   1.1547        0
        0        0   1.1547
Read 2 particles into memory.
Real space cut-off: 10
Real space tolerance: 1e-10

Starting non-linear solver with beta = 0.338719
Solver stopped after 7241 iterations.
Using beta = 0.210558 with value = 9.99998e-11
PME initialization took 1.80837 seconds.

energy: Erecip = 0.0000000000000000 vs. 0.0000000000000000
PME took 0.0454914 seconds.

Real space computation over one shell took 1.7998e-05 seconds.
Ecutoff      = 0.9999999999999993
Etotal (PME) = 1.1244772187883054
Edirect + Erecip + Eextra - Eself = 0.9995941893404763 + 0.0000000000000000 + 0.1253333598041780 - 0.0004503303563489

And your code obtained Etotal (PME) = 1.1244772187883054, which is incorrect. I am trying to figure out how to make your code converge to the correct answer. Here are some possible issues that might be causing it:

1. What units does your code use? I use Hartree atomic units above.
2. What cut off should I use? How do I converge against a cutoff? By increasing it?
3. What tolerance should I use? Do I converge against a tolerance by decreasing it?
4. Why does Erecip say 0?
5. I used 64 PW in each direction. I will converge by increasing PW, that much I understand.

Many thanks for any tips.