蒙特卡罗方法(Ⅲ)

蒙特卡罗方法（三）各种温度下金薄膜外延生长的KMC模拟结构特性：材料生长过程、径向分布函数，原子数密度分布，体系晶体结构，键角分布，健对分布等团簇结构，薄膜生长，合金偏析，高分子材料等热学性质：Lindemann指数，热容量，熔点等Cu19团簇随着温度的降低，体系发生微相分离，左：Micelle右：lamellar研究对象研究性质磁学性质：磁滞回线，矫顽力，居里点等§4蒙特卡罗方法应用---经典粒子体系一个由N个原子构成的系统，通过势函数计算总能，总能E是系统构型R(r)的函数制造随机数，使原子产生随机尝试运动，将原状态R更新为新状态R’通过势能函数计算出新状态与原状态的系统总能量差δE在某一温度下，系统从R构型到R’构型的接受概率为：P=exp(-ΔE/kBT)产生η∈（0，1）之间的随机数，如果Pη，更新状态R→R’如果PQ，保持原状态R,,1()()2ijiiijijiErF通过各种方法产生随机数，经过大量的蒙特卡洛模拟步，系统将达到一个稳定的平衡状态原子随机交换原子随机运动两种尝试运动如果体系较大，可以通过局部能量变化计算整个体系的能量差，以减少计算量若为Bath算法则：P=exp(-ΔE/kBT)/[1+exp(-ΔE/kBT)]§4蒙特卡罗方法应用---经典粒子体系读入系统原子坐标随机交换原子随机移动原子各方向随机移动距离0-0.5Å模拟过程每100MC步输出平均原子能量和最终结构坐标获取二次信息：径向分布函数，结构特征参数计算，原子数分布等4rnrr0.02rnmr----离团簇中心的距离nr----r处薄厚度Δr为处壳层内的原子数Fe和Co原子随机交换位置初始结构文件可以采用XYZ,PDB，BDL等格式能量以eV为单位，最终结构可以为XYZ,PDB，BDL等格式，通过结构可视化程序VMD，ME显示初始结构和最终结构随机运动原子数分布通过随机数控制§4蒙特卡罗方法应用---经典粒子体系,,1()()2ijiiijijiErF2020exp[(1)]exp[(1)]()1()1()eeeerrABrrrrrrr,()ijijjjifr20exp[(1)]()1()eeerfrfrrr1()()()[()()]2()()baabaabbabfrfrrrrfrfr30()(1)0.85inineinFF30()(1)1.15iinoeieFF()[1ln()]()eoeeFFJohnson合金法对势部分嵌入能部分电子密度参考文献：X.W.Zhou，ActaMaterialia49(2005)§4蒙特卡罗方法应用---经典粒子体系势函数参数模拟系综:NTV阶段距离:8.5Å模拟温度:200K~500K模拟步数:5×106边界条件：自由边界条件234-3-2-1012345Fe-CoFe-FeCo-CoGEAMEnergy(eV)r(Angstrom)参考文献：X.W.Zhou，ActaMaterialia49(2005)势函数变化曲线MC模拟参数§4蒙特卡罗方法应用---经典粒子体系02000040000-3.2-3.0-2.8-2.6-2.4-2.2-2.0-1.8Energy(eV/Atom)MCstep平衡结构平均原子势能约为3.01eV最终结构123401020Fe7Co6Nr(0.1/nm)初始结构平均原子势能曲线原子数分布密度正二十面体结构具有幻数（如13，55，147…）个原子的团簇许多都具有正二十面体的低能结构如Cu，Co，Ni等§4蒙特卡罗方法应用---经典粒子体系1234560102030Fe28Co27Nr(0.1/nm)050000100000-3.6-3.3-3.0-2.7-2.4-2.1Energy(eV/Atom)MCstep初始结构最终结构平均原子势能曲线原子数分布密度平衡结构平均原子势能约为3.37eV§4蒙特卡罗方法应用---经典粒子体系逾渗问题§4蒙特卡罗方法应用---经典粒子体系逾渗转变，指的是在庞大无序系统中随着联结程度，或某种密度、占据数、浓度的增加（或减少）到一定程度，系统内突然出现（或消失）某种长程联结性，性质发生突变，我们称发生了逾渗转变，或者说发生了尖锐的相变。统计力学中关于相变的一个著名的标准模型是Ising模型。该模型最初是由W.Lenz提出用来作为铁磁性的一个模型，后来成为他的研究生Ising的博士论文题目。Ising在1925年给出了一维情况下的解，该解显示，在一维下Ising模型没有相变解。Onsager在1944年得到二维Ising模型的准确解，二维时就有了相变。而三维情况下至今还没有严格解，必须依靠数值计算。材料的磁性本质上是量子现象，N.Bohr等人早在量子力学诞生之前就证明了经典力学系统不可能有磁性。现代磁性理论中的关键要素是电子自旋和相应的磁矩。当大量自旋的磁矩取相同方向时，就产生了宏观磁性。物质在外磁场H中的磁化强度M为Ising模型MHχ称为磁化率，§4蒙特卡罗方法应用---经典自旋体系按磁化率的行为，磁性分为：不同磁性物质的磁化率随温度的变化关系§4蒙特卡罗方法应用---经典自旋体系顺磁性铁磁性反铁磁性Heisenberg首先用量子力学方法计算了自发磁化强度，他的计算结果是，量子力学效应来自于电子自旋对中的交换能σi⋅σj。L.Landau和Lifshiz发展铁磁性的磁畴假说，后来得到实验验证,成为铁磁性磁化曲线的理论基础F.Bloch后来又考虑了晶体结构周期性和各种方向的自旋组态，发展了自旋波理论。Langevin给出了顺磁性的经典理论。J.H.vanVleck其后发展了顺磁性的量子理论，由于他对磁学研究的贡献获得1977年的Nobel物理学奖。1928年vanVleck指出铁磁物质的自发磁化可以用电子间的特殊相互作用来解释，该相互作用使得各原子间的电子自旋平行取向。§4蒙特卡罗方法应用---经典自旋体系Ising格子模型：每个自旋位于格点上，自旋取向向上或者向下。自旋值只可能取两个值，故σ=±1。这是一个经典的自旋系统，而要真正用完整的量子力学方法处理时，需要将自旋角动量服从的量子对易规则考虑进来，这样的量子自旋问题仍然是研究的前沿领域。Ising模型Hamilton量为,1IsingijijBiijiHJHIsing模型Ising模型的重要之处在于它的简洁性，但却包含了相变的重要物理内容。其中求和的下标i,j表示近邻自旋对§4蒙特卡罗方法应用---经典自旋体系如果一对自旋方向相同，则能量为−J，相反为J。当J0，体系更趋向于把所有自旋取向排成一致方向以使得能量最低，在没有磁场情况下产生了自发磁化，这是铁磁性。§4蒙特卡罗方法应用---经典自旋体系当J0，自旋对的取向相反才可能使能量最低，宏观不表现磁性，但当加上磁场后逼迫自旋取向相同，产生磁化，这是反铁磁性。当温度升高时，热激发效应使得某些自旋取向随机反转，逐步使系统无序化，在足够高的温度下自发磁化消失成为顺磁性。enforcedferromagnetism，spin-glass-like铁磁性顺磁性反铁磁性亚铁磁性体系的每种自旋构型视为一种状态α，对应于特定的微观态能量Eα。当系统处于温度为T的环境中时，系统的自旋构型并非是一成不变的，它与环境之间不断交换能量，最终达到一个平衡态，在有限的温度下系统将围绕此平衡态进行涨落。在正则系综中系统处于状态α的概率正比于Boltzmann因子:统计力学分布求和是对所有的自旋构型。Boltzmann分布并不要求系统处于基态，只是说系统处于高能态的可能性较小。只有当T→0时，系统才可能处于基态。因此，磁场的作用是使系统有序，而温度使它趋于无序。正是这两种相反的竞争力导致自旋点阵上的有序—无序转变。,()()BBEkTEkTepZTeZT§4蒙特卡罗方法应用---经典自旋体系所有物理量均可通过配分函数解析地求解得到，如1limlnBNFkTZNlnBMkTZH2lnlnln(1)BBUZZkTZkTT==UCT0limHMH平均自由能磁化强度内能比热磁化率§4蒙特卡罗方法应用---经典自旋体系UEpE2221NBiiMp22EpE21NBBiiMNp222()EEE222()MMM222()(),BBEMCkTkT用Boltzmann统计分布将它们表示出来数值模拟求解其中σ为每个格点的平均自旋，它对时间的平均，也为对系综代表点α（系统的自旋构型）的平均。§4蒙特卡罗方法应用---经典自旋体系一维Ising模型的严格解热平衡状态下的配分函数121212231111121111exp[()]2exp[()]NNNNNZKKSSS§4蒙特卡罗方法应用---经典自旋体系11(1,2,1)iiiSiN112121111122(2cosh)NNKSKSKSNSSSZeeeK因子2是因为σiσi+1的取值有4种情形，但Si只有两种。故：新变量为H=0，系统Hamilton量为111NiiiHJBJKkT12i1iNJH≠0，采用环形周期性边界条件，即σN+1=σ1，系统的Hamilton量为111NNiiBiiiEJH1211111exp[()]NNiiiiZKI1limlnln(2cosh)BNFkTZkTKN它是温度的解析函数，除了T=0或∞外，F及其温度的高阶导数没有其它奇点，因此不存在相变。每个自旋的自由能§4蒙特卡罗方法应用---经典自旋体系TkHIBBBJKkTKIKIKIKIeePee配分函数可表示为12112231111111||||||||NNNNNNZPPPPTrP=由本征值方程2220coshsinhKIKIKKKKIKIeeeIeIeee定义一个2×2的矩阵11101(0,1)1(1,0)如果将自旋取值写成矩阵的形式，0111§4蒙特卡罗方法应用---经典自旋体系当H=0时，2coshK=λ+，2sinhK=λ−。低温下，K》1时，λ+≈λ−，自由能222ln{coshsinh}KKKFeIeIe它仍是温度的解析函数，也没有相变。设M∞为T=0时的磁化强度，可得24sinhsinhKMIMIe当去掉外磁场时，H=0，M/M∞=0，即除了T=0点外无自发磁化，没有有序化相变。当H≠0时，可得到自由能与自由边界条件下的结果一致。ln(2cosh)FkTK§4蒙特卡罗方法应用---经典自旋体系二维Ising模型的Onsager解，比热：2224coth2{()()sec2[(2tanh21)()]}2BCkKKhKK122201()1sind12220()1sindK(κ)和E(κ)分别为第一类和第二类椭圆积分K(κ)在κ=1处是发散的，因此相变点是sinh(2)21cosh(2)cccKK可得正方格子的临界温度值（单位：J/kB）平均场下的结果为Tc=2（z=4）。对于三角格子，严格解为Tc=3.641，平均场值为Tc=3