The Poisson-Boltzmann Equation(PBE)is a widely used implicit solvent model for the electrostatic analysis of solvated biomolecules.To address the expo-nential nonlinearity of the PBE,a pseudo-time approach has been developed in the literature,which completely suppresses the nonlinear instability through an analytic integration in a time splitting framework.This work aims to develop novel Finite Element Methods(FEMs)in this pseudo-time framework for solving the PBE.Two treatments to the singular charge sources are investigated,one directly applies the definition of the delta function in the variational formulation and the other avoids nu-merical approximation of the delta function by using a regularization formulation.To apply the proposed FEMs for both PBE and regularized PBE in real protein systems,a new tetrahedral mesh generator based on the minimal molecular surface definition is developed.Numerical experiments of several benchmark examples and free energy calculations of protein systems are conducted to demonstrate the stability,accuracy,and robustness of the proposed PBE solvers.