  A second-order many-body perturbation correction to the relativistic
Dirac-Hartree-Fock energy is evaluated stochastically by integrating
13-dimensional products of four-component spinors and Coulomb potentials. The
integration in the real space of electron coordinates is carried out by the
Monte Carlo (MC) method with the Metropolis sampling, whereas the MC
integration in the imaginary-time domain is performed by the inverse-CDF
(cumulative distribution function) method. The computational cost to reach a
given relative statistical error for spatially compact but heavy molecules is
observed to be no worse than cubic and possibly quadratic with the number of
electrons or basis functions. This is a vast improvement over the quintic
scaling of the conventional, deterministic second-order many-body perturbation
method. The algorithm is also easily and efficiently parallelized with
demonstrated 92% strong scalability going from 64 to 4096 processors for a
fixed job size.
