Monte Carlo simulations in the NpT ensembles have been performed for the structure exploration of aqueous 1,4-dioxane solutions. Three different systems with all-atom dioxane:TIP4P water molar compositions of 2:500 (code:D2), 8:465 (D8), and 17:425 (D17) modeled solutions of 0.22, 0.88, and 1.86 mol/dm(3) concentrations, respectively, at T = 298 K and p = 1 atm. The calculated solution densities increase from 0.992 to 1.002 g/cm(3) with increasing dioxane concentration and approach the experimentally determined densities within 1%. This close agreement was achieved by utilizing RESP charges fitted to the in-solution IEF-PCM/B3LYP/6-31G* electrostatic potential of dioxane taken in its chair conformation and recently developed C, H steric parameters for ethers for calculations with a 12-6-1 all-atom potential. Solution structure analyses pointed out that the dioxane molecules arrange in the solutions with favorable distances of 4-8 angstrom for the ring symmetry centers. Within this range not only pairs of rings but triangular triads and tetrads have also been observed with center-center distances <8 angstrom. For the D8 system, about 25% of the sampled configurations included such a triad. In the case of the D17 model, two simulations starting from different solution configuration predicted different degrees for the dioxane aggregation in aqueous solution. In the more aggregated structure 3-21 triads are consistently maintained and 1-2 tetrads are formed in 58% of the configurations. Each dioxane oxygen forms about one hydrogen bond, on average, to a water molecule in the 0.22-1.86 molar range. The most likely O(dioxane)center dot center dot center dot H(water) hydrogen bond distance is 1.75-1.80 angstrom compared to the optimal distance of 1.72 angstrom in the isolated dimer. The optimal dioxane-water interaction energy of -5.65 kcal/mol indicates a remarkable hydrogen-bond acceptor character for dioxane.