Developers Manual:Open Boundaries/Transport
This site contains a few notes and ideas about the open boundaries implementation.
Contents
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 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:
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 (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.
References
- ↑ 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); http://link.aps.org/abstract/PRB/v73/e205119