# Boundary conditions (Poisson's equation)

RESCU is a real space DFT code, which allows it to implement various
types of boundary conditions. RESCU currently supports three types of
boundary conditions: Dirichlet, periodic and Neumann. The chosen
boundary type is applied in solving the Poisson equation. The boundary
type for other quantities, such as the wavefunctions, is inferred (see
`tab:bctype`{.interpreted-text role="numref"}). In this section, we
describe each type of boundary condition, discuss their purpose and
demonstrate their usage.

The domain is always defined in terms of three lattice vectors
(`domain.latvec`). A boundary condition type may be specified along each
direction using the keyword `domain.boundary`. Note, however, that only
periodic directions are allowed to be nonorthogonal. For instance, in a
2D hexagonal heterostructure, `domain.boundary = [1,1,2]`
(periodic,periodic,Neumann) is valid provided that the third lattice
vector is perpendicular to the first two. The default boundary condition
for every direction is periodic (1). In the following sections, we
provide specific details for each boundary condition type.

+-----------+---------+----------------------+----------------------+
| BC        | type    | Potential $v_H$      | Wavefunctions $\psi$ |
+===========+=========+======================+======================+
| Dirichlet | > \(0\) | $v_H(\par            | $\psi(               |
|           |         | tial\Omega) = f_{inp | \partial\Omega) = 0$ |
|           |         | ut}(\partial\Omega)$ |                      |
+-----------+---------+----------------------+----------------------+
| Periodic  | > \(1\) | $v_H(\p              | $\psi(\pa            |
|           |         | artial\Omega_0) = v_ | rtial\Omega_0) = \ps |
|           |         | H(\partial\Omega_1)$ | i(\partial\Omega_1)$ |
+-----------+---------+----------------------+----------------------+
| Neumann   | > \(2\) | $\parti              | $\psi(               |
|           |         | al v_H(\partial\Omeg | \partial\Omega) = 0$ |
|           |         | a)/\partial x_i = 0$ |                      |
+-----------+---------+----------------------+----------------------+
| Open      | > \(3\) | $v_H(\pa             | $\psi(               |
|           |         | rtial\Omega) = \int  | \partial\Omega) = 0$ |
|           |         | dr' \rho(r')/|r-r'|$ |                      |
+-----------+---------+----------------------+----------------------+
| Dirichlet | > \(4\) | $v_H(                | $\psi(               |
|           |         | \partial\Omega) = \i | \partial\Omega) = 0$ |
|           |         | nt dr' \rho(r')/|r-r |                      |
|           |         | '| \\[-2ex] + f_{inp |                      |
|           |         | ut}(\partial\Omega)$ |                      |
+-----------+---------+----------------------+----------------------+

: Boundary condition types. The value of domain.boundary is indicated in
the parenthesis. The second column indicates the boundary condition
imposed on the Hartree potential (obtained by solving the Poisson
equation). The third column indicates the boundary condition imposed on
the wavefunctions. For the periodic boundary condition, the 0 and 1
subscript indicate the lower and upper boundaries along a given
direction.

## Periodic (1)

Periodic boundary conditions are enforced by default on both the Hartree
potential and the wavefunctions. This boundary condition type should be
used in crystals. It may also be used to simulate open boundary
conditions through the use of supercell geometries, albeit this is
suboptimal in terms of efficiency. When using mixed boundary conditions,
e.g. for a slab simulation, the keyword `domain.boundary` must be set to
1 along the periodic directions (as determined in `domain.latvec`).
Consider the following (partial) input file

``` shell
domain.latvec = [100 7.0 0.0
                 7.0 100 0.0
                 0.0 0.0 10]
domain.boundary = [1 1 2]
```

This tells RESCU that the domain is a flat and large monoclinic cell.
The first two dimensions are nonorthogonal, which is allowed since both
are periodic. The third dimension is nonperiodic (Neumann), which again
is allowed since the third lattice vector is perpendicular to the first
two.

## Neumann (2)

