Reactive molecular dynamics simulations of organometallic compound W(CO)6 fragmentation

TL;DR

This study uses reactive molecular dynamics simulations to accurately model the fragmentation of W(CO)$_6$ under irradiation, reproducing experimental results and revealing mechanisms relevant to radiation chemistry.

Contribution

First detailed reactive MD simulation of W(CO)$_6$ fragmentation that matches experimental data and elucidates energy distribution effects on fragmentation pathways.

Findings

Energy distribution causes CO evaporation and W(CO)$_n^+$ formation.

Localized energy leads to C--O bond cleavage and WC(CO)$_n^+$ fragments.

Calculated fragment energies agree with experiments.

Abstract

Irradiation- and collision-induced fragmentation studies provide information about geometry, electronic properties and interactions between structural units of various molecular systems. Such knowledge brings insights into irradiation-driven chemistry of molecular systems which is exploited in different technological applications. An accurate atomistic-level simulation of irradiation-driven chemistry requires reliable models of molecular fragmentation which can be verified against mass spectrometry experiments. In this work fragmentation of a tungsten hexacarbonyl, W(CO)$_6$, molecule is studied by means of reactive molecular dynamics simulations. The quantitatively correct fragmentation picture including different fragmentation channels is reproduced. We show that distribution of the deposited energy over all degrees of freedom of the parent molecule leads to thermal evaporation of CO groups and the formation of W(CO)$_n^+$ ($n = 0-5$) fragments. Another type of fragments, WC(CO)$_n^+$ ($n = 0-4$), is produced due to cleavage of a C--O bond as a result of the localized energy deposition. Calculated fragment appearance energies are in good agreement with experimental data. These fragmentation mechanisms have a general physical nature and should take place in radiation-induced fragmentation of different molecular and biomolecular systems.

Tables

Table 1. Table 1: Parameters of the reactive force field used in the simulations. In the case of WC(CO) + n superscript subscript absent 𝑛 {}_{n}^{+} , W–C bonds labeled as (1) and (2) refer to the notations of Fig.  1 (c).

	W(CO)${}_{}{}^{+}{}_n{}^{}$		WC(CO)${}_{}{}^{+}{}_n{}^{}$
bond type	W–C	C–O	W–C (1)	W–C (2)	C–O
$r_0$ (Å)	2.11	1.14	1.83	2.31	1.14
$D$ (kcal/mol)	40.7	212.8	143.3	22.6	212.8
$k^{\mathrm{bond}}$(kcal/mol/Å2)	119.9	1493.9	369.4	67.7	1493.9

Equations

$$U^{\mathrm{bond}}(r_{ij})=k_{ij}^{\mathrm{bond}}(r_{ij}-r_{0,ij}{)}^2\text{mbox},$$

$$U^{\mathrm{bond}}(r_{ij})=D_{ij}\left[e^{-2{\beta }_{ij}(r_{ij}-r_{0,ij})}-2e^{-{\beta }_{ij}(r_{ij}-r_{0,ij})}\right]\text{mbox}.$$

$${\vec{v}}_1=\frac{\sqrt{2\mu E_{\mathrm{dep}}}}{m_1}\vec{u},{\vec{v}}_2=-\frac{\sqrt{2\mu E_{\mathrm{dep}}}}{m_2}\vec{u}.$$

$$\sum _i^N\frac{1}{2}{m}_i(\alpha v_i^{\mathrm{eq}}{)}^2=\frac{3N{k}_{\mathrm{B}}T}{2}+E_{\mathrm{dep}},$$

Full text

1 Introduction

Different radiation sources (electron and ion beams, lasers, synchrotrons) are commonly used to explore structure and dynamics of complex molecular systems ISACC_LatestAdvances_2008 ; RadDamage_Gustavo ; SynchrRad_StructBiol . Upon interaction with radiation a molecular target becomes excited or ionized and may relax via different dissociation mechanisms Schlatholter_2005_NIMB.233.62 ; Spezia_2016_FarDisc.195.599 . The analysis of fragments produced in a collision and their kinematic properties provide useful information on the parent molecule, e.g., its geometry and electronic structure, the type of interaction between structural units of the system, etc. Poully_2015_PCCP.17.7172 . Thus, it is important to understand how the energy deposited into the target is distributed among its degrees of freedom and how this distribution leads to the observed fragmentation patterns.

Irradiation of molecular targets and subsequent irradiation-driven chemistry (IDC) are exploited in modern and emerging technologies. For instance, they play an important role in modern radiotherapies such as ion-beam cancer therapy RadDamage_Gustavo ; Nano-IBCT_book that exploits the ability of heavy charged particles to inactivate living cells due to induction of complex DNA damage. Radiation is utilized also for controlling various physical and chemical processes for specific technological needs, e.g. for the formation and controllable growth of metal nanostructures under the exposure to focused electron or ion beams Utke_2008_JVacSciTechnolB.26.1197 ; Utke_book , or for the production of thin films with tailored structural properties.

