Developers Manual:Open Boundaries/Transport

This site contains a few notes and ideas about the open boundaries implementation.

Calculating the Hartree potential for the ground state

In order to do DFT for the transport geometry we have to solve the Poisson equation for a system with broken symmetry. The idea is illustrated by the following figure:

In order to solve the Poisson equation for a system with potential $v'$ and density $n'$ at its left and right and density $n$ in the middle we define an auxiliary system with localized density $\rho$ in the middle and zero density outside. We choose $\rho$ such that $n=\rho+n'$. We then solve the equations

$\nabla^2 v' = -4\pi n'$ (with periodic boundary conditions)

$\nabla^2 u = -4\pi \rho$ (with a method suitable for finite charge distributions)

and obtain the desired potential $v$ as $v=u+v'$ where we have to repeat $v'$ to fill up the central region.

The calculation of $v'$, $n'$ is taken care of by the periodic mode of Octopus. The second equation can be solved by direct summation (in 1D and 2D), FFT with cutoff, ISF, or multigrid (in 3D).

A note for 2D systems

We need $v'$ and $n'$ for a 2D-1D system (in the terminology of ). This is currently not implemented in Octopus but it should be possible to calculate the potential this way:

We define the auxiliary periodic periodic interaction

$\varphi(x, y) = \begin{cases} \frac{1}{\sqrt{x^2+y^2}}, & |y| \le R \\ 0, & |y| > R \end{cases}$

so that we can write the Hartree convolution integral as

$v'(x, y) = \int dx'\int dy' n'(x', y') \varphi(x-x', y-y')$, with $R$ being twice the size of the box in $y$-direction.

Calculating this integral amounts to finding the Fourier-transform of $\varphi(x, y)$ which is

$I(G_x, G_y) = \int_{-\infty}^\infty dx \int_{-R}^{R} dy \frac{1}{x^2+y^2}e^{-i(G_x x + G_y y)}.$

Solving this integral is a bit hard, but Alberto suggested the following way to do it:

From: alberto@physik.fu-berlin.de
To: Florian Lorenzen <lorenzen@physik.fu-berlin.de>
Date: Tue, 22 Jul 2008 12:38:25 +0200 (CEST)
Subject: integral


Hi,

I was trying to figure out the integral, and I came up with the following
solution -- which I cannot integrate further. But they are one-dimensional
integrals, which should be fast:

- Assuming Gx .ne. 0:

  I(Gx,Gy) = 4 \int_0^R cos(Gy*y)*K_0(|Gx|y)

    where K_0 is the modified Bessel function of the second kind, order 0.
    I think that it is the correct solution because, if one takes R ->
    infinity, one retrieves

    I(Gx,Gy) = (2*Pi)/|G|

    which is the expected result, that one obtains easily using polar
    coordinates.

- If Gx .eq. 0:

  I(0, Gy) = -4\int_0^R cos(Gy*y)*log(y)

  In particular, if both Gs are zero:

  I(0, 0) = R(-1+log(R))


Note that, strictly speaking, I(0,Gy) is infinite. But these infinites
should cancel out because of charge neutrality, as explained in the
article of Rozzi.

The 1D integrals that have to be done are singular at zero, which makes
them a bit more tricky, but I think that maybe using the GSL integration
capabilities, they should cause no problem.

Cheers,

Alberto.

Parallelizing the calculation of the memory coefficient and the memory term

One of the bottlenecks of the open boundaries propagator is the calculation of the memory coefficients and the calculation of the memory term which are both of quadratic complexity in the number of timesteps.

It should be possible to parallelize their calculation by domain parallelization with a certain partitioning and the usage of PBLAS instead of BLAS for all the matrix-matrix and matrix-vector multiplication.

The idea is to partition the hrid in such a way that each node gets a share of the interface region, i. e. the part of the wave functions that are affected by the memory term. This is easily obtained by a partitioning in stripes (instead of calling METIS):

Calculating partitions of this shape for the parallepiped by hand is straightforward and a partitioning of the $N\times N$ (with $N$ being the number of interface points) follows naturally which should be compatible with PBLAS requirements.

Of course, this partitioning is not optimal with respect to ghost point communication but should pay off in the end because the dominant terms are calculated faster.