The temporal evolution of nonlinear wave fields of surface gravity waves is studied by large-scale direct numerical simulations of primitive equations in order to verify Hasselmann's theory for nonlinear energy transfer among component gravity waves. In the simulations, all the nonlinear interactions, including both resonant and non-resonant ones, are taken into account up to the four-wave processes. The initial wave field is constructed by combining more than two million component free waves in such a way that it has the JONSWAP or the Pierson–Moskowitz spectrum. The nonlinear energy transfer is evaluated from the rate of change of the spectrum, and is compared with Hasselmann's theory. It is shown that, in spite of apparently insufficient duration of the simulations such as just a few tens of characteristic periods, the energy transfer obtained by the present method shows satisfactory agreement with Hasselmann's theory, at least in their qualitative features.