Kinetic Monte Carlo simulation of shape transition of strained quantum dots中文翻译
蒙特卡罗法
蒙特卡罗法简单介绍和案例蒙特卡罗法历史悠久。
1773年法国G.-L.L.von 布丰曾通过随机投针试验来确定圆周率π的近似值，这就是应用这个方法的最早例子。
蒙特卡罗是摩纳哥著名赌城，1945年 J.von 诺伊曼等人用它来命名此法，沿用至今。
数字计算机的发展为大规模的随机试验提供了有效工具，遂使蒙特卡罗法得到广泛应用。
在连续系统和离散事件系统的仿真中，通常构造一个和系统特性相近似的概率模型，并对它进行随机试验，因此蒙特卡罗法也是系统仿真方法之一。
对于蒙特卡罗技术应用于不可预见费的估算的研究，是对蒙特卡罗技术应用的拓展，能更好地了解尝试其在项目管理方面更多的应用，用其解决项目管理的问题。
用蒙特卡罗技术研究不可预见费，尝试用蒙特卡罗解决一般项目的不可预见费求取问题，避免不可预见费过高过低的问题。
蒙特卡洛方法的基本思想是：将符合一定概率分布的大量随机数作为参数带入数学模型，求出所关注变量的概率分布，从而了解不同参数对目标变量的综合影响以及目标变量最终结果的统计特性。
蒙特卡洛方法的基本原理简单描述如下：假定函数),...,,(21nx x x f y =，蒙特卡洛方法利用一个随机数发生器通过抽样取出每一组随机变量 (ni i i x x x ,...,,21)，然后按),...,,(21n x x x f y =的关系式确定函数的值),...,,(21ni i i i x x x f y =。
反复独立抽样（模拟）多次(ｉ=1，2，…)，便可得到函数的一组抽样数据(n y y y ,...,,21)，当模拟次数足够多时，便可给出与实际情况相近的函数y 的概率分布与其数字特征。
蒙特卡罗法（Monte Carlo Simulation ）也称随机模拟，它主要依据概率分布对随机变量进行抽样，然后将样本带入数学模型进行计算得到应变量。
虽然蒙特卡罗模拟技术只给出的是统计估计而非精确的结果且应用其研究问题需要花费大量的计算时间，但它对问题的维数不敏感，对求解对象是线性问题与否也没有原则性要求，因此在复杂系统的不确定分析中，蒙特卡罗方法成为不可或缺的手段。
蒙特卡罗算法
xi 1 6xi mod11, ui+1 xi 1 /11 （a 6, m 11 ）
如果令 a 3, x0 1 ,得到序列: 1,3,9,5, 4,1,3,9........ 如果令 a 3, x0 2, 得到序列： 2, 6, 7,10,8, 2, 6.......
14
3. 指数分布 指数分布随机变量X的概率密度为
e x , x 0 f ( x) 0, x 0
指数分布常用来描述寿命问题.
15
二.随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数。 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudorandom numbers)。 3.在模拟中，我们需要产生各种概率分布的随机数， 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数。
2
基本思想很早以前就被人们所发现和 利用。17世纪，人们就知道用事件发 生的“频率”来决定事件的“概率”。 19世纪人们用投针试验的方法来决定π。 高速计算机的出现，使得用数学方法 在计算机上大量模拟这样的试验成为 可能。
3
从Buffon(蒲丰)投针问题谈起
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为 l ( l < a) 的 针,此针与任一平行线相交的概率 p。
Monte Carlo Simulation Methods (蒙特卡罗模拟方法)
主要内容: 一. M-C方法概述. 二. 随机数的生成. 三. 模拟训练. 四. 实验题目.
成信院数学与信息科学系 李胜坤
MonteCarlo方法
2. 频率检验
检验每组观测频数 ni与理论频数mi ＝ N １/k之间相差的显著性
Monte Carlo 方法解决实际问题的过程中 ， 主要有以下几个内容 ① 建立简单而又便于实现的概率统计模型 ， 使所求的解正是该模型的某 一事件的概率或数学期望 ， 或该模型能够直接描述实际的随机过程。
② 根据概率统计模型的特点和计算的需求 ， 改进模型 ， 以便减小方差和 减低费用 ， 提高计算效率。
选择递推函数必须注意以下几点 ：
① 随机性好 ； ② 在计算机上容易实现 ； ③ 省时 ； ④ 伪随机数的周期长。
2 伪随机数的产生方法
最基本的伪随机数是均匀分布的伪随机数。 该方法是首先给一个 ２r位的数 ， 取其中间的 r位数码作为第一个伪随机数 ， 然后将这个数平方 ， 构成一个新的 ２ r位的数 ， 再取中间的 r位数作为第 二个伪随机数。 如此循环可得到一个伪随机数序列。 该方法的递推公式为
提高精度一位数 ， 抽样次数要增加１００ 倍 ； 减小随机变量的标准 差 ， 可以减小误差 。
Monte Carlo 方法具有以下四个重要特征 ：
① 由于 Monte Carlo 方法是通过大量简单的重复抽样来实现的 ， 因 此 ， 方法和程序的结构十分简单 。
② 收敛速度比较慢 ， 因此 ， 较适用于求解精度要求不高的问题。
具体试验步骤如下 ：
（1） 产生服从给定分布函数 g（x）的随机变量值 xi （2）检查 xi是否落入积分区域（a≤ x ≤ b） ， 如果满足条件 ， 则记录一次。
Kinetic Monte Carlo simulations of heteroepitaxial growth with an atomistic model
Kinetic Monte Carlo simulations of heteroepitaxial growth with an atomistic model of elasticityPinku Nath,Madhav Ranganathan ⁎Department of Chemistry,Indian Institute of Technology Kanpur,Kanpur 208016,Indiaa b s t r a c ta r t i c l e i n f o Article history:Received 7February 2012Accepted 21May 2012Available online 28May 2012Keywords:Heteroepitaxial systemsKinetic Monte Carlo simulations Submonolayer growth Atomic elasticityWe have implemented Kinetic Monte Carlo (KMC)simulations of growth of heteroepitaxial thin films.A sim-ple cubic Solid-on-Solid (SOS)model is used to describe the atomic con figurations and nearest neighbor bonds are used to describe the energetics.Elastic effects are modeled using harmonic springs between atoms displaced from their lattice positions.The mis fit strain is a consequence of different equilibrium spring lengths for the substrate and film.The consistency of this elastic model with continuum theories for strained surfaces has been shown by performing elastic energy calculations for various morphologies.KMC simula-tions for submonolayer deposition show scaling behavior in the island size distribution.The resulting island shapes are predominantly square and do not show any shape transitions in the physically relevant range of conditions.This method gives a detailed understanding of elastic interactions and their interplay with surface diffusion in heteroepitaxial systems.©2012Elsevier B.V.All rights reserved.1.IntroductionStrained heteroepitaxial systems like Germanium (Ge)on Silicon (Si)are characterized by a small mismatch between the lattice parameters of the two materials,the film and the substrate.This mis fit can affect the morphology of these films.For example,it is well known that a growing heteroepitaxial film develops spontaneous undulations that eventually lead to formation of mounds and,sometimes,quantum dots [1–3].Controlled growth processes like Molecular Beam Epitaxy (MBE)can be used to closely study these phenomena and also grow high quality films and nanostructures [4–15].Heteroepitaxial systems have been of interest for nanoscale optoelectronic devices [16–18];e.g.Field-Effect-Transistors made with Si/Ge [19],GaAsP/InGaAs quantum wells for solar cells [20].The nature of growth of a strained film is in fluenced by the mis fit.For Ge on Si(100),the mis fit is about 4%.The first few layers of the film grow flat over the substrate.As more layers are deposited,the strain accumulated in the film increases and this causes spontaneous mound formation in the film (Fig.2).This type of growth is called Stranski –Krastanov growth [1,21]and is also found in other heteroepitaxial systems like InAs/GaAs [22].The transition from layer-by-layer growth to mound formation after a certain number of layers is referred to as the 2D/3D growth transition.For systems with larger mis fit,it is expected that this transition will be observed at lower coverages.It has been reported that the wetting layer in Ge/Si(100)is about 4–5ML thick [23]where it is about 1–2ML thick forInAs/GaAs [24].The shapes and size distribution of the mounds have also been a subject of many experimental [2,25,7]and theoretical stud-ies [26–31].The behavior during the early stages of mound formation,when the coverage is less than 1ML,has been investigated experimentally [32–34]and theoretically [35–42].During this stage the surface is mostly covered with monolayer islands of different sizes.The shape and size distribution of these islands depend on the surface studied.For example,in the case of Ge/Si(100),these islands are mostly rectan-gular in shape.The relative role of mis fit strain and anisotropic surface energy in the shape and size of these islands is a question of interest to us.KMC simulations have proved to be very useful for the study of surface morphology during growth [43].KMC allows us to study the interplay between deposition and surface diffusion during growth of a certain surface.For heteroepitaxial systems,the elastic strain modi fies the surface diffusion rates.These elastic effects are inherently long range effects.Due to this reason,the resulting KMC scheme is no longer a local scheme.Strain effects have been included in KMC simulations in many different ways.One class of approaches involves calculation of the elastic response of the semi-in finite substrate medium to the forces caused by the morphological features on the surface,like steps,adatoms,etc.The elastic response is typically calculated using continuum approaches involving Green's functions [41,44–46].Alternatively,the effect of elastic strain can be modeled using a reduction in the diffusion barrier for atoms in the film [37],or as an effective interaction between surface atoms in the film [40].In this article we focus on KMC simula-tions in which elastic effects are included through an atomistic model of elasticity with harmonic springs [47–50].Our method is closelySurface Science 606(2012)1450–1457⁎Corresponding author.Tel.:+915122596037;fax:+915122597436.E-mail address:madhavr@iitk.ac.in (M.Ranganathan).0039-6028/$–see front matter ©2012Elsevier B.V.All rights reserved.doi:10.1016/j.susc.2012.05.015Contents lists available at SciVerse ScienceDirectSurface Sciencej o u r na l h o me p a g e :ww w.e l s e v i e r.c o m /l o c a t e /s u s crelated to a prior work in this area,and we apply the method of calcu-lations on a realistic system.The elastic energy calculation involves a multidimensional minimization of the total energy of all the harmonic springs in the system.This makes the elastic computation extremely time-consuming.To implement this model in a KMC simulation,it is therefore necessary to come up with a set of approximations that will reduce the number of elasticity calculations.The rest of this article is organized as follows:In Section 2,we de-scribe the multicomponent Solid-on-Solid model that we use to de-scribe the surface.In Section 3,we show the KMC simulation method and the calculation of the surface hopping probabilities.In Section 4,we explain the model of elasticity used for the heteroepitaxial system.In Section 5,we illustrate how the KMC simulations can be performed with explicit elastic effects.In Section 6,we describe some results of our work,first on elastic energy calculations,and then submonolayer growth simulations.Finally,we make summarize the key findings and give our outlook on future work in this area in Section 7.2.Multicomponent Solid-on-Solid modelWe present here a lattice model for studying the evolution of the surface of a multicomponent system.Our model is based on the sim-ple cubic Solid-on-Solid (SOS)model,in which the crystal lattice is described by a 2D array of columns of various heights [51,43].For the single component SOS model,it is suf ficient to specify the heights of the columns in order to specify the surface.In the multicomponent case,each lattice site can be occupied by one of the various species as shown in Fig.1.In general,the behavior of the surface is dependent on the nature of the species in the columns.Hence it is essential in this model to have a variable to describe the type of the atom in each of the occupied lattice sites.Further,we impose a no-overhang condition,according to which a lattice site can only be occupied if the site one layer below is also occupied.The multicomponent SOS model can then be implemented in a Kinetic Monte Carlo (KMC)scheme in way similar to the standard SOS model using moves to describe the various processes of interest.The crystal is assumed to be growing in the +Z (upward)direction.The substrate is considered to in finitely thick and extend all the way to Z =−∞.We impose periodic boundary conditions in the X and Ydirections with the supercell dimension to be a multiple of the substrate lattice constant a s .Our simulations are designed to model the Ge on Si(100).We note that the systems we are considering do not have a simple cubic lattice,and the true Si(100)surface is reconstructed.The focus of this article is to study the interplay between elasticity and sur-face diffusion and we expect that the qualitative features of this behavior will be reproduced reasonably well in the simple cubic model.In the rest of this article,we will be using “Si ”and “Ge ”to refer to the substrate and film respectively.3.Kinetic Monte Carlo methodWe consider the dynamical processes taking place on the crystal surface during a typical MBE experiment,namely deposition and surface diffusion.We neglect bulk diffusion,defect formation and migration,and exchange of atoms in the bulk of the crystal.Surface diffusion is modeled by thermally activated hopping of surface atoms at random to their neighboring sites.The rate of this process can be written down in the Arrhenius form as r ¼ν0exp −E b =k B T ðÞð1Þwhere E b is energy barrier that has to be overcome in order for the atom to hop,k B is the Boltzmann constant,T is the absolute temperature and ν0is the attempt frequency.The attempt frequency ν0is independent of temperature and its value is around 1013s −1,corresponding to vibra-tional frequency of an atom in the lattice site.We assume the same value of ν0for all surface atoms in the system.The hopping rates satisfy the detailed balance condition,so thatr i →j ðÞ¼e −E j =k B T e −E i =k B ¼P eqj ðÞeq ð2Þwhere r (i →j )is the rate of transition from con figuration i to j involving a surface atom,and E i and E j are the energies of the two con figurations and P eq (i )and P eq (j )are their respective equilibriumprobabilities.Fig.1.Typical con figuration in a two-component simple cubic SOS model (left)and a cross-sectional view of the rectangular section (right).The black and white balls represent different atoms.For example,the black atoms can be Si atoms and the white atoms can beGe.Fig.2.A cartoon showing how strain caused by mismatch between the film (gray circles)and the substrate (black circles)is signi ficantly relaxed at surfaces,and at the edges of steps andmounds.Fig.3.The elastic model showing springs between NN and NNN atoms for a typical con-figuration.The springs in the film (gray)have different parameters than in the substrate (black).In a typical optimization to calculate the elastic energy,we consider numerically op-timize displacements in the top few layers of the substrate (two layers in the above figure),and solve for the displacement field analytically in the rest of the semi-in finite substrate;as indicated by the wavy lines.1451P.Nath,M.Ranganathan /Surface Science 606(2012)1450–1457Detailed balance ensures that when there is no deposition,the system will tend towards equilibrium with a Boltzmann distribution of energy.In the simulation of this model,we carry out a combination of random deposi-tion and hopping moves with the appropriate acceptance rates.The barrier energy E b is the total energy that a surface atom needs to have in order to hop to a neighboring site.In the absence of elastic energy,this barrier is the energy required to break all the bonds of the atom.The value of this barrier depends on the number and strength of these bonds,so E b =E bond =∑j γj where γi is the strength of the bond between the atom (which we attempt to hop)and the i th atom.The sum is performed over all the bonds of the surface atom.The strength of the bond depends on the particular model of interac-tion chosen.For example,if we are describing a Ge atom on the (100)surface of Si using our simple cubic lattice,we use a nearest neighbor (NN)interaction model consisting of two kinds of bonds with energies E h and E v corresponding to bonds in the horizontal and vertical direc-tions respectively.In the absence of elastic energy,the barrier energy is a local quantity which can be stored for all surface atoms.After every move,the barrier energy of a small number of atoms has to be updated.This makes it possible to implement the KMC using a rejection-free method similar to the Bortz –Kalos –Lebowitz algorithm [52].This is not possible when we include elastic effects.Hence,we implement a KMC scheme with ac-ceptance and rejection of atom hopping moves so as to satisfy Eq.(1).4.Elastic modelThe lattice mismatch of a heteroepitaxial film is de fined as =(a f −a s )/a s where a f and a s represent the lattice constants of film and the substrate respectively.For a film of Ge on Si,this value is around 4%.This mismatch can be controlled by depositing a mixture of Si and Ge on the film of varying concentrations.Following previous work [53,54],we implement a model of hetero-epitaxy using springs between atoms.This model is explained in the references but we give a brief review here.Associated with each atom is a 3-dimensional elastic displacement from its lattice position.Each atom is connected by springs to its nearest neighbors (NN)and next nearest neighbors (NNN)(Fig.3).The energy of each spring is given by E spring =(1/2)k (d −l eq )2where d is the distance between the atoms and l eq is the equilibrium length at which the spring energy is zero and k is the spring constant.The distance d between atoms on neigh-boring lattice sites is calculated geometrically using the lattice constants and the elastic displacements.The values of k and l eq depend on the nature of the bond (NN or NNN),and the species of atoms forming the bond.For a pure cubic system,the two independent elastic moduli,C 11and C 12,are related to the spring constants by the relations [55]k L ¼a C 11−2C 12ðÞk D ¼aC 12ð3Þwhere k L and k D are the spring constants for the NN and NNN bonds re-spectively,and a is the equilibrium length of the lateral springs in the simple cubic lattice model.The elastic energy leads to a restoring force in the system that tries to minimize the spring energy of each atom.In order to model heteroepitaxial systems,we assign different equilibrium spring lengths to the film and the substrate atoms corresponding to their lattice con-stants.The equilibrium bond lengths of the NN and the NNN bonds in the film are a f and a f ffiffiffi2p .We assumed that the spring constants k L and k D remain the same for different kinds of atoms,and that the equilibrium length of a bond between the substrate atom and the film atom is equal to that of a bond between two film atoms.We will discuss the implica-tions of these assumptions later.The cell periodicity ensures that for a pure substrate,there are no elas-tic effects.When a single complete monolayer is deposited on the sub-strate,the bonds are compressed in the horizontal (X and Y)directionsas shown in Fig.3.This causes the film to expand in the Z-direction.The equilibrium location of the atoms in the Z direction can be calculated using the minimum energy or zero force condition.Under this condition,the separation of atoms in the Z-direction,a L can be obtained by solving k L a L −a s ðÞþ4k D ðffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffia 2Lþa 2s q −a f ffiffiffi2p Þa L ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffia 2L þa 2sq ¼0:ð4ÞThe elastic energy is a function of the displacements of all the atoms.By minimizing this function with respect to the displacements,we cal-culate the optimum displacements.This is the elastic relaxation that has to be carried out at every step of the Kinetic Monte Carlo simulation.We can do this using any of the standard multidimensional minimization routines like conjugate gradient,steepest descent,etc.We are interested in growth of a thin film on a large substrate.For a semi-in finite substrate with given displacements of the topmost layer,the elastic displacements of the substrate are calculated analytically based on the procedure outlined in Ref.[47].In this procedure,we line-arize the springs in the substrate and expand the displacement vector field {U }in a Fourier basis in the X and Y directions.This Fourier trans-form,~U n o decays exponentially inside the substrate.We expand ~U n o using eigenvalues and eigenvectors and the coef ficients are calculatedusing the displacements of the topmost layer of the substrate.The dis-placement field {U }is used to calculate the contribution to the elastic energy from the substrate.This procedure allows us to greatly reduce the number of atoms involved in the elastic energy calculation.In practice,we treat a few layers of the substrate as discrete layers before treating the remaining as a continuum as shown in Fig.3.Since we treat a few layers of the substrate as discrete layers,we can potentially model intermixing wherein atoms from the substrate migrate into the film.This intermixing is believed to be one of the mechanisms for strain relaxation in these systems.We end this section with a few comments about the elastic model and the KMC simulations.Firstly,our KMC simulations are implemented on a Simple cubic SOS model.We allow for atoms to hop from one site to one of their neighboring sites.The elastic displacements are localized on the lattice sites.These displacements are used in the calculation of the distance between two atoms on neighboring lattice sites.It is,therefore,implicit that these displacements are not too large (not greater than half the spacing between atoms).Secondly,there are some parameters in the model like the spring constants and the equilibrium lengths for all the different types of bonds.The equilibrium lengths of bonds involving only one kind of atom are based on the bulk lattice constants of the atom.But there is some arbitrariness in the choice of the equilibrium length for a bond involving two different atoms.The spring constants of bonds involving only the film or the substrate atoms are chosen so as to reproduce the bulk elastic moduli for these systems.ThereisFig.4.The interaction energy between the ridge and trench features of the inset for D =4plotted as a function of d .The dots show the results of our model and the solid curve is f (d )of Eq.(5)multiplied by 2.2.1452P.Nath,M.Ranganathan /Surface Science 606(2012)1450–1457again some arbitrariness in the choice of the spring constant for a bond involving two different atoms.In particular,the choices of these param-eters influence the wetting layer and the2D/3D growth transition.Fur-ther,real Si or Ge is neither cubic,nor isotropic.In fact the diffusion of Ge atoms on Si(100)surface has been shown to be strongly anisotropic because of the2×1reconstruction.Due to the above reasons,we feel that our model can be used to study the qualitative behavior,but care should be exercised while making quantitative comparisons to other models and experiments.Finally,we note that for an individual atom, the elastic energy is typically a very small part of the entire barrier energy in the KMC simulations discussed in the next section.Despite this,it is be-lieved to play an important role in the structures formed during growth.5.KMC simulations with elasticityThe KMC method we have described in Section3is implemented using the hopping rate given by Eq.(1).The barrier energy is treated as a sum of the energy due to the bonds and the elastic energy due to the displacements.The displacements are calculated for afixed config-uration of atoms.Implicit in this method is the assumption that the elas-tic relaxation of the system takes place on much timescale smaller than that of surface diffusion.The elastic energy causes a reduction of the barrier for surface diffusion,since the springs in our model cannot have energies less than zero.Thus we can write E b=E bond−E el,where E el is the elastic energy contribution to the hopping barrier.To calculate E el,we need to calculate the difference of the total elastic energy of the system with and without the atom.Each elastic energy calculation needs a minimization of the energy function to calculate the atomic dis-placements.The elastic energy computation is the most expensive part of the entire calculation.We therefore consider some approximate schemes that will minimize the number of these calculations.Firstly we assume that the dominant contribution to the elastic en-ergy comes from the elastic springs of the atom to its NN and NNN atoms.Thus we approximate E el as E el¼∑bonds12k bond l bond−l eq bondÀÁ2 where k bond is the spring constant of the bond,l bond is the distance be-tween the atoms forming the bond,and l bondeq is the equilibrium bond length.This approximation allows us to calculate the elastic energy of the atom from its configuration without having to consider the elastic energy of the system without the atom.This reduces the total number of elastic energy calculations by a factor of2.In order to construct more approximations,we need to use some information about the elastic energy.Let us consider an adatom,i.e.a surface atom that is only bonded to one NN and4NNN atoms.Since the bond energy is only due to NN atoms,the total bond energy of an adatom is equal to E v.For an arbitrary surface,the elastic energy of an adatom depends on the detailed configuration of the surface.If the surface is notflat, the energy of adatoms located at different points on the surface will be different.As we shall show in Section1,this difference in elastic energy of adatoms on a given layer is negligible compared to the elastic energy of other surface atoms.We neglect this difference in our KMC simulations.This greatly simplifies the elastic energy calculation for adatoms.The elastic energy plays an important role in the diffusion of other surface atoms and has to be calculated accurately for those atoms.Since a significant fraction of atomic moves on the surface in-volve adatoms,our approximation leads to a significant time-saving. Further,since there are very few atoms that have fewer bonds than adatoms,we will set the time unit in our KMC simulations as the time for one adatom hop.With this approximation,we can write the hopping rate r as r=r ad exp−(E b−E b ad/k B T)where r ad is the hopping rate of an adatom and E b ad is the total barrier energy for an adatom.We choose our timestep so that r ad=1,thereby eliminating rejection of moves in-volving adatoms.This makes it possible to extend our simulations to longer times and study the growth of multiple layers.A similar approx-imation has been used recently in Ref.[50].We mention one more detail about the adatom energies here.The energy of afilm adatom on aflat substrate depends on the bond energy and the elastic interaction parameters involving thefilm and the sub-strate atoms.As we had mentioned earlier,there is someflexibility as regards the various parameters for a bond involving two different atoms.If we assume that the elastic parameters and the bonds are iden-tical for the interface layer and thefilm,then wefind that adatoms from different layers hop at the same rate.Any other choice will lead to an adatom hopping frequency that is dependent on the number of layers. Such dependence plays a role in the transition from2D layer-by-layer growth to3D mound growth.For instance,slower diffusion of adatoms in the upper layers will favor3D growth.A layer dependent adatom energy(elastic or bond)corresponds to a wetting layer in our system.A more careful examination of the wetting layer and the2D/3D transition will be the subject of a subsequent article.6.Results and discussionWe divide our results into two sections.In thefirst,we discuss elasticity calculations for various configurations.These calculations validate the model,and help us in constructing the approximations that were discussed in Section5.In the second,we discuss KMC sim-ulations of submonolayer growth.All calculations were carried out on a128×128system.The various parameters were chosen to mimic the Ge on Si(100).The lattice spacing in Si is2.715Åand that in Ge is 2.832Å.We performed Fourier transforms using the Fast Fourier Transform routine[56].The eigenvalues and eigenvectors for the substrate contribution were calculated using the singular value de-composition routine in LAPACK[57].The multidimensional minimi-zation of the elastic energy to calculate the elastic displacements was accomplished using the conjugate gradient optimizer in the GNU Scientific Library[58].The spring constants are taken as k L=73eV/nm2and k D=96eV/nm2to approximately reproduce the elastic moduli of Si.We have used the same spring constants for all pairs of atoms.The equilibrium bond length of a bond involving an Si and a Ge atom was taken to be equal to that of a Ge bond.To compare the submonolayer growth calculations with earlier work,we have used a temperature of500K,E v=0.8eV and E h=0.3eV.We have used the same value of E v for Ge atoms on bare Si(100)and Ge atoms atop Ge islands.The surface diffusion coefficient D s can be estimated as D s=(1/4)ν0exp(−E v/k B T)a s2which works out to about2×104s−1in length units of a s.Since we ignored all the anisotropic effects in the hop-ping barriers,we do not have any anisotropy in the surface diffusion constant.Besides the Si–Ge model system(4%misfit),we have per-formed calculations for systems with no misfit and8%misfit.Though we do not show this result here,we have verified that3D island forma-tion was observed at low coverages(b0.5ML)when we used8%misfit and spring constants larger by a factor of10.Each KMC simulation tra-jectory was generated on a64-bit Intel X5500series processor;and120 parallel trajectories were simulated to get good statistics.A simulation of1ML on these systems took about144h;but the simulation time in-creases for subsequent layers because of the increase in the total num-ber offilm atoms.6.1.Elastic energy calculationsIt has been proposed that the elastic interaction energy between two adatoms a distance r apart on aflat epitaxial layer is proportional to1/r3. This interaction is due to force dipoles located near each of the adatoms. On the other hand,steps on strainedfilms cause force monopoles and the interaction between steps distance r apart is proportional to log(r).For a periodic array of steps,the dependence is slightly modified.We test this dependence by performing elastic calculations for the configuration shown with d and D in Fig.4.We calculated the elastic energies of a flat surface E f,a surface with a ridge E r,a surface with a trench E t and a surface with both a ridge and a trench E rt.The elastic energy1453P.Nath,M.Ranganathan/Surface Science606(2012)1450–1457of interaction between the ridge and the trench features U int is given by U int =E rt −E r −E t +E f .According to continuum elasticity theory,this interaction energy is proportional tof d ðÞ¼−logsin 2πd þD ðÞ=L ðÞ"#ð5Þwhere L is the periodicity of the system.The constant of proportionality is related to the magnitude of the force monopoles at the step edges.We show that our interaction energy can be fitted to Eq.(5)with a single fitting parameter.The difference at smaller distances is simply a consequence of the difference between the continuum model and our lattice model at short length scales.In general,we expect to see some differences for moder-ate length scales,so it is remarkable that the graph in Fig.4agrees so well at moderate length scales.It is possibly due to the symmetric na-ture of the con figurations considered.Analogous to the usual thermodynamic chemical potential,we de fine an elastic chemical potential μel of an adatom on a given surface as the dif-ference in the elastic energy of the system with and without the adatom.For a given con figuration of a surface,μel changes for different locations asshown in Fig.5.The figure illustrates that the elastic contribution to the hopping barrier is largest at the island edges and is almost constant along a flat surface.Thus the atoms at the edges of islands have a reduced hopping barrier due to elasticity.The adatoms on top of the island have a smaller elastic energy than those on the bare substrate.The reduction in elastic energy causes slower diffusion of adatoms on top of islands.These two effects eventually lead to a transition from layer-by-layer growth to 3D mound formation in the KMC simulations.As a final demonstration of our elastic model,we show the dependence of μel at the island edges on the size of the island in Fig.6.The increase in island energy with island size follows a trend consistent with the continuum elasticity based calcu-lations of Refs.[41,37].As the size of the islands becomes large,the elastic energy increases much more slowly.We note that the chemical potential as calculated by us is signi ficantly smaller than that calculated in Ref.[41].This difference will be re flected again in the results of the KMC simula-tions in the next section.6.2.Submonolayer growthWhen Ge atoms are deposited on an atomically flat Si substrate,the deposited atoms diffuse on the surface and,occasionally,get attached to Ge islands on the surface.The resulting surface shows many islands typ-ically of different sizes as shown in Fig.7for different values of the mis-match.In Fig.7(d),we show a snapshot of a system with 4%mismatch,but a spring constant 5times larger.In this case,we can see some mul-tilayer islands too.If the deposition rate is slow,then the adatoms have suf ficient time to get attached to existing islands and the average size of islands is large.The island size distribution therefore depends on the ratio of the surface diffusion coef ficient to the flux,i.e.the D s /F ratio.We performed KMC simulations of submonolayer growth of Ge on Si.By choosing F =0.1ML/s and F =10ML/s,we get D s /F ratios of 2×105and 2×103(in units of a s 2)respectively.Following the previous section,we take a constant elastic energy for adatoms on any given layer.KMC simulations of submonolayer growth have been performed earlier using continuum elasticity model [41]and using an effective adatom interaction model [40].In both these cases the island size distribution for a given D s /F ratio was consistent with the scaling form [40,38]N s θðÞb s >2θ¼f s =b s >ðÞð6Þwhere N s represent number of islands of size s (s ≥2)when the surface coverage is θand the average island size is b s >.This form was derived for irreversible submonolayer growth conditions in the absence of elas-tic interactions.In Fig.8,we show the scaled island size distribution as a function of s /b s >for θ=0.1−0.29for F =0.1ML/s.The points seem to fit on a master curve indicating that there is very little mound formation.Chemical potential (meV)b)a)102030Elastic Displacementsc)Fig.5.(a)Graph showing μel of an adatom located at various distances on the surface along the line AB.The inset figure shows part of the surface con figuration with the line AB.(b)A contour map of μel on a surface with a 9×9island,and (c)cross-section of the surface showing the elastic displacement of atoms (magni fied by factor 10for better visibility)in the vicinity of an island.102030400200400600µe l (m e V )sFig.6.μel for a Ge adatom next to a square Ge island of size s .1454P.Nath,M.Ranganathan /Surface Science 606(2012)1450–1457。
技术经济学--蒙特卡罗风险分析方法
属于第一类（极低风险）的条件是所有的风险因素均好， 有测试证实，邻近有已开发区。 属于第二类（低风险区）的条件是所有风险因素均较好 或略差，有测试资料证实，但离已开发区较远。 属于第三类（中等风险）的条件是2-3个风险因素（即多 数）较好，1-2个风险因素中等，属于构造带内新区，远 离已开发区域。 属于第四类（高风险）的条件是1-2个风险因素（即少数） 较好，2-3个风险因素（即多数）中等，属于构造带内新 区，邻近无已开发区域。 属于第五类（极高风险）的条件是2-3个风险因素中等 （即多数），1-2个风险因素（即少数）较差或略差，属 于构造带内新区，无邻近已开发区域。
NPV 150 CF ( P / A, IRR, N ) 0
Developed by © Renjin Sun ,School of Business Administration, China University of Petroleum, Beijing
其二， 为各不确定因素选择合适的概率分布函数。
其三， 在0.000-0.999范围内，对每一个不确定因素
随机抽取一个数值，作为累积概率的随机值。 根据抽取的累积概率随机值，并结合累积概率 分布函数，求得相应的不确定性因素的一个值。
Developed by © Renjin Sun ,School of Business Administration, China University of Petroleum, Beijing
二、风险因素及其分布
（一）风险因素 1、勘探投资 2、地面工程建设投资 3、钻井投资 4、操作费用 5、原油价格 6、原油产量风险
（二）风险概率分布类型
蒙特卡罗方法 (Monte Carlo simulation)教材
2019/4/15
Monte Carlo模拟
9
2.Monte Carlo方法简史
Stanislaw Ulam (1909-1984)
S. Ulam is credited as the inventor of Monte Carlo method in 1940s, which solves mathematical problems using statistical sampling.
1
Monte Carlo模拟
第一章 引言 （Introduction)
1. 2. 3. 4. Monte Monte Monte Monte Carlo方法 Carlo方法简史 Carlo模拟的应用 Carlo算法的主要组成部分
2019/4/15
Monte Carlo模拟
2
1.Monte Carlo方法
2019/4/15
Monte Carlo模拟
5
2.Monte Carlo方法简史 Buffon投针实验
1768年，法国数学家Comte de Buffon利用投针实验估 计的值
L
d
p
2019/4/15
d
2L
Monte Carlo模拟
6
2.Monte Carlo方法简史
Problem of Buffon’s needle: If a needle of length l is dropped at random on the middle of a horizontal surface ruled with parallel lines a distance d > l apart, what is the probability that the needle will cross one of the lines?
用Monte+Carlo法模拟纯铜的静态再结晶
５期佟铭明等：用ＭｏｎｔｅＣａｒｌｏ法模拟纯铜的静志再结晶４８７ｌ＿５初始条件ｒ在模拟中采用１００×ｉ００个正六边形单元组成的正方形阵列，不同的灰度代表不同的晶体取向．向系统中一次加入１００个晶胚．温度为４００℃，应变率为ｏ．０５Ｓ～，纯铜（９９．９９５％）单位长度的晶界能（二维）【５１１≈６３ｘ１０～Ｊ／ｍ，初始状态晶粒大小约为１０ｐｒｏ所有计算结果均为计算五次后数值的甲均值．２模拟结果与分析未变形时体系的原始组织是经过１７０ＭＣＳ的正常品粒长大模拟后得到的（图ｌａ），图１中不同的衬度代表不同的晶粒取向（晶界上的锯齿是由于计算机屏幕的分辨率较低造成的，当分辨率足够高时晶界可以变得很光滑）．囝１模拟体系的微观组织ＭｉｃｒｏｓｔｒｕｃｔｕｒｅｏｆｔｈｅｍｏｄｅｌｅｄｓｙｓｔｅｍｗｉｔｈｄｉｆｆｅｒｅｎｔＭＣＳａｎｄｓｔｏｒｅｄｅｎｅｒｇｙ（ａ）＇ＯＭＣＳＦｉｇ．１（ｂ）２ＭＣＳ，Ｈ＝ｌ０ｘ１０—６１．（ｃ）２ＭＣＳ，Ｈ＝Ｉ９５×１０—６１，（ｄ）５ＭＣＳ，Ｈ＝２０ｘ１０６７，（ｅ）５ＭＣＳ，Ｈ＝２．０５×１０～－ｙ，（ｆ）５ＭＣＳ，Ｈ＝３０ｘ１０“？再结晶是一个存储能量的释放过程，所以在一定温度下静态再结晶行为完全由体系的初始变形状态（即初始存储能）决定为了研究初始存储能对静态冉结晶的影响，分别向每个单元中加入由塑性变形导致的初始存储能Ｈ。
为１．ｏ×１０一６１，１．２×１０—６１．１４×１０一６１，１．６×１０一６１，ｌ８×ｌｏ一６’，１．９×１０—６－ｙ，１，９５ｘ１０６－ｙ，２Ｏｘｌ０６１，２．０５×１０—６１，２１Ｘ１０—６１，２２ｘ１０—６～，２４×１０—６１，２．６×１０—６１，２．８×１０“１，３０ｘ１０“１，计算“上１５种条件下的静态再结晶过程．用Monte Carlo法模拟纯铜的静态再结晶作者：佟铭明， 莫春立， 李殿中， 李依依作者单位：中国科学院金属研究所刊名：材料研究学报英文刊名：CHINESE JOURNAL OF MATERIALS RESEARCH年，卷(期)：2002,16(5)被引用次数：8次参考文献(8条)1.A Manonukul;F.P.E.Dunne"Initiation of Dynamic Recrystallization under Inhomogeneous Stress State in Pure Copper " 1999(47)2.B Radhakrishnan;G.B.Sarma;T.Zacharia"Modeling the Kinetics and Microstructural Evolution during static Recrystallization-Mongte Carlo Simulation of Recrystallization " 1998(12)3.G S Srest;D.J.Srolovitz;M.P.Anderson Computer Simulation of Grain Growth- IV. Anisotropic Grain Boundary Energies 1985(03)4.P Peczak;M.J.Luton A Monte Carlo Study of the Influence of Dynamic Recovery on Dynamic Recrystallization 1993(01)5.D J Srolovitz;M.P. Anderson;G.S.Grest Computer Simulation of Grain Growth- III. Influence of A Particle Dispersion 1984(09)6.D J Srolovitz;G.S.Grest;M.P. Anderson Computer Simulation of Grain Growth- V. Abnormal Grain Growth 1985(12)7.M P Anderson;D.J.Srolovitz;G.S.Grest;P.S.Sahni Computer Simulation of Grain Growth- I . Kinetics 1984(05)8.Yoshiyuki SAITO;Masato Enomoto Monte Carlo Simulation of Grain Growth 1992(03)引证文献(8条)1.王狂飞.莫亚林.夏立军.米国发.傅恒志15SiMn钢奥氏体等温相变的元胞自动机模拟[期刊论文]-特种铸造及有色合金 2009(7)2.王狂飞.张锦志.米国发.孙瑞霞晶粒生长的Monte Carlo法模拟研究进展[期刊论文]-热加工工艺 2009(10)3.王狂飞.莫亚林.夏立军.米国发.傅恒志25SiMn钢奥氏体等温相变过程的数值模拟[期刊论文]-铸造 2008(8)4.王狂飞.孙瑞霞.王继红碳含量对奥氏体等温相变影响的数值模拟[期刊论文]-热加工工艺 2008(16)5.方斌.黄传真.许崇海.刘增文.王景海材料微观组织演变的蒙特卡洛仿真原理及应用[期刊论文]-材料导报2007(2)6.郑成武.兰勇军.肖纳敏.李殿中.李依依热变形低碳钢中奥氏体静态再结晶介观尺度模拟[期刊论文]-金属学报2006(5)7.元胞自动机方法模拟再结晶过程的建模[期刊论文]-机械工程材料 2005(10)8.兰勇军低碳钢奥氏体—铁素体相变介观模拟计算[学位论文]博士 2005本文链接：/Periodical_clyjxb200205007.aspx。
蒙特卡罗方法 (Monte Carlo simulation)
[0, )
• The random vector is uniformly distributed on the region [0,d)×[0,). Accordingly, it has probability density function 1/d. • The probability that the needle will cross one of the lines is given by the integral
2014-8-27
Monte Carlo模拟
18
4.Monte Carlo算法的主要组成部分
Monte Carlo算法的主要组成部分 概率密度函数(pdf) 必须给出描述一个物理系统的一组概率密度函数;
随机数产生器 能够产生在区间[0,1]上均匀分布的随机数 抽样规则 如何从在区间[0,1]上均匀分布的随机数出发,随机抽 取服从给定的pdf的随机变量;
p
2014-8-27
0
l sin
0
1 d
2l dAd d
Monte Carlo模拟 8
2.Monte Carlo方法简史 Enrico Fermi
• 1930年,利用Monte Carlo方法研究中子的扩散 • 并设计了一个Monte Carlo机械装置，Fermiac,用于计算核 反应堆的临界状态
2014-8-27
Monte Carlo模拟
12
Monte Carlo模拟
第一章 引言 （Introduction)
1. 2. 3. 4. Monte Monte Monte Monte Carlo方法 Carlo方法简史 Carlo模拟的应用 Carlo算法的主要组成部分
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
金字塔形向圆顶状转变
ห้องสมุดไป่ตู้
arXiv:0912.0053v2 [cond-mat.mtrl-sci] 7 Sep 2010
ergy of the system has to be computed repeatedly in or跳跃 负责 der to simulate the atomic hopping events responsible for 需要 形态演变 the morphological evolution. Applying more advanced The self assembly of quantum dots in heteroepitaxy 算法 algorithms for the solution of the elastic problem and exhibits very interesting physics and has possible appli取样 器件制造 the sampling of atomic hopping events, simulations uscation to device fabrication [1–3]. Growth of Ge or GeSi 适度稳健 纳米材料 标准 ing large and moderate system sizes in 2D [17–19] and 3D alloy nanostructures on Si(100) is the prototype example 分别 沉积 [20–24] respectively became possible. The model was exmost widely studied. Under typical deposition condi晶面 限制边界 tended recently to model (105) facets [25] and atomic intions, pyramidal islands bounded by shallow (105) facets 基片 混杂 浸润层 termixing with substrate atoms [26]. Alternatively, KMC form spontaneously on a wetting layer. They can 支配控制 then simulations can also be performed efficiently using more 经历 undergo transitions into multi-faceted dome islands dom近似 approximate forms of elastic interactions [27–29]. inated by much steeper (113) facets and bounded also by This paper is organized as follows. Section II explains other facets [4]. This shape transition gives rise to a bi解释说明 双峰 our multistate model for elastic solids which can account modal island size distribution路径 with enhanced dome size 均一性 簇 for both a shallow and a steep facet. The KMC simuniformity [5]. An atomic pathway based on step bunch认识 ulation results are presented in Sec. III. In Sec. IV, ising on the pyramids has been identified [6, 7]. land energies and the island transition rate are calculated In this work, we report atomistic dynamic simulations theoretically. In Sec. VI, a scaling relation between the of the pyramid-to-dome transition using a fast kinetic transition island size and the Ge concentration is derived. Monte Carlo (KMC) approach based on a multistate 大量 We conclude in Sec. VI with some further discussions. solid-on-solid model in two-dimensions (2D). Extensive simulations under a wide range of conditions are performed and a simple analytical description is presented. 可归纳的 II. MULTISTATE SOLID ON SOLID MODEL A scaling relation for the transition island size generaliz研究 able to three dimensions (3D) is investigated. Kinetic simulation of a strained film is much more Our model is based on a ball and spring square lat计算 challenging computationally than the unstrained case tice model of elastic solids for Gex Si1−x on Si [17]. The 长程 弹性相互作用 常数 because of the long-range nature of elastic interacsubstrate lattice constant is as = 2.72 ˚ A while the film 失配 节点 tions. First principles calculations [7–11] and molecularmaterial admits a lattice misfit ǫ = 0.04x. Each node dynamics simulations [12, 13] have provided important on the lattice represents an atom and it is connected to 弹簧 information on the energetics of the relevant surfaces and its nearest and next nearest neighbors by elastic springs 静态性能 力系数 steps, but are in general limited to the studies of static with force constants kN = 13.85 eV/a2 s and kN N = kN /2 模量 连续 properties at small system sizes. Continuum simulations respectively. This choice gives the correct modulus c11 正切 剪切模量 有帮助 in contrast are instrumental for investigating large scale of silicon and a shear modulus constant along tangential 对角线 and long time behaviors [14, 15], but lacking atomic disand diagonal directions. 不连续性 成核 creteness, nucleation events associated with island forIn this work, (100), (105) and (113) surfaces must be 水平 surface in the model natmation and shape transition cannot be naturally simeffectively simulated. A leveled ulated. KMC simulations based on lattice models are urally accounts for a (100) surface. However, lattice modhence unique in allowing large scale atomistic studies on els generally lead to islands with a single type of side斜坡 应变异质外延生长 the dynamics of strained heteroepitaxy. walls at 45◦ inclination or steeper [16–24]. A multi-state 弹性固体 弹性 lattice model for elastic solids, Using a ball and spring extension has been introduced recently in Ref. [25] to Orr et. al [16] conducted early KMC simulations of effectively model the much shallower (105) facets of a strained layers in 2D. In the simulations, the elastic enpyramid in 2D. We now further generalize it to simulate