Molecular Dynamics Simulations of Hydrogen Diffusion in Aluminum

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Molecular Dynamics Simulations of Hydrogen Diffusion in Aluminum Xiaowang Zhou, Farid El Gabaly, Vitalie Stavila, and Mark D. Allendorf J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b01802 • Publication Date (Web): 23 Mar 2016 Downloaded from http://pubs.acs.org on March 27, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Dynamics Simulations of Hydrogen Diffusion in Aluminum X. W. Zhou*, F. El Gabaly, V. Stavila, and M. D. Allendorf

Sandia National Laboratories, Livermore, California 94550, USA ABSTRACT Hydrogen diffusion impacts the performance of solid-state hydrogen storage materials, and contributes to the embrittlement of structural materials under hydrogen-containing environments. In atomistic simulations, the diffusion energy barriers are usually calculated using molecular statics simulations where a nudged elastic band method is used to constrain a path connecting the two end points of an atomic jump. This approach requires prior knowledge of the “end points”. For alloy and defective systems, the number of possible atomic jumps with respect to local atomic configurations is tremendous. Even when these jumps can be exhaustively studied, it is still unclear how they can be combined to give an overall diffusion behaviour seen in experiments. Here we describe the use of molecular dynamics simulations to determine the overall diffusion energy barrier from the Arrhenius equation. This method does not require information about atomic jumps, and it has additional advantages, such as the ability to incorporate finite temperature effects and to determine the pre-exponential factor. As a test case for a generic method, we focus on hydrogen diffusion in bulk aluminium. We find that the challenge of this method is the statistical variation of the results. However, highly converged energy barriers can be achieved by an appropriate set of temperatures, output time intervals (for tracking hydrogen positions), and a long total simulation time. Our results help elucidate the inconsistencies of the experimental diffusion data published in literature. The robust approach developed here may also open up future molecular dynamics simulations to rapidly study *

Email: [email protected];

Tel.: 925-294-2851

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32

diffusion properties of complex material systems in multi-dimensional spaces involving composition and defects.

I.

INTRODUCTION As an important rate-limiting mechanism for kinetic processes in solids, mass diffusion is

often the barrier to the development of many new material technologies. For example, replacing fossil fuels with carbon-free hydrogen for transportation energy is widely explored. Some solidstate hydrogen storage materials can provide sufficient hydrogen for vehicles1,2. However, many of the high-capacity materials suffer from impractical hydrogen charging and discharging rates, thought to be due in part to the slow diffusion kinetics2,3. To develop storage materials that can meet DOE target for refuelling time, the fundamental mechanisms and energy barriers of diffusion in hydrogen storage materials must be understood2. The challenge is that the hydrogen charging / recharging cycles are associated with significant changes in both compositions and structures (e.g., phase transformation, interface migration, and defect formation). This material complexity requires studying diffusion in evolving structures, which has inhibited developing new hydrogen storage materials4. In another example, structural materials under hydrogen-containing environments are needed to store and deliver hydrogen. Hydrogen embrittlement is the limiting factor for these materials5. Hydrogen diffusion to defects (e.g., grain boundaries, dislocations, crack tips, etc.), directly contributes to the deterioration of mechanical properties with time. Understanding hydrogen diffusion is equally important to develop new structural materials. Atomistic simulations such as molecular dynamics (MD) and molecular statics (MS) provide a theoretical means to study diffusion. While the accuracy of results depends on the


2

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interatomic potential used in the simulations, the advantages of these methods are that they do not require assumptions, enable complex atomic scale structures to be accurately replicated in the computational cell, and allow diffusion mechanisms (e.g., atom jumps) to be visualized. In the past, diffusion energy barriers were usually calculated using molecular statics (i.e., energy minimization) rather than molecular dynamics simulations. By using a nudged elastic band method6,7,8 to constrain a path connecting the two end points of an atomic jump, a single MS simulation can determine the entire minimum energy path and the associated energy barrier. This approach requires prior knowledge of the “end points” of atomic jumps. In addition, the number of possible atomic jumps with respect to local atomic configurations is tremendous for alloy and defective systems. Even if these diffusion paths can be exhaustively studied, it is still unclear how they can be combined to give an overall diffusion behaviour seen in experiments. Alternatively, molecular dynamics simulations can be used to track the mean square displacement of diffusion atoms over a given period of time at different temperatures. The overall diffusion energy barrier can then be determined from the Arrhenius equation regardless of the number of the atomic jumps. This method has additional advantages such as incorporating finite temperature effects and determining the pre-exponential factor. Unfortunately, this method often has unsatisfactorily large statistical errors, and hence has not been widely used in the past. Using the diffusion of a single hydrogen atom in bulk aluminium as a test case, the objective of the present work is to both develop molecular dynamics models capable of accurately quantifying diffusion energy barriers from the Arrhenius equation, and understand hydrogen diffusion in aluminium. The MD code LAMMPS9,10 is used for this study. To provide researchers with an useful handbook reference missing in literature, we will examine various ensembles including NVT (constant number of atoms, volume, and temperature), NVE (constant


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32

number of atoms, volume, and energy), and NPT (constant number of atoms, pressure, and temperature). Both Nose-Hoover11 and Langevin12 algorithms are used to control temperature and pressure. Effects of various MD parameters (e.g., damping and dragging coefficients of thermostat and barostat, etc.) implemented in LAMMPS are all explored. We will validate the calculations using experimental data available in literature. The Al-H system is ideal for the study because the hydrogen diffusion in aluminium is known to be sensitive to composition and structure13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30 and it is this complexity that needs to be understood for hydrogen storage materials. Theoretical studies of this system also help understand the inconsistent experimental diffusion measurements13,14,15,16,17,18,19,20,21,22,23,24,25,26,27, 20,28

