High efficiency and accuracy phase gratings are of crucial importance for large format heterodyne array receivers at terahertz frequencies. Here, by developing a design approach that can create gratings with arbitrary two-dimensional diffraction distributions, we have realized a reflective metallic phase grating that generates 2×2 diffraction beams at 0.85 THz. The measured total power efficiency of the diffraction beam pattern is 81.9%, which demonstrates at least 17% improvement in efficiency compared with the standard pseudo-2D Fourier phase grating. In addition, we report the realization of up to 10×10 diffraction beam two-dimensional phase grating designs at terahertz wavelengths, using an adaptation of the Gerchberg-Saxton (GS) scheme known as the Mixed-Region-Amplitude-Freedom algorithm. Rigorous full wave simulation proves the efficiency and accuracy of the design, which overcomes the inaccurate intensity of the beam distribution drawbacks originated from the standard GS algorithm. The results pave the way for the development of large-pixel terahertz multi-beam heterodyne receivers.