Development of a Charge-Implicit ReaxFF for C/H/O Systems

Modeling chemical reactions in condensed phases is difficult. Interaction potentials (or force fields) like ReaxFF can perform this modeling with a high overall accuracy, but the disadvantage of ReaxFF is a low simulation speed arising from costly algorithms, in particular charge equilibration. Therefore, we reparametrized ReaxFF to incorporate Coulomb forces into other terms of the force field. Because of this change, our charge-implicit ReaxFF-CHO is >2 times faster than the original parametrization. Despite the lack of explicit electrostatic interactions, our potential can correctly model the reactions and densities of systems containing carbon, hydrogen, and oxygen atoms. We have used the new potential to simulate bombardment of trehalose by water clusters. It has been observed experimentally that these water projectiles can increase the sensitivity of secondary ion mass spectrometry by more than an order of magnitude, but no explanation for this phenomenon was given. Our simulations show that the increase in the intensity of the recorded signal coincides with the emission of trehalose–water complexes.

S2 this stage. The density and potential energy of the system (Eliquid) is averaged during the next step.
After that, the intermolecular interactions are disabled which effectively creates a gaseous system, whose potential energy (Egas) is averaged after a sufficiently long equilibration.
The exact stages of the procedure are given below, along with the time of each step.
1. Create a primitive cubic molecular system with density 10 times lower than the experimental value.
6. Average the density and potential energy of the system (Eliquid) for 100 ps.

Disable intermolecular interactions.
8. Perform an NVT simulation for 1 ps. 9. Average the potential energy of the system (Egas) for 10 ps.
The number of molecules varied between systems and was chosen so that the predicted size of the equilibrated system in every dimension was at least four times higher than the cutoff of the potential. The damping constants of Nosé-Hoover thermostat and barostat were equal to 10 fs -1 and 100 fs -1 , respectively. There were two exceptions from this procedure: the equilibration stage (5) for benzene was three times longer in order to achieve an equilibrium. Additionally, the simulation of freezing water started from the equilibrated water sample for 300 K (steps 1-4 were skipped), which was subsequently cooled to 77 K during a 200 ps simulation with Nosé-Hoover thermostat and barostat and 1 atm target pressure. Subsequently, the procedure was started from step 5.

Oxidation simulations
In order to test the new potential, an oxidation simulation was carried out for o-xylene. The modelling procedure was adapted from [2] and is given below: 1. An initial sample consisted of 100 O2 molecules and a single o-xylene in a periodic cubic box with edge length equal to 2.5 nm.
2. The temperature of the sample was increased to 2500 K with the use of a Berendsen thermostat with the damping constant equal to 0.1 ps during a 10 ps simulation. After that, the system was equilibrated in 2500 K for 100 ps. The C-O and H-O interactions were disabled, so that no reactions would occur at this step.
3. Finally, the temperature was maintained at 2500 K by a Berendsen thermostat with a damping constant equal to 0.5 ps and the simulation was run till the oxidation concluded.

Sputtering simulations
Both water and trehalose samples used in sputtering simulations were hemispherical because the pressure waves induced by a cluster impact have spherical symmetry [3]. The size, density, and the number of atoms in the samples are given in Table S1. The rigid and stochastic layers with thicknesses of 0.7 and 2.5 nm, respectively, were used for preserving the shape of the system, and for suppressing reflections of pressure waves from the sample borders [3]. The simulations were run in 0 K target temperature for 50 ps, which was enough time for achieving saturation of the sputtering yield. Every simulation was performed once as the sputtering yield for cluster projectiles only weakly depends on the impact point [4].  Figure S1 shows the comparison of ci-ReaxFF-CHO and ci-ReaxFF-CHO with tabularized correction to the ZBL potential (a, b). Without the tabulated potential, the ci-ReaxFF-CHO is in agreement with ZBL only for a narrow range of distances. This observation is especially true for the carbon-oxygen interaction. This cannot be avoided, since the ReaxFF formalism allows for specifying short-range repulsion parameters only for elements and uses automatic combination rules for pairs of different elements. Modifying this behavior would break compatibility with existing ReaxFF implementations. The addition of the tabulated potential increases the agreement to interatomic distances equal to about 0.01 nm regardless of the element. The transition from ci-ReaxFF-CHO to ci-ReaxFF-CHO with tabularized potential is smooth, thanks to the use of a splining region. An example of such region is shown in Figure S1c. The cutoffs for splining regions are given in Table S2.  Figure S1 -Comparison of close-range repulsive barriers predicted by charge-implicit ReaxFF and ZBL potential (a,b) and the splining area for the C-C interaction (c). for ReaxFF-2008 [2]. The bond dissociation curves used in the training process are in Figures S2-S4. Figures S5-S7 and Figures S8-S14 contain a comparison of valence and dihedral angle barriers, respectively. Other reactions curves are shown in Figures S15-S17. For QM data, the value in each graph is relative to a structure with minimum potential energy. The ci-ReaxFF-CHO values are relative to the potential energy, predicted by the potential of the same structure.

Simulations of trehalose bombarded by water clusters
Due to the size of the water clusters ( Figure S18a), the diameter of the trehalose sample (Figure S18b