This paper presents a novel method for computing 3D magnetohydrodynamic (MHD) equilibrium states using artificial neural networks. Compared with conventional solvers, we parameterize the Fourier modes as artificial neural networks and minimize the global nonlinear force residual in real space using a first-order optimization technique. As a result, we achieve the same level of residual minimum as that computed by conventional codes at a competitive computational cost, and at the expense of increased computational cost, we obtain a lower residual minimum via neural networks, suggesting a new lower bound on the force residual. With a minimal complexity of neural networks, we expect significant improvements not only in solving single equilibrium states but also in computing neural network models valid for continuous equilibrium distributions.