One of the technological applications of IDC is focused electron beam induced deposition (FEBID) – a novel and actively developing nanofabrication technique that allows controllable creation of metal nanostructures with nanometer resolution Huth_2018_review ; deTeresa_2016_JPD.49.243003 ; Utke_book . FEBID is based on the irradiation of organometallic precursor molecules by keV electron beams whilst they are being deposited on a substrate. Electron-induced decomposition of the molecules releases its metallic component which forms a deposit on the surface with size similar to that of the incident electron beam (typically a few nanometers).

To date, a popular class of precursors for FEBID is metal carbonyls Men(CO)m Huth_2012_BeilsteinJNanotech ; Kumar_2018_Beilstein.9.555 , which are composed of one or several metal atoms (Me) bound to a number of carbon monoxide ligands. Although such precursor molecules (for instance, W(CO)6, Fe(CO)5 or Co2(CO)8) were commonly adopted from the chemical vapor deposition method and thus were not specifically designed to be efficiently and completely dissociated under electron irradiation, recent developments in optimization of the precursor deposition processes as well as the design of novel precursors for FEBID (e.g., HFeCo3(CO)12 Porrati_2015_Nanotechnology.26.475701 ) have made it possible to fabricate fully metallic structures made of Au, Pt, Fe, Co, Pb, and Co3Fe alloy Huth_2018_review .

Metal carbonyls have been widely studied experimentally and a substantial amount of data on thermal decomposition and electron-induced fragmentation has been accumulated over the past decades Beranova_1994_JAMS.5.1093 ; Clements_1976_MetTransB.7.171 ; Cooks_1990_JASMS.1.16 ; Wysocki_1987_IJMSIP.75.181 . This is due to their peculiar structure containing strong C–O bonds and relatively weak Me–C bonds. While the former are difficult to cleave, the latter dissociate easily, what usually happens through a sequential loss of CO groups when sufficient internal energy is available. This allows to determine internal energy distributions after electron-impact ionization, which are commonly used in mass spectrometry. Due to these properties, W(CO)6 and other metal carbonyls are commonly referred to as “thermometer molecules” Susic_1992_JMassSpectrom.27.769 ; Wysocki_1987_IJMSIP.75.181 .

In recent experimental studies of W(CO)6 and Fe(CO)5 molecules Wnorowski_2012_RapCommunMassSpectr.26.2093 ; Allan_2018_PCCP.20.11692 ; Wnorowski_2012_IJMS.314.42 ; Neustetter_2016_JCP.145.054301 ; Lacko_2015_EPJD.69.84 ; Lengyel_2016_JPCC.120.17810 fragmentation mass-spectra upon electron impact were recorded as well as appearance energies for the formation of different molecular fragments were evaluated. These studies considered both low-energy electron collisions (with the typical incident energy of several eV) where fragmentation is deemed to occur via dissociative electron attachment (DEA) Wnorowski_2012_RapCommunMassSpectr.26.2093 ; Allan_2018_PCCP.20.11692 and collisions with more energetic electrons (with the energy of one to several tens of eV) which decompose the parent molecules via dissociative ionization (DI) Wnorowski_2012_IJMS.314.42 ; Neustetter_2016_JCP.145.054301 ; Lacko_2015_EPJD.69.84 ; Lengyel_2016_JPCC.120.17810 .

DI experiments with W(CO)6 Wnorowski_2012_IJMS.314.42 revealed that appearance energy for W(CO)${}_{5}^{+}$ (i.e. loss of one CO group) is about 10 eV while the appearance energy for W+ (loss of all ligands) is about 20 eV. The intermediate fragments, W(CO)${}_{n}^{+}$ ($n=1-4$), were also detected within this energy range. Another type of fragments, WC(CO)${}_{n}^{+}$ ($n=0-3$), appears above 20 eV. Furthermore, the doubly charged fragments, W(CO)${}_{n}^{2+}$ and WC(CO)${}_{n}^{2+}$, were detected above 40 eV, but their abundance is about an order of magnitude smaller than for the respective singly-charged species Wnorowski_2012_IJMS.314.42 .

A detailed atomistic-level understanding of IDC processes (i.e., bond dissociation and further reactivity) in molecular systems is achievable by means of computational modelling. A rigorous quantum-mechanical description of these processes is possible only for relatively small systems containing, at most, a few hundred atoms and evolving on the sub-picosecond time scale. Classical molecular dynamics (MD) is considered as an alternative modelling framework for much larger systems and time scales. However, standard classical MD is unable to simulate the IDC processes as it typically does not account for coupling of the system to the incident radiation nor does it describe the induced quantum transformations.

