A parallel multistate framework for atomistic non-equilibrium reaction dynamics of solutes in strongly interacting organic solvents

We describe a parallel linear-scaling computational framework developed toimplement arbitrarily large multi-state empirical valence bond (MS-EVB)calculations within CHARMM. Forces are obtained using the Hellman-Feynmannrelationship, giving continuous gradients, and excellent energy conservation.Utilizing multi-dimensional Gaussian coupling elements fit to CCSD(T)-F12electronic structure theory, we built a 64-state MS-EVB model designed to studythe F + CD3CN -> DF + CD2CN reaction in CD3CN solvent. This approach allows usto build a reactive potential energy surface (PES) whose balanced accuracy andefficiency considerably surpass what we could achieve otherwise. We use our PESto run MD simulations, and examine a range of transient observables whichfollow in the wake of reaction, including transient spectra of the DFvibrational band, time dependent profiles of vibrationally excited DF in CD3CNsolvent, and relaxation rates for energy flow from DF into the solvent, all ofwhich agree well with experimental observations. Immediately followingdeuterium abstraction, the nascent DF is in a non-equilibrium regime in twodifferent respects: (1) it is highly excited, with ~23 kcal mol-1 localized inthe stretch; and (2) not yet Hydrogen bonded to the CD3CN solvent, itsmicrosolvation environment is intermediate between the non-interactinggas-phase limit and the solution-phase equilibrium limit. Vibrationalrelaxation of the nascent DF results in a spectral blue shift, while relaxationof its microsolvation environment results in a red shift. These two competingeffects result in a post-reaction relaxation profile distinct from thatobserved when DF vibration excitation occurs within an equilibriummicrosolvation environment. The parallel software framework presented in thispaper should be more broadly applicable to a range of complex reactive systems.