The first multimessenger observation of a binary neutron star (BNS) merger in August 2017 demonstrated the huge scientific potential of these extraordinary events. This breaktrough led to a number of discoveries, including the proof that BNS mergers can launch short gamma-ray burst (SGRB) jets and are responsible for a copious production of heavy r-process elements. On the other hand, the details of the merger and post-merger dynamics remain only poorly constrained, leaving behind important open questions. Numerical relativity simulations are a powerful tool to unveil the physical processes at work in a BNS merger and as such they offer the best chance to improve our ability to interpret the corresponding gravitational wave (GW) and electromagnetic emission. Here, we review the current theoretical investigation on BNS mergers based on general relativistic magnetohydrodynamics simulations, paying special attention to the magnetic field as a crucial ingredient. First, we discuss the evolution, amplification, and emerging structure of magnetic fields in BNS mergers. Then, we consider their impact on various critical aspects: (i) jet formation and the connection with SGRBs, (ii) matter ejection, r-process nucleosynthesis, and radiocatively-powered kilonova transients, and (iii) post-merger GW emission.