These deficiencies were overcome recently by means of novel methodologies – reactive CHARMM force field Sushko2015 and Irradiation Driven Molecular Dynamics (IDMD) Sushko2016 ; MBNbook_Springer2017 – allowing, for the first time, high accuracy simulation of IDC in complex molecular systems. These methods were implemented into MBN Explorer Solovyov2012 – a multi-purpose software package for advanced multiscale simulations of structure and dynamics of complex molecular systems with the sizes ranging from the atomic up to the mesoscopic scales MBNbook_Springer2017 . Within the IDMD framework the IDC transformations are treated as random, fast and local processes which can be incorporated locally into classical MD force fields according to the probabilities of the quantum processes that may occur in the system. Such random events result in electronic excitations involving molecular orbitals localized on a specific bond or in a specific part of a molecule. This eventually leads to the cleavage of particular bonds and creation of active species which can undergo further reactivity. In Ref. Sushko2016 IDMD was successfully applied for the simulation of the FEBID process of W(CO)6 precursors. It was demonstrated that IDMD is capable of reproducing experimental observations and making predictions about the morphology, molecular composition and growth rate of tungsten nanostructures emerging on a surface during FEBID Sushko2016 . These methodologies have also been applied recently to study IDC in biomolecular systems in relation to ion-beam cancer therapy deVera_2018_EPJD.72.147 .

An accurate atomistic-level simulation of IDC calls for reliable models of radiation-induced fragmentation of molecular systems. With respect to the FEBID process, electron-beam induced fragmentation and reactivity on a substrate is a complex problem that combines many interlinked phenomena, such as the interaction of precursor molecules with a substrate, primary beam transport and electron multiple scattering on the substrate, ejection of backscattered and secondary electrons of different energies, their interaction with precursors by different mechanisms (e.g., DEA or DI) leading to different fragmentation pathways, reactivity of the fragmentation products among themselves and with other molecules, etc. Mass spectrometry experiments with electron beams represent a convenient means to test and validate the molecular fragmentation models, since the abundance of fragments of different masses and charges can be accurately measured under single electron–molecule collision conditions.

The goal of this work is to develop and validate a model for molecular fragmentation after energy deposition that can be used in IDMD simulations of FEBID Sushko2016 and in other situations where IDC is relevant, e.g., for studying radiation-induced biodamage Nano-IBCT_book ; deVera_2018_EPJD.72.147 . As an illustrative case study we analyze fragmentation of a tungsten hexacarbonyl, W(CO)6, molecule for which mass spectrometry data on electron impact ionization is abundant Wnorowski_2012_IJMS.314.42 . It is demonstrated that the quantitatively correct fragmentation picture including different fragmentation channels can be reproduced by means of classical MD simulations. We show that loss of CO ligands is mainly the result of the thermal evaporation process wherein excess energy is distributed among all degrees of freedom of the parent molecule. This mechanism cannot describe the formation of WC(CO)${}_{n}^{+}$ fragments, which are produced due to the initial cleavage of a C–O bond. This process occurs due to localized energy deposition in the particular bond. Calculated appearance energies for the production of different fragments are in good agreement with experimental data Wnorowski_2012_IJMS.314.42 .

2 Computational methodology

MD simulations were performed in this work using MBN Explorer Solovyov2012 – a software package for advanced multiscale modeling of complex molecular structure and dynamics. Among many options for molecular dynamics available in the software, it includes the reactive CHARMM force field Sushko2015 which was utilized in this work as well as the IDMD module Sushko2016 which will benefit from the fragmentation model developed in this study. The program counts on with the MBN Studio graphical interface MBNStudio_paper_2019 which was used to construct the molecular system, perform the simulations and analyze the results. In the following sections parameters of the reactive force field and the model of W(CO)6 fragmentation are described in detail.

2.1 Reactive CHARMM force field

The reactive CHARMM force field Sushko2015 permits classical MD simulations of a large variety of molecular systems experiencing chemical transformations whilst monitoring their molecular composition and topology changes MBNbook_Springer2017 ; Sushko2016 ; Sushko2015 . It can be applied for studying processes where rupture of chemical bonds plays an essential role, e.g., in irradiation- or collision-induced damage Verkhovtsev_2017_EPJD.71.212 ; deVera_2018_EPJD.72.147 .

