Compare Ewald with NIST reference

[Luthaf]

And fix a bug along the way \o/

We do not reproduce very well the 4th NIST configuration results, I'll check if this is improved by changing the Ewald parameters (alpha & k_max) to use the same values as the NIST reference.

[g-bauer]

Awesome! I am currently running tests with SPC/E water + Ewald for NVT and NVE to compare to the NIST phase equilibria data and I recognized that pressure is doing weird things, your find might explain why.

I was wondering, what is the reasoning behind

alpha: 3.0 * PI / (cutoff * 4.0)

I used the value given by NIST

α = 5.6 / min(Lx,Ly,Lz)

[Luthaf]

alpha: 3.0 * PI / (cutoff * 4.0)

I am afraid I do not remember where this come from. I would like to improve this by allowing the user to only request a precision (say 1e-5), and then compute alpha and kmax to get this precision and enforce the O(n^3/2) behaviour of Ewald.

[g-bauer]

Maybe use kmax and tolerance as inputs with alpha set just like in the NIST reference.

[g-bauer]

Can we get the test to pass with the values from NIST?

I would like to improve this by allowing the user to only request a precision (say 1e-5), and then compute alpha and kmax to get this precision and enforce the O(n^3/2) behaviour of Ewald.

Is that something you want to add in this PR?

[Luthaf]

Can we get the test to pass with the values from NIST?

No, I broke the energy conservation test with these changes, I need to fix it first.

Is that something you want to add in this PR?

No. I'll add it later, in a separated PR =)

[Luthaf]

I finally found the time to finish this. I added more tests, the Ewald code was broken everywhere.

Now it should be better, but please review carefully the tests, as they where not testing what I though before.

[g-bauer]

Still working through your changes. I didn't implement Ewald summation myself yet so I am very slow.
I will very naively put some questions in here, where I don't understand the equations. I hope that's ok.

You mention this in the tests in the comments. Could the large deviation in the 4th sample come from the different alpha values?

• NIST: alpha = 5.6 / 30.0 = 0.1867 (4th example)
• Lumol: alpha = 3*pi / (9.0 * 4.0) = 0.2618 (4th example, cutoff = 9 A)
  The difference is smaller for the smaller boxes of examples 1-3.

[Luthaf]

I hope that's ok.

That's totally ok!

Could the large deviation in the 4th sample come from the different alpha values?

Yes, it does. In the last commit, I added a method to set the value of alpha, and used that for nist-spce tests. This make the relative difference in energy two order of magnitude smaller, at the cost of making the 9A/nist-3 test a bit worst (from 1e-5 to 1e-3).

[Luthaf]

Ping @g-bauer. Again, if you do not have time, I can do a review of this myself. It is old enough that I forgot about how I've done it 😃.

[g-bauer]

Sorry for the delay. I haven't had the time to go through it in one iteration. I think it is better if you do the review. I will continue to go through it, but I don't know when I'll have the time.

The tests passing are a good sign that the implementation is correct. We could further check the different contributions separately, like in the NIST benchmark, section 6).

[Luthaf]

I re-read this and rebased it on top of master. I am waiting for Travis to go green, and then I'll merge the PR.

We could further check the different contributions separately, like in the NIST benchmark, section 6).

We can not do this, because the corresponding functions are private. We could make them public, but this seems very ad-hoc, and I do not see how users of the code would need these functions.

[g-bauer]

That's fine. It was just an idea to further validate the code but if it passes the tests it should be fine.

I am waiting for Travis to go green, and then I'll merge the PR.

🎉