In this paper, we provide an optimal rate convergence analysis and error estimate for a structure-preserving numerical scheme for the Poisson-Nernst-Planck-Cahn-Hilliard (PNPCH) system. The numerical scheme is based on the Energetic Variational Approach of the physical model, which is reformulated as a non-constant mobility gradient flow of a free-energy functional that consists of singular logarithmic energy potentials arising from the PNP theory and the Cahn-Hilliard surface diffusion process. The mobility function is explicitly updated, while the logarithmic and the surface diffusion terms are computed implicitly. The primary challenge in the development of theoretical analysis on optimal error estimate has been associated with the nonlinear parabolic coefficients. To overcome this subtle difficulty, an asymptotic expansion of the numerical solution is performed, so that a higher order consistency order can be obtained. The rough error estimate leads to a bound in maximum norm for concentrations, which plays an essential role in the nonlinear analysis. Finally, the refined error estimate is carried out, and the desired convergence estimate is accomplished. Numerical results are presented to demonstrate the convergence order and performance of the numerical scheme in preserving physical properties and capturing ionic steric effects in concentrated electrolytes.