This reactive force field represents a modification of the standard CHARMM force field MacKerell1998 , which employs the harmonic approximation for the description of interatomic interactions thereby limiting its applicability to small deformations of molecular systems. Contrary to the standard CHARMM force field, its reactive modification goes beyond the harmonic approximation for modeling covalent bonds, angles and dihedral angles. This allows for a more accurate description of the physics of molecular dissociation. In our recent works the reactive CHARMM force field was utilized to study ion-induced radiochemistry deVera_2018_EPJD.72.147 and thermal splitting Sushko2015 of water as well as collision-induced multi-fragmentation of C60 fullerenes Verkhovtsev_2017_EPJD.71.212 . Here it is used to study fragmentation of organometallic compounds.

In the standard CHARMM force field MacKerell1998 all bonded interactions are described by harmonic potentials. For example, the interaction between atoms $i$ and $j$ forming a covalent bond is given by:

[TABLE]

where $k_{ij}^{\rm bond}$ is the force constant of the bond, $r_{ij}$ is interatomic distance and $r_{0,ij}$ is the equilibrium bond length. Angular and dihedral interactions involving triples and quadruples of atoms are defined in a similar manner MacKerell1998 .

In the reactive force field Sushko2015 , the harmonic bond interaction described by Eq. (1) is replaced by a Morse potential:

[TABLE]

Here $D_{ij}$ is dissociation energy for the bond between atoms $i$ and $j$, and the parameter $\beta_{ij}=\sqrt{k_{ij}^{\rm bond}/D_{ij}}$ determines steepness of the potential. An additional parameter can be defined, the cutoff distance, characterizing the distance beyond which the bond is considered as broken and the molecular topology of the system changes. The bond energy calculated by Eq. (2) asymptotically approaches zero at large interatomic distances. Switching functions are defined accordingly to gradually reduce the angular and dihedral interactions as the bond breaks Sushko2015 .

2.2 Determination of the force field parameters

Parameters of the reactive force field for W(CO)6 and its fragments were determined from DFT calculations which were benchmarked against experimental data Wnorowski_2012_IJMS.314.42 . The DFT calculations were performed using Gaussian 09 software Gaussian09 employing the B3LYP exchange-correlation functional and a mixed LanL2DZ/6-31+G(d,p) basis set, wherein the former set described the W atom and the latter was applied to C and O atoms. As we aim at reproducing appearance energies from electron-impact ionization experiments, singly charged parent molecule and molecular fragments were considered. Geometry of each molecule was optimized first (states with spin multiplicity $M=2,4$ and 6 were considered) and a potential energy surface scan was then performed for different W–C and C–O bonds allowing to calculate equilibrium bond lengths, dissociation energies and force constants. Atomic partial charges were obtained through the natural bond orbital analysis Gaussian09 .

To benchmark the methodology, we analyzed equilibrium bond lengths in a neutral W(CO)6 molecule. The calculated values, $r_{0,{\rm{W-C}}}=2.07$ Å, $r_{0,\rm{C-O}}=1.15$ Å, are in good agreement with experimental data Arnesen_1966_ActaChemScand.20.2711 and with the results of earlier DFT calculations Szilagyi_1997_Organometallics.16.4807 . The calculated dissociation energies of W–C bond in the parent cation W(CO)${}_{6}^{+}$ and in different fragments such as W(CO)${}_{n}^{+}$ ($n=1-5$) and WC(CO)${}_{n}^{+}$ ($n=0-5$) are in good agreement with the results of Ref. Wnorowski_2012_IJMS.314.42 with the relative discrepancy of a few kilocalories per mole.

Figure 1(a) shows potential energy curves for W–C and C–O bonds in the parent W(CO)${}_{6}^{+}$ molecule as well as in WC(CO)${}_{5}^{+}$ formed upon removal of an oxygen atom from one of the ligands. Optimized structures of these molecules are shown in panels (b) and (c), respectively. Symbols illustrate the results of DFT calculations while lines show a fit to this data with the Morse potential, Eq. (2). The W–C bond in the parent molecule is significantly weaker ($D=40.7$ kcal/mol) than the C–O bond ($D=212.8$ kcal/mol), which is a common feature of metal carbonyls Diefenbach_2000_JACS.122.6449 . In the case of WC(CO)${}_{5}^{+}$, three different types of W–C bonds, labeled as (1), (2) and (3), can be distinguished, see Fig. 1(c). When an oxygen atom is removed, the remaining carbon atom of the ligand becomes stronger bound to the tungsten atom (bond (1)); the dissociation energy of this W–C bond increases and varies from 119.9 to 143.3 kcal/mol depending on the WC(CO)${}_{n}^{+}$ fragment considered. The opposite CO group becomes weakly bound to the metal atom (W–C bond (2)) and the equilibrium distance between W and C increases from 2.11 Å to 2.31 Å. Four W–C bonds in the perpendicular direction (bond (3)) remain unaffected and their dissociation energy does not change with respect to that in the parent molecule.

