DPPoisson

DPPoisson is an Interactor[1] Module.
Solves the Poisson equation with Gaussian sources at the charges locations in an electroneutral doubly periodic domain. Computes electrostatic force and energy for each particle. Doubly periodic means periodic in XY and unbounded in Z, with the possibility of placing surface charges \sigma_(b/t) at the domain limits and having different permittivities below, inside and above the domain.
A triply periodic solver is also available in UAMMD. See [3] for more info.

Algorithm summary

A complete description of the algorithm can be found in [2] but for completeness a summary is presented here:
This module solves the following equation

Where

With the following boundary conditions:

Where H is the height of the domain.
A similar Ewald splitting technique as employed in the triply periodic solver[3] is used here to separate the algorithm into a near field and far field contribution. This allows to balance the load of the computation between both parts to give an overall complexity that grows linear with the number of charges.
In the doubly periodic case the boundary conditions are taken into account by defining the solution using image charges, reflecting the charges locations across the boundaries (which results in infinite images). A formal description and prove of this method can be found in [2].
Similarly to the triply periodic case the near field interaction kernel decays exponentially, allowing to ignore images beyond the first system reflection.
The far field poses a much more challenging problem since, contrary to the triply periodic case, the problem cannot be solved in Fourier space by simple algebraic multiplication with a Greens function. In our approach the potential is split into two pieces; one with free space boundary conditions and uniform permittivity and another that satisfies the boundary conditions as a correction. Thus

The free space solution is solved by Fourier transforming in the XY directions and interpreting each wave number as a Boundary Value Problem which can be then solved independently in the Chebyshev basis.

The transformation of the charge density in the grid to Chebyshev space is achieved using mirrored 3D FFTs (fast Chebyshev transform) [4].
For the correction potential a Laplace equation for each wavenumber is solved analytically using the jump boundary conditions ensuring that the overall boundary conditions are satisfied in the final potential.

Example

DPPoisson is created as the typical UAMMD module:

#include<uammd.cuh>
#include<Interactor/DoublyPeriodic/DPPoissonSlab.cuh>
using namespace uammd;
...
int main(int argc, char *argv[]){
...
  int N = 1<<14;
  auto sys = make_shared<System>(arc, argv);
  auto pd = make_shared<ParticleData>(N, sys);          
  auto pg = make_shared<ParticleGroup>(pd, sys, "All"); //A group with all the particles
  {
    auto pos = pd->getPos(access::location::cpu, access::mode::write);
    auto charge = pd->getCharge(access::location::cpu, access::mode::write);
   ...
  }
  DPPoissonSlab::Parameters par;
  par.Lxy = {32,32};
  par.H = 32;
  DPPoissonSlab::Permitivity perm;
  perm.inside = 1;
  perm.top = 2;
  perm.bottom = 2;
  par.permitivity = perm;
  par.gw = 1;
  par.Nxy = 72; //Number of Fourier modes for the far field section of the algorithm
  // par.split=1; //Splitting parameter, controls the number of Fourier nodes, choose either this or Nxy directly.
  auto pg = std::make_shared<ParticleGroup>(pd, sys, "All");
  auto dppoisson= std::make_shared<DPPoissonSlab>(pd, pg, sys, par);
...
  myintegrator->addInteractor(dppoisson);
...
return 0;
}

As in [3] the splitting can be tuned via the parameters, in this case either using split or the number of Fourier nodes, Nxy.
The algorithm ensures electroneutrality by placing the opposite of the total system charge as a surface charge in the walls (the same charge to each wall).
The default accuracy parameters, omited in the example (see below), give 3-4 digits of accuracy in the electric field. To ensure this tolerance it is important that charges are kept at least 4*gw away from the walls at all times. For example an steric wall could be placed using ExternalForces[5] to ensure this.

Additional accuracy parameters

These are advanced parameters that control the accuracy of the algorithm. See [2] for additional tested sets of parameters for different accuracies. The default parameters give 3-4 digits of accuracy.

  DPPoissonSlab::Parameters par;
  par.tolerance = 1e-4; //Controls the cut off distance of the near field interaction kernel
  //The far field cell size as: h=sqrt(gw*gw + 1.0/(4.0*split*split))/upsampling;
  //If split is provided upsampling controls the cell size, h. If Nxy is provided upsampling controls the split.
  par.upsampling=1.2; 
  par.support=10; //The number of support cells for the Gaussian charges.
  par.numberStandardDeviations = 4; //Image charges will be considered up to a distance of He=1.25*numberStandardDeviations*sqrt(gw*gw + 1.0/(4.0*split*split));