BUG: DielectricCylinder is sensitive to bin_width

[hejamu]

Summary

The result of DIelectricCylinder for the radial component of the dielectric constant is wrong when using bindwidths larger than 0.1Å

Expected behavior

The quality of the profile should be dependent on the size of the bins, but the result in the bulk region should not be.

Actual behavior

(NOTE: The legend is wrong in this figure, it is Å for the bin width, not nm)
As one can see above, the result for a bin_width of 0.25Å is severely shifted.

Version

0.12

Steps to reproduce

No response

Further information, files, and links

My initial guess as to what is causing this is the following:

I rewrote the DielectricCylinder code some time ago to reduce the effects of larger bins by rewriting the integrals. The calculation for m_r does no longer take bin_width into account. However, the calculation for M_r still does. So larger bins will mean large integration errors.

[hejamu]

Ways to fix this:

• Use a finer spacing to calculate curQ_r_tot internally.
• Find an alternative representation to calculate M_r directly.

[hejamu]

A zoom in on the figure in the original comment shows that even for a bin width of 0.1 Å, which is considered small enough for most applications, the profile still averages out to about -0,25 (!) instead of 1/71....

(NOTE: The legend is wrong in this figure, it is Å for the bin width, not nm)

[hejamu]

I think I figured out the necessary calculations.

We want to compute:

$$M_r = \int_0^R ~ m(r) ~ dr = -\frac{1}{2\pi L} \int_0^R \frac{Q( < r)}{r} ~ dr$$

For charges $q_i$ at (radial) positions $r_i$, the enclosed charge is:

$$Q(< r) = \sum_i q_i , \Theta(r - r_i)$$

Note that this might be the culprit in (at least my) previously failed derivations, since the quantity curQ_r_tot is indeed a sum of Heaviside, not delta functions. (i.e. $$Q(&lt;r)$$ is the total charge enclosed within radius $r$.)

$$M_r = -\frac{1}{2\pi L} \int_0^R \frac{1}{r} \sum_i q_i , \Theta(r - r_i) , dr$$

This step took me a while to realize, but since we have a finite series inside this integral, we can write it out and see that we can take the charge outside of the series of integrals and rewrite them:

$$M_r = -\frac{1}{2\pi L} \sum_i q_i \int_0^R \frac{\Theta(r - r_i)}{r} ~ dr$$

This integral actually has an easy solution:

$$\int_0^R \frac{\Theta(r - r_i)}{r} ~ dr = \int_{r_i}^R \frac{1}{r} ~ dr = \left[\ln r\right]_{r_i}^R = \ln R - \ln r_i = \ln\left(\frac{R}{r_i}\right)$$

Such that we arrive at:

$$M_r = -\frac{1}{2\pi L} \sum_i q_i \ln\left(\frac{R}{r_i}\right)$$

My interpretation of this is that we are weighting every charge by its unique contribution based on its position, charges near the center naturally contribute more the the overall radial dipole (although this integral is not what most people would call a radial dipole...)

[hejamu]

Quick check of 500ns trajectory of the changes in PR #556

Looking good. Should be checked with µs scale data and other bin_widths

[Quewakira]

I compared the results from the integration vs the direct way and they do not match up perfectly. Even for smaller bin 0.001 the result is similar.

bin width
0.01
M_r old integrated
0.0006521030212486064
M_r new dipole like
0.0007952530619446761
mean difference
-0.00014315004069606914
standard error of the difference
1.5896103267912073e-05

I believe the old way (integration) might have been off slightly.

[Quewakira]

I compared the results from the integration vs the direct way and they do not match up perfectly. Even for smaller bin 0.001 the result is similar.

bin width
0.01
M_r old integrated
0.0006521030212486064
M_r new dipole like
0.0007952530619446761
mean difference
-0.00014315004069606914
standard error of the difference
1.5896103267912073e-05

I believe the old way (integration) might have been off slightly.

Nevermind, my system was weird, they both converge