Table 1 summarizes the parameters of the reactive force field used in the calculations. Since the dissociation energy of W–C bond in different W(CO)${}_{n}^{+}$ ($n=1-5$) fragments does not vary significantly Wnorowski_2012_IJMS.314.42 , the value for the formation of W(CO)${}_{5}^{+}$, $D=40.7$ kcal/mol, was used in the simulations.

2.3 Model for energy deposition

Molecular fragmentation as a result of energy deposition involves several stages that take place on different time scales. First, the incident radiation (e.g., an electron) interacts with the molecule and transfers energy to it by means of different mechanisms, e.g., electronic excitation, ionization or DEA. These are fast processes that happen on the sub-femtosecond scale and leave the molecule in an excited electronic state. In the case of ionization, some fraction of deposited energy is spent in overcoming the ionization threshold, another fraction is carried away by the ejected electron, while the remaining part is stored in the target in the form of electronic excitations. The latter can involve different molecular orbitals, being of either bonding or antibonding nature.

An excitation involving an antibonding molecular orbital will quickly evolve through cleavage of a particular bond on the femtosecond timescale. The fragmented parent molecule may still keep some amount of the deposited energy which can lead to the sequential fragmentation of other bonds on the picosecond or even longer timescales.

The excited electronic state may also involve a bonding molecular orbital or may not be localized on a particular bond. In this case, the excess energy can be redistributed over a larger part or even the entire volume of the system and be transferred later into its vibrational degrees of freedom. Relaxation of the deposited energy due to electron-phonon coupling mechanism Gerchikov_2000_JPB.33.4905 leads to an increase in the amplitude of thermal vibrations which, in turn, leads to evaporation of loosely bound CO ligands. As it was shown in the case of small metal clusters Gerchikov_2000_JPB.33.4905 , the electron-phonon coupling is a slow process (as compared to a characteristic time of electron-molecule interaction) that happens on a picosecond time scale. The subsequent evaporation process may last up to microseconds.

In this work, we focus on the events happening from the cleavage of a particular bond to the redistribution of excess energy over the internal degrees of freedom of the molecule, i.e., on the time scales spanning from pico- to microseconds. Although quantum mechanical calculations can simulate the process of energy deposition, excitation of specific electronic orbitals and the initial stage of the fragmentation process, the span of the entire fragmentation process goes far beyond the limits of quantum MD. Thus, classical MD remains the only computational technique allowing exploration of the process within the required time frame at the atomistic level of detail. Quantum mechanical calculations would be useful though to establish the initial conditions for the multiscale simulation of molecular fragmentation, but this analysis goes beyond the scope of the present paper.

Within the framework of classical MD, we propose the following approach to describe the aforementioned fragmentation stages, i.e., fast cleavage of an individual bond (referred hereafter as stage I) and slow energy redistribution over all the molecular degrees of freedom (stage II from now on). Both processes result in an increase of internal energy of the molecule after the energy deposition, which is treated as an initial increase of the kinetic energy of atoms.

Within the approach we have used to model stage I the energy is deposited locally into a specific covalent bond of the target and converted into kinetic energy of the two atoms forming the bond (see Fig. 2(a)). Velocities of these atoms are defined to obey the total energy and momentum conservation laws:

[TABLE]

Here $E_{\rm dep}$ is the amount of deposited energy remaining in the system after ionization (i.e., excess energy over the first ionization potential), $m_{1}$, $m_{2}$ and $\mu=m_{1}m_{2}/(m_{1}+m_{2})$ are, respectively, masses and the reduced mass of the atoms forming the bond, and $\vec{u}$ is a unit vector defining the direction of the relative velocity of these atoms upon bond cleavage. The orientation of $\vec{u}$ may be determined by the field resulting from the local molecular configuration around the bond. In the current model, it has been chosen randomly.

Stage II is governed by the thermal mechanism of fragmentation where the energy is distributed over all degrees of freedom of the target. In this case, equilibrium velocities of atoms corresponding to a given temperature, $v_{i}^{\rm eq}$, are scaled111Note that only the absolute values of the velocities are scaled, while their directions are unaltered; thus, the total momentum is conserved. by a factor $\alpha$ depending on the amount of energy deposited (see Fig. 2(b)). The kinetic energy of $N$ atoms is then given by:

[TABLE]

where the first term on the right hand side corresponds to the initial kinetic energy of the atoms at equilibrium (e.g., $T=300$ K in our simulations), with $k_{\rm B}$ being the Boltzmann’s constant, while the second term is the excess energy deposited in the molecule during the collision.

As it will be discussed in Section 3, each of the above-described mechanisms leads to the formation of a particular group of experimentally observed fragments, while the whole experimental picture is reproduced well when both stages of the fragmentation process are considered. One should note that these stages happen sequentially in the same excited molecule, so that the initial cleavage of a specific bond may be followed by redistribution of the remaining energy over the molecular fragment.

