One of the most important aspects of seismic risk assessment is the characterization of the earthquake hazard through appropriate models that adequately address its variability for different seismicity levels. For applications involving dynamic analysis, especially in the context of performance based earthquake engineering, this requires description of the entire ground motion acceleration time-history. A modeling approach to establish this, gaining increasing interest within the structural engineering community, is the use of stochastic ground motion models. These models are based on modulation of a stochastic sequence through functions that address spectral and temporal characteristics of the ground motion. The parameters of these functions can be related to seismicity (namely moment magnitude and rupture distance) and site (namely soil shear wave velocity) characteristics by appropriate predictive relationships. These predictive relationships provide additional secondary model parameters (such as duration of the excitation and frequency content) that are dependent on the aforementioned seismicity/site characteristics. This work focuses on the selection of these predictive relationships so that compatibility of the resultant ground motions with Ground Motion Prediction Equations (GMPEs) is established. A point-source ground motion model is selected as a stochastic ground motion model and the Next Generation Attenuation (NGA) models are targeted as GMPEs for this study (though the approach is directly extendable to other emerging models). Goal is to offer an app that will provide stochastic ground motions that match specific GMPEs as soon as the end-user defines a set of seismicity characteristics (moment magnitude, rupture distances, local site conditions) to target for this match. The methodology is based on selecting a wide range for the parameters of the predictive relationships, generating ground motions for each of them for different seismicity characteristics and calculating then spectral acceleration for different structural periods. The variability due to the stochastic sequence is addressed by averaging over an ensemble of such sequences. As soon as this database is established and stored, an efficient stochastic search is developed (that also addresses problems with multiple local minima) based on a metamodeling approximations, to select the values for the predictive relationships that optimize the match to GMPEs. This match is defined by comparing the error between the stochastic ground motion predictions and the GMPEs for specific cases, with each case corresponding to a chosen moment magnitude, rupture distance, local site condition and natural period, all of them selected by the end-user. The overall objective function is defined as a weighted least squares error, with the weights chosen based on the importance of having a good match for the respective case. No limitations are placed on the number of cases that can be considered, facilitating a potential match over a great range of seismicity and structural characteristics. Exploiting the computational efficiency of the proposed stochastic search approach a standalone tool is developed to facilitate communication to non-technical audiences, ultimately allowing us to offer a user-friendly applet for generating ground motions compatible with GMPEs.