made on Al-H. Furthermore, aluminum is a light weight structural material suitable for

vehicle applications31. In particular, no hydrogen embrittlement has been reported for aluminum alloys when they are not used under fatigue or corrosive conditions32,33,34. This is in sharp contrast with steels35. As a result, understanding hydrogen diffusion in aluminum can have a high impact on a range of applications including use of this metal as a structural material. Our work will be based on an Al-H bond order potential (BOP)36,37, which we have shown to capture energy and volume trends for a variety of Al, H, and Al-H structures that are the determining factors for hydrogen diffusion. It also captures some important but elusive properties such as the high stacking fault energy of aluminium36, the energy barrier of the H + H2 ↔ H2 + H exchange reaction37, and the Al-H phase diagram36. II.

MOLECULAR DYNAMICS METHODS Our single crystalline, face-centred-cubic (fcc) aluminium cells contain 8 {100} planes in

each of the three coordinate directions. The initial crystals are created based on the room temperature experimental lattice constant a = 4.05 Å38. The system dimension is therefore around


4

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


32 Å3, corresponding to 2048 aluminium atoms. For comparison, we also calculate the theoretical lattice constants at finite temperatures following the published approach39, and find a = 4.05 Å at 300 K and a = 4.06 Å at 700 K. In comparison, the corresponding experimental values are a = 4.10 Å at 300 K and a = 4.15 Å at 700 K40 (note that the experimental values given by references 38 and 40 are different). Bulk crystals are simulated using periodic boundary conditions, and a single hydrogen atom is introduced in each simulated computational cell. The target temperature T is created by assigning velocities to atoms based on the Boltzmann distribution, which is further equilibrated by performing a warm-up MD simulation for more than 0.1 ns. After this, MD simulations of hydrogen diffusion are performed for a total period of tMD, which corresponds to a total time steps of nMD = tMD/dt where dt is the time step size. During MD simulations, the hydrogen location as a function of time, r(t), is tracked on a time interval of ∆t, i.e., at times of t = i∆t, i = 1, 2, …, m (m = tMD/∆t), where ∆t is a multiple of dt. The relative hydrogen displacement per ∆t, measured between (i-1)∆t and i∆t, is then calculated as ∆ri(∆t) = r(i∆t) - r(i∆t - ∆t), i = 1, 2, …, m. Once ∆r per ∆t is known, the relative displacement per larger time intervals of k∆t, measured between (i-1)∆t and (i-1+k)∆t, can be i −1+ k

simply obtained as ∆ri(k∆t) =

∑ ∆r (∆t ) . Clearly, for a given k∆t, we can compute m+1-k values j

j =i

of ∆ri(k∆t): i = 1, 2, …, m+1-k. As a result, the number of ∆ri(k∆t) values reduces to one when k approaches m, but there are many ∆ri(k∆t) values when k 24) seem to deviate from the linear function more significantly than the other points. This temperature effect can be understood because at low temperatures, the number of atom jumps is low resulting in statistical error due to insufficient samples. High temperatures may also pose a problem because an extremely high temperature is equivalent to a barrier-less


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

diffusion (i.e., jump frequency approaches attempt frequency). Then the simulation would be sampling the ballistic motion of the hydrogen atom rather than sensitively reflecting the thermally activated diffusion, especially considering the small hydrogen mass. Hence, it is always safer to apply an intermediate temperature range for calculating the diffusion energy barrier using MD. The temperature range we use for the ∆t = 4.4 ps simulations, T = 500 - 700 K (T/Q = 1220 – 1707 K/eV), appears to be a good choice from Fig. 2(b). Despite different values of ∆t, the statistical errors of both sets of data shown in Fig. 2(b) are too large to yield confident energy barrier calculations. We now explore the convergence of the activation energy as a function of total MD time tMD where the other parameters are held fixed at the better values determined above: dt = 0.001 ps, ∆t = 4.4 ps, and T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K (T/Q = 1220 – 1707 K/eV). The Arrhenius error ζ, and the corresponding activation energy Q, are calculated at different tMD, and the results are shown in Figs. 3(a) and 3(b) respectively. For comparison, the Arrhenius errors obtained from the data shown in Fig. 2(b), and the activation energy obtained from MS simulations, are indicated in Figs. 3(a) and 3(b) respectively using dash lines. Fig. 3(a) more clearly indicates that at tMD = 0.88 ns, increasing ∆t from 0.0088 to 4.4 ps greatly reduces the error of the Arrhenius fit. At the given ∆t = 4.4 ps, increasing tMD beyond 0.88 ns results in a systematic decrease of the Arrhenius error. In particular, the Arrhenius error becomes near zero (on the scale of the figure) when tMD reaches about 7 ns and above. Correspondingly, Fig. 3(b) indicates that the activation energy is extremely sensitive to tMD when tMD is small, but becomes insensitive to tMD at large tMD. In particular, the activation energy nearly saturates when tMD reaches 7 ns. Interestingly, the saturated activation energy shown in Fig. 3(b) is very close to the value obtained from MS simulations. Fundamentally, the MD data


