Multiscale Molecular Dynamics Approach to Energy Transfer in

Subscriber access provided by READING UNIV

Article

A multiscale molecular dynamics approach to energy transfer in nanomaterials John Michael Espinosa-Duran, Yuriy V. Sereda, Andrew Abi-Mansour, and Peter Ortoleva J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00702 • Publication Date (Web): 30 Nov 2017 Downloaded from http://pubs.acs.org on January 3, 2018

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.

Journal of Chemical Theory and Computation 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 37

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

Journal of Chemical Theory and Computation

A multiscale molecular dynamics approach to energy transfer in nanomaterials John M. Espinosa-Duran, Yuriy V. Sereda, Andrew Abi-Mansour† and Peter Ortoleva* Center for Theoretical and Computational Nanoscience, Department of Chemistry, Indiana University, Bloomington, IN 47405, USA. KEYWORDS: Energy transfer, multiscale theory, field variables, nanoscale materials.

ABSTRACT: After local transient fluctuations are dissipated, in an energy transfer process, a system evolves to a state where the energy density field varies slowly in time relative to the dynamics of atomic collisions and vibrations. Furthermore, the energy density field remains strongly coupled to the atomic scale processes (collisions and vibrations), and it can serve as the basis of a multiscale theory of energy transfer. Here, a method is introduced to capture the long scale energy density variations as they coevolve with the atomistic state in a way that yields insights into the basic physics and implies an efficient algorithm for energy transfer simulations. The approach is developed based on the N-atom Liouville equation and an interatomic force field, and avoids the need for conjectured phenomenological equations for energy transfer and other processes. The theory is demonstrated for sodium chloride and silicon dioxide nanoparticles immersed in a water bath via molecular dynamics simulations of the energy transfer between a nanoparticle and its aqueous host fluid. The energy density field is computed

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 37

for different sets of symmetric grid densities, and the multiscale theory holds when slowly varying energy densities at the nodes are obtained. Results strongly depend on grid density and nanoparticle constituent material. A non-uniform temperature distribution, larger thermal fluctuations in the nanoparticle than in the bath and enhancement of fluctuations at the surface, which are expressed due to the atomic nature of the systems, are captured by this method, and not by phenomenological continuum energy transfer models.

1. INTRODUCTION Descriptions of many-particle systems in terms of field variables have been the basis of the classical phenomenological theory of continuous media since XIX century.1 The more traditional statistical mechanical approach is based on deriving continuum theories for manyparticle Newtonian or quantum mechanical systems using techniques such as projection operators2 or asymptotic expansions in a space-time scale ratio.3-5 Field variables have also been used to formulate multiscale theories of molecular systems,6-7 e.g. in the study of liquid crystals where the chemistry of molecules is orientation dependent,6 and served as the starting point of a multiscale approach if the field variables are slow relative to atomic vibrations and collisions.6 In the present and other multiscale methods, field variables are co-evolved with the atomistic microstate avoiding the need for phenomenological laws.7 To simulate energy transfer at nanoscale, a quantum mechanical description can be used,8-10 but with prohibitive computational cost for relatively large systems. However, classical molecular dynamics and Monte Carlo simulations11-13 are an alternative to analyze the classical (non-electronic) component of energy transfer, but are computationally too expensive for long time scale processes. To simulate large systems for long periods of time, multiscale methods14-16 can be very useful. A coarser approximation is developed using models based on quantum17-18 or


2

Page 3 of 37

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


classical macroscopic19-21 phenomenological laws. These models work under certain conditions but do not provide detailed local or atom-resolved information, nor do they fully account for the highly fluctuating nature of nanosystems. These challenges have been addressed limitedly by calibrating mathematical models using data from classical MD22 and quantum simulations,23 but the computational cost of these simulations is high and can lead to inaccuracies due to the lack of recalibration or feedback to readjust the model parameters. Some approaches make a complete description of the system and the energy exchange in a framework where a range of theories is included to describe them on different time and length scales, going from all-atom models to classical phenomenological ones.24 Unfortunately, the framework was found to be developed only for one type of problem: welding. A top-down approach presented elsewhere,25 couples a continuum description of the system in a Finite Element model with classical MD, using the temperature defined at the nodes of the mesh as the coarse-grained (CG) variable. In the latter work, the authors did not keep track of the potential energy changes, and temperature was the unique CG variable, thereby missing a very relevant portion of the information. The development of continuum-microstate coevolution approaches has been placed within the framework of multiscale perturbation theory; the result was Langevin equations for CG continuum variables that co-evolve with the many-particle microstate.26-27 In previous studies, multiscale factorization (MF)28 was developed and used to coevolve the microstate and the CG state.7, 28-29 MF carries out the microstate and CG coevolution taking into account the timescale separation between atomistic and CG dynamics, as well as the stationarity principle of atomistic fluctuations.28-29 Most closely related to the energy transfer approach presented here is the previous work where the internal energy of a system is considered to be a slowly-varying CG variable that coevolves with the rapidly fluctuating atomistic degrees of freedom.29 The analysis


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 37

