I mean that most of the potentials (Lennard-Jones, DLVO etc..) have a spherical symmetry thus are not appropriate when getting close to a surface/interface. Indeed, the physics of a surface is not the same as the bulk.
It depends in the system and the level of approximation you want. Some potentials that are used in classical md simulations include:
The Tersoff potentials, used for simulations of carbon, silicon, germanium and a range of other materials, which include three atoms interactions.
The (and modified) embedded-atom method and Tight-Binding second moment approximation potentials, which includes the calculation of electron density around an atom from contributions of all surrounding atoms.
Some references:
J. Tersoff, phys. Rev. B, 39, 5566, 1989
Daw et al, Mat. Sci and Engr Rep, 9, 251, 1993
Cleri and Rosato, phys rev B, 48, 1993
Baskes, phys rev B, 46, 2727, 1992
ab-initio md simulations could be another choice, but of course it depends on the system size, since may be too expensive computationally.