We report the results of Molecular Dynamics (MD) simulations and formal modeling of the free energy surfaces and reaction rates of primary charge separation in the reaction center of Rhodobacter sphaeroides. Two simulation protocols were used to produce MD trajectories. Standard force field potentials were employed in the first protocol. In the second protocol, the special pair was made polarizable to reproduce a high polarizability of its photoexcited state observed by Stark spectroscopy. The charge distribution between covalent and charge-transfer states of the special pair was dynamically adjusted during the simulation run. We found from both protocols that the breadth of electrostatic fluctuations of the protein/water environment far exceeds previous estimates resulting in about 1.6 eV reorganization energy of electron transfer in the first protocol and 2.5 eV in the second protocol. Most of these electrostatic fluctuations become dynamically frozen on the time-scale of primary charge separation resulting in much smaller solvation contributions to the activation barrier. While water dominates solvation thermodynamics on long observation times, protein emerges as the major thermal bath coupled to electron transfer on the picosecond time of the reaction. Marcus parabolas were obtained for the free energy surfaces of electron transfer by using the first protocol while a highly asymmetric surface was obtained in the second protocol. A non- ergodic formulation of the diffusion-reaction electron transfer kinetics has allowed us to reproduce the experimental results for both the temperature dependence of the rate and the non-exponential decay of the population of the photoexcited special pair.