2.4 Simulation details

The MD simulations of W(CO)6 fragmentation were performed using the MBN Explorer software Solovyov2012 . First, structure of the molecule was optimized using the parameters obtained from the DFT calculations discussed above. Then the molecule was equilibrated at $T=300$ K for 100 ns. The equilibration simulation was performed using the Langevin thermostat with damping time of 2 ps. Atomic coordinates and velocities were recorded every 100 ps. The equilibrated trajectory was sampled to generate random initial geometries and velocity distributions for the simulation of fragmentation. About 1000 constant-energy simulations each of 1 $\mu$s duration were conducted at different values of $E_{\rm dep}$ ranging from 0 to 475 kcal/mol. The upper limit is several times larger than the energy needed to break one W–C bond (see Fig. 1) which enables simulation of evaporation of several CO ligands. Fragments produced after 1 $\mu$s of simulation were analyzed, and the corresponding appearance energies were evaluated from this analysis and compared to experimental data.

3 Results and discussion

In the following, the fragmentation pathways of a W(CO)6 molecule upon electron impact ionization are discussed using the energy deposition model described in Section 2.3. Simulated appearance energies are compared to the most recent set of experimental data obtained from mass spectrometry for positively charged fragments produced by electron beams of energy $\leq 140$ eV Wnorowski_2012_IJMS.314.42 , as well as to the previously reported experimental values Michels_1980_InorgChem.19.479 ; Winters_1965_InorgChem.4.157 ; Bidinosti_1967_CanJChem.45.641 ; Foffani_1965_ZPhysChem.45.79 ; Qi_1997_JCP.107.10391 .

Figure 3 shows the number of CO fragments, $n$, produced upon dissociation of the parent W(CO)${}_{6}^{+}$ molecule for a given amount of excess energy $E_{\rm dep}$. Open blue squares show the results of simulations of stage I. In this case, a given amount of energy was deposited into a W–C bond resulting in a prompt release of one CO group. However, no further fragmentation has been observed even at high values of $E_{\rm dep}$. Due to a large difference in masses of a carbon and a tungsten atoms more than 90% of deposited energy is transferred into kinetic energy of the C atom and carried away by the released CO group. Therefore, independently of the value of $E_{\rm dep}$ given to the molecule through cleavage of a W–C bond, only a small amount of energy is transferred to the remaining W(CO)5 fragment, which is not sufficient to observe further fragmentation events.

As discussed in more detail further in this section, a single atom or a small fragment produced after cleavage of a specific bond can hit the remaining large fragment upon its escape, redepositing some amount of energy into the large fragment and thus triggering further fragmentation at stage II. Filled blue squares describe the situation when the CO group released due to cleavage of a W–C bond collided with the remaining W(CO)5 molecule, which led to the loss of another CO group. However, not more than two CO ligands have escaped the parent molecule in this case.

The results of simulations describing stage II are shown in Fig. 3 by filled red circles. Distribution of deposited energy over all degrees of freedom of the target leads to evaporation of multiple CO fragments. These results are in good agreement with the appearance energies reported in a recent experimental study by Wnorowski et al. Wnorowski_2012_IJMS.314.42 (black dashed lines) as well as with appearance energies obtained in earlier experiments Michels_1980_InorgChem.19.479 ; Winters_1965_InorgChem.4.157 ; Bidinosti_1967_CanJChem.45.641 ; Foffani_1965_ZPhysChem.45.79 ; Qi_1997_JCP.107.10391 (shaded areas). Note that the first ionization potential was subtracted from the experimental appearance energies to convert them into the excess deposited energy. Due to statistical nature of the fragmentation process, emission of a given number of CO groups was observed at different values of $E_{\rm dep}$. The lowest values at which five or less CO molecules were recorded after 1 $\mu$s agree nicely with the experimental results. However, complete fragmentation (i.e, loss of 6 CO molecules) has been observed only in a few trajectories at $E_{\rm dep}=375$ and 400 kcal/mol, which are significantly higher than the experimental appearance energy for W+. This indicates that the complete fragmentation takes place on a larger time scale and longer simulations are needed to observe the complete fragmentation at smaller energies.

The results of simulations describing stage II can be used to evaluate branching ratios for the production of different fragments for a given amount of excess energy $E_{\rm dep}$. This analysis provides information on how many carbonyl groups will be most likely evaporated at a given $E_{\rm dep}$ after 1 $\mu$s. These results are presented in Figure 4. It shows that emission of one, two and three CO groups ($n=1,2,3$) takes place in rather well separated energy “windows”. Fragments corresponding to $n=1$ and 2 were recorded in these energy ranges with the maximal probability corresponding to the branching ratio of 1. For the larger numbers of emitted CO the maximal branching ratios drop down to 0.8 suggesting an increased probability of observing different fragments. Note also that the characteristic energy ranges for the emission of $n$ CO groups increase with $n$. This analysis allows for evaluation of a typical amount of energy that should be deposited into the W(CO)6 molecule to observe a specific fragment. These values will be used in our future work as input for IDMD simulations in order to define the probabilities for different fragmentation pathways.