The 2D periodic + 1D open geometry is frequently encountered (in
suspended two-dimensional crystals for example). This geometry often
result in the system having a significant global dipole moment, an
artifact of the infinite in-plane dimensions. The Hartree potential then
tends to different values on each side of the system in the
perpendicular dimension. If periodic boundary conditions are enforced,
the Hartree potential is linear in the vacuum region, with a slope that
is inversely proportional to the vacuum size. The convergence of the
Hartree potential is thus very slow when adopting the supercell method.

<figure>
<img src="data/grbn_bc1-2.png" id="fig:grbn_bc1_2" style="width:75.0%"
alt="data/grbn_bc1-2.png" />
<figcaption>Average Hartree potential of a graphene-boron-nitride
heterostructure with respect to the z-coordinate. The Hartree potential
is average in the in-plane directions. The blue curve represents the
Hartree potential as calculated in a periodic supercell. The periodic
boundary conditions force a linear potential in the vacuum region. The
orange curve represents the Hartree potential as calculated with
homogeneous Neumann boundary conditions. The potential then quickly
settles to constant values on either side of the heterostructure. The
yellow curve represents the Hartree potential as calculated with open
boundary conditions.</figcaption>
</figure>

The situation is remedied by adopting homogeneous Neumann boundary
conditions along the perpendicular dimension when solving the Poisson
equation. Simply put, this means that the Hartree potential *derivative*
is set to zero on the domain boundary. Indeed, at a sufficient distance,
we expect the Hartree potential (of a neutral system) to reach a
constant. In this way, the Poisson equation is well defined and the
constant on either side of the structure is allowed to differ, thereby
removing the artificial linear potential between the images. In several
planewave codes, this is referred as *dipole correction*. Apart from
dealing with the global dipole moment, Neumann boundary conditions
usually allow reducing the size of the vacuum used in a simulation,
which may lead to significant computational gains. An example is found
in `fig:grbn_bc1_2`{.interpreted-text role="numref"} where is depicted
the Hartree potential of a graphene-boron-nitride heterostructure as
calculated by periodic and Neumann boundary conditions.

We stress that Neumann boundary condition are only enforced upon the
Hartree potential. The wavefunctions should be confined to the structure
along the nonperiodic dimension, they should decay to zero at the domain
boundary. Mathematically, this means that homogeneous Dirichlet boundary
conditions are imposed upon the wavefunctions (see
`tab:bctype`{.interpreted-text role="numref"}). This in turns result in
homogeneous Dirichlet boundary conditions for the density.

Another point of caution is that the absolute value of the electrostatic
potential is (conventionally) set to zero. This is consistent with the
case of fully periodic boundary conditions. Practically, this means that
the relative vacuum values of the Hartree potential are independent of
the vacuum thickness (beyond a certain vacuum thickness), but their
absolute values vary such as to maintain an average electrostatic
potential of zero.

Finally, as mentioned in the previous section, Neumann boundary
conditions can be imposed along a direction by setting the keyword
`domain.boundary` to 2.

## Induced Dirichlet (0)

Certain simulations demand imposing a specific boundary value to the
Hartree potential. The corresponding type of boundary condition is
called Dirichlet. This can be useful in simulating a gate
nonatomistically, for instance. Note that specifying a value on the
boundary is tantamount to enforcing an induced field, not an external
field. As in the case of Neumann boundary conditions, the wavefunctions
should be confined to the structure along the nonperiodic dimensions,
and hence they satisfy homogeneous Dirichlet boundary conditions.

Dirichlet boundary conditions are selected by setting the keyword
`domain.boundary` to 0. By default, the boundary values are set to zero.
In general, the boundary values along the first, second and third
lattice vectors are specified using the keywords `domain.bvalx`,
`domain.bvaly` and `domain.bvalz` respectively. Without loss of
generality, let's consider `domain.bvalz`. The values on the lower and
upper boundaries are passed via a *cell array*, denote by curly braces.
The value on a boundary can be specified in three ways:

-   If the value on a boundary is constant, it can be specified as a
    simple number;
-   If the value on a boundary is function with an analytical form, it
    can be specified as a function handle which evaluates the function.
    Note that function handle should evaluate a function of two
    variables on the domain $[0,1]\times[0,1]$ in general;
