  Relativistic corrections to the evolution of structure can be used to test
general relativity on cosmological scales. They are also a well-known
systematic contamination in the search for a primordial non-Gaussian signal. We
present a numerical framework to generate RELativistic second-order Initial
Conditions ($\texttt{RELIC}$) based on a generic (not necessarily separable)
second-order kernel for the density perturbations. In order to keep the time
complexity manageable we introduce a scale cut that separates long and short
scales, and neglect the "short-short" coupling that will eventually be swamped
by uncontrollable higher-order effects. To test our approach, we use the
second-order Einstein-Boltzmann code $\texttt{SONG}$ to provide the numerical
second-order kernel in a $\Lambda$CDM model, and we demonstrate that the
realisations generated by $\texttt{RELIC}$ reproduce the bispectra well
whenever at least one of the scales is a "long" mode. We then present a generic
algorithm that takes a perturbed density field as an inputand provides particle
initial data that matches this input to arbitrary order in perturbations for a
given particle-mesh scheme. We implement this algorithm in the relativistic
N-body code $\texttt{gevolution}$ to demonstrate how our framework can be used
to set precise initial conditions for cosmological simulations of large-scale
structure.