10

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


should be taken as the standard if the statistical error is not large. This is because MD incorporates the temperature effect (also the mass effect as will be discussed below), and reflects the overall behaviour of all possible diffusion paths (e.g., O-T and T-O two-way jumps).

Fig. 3. (a) Arrhenius error ζ and (b) activation energy Q as a function of MD time tMD, obtained from MD simulations using the NVT ensemble at a time step size of dt = 0.001 ps, an output time interval of ∆t = 4.4 ps, and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K (T/Q = 1220 – 1707 K/eV). The discussions above indicate that a very low Arrhenius error can be achieved from MD simulations using the NVT ensemble at dt = 0.001 ps, ∆t = 4.4 ps, tMD = 13 ns, and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K. We now directly visualize the mean square displacement and the Arrhenius fits for such simulations in Fig. 4. Fig. 4(a) confirms that the mean square displacement has very small statistical errors and its median points satisfy very well the linear relationship with time for all the nine temperatures. More importantly, Fig. 4(b) indicates that the MD diffusivities satisfactorily fit the Arrhenius equation. This is an important result indicating that indeed MD can be a robust method to study diffusion, which is not limited by the complexity of the system.


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32

Fig. 4. (a) Mean square displacement as a function of time and (b) the Arrhenius plot for MD simulations conducted using the NVT ensemble at a time step size of dt = 0.001 ps, an output time interval of ∆t = 4.4 ps, a total MD time of tMD = 13.2 ns, and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K. C. Parameter Study Having established MD as a viable method for diffusion calculations, questions arise as how results depend on the many MD simulation parameters such as ensembles, thermostat methods, time step size, damping and dragging parameters for the thermostat and barostat, etc. To understand these, extensive MD simulations are performed to calculate diffusion activation energies at different ensembles and MD parameters. According to discussion above, all of our calculations use appropriately large values of ∆t (≥ 2.2 ps) and tMD (≥ 6.6 ns), and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K. By solving atom positions entirely from Newton’s equations, the NVE ensemble does not introduce additional parameters other than dt, ∆t, tMD and T. The NVT ensemble, based on either the Nose-Hoover11 or the Langevin12 thermostats, requires additional thermostat parameters9. For the Nose-Hoover thermostat, these include the temperature damping parameter (Tdamp), number of sub-cycles to perform thermostat within each time step (Tloop), and the drag factor (drag). For


12

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the Langevin thermostat, the parameters include Tdamp, the switch parameter for the GronbechJensen/Farago (gjf = yes or no) algorithm42,43, and the switch parameter to set the total random force to zero (zero = yes or no). The NPT ensemble is based on the Nose-Hoover barostat thermostat. In addition to Tdamp, Tloop, and drag mentioned above, it also involves pressure damping parameter (Pdamp), number of sub-cycles to perform barostat within each time step (Ploop), the number of steps for resetting the reference cell (nreset), and the switch parameter to use Martyna-Tuckerman- Klein correction44,45 (mtk = yes or no). For NVT or NPT ensembles, there are three additional conditions that can be specified: A, the thermostat is only applied to aluminium (i.e., leaving the hydrogen to the NVE ensemble) and the initial velocity of the hydrogen atom is set according to the simulated temperature T; B, the thermostat is only applied to aluminium and the hydrogen initial velocity is set to zero; and C, the thermostat is applied to both aluminium and hydrogen but the hydrogen initial velocity is set to T. Finally, the mass of the diffusion species can also be treated as a parameter so that we can explore the effects of mass (for example, effects of hydrogen isotopes). With the MD parameters defined, we now point out that the results discussed in Figs. 2-4 are obtained from an NVT ensemble with a Nose-Hoover thermostat applied only to aluminium (using condition B) at Tdamp = 100dt, Tloop = 1, and drag = 0.0. NVE does not have a thermostat or a barostat so that care is made to ensure correct NVE simulations. First, due to the equal partition of kinetic and potential energies, the initial temperature of the system is created at the value twice of the target temperature. After ignore the warm-up MD simulation, the time averaged temperature during the entire tMD simulation period is taken as the true temperature. The true temperatures slightly differ from the target temperatures (e.g., 485 vs. 500 K). We can of course repeat the simulation with the initial


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32

