Accessing the structure of silicate glasses by atomistic simulations remains challenging. Indeed, traditional molecular dynamics simulations depending on the melt-quench method are limited to very high cooling rate and, hence, tend to overestimate the degree of disorder in glasses. Also, conventional reverse Monte Carlo simulations can yield unrealistic structures since a virtually infinite number of configurations can have the same pair distribution function. Here, to overcome these limitations, we report a new method that combines energy minimization with reverse Monte Carlo. Through studying the example of sodium silicate glasses, we show that this method provides some glass configurations that exhibit higher stability than those produced by melt-quench or reverse Monte Carlo. Also, we discuss how this increased degree of stability manifests itself in the atomic structure (with a focus on the medium-range order). This paves the way toward an increased ability to accurately model the structure of silicate glasses.