This paper develops a numerical discretization for solving the Vlasov–Poisson–Fokker–Planck system, which is a model for describing collisional plasma in the presence of a self-consistent electric field. This numerical discretization is established in the particle-in-cell framework in which the dynamical behavior of charged particles follows a stochastic differential equation. We study the geometry structure of the corresponding stochastic differential equation, and construct contact methods to numerically solve it. We also analyze the time accuracy order of the proposed numerical schemes. Numerical experiments are conducted to validate our methods.