As the semiconductor technology node continues to shrink, achieving smaller critical dimension in lithography becomes increasingly challenging. Negative tone development (NTD) process is widely employed in advanced node due to their large process window. However, the unique characteristics of NTD, such as shrinkage effect, make the NTD resist model calibration more complex. Gradient descent (GD) and heuristic methods have been applied for calibration of NTD resist model. Nevertheless, these methods depend on initial parameter selection and tend to fall into local optima, resulting in poor accuracy of the NTD model and massive computational time. In this paper, we propose cluster sampling and scalable Bayesian optimization (BO) with constraints method for NTD resist model calibration. This approach utilizes cluster sampling strategy to enhance the capability for global initial sampling and employs scalable BO with constraints for global optimization of high-dimensional parameter space. With this approach, the calibration accuracy is significantly enhanced in comparison with results from GD and heuristic methods, and the computational efficiency is substantially improved compared with GD. By gearing up cluster sampling strategy and scalable BO with constraints, this method offers a new and efficient resist model calibration.