We have developed a 3D magnetohydrodynamics simulation code for applications in the solar convection zone and photosphere. The code includes a non-local and non-grey radiative transfer module and takes into account the effects of partial ionization. Its parallel design is based on domain decomposition, which makes it suited for use on parallel computers with distributed memory architecture. We give a description of the equations and numerical methods and present the results of the simulation of a solar plage region. Starting with a uniform vertical field of 200 G, the processes of flux expulsion and convective field amplification lead to a dichotomy of strong, mainly vertical fields embedded in the granular downflow network and weak, randomly oriented fields filling the hot granular upflows. The strong fields form a magnetic network with thin, sheetlike structures extending along downflow lanes and micropores with diameters of up to 1000 km which form occasionally at vertices where several downflow lanes merge. At the visible surface around optical depth unity, the strong field concentrations are in pressure balance with their weakly magnetized surroundings and reach field strengths of up to 2 kG, strongly exceeding the values corresponding to equipartition with the kinetic energy density of the convective motions. As a result of the channelling of radiation, small flux concentrations stand out as bright features, while the larger micropores appear dark in brightness maps owing to the suppression of the convective energy transport. The overall shape of the magnetic network changes slowly on a timescale much larger than the convective turnover time, while the magnetic flux is constantly redistributed within the network leading to continuous formation and dissolution of flux concentrations.