数值计算方法讲座3

MonteCarlo方法的应用磁性材料与相变磁性材料的重要特征是高温处于无序态〔无磁性〕，低温处于有序态〔有磁性〕。在物理上，同一体系的这样的不同的宏观状态，称为‘相’。自然界的法那么是，体系如果存在不同的相，随着控制参数的改变，体系会发生相变----宏观状态的突变，这一现象称为相变。材料的宏观磁性，来自微观的分子或原子的磁子的有序状态。支持有序态的力量是磁子之间的相互作用。破坏有序态的力量是来自环境或杂质的无规相互作用，常常用温度描述。这样两种力量的竞争是相和相变的根源。磁性材料通常是固体问题：为什么？因为液体一般无法保证磁子的有序状态。设体系已到达平衡态，统计物理学很好地描述磁性材料的性质。不过，我们首先需要建立体系的微观模型。问题：什么是平衡态？物理学家的重要思维方法之一是简化问题，寻找现象的本质。理想晶体、无杂质代表晶格上的磁子，最简单情形假设只有两种取向，用描述假设只有最近邻的磁子存在相互作用，但磁子可以与外磁场相互作用例如，二维Ising模型的Hamiltonian〔能量〕<ij>记最紧邻相互作用。Hamiltonian，在理论力学中，我们可以写下体系的根本运动方程，即Hamiltonian方程。但是，由于Ising模型的特殊性，Hamiltonian方程不存在。不过，统计物理学的出发点不是运动方程，而是假设每个微观状态出现的概率遵从一个统计分布。对于与温度为T的大热源接触的体系，在平衡态时，这一分布由体系的Hamiltonian决定，热源热源体系体系LudwigBoltzmann,1844-1906宏观物理量是对所有微观状态平均的结果。例如，(单位体积)磁化强度MTCT对空间维数d=1和2，存在准确解。但对d=2,准确解极复杂。设h＝0，磁化的曲线如图，对d=1,Tc＝0，没有真正的有序态。对d>1,在T=Tc处，发生二级相变。M在Tc处连续，但其一阶导数不连续必须强调从理论上说，真正意义的相变只对粒子数无穷多的体系存在如果是一级相变，M在Tc处不连续MTCT在远离相变点的区域，数值模拟和物理量的测量一般比拟容易。困难在于如何模拟系统在相变点附近的特征行为。二级相变的特征：物理量遵从幂次行为具有普适性！例如，磁化其中约化温度d＝2，这里称为‘临界指数’，只与几个重要物理条件有关：例如，对称性空间维数相互作用的力程临界指数把自然界的二级相变分为假设干普适类阅读材料：-----------------------------------------------------------------------------------------关联长度按幂次发散或者説关联函数取幂次形式关联函数〔或两点关联函数〕…，当足够大称为关联长度。在二级相变点附近，设当，发散。这时关联函数取幂次形式二级相变这两个特征高度非平庸，其关键是关联长度发散。因为关联长度在相变点附近发散，相应地统计涨落无穷大，所以微观细节不影响体系的大范围性质，如幂次行为、临界指数等，这便是普适性。-----------------------------------------------------------上述的行为特征可以一般表达为所谓的标度形式例如，对磁化强度，这里b是任意的标度因子。问题：如何从标度形式导出磁化的幂次行为？进一步，标度形式可以推广到有限体系，其中L是体系的尺寸，必须足够大。换句话说，当L足够大，但不是无穷大时，体系具有类相变行为。我们研究什么？体系是否存在相和相变？相变点在哪里？一般说来，如果存在不同的相，便会发生相变当然，对于理论家或计算理论家，必须先建立模型确定相变的阶如果是一级相变，测量不连续性等如果是二级相变，测量物理量的幂次行为和临界指数等当然，物理学家和材料学家所关心的问题侧重点会有所不同，甚至相当不同。例如，对高温超导材料，物理学家更关心超导的机制，即如何建模，得到高温超导相相变点在哪？与实验是否符合？相变的阶，相变的特征行为如果根本的物理机制一时不清楚，如何构造唯象的理论等材料学家可能更关心材料的性能，实用价值等这里我们必须指出，真实的磁性材料通常没有Ising模型描述的那么简单。所以，我们需要不断改良模型，以其更好地描述实际材料的特性。例如，磁子可以有多个取向，甚至连续取向，甚至需要量子化Potts模型XY模型－>量子自旋模型是量子力学中的算符相互作用次近邻作用，在次近邻作用，……长程作用等一般地，短程作用同属一个普适类。无序和掺杂根据无序和掺杂的不同，物理系统的相、相变和物理性能等会非常不同。例如，Tc这是所谓的相图不同的晶格，晶格的缺陷，建模的问题不是纯粹的数值模拟问题，解析的理论研究〔准确方法或近似方法〕也需要建模。但是，有些模型会更适合特定的方法。所以，我们应中选择最适合数值模拟的模型。归纳起来，数值模拟的根本思路如下：科学问题的表述选定或开展一种理论，模型的构造研究内容和研究方案确实定选定计算方法，克服可能的困难计算，分析结果，撰写论文数值模拟的关键是如何‘模拟’微观世界，和选定或提出有效的计算方法。第二节Ising模型的MonteCarlo模拟现在的问题包括两方面：如何计算物理量的平均值 ， ， 如何准确测量相变温度〔或者说〕，和相变点附近的物理量的特征行为，如幂次行为和临界指数等。〔不在相变点附近的物理量的测量一般较容易。〕对多体问题，直接数值计算物理量的平均值是不可能的。例如，对二维格点，当L＝100时，格点数为，那么需要求和的项数为，这是天文数字，根本无法完成。对这样的多自由度求和〔或积分〕问题，MonteCarlo模拟方法十分有效。我们应用重要抽样方法，给予概率分布的意义，引入恰当随机过程，产生一系列自旋构形 当足够大时， 遵从分布，那么有 例如 ，L是格点尺度。思考题：如果不能给予概率分布的意义怎么办？关于随机过程的平衡态我们扔掉前面个自旋构型，是为了等待动力学系统到达所谓的平衡态。对单自由度或少自由度系统，通常不大，但对多自由度系统，可以很大。还与初始状态有关，低温的初始状态较大。我们说系统到达平衡态，并非指自旋构型不再作动力学演化，处于静止状态。而是指自旋构型的概率分布不再随时间演化。换句话说，如果我们在足够长的时间里对自旋构型作物理量的平均值，物理量不随时间改变。构造随机过程，关键是选择一个恰当的转移概率，满足各态历经和细致平衡条件。仔细点说，这里的问题有两方面迭代方式即每次迭代更新多少个自旋。例如，我们可以每次更新一个、二个、三个，或N个自旋。一般而言，每次更新一个自旋比拟有效，称单自旋迭代法〔或翻转法〕。因为一次更新太多自旋，他们的能量改变会互相抵消。一旦决定采用单自旋翻转法，还存在扫描方式的选择。例如，我们可以选择随机扫描法，即随机找一个自旋进行迭代；也可以选择顺序扫描法，即按一定次序扫描，如一行一行地扫描。对一般问题而言，扫描方式不影响结果。但顺序扫描法比拟省时间。MonteCarlo算法例如，Heat-bath算法 选定，取 注意：这一算法的跃迁概率与的值无关！这与Metropolis的方法不同。 的能量 的能量由于每次只迭代一个自旋，与无关的自旋的能量不必计算。设 各态历经是显然的。 细致平衡 边界条件：自然边界周期边界练习： 构造Metropolis算法 构造二自旋迭代的Heat-bath和Metropolis算法在计算机上实现Heat-bath的算法 选定 计算 产生均匀分布的随机数 ，如果 否那么 0 1 概率 计算机程序DO20I=1,NIR(I)=I20CONTINUEIR(0)=NIR(N+1)=1DO40I=-4,4,1XX1=DEXP(-AT1*I)XX2=XX1+1.D+00/XX1PABHM1(I)=XX1/XX240CONTINUEDO120I=1,NDO120J=1,NCthisisHeat-bathalgorithmIEM1=S1(I,IR(J-1))+S1(I,IR(J+1))c+S1(IR(I-1),J)+S1(IR(I+1),J)callranecu(iseed1,iseed2,ranec1)IF(ranec1.LE.PABHM1(IEM1))THENS1(I,J)=-1ELSES1(I,J)=1ENDIF120 CONTINUE在远离相变点的区域计算物理量的平均值一般不太难，因为统计涨落较小。关键是如何较准确测量相变点和物理量的幂次行为以及临界指数。这里的问题是，严格的相变现象只在无穷多粒子体系存在。Naïve的方法比方，我们可以测量磁化强度M及其k次矩 M 当 这里的要点是M对温度的导数不连续。事实上，因为，M在低温一侧是发散的，。在数值模拟中，因为L总是有限，所以，无法观测到这样的不连续性。一般只能逐步增加L，看曲线的趋势。这样的方法难以准确定出相变点，也难以测量临界指数，已经较少应用。有限尺度标度行为方法相变点附近系统具有标度行为 对有限尺度体系，重整化群方法可以论证，具有有限尺度标度行为 利用这一标度行为，可以较简单测量引入Bindercumulant 由标度行为形式 当 从这一行为，可以测定求U的温度的导数 对二次矩 小结一下，数值测量的要点是构造 恰当的物理量，减少待测量的参数Ising模型的另一种应用假设代表粒子A，代表粒子B，而两种粒子的数目是守恒的。换句话说，Ising模型的磁化M是守恒的。对这一情形，Markov过程的自旋不能任意翻转，必须保证磁化M守恒。例如，单自旋翻转是不允许的。保证磁化守恒的最简单的实现，是交换两个自旋的值。习题：构造交换两个自旋的Metropolis和heat-bath算法。磁化不守恒的Markov过程称模型A，而磁化守恒的Markov过程称模型B。模型B的演化比模型A慢得多。第三节动力学慢化以Markov过程为根底的MonteCarlo方法----目前这是几乎唯一的方法，的显而易见的弱点是，产生的状态在时间方向上会有‘关联’。什么是状态的时间关联？在一定的演化时间内，状态的大局部自旋没有改变。这些大局部自旋相同的状态，几乎可以看成同一状态。假设在时间内状态没有改变，那么称为关联时间。显然，当存在时间关联时，MonteCarlo模拟的误差应当重新估计，修正为这里M是抽样总数。当不大时，这对数值模拟的影响并不大。但是，在一些重要而特殊的情形，例如，相变点附近，发散。这会极大地降低MonteCarlo方法的效率，称之为dynamicslowingdown。如果是在二级相变点附近，便是所谓的criticalslowingdown----临界慢化。为了定量描述状态的时间关联，可以引入时间关联函数 包括统计平均和对平均。当t足够大，一般地，称之为关联时间注意：当我们对平均，便假设系统已经处于平衡态，我们考察系统在平衡态的动力学涨落。思考题：为什么当足够大，？对二级相变， * ，这里z是所谓的动力学临界指数，用于描述关联时间的发散。 当 ，这是动力学标度行为的物理根底 * ， 当， 这便是通常说的临界慢化 为什么会出现临界慢化现象？设，如果每次迭代一个自旋， 在温度较低时，相同取向的自旋倾向于联成一片，也即所谓的‘cluster’。例如，当的所有或大多数近邻自旋，是难以取值－1的。当温度远低于相变温度时，系统处于有序态，大多数自旋形成cluster，并取确定的值，＋1或－1。换句话说，统计涨落较小，这不会导致什么数值模拟上的困难。当温度在相变温度附近时，系统处于临界状态，序参数为零，但大多数自旋形成cluster，不取确定的值，时而＋1时而－1。在数值模拟中，对于较大的cluster，难以从＋1翻转到－1。问题：如果每次迭代多个自旋呢？阅读材料：--------------------------------------------------------------------------------------理论上说，criticalslowingdown与动力学标度行为密切相关。设 ,当t足够大， 假设 当然 显然，指数上的b的因子必须自身抵消掉〔为什么?〕，即与b无关 验证：思考题：为什么 ？提示：当t足够大，----------------------------------------------------------------------对于实际的数值模拟，动力学慢化意味这什么？必须等待非常长的时间〔即〕，才能到达平衡态。即便已到达平衡态，必须等待非常长的时间才能得到一个新的〔微观〕自旋状态。进一步，动力学慢化还影响动力学测量，例如，传统的测量的方法 两难境地：要测量准确，需要较大 但当大，遇到临界慢化困难。Cluster算法过去一、二十年，如何改良MonteCarlo算法，克服临界慢化困难，一直是物理学家和计算科学家的一个重要课题。例如，Cluster算法Multi-canonical算法reweighting算法Hybrid算法MonteCarlo重整化群方法……Cluster算法的重要代表算法：Swendsen-Wang算法和Wolff算法 仍然以二维Ising模型为例,正那么分布写为 Wolff算法 从某一初始态出发随机取一格点i，翻转的值，即取如果近邻的自旋与不同值〔即与同值〕，以概率翻转该自旋，这样得到的与相同值的自旋称之为一个Cluster,如果近邻的自旋和同值，不做任何事情。 〔iii〕对该Cluster的近邻自旋操作〔ii〕直至停止。证明：各态历经自然满足，因为不管周围的自旋构形如何，随机选中的自旋总被翻转，所以，关键是细致平衡条件。→

→连起来的键 —— 对断开的键 对不做任何事情的键 对 随机抽到和的概率一