In Poisson-Boltzmann(PB)model,solute and solvent regions are separated by a dielectric interface and we are interested in tracing its location and conformational change.We numerically solved the PB model under a well-posed boundary integral formulation.The fast treecode algorithms and GPU computing combined make it realistic to apply the PB model to biomolecular dynamics simulation.We derived and implemented the numerical schemes for computing the dominant and computationally demanding electrostatic solvation force.