Since 2015, when the LIGO-Virgo collaboration announced the first simultaneous detections of gravitational waves (GW150914) until now, more than 66 gravitational wave detections were reported, but only two signals came from a binary neutron star collision, the GW170817 and GW190425 events. GW170817 was accompanied by an electromagnetic outburst manifested as a kilonova and an off-axis jet. However, no conclusive electromagnetic signature was found to come from the GW190425. Indeed, nature proves again more complicated than our models, and it is still a big question how to model kilonovae, or to understand the mechanisms driving astrophysical jets and gamma ray bursts. Not only the electromagnetic counterpart of such collisions, but also the nature of the remnant is also mostly unknown, and its investigation can give glimpses of the internal structure of neutron stars. This talk will focus on the challenges encountered endeavoring to perform numerical relativity simulations of realistic binary neutron star collisions, using open-source General-Relativistic Magnetohydrodynamic codes such as GRHydro, Spritz and IllinoisGRMHD, that run within the Einstein Toolkit software platform. We construct initial data for 6 parametrized hybrid equations of state, and add both interior magnetic field, and exterior magnetosphere, then perform simulations of BNS mergers within the mass range of the two detections. Besides addressing the problem of the reproducibility of the results, the aim of this talk is to forge a path through this complicated procedure, and make it approachable to graduate students, or researchers that just enter the field.