The electrostatic potential (EP) generated by the protein α-lactoalbumin in the presence of added salt is computed as a thermal average at a given point in space. With this aim, constant pH Monte Carlo (MC) simulations are performed within the primitive model, namely, the solvent is treated as a continuum dielectric. The study of the thermal and spatial fluctuations of the EP reveals that they are remarkably high inside the protein. The calculations indicate that fluctuations inside the protein are mainly due to the asymmetric distribution of the charge groups, while the charge fluctuations of the titratable groups play a minor role. The computed EP matches very well with the one obtained from the Poisson equation for the average charge density in spherical symmetry. The Tanford-Kirkwood multipole expansion reproduces the simulated angular-averaged potential rather accurately. Surprisingly, two of the simplest mean-field models, the linear Poisson-Boltzmann (PB) equation and Donnan potential, provide good estimations of the average EP in the effective protein surface (surface EP). The linear PB equation predicts a linear relationship between charge and surface EP, which is numerically reproduced only if the small ions within the protein are taken into account. On the other hand, the partition coefficients of the small ions inside and outside the protein predicted by Donnan theory reproduce reasonably well the simulation results.