2
$\begingroup$

The 2D Poisson equation is:

(1)$$\frac{d^2\varphi(x,y)}{dx^2}+\frac{d^2\varphi(x,y)}{dy^2}=-\frac{\varrho(x,y)}{\epsilon_0\epsilon}$$

And in $k$-space it is in form of:

(2)$$(k_x^2+k_y^2) \phi(k_x,k_y)=-\frac{\rho(k_x,k_y)}{\epsilon_0\epsilon}.$$

To solve numerically I use complex FFT (FFTW library in C). For area of physical size L and grid of size N (grid const. h=L/N), discrete coordinates and periodic boundaries I have:

(3)$$\rho[k_x,k_y]=\frac{1}{N^2}\sum_{x=0}^{N-1}\sum_{y=0}^{N-1}\varrho[x,y]e^{-j2\pi(xk_x+yk_y)/N}$$

I can multiply both sides by $-\frac{1}{\epsilon\epsilon_0}$, then dividing $\rho[k_x,k_y]$ by $(k_x^2+k_y^2)$ at every point (taking into account that FFT output is symmetric over $k_{x}$ = N/2 and $k_y$ = N/2) I get $\phi[k_x,k_y]$. By inverse FFT:

(4)$$\varphi[x,y]=\sum_{k_x=0}^{N-1}\sum_{k_y=0}^{N-1}\phi[k_x,k_y]e^{j2\pi(xk_x+yk_y)/N}$$ I get 2D potential coming from density of charge $\varrho$. But I am missing some scaling factor in above definitions and I can not figured this out. In short what to do to match the SI units in above definitions so they fit, both in terms of Fourier transform and real result of $\varphi$. My main concern is that if I input the test charge of density $\varrho(x,y)=\frac{e}{h^2}\delta[x,y]$ in units of $C/m^2$ in the middle of area with $L=10^{-6}m$ I expect output similar to Coulomb potential $\frac{1}{r}\frac{e}{4\pi\epsilon\epsilon_0}$, but the solution differs in many orders. Why?

So I think it may be a problem with dividing by k. What should be k for this definition of discrete transform? For example $k=2\pi h/L$ or $k=2\pi x/N$

||cite|improve this question
$\endgroup$
  • $\begingroup$ How did you get the k-space equation? Unless I am terribly mistaken, when you Fourier transform the first equation, the $\epsilon$s are not supposed to drop, and the left-hand side is supposed to get a minus sign. $\endgroup$ – alarge Dec 9 '14 at 13:01
  • $\begingroup$ My mistake I have changed it and added some additional comments. $\endgroup$ – KTKT Dec 10 '14 at 2:06
  • $\begingroup$ I still think the sign of eq. (2) might be wrong. Also, do note that FFTW does not normalize the inverse transform (most FFTs do, but FFTW is really optimized for performance, so this is something that must be manually done). $\endgroup$ – alarge Dec 10 '14 at 2:10
  • $\begingroup$ Yes I normalize manually in (3), about the sign probably you are right but will this change the order of solution? I think only the sign. The sign should be a matter of polarity for positive and negative e charge I assumed on the beginning of calculation. $\endgroup$ – KTKT Dec 10 '14 at 2:33
  • $\begingroup$ I have done figuring out FFTW things once, and damn, it was tedious. Note that, r2c and c2r is used here. gitlab.com/askhl/libvdwxc/blob/master/src/vdw_include.c#L85 ^Lines 85-101 figure out what is k at particular FFT point. Reciprocal cell is calculated here: gitlab.com/askhl/libvdwxc/blob/master/src/vdw_core.c#L182 $\endgroup$ – Mikael Kuisma Nov 24 '15 at 12:55
1
$\begingroup$

Ok I think I solved the problem. So to divide the density FFT by $k^2$ I need actual values of $k$ in $k$-space for my system. FFTW orders the result of transformation in so called "in-order" output, that means in first quadrant of FFT the first pixel corresponds to DC frequency and next to $k/L$ frequency ($k$ is from $0$ to $N-1$) where $L$ is length of whole system. Smallest wave length of the system is $h=L/N$, then $1/h$ is highest frequency. So $k_x$ in above equations should be changed to $2 \pi*k/L$. The coordinate $k$ should also depend on which quadrant of output fft we are, for symmetric quadrant we should use $N-k$ instead of $k$. Probably that is it.

||cite|improve this answer
$\endgroup$
0
$\begingroup$

Substitute $$ \varphi(x,y) = \int dx dy \phi(k_x, k_y) e^{i k_x x + i k_y y},~~\varrho(x,y) = \int dx dy \rho(k_x, k_y) e^{i k_x x + i k_y y} $$ in the first equation and this should immediately give the equation you desire. Also, just for potential future purposes, note $$ \int dx e^{i k x } = 2\pi \delta (k) $$ (See here for discussion.)

||cite|improve this answer
$\endgroup$
  • $\begingroup$ I thought this leads to my equations, isn't it? $\endgroup$ – KTKT Dec 10 '14 at 3:06
  • $\begingroup$ I understood v1 of your question to mean how you obtain the second equation and if you had gotten all the factors of $2\pi$ and what not correct. If that was not what you wanted, then disregard this answer. $\endgroup$ – Prahar Dec 10 '14 at 3:36
  • $\begingroup$ Maybe i was not precise enough in 1st version of the question, I will add the 4pi^2 factor to equation. But my problem is not solved by this. Thank you for answer anyway. $\endgroup$ – KTKT Dec 10 '14 at 3:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.