Designing efficient and high-accuracy numerical methods for complex dynamic incompressible Magnetohydrodynamics (MHD) equations remains a challenging problem in various analysis and design tasks. This is mainly due to the nonlinear coupling of the magnetic and velocity fields occurring with convection and Lorentz forces, and multiple physical constraints, which will lead to the limitations of numerical computation. In this paper, we develop the MHDnet as a physics-preserving learning approach to solve MHD problems, where three different mathematical formulations are considered and named $B$ formulation, $A_1$ formulation, and $A_2$ formulation. Then the formulations are embedded into the MHDnet that can preserve the underlying physical properties and divergence-free condition. Moreover, MHDnet is designed by the multi-modes feature merging with multiscale neural network architecture, which can accelerate the convergence of the neural networks (NN) by alleviating the interaction of magnetic fluid coupling across different frequency modes. Furthermore, the pressure fields of three formulations, as the hidden state, can be obtained without extra data and computational cost. Several numerical experiments are presented to demonstrate the performance of the proposed MHDnet compared with different NN architectures and numerical formulations. In future work, we will develop possible applications of inverse problems for coupled equation systems based on the framework proposed in this paper.