andysim / helpme

helPME: an efficient library for particle mesh Ewald
BSD 3-Clause "New" or "Revised" License
27 stars 13 forks source link

mu-q and mu-mu energy #57

Closed tgakay closed 2 years ago

tgakay commented 3 years ago

Hi Andy,

I'm running into energy dependence on the kappa parameter when dipoles are present. I've noticed that helPME does not have kappasweep tests for this case so I'm wondering if you've verified it. When testing, I did add the E_self(mu) but it didn't help much. I also tried with a wide range of k-space grids and with a very large real-space cutoff, and both with zero and non-zero total dipole moment to no avail. Below is a table of one of the test runs for a water dimer from your unit tests (including coordinates, charges and dipole moments) with different kappas (called alpha). All quantities are in atomic units (Hartree and Bohr).

grid dim alpha [a.u.] GridSp [a.u.] range [a.u.] E_self [a.u.] E_real [a.u.] E_rec [a.u.] E_total [a.u.] T(real) [s] T(rec) [s]
20 0.02 2 283.936361 -0.023545708 -0.393121622 1E-10 -0.41666733 0.005952 0.003026
30 0.03 1.333333 189.290907 -0.035318903 -0.381349666 1.3161E-06 -0.416667253 0.002008 0.001949
48 0.045 0.888889 126.193938 -0.052979502 -0.363781306 8.88422E-05 -0.416671965 0.001029 0.004987
64 0.0675 0.592593 84.1292922 -0.079473126 -0.338054364 0.000779148 -0.416748342 0.001016 0.011076
96 0.10125 0.395062 56.0861948 -0.119222761 -0.300991047 0.003060132 -0.417153676 0 0.039894
144 0.151875 0.263374 37.3907965 -0.178878262 -0.248487195 0.008947323 -0.418418134 0 0.135686
216 0.2278125 0.175583 24.9271977 -0.268466299 -0.176956275 0.023667676 -0.421754898 0 0.420907
324 0.34171875 0.117055 16.6181318 -0.403202007 -0.093209204 0.065818528 -0.430592682 0 1.399255
486 0.512578125 0.078037 11.0787545 -0.606499142 -0.026404418 0.182708722 -0.450194838 0 4.846253
750 0.768867188 0.052025 7.38583635 -0.91547316 -0.001493128 0.440175597 -0.476790692 0 18.04084

I would appreciate if you'd let me know if I'm missing something and maybe add a kappasweep test with dipoles.

Kind regards, Alexei

Edit: the GridSp is approximate. It's an input parameter used to compute the grid dimensions. Edit2: the splineOrder was 7 in all tests but I've also run with higher orders with the same results.

tgakay commented 3 years ago

A small update. I tried to reproduce numbers for the "Two dipoles" case in table 2 from https://doi.org/10.1063/1.481216 (E(eps_sur=Inf)=-2.0087525) and could only do that for a very small kappa=0.05. So it's obviously not the surface term. Alexei

andysim commented 3 years ago

Hi, Alexei,

Thanks for the very thorough debugging effort. I don't directly have a kappa (alpha, also beta in the PME literature!) sweep test for the dipole case, because I don't support the real space term in the code; I wanted to avoid the slippery slope of becoming a full MD engine, rather than just a reciprocal space helper. I can mimic a dipole with a set of equal and opposite charges and make a test case out of that (I already made something similar to test the electric field last week) so I'll let you know what I find from that. The dipole code itself was tested via a kappa sweep test in the code that inspired me to add multipoles, but I will check it again here. I'm out of town for a few days, so apologies in advance if you don't hear back from me immediately. I am very eager to get this resolved though, so that I can be sure the code is doing what I expect and so that you can use it for what your needs.

Thanks again,

Andy.

andysim commented 3 years ago

Sorry for the long delay in addressing this - it completely fell off my radar for a while. The dipoles should hopefully be fixed by #58 now, but please ping back if there are any further problems. I'm still a bit confused by the test that compares fields to finite differences of potentials; the test is great for cubic boxes, but not otherwise, so I still need to do some reading. Possible h=0 terms aside, the "kappa sweep" test and comparison to the literature values you cited above make me think that the current implementation is correct.

tgakay commented 3 years ago

Hi Andy,

Thank you for the fix and for the heads up. I finally came around to testing it and can confirm that it fixes the issues I reported. I've now incorporated it into the AMS ForceField engine.

Thank you again and kind regards,

Alexei

On 25/06/2021 23:17, Andy Simmonett wrote:

Sorry for the long delay in addressing this - it completely fell off my radar for a while. The dipoles should hopefully be fixed by #58 https://github.com/andysim/helpme/pull/58 now, but please ping back if there are any further problems. I'm still a bit confused by the test that compares fields to finite differences of potentials; the test is great for cubic boxes, but not otherwise, so I still need to do some reading. Possible h=0 terms aside, the "kappa sweep" test and comparison to the literature values you cited above make me think that the current implementation is correct.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/andysim/helpme/issues/57#issuecomment-868837890, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG3TYKHRVCXJXAMOMX2U2UTTUTW5DANCNFSM4WT6CAJQ.