6
$\begingroup$

I know that this question might be very specific and maybe nobody will know the answer, but this is probably the only community where I could find an answer:

So, as part of my master's project, I am currently writing a multigrid code for Poisson's equation in order to calculate particle interactions in a molecular dynamcis setup. I'd like to use periodic BCs, but here arises one problem:

According to the Multigrid book by Trottenberg, periodic BCs only have a solution if $\sum_{x, y, z} f_h(x,y,z) = 0$, where $f_h$ is the right hand side and the sum is over all gridpoints.

The algorithm I found in some publications samples the charged particles via a density function to the grid. So in general, even if the total charge of the system is neutral, this condition will not be fulfilled.

Unfortunately, all papers I found only deal with this PBC problem on a quick sidenote and I wasn't able to get any valuable information out of it.

The code I wrote so far produced fine results when using different functions which fulfill this condition, but fails when going to the real data, i.e., sampled particle charge densities. My question now: Does anybody know this problem and the solution, or any idea which could help me?

$\endgroup$
  • $\begingroup$ Hi @chris, and welcome to Scicomp! I'm curious... Are you using purely periodic boundary conditions? $\endgroup$ – Paul Mar 3 '13 at 19:46
  • $\begingroup$ Hey and thank you. Yes I am. It is a conventional way when using FFT based methods such as SPME or P3M, so I also want to have periodic BCs to be able to compare the results $\endgroup$ – chris Mar 3 '13 at 20:27
  • $\begingroup$ What kind of function are you using to map your particles to the grid? Do they preserve the sum of the quantity you are mapping? $\endgroup$ – Pedro Mar 3 '13 at 21:07
6
$\begingroup$

This isn't a multigrid issue, it's a problem formulation issue. Consider the symmetric system $A x = b$ and suppose that $A e = 0$ for some nonzero vector $e$ (the constant when $A$ is the Laplacian with periodic or Neumann boundary conditions). This system is singular and has no solution if you choose $b$ such that $e^T b \ne 0$. Indeed, the exact periodic problem is ill-posed because the integral of the Coulomb potential with some periodic change density (not averaging to zero) over all space is infinite. That is, if $\sigma(r')$ is periodic with nonzero integral, then

$$\Phi (r) = \int_{r'\in R^3} \frac{\sigma(r')}{4 \pi \vert r - r'\rvert} \mathrm{d} r' = \pm \infty .$$

Fortunately, all that matters is the field, not the potential, so you can renormalize. This amounts to projecting out the null space by solving

$$ A x = (I - e e^T) b $$

where $e$ is the normalized null space vector.

$\endgroup$
  • 2
    $\begingroup$ Thank you, this works. Just for my personal understanding: Why can we introduce this shift of $b$ without affecting the field? $\endgroup$ – chris Mar 3 '13 at 23:13
  • 2
    $\begingroup$ The field $F = \nabla \Phi = \nabla (\Phi + C)$ for any constant $C$. $\endgroup$ – Jed Brown Mar 3 '13 at 23:19

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.