The Problem
In multislice simulations of transmission electron microscope experiments, it is common to make the approximation that atoms can be thought of lying within a single plane instead of having a finite extent along the $z$ direction. The projected potential of an atom located at $(x_a, y_a)$ will be some 2D radial function $v_z(r)$ where $r = \sqrt{(x - x_a)^2 + (y - y_a)^2}$. Also, these simulations need to sample atomic potentials over a grid of rectangular pixels each with size $\Delta x \times \Delta y$ ($\Delta x$ not necessarily equal to $\Delta y$).
There is a tricky little detail: a common way of parametrizing projected atomic potentials writes $v_z(r)$ as a sum of Gaussian terms like $a \exp(-b r^2)$ and terms like $c K_0(dr)$ ($a,b,c,d$ are constants). Annoyingly, the modified Bessel function $K_0(r)$ goes to infinity as $r \to 0$. However, $\int_0^{r_c} K_0(r) 2\pi r dr$ is finite for any $r_c > 0$, making it still useful for describing sharply peaked atomic potentials.
So this raises a more general question: Consider any 2D radial function $f(r)$ such that
- $\lim_{r \to 0^+} f(r) = +\infty$ and
- $\int_0^{r_c} f(r) 2\pi r dr$ is finite for any $r_c > 0$.
Given a grid of rectangular pixels with locations $\{(x_i, y_j)\}$, how would you sample $f(r)$ over those pixels without serious artifacts? Assume the most general case - that there might or might not be a pixel for which $r = 0$.
Then, if we evaluate the function at each grid point, the pixels nearest to the singularity could be assigned large values that are sensitive to the exact positioning of the grid.
Simple solutions
A common way to mitigate those artifacts is to define some minimum distance $r_{\rm min}$. Then, replace the function $f(r)$ with a better behaved one: $$ f_b(r) = \begin{cases} f(r_{\rm min}), & r < r_{\rm min} \\ f(r), & r \geq r_{\rm min} \end{cases}.$$
In addition to that, I know that the prismatic package also oversamples each pixel.
You can also simply choose a finer sampling, but if you have too much grid points, your simulations will be very computationally expensive.
Yet another way is to tackle this in reciprocal space instead of real space – sample the Fourier transform of $f(r)$ and then take the inverse discrete Fourier transform. At least for $f(r) = K_0(r)$, there will be a nice Fourier transform that does not diverge.
My approach
Back when I was working on this problem in December 2021, I did not know all those things, and I decided to come up with a way in real space to avoid artifacts of sampling $f(r)$. I knew that Michael Cao had a partial solution to this probem in the appendix of his PhD thesis, where he sought to find sampled values $f_s(x_i, y_j)$ such that the integral will be preserved: $$ \sum_{i,j} f_s(x_i, y_j) = \int f(r) 2\pi r dr.$$ So he needed to get $$\int_{\rm pixel} f(r) dA.$$ But he only worked out an approximate analytical solution for the value of a pixel located at $r = 0$. I wanted to finish the job.
I despaired of finding an analytical solution for integrating $f(r)$ over rectangular regions, and turned to numerical integration. (At the time I was making tweaks to Earl Kirkland’s C++ code for multislice simulations, and I just found out about a GSL, a C++ package that can do numerical integration.)
I decided to integrate in polar coordinates, not rectangular ones. This way, I could count on the extra factor of $r$ in the integrand to take care of the singularity at $f(r = 0)$. And then, it was a matter of figuring out how to divide the rectangle into regions to integrate over, being careful to enumerate and handle all possible cases. Here, I will not describe all the cases, but just a couple of illustrative examples.
Examples of integrating over rectangles in polar coordinates
Let the rectangle be bounded by $x_1 < x < x_2$ and $y_1 < y < y_2$.
Define $$ r_{11} = \sqrt{x_1^2 + y_1^2}, $$ $$ r_{12} = \sqrt{x_1^2 + y_2^2}, $$ $$ r_{21} = \sqrt{x_2^2 + y_1^2}, $$ $$ r_{22} = \sqrt{x_2^2 + y_2^2}. $$
Example 1: $x_1 > 0$, $y_1 = 0$, $r_{12} < x_2$
We can break this into three regions:
| Region | Integral |
|---|---|
|
$$+\int_{x_1}^{r_{12}} r \arccos \left( \frac{x_1}{r} \right) f(r) dr$$ |
|
$$+\int_{r_{12}}^{r_{22}} r \arcsin \left( \frac{y_2}{r} \right) f(r) dr$$ |
|
$$-\int_{x_2}^{r_{22}} r \arccos \left( \frac{x_2}{r} \right) f(r) dr$$ |
Example 2: $x_1 > 0$, $y_1 > 0$, $r_{21} < r_{12}$
From the previous example, we see that we are replacing the $2\pi$ in the integral $\int f(r) 2\pi rdr$ with some other angle $\theta(r)$, which is determined by the lines that bound that region. For this example we also break it into three regions.
| Region | Limits | $\theta(r)$ |
|---|---|---|
| 1 | $r_{11}, r_{21}$ | $$\arccos \left(\frac{x_1}{r}\right) - \arcsin \left(\frac{y_1}{r}\right)$$ |
| 2 | $r_{21}, r_{12}$ | $$\arccos \left(\frac{x_1}{r}\right) - \arccos \left(\frac{x_2}{r}\right)$$ |
| 3 | $r_{12}, r_{22}$ | $$\arcsin \left(\frac{y_2}{r}\right) - \arccos \left(\frac{x_2}{r}\right)$$ |
Conclusion
I modified Earl Kirkland’s computem/temsim code to sample atomic potentials in this way. You can find the relevant function here. This way of calculating atomic potentials, of course, takes significantly longer, so I ended up having to store some precomputed values in a table and then bilinearly interpolate from that table.
All this might not be worth the trouble for improving multislice simulations, since this sampling problem is not always that severe and can be mitigated or addressed in other ways. But this gave me an interesting math/programming problem to work on for a month or so.