A new approximate method for calculating the thermodynamic functions of multispin systems is considered for the example of the simplest Ising lattices. The idea of the method is based on the hypothesis that the Helmholtz energy of a system is obtained by averaging the Helmholtz energies of clusters over the lattice. The partition function of a cluster is found through usual spur calculations over the nodes constituting its nucleus followed by weighted averaging over the nearest-neighbor nodes. Some random correlated functions are used as weight functions. It is assumed that precisely these functions contain correlations between spins, that is, take into account some long-range order elements. Two clusters (mono- and binuclear) are considered, and two approaches are formulated. In the first approach, unknown parameters are determined using the variational principle. In the second one, temporal spin correlation functions are introduced, and the equations for parameter calculations are found from the boundary conditions for the temporal correlation functions. Both approaches are already equivalent to the Bete-Peierls approximation in the zeroth approximation.