In regression applications, the presence of nonlinearity and correlation among observations offer computational challenges not only in traditional settings such as least squares regression, but also (and especially) when the objective function is nonsmooth as in the case of quantile regression. Methods are developed for the modelling and estimation of nonlinear conditional quantile functions when data are clustered within two-level nested designs. The proposed estimation algorithm is a blend of a smoothing algorithm for quantile regression and a second order Laplacian approximation for nonlinear mixed models. This optimization approach has the appealing advantage of reducing the original nonsmooth problem to an approximated L2 problem. While the estimation algorithm is iterative, the objective function to be optimized has a simple analytic form. The proposed methods are assessed through a simulation study and two applications, one in pharmacokinetics and one related to growth curve modelling in agriculture.