1、陈效双 中科院上海技术物理研究所、上 海科技 大学 计 算 材料学 计算材料学 第 四 章 原子模拟方 法本章提纲 1. 原子模拟基础 2. 结构优化方法 3. 分子动力学方法 4. Monte-Carlo方法 2计算材料学第四章 4-4. Monte-Carlo 方法 3Monte-Carlo方法 4Dice of Casino: random number Monte Carlo hotel, Las Vegas Monte Carlo方法的来由 随机 数 5 Monte Carlo 的思想起源于von Neumann等人对裂 变材料的中 子扩散问题 研究。 在Metropolis 等人
2、建 立了计算机模拟的Monte Carlo 方法以后, 这一方 法在解决多粒子体系的相关物理问题的研究中被 广泛使用 。 最早利用计算机模拟研究统计力学体系以及相关 物理问题的是Metropolis 等人于1953 年在美国Los Alamos 国立 实验室的第一代电子计算机上完成的, 并由此建立了计算机模拟的Monte Carlo方法。 Monte Carlo方法的发展 历史 6Monte Carlo方法的发展 历史 1968 年, Wood建 立了NPT 正则系综 的Monte Carlo 方法 。 1969 年Norman 和Filinov 建立了巨正则系综的Monte Carlo抽样方
3、法 。 1987 年 Panagiotopoulos 把 Monte Carlo 方法应用于 Gibbs 系综 。 1986 年, Voter 在点 阵气体模 型基础上 提出了描 述表面 原子运动Kinetic Monte Carlo方法 ， 被迅速应用于 薄 膜生长模拟 。 7 MC方法在数学上成为 随机模拟方法 、 随机 抽 样技术等； MC方法的基本思 想是： 针对某 一具体 的问题 ， 通过建立一个概率模型或随机过 程模型 ，使 其参数等于实际问题的解。 Monte Carlo方法的基本 思想 8What Is MC and What is it for? MC 随机地探索一个系统的具
4、有几率状态 ， 与物 理 上期 待的 量相匹配 Stochastic 意 味着 涉 及或含 有一 个随机 变量 或多个 随机 变 量 , 在实际中意味着那个方法基于随机数的值做事情 MC 用于得到热力学平均值，热力学势能 （ 来自 平 均值 ） ， 和研究相变 MC 有 许多材料科学以外的 其 他应用, 覆盖了大量的用随机 数的方法 Called Monte Carlo 正是因为那个赌博发生 很多机会 ! 9Monte Carlo积分圆周 率 的计算 4 circle hit square shot St St = = 102 1 2 1 d () () d () () x x x x F x
5、f x fx Fx x x = = Monte Carlo积分 max 21 21 1 max 1 () () x xx xx Ff = = x y x 1 x 2 f(x) 11分子动力学方法的局限 12统计力学与热力学简单回顾 13 微观系统的几率分布和它的 Ha mil ton ian 是与宏观 系统的自由能相联系统计力学与热力学简单回顾 14 一系列的微观状态的平均可以获 得热力 学量系综平均 15 平均是很多微观状态的平均 ， 而不是动力学系统 时间上的平均系综平均 16 简 单 抽样方 法在随 机选取 事件（ 如：空 间位型 ）时， 只 考虑事件 存在的可 能性，而 不考虑事 件本
6、身存 在的概率 。 在 简 单抽样 过程中 ，统计 平均量 是通过 已出现 事件及 其 存 在 的概率 进行统 计平均 而给出 。 抽样满足 波尔兹曼 统计，物 理量A 的期 望值由下 式给出： 简单抽样 = = M l M l kT H kT H A A 1 1 ) / ) ( exp( ) / ) ( exp( ) ( ) ( 17 重要抽样方法在随机选取事件 （ 如：空间位型 ） 时 ， 同 时考虑事件存在的可能性 和 事件 本身存 在的概 率 。 事件本身存在的概率通过Markov 链中的转移矩阵而确定 。 在重要抽样过程中 ， 统计平均量是通过已出现事件直接 进行统计平均而给出 。
7、对于有限温度下的统计 ， 不再是按照简单抽样那样随机 生成样本 ， 而是 对温度T起 重要 作用 的那 些样本进 行择 优 抽样 。 即采取与分布函数相同的 概率形 式 抽取 样本 。 重要抽样 18 An algorithm to pick a series of states so that they asymptotically appear with probability (s)=exp(-E(s) 1. Assume we are in state s i 2. Choose a new state s*, and define E=E(s*)-E(s) 3. If E0 then
8、 accept s* with probability exp(-E) 5. If we accept s* then set s i+1 =s* and return to 1. in state s i+1 6. If we reject s* then return to 1. in state s i This is a Markov process, which means that the next state depends only on the previous one and none before than (not like MD). The Metropolis Al
9、gorithm (General) 19 We only need to consider atomic positions (since momenta do not enter into (s) 1. Assume the atoms have positions (r 1 , r j , r N ) 2. Choose a new set of positions r j * = r j + (2-1)r, where j is chosen at random and r is a fixed step length and is a random number between 0 a
10、nd 1. 3. Find E=E(r 1 , r j * , r N )-E(r 1 , r j , r N ) (note that this can be done quickly be only recalculating the energy contribution of atom j and its neighbors) 4. If E0 then accept r j * with probability exp(-E) 6. If we reject r j * then change nothing and return to 1 The probability of se
11、eing any set of coordinates r will tend asymptotically to Metropolis Algorithm for Interacting Atoms 20 ( ) ( ) ( ) exp VZ = r rr1. Start with some configuration 2. Choose perturbation of the system, e.g., move one atom from R to R 3. Compute energy change due to that perturbation 4. If E0 ，accept p
12、erturbation with probability exp-E/kT 5. Choose next perturbation 6. 7. Run sufficient MC step and average over all these states R R Metropolis Monte Carlo 样本的分布为正则分布， 可以很好模拟NVT系综 Metropolis et al., J. Chem. Phys. 21, 1087 (1958) 21随机数 22Monte Carlo 轨迹 23Thermal Average in MCMC Step Energy Equilibra
13、tion period: Not equilibrated, thermal averages will be wrong Converged: Equilibrated, thermal averages will be right24Obtaining Thermal Averages From MC The MC algorithm will converge to sample states with probability (s) = probability expected from statistical mechanics. So a thermal average is gi
14、ven by Note that N mcs should be taken after the system equilibrates Fluctuations are also just thermal averages and calculated the same way ( ) ( ) 1 1 mcs N i si mcs A As s A N = = = ( ) ( ) ( ) 2 1 1 mcs N i si mcs A A A As s A N = = = 25Monte Carlo方法模拟系 综性质 的精 度 26Finite Size Effect 27小结：MC 方法是对
15、系综的有 效抽样 28MC on Spin Systems MC is often used for spin systems, used to model, e.g., magnetism and alloys. The Ising Model each site on a lattice has a spin that can be up or down. Simple model for magnetic moments interacting in a crystal. ( ) 1 for up spin, -1 for down spin J0 ferromagnetic order
16、ing tendency ij i j ij i EJ = =+ 29Metropolis Algorithm for Ising Model We only need to consider spin states 1. Assume the spins have value ( 1 , j , N ) 2. Choose a new set of spins by flipping, j * = - j , where j is chosen at random 3. Find E=E( 1 , - j , N )-E( 1 , j , N ) (note that this can be
17、 done quickly be only recalculating the energy contribution of spin j and its neighbors) 4. If E0 then accept spin flip with probability exp(-E) 6. If we reject spin flip then change nothing and return to 1 The probability of seeing any set of spins will tend asymptotically to ( ) ( ) ( ) exp EZ = 3
18、0Ising Model MC: Magnetization vs. T T T c M M=1/site 0 0 31Monte Carlo 方法的利弊 Conceptually simple Easy to implement Can Equilibrate any degree of freedom /No Dynamics needed Accurate Statistical Mechanics Advantages 32 No Kinetic Information. Cannot do explicit time evolution (except under some circ
19、umstances with kinetic MC), so no transport coefficients May be slow to explore large phase space since small steps can take to long to move the system through bottlenecks (e.g., large molecules). This can be improved with force-bias MC. Requires many Energy Evaluations Stochastic nature gives noise
20、 in data Not easy to get entropy/free energy Monte Carlo 方法的利弊 Disadvantages 33 Explicit time evolution: Easy in MD, e.g., transport coefficients from time correlation functions Hard in MC, must use kinetic MC formalism NPT ensemble: Hard in MD, since must use “tricks” Easy in MC, since Metropolis a
21、lgorithm applies naturally Discrete models with not explicit equations of motion (e.g., lattice of spins to model magnetism) Hard in MD since you cannot naturally define update steps Easy in MC since the states are straightforward to step through (e.g., spin flips) MC vs. MD 34 Exploring complex pha
22、se space with many local minima Hard in MD since barriers are difficult to overcome and you can get stuck in one minimum but at least forces guide you Easy in MC since you can jump between minima through unphysical steps Hard in MC when you need to explore space without clear local minima, since ste
23、ps may need to be small and are may not be clearly directed (e.g., large molecules) MC + MD together: This can be done to get the best of both methods. MC vs. MD 35Kinetic MC (KMC) Kinetic Monte Carlo (KMC, sometimes called Dynamic MC) is meant to address materials-related problem whose physical out
24、comes are very much governed by the energy barriers between various possible states, or local configurations of the system. The rates for these physics-based transitions control the dynamics. 36Limitation of Metropolis walk Let W( i f ) be the transition probability per unit time that the system wil
25、l undergo a transition from state i to f and let be the energy change during the transition. Standard Metropolis walk: is the usual MC time step. Metropolis MC samples phase randomly jumping here regardless of barriers required to reach a particular lower- energy minimum, thus can not reproduce the
26、real dynamics of system. Example for barrier between two sites i and j with energy E between minima and additional barrier. Hop from i to j would cost more energy than j to i. 37Thermally-excited process in KMC the barrier depends on initial and final position, that is, it is not symmetric: 38Typica
27、l application of KMC Surface diffusion Molecular beam epitaxy (MBE) growth Chemical vapor deposition (CVD) growth Vacancy diffusion in alloys Dislocation motion Compositional patterning in alloys under irraditation Topological constraints vs. entanglements of polymers Coarsening of domain evolution
28、and more 39薄膜生长的KMC模拟 Adatom diffusion h = exp(E b /kT) Large scale Long time scale Temperature effects Bond energy Diffusion barrier 40薄膜生长的KMC模拟 100K 300K 400K 450K Morphology of Au epitaxy films at various temperature, The coverage is 0.2 (a) 100 K (b) 300 K (c) 400 K (d) 450 K 41薄膜生长的KMC模拟 Homoepitaxy on Fe(110) Kinetic MC simulation 42薄膜生长的KMC模拟 Zr 65 Al 7.5 Cu 27.5 film deposited at RT on SiO 2 Kinetic MC simulation 43