of nanosystems using density-like field variables, within a CG-microstate coevolution scheme that follows from the all-atom physics, accounts for dynamical phenomena that cannot be described by classical phenomenological continuum models which miss the intrinsic random atomic vibrations and Boltzmann distribution of velocities. For example, the formation of waves on the surface of a nanoscale material which has been excited by laser radiation30 cannot be described by current continuum phenomenological methods, since they are unable to describe atomic fluctuations which give rise to this phenomena. However, density-like variables can provide a framework to analyze surface-enhanced fluctuations as will be shown here. Here, an efficient algorithm for simulating energy transfer based on the coevolution of the energy density field and the microstate using the MF approach28 is presented. Since the energy can naturally be described as a field or density-like variable, the energy density (field) is introduced in a manner that is consistent with Newtonian mechanics (via the Liouville equation). The algorithm is developed for all-atom models, but it can be used with CG structural models like Martini,15 and does not require energy transfer phenomenological laws. The algorithm is proposed in a way that preserves the all-atom state of the system, but guides its evolution with the slowly varying density field. The energy density field is computed at the nodes of a grid, whose positions have no direct relations with the atomic positions. A Gaussian function centered at the nodes is used to weight the individual atomic contributions to the energy density, thereby obtaining a slow CG variable dynamics. The algorithm is cast in terms of classical all-atom force fields and therefore does not include any electron-related dynamic phenomena that could contribute to the energy transfer. Therefore, it works for classical systems, mainly electrical insulators, where the energy transfer is mostly due to phonons, and not electrons. The algorithm


4

Page 5 of 37

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


proposed in this work is based on the dynamics of the total energy density (i.e., potential plus kinetic energy), and not only the temperature (kinetic energy) as others have done.25 Here, the coevolution is achieved via the alternation of microstate dynamics and CG evolution. The microstate dynamics (achieved computationally via MD) is used to obtain information about the local energy exchange between system and bath occurring on a short time scale. The collected information is used to evolve the CG state (i.e., the energy density field), on a longer time scale. At the end of this CG state evolution, an updated energy density field is obtained and the all-atom state is recreated in a manner that is consistent with it. This is done in a way that the recreated or new microstate is free of unlikely configurations (i.e., does not have extremely stretched or compressed bonds). The efficiency of the method is related to the stationarity property of the atomistic fluctuations by which the system explores a rich ensemble of CG velocities (rates of change of the CG variable) in a short time scale compared to the characteristic time of energy evolution.28-29 Due to this, one can take longer time steps in this discretized evolution algorithm. This paper is organized as follows. The energy density is introduced as a field variable (Section 2), and the Liouville operator is unfolded to allow the coevolution (Section 3) via the multiscale Lie-Trotter factorization (Section 4). This approach is evaluated for a set of nanostructures immersed in an aqueous bath (Section 5), and conclusions are drawn (Section 6). 2. THEORY 2.1 Field variable formulation Consider an -atom system evolving via Newtonian mechanics. The energy density field

Φ () is introduced such that the total energy of a system is


5


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 6 of 37

(1)

Due to the incapacity of sampling all the CG space with infinite resolution, only discrete

points in the CG space are considered Φ → {Φ }. The CG Space is a scale or

projected version of real space such that any coordinate in real space can be associated with a coordinated in CG space, i.e., = . ≤ 1 is a scaling factor that expresses the multiscale

nature of the problem and projects the atomic coordinates from real space to CG space . The

