The equations governing the motion of N magnetic particles suspended in a viscous fluid at low Reynolds and finite Stokes numbers are solved by direct numerical simulations for different Peclet numbers. The Langevin dynamics simulations include all dipole-dipole magnetic interactions for force and torque. An external applied magnetic field and near field interactions represented by contact and repulsion forces are also considered. Repulsive forces are modeled using a variation of the screened-Coulomb type potential. The initial particle distribution is an ergodic ensemble in which each member consists of N mutually impenetrable spheres whose centers are randomly distributed in a prismatic cell of volume V with wall boundaries. The stability of the proposed numerical method and its convergence in calculating some relevant macroscopic properties of the magnetic suspension are carefully examined. The simulations are used to investigate structure transition from an isotropic random distribution of particles to other structures in the presence of an external magnetic field and magnetic particle-particle interactions. The simulations show dimmers and short chain formation in the suspension space at low volume fraction. When the volume fraction is increased long chains and thin anisotropic structures may be observed along the magnetic field direction. The numerical method is also used to calculate the steady state equilibrium magnetization, and accurate results are obtained for different particle volume fractions phi in agreement with O(phi(3)) asymptotic theories. The method presented is able to consider up to 3000 particles with accuracy and computational efficiency. Typical configurations and particle trajectories are also shown and discussed from a physical point of view. (C) 2015 Elsevier B.V. All rights reserved.