The extension of molecular mechanics to reactive systems, metals, and covalently bonded clusters with variable coordination numbers requires new functional forms beyond those popular for organic chemistry and biomolecules. Here we present a new scheme for reactive molecular mechanics, which is denoted as the valence-bond order model, for approximating reactive potential energy surfaces in large molecules, clusters, nanoparticles, solids, and other condensed-phase materials, especially those containing metals. The model is motivated by a moment approximation to tight binding molecular orbital theory, and we test how well one can approximate potential energy surfaces with a very simple functional form involving only interatomic distances with no explicit dependence on bond angles or dihedral angles. For large systems the computational requirements scale linearly with system size, and no diagonalizations or iterations are required; thus the method is well suited to large-scale simulations. The method is illustrated here by developing a force field for particles and solids composed of aluminum and hydrogen. The parameters were optimized against both interaction energies and relative interaction energies. The method performs well for pure aluminum clusters, nanoparticles, and bulk lattices and reasonably well for pure hydrogen clusters; the mean unsigned error per atom for the aluminum-hydrogen clusters is 0.1 eV/atom.