positions of these discrete points have no direct relation with the coordinates of the atoms in the system and therefore can be located anywhere in the simulation box. In this case, for

simplicity, it is assumed that those discrete points are uniformly distributed in a three dimensional grid, therefore called grid points or nodes (as those in Smoothed Particle Hydrodynamics31 or Finite Element Method32). Then, the integral in eq. (1) becomes a summation over the grid nodes (considering a differential volume element ∆ ), ≈ Φ ( )∆ .

(2)

To account for the interatomic interactions around every grid node, a Gaussian function ( − ) centered at the grid node with coordinates is introduced such that. − = 1,

− =

1

(2") #

' %()* )+ %(& $ ,-+ ,

(3)

(4)

where is the coordinate of the .-th atom whose energy contributes to the energy density computed for the grid node located at and # is the width of the Gaussian. The parameter in

general is taken to be the average nearest-neighbor atom distance divided by a long characteristic


6

Page 7 of 37

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


length. The latter might be, for example, the size of the region wherein an initial input of energy from a laser pulse is injected during a heating experiment.33 Then, as earlier,29 via a simplified approach one could follow the time course of Φ ( ) as energy is transferred from an excited nanostructure to the host medium, or in the opposite direction. After a transient period, Φ ( )

is expected to evolve slowly through energy redistribution within the nanostructure, and in a longer time scale, across the host medium. The energy density Φ ( ) must be related to the all-atom state of the system, Γ (the microstate which corresponds to the 6N atomic positions and momenta in real space, with N equal to the number of atoms in the system), such that Φ ( ) → Φ ( , Γ) . Then, Φ ( , Γ) is

expressed as the summation of the energy experienced by every atom i, (Γ ), weighted by the Gaussian function introduced above, Φ , Γ = − ⋅ (Γ )

(5)

For computational efficiency, only atoms in the vicinity of the grid node located at could be

considered in this summation. The size of the vicinity around can be chosen as a function of #

and the internode distance, considering that the integral of a Gaussian function converges within 3 # to the 99.73% of its value (0.27 % error). Then, the energy experienced by an atom can be written as the sum of the atomic kinetic and potential energies, K and U, respectively, (Γ ) = 2(3 ) + 5( )

(6)

where Γ = { , 3 } with and 3 being the atomic positions and momenta, respectively. The kinetic energy of the .-th atom is expressed as a function of its momentum (3 ) and mass (6 ),


7


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 2(3 ) = 26 ,

Page 8 of 37

(7)

The potential energy experienced by the .-th atom 5( ) depends on its position and the

positions of all the atoms in the system (which is a subset of Γ). It can be expressed as the summation of single, pair, and more particle interactions. Considering that in this work standard interatomic force fields are used, up to 4 particle interactions are accounted for: 1 1 1 5( ) = 5 (7) ( ) + 5 (,) , 8 + 5 ( ) , 8 , : + 5 (5 ps for all the initial integration times tested. However, Figure 4.a shows that Φ has anomalous very large fluctuations during the first picoseconds. These results suggest that for this system under the given conditions, the length of the MD simulation to evolve the microstate (] JK)L ), and during which Φ is calculated, should have a duration δ>5 ps, for t>5 ps. In the 3x3x3 case (Figure 4.b), for t=0 ps, which implies taking into account the very rapid initial change, the integral does not converge within the plotted time (60 ps). For t=5 ps, and t=15 ps, the integral starts reaching convergence by the end of the plotted time. For t=25 ps and t=35 ps, the integral converges for τ>5 ps, which is in agreement with the results presented in Figure 2.b which indicate that for the 3x3x3 grid, Φ significantly slows after 20 ps. These results suggest that for this system under the given conditions, the microstate (] JK)L ) simulations should have a duration δ>5 ps, for t>15 ps.

a)

b)

Figure 4. Stationarity integral I (t,τ) (eq. (24)) for the cooling of an NaCl nanocrystal immersed

in a water bath, for Φ located at: a) the bottom back corner of a 2x2x2 grid with σ=29 Å; b) the central node of a 3x3x3 grid with σ=19.5 Å.

Since the lower limit for δ is known for slow and smooth Φ , the CG evolution

(extrapolation) of Φ was studied for different combinations of δ (MD time) and ∆ (CG time), ACS Paragon Plus Environment

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

