The electrostatic osmotic pressure in a system of hexagonally packed DNA molecules has been calculated with the Monte Carlo (MC) simulation method. The DNA molecules were modeled as hard cylinders with charged groups located at the sites corresponding to B-DNA, with the ions considered as point charges with repulsive r(-12) potentials and the solvent treated as a dielectric medium. Periodic boundary conditions for a hexagonal cell were used with Ewald summation of the electrostatic interactions. The pressure was calculated from the relation P = -Delta F/Delta V, where differences of free energies F were obtained with the expanded ensemble method. The calculations were carried but both for salt-free solutions and for solutions containing added salt; in the latter case the simulations were performed within the grand canonical ensemble. In the system with only monovalent ions, the forces between DNA were found to be always repulsive. In the case of divalent :counterions, an effective attraction between DNA molecules may appear for distances of 5-15 Angstrom between the surfaces, depending on the ion size and salt concentration. The results of the simulations showed that a correct statistical-mechanical treatment of the electrostatic interactions in the frame of a continuum dielectric model can reproduce the essential features of available experimental data, indicating that this contribution to the pressure is an important contributor to the experimentally observed pressure versus distance curves for ordered DNA.