The experimental fragmentation mass-spectra Wnorowski_2012_IJMS.314.42 ; Michels_1980_InorgChem.19.479 ; Winters_1965_InorgChem.4.157 ; Bidinosti_1967_CanJChem.45.641 ; Foffani_1965_ZPhysChem.45.79 ; Qi_1997_JCP.107.10391 revealed the formation of not only W(CO)${}_{6-n}^{+}$ ($n=0-6$) but also WC(CO)${}_{5-n}^{+}$ ($n=2-5$) molecules, which, however, have not been observed in the simulations of stage II even after 1 $\mu$s of simulation. This is due to the very low probability for observing a statistical cleavage of a C–O bond, owing to the much lower strength of the W–C bond. We therefore assumed that a C–O bond can break after a localized energy deposition into it. As described in Section 2.2, loss of an oxygen atom from a CO group makes the opposite ligand weaker bound to the W atom (bond (2)), while the W–C bond corresponding to the cleaved C–O (bond (1)) becomes much stronger. Therefore, localized energy deposition into one C–O bond leads to a prompt release of an O atom, together with the recoil of the C atom to the parent molecule, which ends up being vibrationally excited and subsequently releases a CO fragment. This scenario explains also the production of smaller fragments due to subsequent loss of several CO groups from WC(CO)5 which were observed experimentally and in the simulations at $E_{\rm dep}>180$ kcal/mol, as shown in Fig. 5.

The fragmentation mechanism in this case is more complex than the simple evaporation of CO molecules, and thus deserves a more detailed analysis. An exemplary simulated trajectory for $E_{\rm dep}=375$ kcal/mol is shown in Fig. 6 (a movie is provided as supplementary material). Initial velocities of C and O atoms are shown by arrows in panel (a). As follows from the analysis of this trajectory, the initial cleavage of the C–O bond is a fast event happening within the first 50 fs of the simulation (panels (b-c)). It leads to the release of an energetic O atom as a result of energy redistribution according to Eq. (3) (stage I). This atom can directly escape from the parent molecule or collide with it in its way out, depending on the orientation of the vector $\vec{u}$. In the former case, less extensive fragmentation is observed as the O atom is ejected with a notable amount of kinetic energy. In the latter case, a fraction of its kinetic energy is redeposited into the parent molecule leading to further fragmentation (stage II).

The collision of the escaping O atom with the parent molecule is depicted in the next panels of Fig. 6. After about 50 fs (panel (c)) the ejected O atom collides with a neighbouring CO ligand. The latter is ejected from the parent molecule after approx. 130 fs (panel (d)) while the remaining C atom, becoming strongly bound to W, does not detach from the parent molecule (panel (e)). The remaining molecular fragment is still excited, and the evaporation process follows with the ejection of a CO group by 6 ps (panel (f)), another one by $\sim 15$ ps (panel (g)) and the last one by $\sim 230$ ps (panel (h)).

It should be noted that in this case the reactive force field is rather sensitive to the parameters used, particularly to the dissociation energies, and variation of these parameters can impact the results of simulations. In Fig. 5, black squares depict the results of simulations using the dissociation energies calculated by DFT (see Section 2.2), specifically $D_{\rm C-O}=212.8$ kcal/mol. While the appearance energy for WC(CO)4 is well reproduced,222Although this fragment was not observed experimentally, it appearance energy (gray dotted line in Fig. 5) can be estimated from the data reported in Ref. Wnorowski_2012_IJMS.314.42 . The calculated first appearance energy is smaller than the expected experimental value: this is due to the thermal energy of the molecule which also contributes to fragmentation together with the deposited energy. the rest are shifted towards larger energies. However, by analyzing the appearance energies reported by Wnorowski et al. Wnorowski_2012_IJMS.314.42 , one can estimate the dissociation energy of C–O bond of about 180 kcal/mol. It should be noted that the C–O bond in metal carbonyls is rather complex and can be characterized as a mixture of a triple and a double bond; $D_{\rm{C-O}}$ varies from 126 kcal/mol for a double bond to 257 kcal/mol for a triple bond Darwent1970 ; Dean1972 . Thus, values within this range are meaningful. Simulation results using $D_{\rm{C-O}}=180$ kcal/mol are shown in Fig. 5 by open red circles. The calculated appearance energies are shifted to lower energies with respect to the results employing the dissociation energy obtained from DFT, being closer to the experimental results.