temperature scaled up correspondingly so that the true temperature matches exactly the target temperature. This, however, is not necessary because the analyses use only the true temperatures. Second, we do verify that the total energy is conserved extremely well (e.g., statistically changes from -6618.68 to -6618.50 eV, which can be attributed to numerical noises). 18 series of MD simulations are performed to examine effects of each MD parameter on activation energy barrier Q, pre-exponential factor D0, and the Arrhenius error ζ. In each series, one or several related MD parameters are varied whereas the other parameters are kept constant. Some parameters are changed together. For instance, when the number of steps are the same, increasing dt would cause a corresponding increase in ∆t and tMD. Under the conditions Tdamp = 100dt and Pdamp = 1000dt (which are normally used in MD simulations9), increasing dt would also cause a corresponding increase in Tdamp and Pdamp. The results of all of the 18 series of simulations are summarized in Table 1, where different ensembles are separated by double lines, and different series are separated by single lines and are marked by S#. We emphasize that the results in Table 1 are not trivial, involving nearly 2000 simulations. This is because each activation energy value is derived from nine temperatures, and each temperature consists of three simulations with different random number seeds to enable us to achieve the total tMD within 1/3 of the wall clock time. Table 1. Numerical values of pre-exponential factor D0, activation energy barrier Q, and error of the Arrhenius fit as a function of MD parameters dt, ∆t, tMD, Tdamp, Tloop, drag, mH, gjf, zero, Pdamp, Ploop, nreset, mtk, and “other”. Here other = A means that thermostat applies only to aluminum with hydrogen initial velocity not zero; other = B means that thermostat applies only to aluminum with hydrogen initial velocity zero; and other = C means that thermostat applies to all atoms with hydrogen initial velocity not zero. Note that in the table, different ensembles are separated by double lines, and different series are separated by single lines and are marked by S#. Here a series explores effects of targeted MD parameters while keeping the other parameters constant. S#

NVE parameters (mH=1amu) dt (ps) ∆t (ps)

tMD (ns)


D0 (Å2/ps)

Q (eV)

ζ

14

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

0.00075 3.3 0.00100 4.4 0.00125 5.5 NVT parameters (Nose-Hoover) S# dt (ps) tMD Tdamp Tloop drag ∆t (ps) (ps) (ns) 2 0.00050 2.2 6.6 0.050 1 0.0 0.00075 3.3 9.9 0.075 1 0.0 0.00100 4.4 13.2 0.100 1 0.0 0.00125 5.5 15.5 0.125 1 0.0 3 0.00075 3.3 9.9 0.075 1 0.0 0.00100 4.4 13.2 0.100 1 0.0 0.00125 5.5 15.5 0.125 1 0.0 4 0.00075 3.3 9.9 0.075 1 0.0 0.00100 4.4 13.2 0.100 1 0.0 0.00125 5.5 15.5 0.125 1 0.0 5 0.00075 3.3 9.9 0.010 1 0.0 0.00100 4.4 13.2 0.010 1 0.0 0.00125 5.5 15.5 0.010 1 0.0 6 0.00050 2.2 6.6 0.005 1 0.0 0.00050 2.2 6.6 0.010 1 0.0 0.00050 2.2 6.6 0.030 1 0.0 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.070 1 0.0 0.00050 2.2 6.6 0.090 1 0.0 0.00050 2.2 6.6 0.110 1 0.0 7 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 2 0.0 0.00050 2.2 6.6 0.050 3 0.0 0.00050 2.2 6.6 0.050 4 0.0 0.00050 2.2 6.6 0.050 5 0.0 8 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.2 0.00050 2.2 6.6 0.050 1 0.4 0.00050 2.2 6.6 0.050 1 0.6 0.00050 2.2 6.6 0.050 1 0.8 0.00050 2.2 6.6 0.050 1 1.0 0.00050 2.2 6.6 0.050 1 1.2 0.00050 2.2 6.6 0.050 1 1.4 9 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.0 0.00050 2.2 6.6 0.050 1 0.0

9.9 13.2 15.5 mH (amu) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 4 6 8 10

0.38 0.39 0.39 Q (eV)

0.0044 0.0124 0.0093 ζ

other*

73 81 83 D0 (Å2/ps)

A A A A B B B C C C A A A A A A A A A A A A A A A A A A A A A A A A A A A A A

203 161 149 107 160 180 121 146 137 138 147 80 68 92 225 151 203 121 140 371 203 60 114 70 389 203 87 220 188 79 179 101 64 203 103 94 26 24 26

0.44 0.42 0.42 0.40 0.42 0.43 0.41 0.42 0.42 0.42 0.42 0.39 0.38 0.39 0.44 0.42 0.44 0.41 0.42 0.47 0.44 0.37 0.41 0.38 0.47 0.44 0.40 0.44 0.43 0.39 0.43 0.40 0.37 0.44 0.41 0.42 0.36 0.36 0.37

0.0110 0.0145 0.0121 0.0095 0.0096 0.0040 0.0133 0.0085 0.0064 0.0166 0.0151 0.0138 0.0023 0.0145 0.0107 0.0215 0.0110 0.0061 0.0082 0.0187 0.0110 0.0296 0.0108 0.0043 0.0239 0.0110 0.0330 0.0258 0.0073 0.0134 0.0202 0.0211 0.0133 0.0110 0.0031 0.0083 0.0064 0.0044 0.0046


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

S# 10

11

12 13

S# 14

15

16

17

18

