  Electrostatic interactions play crucial roles in biophysical processes such
as protein folding and molecular recognition. Poisson-Boltzmann equation
(PBE)-based models have emerged as widely used in modeling these important
processes. Though great efforts have been put into developing efficient PBE
numerical models, challenges still remain due to the high dimensionality of
typical biomolecular systems. In this study, we implemented and analyzed
commonly used linear PBE solvers for the ever-improving graphics processing
units (GPU) for biomolecular simulations, including both standard and
preconditioned conjugate gradient (CG) solvers with several alternative
preconditioners. Our implementation utilizes standard Nvidia CUDA libraries
cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy
can be achieved given that the single precision is often used for numerical
applications on GPU platforms. The optimal GPU performance was observed with
the Jacobi-preconditioned CG solver, with a significant speedup over standard
CG solver on CPU in our diversified test cases. Our analysis further shows that
different matrix storage formats also considerably affect the efficiency of
different linear PBE solvers on GPU, with the diagonal format best suited for
our standard finite-difference linear systems. Further efficiency may be
possible with matrix-free operations and integrated grid stencil setup
specifically tailored for the banded matrices in PBE-specific linear systems.
