The selfconsistent solutions of the Dyson equation are obtained using a plane wave basis set for 7 small molecules. Such selfconsistent solutions can help to unify the different GW selfconsistent schemes, reduce the scatter of results in current GW calculations, and shed light on the true effects of GW selfconsistency. Unlike other works of selfconsistent GW calculations, in the present work the Green's function is expressed as a matrix under the plane wave basis set. The algorithmic details

more »

... ch enable such calculations are presented. The ability to solve the full Greens function using a plane wave basis set may open the door for future beyond-GW many-body perturbation theory calculations. GW method has been used as one of the most accurate methods to calculate the electronic structures of materials from bulk crystals to molecules 1 . However, there is a strong dependence of the traditional GW results on the initial input single particle electron wave functions {ψ i } and eigen energies {ǫ i }, especially for G0W0 calculations 2-5 (e.g., up to 1 eV for some oxide band gaps 6 ) where the Green's function G and screened electron-electron interaction W are not updated. One way to solve this problem is to introduce selfconsistency in the solution. However, there are different approaches to solve the GW problem selfconsistently. One can divide these approaches into two major categories. In the first category, the Green's function is still described by an noninteracting Green's function using the eigen functions and eigen energies of the single particle orbital,e.g., G 0 (ω) = i ψ i (r 1 )ψ * i (r 2 )/(ω + µ − ǫ i ± iδ i ), although ψ i and ǫ i will be updated during the selfconsistent iterations 7-11 . However, some of the selfconsistence conditions could be a bit arbitrary since one can propose different selfconsistent schemes 6 . Furthermore, in some cases, there can still be initial wave function dependence 12 . In the second category, for which we will focus on, the Dyson equation is solved selfconsistently, and the Green's function G can no longer be described by a noninteracting single particle Green's function G 0 . Schöne and Eguiluz 13 have solved the Dyson equation by expanding the Green's function G with the input single particle orbital and used a truncation on the number of these orbital. They used the Baym-Kadanoff formalism 14 to express the Dyson equation. Similarly, Kutepov, Savrasov, and Kotliar 15,16 also used the band states to expand the Green's function and iteratively solved the Dyson's equation in a bulk. The Matsubara's frequency mesh 16 on imaginary time and frequency axis is used for frequency integrations. Caruso, Scheffler et al. have used atomic orbital to solve the Dyson equation 17,18 , so is the group of Thygesen 19,20 , they have solved the Dyson equation using localized atomic orbital for 34 different molecules. The atomic orbital has also been used by van Leeuwen et.al 21,22 and Koval et al. 23 to solve the Dyson equation. Computationally, the atomic orbital has the advantage for being able to significantly reduce the dimension of the problem while still be able to cover the important high energy single particle excitations. They are thus particularly suitable for isolated molecule systems, while the periodic crystals are more traditionally solved with plane waves or full potential linearized augmented plane wave (FLAPW) methods. We notice that, for the methods starting with plane wave, or FLAPW basis sets, the band states are often used to expand the Green's function G 13 . However, this could lead to issues related to the truncations of these band states 24,25 . In this work, we will use a plane wave basis set to directly represent the Green's function matrix G without any further truncation. Although doing so will significantly increase the computational cost, as we will show in this work, with the help of modern supercomputers, it is now possible to carry out such calculations. There could be other advantages for adopting this approach. For example, by representing the Green's function matrix in reciprocal and real space, the formalism becomes simpler. This might ease the step to adopt other formalisms beyond GW approximation in the future. One might ask why one should choose the Dyson equation as the selfconsistent solution of the GW problem, given all the possible ways for the GW selfconsistent calculations (i.e., an input equaling output criterion in an iterative procedure). Baym and Kadanoff 14,26 have shown that many conservation laws are preserved under the selfconsistent solution of the Dyson equation. The same is true for the charge conservation law 27 . The Dyson equation is the variational minimum (or stationary point) solution of Klein's total energy expression 28 under the random phase approximation (RPA). This is like the Kohn-Sham equation is the variational minimum solution of the density functional theory (DFT) total energy. Furthermore, it has been shown that 29,30 , under such a variational solution, the differences of the RPA total energies after adding or subtracting one electron equal the GW quasi-particle eigen energies. Recently, there is a surge of using RPA for total energy calculations 29,31-33 . But many such calculations are based on the input (e.g., DFT) noninteracting single particle Green's function G 0 . The selfconsistent solution of the Dyson equation is to find the electronic ground state of the total energy expression. As an result, one can, for example, use the Hellmann-Feynman theory to calculate the atomic force under the RPA total energy. Considering all these factors, it is not difficult to conclude that the Dyson equation as derived from the original GW formalism 34 is the most natural choice for the GW selfconsistent calculations. However, the cost is to represent the Green's function G as a full matrix. It can non longer be represented by a set of single particle wave functions and eigen energies. The effects of self-consistency for GW calculations for homogeneous electron gas has been studied by Holm and Barth 35 . They found an over estimation of the free electron bandwidth and the disappearance of the plasmon satellite structures in the spectral function due to selfconsistent calculations. They thus concluded that the nonselfconsistent G0W0 calculation is preferred unless the vortex correction is included, despite the fact that the selfconsistent GW total energy is found to be rather accurate. Nevertheless, their conclusion is based on model metallic systems. It is thus interesting to test the selfconsistent GW results on real and nonmetallic systems. Ku and Eguiluz 36 have calculated bulk Si and Ge using selfconsistent GW (sc-GW) method, and concluded that the selfconsistency and core level should be used together as their effects can cancel out each other, although this conclusion has been contested by some later studies 24,25 due to the band state truncation issue. In terms of molecules, Caruso et al. 18 found that the selfconsistency does not necessarily make the spectrum result worse than the nonselfconsistent results. Thus, this