Tribochemical Reactions of MoDTC Lubricant Additives with Iron by Quantum Mechanics/Molecular Mechanics Simulations

The remarkable lubricant properties of molybdenum dithiocarbamates (MoDTCs) make this class of oil additives well-known in the automotive industry. However, the mechanism of function of these compounds is still not completely understood at the atomistic level. We provide new insights into the dissociation of MoDTCs in tribological conditions, which are the key to describe the debated mechanism to form MoS2. Quantum mechanics/molecular mechanics (QM/MM) dynamic simulations allowed us to monitor in real time the tribochemical reactions occurring at the iron interface and revealed that the presence of the iron substrate and the mechanical stresses alter the dissociation path with respect to what is expected for the isolated MoDTC molecules. Moreover, they uncovered the important role of molecular oxidation on the dissociation pattern: the presence of oxygen atoms in the ligand position of MoDTCs favors the release of the central units of the molecules, containing just Mo and S atoms with the correct stoichiometry to form MoS2. This work demonstrates how the predictive power of ab initio simulations can be very valuable to design new lubricant additives.


■ INTRODUCTION
In the pursuit of more sustainable solutions to the environmental challenges, friction modifiers play a fundamental role in attempting to minimize energy losses. Organo-molybdenum compounds are among the most successful friction modifiers in the boundary lubrication regime, where the asperity contact at the steel/steel interface generates extreme tribological conditions. 1,2 The success of molybdenum dithiocarbamate (MoDTC) additives in reducing the friction in automotive applications is, in fact, well-established. It is known that the beneficial effect of this class of organo-molybdenum compounds resides in the possibility to generate layers of molybdenum disulfide (MoS 2 ) in tribological conditions. 3 Recent works regarding MoDTC complexes have demonstrated their capability to reduce friction also at metal/ diamond-like−carbon and metal/rubber contacts, 4−6 and promising applications have been described in combination with textured metallic surfaces 7 and nanoparticles. 8,9 While wide experimental research has been carried out to describe the chemical and tribological behavior of MoDTC in engine oils and its interaction with ZnDTP, 10−13 little attention has been paid to the fundamental steps of the tribochemical reaction that forms the beneficial MoS 2 tribolayers. The first elementary step of this reaction is the dissociation of the molecule into two separate fragments. Grossiord et al. characterized the formation of the MoS 2 layers by means of analytical tribology and proposed a mechanism for the dissociation of MoDTC. 14 They suggested that these compounds undergo a homolytic breaking between the sulfur atoms of the dithiocarbamate (DTC) units and the molybdenum atoms. Khaemba et al., on the contrary, proposed that the bond rupture occurs between the carbon and sulfur atoms of the DTC, leaving the S atoms attached to molybdenum. 15 The actual mechanism of the dissociation of MoDTC is still debated to this day. In a previous study, our simulations attempted to quantify the strength of the bonds broken at the early stages of the reaction. We found that for isolated complexes the most favorable way of cutting the molecule in two fragments is between S and Mo atoms, in agreement with Grossiord's proposition. 16 However, the presence of a metallic substrate may modify this picture because the stabilization effect provided by the iron atoms can vary among the different fragments obtained in the dissociation. Onodera et al. investigated the interaction between one of the possible structures of MoDTC and a (001) iron surface by means of a hybrid quantum chemical/ classical molecular dynamics approach based on the tightbinding approximation. 17 They observed the elongation of the bond between Mo and the O atoms in ligand position, i.e., between the C atom of the DTC and molybdenum, according to the notation of Khaemba,15 and they identified the role of iron as a catalyst in the dissociation of MoDTC. However, this bond stretching was observed in a few hundreds of femtoseconds, which is a short time interval for drawing conclusions on the fundamental steps of the molecular decomposition, as we show in detail in the following paragraphs.
In this work, we investigated the dissociation of MoDTC by considering three dimeric structures, in which the position of O atoms and the ratio between S and O atoms vary. In a previous study, our simulations showed that the ligand position, between C and Mo atoms, is the most favorable position for the O atoms, and we found that the oxidation of the complexes is generally favorable. 16 These results were in agreement with the experimental observations of De Feo et al., 18 who demonstrated that MoDTC structures become spontaneously rich in oxygen, which eventually occupies the ligand position, in addition to the terminal position, i.e., on top of Mo, due to a multistep reaction. Therefore, it is worth comparing complexes with different positions and ratios between S and O atoms to better understand the role of oxygen in the dissociation process of MoDTC.
Here, the interaction of these complexes with a metallic substrate was investigated first by performing dynamic simulations following a quantum mechanics/molecular mechanics (QM/MM) approach 19,20 based on the density functional theory (DFT). Such a computational scheme allows an accurate description of the chemical reactions while maintaining the calculations computationally affordable. The region in which the chemical reactions take place is treated at the quantum mechanical level by explicitly taking into account the electrons of the reacting chemical species. The electronic structure of the atoms far away from the interface can be neglected, as they do not take part in any chemical reaction. However, the presence of these atoms is necessary to absorb part of the atomic displacements and oscillations generated by the reactions at the interface. Therefore, they are assigned to the classical region, as they are treated by means of classical mechanics. Our recent studies showed that the QM/MM approach can be very promising in tribochemistry, as it properly describes the chemical reactions occurring at the sliding interface with a significant speedup with respect to simulations where all the atoms are treated at the quantum level. Here we extend the field of application of the QM/MM approach to investigate the mechanism of function of lubricant additives, which is rather challenging because of the complexity of the systems and the processes involved. 21,22 After describing the tribological behavior of the three MoDTC complexes by means of the QM/MM approach, we carried out additional static simulations to confirm and quantify the results obtained from the dynamic simulations. We calculated the adsorption energy of the complexes and their dissociation energy on the clean iron surface. We considered iron, which is different from steel, as it can be an appropriate model for the native metallic surfaces exposed during scratching. 23 Our aim is to clarify the dissociation mechanism of MoDTC, describe quantitatively the interaction of these complexes with an iron surface, and understand its catalytic effect to shed some new light on the tribochemical reactivity of organo-molybdenum friction modifiers.

