SAFIRES - A flexible boundary method for QMMM calculations with ASE

This MR is for a new feature we would like to introduce to ASE: the Scattering-Adapted Flexible Inner Region Ensemble Separator (SAFIRES) as presented in our recent paper.[1] SAFIRES is a restrictive boundary method for hybrid calculations (e.g. QMMM) of molecular dynamics (MD) of a solute in a solvent environment. The goal of SAFIRES is to enable us to study the properties of the model system in contact with a large number of explicit solvent molecules, i.e. in a realistic solvation environment, see illustration in Figure 1.

Figure 1: Illustration of the QMMM boundary in a simulation using a bulk water model (left) and of water in contact with a periodic surface model (right).

When performing QMMM-MD simulations of a solute ion or molecule, or a periodic model surface, in contact with water, a key issue is the treatment of the boundary between QM and MM regions (we assume QM and MM instead of a general hybrid simulation approach from hereon for simplicity). Often it is a good choice to include the solute as well as the first solvation shell in the QM region; this way, electronic interactions across the boundary only occur between solvent molecules, thus simplifying the interaction terms. However, since solvent molecules are mobile, the question then becomes how to treat boundary events, i.e. events where a solvent molecule from the QM region enters the MM region or vice versa.

One possibility is to switch the description between QM and MM as molecules travel into the opposite region. This approach is commonly referred to as adaptive QMMM. Ad-hoc switching of the descriptions however introduces discontinuities to the energy and forces. To smooth over these discontinuities, a large number of additional energy calculations is required, thereby increasing the computational cost by orders of magnitude.[2]

Another possibility is to restrict exchange of particles across the boundary (restrictive boundary methods). Unlike in adaptive QMMM methods, QM and MM particles cannot change description in restrictive models and no additional calculations are required to conserve the energy and forces. However, dynamics can only be studied deep inside the QM region in this approach due to the redirection of particles at the boundary. In 2012, Rowley and Roux presented the FIRES boundary method which uses a harmonic potential to contain QM and MM particles within their respective regions.[3] Although Rowley and Roux show mathematical proof that correct average thermodynamic properties can in principle be obtained using this restrictive method despite trajectories of molecules being artificially altered, the FIRES method was still found to introduce artifacts around the location of the boundary.[4]

Our group recently developed SAFIRES, a scattering-adapted implementation of FIRES. Here, boundary events are restricted not by means of a harmonic potential but by performing an elastic collision between the QM and MM particles which is mediated by the boundary. The scattering approach required implementation of a variable-time-step propagator which can identify and propagate to the exact moment in the dynamics where a collision would occur. SAFIRES is implemented as a dynamics class derived from the Langevin class. The propagator reduces to the Langevin propagator for constant time steps and further to velocity Verlet dynamics for constant time steps and zero friction. SAFIRES conserves energy and forces, and forces are only updated at full default time step intervals to ensure temporal consistency. In benchmark calculations using the same Lennard-Jones potential for both regions we have shown that SAFIRES is able to recover the radial distribution function (RDF) of a simulation without boundary while FIRES introduces artifacts at the boundary, see Figure 2.

Figure 2: RDFs of an Ar Lennard-Jones liquid from a simulation without any boundary (black dotted line) as well as from simulations using the FIRES (red) and SAFIRES (blue) boundary methods. RDFs are sampled between a central "solute" Ar particle and all other Ar particles to allow us to give distributions for the location of the FIRES (grey solid line) and SAFIRES boundary (grey dashed-dotted line).

Here, we would like to make the code available for all ASE users. This merge request contains:

• The SAFIRES dynamics class, which constitutes a flexible-time-step Langevin propagator and has been derived from the existing Langevin class, located in /ase/md/safires.py.
• Documentation of the SAFIRES class and its parameters, located in /doc/ase/md.rst.
• Tutorial documentation and executable tutorial script file on how to use the SAFIRES class, located in /doc/tutorials/safires.
• Unit test for the SAFIRES class, located in /ase/test/test_safires.py.
• Changes to /ase/md/logger.py to log SAFIRES-specific output (location of boundary, number of resolved boundary events).

We're looking forward to feedback on the MR and improving the code along the way.

Cheers,
Björn Kirchhoff

[1] B. Kirchhoff, E. Ö. Jónsson, A. O. Dohn, T. Jacob, H. Jónsson, J. Chem. Theory Comput. 2021, 17, 5863–5875.
[2] H. C. Watanabe, Q. Cui, J. Chem. Theory. Comput. 2019, 15, 3917–3928.
[3] C. N. Rowley, B. Roux, J. Chem. Theory. Comput. 2012, 8, 3526–3535.
[4] R. E. Bulo, C. Michel, P. Fleurat-Lessard, P. Sautet, J. Chem. Theory Comput. 2013, 9, 5567–5577.