including structural properties (e.g., polar nanoregions or PNRs) using neutron scattering [@PhononLocalization], dielectric responses [@Bokov2006a; @Nuzhnyy2012; @Petzelt2014; @Kleeman2014; @Wang2014; @Wang2016], the crossover from ferroelectrics to relaxor [@Kleeman2014]. To understand such systems, many theories have been proposed [@Pirc1999; @Bokov2012; @Akbarzadeh2012; @Sherrington2013a; @Uchino2014; @Sherrington2014; @PhononLocalization; @Kleemann2015]. Relaxors are complex systems, to some extent similar to spin glasses [@Sherrington2013a; @Sherrington2014a; @Sherrington2014], in that their compositions are, without exception, made of complex oxides containing different ions and inevitably inhomogeneous. For instance, the B-site ions of typical relaxors Ba(Zr,Ti)O$_{3}$ (BZT) and Pb(Mg$_{1/3}$Nb$_{2/3}$)O$_{3}$(PMN) are randomly distributed.

The dielectric response of ferroelectric relaxors is the defining feature that differentiates them from normal ferroelectrics: (i) large susceptibilities at low frequency (GHz or lower) ; (ii) even more unusual, the characteristic temperature $T_{m}$, at which the susceptibility peaks strongly depends on the frequency of the probing *ac* electric field. In other words, the susceptibility, $\chi$, depends on both temperature, $T$, and the probing frequency, $\nu$. While such phenomena are well known experimentally [@Bokov2007BZT; @Bokov2012; @Gridnev2004; @Bokovo2006; @BZT-susceptibility; @Tagantsev1999; @PMN-BZT], numerical generation of relaxor’s dielectric response from model-based simulations has been a challenging work. For instance, the shift of $T_{m}$ of the lead free relaxor BZT was only numerically achieved 