Page 24 of 37

and therefore different speed-up factors k=∆/δ. δ was varied from 5 ps to 10 ps, and ∆ from 10 ps to 90 ps, (and k from 2 to 10). As shown in Figure 3, the CG evolution is coherent, i.e., readily extrapolatable over long time intervals of duration ∆ after the initial rapidly-fluctuating transient t=10 ps for the 2x2x2 grid, and from t=20 ps for the 3x3x3 grid (obtained from Figure 4 analysis). Figure 5 shows cases for which the multiscale theory (extrapolation of Φ ) was applied achieving speed-up over MD and relatively low error (as compared with the continuous Φ from the MD trajectory). To do the extrapolation in each ∆ time-step, in this case, the change of Φ (∆Φ ) during δ was computed as the difference between the average first and last 10 saved frames in δ. No energy redistribution was applied after each ∆ time step.

a)

b)

Figure 5. Extrapolation of Φ computed at all the nodes in different grid densities given different theoretical parameters, for the cooling of an NaCl nanocrystal in water. The

extrapolated Φ are in black (linear segments of duration ∆ each), and the exact values from MD are in different colors. a) For the 2x2x2 grid, with CG evolution after t=10 ps for δ=9 ps and ∆=27 ps (k=3). b) For the 3x3x3 grid, with CG evolution after t=20 ps for δ=8 ps and ∆=56 ps (k=7). For the 2x2x2 grid, Figure 5.a shows the extrapolated Φ in black, for δ=9 ps and ∆=27 ps (k=3). The maximum speed-up that maintained accuracy for this grid was k=8, for δ=8 ps and


24

Page 25 of 37

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


∆=64 ps. For the 3x3x3 grid, Figure 5.b shows the extrapolated Φ in black for δ=8 ps and ∆=56 ps (k=7). The latter correspond to the maximum speed-up that maintained accuracy for this grid. Other cases evaluated are presented in Figure SI.3. From these results it is observed that the higher the speed-up, the greater the extrapolation error, so it is crucial find a tradeoff between error and speed-up over MD for each problem and system.

Figure 6. Redistribution of increased kinetic energy into potential energy to recover a microstate consistent with the updated energy density for small NaCl nanocrystal. Considering eq. (32), to obtain a new microstate consistent with the updated energy density Φ (G + ∆), it is necessary to scale the atomic velocities according to eq. (33), and then run an NVE simulation for the redistribution of kinetic energy into potential. Figure 6 shows the redistribution of kinetic energy into potential energy during an NVE simulation, after scaling up the atomic velocities of the NaCl nanocrystal from 50 K to 300 K. The nanoparticle redistributes locally (intramolecular mainly) the increased kinetic energy into potential energy very quickly, i.e. between 200 fs and 400 fs. This process is coupled to the energy transfer between the nanocrystal and the heat bath, which occurs in a much longer time scale. After this short energy redistribution, the algorithm repeats: an NVE simulation of duration δ is used to accumulate the data needed to compute Φ by eq. (5); then, Φ is evolved from time G to G + ∆ via the MF ACS Paragon Plus Environment

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 37

algorithm using eq. (30). Finally, atomic velocities are scaled according to Φ (G + ∆) and then the microstate is reconstructed. 4.2 Heating of a cool SiO2 nanoparticle in room temperature water The other problem analyzed here is the heating of a cooled SiO2 nanoparticle immersed in room temperature water. Applicability of the theory for this system was tested at the same grid densities used for the NaCl system (see Figure 7.a and Figure SI.4). Since both systems have the same number of atoms and roughly the same dimensions, the location of the nodes relative to the nanoparticle was considered to be the same. The two nanoparticles are made of a different material, thus differences arose from the different set of parameters used to represent the atomic interactions. Figure 7.b shows the temporal evolution of the temperature in the composite system (total), nanoparticle (SiO2), and the bath (water) for this simulation. In this case, the system reached thermal equilibrium in approximately 50 ps. The thermal fluctuations in the nanoparticle are much larger than in the bath, as thermodynamically expected, since the bath has 10 times more atoms.

a)

b)

Figure 7. a) Small SiO2 nanoparticle immersed in a water bath (for water, only the oxygen atoms are shown as red dots) for a 4x4x4 grid (light green spheres). b) Temperature profiles for the first 100 ps of the nanoparticle heating simulation.


