  We present a calculation of the matter power spectrum covariance matrix ${\rm
Cov}(\bf{k}_1,\bf{k}_2)$ that uses power spectrum responses to accurately
describe the coupling between large- and small-scale modes beyond the
perturbative regime. These response functions can be measured with
(small-volume) N-body simulations, which is why the response contributions to
the covariance remain valid and predictive at all orders in perturbation
theory. A novel and key step presented here is the use of responses to compute
loop contributions with soft loop momenta, which extends the application of the
response approach beyond that of previously considered squeezed $n$-point
functions. The calculation presented here does not involve any fitting
parameters. When including response-type terms up to 1-loop order in
perturbation theory, we find that our calculation captures the bulk of the
total covariance as estimated from simulations up to values of $k_1, k_2 \sim
1\ h/{\rm Mpc}$. Moreover, the prediction is guaranteed to be accurate whenever
the softer mode is sufficiently linear, ${\rm min}\{k_1,k_2\} \lesssim 0.08\
h/{\rm Mpc}$. We identify and discuss straightforward improvements in the
context of the response approach, which are expected to further increase the
accuracy of the calculation presented here.
