SLIDE 1 Report
Francesco Marino
Contents
1 Introduction 1 2 Theoretical backgroung 2 2.1 Quasi-elastic neutron scattering . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 Molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Simulation details 6 4 Results and discussion 7 5 Summary and conclusion 9 References 9
Abstract This report details the work I have undertaken at Institut Laue Langevin (ILL) as a participant of the XRay and Neutron Science Summer Programme 2017 over the period September 4 - 29. I worked under the supervision of Peter Fouquet, helping him set up a molecu- lar dynamics (MD) simulation of the diffusion of ammonia molecules on a graphite substrate. MD results show that adsorbate molecules follow a very fast jump diffusion mech- anism, characterized by a diffusion coefficient of the order of 10−8m2/s, in agreement with time-of-flight neutron spectroscopy measurements. These encouraging preliminary studies suggest tha MD can be a valid tool for future investigation of surface diffusion.
1 Introduction
The diffusion of ammonia on graphite is particularly interesting for potential applications
- f graphene and graphitic material surfaces. These include the use of graphene as sensor
for ammonia groups and the n-doping of graphene by means of ammonia gas. 1
SLIDE 2
In this project, the properties of the adsorption of ammonia on an exfoliated graphite substrate were studied combining experimental and computational methods [2]. Experi- ments were carried out with neutron time-of-flight (TOF) technique at IN6 beamline at Institut Laue Langevin (ILL). Quasi-elastic neutron scattering (QENS) is well suited for studying fast diffusion processes as it can resolve the motion of atoms and molecules on a timescale from ns to ps. Furthermore, models used for describing different types of diffusion make distinct predictions on the properties of the dynamical scattering function, namely about the shape and momentum transfer dependence of the quasi-elastic broadening; this makes it possible to compare theory and QENS data. As for the computational side, we set up a molecular dynamics (MD) simulation. MD allows to follow the evolution of the system in time and is therefore useful to have a more intuitive feeling of the mechanism that rules the diffusion process. Moreover, since MD is an ab initio method, no additional assumptions are needed: the observed behaviour is the consequence of the interaction between atomic constituents. COMPASS II was chosen as a force field for our simulations. Whilst it has been extensively validated for bulk studies, there is few evidence concerning its accuracy as far as surface processes are concerned. Thus, one further motivation for our work on MD is carrying a preliminary test on whether this kind of simulations may be reliable tools for studying adsorption. Experimental and molecular dynamics results show a good agreement, at least in the general features of the adsorption mechanism. In fact, they both suggest that ammonia molecules follow a very fast jump motion.
2 Theoretical backgroung
2.1 Quasi-elastic neutron scattering
The time-of-flight (TOF) spectroscopy technique is used to measure the dynamic structure factor, or dynamic scattering function, S(Q, ∆E) of neutrons scattered by a sample. ∆E = Ef −Ei and Q = kf −ki are the energy transfer and the momentum transfer, respectively. The basic experimental procedure is the following (Figure 1). A beam of neutrons is sent through a chopper or a velocity selector. Only neutrons with a defined velocity are allowed to pass: as a result, a monochromatic beam of known wavelength and kinetic energy is produced. The selected neutrons are then directed against the sample, where they are scattered by nuclei. The interaction leads to a variation in the energy as well as the momentum of probes. Scattered neutrons are collected by the detectors, which cover a wide angle around the sample. The time interval between scattering and detection is recorded. As the flight distance is known, this allows to calculate the kinetic energy of scattered neutrons. Thus, the energy transfer ∆E can be easily obtained. 2
SLIDE 3 The momentum transfer Q, instead, is related to the scattering angle by simple formula: Q = 4π
λ sin(θ), where λ is the wavelength and θ is the angle between the incident and the
scattered beams. Figure 1: Schematic representation of IN6 beamline at ILL, where TOF measurements have been performed In our experiment, TOF was used to carry out quasi-elastic neutron scattering (QENS)
- measurements. QENS [3] is characterized by small energy transfers (of the order of neV to
meV ) between neutrons and target particles, leading to a broadening of the elastic peak in the energy distribution. This so called quasi-elastic broadening is the signature of diffusive processes. Theoretical models describing diffusion mechanisms make prediction on the quasi-elastic broadening, namely about the shape of the quasi-elastic profile and the dependence of its width on the momentum transfer . We will refer to the FWHM of the energy distribution as Γ(Q). Generally, there are three standard cases for diffusion: ballistic diffusion, continuous diffusion and jump diffusion. Ballistic diffusion describes processes when there is no sig- nificant friction between adsorbate and substrate particles. The shape of the quasi-elastic profile as a function of energy transfer is a Gaussian curve. The quasi-elastic broadening is proportional to the magnitude of the momentum transfer: Γ(Q) = 2
m Q (1) Continuous diffusion is described by the Brownian model. Adsorbate particles are 3
SLIDE 4 assumed to follow a random walk, i.e. molecules move uniformely and, when two of them collide, their velocities change randomly. Brownian motion is charaterized by a high friction coefficient between substrate and adsorbates. The quasi elastic profile, in this case, is a Lorentzian distribution, whose FWHM is proportional to Q2. Γ(Q) = 2DQ2 (2) Here, D is the diffusion coefficient, related to the temperature T and the friction coef- fcient η by the Einstein’s relation: D = kBT mη (3) Finally, jump diffusion model describes processes where there is little friction, but activation energy is high. Adsorbate molecules follow a stepwise motion: they spend on average a time τ on an equilibrium position, then make an istantaneous jump to another (random) site, as a consequence of thermal fluctuations. The quasi-elastic profile is still Lorentzian-shaped, but the dependence on the momentum transfer is different. In fact, if we call l the jump vector connecting two adiacent sites, the expression of the quasi-elastic broadening is a sinusoidal function of Q · l, namely: Γ(Q) = τ
sin2(l · Q 2 ) (4) Figure 2: 2D plot of the dynamic scattering function S(Q, δE) extracted from the neutron TOF data relative to a 0.9 ML ammonia coverage sample at 94 K. 4
SLIDE 5 2.2 Molecular dynamics
We briefly present the fundamental ideas of molecular dynamics. For a more detailed introduction to the subject, see [5] and [4]. For COMPASS force field, read [1], [6] and [7]. Molecular dynamics (MD) is a simulation method in which a system made of atoms and molecules is allowed to evolve from an initial configuration according to the laws of classical physics. Positions and velocities of the constituents of the system at each time step are calculated numerically solving Newton’s equation: mi ¨ ri = −∇iV (r1, ..., rN) (5) Here, ri denotes the position vector of the i-th particle of the system, V is the potential energy and ∇i = ( ∂
∂xi , ∂ ∂yi , ∂ ∂zi ).
The potential energy function has to be carefully designed in order to achieve a realistic
- utput. In COMPASS II force field, V is splitted in bonded, non bonded and cross- terms:
V = Vbond + Vnon−bond + Vcross (6) Bonded terms deal with strong chemical bonds (covalent, ionic). Bonds are modelled as springs (with additional degrees of freedom) connecting atoms. Vbond is the sum of polinomial (2nd to 4th order) and sinusoidal contributions used to describe stretch, bending and torsion of the bonds as perturbations around equilibrium. Non-bonded terms, instead, concern weaker inter-atomic pair interactions, i.e. electro- static and van der Waals interactions: Vnon−bond = VvdW + Velec (7) As for the former, the Coulomb potential is used: Velec =
qiqj rij (8) The sum runs over all pairs of atoms in the system. qi and qj denote the partial electric charges on the i-th and j-th atoms and rij is the distance between nuclei. Van der Waals forces are described by the so called Lennard-Jones pair potential, whose general expression is: VvdW =
ǫ(Aij rn
ij
− Bij r6
ij
) (n > 6) (9) While the term proportional to −r−6
ij
is motivated by theoretical calculations, the ex- ponent n and the values of the coefficients are empirically chosen. In COMPASS II, n = 9. 5
SLIDE 6 Figure 3: Plot of Lennard-Jones pair potential as a function of the distance between nuclei. Lennard-Jones potential captures the most relevant features of van der Waals interac- tions (Figure 3), which are repulsive at short distance and attractive at long distance. In particular VvdW goes asymptotically to zero for rij− > ∞, i.e. the force between atoms becomes negligible when they are far from each other. Also, the potentials has a global minimum which corresponds to the only equilibrium state of the two-particles system. In order to make simulations more accurate, force fields like COMPASS II include in the potential energy function further terms, describing the coupling between stretching, bending, and torsion. An example of cross term is: Vbond−bond = 1 2(r − r0)(r′ − r′
0)
(10) Parameters of the model are determined fitting a large data set. As experimental values are sometimes unreliable, input data are often obtained from ab initio quantum calculations.
3 Simulation details
The program package Materials Studio, and in particular the force field engine Forcite, was used throughout the project. Simulations were performed with a time step of 1 fs, writing the atomic positions to the trajectory files every 20 fs. The NVT (constant number of particles, volume and temperature) thermodynamic ensemble was used, together with a Nos´ e thermostat to keep the system at a steady temperature. The total simulation time 6
SLIDE 7 was set to 400 ps (400,000 steps). The system was allowed the first 50 ps to equilibrate, before using data for analysis. We carried tests with COMPASS II and Universal force fields. COMPASS II proved to be much more accurate than Universal, although more computationally expensive. There- fore, all the subsequent discussion refers to data gathered with COMPASS II. Atom types and charges were automatically assigned by the force field. The atom type h1n was assigned to the hydrogen atoms, the atom type n3* to the nitrogen atoms and the atom type c3 to the carbon atoms in the graphite, respectively. The van der Waals forces were calculated based on single atoms and truncated at a cut-off distance of 15.5 ˚ A by a cubic spline with a spline width of 1 ˚
- A. The atomic position before starting the simulations
were determined by converging the total energy with a convergence limit of 10−4 kcal/mol (4.3 × 10−3 meV). The model system was composed of a graphite substrate and of ammonia molecules. To describe high and low adsorbate coverage, we used 1 and 8 NH3 molecules respectively. As for the graphite substrate, a geometry of 5 × 5 unit cells on the surface plane and 2 unit cells on the vertical direction (c-axis) was chosen. Since the graphite unit cell is made
- f two parallel graphene cells, this setup yields a total of 4 sheets, the two inner of whom
were constrained. To avoid unwanted interactions between the ammonia molecules on the upper and lower sheets, the vacuum space above the graphene surface was increased to 35 ˚ A.
4 Results and discussion
The simulation film can give hints on the qualitative features of the dynamics of the system. In our case, the simulation output suggests that ammonia diffusion follows a jump mechanism. In fact, the general trend is made of molecules oscillating around an equilibrium position for a certain time interval, then suddenly jumping to another point. The adsorbate molecules’ motion is characterized by an unusual rolling pattern that may be worth further investigations. It is also interesting to point out that in our simulation NH3 diffuses not only on the substrate surface, but also in the vertical direction, although in the latter jumps are of limited length. To quantitatevely characterize the process, we calculated the diffusion coefficient D. For this purpose, we made use of the Einstein relation which relates it to the mean square displacement (MSD). MSD is defined as MSD(t) =< ∆r(t)2 >=< |r(t) − r(t0)|2 > (11) Here, r is the position vector at time t. In MD simulations, the average is performed
- n both time and system’s particles. In the long time range, a linear dependence of MSD
7
SLIDE 8 with respect of t is expected. In formulas: MSD(t) = 2fDt + C (12) f refers to the number of degrees of freedom of the system, while C is a constant. From the previous discussion follows that the diffusion coefficient is easily calculated once the coefficient of proportionality is obtained by means of a linear fitting of MSD(t). Table 1 shows values for both the constant of proportionality and the diffusion coeffi- cient for 8 different runs at 94 K at high ammonia coverage. MSD was calculated starting from 100 ps, for a duration of 100 ps. We assumed 3 degrees of freedom, because, as said before, adsorbate molecules show motion along the c-axis too. The limited number of data and the large discrepancies between measures make it difficult to work out an accurate statistical analysis. Nevertheless, the mean value may be used to have a rough estimate of the diffusion coefficient at 94 K, to be compared with the experimental evidence. The results are DMD = (5.7 ± 2.5)10−8m2/s and Dexp = (3.9 ± 0.4)10−8m2/s for the molecular dynamics and experimental data respectively. Although the uncertainty on the simulation-derived diffusion coefficient is very big, the orders of magnitude of DMD and Dexp coincide. This is a remarkable fact, as it suggests that molecular dynamics simulations lead to a diffusion mechanism which is, at least in its general features, similar to the one experimentally observed. In particular, both MD and experiment show that the adsorption
- f ammonia on graphite is a very fast process. Neutron spin echo measurements, also, do
not highlight any motion at longer time scales and thus further validate the previous conclusions. Coeff(10−8m2/s) D(10−8m2/s) 14.78 2.46 21.98 3.66 22.63 3.77 23.19 3.87 39.28 6.55 44.95 7.49 48.60 8.10 51.75 8.63 Table 1: The table shows coefficient of proportionality between MSD and t (Coeff) and the diffusion coefficient D relative to molecular dynamics simulations of the high coverage system at 94K. 8
SLIDE 9 5 Summary and conclusion
We have studied the adsorption of ammonia on exfoliated graphite using quasi-elastic scattering and molecular dynamics simulations. Time-of-flight measurements suggest that ammonia molecules follow a very fast jump motion. At 94 K, the 0.9 ML ammonia coverage sample is characterized by a diffusion coefficient of Dexp = (3.9±0.4)10−8m2/s. At the same temperature, MD gives for the high coverage system a value of DMD = (5.7±2.5)10−8m2/s. Molecular dynamics agrees with neutron scattering experiments on the general features
- f the motion, i.e. both suggest a fast jump mechanism with a high diffusion coefficient
- f the order of 10−8m2/s. Additional spin-echo measurements do not show any relevant
longer time-scale order, validating the fast dynamics hypothesis. In conclusion, MD simulations with COMPASS II force field proved quite accurate in reproducing the adsorption of ammonia on graphite. The next step after these encouraging preliminary results would be to make a full scale detailed simulation to further investigate the properties of this process.
References
[1] Accelrys, San Diego. Cerius2 4.8, Forcefield-Based Simulations, 2003. [2] I. Calvo-Almazan M. Zbiri M.M. Koza W.E. Ernst P. Fouquet Anton Tamtoegl, M. Sac-
- chi. Ultrafast molecular transport on carbon surfaces: The diffusion of ammonia on
- graphite. Carbon, 126:23–30, 2018.
[3] I. Calvo-Almaz` an and P. Fouquet. The application of quasi-elastic neutron scattering techniques (qens) in surface diffusion studies. Eur. Phys. J. Special Topics, 213, 2012. [4] M. A. Gonz´
- alez. Force fields and molecular dynamics simulations. Collection SFN,
12:169–200, 2011. [5] Ashlie Martini. Short course on molecular dynamics simulation, 2009. [6] H. et al. Sun. J Phys Chem B, 102:7338–7364, 1998. [7] Jin Z. Yang C. et al. Sun, H. . J Mol Model, pages 22–47, 2016. 9