NVT parameters (Langevin, other=A, mH=1amu) dt (ps) tMD Tdamp gjf zero ∆t (ps) (ps) (ns) 0.00050 2.2 6.6 0.050 yes yes 0.00075 3.3 9.9 0.075 yes yes 0.00100 4.4 13.2 0.100 yes yes 0.00125 5.5 15.5 0.125 yes yes 0.00050 2.2 6.6 0.030 yes yes 0.00050 2.2 6.6 0.040 yes yes 0.00050 2.2 6.6 0.050 yes yes 0.00050 2.2 6.6 0.060 yes yes 0.00050 2.2 6.6 0.070 yes yes 0.00050 2.2 6.6 0.050 yes yes 0.00050 2.2 6.6 0.050 no yes 0.00050 2.2 6.6 0.050 yes yes 0.00050 2.2 6.6 0.050 yes no NPT parameters (Tdamp=100dt, Tloop=1ps, drag=0, other=A, mH=1amu) dt (ps) Pdamp Ploop nreset mtk tMD ∆t (ns) (ps) (ps) 0.00050 2.2 6.6 0.50 1 0 yes 0.00075 3.3 9.9 0.75 1 0 yes 0.00100 4.4 13.2 1.00 1 0 yes 0.00125 5.5 15.5 1.25 1 0 yes 0.00050 2.2 6.6 0.30 1 0 yes 0.00050 2.2 6.6 0.40 1 0 yes 0.00050 2.2 6.6 0.50 1 0 yes 0.00050 2.2 6.6 0.60 1 0 yes 0.00050 2.2 6.6 0.70 1 0 yes 0.00050 2.2 6.6 0.50 1 0 yes 0.00050 2.2 6.6 0.50 3 0 yes 0.00050 2.2 6.6 0.50 5 0 yes 0.00050 2.2 6.6 0.50 7 0 yes 0.00050 2.2 6.6 0.50 9 0 yes 0.00050 2.2 6.6 0.50 1 0 yes 0.00050 2.2 6.6 0.50 1 2200 yes 0.00050 2.2 6.6 0.50 1 4400 yes 0.00050 2.2 6.6 0.50 1 6600 yes 0.00050 2.2 6.6 0.50 1 8800 yes 0.00050 2.2 6.6 0.50 1 11000 yes 0.00050 2.2 6.6 0.50 1 0 yes 0.00050 2.2 6.6 0.50 1 0 no


Page 16 of 32

D0 (Å2/ps)

Q (eV)

ζ

75 124 116 166 84 111 75 90 82 75 42 75 fail D0 (Å2/ps)

0.38 0.41 0.41 0.42 0.39 0.40 0.38 0.40 0.39 0.38 0.36 0.38 fail Q (eV)

0.0033 0.0029 0.0017 0.0030 0.0019 0.0053 0.0033 0.0077 0.0079 0.0033 0.0040 0.0033 fail ζ

72 98 137 103 77 44 72 192 190 72 75 79 190 181 72 81 108 59 67 115 72 105

0.39 0.40 0.42 0.40 0.39 0.36 0.39 0.44 0.44 0.39 0.39 0.39 0.44 0.43 0.39 0.39 0.41 0.37 0.38 0.41 0.39 0.40

0.0129 0.0120 0.0057 0.0081 0.0448 0.0557 0.0129 0.0134 0.0172 0.0129 0.0131 0.0316 0.0250 0.0355 0.0129 0.0242 0.0105 0.0163 0.0214 0.0101 0.0129 0.0081

16

Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1 provides significant insights on MD simulations of diffusion. First, the error of the Arrhenius fit is significantly reduced in all of the MD simulations presented in Table 1 as compared to the ζ = 0.5188 and ζ = 0.0759 of the shorter time MD simulations shown in Fig. 2. As ζ is reduced, the activation energy Q increasingly converges towards ~0.41 eV (except for the mass effect). According to Fig. 3, ζ can always be reduced by increasing tMD. Hence, Table 1 strongly indicates that MD is a powerful and accurate method to calculate activation energy barrier of diffusion where converged results independent of the MD parameters can be achieved with a sufficiently long tMD. S#9 indicates that increasing the mass of the diffusion species causes a systematic decrease in both the pre-exponential factor D0 and activation energy barrier Q. The effects of mass on D0 can be understood because D0 scales with vibration frequency and vibration frequency reduces with increasing mass. On the other hand, within the mass range of the three hydrogen isotopes (hydrogen = 1 amu, deuterium = 2 amu, and tritium = 3 amu), the diffusion barriers are not too different. Hence, one strategy to calculate diffusion barriers for hydrogen is to artificially increase its mass to 2 or 3 amu because unlike mH = 1 amu which produces an error of ζ = 0.0110, increasing the mass above 2 seem to consistently lead to errors below the 0.01 margin at the same tMD. Note that MD simulations of the light hydrogen atoms must use extremely small time step size dt (e.g., dt 25 mm 15 mm 4 mm 2 mm 3 mm 3 mm 4 mm 99.999 % 99.999 % + oxides (water-treated) 99.8 % -----------------------------------------

D0 (Å2/ps) Q (eV) 9.530×102 0.46 4.580×102 0.38 1.520×103 0.55 2.110×104 0.69 1.540×105 0.84 4.220×106 1.11 6.000×105 0.92 1.900×103 0.41 1.100×105 0.74 2.510×106 1.100×103 2.100×103 1.101×103 2.000×102 9.200×103 6.100×103 2.600×103 1.300×105 1.750×100 1.200×109

