The discrete element method (DEM) is the most popular numerical method for simulating granular materials; however, the high computational costs of DEM limit the complexity and size of engineering problems that can be solved. Simulation of granular soil is particularly challenging because soil particles have complex shapes that cannot be simplified to spheres or ellipsoids without great care. Some DEM codes have introduced artificial rotational moment of inertia or clumping of spherical particles to account for particle shape effects indirectly. An alternative to particle shape simplification is explicit representation of soil particle shapes with polyhedra. However, detecting contacts of polyhedral shapes is computationally expensive and small integration time steps are required for numerical stability. Here, we summarize the functionality and applications of BLOKS3D, a research DEM code focused on efficient simulation of soil using polyhedral shapes. By using direct representation of polyhedral shapes, the contractive and dilative behavior of soils can be simulated using a simple spring-dashpot-slider collision model. Conventional DEM uses an explicit scheme to double integrate acceleration (forces) with a maximum time step limited to ensure numerical stability. In computer graphics, particle simulations are commonly performed using impulse calculations without computing contact forces. BLOKS3D introduces an impulse-based DEM (iDEM) scheme that can also retrieve contact forces. This method allows engineers to harness the speed of the computer graphics algorithm while maintaining access to the forces and stresses needed for design. The iDEM formulation has a less restrictive requirement on maximum time step and prior work has shown that an iDEM time step of approximately 100 times that of conventional DEM produces similar results. This translates to a runtime speedup of two orders of magnitude. BLOKS3D was parallelized to further speedup simulations. Contact detection constitutes roughly 70% of the serial computational cost, and therefore, the contact detection steps were the focus of the parallelization. A two-level contact detection algorithm was implemented to improve contact detection efficiency. First, a rapid algorithm identifies possible contacts by discretizing the space domain into cubes. Only particles in the same or neighboring cubes are passed to the second algorithm. The second algorithm then uses the rigorous shortest link method to identify final contacts. The parallelized two-level contact detection algorithm for shared memory systems can be implemented on systems ranging from HPC platforms to engineering workstations using the OpenMP framework. The speedup obtained from parallelization is dependent on the number of computing threads and the number of particles. For example, simulations performed with 48 threads were approximately 30 times faster than simulations in series for particle packing simulations of 100,000 polyhedral particles with the conventional DEM formulation. The iDEM formulation and the parallelized two-level contact detection can be employed together. Combining these improvements led to the first simulation directly modeling more than 50 million explicit polyhedral DEM particles.