This paper is concerned with the numerical integration in time of nonlinear Schrödinger equations using different methods preserving the energy or a discrete analogue of it. The Crank–Nicolson method is a well-known method of order $2$ but is fully implicit and one may prefer a linearly implicit method like the relaxation method introduced in Besse (1998, Analyse numérique des systèmes de Davey-Stewartson. Ph.D. Thesis, Université Bordeaux) for the cubic nonlinear Schrödinger equation. This method is also an energy-preserving method and numerical simulations have shown that its order is $2$. In this paper we give a rigorous proof of the order of this relaxation method and propose a generalized version that allows one to deal with general power law nonlinearites. Numerical simulations for different physical models show the efficiency of these methods.