where f I is the Lagrange multiplier associated with the expectation position constraint for the Ith quantum nucleus. Making the Lagrangian stationary with respect to electronic and quantum nuclear orbital variations will yield a new set of coupled multicomponent Kohn-Sham equations that can be solved self-consistently together with the nuclear expectation constraints. The energy of the final quantum states is a function of both quantum nuclear expectation positions and classical nuclear positions. Thus, it can be evaluated at any molecular geometry, resulting in the CNEO energy surface that naturally incorporates NQEs of the quantum nuclei. It can be proved that the Lagrange multiplier f I acts as the classical force on the Ith quantum nucleus in the complete basis limit, although in practical calculations with a finite basis set, the Pulay force for quantum nuclei will be included in gradient calculations. The CNEO energy surface can be employed for classical molecular simulations. The underlying theoretical justification is the constrained minimized energy surface molecular dynamics (CMES-MD), a MD framework that our group developed to incorporate nuclear quantum effects in classical simulations. The equations of motion for CMES-MD are