Path integral molecular dynamics framework for spin-phonon systems

We build on the framework for studying systems containing spin-phonon coupling using molecular dynamics as the method for sampling the positions of atoms carrying spins, whose interaction is described by the Heisenberg model, and which we presented in a previous paper Srpak et al. [Phys. Rev. B 112, 205132 (2025)]. In applying the real-space path integral framework to the spin-phonon problem, neighboring imaginary time slices become coupled via overlap of eigenvector matrices propagating along the path integral, leading to interactions that go beyond the pure harmonic interactions between neighboring beads of the ring-polymers in Born-Oppenheimer-based path integral molecular dynamics (PIMD) simulations of quantum nuclei. We show how the use of a spin-force, defined as the negative gradient of the spin free energy, accounts for this type of coupling, valid for systems in thermal equilibrium. The spin force is computed using either ring-polymer exact diagonalization (for small systems) or stochastic series expansion (for large systems). We provide the technical details for these implementations. We present a neural network approach for this type of PIMD simulation, where we compute the spin force based on training data obtained from the ring-polymer stochastic series expansion method. Utilizing these methods, we compare them with the literature and demonstrate how one can obtain phase diagrams clearly showing how the thermal and nuclear quantum effects influence the stability of the spin-Peierls dimerized phase. We also show how machine learning makes it possible to apply the method to two dimensional systems.