26

Page 27 of 37

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)

b)

c)

d)

Figure 8. Temporal evolution of Φ during the heating of a SiO2 nanoparticle in water. The

energy density Φ at each node was plotted for different grid densities. a) 2x2x2, b) 3x3x3, c) 4x4x4, d) 5x5x5.

It was concluded from the results presented in the previous section that σ should be approximately half the internode distance, therefore this is taken into account in the following analysis since the systems have similar size. Figure 8 shows Φ for the different grid densities;

each one of the curves in the plots corresponds to Φ at a given node. For the 2x2x2 and 4x4x4 grids, fluctuations in Φ were large, so they are not considered slowly varying CG variables. Therefore, it was not possible to extrapolate their value using the proposed theory or to obtain a speed-up over conventional MD. However, the formalism for these two grid densities captured


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

Page 28 of 37

surface-enhanced fluctuations which are key to understand solid-liquid interfaces at nanoscale. This is notable, given the fact that they are not easily accessible experimentally nowadays and are neglected in most phenomenological continuum models. For the 3x3x3 and 5x5x5 grids, the fluctuations were much smaller and the CG variables were much slower and smoother than in the other cases; therefore, the theory was applied to extrapolate the CG variables and a speed-up over MD could be obtained. Regarding identifying the general location of the grid nodes using their energy density dynamics observed from the plots, one can say that they follows the same trend described for the cooling experiment in the Section SI.2, but from top to bottom of each plot.

Figure 9. Stationarity integral I (t,τ) given Φ at the central node of the SiO2 nanoparticle for a

3x3x3 grid with σ=19.2 Å. It shows that the time required for the integral to converge is τ >10ps, for t > 30ps.

Since the thermal properties of this nanoparticle are different from those of the previous one, the stationarity time is different. The following analysis was carried out only for the 3x3x3 grid since slowly varying CG variables were obtained and it had a treatable amount of nodes. Figure 9 shows the stationarity integral (eq. (24)) solved for different initial times t, for the central node in a 3x3x3 grid for σ=19.2 Å. The results indicate that the integral converges for


28

Page 29 of 37

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


τ>10 ps for t> 40 ps. These suggest that for this system under the given conditions, the length of the MD simulation to compute ] JK)L should be δ> 10 ps. Given the lower limit for δ from the stationarity integral, for the grid densities for which Φ was slow and smooth, the evolution of the CG variables was studied for different combinations of δ and ∆ (and k). The parameter δ was varied from 1 ps to 20 ps and ∆ from 2 ps to 60 ps (and k from 2 to 10). Figure 10 shows cases for which the multiscale theory (extrapolation) was successfully applied and a speed-up with relatively low error with respect to the results derived from MD was obtained. In these cases, the CG evolution was applied after t=40 ps (to avoid the initial rapidly-changing transient) as was found from Figure 9. Using the MD trajectory to compute Φ for the 3x3x3 grid, ∆Φ was computed as k times the difference between the average of the first and last 5 saved frames in δ. Figure 10 shows the extrapolated Φ in black for δ=1 ps and ∆=2 ps (k=2), and for δ=10 ps and ∆=30 ps (k=3). The first case shows an exceptional situation where δ is smaller than the suggested from the convergence of the stationarity integral (Figure 9). This is possible since the CG variables fluctuate in a time scale comparable to τ (the stationarity time). Therefore to properly capture and predict its fluctuations, it is required to use a smaller time δ. Then, for δ=7 ps and ∆=35 ps (k=5) (Figure SI.5), the maximum speed-up over MD that retained accuracy was achieved for this system. In this case, as in the previous section, it is observed that the higher the speed-up, the higher the error in the extrapolation.


29


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)

Page 30 of 37

b)

Figure 10. Extrapolation of Φ computed at all the nodes of different grid densities given different theoretical parameters, for the heating of an SiO2 nanoparticle in water. The

extrapolated Φ values are in black, and the dynamic values from MD are in different colors. For a 3x3x3 grid (after t= 40 ps): a) δ=1 ps and ∆=2 ps (k=2); b) δ=10 ps and ∆=30 ps (k=3).

