Semiclassical molecular dynamics simulations are developed as a tool for studying anharmonic clusters and solids at energies near the zero point. The method employs the time-dependent self-consistent-held approximation, that describes each mode as moving in the mean dynamical field of all other modes. The method further describes each mode by a semiclassical Gaussian wave packet; The scheme is carried out in normal modes. The method is restricted to systems of moderate anharmonicity at low temperatures. It is, however, computationally efficient and practically applicable to large systems. It can be used for the dynamics of nonstationary states as well as for stationary ones. Structural, dynamical and a variety of spectroscopic properties can easily be evaluated. The method is tested for thermal equilibrium states of (Ne)(13), (Ar)(13) against "numerically exact" quantum Feynman path integral simulations. Excellent quantitative agreement is found for the atom-atom pair distribution functions. The method is also applied to (H2O)(n) clusters. Good agreement is found with experimentally available fundamental stretch-mode frequencies.