Finally, I complement the ndings on the CT model of Man1 and Sec. 4.2.1 by one additional observation on the total energy of this system when static external elds are applied. If one plots as in Fig. A.3 the total energy of the model system versus the eld strength of the external eld, one observes almost straight line segments of the total energy that are interrupted by discontinuous changes of the slope. Slight deviations from these straight lines can be observed in the intermediate segment. Note that GOEP-SIC calculations for eld strengths up to3.5×109 V/m are in almost perfect agreement with the GKLI-SIC ndings.
The eld strengths where the slope changes discontinuously equal the eld strengths where
3The B3LYP polarizabilities were kindly provided by Andreas Karolewski.
A.3. STATIC CHARGE TRANSFER IN A TRANSPARENT MODEL SYSTEM 61
Figure A.3: Dependence of the total energy of the CT model system of the left part of Fig. 4.2 on the strength of the external eld.
I present the results computed with LSDA, GKLI-SIC, and GOEP-SIC. For this illus-tration, I subsumed the degenerate realiza-tions of GKLI-SIC (see Sec. 4.2.1). As a guide for the eye, I tted straight lines with a least-squares linear regression to the three segments of the GKLI-SIC data points from 0.0V/m to2.5×10−9 V/m, from2.5×10−9 V/m to5.5×10−9 V/m, and from6.0×10−9 V/m to 8.0×10−9 V/m. Theses segments correspond to the situation where no electron is transferred, one electron has jumped from the donor to the acceptor chain, and two elec-trons are transferred.
the electrons jump from the donor to the acceptor chain according to Man1 and the right part of Fig. 4.2. Thus, the discontinuous changes of the slope of the total energy can be attributed to these electron jumps. At the corresponding eld strengths, the density reshues notably so that the eld acts on very dierent density distributions before and after the electrons transfer. In Fig. A.3, this manifests in the dierent slopes of the total energy with respect to the electric eld strength. Details of the latter dependence remain to be investigated.
At rst sight, the straight line segments of the total energy of Fig. A.3 appear reminiscent of the straight line behavior of the total energy with fractional changes of the electron number that Perdew et al. reported in Ref. [PPLB82] (for an overview, see Sec. 2.4). Yet, the picture that emerges here is a dierent one as the calculations do not involve fractional charges of the entire system and Fig. A.3 shows the total energy versus the eld strength in contrast to fractional charges that are depicted in Ref. [PPLB82]: Although some values of the external eld strength can be related to specic integer occupations of D and A in GKLI-SIC calculations, none of the intermediate points corresponds to any fractional occupation of the two hydrogen chains. Still, the abrupt changes of the slope at the eld strengths where electrons jump are a manifestation of the integer electron transfer behavior.
In principle, a relation between the external eld strength and the (fractional) charge of the D or of the A chain may be set up in the case of LSDA using the fractional occupations of the single moieties in the LSDA data of the right part of Fig. 4.2. Yet, this picture does not coincide with the fractional occupation picture of the entire system of Perdew et al. [PPLB82] despite the striking analogy, because here the single moieties are not isolated and interactions between D and A need to be taken into account. In Fig. A.3, the LSDA functional exhibits a continuous concave shape of the total energy with respect to the external eld strength. This nding reects the fractional electron transfer behavior of LSDA where no abrupt density reshuing occurs.
Appendix B
Multigrid Poisson solver
The solution of Poisson's problem is one of the key features of the PARSEC code as it is exploited to determine the Hartree potential instead of computing the Hartree potential in-tegral directly [KMT+06]. The standard method to solve Poisson's equation in PARSEC is the conjugate gradient (CG) method [Saa03]. In a standard (semi)local density functional PARSEC run, the solution of Poisson's equation is not of central limiting character to the overall performance of the code as it is required only once per GS self-consistency iteration and once per time step if one assumes that one needs to perform a single potential evaluation per time step. Yet, for orbital-dependent functionals as for example SIC, EXX, and espe-cially energy-minimizing GSIC, the eciency of the Poisson solver is crucial for the overall performance. For instance, in SIC the Poisson solver needs to be called N + 1 times both per self-consistency iteration and per time step, whereN is the number of orbitals involved.
The situation is even more complicated in energy-minimizing GSIC calculations. Here, the Poisson solver is needed at least N times in each step that needs to be taken during the iterative unitary optimization of the unitary transformation in a GS calculation (see Ap-pendix C). Most Poisson calls, however, are performed during TDGSIC time propagation.
For instance, propagation of the 28 occupied orbitals of the DMABN molecule of Sec. 4.2.2 for 40 fs requires5×104 time steps, where in the example per time step and orbital on aver-age4.2Poisson calls are performed. In total this makes about6×106 Poisson calls using the most ecient propagation scheme and algorithms available in PARSEC at current terms.
This example clearly demonstrates that the solution of Poisson's equation is a performance limiting factor if one uses orbital functionals in PARSEC. For the feasibility of many calcu-lations presented in this thesis, I addressed this performance bottleneck and implemented a new method based on multigrid (MG) and defect correction ideas to the solution of Poisson's equation in the PARSEC program package.
In this appendix, I give an elementary introduction to the basic MG and defect correction concepts and explain the implementation of the MG solver in the PARSEC code. The new algorithm turns out to be much more ecient than the CG Poisson solver. Moreover, the parallelization of the MG Poisson solver is adjusted to the requirements of SIC and EXX calculations, where orbital density dependent Poisson problems can be performed in parallel.
For a more detailed introduction to the MG method, I recommend the book of Trottenberg, Oosterlee, and Schüller [TOS01], the book of Hackbusch [Hac85], and the introduction of MG methods in the Numerical Recipes [PTVF92]. Parts of the presentation in the following are based on these books.
63
B.1 Multigrid and defect correction
B.1.1 The multigrid idea
Typically, one uses the MG method for solving second-order partial dierential equations.
Here, for a clear-cut presentation in the PARSEC context, I focus on Poisson's equation in three dimensions
∇2v(r) =−4πn(r) =f(r), (B.1)
wherev(r) is the Hartree potential andn(r)the density. I abbreviate the inhomogeneity on the right hand side of Eq. (B.1) byf(r). For a numerical solution, the Laplacian operator as well as the functionsv(r) andf(r)are discretized on a grid. The discretized equation reads
Lhvh =fh, (B.2)
where the discretized functions carry the index h corresponding to the grid spacing h and Lh denotes the discretized version of the Laplacian operator. If one introduces v˜h to be an approximate solution of Eq. (B.2), the error to obtain vh is
uh =vh−v˜h. (B.3)
The residual between the approximate result of Lhv˜h and the true result of Lhvh reads
rh =Lh˜vh−fh. (B.4)
Based on these denitions, Poisson's equation can be written in terms ofrh and uh as
Lhuh=−rh. (B.5)
One naturally arrives at an iterative procedure for new potentials
˜
v(n+1)h = ˜v(n)h +u(n)h (B.6) that can be obtained from previous potentials ˜v(n)h and solutionsu(n)h of Eq. (B.5), whererh is computed at each step nof the iteration according to Eq. (B.4).
The MG method is used to solve the central Eq. (B.5) of the previous iterative scheme.
Understanding the main idea of the MG method is related to two basic principles [TOS01]:
the smoothing principle and the coarse grid principle. I illustrate both principles using the expansion
uh(x, y, z) =
NXx−1 k=1
NXy−1 l=1
NXz−1 m=1
aklmsin kπ
Nxhx
sin lπ
Nyhy
sin mπ
Nzhz
(B.7) of the error [TOS01], where Nx, Ny, and Nz are the number of grid points in x, y, and z direction, respectively, and uh(x, y, z) vanishes at the boundaries of the grid. The highest frequency components (large k, l, or m) contributing to this series are given by the grid spacing h: Only functions that vary on a length scale that is larger than the grid spacing can be represented after discretization. The maximum extension of the grid determines the lowest frequency components (k,l, andmare small) of the Fourier series. Classical iterative
B.1. MULTIGRID AND DEFECT CORRECTION 65
Figure B.1: The left panel illustrates the MG grid-level management with a one-dimensional example and four grid levels. At each coarser grid level the grid spacing doubles in comparison to the previous ner grid level. The grid points are depicted by full circles . Circles with smaller size are used only on ner grid levels. The right panels show the cycle succession of the three most prominent MG cycling schemes, the V-, W-, and F-cycle. On the very coarsest grid level an exact solution should be computed, whereas on the ner grid levels a couple of iterations of the iterative relaxation method are performed.
relaxation methods, for example Jacobi or Gauss-Seidel relaxation, eciently solve high frequency components of Eq. (B.5). Therefore, application of these techniques amounts to a strong smoothing of the error of any approximation (smoothing principle) [TOS01]. Yet, determining lower frequency components is a tedious task using classical relaxation methods.
This is where the coarse grid principle comes into play [TOS01]. For an explanation, consider two grids: a coarse and a ne one. The highest frequency components that can be sampled on the coarse grid are smaller than the highest frequency components that can be sampled on the ne one. Therefore, only functions that are smooth on the ne grid can be approximated reliably on the coarse grid. If one is sure that the functions that are involved in Poisson's problem are smooth on the ne grid, the problem can be transferred to the coarse grid where nding solutions is more ecient.
The main idea of the multigrid method is to exploit those two principles for a fast solution of Eq. (B.5): High frequency components of the error are relaxed on ne grids and lower frequency components are solved on successively coarser grids. Moreover, this idea can equally be applied starting on the coarse grid side and proceeding to ner grids.
Typically, one arranges the grid levels such that the grid spacing always changes by a factor of two. This type of grid management is illustrated with a one-dimensional example in the left part of Fig. B.1. In this scheme, only a few relaxation steps are needed per grid level.
Subsequently, the solution is transferred either to a coarser grid by a so-called restriction operation or to a ner grid by prolongation. On this new grid level further relaxations are performed. As the classical relaxation methods can be conducted eciently especially on the coarser grid levels, the solution of Eq. (B.5) can be obtained with comparably low numerical eort. For instance, the CG method is reported to need numerical operations on the order of N3/2logεto solve the two-dimensional Poisson problem, while the iterative MG method yields a solution after about Nlogε operations [TOS01]. Here, N denotes the number of unknowns and logε reects the fact that the numerical accuracy and, therefore, the stop criterion are assumed to be in the range of the discretization accuracy.
Based on the central MG idea, dierent cycling schemes have emerged to conduct a MG
cycle. In the following, I introduce the most prominent ones, the so-called V-, W-, and F-cycles [TOS01] as illustrated in the right part of Fig. B.1. All cycling schemes have in common that on the coarsest grid level a high quality solution of Poisson's equation needs to be obtained. In Fig. B.1, this aspect is depicted by the symbol 2. On all other grid levels a small number of standard relaxations are performed. The most basic cycling scheme is the V-cycle. Here, one goes from the nest grid successively down to the coarsest grid and straight back. In the W-cycling scheme, when going from grid level to grid level one always performs one additional V-cycle compared to the V-cycling scheme (see Fig. B.1).
This amounts to a higher number of coarse grid solutions in the W- than in the V-cycle.
One further alternative is the so-called F-cycle. Here, during the back transfer from coarse to ne grids, one consecutively goes nsteps up and down starting withn= 1 and increases n after each round until the nest grid level is reached. The numerical eort related to a full cycle diers [TOS01] for each cycling scheme, as every time a dierent number of steps and operations are involved. Yet, the dierent methods provide means to adapt the MG method to the special requirements of particular problems as one is able to tune the number of relaxations before and after transferring the error between grid levels.
B.1.2 An introduction to defect correction
The numerical load involved in a classical relaxation step strongly depends on the number of neighboring points that are used to compute the Laplacian operator on the grid. One strategy of the PARSEC code to gain numerical eciency is to treat the Laplacian operator of the kinetic energy with a higher-order expansion [KMT+06] as this allows for larger grid spacings. Therefore, also the solution of Poisson's equation needs to be performed by a high-order scheme. Yet, the high numerical performance of MG solvers is obtained by low order discretization using only immediate neighbors in the Laplacian operator as well as the ne to coarse grid restriction and the coarse to ne grid interpolation operations. A solution to this conict is the so-called defect correction [Ste78, Hac81, Sch84, Hac85, TOS01].
Using defect correction one does not solve Poisson's equation
Lˆhvh =fh, (B.8)
with a higher-order discretization Lˆh of the Laplacian operator. Instead, the central idea is to determine the defect
dh =fh−Lˆhvh+Lhvh (B.9) related to dierent discretizations and solve Poissons's equation for this defect
Lhvh =dh (B.10)
with a lower order Laplacian discretization Lh by the MG method. The defect correction idea can be combined with the MG cycling schemes. In this case one needs to perform the following steps:
(i) start from an initial guessv˜h(0), setk= 1
(ii) compute one MG cycle with lower order discretization to obtain ˜v(1)h
(iii) determine the defectd(k)h to the higher-order discretization using the actual approximate solutionv˜h(k) in Eq. (B.9)
B.2. FEATURES OF THE PARSEC CODE IN THE MULTIGRID CONTEXT 67