In this manuscript, we describe a new systematic workflow for generating optimized geometries along grids of torsion constraints by wavefront propagation, which addresses the drawbacks of serial relaxed scanning described above. The method, called TorsionDrive, generates results that are independent of scan direction and naturally incorporates multiple initial guesses into a single grid of constrained minimized structures. A predecessor to the present method was used to scan the 2-dimensional torsion angles of blocked amino acid dipeptides in Ref. 8, but due to limitations of implementation, the method was limited to 2-D grids and was not easily applicable to other molecules. The present method can be applied to any molecule subject only to the limits of the underlying energy and gradient method. Furthermore, it is capable of driving an arbitrary number of torsions to generate N -dimensional grids of optimized structures and energies where N ≥ 1. The current workflow is implemented as a Python package that interfaces with energy minimization routines in a modular way, including the open-source geomeTRIC optimization package that uses externally obtained energy and gradient information, as well as "native" optimization methods implemented in many quantum chemistry codes. The method is released as an open source package 32 and includes a number of useful features such as including energy upper limits, extra constraints, and limited scan ranges. In addition to standalone operation mode, TorsionDrive is also implemented as a service in the MolSSI QCArchive ecosystem, and it is available to compute results for any implemented gradient method in QCArchive; this includes not only quantum chemistry methods, but also some MM force fields and recently developed neural network potentials parameterized by machine learning methods. The main idea of TorsionDrive can be described conceptually as scanning the torsion angles with wave-front propagation. The details are illustrated by walking through a complete scan procedure, shown in Figure . Before the start of the scan, we specify the dihedral angles in the molecule of interest using quartets of atomic indices, and the spacing (resolution) of the scan. A grid of constraint values is created which has the same dimension as the number of dihedral angles provided. In the illustration, we perform a 1-D scan with a 60 • spacing. The data associated with each grid point is represented by a circle, which consists of one or more "optimization datasets" (i.e. the Cartesian coordinates of a constrained local minimum, and corresponding QM energy and gradient). Importantly, each grid point is able to contain multiple constrained optimization datasets that may correspond to local minima with different energies and values of the orthogonal degrees of freedom.