The next-to-next-to-next-to-leading order ((NLO)-L-3) Hamiltonian of potential nonrelativistic QCD is derived. The complete matching of the Hamiltonian and the contribution from the ultrasoft dynamical gluons relevant for perturbative bound-state calculations is performed including one-, two-, and three-loop contributions. The threshold expansion is used to disentangle and match contributions of different scales in the effective-theory calculations. As a physical application, the heavy-quarkonium spectrum is obtained at (NLO)-L-3 for the case of vanishing QCD beta function. Our results set the stage for a full (NLO)-L-3 analysis of the heavy-quarkonium system. (C) 2002 Elsevier Science B.V. All rights reserved.