The relative dispersion of one fluid particle with respect to another is fundamentally related to the transport and mixing of contaminant species in turbulent flows. The most basic consequence of Kolmogorov's 1941 similarity hypotheses for relative dispersion, the Richardson-Obukhov law that mean-square pair separation distance grows with the cube of time at intermediate times in the inertial subrange, is notoriously difficult to observe in the environment, laboratory, and direct numerical simulations (DNS). Inertial subrange scaling in size parameters like the mean-square pair separation requires careful adjustment for the initial conditions of the dispersion process as well as a very wide range of scales (high Reynolds number) in the flow being studied. However, the statistical evolution of the shapes of clusters of more than two particles has already exhibited statistical invariance at intermediate times in existing DNS. This invariance is identified with inertial-subrange scaling and is more readily observed than inertial-subrange scaling for seemingly simpler quantities such as the mean-square pair separation Results from dispersion of clusters of four particles (called tetrads) in large-scale DNS at grid resolutions up to 4096 points in each of three directions and Taylor-scale Reynolds numbers from 140 to 1000 are used to explore the question of statistical universality in measures of the size and shape of tetrahedra in homogeneous isotropic turbulence in distinct scaling regimes at very small times (ballistic), intermediate times (inertial) and very late times (diffusive). Derivatives of fractional powers of the mean-square pair separation with respect to time normalized by the characteristic time scale at the initial tetrad size constitute a powerful technique in isolating cubic time scaling in the mean-square pair separation. This technique is applied to the eigenvalues of a moment-of-inertia-like tensor formed from the separation vectors between particles in the tetrad. Estimates of the proportionality constant "g" in the Richardson-Obukhov law from DNS at a Taylor-scale Reynolds number of 1000 converge towards the value g=0.56 reported in previous studies. The exit time taken by a particle pair to first reach successively larger thresholds of fixed separation distance is also briefly discussed and found to have unexplained dependence on initial separation distance for negative moments, but good inertial range scaling for positive moments. The use of diffusion models of relative dispersion in the inertial subrange to connect mean exit time to "g" is also tested and briefly discussed in these simulations. Mean values and probability density functions of shape parameters including the triangle aspect ratio "w," tetrahedron volume-to-gyration radius ratio, and normalized moment-of-inertia eigenvalues are all found to approach invariant forms in the inertial subrange for a wider range of initial separations than size parameters such as mean-square gyration radius. These results constitute the clearest evidence to date that turbulence has a tendency to distort and elongate multiparticle configurations more severely in the inertial subrange than it does in the diffusive regime at asymptotically late time. Triangle statistics are found to be independent of initial shape for all time beyond the ballistic regime. The development and testing of different schemes for parallelizing the cubic spline interpolation procedure for particle velocities needed to track particles in DNS is also covered. A "pipeline" method of moving batches of particles from processor to processor is adopted due to its low memory overhead, but there are challenges in achieving good performance scaling.