  We present an implementation of the quantum mechanics/molecular mechanics
(QM/MM) method for periodic systems using GPU accelerated QM methods, a
distributed multipole formulation of the electrostatics, and a pseudo-bond
treatment of the QM/MM boundary. We demonstrate that our method has
well-controlled errors, stable self-consistent QM convergence, and
energy-conserving dynamics. We further describe an application to the catalytic
kinetics of chorismate mutase. Using an accurate hybrid functional
reparametrized to coupled cluster energetics, our QM/MM simulations highlight
the sensitivity in the calculated rate to the choice of quantum method, quantum
region selection, and local protein conformation. Our work is provided through
the open-source \textsc{PySCF} package using acceleration from the
\textsc{GPU4PySCF} module.