0.93 0.42 0.47 0.49 0.52 0.57 0.57 0.61 0.69 0.17 1.45

literature Ichimura et al13,14

Papp et al15,16

Eichenauer et al17,18 Outlaw et al19 Matsuo et al20 Ishikawa et al21 Saitoh et al22 Hashimoto et al23 Csanady et al24 Young et al25 Ransley et al26

Table 2 confirms that the literature data for hydrogen diffusion in aluminium is extremely scattered. However, Ichimura et al13,14 studied the effects of grain boundaries and Papp et al15,16 studied the effects of aluminium purity. Their results indicate that for pristine materials with increasing grain size and purity, the hydrogen diffusion energy barrier increasingly approaches the 0.41 – 0.46 eV margin, in good agreement with our results. Interestingly, Table 2 also clearly shows a “compensation” effect13, that is, large activation energy barriers are always associated with large pre-exponential factors so that the overall diffusion coefficients do not vary much. This suggests that the measured large energy barriers pertain to traps rather than lattice sites. This is because to have a large overall diffusivity when the barrier is high, the activated hydrogen atom must diffuse a long distance before being de-


20

Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


activated again. Consistent with the this notion, the increase in energy barrier due to grain boundaries and impurities has been interpreted by the trapping of hydrogen atoms to these structural defects13,25,27,28. In this regard, it is the low end of the measured activation energies that pertain to the lattice barrier calculated in our simulations. The lowest energy barrier is 0.17 eV measured by Young et al25. This data point is much lower than all the other measurements. If this point is excluded, the lower end of the measured energy barriers are again in the 0.38 – 0.46 eV range, in agreement with our simulations. Confidence in the simulations may also stem from the interatomic potential that captures a variety of Al-H structures, as mentioned above. We are currently extending our methodologies to include various defects. A complete set of data will allow us to explore the interplay between diffusion barrier, defects, and interatomic potential. Such interplay is needed to better connect our results with the extremely low diffusion barrier measured by Young et al25, and other quantum mechanical calculations46. The main goal here is to identify if the lowest energy barrier seen in experiments corresponds to the lattice diffusion barrier. For instance, the energy barrier for surface H atoms to recombine to H2 molecules may be lower than the diffusion barrier. In addition, while grain boundaries can trap diffusion atoms from the bulk, the trapped atoms can still diffuse along grain boundaries where the diffusion barrier is in fact lower than in the bulk. Grain boundaries only increase diffusion barriers when there are diffusion-stopping traps on the grain boundaries, such as the junction points of the grain network13. Likewise, dislocation network can trap hydrogen atoms, but hydrogen diffusion along the open channel of straight edge dislocations can be faster. By performing the MD simulations (on both H2 recombination and diffusion) using systems containing the relevant defects, our methods will allow us to confidently study the effects of


21


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32

recombination, vacancies, grain boundaries, and dislocations etc. in the future so that this issue can be directly addressed. It is important to note that the experimental pre-exponential factors shown in Table 2 are in general larger than the calculated values shown in Table 1. Unlike diffusion energy barrier that can be directly tied to the energetics of the interatomic potential, the discrepancy in preexponential factor is less clear. Because experimental D0 is larger, there might be ballistic motion of hydrogen atoms (i.e., jump length is longer than the lattice size) in experiments that is not captured by simulations. Our highly converged results should establish an interesting topic for future studies to clarify this issue. V.

CONCLUSIONS We have demonstrated that robust MD methodologies can produce reliable diffusion

energy barriers from the Arrhenius equation and temperature-dependent diffusivity values. Most importantly, such methods can be generally applied to complex materials including defective metal alloys, hydrides, and oxides. Simulations using these methods lead to the following insights: 1. The time interval ∆t for tracking atom positions can be comparable or above the residence time of diffusion species without impacting the convergence of the results, and highly converged results can be achieved by extending the total simulated time tMD; 2. Parameters drag, Tdamp, Tloop, Pdamp, Ploop, nreset, and mtk do not impact results in a significant way, whereas gjf = yes is desirable and zero = yes is mandatory for the Langevin thermostat. Normally the default values of drag = 0, Tdamp = 100dt, Tloop = 1, Pdamp 1000dt, Ploop = 1, nreset = 0, mtk = yes, gjf = yes, and zero = yes can be used, although researchers can always make judicious choices based on Table 1;


22

Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3. Increasing the mass of the diffusion atoms reduces both pre-exponential factor and diffusion energy barrier, although the change of diffusion barrier is not significant for the mass range of the three hydrogen isotopes (hydrogen, deuterium, and tritium); 4. NVE ensemble produces highly consistent results. When only simulating a single diffusion atom, the NVT ensemble with a Langevin thermostat produces the lowest Arrhenius error as compared with other ensembles. However, we find that the Arrhenius error of NPT ensemble can be even lower if the system contains many diffusion atoms; 5. Our MD prediction of diffusion energy barrier of Q = 0.41 eV is in good agreement with the 0.38 – 0.46 eV experimental range for pristine samples with large grains and high purities. The MD value also agrees with the MS calculations. Our study helps to elucidate the inconsistence of the experimental results: the high measured energy barriers likely pertain to traps, rather than lattice sites; VI. ACKNOWLEDGEMENTS Sandia National Laboratories is a multiprogram laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the US Department of Energy's National Nuclear Security Administration under contract DE-AC0494AL85000. The authors gratefully acknowledge research support from the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, Fuel Cell Technologies Office, under Contract Number DE-AC04-94AL85000.

