A variational, self-consistent implementation of the Perdew-Zunger self-interaction correction (PZ-SIC), based on a unified Hamiltonian and complex optimal orbitals, is presented for finite systems and atom-centered basis sets. A simplifying approximation allowing the use of real canonical orbitals is proposed. The algorithm is based on two-step self-consistent field iterations, where the updates of the canonical orbitals and the optimal orbitals are done separately. Calculations of the energy of atoms ranging from H to Ar are presented, using various generalized gradient functionals (PBE, APBE, PBEsol) and a meta-generalized gradient functional (TPSS). While the energy of atoms is poorly described by PBEsol, which is a functional optimized to reproduce properties of solids, the PZ-SIC brings the calculations into good agreement with the best ab initio estimates. The importance of using complex optimal orbitals becomes particularly clear in calculations using the TPSS functional, where the original functional gives good results while the application of PZ-SIC with real orbitals gives highly inaccurate results. With complex optimal orbitals, PZ-SIC slightly improves the accuracy of the TPSS functional. The charge localization problem that plagues Kohn-Sham DFT functionals, including hybrid functionals, is illustrated by calculations on the CH3 + F- complex, where even PBEsol with PZ-SIC is found to give estimates of both energy and charge with accuracy comparable to that of coupled cluster calculations. (Figure Presented).