The computational expense of the LSC nonadiabatic dynamics stems from the computation of the electronic structure variables, E n (R), g n = ∂En(R) ∂R , and d nm (R), for all electronic states under consideration at a given nuclear geometry, R. Here, we describe how the workload of the TeraChem electronic structure calculations are distributed amongst the available GPU compute nodes used in a simulation and how they communicate with our LSC code to propagate nuclear motion. The PySCES code is available open-source on Github along with appropriate documentation and input files to test the code. Figure depicts the workflow of running an LSC simulation with the PySCES code, starting with an optimized geometry. The user provides PySCES with the location of a TeraChem normal mode analysis of the optimized structure, for which initial nuclear and electronic conditions are assigned from the Wigner distributions in Eqs. 8 and 11, respectively. Starting from this new randomized geometry, and for all subsequent geometries throughout the simulation, both the gradients g i of each adiabatic energy surface and the NAC vectors d ij between two surfaces need to be computed. The sign of the NACs are then checked and, if necessary, are corrected for random sign-flips between time steps that occur due to a change in phase of the wave function(see Section III A below). We then obtain updated nuclear and electronic variables (position and momenta) via a Runge-Kutta 4 th -order integrator as implemented in SciPy. The updated nuclear geometries are then used once again as inputs for electronic structure calculations, and the process continues. All relevant information, including the nuclear and electronic variables, gradients, NAC vectors, timings of the individual electronic structure calculations, and checkpoint files needed to restart the trajectory, are logged throughout the simulation as well.