The evaluation of reduction potentials of proteins by ab initio approaches presents a major challenge for computational chemistry. This is addressed in the present investigation by reporting detailed calculations of the reduction potentials of the blue copper proteins plastocyanin and rusticyanin using the QM/MM all-atom frozen density functional theory, FDFT, method. The relevant ab initio free energies are evaluated by using a classical reference potential. This approach appears to provide a general consistent and effective way for reproducing the configurational ensemble needed for consistent ab initio free energy calculations. The FDFT formulation allows us to treat a large part of the protein quantum mechanically by a consistently coupled QM/QM/MM embedding method while still retaining a proper configurational sampling. To establish the importance of proper configurational sampling and the need for a complete representation of the protein+solvent environment, we also consider several classical approaches. These include the semi-macroscopic PDLD/S-LRA method and classical all-atom simulations with and without a polarizable force field. The difference between the reduction potentials of the two blue copper proteins is reproduced in a reasonable way, and its origin is deduced from the different calculations. It is found that the protein permanent dipole tunes down the reduction potential for plastocyanin compared to the active site in regular water solvent, whereas in rusticyanin it is instead tuned up. This electrostatic environment, which is the major effect determining the reduction potential, is a property of the entire protein and solvent system and cannot be ascribed to any particular single interaction.