We propose a novel bulk phase Monte Carlo simulation technique, in which the energy is calculated by quantum mechanical methods. The semiempirical fragment self-consistent field technique applied divides the periodic simulation cell into two parts. The first is the subsystem where the important change (the random movement of an atom or molecule) takes place and the second is the environment exerting only secondary effects on the former. Expanding the electronic wave function on the basis of strictly localized molecular orbitals and/or atomic hybrid orbitals the wave function of the environment is obtained from simple coupled 2X2 secular equations. The conventional self-consistent field equations, with a perturbation term in the Fockian, have to be solved only for the subsystem. In this way the computational efforts are decreased drastically, as the dependence on the number of atoms in the environment reduces to quadratic instead of cubic or quartic as in conventional semiempirical or ab initio methods, respectively. We wrote a computer code and applied our method to amorphous silicon. Starting from a distorted tetrahedrally bonded random network model we performed Monte Carlo simulations using the fragment self-consistent field energy calculation. After equilibration we obtained distribution functions almost identical to the ones corresponding to the distortion free tetrahedrally bonded network. This finding confirms the adequacy of our method for this specific case.