  Numerical treatment of two dimensional strongly-correlated systems is both
extremely challenging and of fundamental importance. Infinite projected
entangled-pair states (PEPS), a class of tensor networks, have demonstrated
cutting-edge performance for ground state calculations, working directly in the
thermodynamic limit. Furthermore, in recent years the application of PEPS has
been extended to also low-lying excited states, using an ansatz that targets
quasiparticle states above the ground state with high accuracy. A major
technical challenge for those simulations is the accurate evaluation of
summations of two- and three-point correlation functions with reasonable
computational cost. In this work, we show how a reformulation of $n$-point
functions in the context of PEPS leads to extra contributions to the results
that prove to play an important role. Benchmarks for the frustrated $J_1-J_2$
Heisenberg model illustrate the improved precision, efficiency and stability of
the simulations compared to previous approaches. Leveraging automatic
differentiation to generate the most tedious and error-prone parts of the
computation, the straightforward implementation presented here is a step
towards broader adoption of the PEPS excitation ansatz in future applications.