■ METHODS AND SYSTEMS
In this work, we considered the three MoDTC complexes shown in Figure 1, which will be referred to as standard MoDTC (sMoDTC), Isomer H, and Substitution L, by following the nomenclature adopted in our previous publication, where these complexes were extensively described. 16 The metallic substrate was simulated by considering iron slabs, with the molecules interacting with the (110) surface, as it is the most stable surface for iron. 24 The QM/MM calculations contained two sMoDTC (Simulation A) or two Substitution L (Simulation B) complexes per 7 × 3√2 simulation cell, intercalated between the (110) iron slabs, each composed by eight atomic layers. The two molecules and two iron layers, one for the top and one for the bottom iron slabs, belonged to the quantum region, while the classical region was composed of the remaining 14 iron layers, seven for each iron slab. Figure 2 shows a schematic representation of the QM/MM system, which was prepared by following a procedure explained in the Supporting Information. Inside the chemically reactive region, the interactions among the atoms were treated at the DFT level. In this region, the electronic structure of the atoms was optimized at each step of the simulation by Quantum ESPRESSO package, 25 effectively making this computational approach a Born−Oppenheimer molecular dynamics technique. 26 The interactions among the atoms in the classical region and at the border between the quantum and the classical regions were described by the embedded-atom method, 27,28 with the overall dynamics handled by the LAMMPS suite. 29 The total energy of the system E tot was obtained through a subtractive scheme 30 = where E tot,MM is the energy of the whole system computed classically and E crr,MM and E crr,QM are the energies of the chemically reactive region computed classically and at the DFT level, respectively. By following this scheme, the forces of the whole system are first computed by means of classical force fields. Then, the forces of the chemically reactive region are calculated at the quantum level and substituted within the same region. In this way, it is possible to maintain an easy implementation, which does not require any further assumption to couple the quantum and the classical regions, as the forces at the boundary are provided by the force fields.
In the simulations, a thermalization process, lasting 1 ps, allowed the nuclei to globally reach 300 K thanks to a Nose− Hoover thermostat with a damping parameter of 33.3 fs. An initial ensemble of velocities for all of the nuclei was generated from a Gaussian distribution with a standard deviation scaled to produce a temperature of 300 K. After 1 ps simulation, the temperature of the quantum region was left to evolve freely, while the classical atoms continued to be thermostatted, absorbing the heat generated at the interface. This approach can be considered as an intermediate scheme between a base and a full thermostatting. 31 Lateral motion of the iron slabs along the x-axis was introduced to simulate the sliding of the two surfaces. A velocity kick of 300 m/s in two opposite directions for the top and the bottom slabs was provided to the outermost iron layers at each picosecond of the simulation to counteract the interfacial friction. This value for the sliding velocity was chosen to accelerate the reaction rates to observe the desired chemical events in reasonable simulation time. The outermost iron layers also experienced vertical forces that impose a load of 4 GPa to the system for the whole duration of the simulation. Such a value for the load is quite large compared to the usual experimental conditions. However, the ideal surface considered in the calculations is meant to represent a nanoasperity contact, where the contact pressures can reach values way above the parameters imposed in the experimental setups. 22,32,33 Furthermore, the chosen value for the load is still much lower than the Young's modulus of iron, ensuring no failure of the material in these conditions. The time step of the simulation was 1 fs, and periodic boundary conditions in all directions were employed for the supercells, as the current implementation of the QM/MM scheme requires the same simulation cell and periodicity conditions for the quantum and classical subsystems. A cell height of approximately 48 Å allowed us to obtain a separation of about 15 Å between the vertical periodic replicas of the whole QM/MM system. All of the DFT calculations in this work were carried out by employing the Perdew−Burke−Ernzerhof (PBE) 34 approximation for the exchange-correlation functional. The pseudopotential/plane wave approach was adopted. Ultrasoft pseudopotentials were employed, and the plane wave expansion of the electronic wave function (charge density) was truncated using a 30 Ry (240 Ry) cutoff for the kinetic energy. The value of 30 Ry for the kinetic energy cutoff of the wave functions was chosen to obtain an accurate description of the iron bulk, as shown in previous work from our group, 24 while the energy difference between two isolated MoDTC structures was at convergence below the threshold of 1 meV per atom, as explained in the Supporting information. Spin polarization was taken into account. We added Gaussian smearing of 0.02 Ry to better describe the occupation of the electronic states of the metal around the Fermi energy.  Figure 3, in a top view representation of the interface, with the iron atoms hidden for clarity. Both sMoDTC and Substitution L are represented. Panels A and B are directly comparable with panels D and E, respectively, as they correspond to the configurations assumed by the two structures at the beginning of the simulation and after 1 ps. Panels C and F show the instant in which full dissociation of the two compounds was achieved. Initially, Substitution L appeared to experience a quicker dissociation than sMoDTC, as two carbamate units were already separated from the central unit after 1 ps. However, full dissociation of the DTC units was achieved more quickly for sMoDTC (1.8 ps) rather than Substitution L (3.  Adsorption and Dissociation on Iron. To verify and quantify the results of the QM/MM simulations, we calculated the adsorption and fragmentation energies of sMoDTC, Isomer H, and Substitution L on iron. First, the complexes were put close to the surface, with the terminal atoms pointing to the iron atoms. The result of the geometry optimizations, visible in Figure 4A−C, shows that the position of oxygen in the molecule plays a role in the chemisorption process. The sMoDTC molecule, with oxygen atoms in terminal position, is not able to fully lie on the surface, as three sulfur atoms on one side of the central unit remain distant from the iron surface. For the other two complexes, there is no such hindrance, and five (for the isomer) and four (for the substitution) sulfur atoms reach the surface. In particular, two sulfur atoms of Substitution L detach from the central unit in an attempt to occupy the hollow sites of the iron surface. This phenomenon is reflected upon the adsorption energies E ads , calculated as where E tot , E mol , and E iron are the total energies of the systems composed of the molecule and the iron slab at the end of the optimization process, the molecule alone, and the iron slab alone, respectively. The adsorption energies turned out to be −4.33, −7.84, and −7.17 eV for sMoDTC, Isomer H, and Substitution L, respectively. All of the calculated values of the adsorption energies are compatible with chemisorption, with the isomer and the substitutional molecule achieving the highest energy gain upon adsorption.
To verify which is the most probable dissociation pattern for sMoDTC, Isomer H, and Substitution L on iron, we calculated the total energies, after geometry optimization, of their fragments adsorbed on iron, in the same supercell employed for the adsorption calculations. We considered 15 different fragments, visible in Figure 4, to take into account two steps of the fragmentation, i.e., the dissociation of the first and then the second DTC unit. We considered two possible ways of dissociating the complexes by following the same nomenclature employed in our previous work: Cut 1, in which both O or S atoms in ligand position remain on the ligand after the dissociation, and Cut 3, in which the O or S atoms remain attached to Mo. Cut 2, which is not considered here, is the intermediate situation in which one of the two atoms in the ligand position remains with the ligand unit, while the other atom remains attached to the Mo atom. Due to the asymmetry of Isomer H, there are two possibilities to perform the two dissociative steps: first the carbamate unit and then the dithiocarbamate unit are detached (O→S) or vice versa (S→ O). We calculated the dissociation energy as the difference between the total energies of the fragments adsorbed on iron and the entire molecule adsorbed on iron, as explained more extensively in the Supporting Information. We summarize the results of this analysis in Table 1.
For sMoDTC, both Cut 1 and Cut 3 are favorable on the iron surface, with Cut 3 appearing to be the most favorable fragmentation pattern. In fact, while the first step of the dissociation is energetically very similar for Cut 1 and Cut 3, the products obtained from Cut 3 are much more stable than those from Cut 1 after the second step. This is in contrast with the situation of an isolated complex, for which Cut 1 is the most favorable dissociative pattern. 16 The iron surface plays a key role in stabilizing the products of the dissociation and favors a fragmentation mechanism that would be the least favorable for the isolated compound. Therefore, we expect to  The Journal of Physical Chemistry C pubs.acs.org/JPCC Article observe the C atom of the dithiocarbamate unit, which is usually attached to the S atoms, directly chemisorbed on the iron surface. The two S atoms belonging to the ligand position should remain close to the Mo atom, when the central unit is adsorbed on iron.
The most favorable fragmentation appears to be Cut 3 for Isomer H as well. In particular, when the first unit to detach is a carbamate unit, Cut 3 is just barely favorable in the first step, while the second step is about 10 times more favorable than the first one. When the dithiocarbamate unit is detached first by following the pattern of Cut 3, the energy values are simply reversed in order with respect to the previous case. For Cut 1, the situation appears to be more complex. The dissociation of a (CH 3 ) 2 NC unit, leaving two oxygen atoms attached to the Mo atom, is not favorable, while the subsequent step is favorable. When the dithiocarbamate unit is detached first, the energy values are slightly different from the latter case due to the different interactions arising between the atoms of the central units with dangling bonds and the metallic substrate.
Cut 1 is the most favorable way of detaching the first carbamate unit of Substitution L, in contrast with sMoDTC and Isomer H. It is preferable, however, to detach the second carbamate unit with Cut 3, as the energy cost is lower. This reveals that oxygen atoms in MoDTC can play an important role in governing the dissociation process of the (dithio)carbamate units. The presence of O atoms in the ligand position makes Cut 1 more favorable than Cut 3 for the dissociation of the first carbamate unit of Substitution L on iron. As expected, the metallic substrate has a deep impact on the stabilization of the resulting fragments. While the most favorable fragmentation of isolated sMoDTC is Cut 1, the same complex on iron most likely undergoes dissociation by following Cut 3. These results fully support the predictions offered by the dynamic simulations and suggest that the two proposed mechanisms for the dissociation of MoDTC are both valid, as the most probable fragmentation pattern varies by changing the ratio between S and O atoms in ligand position. The mechanism proposed by Grossiord is in fact most likely to occur when O atoms are in ligand position, while the mechanism proposed by Khaemba is most probable when S atoms are in ligand position.
Discussion. The QM/MM calculations showed that the dissociation of MoDTC complexes can be completed in the first picoseconds of the simulation. Such kinetics undoubtedly depends on the tribological conditions of the system, as the rate of the chemical reactions is strongly influenced, for instance, by temperature and the applied load. Sufficient simulation time must be considered, though, to capture the whole dissociation process and draw conclusions on the possible products of the tribochemical reactions. Onodera et al. observed a Mo−O bond elongation in 100 fs of a tightbinding-based hybrid quantum chemical/classical molecular dynamics simulation. 17 However, we confirmed that Mo−O stretching vibration is found at 1000 cm −1 in the vibrational spectrum, 16 corresponding to a period of around 33 fs. This means that Mo−O bond vibrates only a few times in such a short timescale, even in tribological conditions. Therefore, our simulation time seems to be more adequate to describe the bond breaking between ligand and central units of MoDTC. Furthermore, the hybrid approach employed by Onodera and co-workers relies on the calculation of the properties of the quantum region by means of a tight-binding technique, which, in turn, requires several parameters. Such parameters are determined by DFT calculations. The advantage of our approach is the direct computation of all of the chemical properties of the tribological system at the DFT level without further parametrization.
Our QM/MM simulations provided reliable results even using one layer per iron slab in the quantum region, which is a rather severe approximation. This choice was also made by Onodera et al. in their study, and we additionally estimated the extent of such approximation on the adsorption of sMoDTC. It turns out that by adsorbing sMoDTC on top of an iron bilayer, with the bottom layer fixed, an additional sulfur atom in the bridging position of the complex adsorbs on the surface, while the geometry of the rest of the system remains almost identical. The adsorption of sMoDTC on an iron trilayer, with the bottom layer fixed, leads to a molecular geometry identical to that in the case of the bilayer, with adsorption energy that is 1 eV higher, hence less favorable, in the case of the trilayer. With fewer iron layers, the surfaces become more reactive, hence stickier. However, we reasonably believe that the trajectories obtained by the QM/MM simulations are qualitatively consistent with the predictions offered by the static calculations, and no spurious products were generated due to the approximation of the iron slab to a combination of a single quantum and seven classical atomic layers. In fact, all of the reactants and the products of dissociation can bind with increased adsorption energy to the stickier iron surface, and the resulting trajectory is simply accelerated with respect to the case of a thicker iron slab.
Our QM/MM simulations revealed the important role of oxygen in governing the dissociation mechanism of MoDTC. De Barros Bouchet and co-workers described the decrease in the lubricating performance due to the oxidation of MoDTC, which has an impact on the number, the shape, and the composition of the MoS 2 sheets. 35 They found that the presence of ZnDTP allows MoDTC to maintain sufficient friction reduction properties even when 50% of the additive is degraded by oxidation because of the possibility to replenish the S atoms needed to build MoS 2 . Our simulations support those findings, as Substitution L, which is a slightly oxidized MoDTC compound, appears to be a better candidate in producing small units of molybdenum disulfide with proper stoichiometry. Some authors pointed out that the presence of oxygen atoms in MoS 2 sheets is unavoidable in the synthetic process and by mere exposure to ambient conditions. 36,37 The lubricant properties of pristine MoS 2 , however, is unrivaled with respect to layers of MoS 2 containing oxygen also. Therefore, we believe that the controlled oxidation of This study demonstrates the merits of ab initio calculations and, particularly, the QM/MM approach in tribochemistry to identify efficient lubricant additives. The insights provided by these simulations can be valuable to predict new solutions for the general problem of friction reduction.

