To realize the numerical calculation of the unified hardening (UH) model for soils, the closest point projection method (CPPM) is used to do the stress integration in the finite element analysis. A package of approaches is proposed to improve the convergence and efficiency of the algorithm. An implicit differentiation method, which can greatly save the computational cost and storage space, is employed to calculate the second derivatives of the 3D yield function generalized by the transformed stress (TS) method. In addition, the denominator of the UH law is equal to 0 at the characteristic state and critical state, while the TS equation contains a square root and may be invalid for tensile stress. These singular problems are solved by introducing subtle treatments without changing the calculation results of the UH model. The consistent tangent operator corresponding to the stress integration can be expressed in a concise form. Finally, a series of simple loading tests and finite element simulations of the bearing capacity of a rigid circular foundation are conducted to evaluate the accuracy, convergence and efficiency of the proposed algorithm.