Figure 11 shows redistribution of increased kinetic energy (from 50 K to 300 K) into potential to recover a microstate consistent with the updated energy density for the SiO2 nanoparticle. The nanoparticle redistributes (mainly intramolecularly) the kinetic energy into potential in between 400 fs and 600 fs, almost two times slower than in the NaCl nanocrystal. This time difference could be due to: a) the different structural symmetry, and b) different strength and nature in the intramolecular interactions. In the first case, the NaCl nanoparticle is a crystal and the SiO2 nanoparticle is an amorphous solid, and in a crystal, the phonons responsible for the energy transfer are able to transfer their momentum (kinetic energy) to the surroundings faster than amorphous solids, due to the higher atomic order. In the second case, the NaCl nanoparticle is held together by Coulomb interactions and the SiO2 nanoparticle is held by covalent bonds.


30

Page 31 of 37

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


Figure 11. Redistribution of increased kinetic energy into potential energy to recover a microstate consistent with the updated energy density for a small SiO2 nanoparticle. 5. CONCLUSIONS A multiscale theory to study energy transfer in nanoscale systems is proposed based on the energy density as a coarse-grained field variable and on the multiscale factorization method to evolve the system in time. The theory is verified for cubic nanostructures individually immersed in a water bath (with 10 times more atoms). The nanostructures considered were an NaCl nanocrystal and an SiO2 amorphous nanoparticle. Both directions of energy transfer, from the bath to the nanoparticle and vice-versa, were studied. The theory uses the energy density as a field variable assuming it is slowly varying in time. For the volume of the composite systems, different 3D symmetric grids were tested and the CG variables were computed at the grid nodes; each grid node position is not related to positions of the atoms. The stochastic nature of nanoscale systems for which the surface to volume ratio is high and atomic processes are very relevant requires the use of all-atom multiscale approaches for computational studies. Phenomena like Boltzmann distribution of velocities, surface-enhanced fluctuations, and size-dependent thermal fluctuation are relevant for nanoscale systems and are


31


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 32 of 37

captured by the multiscale formalism proposed here. Note that these phenomena are not included usually in continuum phenomenological approaches. Based on the present simulations, the theory can be applied to bulky nanostructures, with ionic and covalent interactions. The ideal Gaussian width σ was suggested to be around half the internode distance for which CG variables were slowly varying and energy profiles could be captured by the theory. The CG evolution should be started after the initial rapidly fluctuating transient has ceased, since the convergence of stationarity integral strongly depends on the time G chosen to start the integration. The maximum speed-up achieved for the systems studied here was k=8. Future studies should assess the applicability of the method to phenomena such as surface-enhanced fluctuations, as well as for larger and more slowly varying systems. The process of local redistribution of kinetic energy into potential was shown to happen quickly, i.e., in a few hundred MD time steps. Therefore, rescaling the velocities to scale the kinetic energy assuming that any energy change is a change in kinetic energy likely will not affect the overall evolution. Supporting Information. Images of the systems simulated in this paper showing the different grids tested are included. A description of how the energy density dynamics at a given node can be related to the node location in the cooling problem is provided. Additional plots of the energy density per node for the cooling and heating simulations are included for different grid densities and speed-up factors. The supporting information file is available free of charge at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author


32

Page 33 of 37

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


*

E-mail: [email protected]

Present Address †

Center for Materials Science and Engineering, Merck & Co., Inc., West Point, PA 19486.

