Optimal prediction in molecular dynamics
Optimal prediction approximates the average solution of a large system of ordinary differential equations by a smaller system. We present how optimal prediction can be applied to a typical problem in the field of molecular dynamics, in order to reduc…
Authors: Benjamin Seibold
Optimal Predition in Moleular Dynamis Benjamin Seib old Departmen t of Mathematis T e hnial Univ ersit y of Kaiserslautern Gottlieb-Daimler-Straÿe, 67653 Kaiserslautern, German y Email: seiboldmathematik.uni-kl. de Abstrat Optimal predition appro ximates the a v erage solution of a large system of ordinary dieren tial equations b y a smaller system. W e presen t ho w optimal predition an b e applied to a t ypial problem in the eld of moleular dynamis, in order to redue the n um b er of partiles to b e tra k ed in the omputations. W e onsider a mo del problem, whi h desrib es a surfae oating pro ess, and sho w ho w asymptoti metho ds an b e emplo y ed to appro ximate the high dimensional onditional exp etations, whi h arise in optimal predition. The th us deriv ed smaller system is ompared to the original system in terms of statistial quan tities, su h as diusion onstan ts. The omparison is arried out b y Mon te-Carlo sim ulations, and it is sho wn under whi h onditions optimal predition yields a v alid appro ximation to the original system. Keyw ords: Optimal predition, moleular dynamis, surfae oating, hopping, Laplae's metho d, lo w temp erature asymptotis, Mon te-Carlo 1 In tro dution Computations in the eld of moleular dynamis t ypially require a large omputational eort due to t w o fators: 1. Small time steps are required to resolv e the fast atomi osillations. 2. Large systems are obtained due to the large amoun t of atoms whi h ha v e to b e omputed. A wide v ariet y of metho ds has b een dev elop ed to remedy these problems. Larger time steps are admitted e.g. b y smo othing algorithms, whi h a v erage in time o v er the fast osillations. V arious other metho ds redue the degrees of freedom, e.g. m ultip ole metho ds [14℄ in the on text of long range partile in terations. In this pap er w e in v estigate whether the metho d of optimal predition, as presen ted and analyzed in [1, 4 11, 15, 19℄, an in priniple b e applied to problems in moleular dynamis in order to redue the n um b er of atoms to b e omputed. As a rst step in this in v estigation w e onne to a one dimensional mo del problem whi h inherits partiular prop erties from a real moleular dynamis problem. In Setion 2 w e presen t the real problem as it arises in appliations, and deriv e the simplied mo del problem. The onsidered problem is Hamiltonian, hene in Setion 3 w e presen t the metho d of optimal predition in the sp eial ase of Hamiltonian systems. In Setion 4 1 2 PR OBLEM DESCRIPTION 2 w e apply the metho d of optimal predition to the mo del problem. This yields expressions in v olving high dimensional in tegrals. W e ev aluate these in tegrals b y asymptoti metho ds, emplo ying the fat that the pro ess is running at a lo w temp erature, whi h yields a new and smaller system. Setion 5 deals with the n umerial sp eed-up. In Setion 6 w e dene riteria, ho w to in v estigate whether the new system is a v alid appro ximation to the original system. W e sho w ho w imp ortan t statistial quan tities, su h as diusion onstan ts and energy utuations, an b e obtained b y n umerial exp erimen ts. These are presen ted in Setion 7, and w e in v estigate under whi h onditions optimal predition preserv es the relev an t statistial quan tities. 2 Problem Desription 2.1 The Ph ysial Problem In the pro dution of semiondutors a thin la y er (a few atomi monola y ers) of opp er has to b e oated (sputtered) on to a silion rystal. A te hnial desription of the pro ess an b e found in [18℄. Imp ortan t for the qualit y of the pro dut is that opp er atoms m ust not p enetrate to o deeply in to the silion rystal. The opp er atoms p enetrate rstly b y their kineti energy when hitting the rystal surfae, seondly b y the pro ess of atomi hopping , whi h will b e desrib ed in the follo wing. In order to obtain sp ei kno wledge ab out these pro esses, moleular dynamis sim ulations ha v e to b e arried out, as desrib ed in [13℄. One imp ortan t asp et of the oating pro ess is that the system is out of its thermo- dynamial equilibrium only for v ery short times, namely for ab out 10 − 11 seonds after one opp er atom has hit the surfae of the silion rystal. During this time the opp er atom p enetrates in to the rystal and soni w a v es transp ort the impat energy a w a y to the b ottom of the rystal, whi h is onstan tly b eing o oled. Hene, after 10 − 11 seonds the whole rystal is in equilibrium again. On the other hand, the time b et w een t w o opp er atoms hitting the rystal surfae is on a time sale of 10 − 4 seonds, i.e. the system is in equilibrium nearly all the time, in partiular the temp erature is onstan t with resp et to spae and time. Ho w ev er, ev en in the state of thermo dynamial equilibrium single opp er atoms an hange their p osition in the silion rystal b y hopping ev en ts, i.e. a opp er atom gains b y aiden t enough energy to o v erome the p oten tial barrier b et w een t w o la y ers in the silion rystal and th us hops to a neigh b oring ell. By atomi hopping opp er atoms an p enetrate m u h deep er in to the silion rystal as their impat energy w ould allo w, hene the pro ess in equilibrium annot b e omitted from the omputation. The a v erage time b et w een t w o hopping ev en ts is on a time sale of 10 − 10 seonds, while the fast atomi osillations happ en on a time sale of 10 − 14 seonds. In this pap er w e sho w ho w the metho d of optimal predition an b e applied to the system in equilibrium. Only the atoms at the top la y ers of the rystal, where opp er atoms an b e found, are omputed exatly , while the silion atoms in the lo w er la y ers are k ept tra k of only in an a v erage sense. In order to k eep the te hnial diulties at a minim um, w e set up a one dimensional mo del problem whi h sim ulates atomi hopping. 2.2 The Mo del Problem In the mo del problem, w e assume t w o ma jor simpliations: 2 PR OBLEM DESCRIPTION 3 • F o us on a one dimensional problem, i.e. w e onsider n atoms lined up lik e b eads on a ord. A single opp er atom is inserted in to a line of n − 1 silion atoms. • The p oten tial V ( q ) dep ends only on the pairwise distanes of the atoms, i.e. V ( q 1 , . . . , q n ) = n X i,j =1 i k 6 = 0 else (4.21) In blo k form it an b e written as ∂ 2 V ∂ ˆ q ∂ ˜ q ( ˆ q , r ( ˆ q )) = 0 B 12 0 0 , (4.22) where B 12 is an upp er triangular k × k matrix. Substituting these sp eial v etors and matries in to the equation for the virtual atoms in (4.17 ) yields the follo wing equation of motion in blo k form A 11 A 12 A 21 A 22 · ˙ r V ˙ r m + k e = 0 B 12 0 0 · M − 1 · ˆ p , (4.23) 5 COMPUT A TIONAL SPEED-UP 11 whi h implies the k -dimensional system for r V ¯ A 11 · ˙ r V = B 12 · M − 1 l · ˆ p l . (4.24) Here ¯ A 11 = A 11 + ( A 12 · e ) · e T k , where e k = (0 , . . . , 0 , 1) T . The momen ta of the last k real atoms are denoted b y ˆ p l = ( p m − k + 1 , . . . , p m ) T , and M l is the diagonal matrix on taining the orresp onding masses. This relation an b e in terpreted as a b oundary la y er ondition whi h ats as if the rystal of silion atoms w as on tin ued to innit y , although it is atually ut o after the m -th atom. With the ab o v e mo diations the zero temp erature limit optimal predition system (4.17 ) b eomes an expliit (2 m + k ) -dimensional system of equations. Hene, one an exp et this system to yield a omputational sp eed-up, dep ending on the v alues of n , m and k . The question of sp eed-up will b e onsidered in Setion 5. Due to the v arious appro ximations in a hieving the smaller system it is not at all lear whether the smaller system yields the same dynamis as the original system. In Setion 6 w e will ompare the new system to the original system and in v estigate under whi h onditions the new system reets the orret dynamis, and under whi h onditions it do es not. 5 Computational Sp eed-Up Besides the ab o v e ph ysial in terpretations, the atual in ten tion w as to use optimal pre- dition as a metho d to redue the omputational eort. In this setion, w e onsider b oth the v ersion with l virtual atoms and the b oundary la y er v ersion. W e ompare the CPU times for omputations of the t w o optimal predition systems with the CPU time for the orresp onding omputations of the original system. The omparison is arried out in de- p endene of the sizes n and m . Sine the original v ersion of optimal predition do es only replae real atoms b y virtual ones, while on the other hand the b oundary la y er ondition v ersion allo ws to really omit atoms, a signian t sp eed-up an only b e exp eted from the b oundary la y er ondition v ersion. The omputations w ere p erformed on a net w ork of AMD A thlon-6 1.4 GHz omputers. The b oundary la y er ondition v ersion w as omputed with k = 10 virtual atoms. Figure 3 sho ws the sp eed-up fators with resp et to the original system for the t w o v ersions of opti- mal predition. The original v ersion of optimal predition do es not yield an y aeleration (the sp eed-up fators are less than 1). Apparen tly , for our mo del problem setting up and solving the linear system in (4.17 ) is more exp ensiv e than omputing the full system of equations. The b oundary la y er ondition v ersion, on the other hand, yields signian t sp eed-up fators for small m . In priniple, one an a hiev e arbitrarily high sp eed-up fa- tors b y k eeping m and k xed and inreasing n . In most ases, ho w ev er, the original rystal size n is giv en a priori, and suitable v alues for m are giv en b y the requiremen t that the new system m ust ha v e the same dynamis as the original one, as disussed in Setion 6. A hoie of m whi h is reasonable in man y ases, is utting the whole rystal of n atoms in to t w o halv es, i.e. m = ⌊ n 2 ⌋ . The sp eed-up fators for the t w o optimal predition systems are sho wn in Figure 4. While the original v ersion do es not redue the omputational eort, the b oundary la y er ondition v ersion yields a go o d sp eed-up for larger n . Care should b e tak en when generalizing these results. In our mo del problem the p oten tial funtions are fairly heap to ev aluate. In more ompliated systems it ould 6 COMP ARING TW O MOLECULAR D YNAMICS SYSTEMS 12 0 20 40 60 0 20 40 60 0 0.5 1 m optimal prediction speed−up n speed−up factor 0 20 40 60 0 20 40 60 0 5 10 m boundary layer condition speed−up n speed−up factor Figure 3: Sp eed-up fators dep ending on n and m . 0 50 100 150 0 1 2 3 4 5 6 7 8 9 10 n = 2m speed−up factor optimal prediction boundary layer condition Figure 4: Sp eed-up fators for n = 2 m v ery w ell pa y o to solv e the linear system in ( 4.17 ) instead, or to solv e the minimization problem diretly . On the other hand, the b oundary la y er ondition v ersion of optimal predition ould p ossibly fail to w ork in other appliations, e.g. in three spae dimensions. Su h questions will ha v e to b e in v estigated when applying optimal predition to a new problem. 6 Comparing T w o Moleular Dynamis Systems In moleular dynamis tra jetories in high dimensional phase spae are no appropriate means for omparing t w o systems, sine initial p ositions and momen ta an nev er b e kno wn exatly , and moleular dynamis is t ypially haoti. Instead, omparing means to test whether b oth systems ha v e similar dynamis. This is represen ted b y statistial quan tities su h as time orrelation funtions, diusion onstan ts, utuations of energy , et. W e onsider the follo wing statistial pro esses in order to ompare the t w o systems: • The distribution of the p osition the opp er atom. A opp er atom, whi h is initially lo ated alw a ys at the same p osition, is in the ensem ble of man y exp erimen ts desrib ed b y a diusion pro ess, whose distribution an b e appro ximated b y Mon te- Carlo sampling. • The n um b er of hopping ev en ts up to a xed time. 6 COMP ARING TW O MOLECULAR D YNAMICS SYSTEMS 13 • The utuation of energy of the rst m atoms. The energy of the rst m atoms utuates around some xed a v erage. W e onsider the v ariane of the energy o v er a xed time in terv al. The rst t w o quan tities are related to the diusion of a single opp er atom in the silion rystal due to hopping ev en ts. Sine opp er hopping has b een the eet of in terest in the rst plae, it is a natural riterion for omparison. In moleular dynamis, statistial quan tities of ergo di systems an b e omputed b y long time a v eraging or Mon te-Carlo sampling. Long time a v eraging means running a single omputation and using limiting pro esses to appro ximate statistial quan tities. Im- p ortan t examples, e.g. for appro ximating self-diusion onstan ts [2℄, are the Gr e en-Kub o formula [17, 23℄ or the Einstein r elation , whi h b oth appro ximate ensem ble a v erages for ergo di systems b y long time a v eraging. In our appliation, ho w ev er, long-time omputa- tions are problemati, sine rstly the opp er atom ma y tra v el to the b oundaries of the silion rystal and seondly soni w a v es will ome in eet, as sho wn in Setion 7 . Hene, w e obtain the diusion onstan t b y Mon te-Carlo sampling instead, i.e. w e solv e the same system o v er and o v er again with short omputation times, where the initial onditions are sampled from the anonial measure. In other w ords: W e obtain the a v eraged quan tities not b y long-time a v eraging, but b y a v eraging o v er man y samples. Sampling b oth the initial p ositions q i and momen ta p i from the anonial measure Z − 1 e − H ( q,p ) w ould require exp ensiv e metho ds lik e a eptan e-r eje tion metho ds or Metr op o- lis sampling [12, 16℄ due to the struture of the p oten tial V ( q ) . W e irum v en t su h prob- lems b y setting the initial p ositions q i in to the p oten tial minim um and sampling the initial momen ta p i from ˜ Z − 1 e − T ( p ) , whi h is just sampling indep enden tly from Gaussian distri- butions. In our sim ulations after ab out 5 · 10 − 14 s the Hamiltonian dynamis has driv en the system in to equilibrium. A dditionally , k eeping the initial p ositions xed automatially guaran tees to remain in the orret domain M L as giv en b y (4.2 ). 6.1 A Random W alk Mo del for the Copp er Diusion The onsidered opp er diusion is due to hopping ev en ts, while displaemen ts due to short osillations b et w een the same t w o silion atoms are not tak en in to aoun t. Hene, the or- resp onding diusion pro ess is disrete in spae on the spatial grid {− ν d 0 , . . . , 0 , . . . , ν d 0 } . Let us assume for the momen t, that the diusion pro ess is linear and hene desrib ed b y a (2 ν + 1) -dimensional omp artment mo del ˙ u ( t ) = A · u ( t ) , u (0) = e m +1 = (0 , . . . , 0 , 1 , 0 , . . . , 0) T , (6.1) where A ∈ R (2 ν +1) × (2 ν +1) has olumn sums equal to zero to ensure mass onserv ation, i.e. d dt P i u i ( t ) = 0 . The analytial solution to (6.1 ) is u ( t ) = exp( tA ) · e m +1 . (6.2) One example for su h a pro ess is giv en b y the tridiagonal T o eplitz (up to the b oundary en tries) matrix with the stenil κ d 2 0 (1 , − 2 , 1) , whi h is a nite dierene appro ximation to the heat equation. Hene, the ompartmen t mo del (6.1 ) with this matrix on v erges to the heat equation as d 0 → 0 . Still, there are man y other matries A , whose disrete diusion pro esses all on v erge to the heat equation as d 0 → 0 . 6 COMP ARING TW O MOLECULAR D YNAMICS SYSTEMS 14 The hopping b eha vior of the opp er atom an b e mo deled as a sp ei random w alk. One imp ortan t asp et of the opp er hopping is that hopping ev en ts are orr elate d , in the sense that giv en the opp er atom has just hopp ed, it is m u h more lik ely than normally that a seond hopping ev en t to the same diretion follo ws, sine the kineti energy is not ompletely lost in one single jump. Hene, the opp er diusion is formally non-Mark o vian. Still, the hopping an b e desrib ed as a Mark o vian random w alk b y assuming hopping ev en ts to b e unorrelated, but to allo w the opp er atom to hop o v er more than one silion atom in one single hopping ev en t. These assumptions lead to a mo del whi h is t ypial in the on text of sto hasti pro- esses, see [22℄. Let X t denote the p osition of the opp er atom at time t . F or the n th hopping ev en t let T n b e the hopping time, ∆ T n = T n − T n − 1 the time sine the previous hopping ev en t, and ∆ n ∈ Z \ { 0 } the n um b er of silion atoms whi h the opp er atoms hops o v er to the righ t. Here ∆ n < 0 means hopping to the left. Assume no w that b oth the ∆ T n and the ∆ n are indep enden t and distributed aording to P (∆ T n ∈ [ s, s + d s )) = α · exp( − αs ) d s, (6.3) P ( | ∆ n | = i ) = p i . (6.4) Here α is a parameter on trolling the hopping rate (see [21, 22℄ for a deriv ation), and ( p i ) i ∈ N is a non-negativ e sequene satisfying P i ∈ N p i = 1 2 and the onstan t γ = P k i =1 ( d 0 i ) 2 p i is nite. F or our mo del problem w e simply assume p = ( p 1 , . . . , p k ) . A dditionally , w e restrit to symmetri random w alks, whi h is reasonable as long as the opp er atom do es not approa h the rystal b oundaries. An analysis of the desrib ed random w alk, follo wing the analysis in [22℄, yields that the v ariane of the pro ess at time t equals V ar ( X t ) = 2 γ αt . A omparison with the v ariane giv en b y the heat equation yields the relation κ = γ α, (6.5) whi h allo ws us to sp eak onsisten tly of a diusion onstan t κ also for the disrete diusion pro ess (6.1). As deriv ed in [22℄ the probabilit y distribution of the random w alk X t is desrib ed b y (6.1), where the so alled innitesimal gener ator A is a band matrix with the stenil α · ( p k , . . . , p 1 , − 1 , p 1 , . . . , p k ) . Dene ˜ A as a orresp onding band matrix with the stenil γ − 1 · ( p k , . . . , p 1 , − 1 , p 1 , . . . , p k ) , hene A = κ ˜ A . Both matries ha v e to b e hanged in the upp er and lo w er ro ws aording to the appropriate b oundary onditions. Assume no w, that the v alues p 1 , . . . , p k are kno wn. Hene, the matrix ˜ A is ompletely determined. Let v ( t ) ∈ R 2 ν +1 denote the n umerially obtained distribution v etors for the p osition of the opp er atom for all times t . Initially v (0) = e m +1 , as in (6.1 ). Sine it is unlear, whether the diusion parameter κ is onstan t in time, w e let it b e a time dep enden t parameter κ ( t ) , whi h w e ompute at times t 1 , . . . , t µ . F or obtaining the v alue κ j = κ ( t j ) w e onsider the diusion on the time in terv al I j = [ t j , t j ] , where t j = t j − ∆ t 2 and t j = t j + ∆ t 2 . Giv en ∆ t is not to o large, w e an appro ximately assume the diusion parameter to b e onstan t on the giv en in terv al: κ ( τ ) = κ j ∀ τ ∈ I j . On I j w e use the data pro vided b y the real ev olution v ( t ) to ompute an L 2 -t of the diusion parameter κ j with resp et to error funtional F ( κ j ) = Z t j t j k e ( τ − t j ) κ ( t j ) ˜ A · v ( t j ) − v ( τ ) k 2 2 d τ , (6.6) 7 NUMERICAL EXPERIMENTS 15 whi h is partiularly stable with resp et to errors in v due to the Mon te-Carlo sampling. The time ∆ t m ust on the one hand b e small enough to justify the appro ximation that the diusion parameter is onstan t on the in terv als I j , on the other hand it m ust b e large enough to ha v e already some diusion tak en plae. In Subsetion 7.1 w e ompute diusion parameters for the n umerial data obtained for our mo del problem. 7 Numerial Exp erimen ts W e sim ulate a rystal of n = 70 atoms with 69 silion atoms and 1 opp er atom, whi h starts at the 22 nd p osition. F or the optimal predition system, w e ho ose m = 50 . These are enough atoms su h that a v eraged quan tities mak e sense and one an sp eak of a ther- mo dynamial equilibrium. A dditionally , the rystal has a reasonable in terior region whi h is not aeted b y an y b oundary eets. The in tegration is p erformed b y the lassial ex- pliit fourth order Runge-Kutta metho d, whi h w e prefer o v er energy preserving metho ds in this on text, as it allo ws signian tly larger time steps, while still resolving the hopping ev en ts aurately . The in tegration time step is ∆ t = 2 . 5 · 10 − 15 s, and the omputation time is t ∗ = 4 . 0 · 10 − 13 s, whi h is short enough, su h that the hange in total energy due to the in tegrator is insignian t. W e solv e the system N = 25000 times at a temp erature T = 7000 K , with the initial data sampled as desrib ed in Setion 6. F or our exp erimen ts, Mon te-Carlo estimates yield an exp eted error in the distribution of ab out 4 · 10 − 3 , whi h is signian tly smaller than the dierene of an y t w o quan tities in omparison. Note that the short omputation time t ∗ allo ws to exlude an y eets aused b y soni w a v es tra v eling through the rystal. The v elo it y of sound an b e estimated as deriv ed in [17℄: The equilibrium distane inside the silion rystal is appro ximately d 0 = 1 . 87 · 10 − 10 m . Using the Y oung mo dulus Y = d 0 · ∂ 2 H lo ∂ d 2 0 ( d 0 ) = 5 . 05 · 10 − 9 kg · m s 2 w e obtain a v elo it y of sound of c = q Y · d 0 m S i = 4 . 50 · 10 3 m s . The shortest distane of the opp er atom to an y b oundary in the rystal is s left = q 22 − q 1 = 2 . 01 · 10 − 9 m in the original system and ˜ s righ t = q 50 − q 22 = 1 . 90 · 10 − 9 m in the optimal predition system. Hene, the minim um time a soni w a v e tak es to tra v el from a b oundary to the opp er atom is appro ximately t min = ˜ s righ t /c = 4 . 22 · 10 − 13 s, whi h is longer than the omputation time t ∗ = 4 . 0 · 10 − 13 s. In Subsetion 7.4 w e will deal with the ase when soni w a v es are presen t. Figure 5 sho ws the distribution of the p osition of the opp er atom in the silion rystal when solving the original system. The analogous distribution for the optimal predition system lo oks indistinguishably similar in this kind of plot. The pro ess is apparen tly of a diusiv e nature, but the diusion parameter hanges with time, as already the ev olution of the maxim um indiates. The follo wing analysis will onrm this observ ation. The time-dep enden t relativ e error b et w een the distributions for the original and the optimal predition system e ( t ) = max x | v ( x, t ) − ˜ v ( x, t ) | max x | v ( x, t ) | (7.1) is less than 3% up to the time t = 3 . 0 · 10 − 13 s and still less than 9% o v er the omplete time in terv al, whi h is not o v erwhelmingly small, but do es indiate denite similarities b et w een the t w o distributions. 7 NUMERICAL EXPERIMENTS 16 0 1 2 3 4 x 10 −13 15 20 25 30 0 0.2 0.4 0.6 0.8 1 position in crystal t Figure 5: Distribution of the opp er atom's p osition 7.1 Diusion P arameters W e ompute the diusion parameters κ ( t ) for the original and the optimal predition system, using the metho d desrib ed in Subsetion 6.1 . The parameters ( p 1 , . . . , p k ) of the orresp onding random w alk mo del are obtained b y Mon te-Carlo sampling. W e use the exp erimen ts for the original system, whi h already yielded the distribution sho wn in Figure 5. In ea h Mon te-Carlo exp erimen t w e follo w the path of the opp er atom and luster onseutiv e hopping ev en ts whi h happ en in an in terv al of ∆ t 1 = 6 . 0 · 10 − 14 s to a single one. Prior to this lustering, fast double hopping ev en ts in to opp osing diretions, i.e. those happ ening inside of ∆ t 2 = 2 . 0 · 10 − 14 s , whi h are solely due to osillations of silion atoms, m ust b e exluded. The times ∆ t 1 and ∆ t 2 are suitably hosen for our mo del problem. The results of the ab o v e desrib ed Mon te-Carlo exp erimen ts are sho wn in Figure 6. The hopping probabilities are giv en b y the b o x histogram. Note that due to Mon te-Carlo errors and b oundary eets the resulting v alues are not exatly symmetri. Sine w e wish to onsider a symmetri random w alk, w e ho ose as p 1 , . . . , p k the a v erage v alues of the obtained results. Here, only p 1 , . . . , p 11 are greater than zero. Hene, w e ho ose k = 11 for the omputation of the diusion parameters. The urv e denotes the v alues i 2 p i (saled to t in to the same plot), whi h are relev an t, sine the v ariane P k i =1 i 2 p i = γ d 2 0 is prop ortional to κ (see relation (6.5)). Figure 7 sho ws the time-dep enden t diusion parameter κ ( t ) for the original system and the optimal predition appro ximation, omputed as desrib ed in Subsetion 6.1 . In b oth ases w e ho ose ∆ t = 2 . 5 · 10 − 14 s . Imp ortan t observ ations and estimates are: • The diusion parameters start with a strong p eak, and after 1 · 10 − 13 s they utuate around a xed v alue of κ = 4 . 4 · 10 − 10 m 2 s . This b eha vior is to o pronouned to b e only due to Mon te-Carlo and appro ximation errors. The initial b eha vior is most lik ely due to the fat that the system do es not exatly start in its thermo dynamial equilibrium. 7 NUMERICAL EXPERIMENTS 17 −10 −5 0 5 10 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 hopping width probability p i p i i 2 Figure 6: Probabilities of hopping o v er more than a single atom 0 0.5 1 1.5 2 2.5 3 3.5 4 x 10 −13 0 1 2 3 4 5 6 7 8 x 10 −8 t / s κ / m 2 /s original system optimal prediction Figure 7: Diusion parameters o v er time • Considering that κ ( t ) sho ws su h pronouned b eha vior, the t w o funtions κ ( t ) for the original system, and ˜ κ ( t ) for the optimal predition appro ximation are remark ably lose to ea h other. Up to the time t = 1 . 4 · 10 − 13 s the the distane b et w een the t w o urv es is quite small. After that time the urv es dier more, but sho w the same features. After 3 . 0 · 10 − 13 s the error tak es its maxim um, whi h oinides with the errors observ ed in the distributions sho wn in Figure 5 . • In order to roughly estimate whether the obtained diusion parameter is reasonable, w e ompare it with a diusion o eien t measured in a real material. In [3℄ the diusion onstan t of opp er in a silion rystal at a temp erature of T 0 = 1273 K is giv en as κ 0 = 4 . 4 · 10 − 10 m 2 s . Assuming that the diusion onstan t dep ends linearly b oth on temp erature and the p oten tial barrier, as e.g. with the Einstein form ula [17℄, w e obtain an estimate κ ≈ κ 0 · T T 0 · ∆ E 0 ∆ E = 1 . 7 · 10 − 8 m 2 s for T = 7000 K , using the v alues ∆ E 0 = 3 e V and ∆ E = 0 . 43 e V . A lo ok at Figure 7 indiates that our exp erimen tally obtained diusion parameters are indeed in this region. 7.2 The Num b er of Hopping Ev en ts Besides diusion parameters, also the pure n um b er of hopping ev en ts whi h happ en up to a giv en time t ∗ should b e preserv ed b y optimal predition. Figure 8 sho ws the n um b er of hopping ev en ts for the ab o v e omputation plotted as histograms. The solid line represen ts the original system, and the dashed line stands for the optimal predition system. The t w o graphs dier only for the probabilities of zero and one hopping ev en t. Apart from that, one an sp eak of the same hopping b eha vior. 7 NUMERICAL EXPERIMENTS 18 0 2 4 6 8 10 12 14 16 18 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 number of hopping events N fraction of experiments with N hopping events original system optimal prediction Figure 8: Num b er of hopping ev en ts 0 0.5 1 1.5 2 x 10 −51 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 energy variance / kg 2 m 4 /s 3 fraction of experiments original system optimal prediction Figure 9: T otal energy utuation 7.3 Energy Flutuations While the total energy is onstan t for the original system as w ell as for the optimal predition appro ximation, the energy of the rst m atoms E left ( t ) = 1 2 m X i =1 p 2 i ( t ) m i + m X i,j =1 i
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment