Developers Manual:Open Boundaries/Transport

From OctopusWiki
Jump to navigation Jump to search

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:

Open Boundaries ground state hartree.png

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

(with periodic boundary conditions)

(with a method suitable for finite charge distributions)

and obtain the desired potential as where we have to repeat to fill up the central region.

The calculation of , 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 and for a 2D-1D system (in the terminology of [1]). 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

so that we can write the Hartree convolution integral as

, with being twice the size of the box in -direction.

Calculating this integral amounts to finding the Fourier-transform of which is

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

To: Florian Lorenzen <>
Date: Tue, 22 Jul 2008 12:38:25 +0200 (CEST)
Subject: integral


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

- 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.



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):

Open boudaries partition.png

Calculating partitions of this shape for the parallepiped by hand is straightforward and a partitioning of the (with 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.


  1. C. A. Rozzi, D. Varsano, A. Marini, E. K. U. Gross, and A. Rubio, "Exact Coulomb cutoff technique for supercell calculations", Phys. Rev. B 73, 205119, (2006);