ACKNOWLEDGMENT The authors acknowledge the National Science Foundation (NSF, DMREF-1533988) for support. The authors also acknowledge support by the IU College of Arts and Sciences via the Center for Theoretical and Computational Nanoscience, and Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute, and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at IU is also supported in part by Lilly Endowment, Inc. The work done by J.M.E.D. was supported by the Colombian Administrative Department of Science, Technology, and Innovation – COLCIENCIAS – through the Doctorates Abroad Program 568-2012. REFERENCES 1. Maxwell, J. C., A dynamical theory of the electromagnetic field. Philosophical Transactions of the Royal Society of London 1865, 155, 459-512. 2. Zwanzig, R., Nonequilibrium Statistical Mechanics. Oxford University Press: 2001. 3. Deutch, J. M.; Oppenheim, I., The Lennard-Jones Lecture. The concept of Brownian motion in modern statistical mechanics. Faraday Discuss. Chem. Soc. 1987, 83, 1-20. 4. Shea, J. E.; Oppenheim, I., Fokker-Planck equation and non-linear hydrodynamic equations of a system of several Brownian particles in a non-equilibrium bath. Physica A 1997, 247, 417-443. 5. Ortoleva, P. J., Nanoparticle dynamics: A multiscale analysis of the Liouville equation. J. Phys. Chem. B 2005, 109, 21258-21266. 6. Shreif, Z.; Pankavich, S.; Ortoleva, P. J., Liquid-crystal transitions: A first-principles multiscale approach. Phys. Rev. E 2009, 80, 031703. 7. Abi Mansour, A.; Ortoleva, P. J., Implicit Time Integration for Multiscale Molecular Dynamics Using Transcendental Padé Approximants. Journal of Chemical Theory and Computation 2016, 12, 1965-1971. 8. Gerber, R. B.; Buch, V.; Ratner, M. A., Time‐dependent self‐consistent field approximation for intramolecular energy transfer. I. Formulation and application to


33


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 34 of 37

dissociation of van der Waals molecules. The Journal of Chemical Physics 1982, 77, 30223030. 9. High, J. S.; Rego, L. G. C.; Jakubikova, E., Quantum Dynamics Simulations of Excited State Energy Transfer in a Zinc–Free-Base Porphyrin Dyad. The Journal of Physical Chemistry A 2016, 120, 8075-8084. 10. Arango, C. A.; Kennerly, W. W.; Ezra, G. S., Semiclassical IVR approach to rotational excitation of non-polar diatomic molecules by non-resonant laser pulses. Chemical Physics Letters 2006, 420, 296-303. 11. Frenkel, D.; Smit, B., Understanding Molecular Simulation: From Algorithms to Applications. Academic Press 2001. 12. Gould, H.; Tobochnik, J., Statistical and Thermal Physics: With Computer Applications. 9780691137445: 2010. 13. Tuckerman, M. E., Statistical Mechanics: Theory and Molecular Simulation. Oxford University Press, USA: New York, 2010. 14. Jaramillo-Botero, A.; Nielsen, R.; Abrol, R.; Su, J.; Tod Pasca; Mueller, J.; III, W. A. G., First-Principles-Based Multiscale, Multiparadigm Molecular Mechanics and Dynamics Methods for Describing Complex Chemical Processes. Top Curr Chem 2012, 307, 42. 15. Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H., The MARTINI forcefield: coarse grained model for biomolecular simulations. J. Phys. Chem. B 2007, 111, 7812-7824. 16. Tuckerman, M.; Berne, B. J.; Martyna, G. J., Reversible multiple time scale molecular dynamics. Journal of Chemical Physics 1992, 97, 1990-2001. 17. Lake, R.; Datta, S., Energy balance and heat exchange in mesoscopic systems. Physical Review B 1992, 46, 4757-4763. 18. Galperin, M.; Nitzan, A.; Ratner, M. A., Heat conduction in molecular transport junctions. Physical Review B 2007, 75, 155312. 19. Cao, B.-Y.; Guo, Z.-Y., Equation of motion of a phonon gas and non-Fourier heat conduction. Journal of Applied Physics 2007, 102, 053503. 20. Xing, Z.; Koji, T.; Motoo, F., Charge and Heat Transport in Polycrystalline Metallic Nanostructures Chinese Physics Letters 2008, 25, 3360-3363. 21. Luo, T.; Chen, G., Nanoscale heat transfer - from computation to experiment. Physical Chemistry Chemical Physics 2013, 15, 3389-3412. 22. Xu, X.; Cheng, C.; Chowdhury, I. H., Molecular Dynamics Study of Phase Change Mechanisms During Femtosecond Laser Ablation. Journal of Heat Transfer 2004, 126, 727734. 23. Esfarjani, K.; Chen, G.; Stokes, H. T., Heat transport in silicon from first-principles calculations. Physical Review B 2011, 84, 085204. 24. Tong, M.; Liu, J.; Yu, X.; Dong, H. B.; Davidchack, R. L.; Dantzig, J.; Ceresoli, D.; Marzari, N.; Cocks, A.; Zhao, C.; Richardson, I.; Kidess, A.; Kleijn, C.; Hoglund, L.; Wen, S. W.; Barnett, R.; Browne, D. J., An integrated framework for multi-scale multi-physics numerical modelling of interface evolution in welding. IOP Conference Series: Materials Science and Engineering 2012, 33, 012029. 25. Wagner, G. J.; Jones, R. E.; Templeton, J. A.; Parks, M. L., An atomistic-to-continuum coupling method for heat transfer in solids. Computer Methods in Applied Mechanics and Engineering 2008, 197, 3351-3365.