-   If the value on a boundary is numeric-valued function, the function
    can be read from a MAT-file. The path of the MAT-file containing the
    function values should be passed as a string. The function values
    should be arranged in an array that has the exact shape of the
    boundary in question. The exact name of the variable is unimportant,
    but the variable should be lone in the MAT-file. For instance, the
    lower y-boundary of a \[20,30,40\] points domain (`domain.cgridn`)
    should be an array of size \[20,1,40\].

We close this section with the following example

``` shell
domain.latvec = [100 7.0 0.0
                 7.0 100 0.0
                 0.0 0.0 5]
domain.boundary = [1,1,0]
domain.bvalz = {'bvalz0.mat' @(u,v) exp(-((u-0.5).^2 + (v-0.5).^2)/0.05)}
domain.cgridn = [200,200,10]
```

Dirichlet boundary conditions are imposed along the z-direction of this
monoclinic unit cell. The value on the lower boundary is specified by an
array stored in `bvalz0.mat`. Here the array is $200\times 200$, so it
is convenient to define the grid by hand using the keyword
`domain.cgridn`. Adjusting `domain.lowres` to yield a
$200\times 200\times 10$ grid may be difficult. The value on the upper
boundary, of cylindrical Gaussian shape, is specified by a function
handle. A word of caution, the boundary values of a given direction
should respect the boundary conditions along the remaining directions as
possible. Here, the Gaussian function violates the periodic boundary
conditions along the xy-directions. However, its width is sufficiently
small for it to approximately respect periodic boundary conditions.

## Open (3)

If one of the dimensions is open while the other two are periodic, then
the potential on the open boundary can be calculated exactly
`eq:potBC2D`{.interpreted-text role="eq"}.

[$$v_H(\mathbf{r},z) = \sum\limits_{z_i\in \mathcal{H}_i} -2\pi\rho_\mathbf{0}(z_i)|z-z_i| + 2\pi \sum\limits_{\mathbf{k}}^{'}e^{i\mathbf{k}\cdot\mathbf{r}} \frac{\rho_\mathbf{k}(z_i)}{k} e^{-|k||z-z_i|}$$]{label="eq:potBC2D"}

If two of the dimensions are open while the other is periodic, then the
potential on the open boundary can be calculated using
`eq:potBC1D`{.interpreted-text role="eq"}.

[$$v_H(x,y,z) = \sum\limits_{x_i,y_i\in \mathcal{H}_i} -2\log(d_{i}(x,y)) \rho_\mathbf{0}(x_i,y_i) + 2 \sum\limits_{k}^{'}e^{ikz} K_0(d_{i}(x,y)k) \rho_k(x_i,y_i)$$]{label="eq:potBC1D"}

where $d_{i}(x,y) = \sqrt{(x-x_i)^2+(y-y_i)^2}$ and $K_0$ is the
decaying modified Bessel function. After calculating the Hartree
potential on the boundary, the equation can be solved using it as a
regular Dirichlet boundary condition. This is done by setting the value
of the keyword `domain.boundary` equal to 3 for the corresponding
dimension. As can be seen in `fig:grbn_bc1_2`{.interpreted-text
role="numref"}, the transverse-averaged Hartree potential calculated
using open boundary conditions matches the one calculated with Neumann
boundary conditions, albeit with a global shift. This is because the
Neumann absolute level is chosen such that the average potential is zero
while it is calculated from `eq:potBC2D`{.interpreted-text role="eq"} in
the former case. Note that unlike homogeneous Neumann boundary
conditions, open boundary conditions can capture all multipole moments
of a slab sytems and not only the perpendicular dipole moment.

## External Dirichlet (4)

Often, the external field should be specified instead of the induced
field. In this case, `domain.boundary` should be set to 4 instead of 0.
This calculation requires calculating the potential on the open
boundaries, and hence the same restrictions as open boundary
calculations apply (one of the dimensions is open while the other two
are periodic). The external field can be specified using `domain.bvalx`,
`domain.bvaly` or `domain.bvalz` just like the type-0 Dirichlet boundary
conditions.