■ CONCLUSIONS
In this work, we presented the application of a recently developed QM/MM technique to the dissociation of MoDTC, which is a problem of key relevance in understanding the mechanism of function of Mo-based lubricant additives. The QM/MM simulations demonstrated that the dissociation mechanism of MoDTC depends on the level of oxidation of the complex, as compounds with oxygen atoms in ligand position tend to generate central units with the same chemical composition and stoichiometry of MoS 2 . Static calculations based on DFT confirmed the predictions of the dynamic calculations, as the different dissociation mechanisms reflected the fragmentation energies calculated for the MoDTC compounds, and revealed that such a mechanism strongly depends on the presence of the metallic substrate. A further step of the analysis, which will be described in a forthcoming article, is to consider also passivated metallic surfaces and the effect of oxygen in the substrate to better describe the first tribochemical events, when the native surface is yet to be reached.
The QM/MM approach was proven to be extremely valuable for tribochemistry, as it provided the necessary insight to clarify the debated dissociation mechanism of MoDTC.
Such an approach appears to be very promising for the design of new materials to reduce friction, where the chemistry of the species active in a tribological environment must be accurately taken into account.
Preparation of the QM/MM systems; the cutoff for the kinetic energy of the wave functions; and details on the dissociation of MoDTC complexes on iron (PDF)