Abstract Shear transformation zone dynamics models of metallic glass deformation access experimentally-relevant time scales by using the kinetic Monte Carlo method to simulate small, fast, often discrete events, while the finite element method calculates macroscopic shape change and continuum-level interactions within samples. The most time-consuming portion of these models is the finite element method calculation on each step. However, in cases where the finite element mesh geometry and element elastic properties do not change from step to step, the finite element stiffness matrix (and its Cholesky factors) from previous steps can be reused. This strategy improves the asymptotic complexity of these models and in practice accelerates their execution by nearly 200 $$\times $$ × . This enables simulation of larger samples in more reasonable time. A set of three-dimensional shear transformation zone dynamics simulations, with larger length scales than any currently in the literature, illustrates the utility of this approach.