We develop scalable parallel domain decomposition algorithms for nonlinear complementarity problems including, for example, obstacle problems and free boundary value problems. Semismooth Newton is a popular approach for such problems, however, the method is not suitable for large scale calculations because the number of Newton iterations is not scalable with respect to the grid size; i.e., when the grid is refined, the number of Newton iterations often increases drastically. In this paper, we introduce a family of NewtonKrylov-Schwarz methods based on a smoothed grid sequencing method, a semismooth inexact Newton method, and a two-grid restricted overlapping Schwarz preconditioner. We show numerically that such an approach is totally scalable in the sense that the number of Newton iterations and the number of linear iterations are both nearly independent of the grid size and the number of processors. In addition, the method is not sensitive to the sharp discontinuity often associated with obstacle problems. We present numerical results for several large scale calculations obtained on machines with hundreds of processors.