The Heisenberg {\it ab initio} theory of magnetization is developed to apply for multilayer nanoparticles. The theory is based on distribution and partition functions modification with account the difference between exchange integral and closest neighbour numbers, that change the system of resulting transcendental equation for magnetization and its reversal to form either a paramagnetic type curve or hysteresis loops patterns. The equations are obtained within the Heisenberg partition function construction by Heitler diagonalization of energy matrix via irreducible representations of permutation symmetry group. A combination with the Gauss distribution gives the explicit expression for the partition function in the asymptotic limit] at large spin range in terms of transcendent function. The exchange integral, as a parameter of the equation of state (material equation) is evaluated from Curie temperature value by means of a formula derived within the presented theory. Methods of data processing from the simultaneous solution of the material equation system are proposed. The multi-valued function of hysteresis loop is found by combination of graphical approach and special procedure for elimination of mistaken peaks and prolapses of the patterns. The theory and computation methods are applied to spherical particles with separate surface layers consideration. The contribution of the surface layers, that are specified by number of closest neighbors and exchange integrals into overall magnetization, is studied for two-layer and three-layer models, that are discussed and compared graphically.