Filled blue triangles show the results of simulations where the collision of the escaping O atom with the parent molecule has been taken into account. In this case, the first four appearance energies (WC(CO)4 to WC(CO)) are well reproduced, but the last one (WC) is not observed. This happens because the deposition of large amounts of $E_{\rm dep}$ larger than 390 kcal/mol leads to the cleavage of not only C–O bond but also the W–C bond labelled as (1) in Fig. 1(c). One may expect that this bond becomes stronger as the coordination number of W decreases (i.e. when smaller fragments are formed) that happens at large values of $E_{\rm dep}$. Thus, $D_{\rm{W-C}}$ could be larger than 143 kcal/mol, the value obtained from the DFT calculations for the WC(CO)${}_{5}^{+}$ molecule. We have therefore performed a set of simulations using $D_{\rm{W-C}}=180$ kcal/mol, that is similar to the dissociation energy of the C–O bond. These results are shown by filled green diamonds in Fig. 5.

Under these conditions, the experimental appearance energies are well reproduced. Although the reactive force field is rather sensitive to the parameters used, an appropriate choice of the parameters leads to a quantitative agreement with experiments. This shows that reactive classical MD simulations are appropriate to simulate the fragmentation patterns of W(CO)6. The method, being general, can be applied to other organometallic precursors for FEBID as well as to many other inorganic, organic and biological molecules.

4 Conclusion and outlook

We have developed and successfully validated a model for molecular fragmentation upon interaction with ionizing radiation. Two stages of the fragmentation process were considered, namely a localized energy deposition into a specific covalent bond and the redistribution of deposited energy among all the bonds of a target molecule. These situations represent limiting cases of the possible pathways for energy redistribution in a molecule, which involve the excitation to specific antibonding electronic orbitals or the energy transfer to vibrational modes of the molecule.

As a case study we analyzed fragmentation of a tungsten hexacarbonyl, W(CO)6, molecule which is widely used as a precursor for focused electron beam induced deposition (FEBID) and for which mass spectrometry data on electron impact ionization are abundant. We demonstrated that the quantitatively correct fragmentation picture including different fragmentation channels is reproduced well by means of classical MD simulations with reactive force fields.

The analysis revealed that loss of CO ligands and the formation of W(CO)${}_{n}^{+}$ ($n=0-5$) fragments is mainly the result of the thermal evaporation process where excess energy is distributed among all degrees of freedom of the parent molecule. Another type of fragments, WC(CO)${}_{n}^{+}$ ($n=0-4$), is produced mainly due to cleavage of a C–O bond as a result of the localized energy, followed by redistribution of excess energy in the parent molecule resulting in further evaporation of CO ligands. Calculated appearance energies for different fragments are in good agreement with experimental data Wnorowski_2012_IJMS.314.42 . We also analyzed the probability for the formation of different W(CO)${}_{n}^{+}$ fragments as a function of deposited energy and evaluated characteristic values of deposited energy needed to release a given number of CO ligands. This information will be used in our future work to model FEBID of metal carbonyl precursors and the electron-induced formation and growth of metal nanostructures within the Irradiation Driven Molecular Dynamics approach Sushko2016 .

The fragmentation model considered in this study has a clear physical explanation and describes the different stages of electron-induced fragmentation of a molecular target. This approach can therefore be utilized to study radiation-induced fragmentation of other FEBID precursors (other metal carbonyls in particular) as well as biomolecular systems whose fragmentation pathways are of great relevance for radiotherapy applications. Results of this analysis may bring insights into irradiation-driven chemistry of various molecular systems which is exploited in different modern and emerging technologies. Additionally, the model developed represents a first step towards the modelling of fragmentation mass spectra by means of classical MD simulations. This objective, which goes beyond the scope of the present paper and might be accomplished in subsequent works, would require the development of appropriate quantum mechanical models for predicting the probabilities of energy deposition after an inelastic collision, as well as the probability of energy transfer through the localized or thermal mechanisms.

Acknowledgements

PdV acknowledges the Alexander von Humboldt Foundation/Stiftung for its financial support by means of a postdoctoral fellowship. AV is grateful for the financial support from the DKFZ Postdoctoral Fellowship. This work was also supported in part by the Deutsche Forschungsgemeinschaft (Project no. 415716638), the Spanish Ministerio de Ciencia, Innovación y Universidades and the European Regional Development Fund (Project no. PGC2018-096788-B-I00) and by the COST Action CA17126 “Towards understanding and modelling intense electronic excitation” (TUMIEE). The possibility to perform computer simulations at Goethe-HLR and FUCHS clusters of the Frankfurt Center for Scientific Computing is gratefully acknowledged.