This paper presents a parallel finite element algorithm based on the toolbox PHG for solving the nonlinear Poisson-Boltzmann equations for biomolecular systems. Previously TMSmesh, TransforMesh and ISO2Mesh were used to generate high resolution meshes. In this work a new method which combines mesh generation and adaptive mesh refinement is introduced to avoid complicated mesh generation steps. The new approach makes use of the capabilities of PHG in local mesh refinement and dynamic load balancing. The efficiency of the new method is demonstrated with a sphere model and the AChE system through both the decay of the a posteriori error estimates and the convergence of solvation-energy. In the meantime the method is also successfully applied to the Gramicidin Aion channel system.