The PWFs obtained in this way provide an exact rep- resentation of the occupied PBE eigenstates. However, subject to the condition span { φ$_{i}$ } ⊇ span { ψ$_{n}$ } occ as described in Ref. 44. The target NAOs are given by Φ$_{Alm}$( r ) = ζ$_{Al}$ ( r ) Y$_{lm}$ ( r ) where ζ$_{Al}$ is a modified Gaus- sian which vanish outside a specified cut-off radius, and Y$_{lm}$ are the spherical harmonics corresponding to the va- lence of atom A . The number of PWFs equals the num- ber of target NAOs. For example we obtain one PWF for H ( l$_{max}$ = 0), and four PWFs for C ( l$_{max}$ = 1). The PWFs mimick the target atomic orbitals but in addition they allow for an exact representation of all the occupied molecular eigenstates. The latter are obtained from an accurate real-space PAW-PBE calculation $^{45,46}$. In this equation all quantities are matrices in the aug- mented Wannier basis, e.g. H$_{KS}$$_{,ij}$ = 〈 φ$_{i}$ | ˆ H$_{KS}$ | φ$_{j}$ 〉 is the KS Hamiltonian matrix and S$_{ij}$ = 〈 φ$_{i}$ | φ$_{j}$ 〉 is an overlap matrix. The term Δ v$_{H}$ represents the change in the Hartree potential relative to the DFT Hartree potential already contained in H$_{KS}$ , see Appendix A. The central object is the retarded Green function, G $^{r}$, In Ref. 47 the GW-transport scheme was described for the case of an orthogonal basis set and for a truncated, two-index Coulomb interaction. Below we generalize the relevant equations to the case of a non-orthogonal basis and a full four-index Coulomb interaction. Some cen- tral results of many-body perturbation theory in a non- orthogonal basis can be found in Ref. 48. For the GW calculations we apply a basis set consist- ing of projected Wannier functions (PWF) augmented by numerical atomic orbitals (NAO). The PWFs, φ$_{i}$ , are obtained by maximizing their projections onto a set of target NAOs, Φ$_{Alm}$, subject to the condition that they span the set of occupied eigenstates, ψ$_{n}$ . Thus we maxi- mize the functional The HF and GW calculations for isolated molecules are performed using a Green function code developed for quantum transport. 47 In principle, this scheme is de- signed for a molecule connected to two electrodes with different chemical potentials µ$_{L}$ and µ$_{R}$ . However, the case of an isolated molecule can be treated as a special case by setting µ$_{L}$ = µ$_{R}$ = µ and modelling the cou- pling to electrodes by a constant imaginary self-energy, Σ$_{L/R}$ = iη . The chemical potential µ is chosen to lie in the HOMO-LUMO gap of the molecule and the size of η , which provides an artifical broadening of the discrete levels, is reduced until the results have converged. In this limit of small η the result of the GW calculation becomes independent of the precise position of µ inside the gap. The GW and HF calculations presented in Sec. III were performed using a DZP augmented Wannier basis. This gives a total of 5 basis functions per H, Li, and Na, and 13 basis functions for all other chemical elements con- sidered. In Sec. III C we discuss convergence of the GW calculations with respect to the size of the augmented Wannier basis. The paper is organized as follows. In Sec. II we de- scribe the theoretical and numerical details behind the GW calculations, including the augmented Wannier func- tion basis set, the self-consistent solution of the Dyson equation, and the evaluation of valence-core exchange within PAW. In Sec. III we discuss and compare the results of G$_{0}$W$_{0}$, GW, HF, PBE0, and PBE calculations. We analyze the role of dynamical screening, and discuss the effect of self-consistency in the GW self-energy. We conclude in Sec. IV. Here we report on single-shot G$_{0}$W$_{0}$ and fully self- consistent GW calculations of QP energies for a set of 33 molecules. The calculated IPs are compared with ex- perimental values as well as single-particle eigenvalues obtained from Hartree-Fock and DFT-PBE/PBE0 theo- ries. As additional benchmarks we compare to second- order M¨ oller-Plesset (MP2) and DFT-PBE total energy differences between the neutral and cation species. Spe- cial attention is paid to the effect of selfconsistency in the GW self-energy and the role of the initial Green func- tion, G$_{0}$ , used in one-shot G$_{0}$W$_{0}$ calculations. The use of PAW rather than pseudopotentials facilitate the in- clusion of core-valence exchange, which we find can con- tribute significantly to the HF and GW energies. Our results show that the GW approximation yields accurate single-particle excitation energies for small molecules im- proving both hybrid DFT and full Hartree-Fock results. this does not suffice for GW calculations because the po- larizability, P , and the screened interaction, W , do not live in this subspace. Hence we augment the PWFs by additional NAOs including so-called polarization func- tions which have l = l$_{max}$ + 1 and/or extra radial func- tions (zeta functions) for the valence atomic orbitals. For more details on the definition of polarization- and higher zeta functions we refer to Ref. 46. To give an exam- ple, a double-zeta-polarized (DZP) basis consists of the PWFs augmented by one set of NAOs corresponding to l = 0 , ..., l$_{max}$ and one set of polarization orbitals. Note that the notation, SZ, SZP, DZ, DZP, etc., is normally used for pure NAO basis sets, but here we use it to de- note our augmented Wannier basis set. We find that the augmented Wannier basis is significantly better for HF and GW calculations than the corresponding pure NAO basis. assessment of the GW approximation itself. 2 B. GW calculations II. METHOD A. Augmented Wannier function basis Ω = ∑ i ∑ $_{A,l,m}$| 〈 φ$_{i}$ | Φ$_{Alm}$ 〉| 2 (1) G $^{r}$( ε ) = [( ε + iη ) S − H$_{KS}$ + v$_{xc}$ − Δ v$_{H}$ − Σ r $_{xc}$[ G ]( ε )] − 1 (2)