Large-scale modeling of ultra-thin SiO₂ films on Si(001) surfaces has been performed by means of molecular dynamics utilizing our original inter-atomic potential energy function for Si, O mixed systems. The SiO₂ film is formed by layer-by-layer insertion of O atoms into Si-Si bonds in a Si wafer from the surface. The obtained models reproduce quantitatively the structural transition layers near the interface. Through a modeling of vicinal SiO₂/Si(001) model including atomic steps, it has been found that oxide film near the step-edge is preferentially amorphized. For a more advanced modeling method, we propose a new simulation procedure where O atoms are introduced into the substrate in one-by-one manner. In the calculation, the oxidation is started from the surface and abrupt change in composition at the SiO₂/Si interface is reproduced. Thus, the classical molecular dynamics is a powerful method together with a simplified inter-atomic potential function applicable to mixed systems.