We present distribution functions and mark correlations of the shapes of massive dark matter halos derived from Hubble volume simulations of a ?CDM universe. We measure both position and velocity shapes within spheres that encompass a mean density 200 times the critical value and calibrate small-N systematic errors using Poisson realizations of isothermal spheres and higher resolution simulations. For halos more massive than 3 ? 1014 h-1 M?, the shape distribution function peaks at (minor/major, intermediate/major) axial ratios of (0.64, 0.76) in position and is rounder in velocity, peaking at (0.72, 0.82). Halo shapes are rounder at lower mass and/or redshift; the mean minor-axis ratio in position follows c/a(M,z) = c15,0[1 - ? ln(M/1015 h-1 M?)](1 + z)-, with c15,0 = 0.631 ? 0.001, ? = 0.023 ? 0.002, and = 0.086 ? 0.004. Position and velocity principal axes are well aligned in direction, with median alignment angle 22?, and the axial ratios in these spaces are correlated in magnitude. We investigate mark correlations of halo pair orientations using two measures: a simple scalar product shows ?1% alignment extending to 30 h-1 Mpc, while a filamentary statistic exhibits nonrandom alignment extending to scales ~200 h-1 Mpc, 10 times the sample two-point correlation length and well into the regime of negative two-point correlation. Shapes of cluster halos are little affected by the large-scale environment; the distribution of supercluster member minor-axis ratios differs from that of the general population at only the few percent level.