where r i indicates the Cartesian coordinates of the ith restrained atom, whose anchor position is r 0 i . The same force constant is used for each restrained atom, and its value (k = 450 N/m) reflects a typical C-C single bond. In the present work, all restrained atoms are C α carbons of the peptide backbone. Analytic derivatives of V conf are incorporated into the gradient and Hessian calculations, and no gradient or Hessian elements are set to zero or otherwise manipulated within this approach. Electronic structure calculations were performed at the B3LYP+D3(BJ)/6-31G(d,p) level, where D3 indicates Grimme's empirical dispersion correction, 64 evaluated using the Becke-Johnson (BJ) damping function. The SG-1 quadrature grid is used for the exchangecorrelation integration. The self-consistent field (SCF) convergence criterion was set to 10 -8 Ha, with integral and shell-pair drop tolerances set to 10 -11 a.u.. All calculations were performed using the Q-Chem program, in which we have implemented both harmonic confining potentials and the zeroed-out Hessian approach. (Both methods are available starting from Q-Chem v. 6.0.) For some calculations, the conductor-like polarizable continuum model (C-PCM) was used to incorporate dielectric boundary conditions. For these calculations, the solute cavity was constructed from atomic spheres with radii 1.2 times larger than those in the modified Bondi set. This cavity was discretized using the switching/Gaussian algorithm, with 110 Lebedev grid points for each hydrogen and 194 Lebedev points for other atomic spheres. Note that low-dielectric boundary conditions (with ε = 2-4) can be critically important for converging SCF on large protein models with ionic side chains.