References 1 2

Schlapbach, L.; Zuttel, A. Hydrogen-Storage Materials for Mobile Applications. Nature 2001, 414, 353-358. Bhattacharyya, R.; Mohan, S. Solid State Storage of Hydrogen and Its Isotopes: An Engineering Overview. Ren. Sus. Ener. Rev. 2015, 41, 872-883.


23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3

4 5 6

7

8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

Page 24 of 32

Crosby, K.; Wan, X.; Shaw, L. L. Improving Solid-State Hydriding and Dehydriding Properties of the LiBH4 Plus MgH2 System with the Addition of Mn and V Dopants. J. Power Sources 2010, 195, 7380-7385. Chen, X.; Li, C.; Grätzel, M.; Kostecki, R.; Mao, S. S. Nanomaterials for Renewable Energy Production and Storage. Chem. Soc. Rev. 2012, 41, 7909-7937. Zheng, J.; Liu, X.; Xu, P.; Liu, P.; Zhao, Y.; Yang, J. Development of High Pressure Gaseous Hydrogen Storage Technologies. Inter. J. Hydro. Ener. 2012, 37, 1048-1057. Henkelman, G.; H. Jonsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978-9985. Henkelman, G.; Uberuaga, B. P.; Jonsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901-9904. Nakano, A. A Space-Time-Ensemble Parallel Nudged Elastic Band Algorithm for Molecular Kinetics Simulation. Comp. Phys. Comm. 2008, 178, 280-289. LAMMPS download site: lammps.sandia.gov. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comp. Phys. 1995, 117, 1-19. Hoover, W. G. Phys. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. B 1985, 31, 1695-1697. Schneider, T.; Stoll, E. Molecular-Dynamics Study of a Three-Dimensional OneComponent Model for Distortive Phase Transitions. Phys. Rev. B 1978, 17, 1302-1322. Ichimura, M.; Sasajima, Y.; Imabayashi, M. Grain-Boundary Effect on Diffusion of Hydrogen in Pure Aluminum. Mater. Trans., JIM 1991, 32, 1109-1114. Ichimura, M.; Imabayashi, M. Measurement of the Diffusion-Coefficient and Solubility of Hydrogen in Solid Aluminum. J. Jpn Inst. Metals 1979, 43, 876-883. Papp, K.; Kovacs-Csetenyi, E. Diffusion of Hydrogen in High-Purity Aluminum. Scr. Metall. 1981, 15, 161-164. Papp, K.; Kovacs-Csetenyi, E. Diffusion of Hydrogen in Solid Aluminum. Scr. Metall. 1977, 11, 921-923. Eichenauer, W.; Hattenbach, K.; Pebler, A. Die Loslichkeit won Wasserstoff in Festem und Flussigem Aluminium. Z. Metallk. 1961, 52, 682-684. Eichenauer, W.; Pebler, A. Messung des Diffusionskoeffizienten und der Loslichkeit von Wasserstoff in Aluminium und Kupfer. Z. Metallk. 1957, 48, 373-378. Outlaw, R. A.; Peterson, D. T.; Schmidt, F. A. Diffusion of Hydrogen in Pure Large Grain Aluminum. Scr. Metall. 1982, 16, 287-292. Matuso, S.; Hirata, T. Diffusion of Hydrogen in Aluminium. Trans. Nat. Res. Inst. Met. 1969, 11, 88-92. Ishikawa, T.; McLellan, R. B. The Diffusivity of Hydrogen in Aluminum. Acta Metall. 1986, 34, 1091-1095. Saitoh, H.; Iijima, Y.; Tanaka, H. Hydrogen Diffusivity in Aluminum Measured by a GlowDischarge Permeation Method. Acta Metall. Mater. 1994, 42, 2493-2498. Hashimoto, E.; Kino, T. Hydrogen Diffusion in Aluminum at High-Temperatures. J. Phys. F: Met. Phys. 1983, 13, 1157-1165.


24

Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


