Although large ground motion databases are today widely available, the selection of records still hampers the use of nonlinear response history analysis in engineering practice. We propose a method that is based on a novel optimization algorithm in order to create subset of records chosen from a large database. Existing heuristic methods select and/or scale ground motion records so that their mean spectrum is close to a target/design design spectrum. The records obtained following this practice offer good estimates of the mean response but considerably underestimate the inherent response variability, thus providing no insight regarding the dispersion around the mean. The proposed method selects and scales the ground motion records so that their mean spectrum and the (period-depended) dispersion fit best a target spectrum and its dispersion. The problem is formulated as a two-objective optimization problem, where the record selection considers both the mean spectrum and its dispersion at the range of periods of interest. The problem is solved with the aid of the Differential Evolution algorithm which searches for potential record combinations whose mean and variance match the target spectral values (median and percentiles) obtained either from a code-compatible smooth spectrum or from a ground motion prediction equation (GMPE). The proposed procedure is efficient, easy to implement and able to quickly search a large pool of ground motion records, identifying record subsets that provide estimates of the response quantities of interest with a minimum number of ground motions. A three-storey moment frame building is considered as a benchmark problem in order to demonstrate the efficiency of the proposed optimization scheme.