Frozen density embedding (FDE) has become a popular subsystem density functional theory (DFT) method for systems with weakly overlapping charge densities. The failure of this method for strongly interacting and covalent systems is due to the approximate kinetic energy density functional (KEDF), although the need for approximate KEDFs may be eliminated if each subsystem's Kohn-Sham (KS) orbitals are orthogonal to the other, termed external orthogonality (EO). We present an implementation of EO into the FDE framework within the Amsterdam density functional program package, using the level-shift projection operator method. We generalize this method to remove the need for orbital localization schemes and to include multiple subsystems, and we show that the exact KS-DFT energies and densities may be reproduced through iterative freeze-and-thaw cycles for a number of systems, including a charge delocalized benzene molecule starting from atomic subsystems. Finally, we examine the possibility of a truncated basis for systems with and without charge delocalization, and found that subsystems require a basis that allows them to correctly describe the supermolecular delocalized orbitals.