24 Csanady, A.; Papp, K.; Pasztor, E. The Relation between Hydrogen Desorption and the Surface Conditions of High-Purity Aluminum. Mater. Sci. Eng. 1981, 48, 35-39. 25 Young, G. A.; Scully, J. R. The Diffusion and Trapping of Hydrogen in High Purity Aluminum. Acta Mater. 1998, 46, 6337-6349. 26 Ransley, C. E.; Talbot, D. E. J. Wasserstoff-Porositat in Metallen unter Besonderer Berucksichtigung des Aluminiums und Seiner Legierungen. Z. Metallk. 1955, 46, 328-337. 27 McLellan, R. B. Hydrogen Diffusion in Aluminum. Scr. Metall. 1983, 17, 1237-1240. 28 Edwards, R. A. H.; Eichenauer, W. Reversible Hydrogen Trapping at Grain-Boundaries in Superpure Aluminum. Scr. Metall. 1980, 14, 971-973. 29 Gunaydin, H.; Barabash, S. V.; Houk, K. N.; Ozolins, V. First-Principles Theory of Hydrogen Diffusion in Aluminum. Phys. Rev. Lett. 2008, 101, 075901-4. 30 Wolverton, C.; Ozoliņš, V.; Asta, M. Hydrogen in Aluminum: First-Principles Calculations of Structure and Thermodynamics. Phys. Rev. B 2004, 69, 144109-16. 31 Hirsch, J.; Al-Samman, T. Superior Light Metals by Texture Engineering: Optimized Aluminum and Magnesium Alloys for Automotive Applications. Acta Mater. 2013, 61, 818843. 32 Marlaud, T.; Malki, B.; Henon, C.; Deschamps, A.; Baroux, B. Relationship between Alloy Composition, Microstructure and Exfoliation Corrosion in Al-Zn-Mg-Cu Alloys. Corr. Sci. 2011, 53, 3139-3149. 33 Kamoutsi, H.; Haidemenopoulos, G. N.; Bontozoglou, V.; Pantelakis, S. Corrosion-Induced Hydrogen Embrittlement in Aluminum Alloy 2024. Corros. Sci. 2006, 48, 1209-1224. 34 Marlaud, T.; Malki, B.; Deschamps, A.; Baroux, B. Electrochemical Aspects of Exfoliation Corrosion of Aluminium Alloys: the Effects of Heat Treatment. Corros. Sci. 2011, 53, 1394-1400. 35 Michler, T.; San Marchi, C.; Naumann, J.; Weber, S.; Martin, M. Hydrogen Environment Embrittlement of Stable Austenitic Steels. Inter. J. Hydro. Ener. 2012, 37, 16231-16246. 36 Zhou, X. W.; Ward, D. K.; Foster, M. E. An Analytical Bond-Order Potential for the Aluminum Copper Binary System. submitted. 37 Zhou, X. W.; Ward, D. K.; Foster, M.; Zimmerman, J. A. An Analytical Bond-Order Potential for the Copper-Hydrogen Binary System. J. Mater. Sci. 2015, 50, 2859-2875. 38 Donnay, J. D. H.; Ondik, H. M. Crystal Data, Determinative Tables, 3rd ed., Vol. 2 (Inorganic Compounds); U. S. Department of Commerce, National Bureau of Standards, and Joint Committee on Power Diffraction Standards, U.S.A., 1973. 39 Zhou, X. W.; Ward, D. K.; Zimmerman, J. A.; Cruz-Campa, J. L.; Zubia, D.; Martin, J. E.; van Swol, F. An Atomistically Validated Continuum Model for Strain Relaxation and Misfit Dislocation Formation. J. Mech. Phys. Solids. 2016, accepted. 40 Royset, J.; Ryum, N. Some Comments on the Misfit and Coherency Loss of Al3Sc Particles in Al-Sc Alloys. Scr. Mater. 2005, 52, 1275-1279. 41 Reif, F. Fundamentals of Statistical and Thermal Physics; McGraw Hill: New York, U.S.A., 1965. 42 Gronbech-Jensen, N.; Farago, O. A Simple and Effective Verlet-Type Algorithm for Simulating Langevin Dynamics. Mol. Phys. 2013, 111, 983-991. 43 Gronbech-Jensen, N.; Hayre, N. R.; O. Farago, O. Application of the G-JF Discrete-Time Thermostat for Fast and Accurate Molecular Simulations. Comp. Phys. Comm. 2014, 185, 524-527.


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 32

44 Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nose-Hoover Chains: the Canonical Ensemble via Continuous Dynamics. J. Chem. Phys. 1992, 97, 2635-2643. 45 Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117-1157. 46 Verners, O.; Psofogiannakis, G.; van Duin, A. C. T. Comparative Molecular Dynamics Study of fcc-Al Hydrogen Embrittlement. Corr. Sci. 2015, 98, 40-49.


26

Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Figure


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Illustration of the jump between the nearest octahedral and tetrahedral sites in fcc Al. 101x106mm (300 x 300 DPI)


Page 28 of 32

Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) Mean square displacement as a function of time and (b) the Arrhenius plot for two sets of MD simulations conducted using the NVT ensemble at two different temperature ranges and ∆t values, but the same time step size of dt = 0.001 ps and total MD time of tMD = 0.88 ns. 152x66mm (300 x 300 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a) Arrhenius error ζ and (b) activation energy Q as a function of MD time tMD, obtained from MD simulations using the NVT ensemble at a time step size of dt = 0.001 ps, an output time interval of ∆t = 4.4 ps, and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K (T/Q = 1220 – 1707 K/eV). 67x29mm (600 x 600 DPI)


Page 30 of 32

Page 31 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a) Mean square displacement as a function of time and (b) the Arrhenius plot for MD simulations conducted using the NVT ensemble at a time step size of dt = 0.001 ps, an output time interval of ∆t = 4.4 ps, a total MD time of tMD = 13.2 ns, and nine temperatures T = 500, 525, 550, 575, 600, 625, 650, 675, and 700 K. 67x29mm (600 x 600 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table of Content Figure 44x24mm (300 x 300 DPI)


Page 32 of 32

Molecular Dynamics Simulations of Hydrogen Diffusion in Aluminum

Recommend Documents