A new version of the Delphi program, which provides numerical solutions to the nonlinear Poisson-Boltzmann (PB) equation, is reported. The program can divide space into multiple regions containing different dielectric constants and can treat systems containing mixed salt solutions where the valence and concentration of each ion is different. The electrostatic free energy is calculated by decomposing the various energy terms into Coulombic interactions so that that the calculated free energies are independent of the lattice used to solve the PB equation. This, together with algorithms that optimally position polarization charges on the molecular surface, leads to a significant decrease in the dependence of the electrostatic free energy on the resolution of the lattice used to solve the PB equation and, hence, to a remarkable improvement in the precision of the calculated values. The Gauss-Seidel algorithm used in the current version of DelPhi is retained so that the new program retains many of the optimization features of the old one. The program uses dynamic memory allocation and can easily handle systems requiring large and dimensions-for example a 300(3) system can be conveniently treated on a single SGI R12000 processor. An algorithm that estimates the best relaxation parameter to solve the nonlinear equation for a given system is described, and is implemented in the program at run time. A number of applications of the program are presented.