An efficient computational method to evaluate the binding energies of many protons in large systems was developed. Proton binding energy is calculated as a corrected nuclear orbital energy using the second-order proton propagator method, which is based on nuclear orbital plus molecular orbital theory. In the present scheme, the divide-and-conquer technique was applied to utilize local molecular orbitals. This use relies on the locality of electronic relaxation after deprotonation and the electron-nucleus correlation. Numerical assessment showed reduction in computational cost without the loss of accuracy. An initial application to model a protein resulted in reasonable binding energies that were in accordance with the electrostatic environment and solvent effects.