This work presents a few variational multiscale models for charge transport in complex physical, chemical, and biological systems and engineering devices, such as fuel cells, solar cells, battery cells, nanofluidics, transistors, and ion channels. An essential ingredient of the present models, introduced in an earlier paper [Bull. Math. Biol., 72 (2010), pp. 1562--1622], is the use of the differential geometry theory of surfaces as a natural means to geometrically separate the macroscopic domain from the microscopic domain, while dynamically coupling discrete and continuum descriptions. Our main strategy is to construct the total energy functional of a charge transport system to encompass the polar and nonpolar free energies of solvation and chemical potential related energy. By using the Euler--Lagrange variation, coupled Laplace--Beltrami and Poisson--Nernst--Planck (LB-PNP) equations are derived. The solution of the LB-PNP equations leads to the minimization of the total free energy and explicit profiles of electrostatic potential and densities of charge species. To further reduce the computational complexity, the Boltzmann distribution obtained from the Poisson--Boltzmann (PB) equation is utilized to represent the densities of certain charge species so as to avoid the computationally expensive solution of some Nernst--Planck (NP) equations. Consequently, the coupled Laplace--Beltrami and Poisson--Boltzmann--Nernst--Planck (LB-PBNP) equations are proposed for charge transport in heterogeneous systems. A major emphasis of the present formulation is the consistency between equilibrium Laplace--Beltrami and PB (LB-PB) theory and nonequilibrium LB-PNP theory at equilibrium. Another major emphasis is the capability of the reduced LB-PBNP model to fully recover the prediction of the LB-PNP model at nonequilibrium settings. To account for the fluid impact on the charge transport, we derive coupled Laplace--Beltrami, Poisson--Nernst--Planck, and Navier--Stokes equations from the variational principle for chemo-electro-fluid systems. A number of computational algorithms are developed to implement the proposed new variational multiscale models in an efficient manner. A set of ten protein molecules and a realistic ion channel, Gramicidin A, are employed to confirm the consistency and verify the capability of the algorithms. Extensive numerical experiments are designed to validate the proposed variational multiscale models. A good quantitative agreement between our model prediction and the experimental measurement of current-voltage curves is observed for the Gramicidin A channel transport. This paper also provides a brief review of the field.