34

Page 35 of 37

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


26. Pankavich, S.; Ortoleva, P., A Review of Two Multiscale Methods for the Simulation of Macromolecular Assemblies: Multiscale Perturbation and Multiscale Factorization. Computation 2015, 3, 29-57. 27. Singharoy, A.; Joshi, H.; Miao, Y.; Ortoleva, P. J., Space warping order parameters and symmetry: Application to multiscale simulation of macromolecular assemblies. J. Phys. Chem. B 2012, 116, 8423-8434. 28. Abi Mansour, A.; Ortoleva, P. J., Multiscale factorization method for simulating mesoscopic systems with atomic precision. J. Chem. Theory Comput. 2014, 10, 518-523. 29. Sereda, Y. V.; Espinosa-Duran, J. M.; Ortoleva, P. J., Energy transfer between a nanosystem and its host fluid: A multiscale factorization approach. J. Chem. Phys 2014, 140. 30. Cremons, D. R.; Plemmons, D. A.; Flannigan, D. J., Femtosecond electron imaging of defect-modulated phonon dynamics. Nature Communications 2016, 7, 11230. 31. Liu, M. B.; Liu, G. R., Smoothed particle hydrodynamics (SPH): An overview and recent developments. Arch. Comput. Meth. Eng. 2010, 17, 25-76. 32. Akin, J. E., Application and implementation of finite element methods. Academic Press: London, 1982. 33. Smith, A. F.; Weiner, R. G.; Skrabalak, S. E., Symmetry-Dependent Optical Properties of Stellated Nanocrystals. The Journal of Physical Chemistry C 2016, 120, 20563-20571. 34. Verlet, L., Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Physical Review 1967, 159, 98-103. 35. Birdsall, C. K.; Langdon, A. B., Plasma Physics via Computer Simulations. McGraw-Hill: 1985. 36. Trotter, H. F., On the product of semi-groups of operators. Proceedings of the American Mathematical Society 1959, 10, 545-551. 37. Hall, B., Lie Groups, Lie Algebras, and Representations: An Elementary Introduction. Springer: 2004; p 354. 38. Strang, G., On the Construction and Comparison of Difference Schemes. SIAM Journal on Numerical Analysis 1968, 5, 506-517. 39. Lax, P. D.; Wendroff, B., Difference schemes for hyperbolic equations with high order of accuracy. Communications on Pure and Applied Mathematics 1964, 17, 381-398. 40. Courant, R.; Friedrichs, K.; Lewy, H., On the partial difference equations of mathematical physics. IBM J. Res. Dev. 1967, 11, 215-234. 41. Crowley, W. P., Second-order numerical advection. Journal of Computational Physics 1967, 1, 471-484. 42. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79. 43. Foloppe, N.; MacKerell, J. A. D., All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. Journal of Computational Chemistry 2000, 21, 86-104. 44. Lopes, P. E. M.; Murashov, V.; Tazi, M.; Demchuk, E.; MacKerell, A. D., Development of an Empirical Force Field for Silica. Application to the Quartz−Water Interface. The Journal of Physical Chemistry B 2006, 110, 2782-2792. 45. Plimpton, S., Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of Computational Physics 1995, 117, 1-19. 46. Hoover, W. G., Canonical dynamics: Equilibrium phase-space distributions. Physical Review A 1985, 31, 1695-1697.


35


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 36 of 37

47. Fennell, C. J.; Gezelter, J. D., Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics. The Journal of Chemical Physics 2006, 124, 234104. 48. Hockney, R. W.; Eastwood, J. W., Computer Simulation Using Particles. 1989. 49. Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O., MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. Journal of Computational Chemistry 2011, 32, 2319-2327.


36

Page 37 of 37

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


37

Multiscale Molecular Dynamics Approach to Energy Transfer in

Recommend Documents