follows $$m(x,\delta) =
  \left\{ \begin{array}{ll}
          \displaystyle\frac{\delta^2}{2\delta-x} ,&\mbox{ if } x < \delta\\
          x,&\mbox{ otherwise},
            \end{array}\right.$$ If $\delta \to 0$ we recover the classical Macaulay brackets: $m(x,\delta) \to \langle x\rangle$. This approximation creates small repulsive forces between the sphere and asperities, but they decay quite rapidly (as $z^{-3/2}$), the parameter $\delta$ reduces to zero during the convergence loop enabling to ensure a better convergence for Newton-Raphson method, initial value of $\delta$ is chosen to be $\sigma/100$. With this smoothing the system of equations  reduces to $${\boldsymbol{F}} - {\boldsymbol{C}} m( {\boldsymbol{G}}_0 - {\boldsymbol{B}} {\boldsymbol{F}},\delta)^{3/2} = 0,
 \label{eq:sys_approx}$$ whose solution tends to the solution of  when $\delta\to0$. The initial guess for the force vector ${\boldsymbol{F}}_0$ is computed from a model without interaction (${\boldsymbol{B}} = 0$), then from Eq. : $${\boldsymbol{F}}_0 = {\boldsymbol{C}} \langle {\boldsymbol{G}}_0 \rangle^{3/2}.$$ A classical choice to take the initial guess from previous converged load step is not used since the possible set of active asperities changes from step to step (see next Section) and even if not, such initial guess often leads to divergence of the computational step. From the initial guess, we obtain the increment of $\Delta {\boldsymbol{F}}$ as: $$\Delta{\boldsymbol{F}} = -{\boldsymbol{K}}^{-1}{\boldsymbol{R}},$$ where the free term ${\boldsymbol{R}}$ and the 