CN109190247B - 优化的量子蒙特-卡洛模拟方法在研究复杂磁系统中的应用 - Google Patents

Abstract

本发明提出的优化的量子蒙特‑卡洛方法，用于研究磁性材料，属于计算物理和计算材料学领域。它包括：S1、系统哈密顿量中的自旋或磁矩为量子力学算符，所有物理量都按照量子理论算出；S2、模拟初始时令全部自旋随机取向；S3、每一步中都随机选取一个自旋，令其旋转一个随机的立体角度；S4、根据Metropolis算法判断自旋的新取向是否被接受；S5、如被接受，则更新近邻的能态；S6、判断本次循环是否结束，否则返回S3；S7、判断本次循环是否满足收敛条件或循环数大于某个整数，否则返回S3；S8、计算和输出磁结构和其他物理量。因量子理论的应用和算法的优化，大大提高了计算速度，且克服了经典模拟方法无法解决的困难。

Description

优化的量子蒙特-卡洛模拟方法在研究复杂磁系统中的应用

技术领域

本发明属于计算物理和计算材料学领域，其理论、模型、方法及开发出的相关软件将用于模拟磁系统的微观磁结构，计算其磁化强度、磁化率、磁滞回线、热容量等宏观物理量。

背景技术

数十年来，国内外研究者普遍运用经典蒙特-卡洛(Classical Monte Carlo)[1-5]和微磁学(Micromagnetics)[6-9]两种数值方法模拟磁性材料的微观磁结构，研究它们的宏观物理性质。

然而，这二种方法都建立于经典物理基础之上，即在模拟中，磁系统中的自旋和磁矩被当做长度不变、但可在空间旋转的经典矢量。显然，如此简单的处理方法，将有碍于对磁系统的精确描述，也将对模拟的运算速度和结果产生不良的影响。

在某温度下，为使模拟收敛于系统的平衡态，经典蒙特-卡洛法(CMC)在运算中不断旋转各自旋的空间取向，比较旋转前后系统的能量变化，根据Metropolis算法，确定自旋的新状态是否被接受，以降低系统的总能量；在运行数万甚至数百万次循环后，在最后的数万次循环中，对系统中每个自旋的矢量值求平均，作为平衡态中各个自旋的矢量值；然后再计算整个系统的磁化强度、磁化率、热容量等宏观物理量。

微磁学模拟法通常把磁体分成许多网格，每个网格的总磁矩为

开假设在整个磁体中，所有

的大小M_S处处相等，求解耦合的Landau-Lifshitz-Gilbert微分方程，以确定

随时间的变化规律[6-9]。但是，微磁学通常不考虑温度效应。为了克服这一缺陷，Skomski等人把磁矩的大小M_S以及系统参量K₁等作为与温度相关的量[8]。但是，如何确定M_S(T)和K₁(T)二函数又成了新的问题。而且，如纳米等有限大小的系统，各处的磁矩大小显然不同，但微磁学却忽略了这种差异。

此外，微磁学建立在经典的连续模型之上，它成立的条件是，系统的磁性质在空间缓慢变化，可以近似视作连续变化[10]。当斯格明子的直径大于数十、甚至数百纳米时，连续模型近似成立，微磁学能够算出较为满意的磁结构。但是近年来，欧洲的科学家们发现，在多层磁材料的界面上形成的斯格明子的直径仅为几个纳米[11-14]。这些小尺度的斯格明子具有极高的数据存储价值。在此情况下，经典的连续模型不再适用，微磁学的描述不再精确和可靠[10]。

经典蒙特-卡洛和微磁学二方法的上述不足可能会大大影响收敛速度和模拟结果的正确性，可以发现一些文献中用它们算出的磁结构并不具有系统的对称性。例如，考虑圆型纳米盘上磁偶极矩之间的长程相互作用，用二经典方法算出的磁结构都不对称；再如，采用经典蒙特-卡洛法无法模拟出非零温度下无限二维正方结构上的反铁磁斯格明子晶体。

根据量子理论，A.W.Sandvik等人提出了模拟简单自旋系统的量子蒙特-卡洛方法[15，16]。运用此法需计算系统的配分函数的路径积分，理论上十分缜密，也因此甚为繁杂，故仅被用于模拟自旋S＝1/2的、低维的、简单的自旋系统。然而，实际的磁系统往往由多种原子构成，它们的自旋可取不同的、更大的量值，磁系统内部也会存在多种复杂的相互作用。所以这一量子蒙特-卡洛法与实际运用还有较远的距离。

为克服上述方法的不足，本人近年来研发出二种量子模拟方法。在此二种模型中，磁系统哈密顿量中的角动量或磁矩是量子力学算符，不再是经典矢量；任意温度下各种物理量，如磁矩、磁化强度、系统的总能量、总自由能、热容量等，都严格按照量子理论算出。

第一种量子方法运用自洽算法(Self-Consistent Algorithm)，故简记为SCA法。由于量子理论的引入，程序能自发地收敛于系统的平衡态，而不必像CMC那样每一步都对系统的总能量求极小，因而显得十分简捷。运用SCA方法，笔者已经成功地模拟了纳米颗粒、纳米线、纳米盘等系统的复杂的磁结构，如纳米盘上的铁磁和反铁磁涡旋[17-20]。

第二种为改进的量子蒙特-卡洛模拟方法。它将量子磁学与Metropolis算法相结合，简记为QMC方法。QMC的要点为：系统哈密顿中的自旋为量子力学算符；所有物理量都严格按照量子理论算出；每一步都如经典蒙特-卡洛方法那样，根据Metropolis法则确定所选自旋的新状态是否被接受；经过数千次循环后，程序即可能收敛于平衡态附近。至今，笔者已用此量子方法(QMC)模拟了纳米颗粒和纳米盘等，其模拟结果与前一量子方法(SCA)的结果相符[21,22]。

在CMC和QMC模拟中，为确定自旋的新状态是否被接受，需要计算所选自旋的旋转引起的系统总能量的变化[4]。如果所研究的磁系统甚小，或者需计及磁偶极之间的长程相互作用，直接计算系统总能量的变化似乎可行。但是，如磁系统中包含大量的自旋，如每一步都计算系统总能量的变化，将浪费大量时间。实际上，所选自旋的旋转引起的系统总能量的变化仅限于局域。因受Ising模型的影响，人们在程序设计中往往仅仅注意到了这一点，但却忽视了另一重要的事实，即该自旋的转动也引起局域内其它自旋的能态的变化。这一严重的疏忽将导致计算收敛缓慢、无法收敛、甚至收敛于不正确的系统状态。

对于简单的磁系统，如忽视局域自旋状态的更新，在经过数千次循环后，对最后的数百次循环求平均，来确定系统的“平衡状态”，或采用其它的补救措施，也可能模拟出较为合理的磁结构，如笔者运用改进的QMC模拟纳米球和纳米盘所做的那样[21,22]。但这些权宜之计大大减缓了计算速度，而且，笔者在近期的研究中发现，对于含有复杂相互作用(如Dzyaloshinsky-Moriya作用)的磁系统，权宜的方法却未必能保证模拟出正确的磁结构。

本发明就是基于上述分析，对笔者提出的量子蒙特-卡洛方法作进一步的优化，并应用于模拟铁磁和反铁磁的二维Bloch型和Néel型斯格明子晶体，以论证新法的有效性和正确性，为其广泛的应用奠定坚实的基础。

参考文献

[1]KurtBinder，Dieter Heermann，Monte Carlo Simulation in StatisticalPhysics，An Introduction，5-th edition.Springer，2010.

[2]H.D.Rosales，D.C.Cabra，and P.Pujol，Three-sublatticeskyrmion crystalin the antiferromagnetic triangular lattice，Phys.Rev.B 92(2015)214439.

[3]R.Keesman，M.Raaijmakers，A.E.Baerends，G.T.Barkema，and R.A.Duine，Skyrmions in square-lattice antiferromagnets，Phys.Rev.B 94(2016)054402.

[4]D.Hinzke，U.Newak，Monte Carlo Simulation of the MagnetizationSwitching in a Heisenberg Model for Small Ferromagnetic Particles，ComputerPhysics Communications，121-122(1999)334.

[5]E.Y.Vedmedenko，H.P.Oepen，A.Ghazali，J.-C.S.Lévy，and J.Kirschner，Magnetic Microstructure of the Spin Reorientation Transition：A ComputerExperiment，Phys.Rev.Lett.84(2000)5884.

[6]W.F.Brown，Micromagnetics，New York：Interscience，1963.

[7]R.H.Kodama，A E Beikwitz，Atomic Scale Magnetic Modelling of OxideNanoparticles，Phys.Rev.B 59(1999)6321.

[8]R.Skomski，P.Kumar，George C.Hadjipanayis，and D.J.Sellmyer，Finite-Temperature Micromagnetism，IEEE Transactions on Magnetics 49(2013)3229.

[9]Ph.Depondt，J.-C.S.Lévy，F.G.Mertens，Vortex polarity in 2-D magneticdots by Langevin dynamics simulations，Phys.Lett.A 375(2011)628-632.

[10]Magnetic skyrmion，https：//en.wikipedia.org/wiki/Magnetic skyrmion

[11]S.Heinze，K.von Bergmann，M.Menzel，J.Brede，A.Kubetzka，R.Wiesendanger，G.Bihlmayer，and S.Blügel，Nat.Phys.7(2011)713-718.

[12]B.Dupé，G.Bihlmayer，S.Blügel，and S.Heinze，Nat.Commun.7(2016)11779。

[13]J.Hagemeister，D.Iaia，E.Y.Vedmedenko，K.von Bergmann，A.Kubetzka，andR.Wiesendanger，Phys.Rev.Lett.117(2016)207202.

[14]J.Hagemeister，N.Romming，K.von Bergmann，E.Y.Vedmedenko，R.Wiesendanger，Nat.Commun.6(2015)8455.

[15]A.W.Sandvik，Phys.Rev.B 56(1997)11678.

[16]A.W.Sandvik，Computational Studies of Quantum Spin Systems，arXiv：1101.3281v1[cond-mat.str-el]17Jan 2011.

[17]Z.-S.Liu，M.

and V.

Magnetism of PrAl₂nanoparticleinvestigated with a quantum simulation model，J.Phys.：Condens.Matter 23(2011)016002.

[18]Z.-S.Liu，V.

and M.

Magnetism of DyNi₂B₂Cnano-particle investigated with a quantum simulation model，Phys.Status Solidi B，249(2012)202-208.

[19]Z.-S.Liu，H.Ian，Effects of Dzyaloshinsky-Moriya Interaction onMagnetism in Nanodisks from a Self-Cohsistent Approach，J.Nanopart.Res.(2016)18：9.

[20]Z.-S.Liu，H.Ian，Numerical Studies on Antiferromagnetic Skyrmionsin Nanodisks by Means of a New Quantum Simulation Approach，Chem.Phys.Lett.649(2016)135.

[21]Z.-S.Liu，V.

and M.

Mutual Verification of Two NewQuantum Simulation Approaches for Nanomagnets，Physica E 62(2014)123-127.

[22]Z.-S.Liu，O.Ciftja，X.-C.Zhang，Y.Zhou，H.Ian，Vortical Structures forNanomagnetic Memory Induced by Dipole-DipoleInteraction in Monolayer Disks，Superlattices & Microstructures，

[23]X.-Z.Yu，Y.Onose，N.Kanazawa，J.H.Park，J.H.Han，Y.Matsui，N.Nagaosa，Y.Tokura，Nature，465(2010)901

[24]A.O.Leonov，M.Mostovo，Multiply periodic states and isolatedskyrmions in an anisotropic frustrated magnet，Nat.Commun.6(2015)8275

发明内容

本发明的目的在于提供一种优化的量子蒙特-卡洛方法，用于模拟复杂磁系统的微观磁结构，研究其宏观物理性质。它克服了现存经典模拟方法的收敛缓慢、不正确收敛的严重困难。

本发明的实现包括以下步骤：

S1、对磁系统哈密度量进行量子化处理，即其中的自旋或磁矩为量子力学算符，而非经典矢量；

S2、模拟始于温度T＝T₀，T₀高于磁临界温度，故令所有的自旋量随机取向；

S3、模拟的每步中，都随机地选取一个自旋量，令其旋转一个随机的立体角度；

S4、根据Metropolis算法，判断该自旋的新取向是否被接受，如被接受，则更新局域中其他自旋的能态，并执行下一步，如不被接受，则直接执行下一步；

S5、判断本次循环是否结束，即模拟的步数是否等于自旋的总数，如是，则执行下一步，如否，则返回步骤S3；

S6、判断本次循环是否满足收敛条件或循环次数大于某个整数，如是，则执行下一步，否则，返回步骤S3；

S7、根据量子理论计算和输出微观磁结构和其他宏观物理量。

本发明的进一步技术方案还包括以下步骤：

S8、令T＝T-△T，即降低温度，△T为温度循环步长，如T低于预先设定的最低温度T_f，则计算结束，否则返回S3。

本发明的进一步技术方案是：所述步骤S4中还采取如下措施：

S4-1、如涉及边界自旋，则需运用周期性边界条件；

S4-2、运用量子理论计算自旋的平均值及能量等。

本发明的进一步技术方案是：所述优化的量子蒙特卡洛方法通过更新局域中其他自旋的能态，使得模拟更迅速地收敛于系统的平衡态。

本发明的又一技术特征是，最后一次循环计算得到的状态，即可精确地表示系统在此温度下的平衡态，无需像经典蒙特-卡洛方法那样，需要对最后的数万次循环求平均。

本发明的又一技术特征是：所述优化的量子蒙特卡洛方法，在任意温度下，都能细致地计算出材料中各点处磁矩的不同幅值及方向。

本发明的目的在于提供一种优化的量子蒙特-卡洛方法，用于模拟和研究磁性材料，比如实例中的铁磁和反铁磁的二维Bloch型和Néel型斯格明子晶体。

本发明的有益效果是：通过量子理论的应用和对近邻自旋态的及时更新，不仅大大提高了计算速度，而且克服了现存经典模拟方法的收敛缓慢、无法正确收敛的严重困难。

附图说明

图1是本发明实施例提供的优化的量子蒙特-卡洛模拟方法的流程图。

图2是本发明实施例提供的Bloch型二维铁磁斯格明子晶体示意图。

图3是本发明实施例提供的Bloch型二维反铁磁斯格明子晶体示意图一。

图4是本发明实施例提供的Bloch型二维反铁磁斯格明子晶体示意图二。

图5是本发明实施例提供的Néel型的二维铁磁斯格明子晶体示意图。

图6是本发明实施例提供的Néel型的二维反铁磁斯格明子晶体示意图。

具体实施方式

一、优化的Metropolis算法

假设磁系统中有N个自旋，模拟从高温T₀开始。在某一温度T下，循环的步骤如下：

1、随机地选取某一自旋S_i；

2、使S_i转过一随机的空间立体角度；

3、计算由此转动引起的系统能量变化ΔE_i；

4、计算p_i＝exp(-ΔE_i/k_BT)；

5、产生一个随机数r_i；

6、如果r_i小于p_i，则此操作被接收，且更新局域的自旋能态；

7、如果r_i大于p_i，则此操作被抛弃；

8、如上述步骤数小于N次，则返回(1)；

9、如不满足收敛条件及循环次数小于给定的整数，返回(1)；

10、根据量子理论，计算和输出微观磁结构及宏观物理性质。

在上述过程中，所有S_i和ΔE_i都需按照量子力学公式算出。特别是，笔者在运用数种不同方案，比较多次模拟结果后发现，步骤6中局域自旋状态的及时更新是实现优化、高效的重要措施。

二、存在Dzyaloshinsky-Moriya作用的二维磁系统

对于此类二维磁系统，如施加垂直外加磁场，则在表面或界面上形成斯格明子晶体。磁系统的哈密顿量可写为

其中第一项表示海森堡(Heisenberg)交换作用，第二项表示Dzyaloshinsky-Moriya作用(DMI)，第三项表示垂直于二维表面或界面的单轴各向异性，最后一项则表示系统与外加磁场的相互作用能量。J_ij和

分别表示相邻的第i个和第j个自旋间的海森堡相互作用和DMI相互作用的强度。如J_ij＞0，系统为铁磁耦合；如J_ij＜0，系统为反铁磁耦合。如果矢量

(

表示第i个自旋到第j个自旋的空间矢量)，则斯格明子为涡旋状的Bloch型；如矢量

(

为垂直于二维系统表面或界面的单位矢量)，则斯格明子为发散或者汇聚形状的Néel型。

三、存在DMI及Compass型各向异性作用的二维磁系统

此类二维磁系统内除了Heisenberg作用和DMI之外，还存在Compass型各向异性作用，其哈密顿量记为：

其中的最后二项即表示Compass型各向异性作用。此种系统不需外加磁场，斯格明子晶体即可在临界温度下自发形成。

四、优化的量子蒙特-卡洛方法(OQMC)的其它技术要素

本人在过去的研究工作已采取了以下措施：

(1)引入量子理论，磁系统哈密顿量中的自旋不再是经典矢量，而是量子力学算符，所有物理量都严格按照量子力学公式算出；

(2)模拟的每一步都采用Metropolis算法，以确定自旋的新状态是否被接受；

(3)对最后的数百次循环中算出的各自旋矢量求平均。

因(3)需要计算各自旋平均值，模拟不仅相当费时，且效果也不甚理想。故在此发明中，笔者采取了以下重要的优化措施：

(4)如所选自旋的新状态被接受，则立即更新局域中其它自旋的能态。

为了模拟二维系统上的斯格明子晶体，还需要考虑：

(5)周期性边界条件；

(6)垂直外加磁场；

(7)对于存在Compass型各向异性的磁系统，则不需要考虑外加磁场，系统便会在临界温度下自发形成斯格明子晶体。

采取(4)这一优化措施后，程序便能迅速地收敛于平衡态，即最后一次循环算得的状态就可精确地表示系统的平衡态，不必再对最后的许多次循环求平均，从而大大提高了运算速度。

为证明这一措施(4)的必要和正确性，笔者比较了以下数种方案：

(a)采用未优化的Metropolis算法；

(b)采用未优化的Metropolis算法，但在每一次循环结束后再更新所有自旋的能态。

(c)采用优化的Metropolis算法，即当所选自旋被旋转后，便立即更新局域内自旋的能态。

考虑这三个方案及其组合，笔者模拟了Bloch型的二维斯格明子晶格(设J＝1K，D＝1.02733K，B＝0.12Tesla)，结果是：

如仅考虑(a)方案，模拟不出斯格明子晶格。

如运用(b)方案，程序在温度T＞0.3K的温区内理想收敛，在低温区收敛甚差；在温区3.1≥T≥0.3K内模拟出较好的斯格明子晶格；但当T＝0.1K时，磁结构却为Skyrmions+Bimerons。

如运用(c)方案，程序在整个温区快速收敛，且在温度T≤3.1K的整个温区内模拟出十分对称、周期分布的斯格明子晶格。

如采用(b，c)方案的组合，模拟结果与单独采用(c)方案相同，仅仅斯格明子阵列相对地位移，且算出的系统总自由能等宏观热力学量，在整个温区相重合。

因此可以推断，方案(c)对于快速、正确的模拟是必不可少的。

图1示出本发明提供的优化的量子蒙特-卡洛方法的流程图，详述如下：

步骤S1，对磁系统的哈密顿量进行量子化处理，即其中的自旋或磁矩为量子力学算符；

步骤S2，模拟始于T＝T₀，T₀需高于磁临界温度，令系统中所有自旋的取向随机分布；

步骤S3，模拟的每一步中，都随机选取一个自旋量，令其旋转一个随机的空间立体角度；

步骤S4，根据Metropolis算法判断自旋的新取向是否被接受，如被接受，则更新局域中其他自旋的能态并执行下一步，如不被接受，则直接执行下一步；

步骤S5，判断本次循环是否结束，即模拟的步数是否等于系统中自旋的总数N，如是，则执行下一步，如否，则返回步骤S3；

步骤S6，判断本次循环是否满足收敛条件，即前后两次循环算出的所有自旋的矢量差值的相对变化是否小于某个给定的小量δ，或循环数是否大于某个给定的上限整数，如是，则执行下一步，如否，则返回步骤S3；

步骤S7，根据量子理论，计算和输出微观磁结构及其他宏观物理量。

步骤S8，降低温度，即令T＝T-△T，△T为温度的循环步长，如T低于预先设定的最低温度T_f，则计算结束，否则返回S3。

五、本发明的有益效果

1、本发明引入了量子理论，克服了数十年来国内外广泛应用的微磁学和蒙特-卡洛方法的经典局限性，是理论上的一大进步。

2、磁材料的性质缘于原子或离子之间的微观相互作用，而这些相互作用，又是由局域电子及传导电子的分布决定的，量子理论才能给出磁结构和磁性质的精确描述，因此，量子理论的引入和应用是非常必要的。

3、量子理论的应用，使得程序设计变得简单、易行，模拟十分快速。例如，当研究稀土磁性系统时，需要考虑复杂的晶体场作用(CEF)。如采用微磁学或者CMC，用经典矢量表示晶体场作用和描述磁结构的变化，就甚为不便。相反，如果用量子力学算符表示CEF，根据量子理论计算各物理量，程序设计就十分容易[18]。

再如，对于二维磁系统，国外学者仅能模拟出有限大小的正方磁系统中单个反铁磁斯格明子，却未能模拟出无限二维正方系统中的反铁磁斯格明子晶体。然而，笔者运用量子的OQMC方法，就能十分容易地算出非零温度下无限大二维系统的反铁磁Néel型和Bloch型斯格明子晶体。有关结果将在实例中给出。

此外，国外学者运用CMC和微磁学方法，难以模拟出有限大小纳米盘上强磁偶极作用产生的对称涡旋型磁结构，而用笔者的两种量子方法，都算出了十分对称的涡旋磁结构，且此二法的结果一致[22]。

4、本发明的OQMC法，保证了模拟的正确收敛。在实例四中，如在每次循环完成之后才对所有自旋的能态更新，则仅能在0.8K≤T≤1.9K的温区内模拟出Néel型的二维反铁磁斯格明子晶体，而在更低温区却不能算出正确的磁结构。

5、本发明的OQMC法，大大提高了计算速度。举例来说，如模拟28×28个自旋构成的二维正方系统上的反铁磁斯格明子晶体，采用了周期性边界条件，从T＝2.5K到0.1K共25个温度点，T表示的是折合的温度，算出临界温度下18幅反铁磁斯格明子晶体图，使用ThinkPad T470P笔记本，仅需要约11分钟。为了确保计算结果的正确性，模拟的计算精度要求很高(取δ＝10^-6)；如取稍低的精度，运算速度会更快。

6、本发明OQMC法的优点还表现为，每一温度下的最后一次循环，就能精确地给出系统的平衡态，而不必像经典蒙特-卡洛那样，需对最后的数万次循环求平均。

六、OQMC法与SCA法的区别

1、二法的共同之处

(1)两种量子模拟方法都采用量子理论，即都把系统哈密顿中的自旋或角动量当做量子力学算符，而不是经典矢量。

(2)所有物理量都严格按照量子理论算出，而不是像经典蒙特-卡洛法那样通过简单的代数平均算出。

(3)二法可分别应用于模拟同一磁系统，并给出一致的结果。比如四个应用实例，尽管二法考虑的温度不同，但结果非常相似。

2、二法的不同之处

(1)SCA法采用自洽算法，它不必对磁系统内的自旋进行人为的干预，程序就可以自发收敛于平衡态。

(2)OQMC法采用改进的Metropolis算法，每一步都需对随机选取的自旋进行旋转操作，判断它是否被接受，更新近邻能态等等，所以二者的算法完全不同。

实施例

以下通过模拟Bloch型和Néel型的二维铁磁和反铁磁斯格明子晶体，论述OQMC法的具体应用，证明其正确性和有效性。如前所述，这类二维系统的哈密顿量可写为：

式中，当J_ij＞0时，系统为铁磁耦合；当J_ij＜0时，系统为反铁磁耦合。如果

可模拟出Bloch型斯格明子；如果

则可模拟出Néel型斯格明子。为了简化模型，以下的四个应用中仅考虑最近邻自旋之间的相互作用，且假设它们处处相等，即J_ij＝J，D_ij＝D。

实施例一：Bloch型的二维铁磁斯格明子晶体(Ferromagnetic Skymion Crystalof Bloch-Type)

选用30×30的正方格子，每个格点上都有一个S＝1的自旋。为了模拟无穷二维系统，采用了周期性边界条件。此处未考虑单轴各向异性。再设J＝1K，令D＝1.02733K，于是根据相关理论得知，在弱的垂直外加磁场中斯格明子之间的周期距离λ＝10。

图2绘出B＝0.12Tesla，T＝1K时正方格子上的六角密排结构的铁磁斯格明子阵列，它具有完美的几何对称性，空间周期为10，与理论和实验的结果相符[23]。

这里，外加磁场沿着z方向，每个斯格明子中心区域的自旋磁矩的z分量都沿着-z方向，而其周边区域的自旋磁矩的z分量取向相反，所以都是严格的铁磁斯格明子。此外，每个斯格明子磁涡旋都是顺时针的，都被四个逆时针的磁涡旋包围着，从而降低系统的总能量，磁结构更加稳定。

实施例二：Bloch型的二维反铁磁斯格明子晶体(Antiferromagnetic SkymionicCrystal of Bloch-Type)

模拟中采用35×35的正方格子，每个格点上有一个S＝1的自旋。取J＝-1K，D＝1K(于是D/J＝-1)，同样忽略了单轴垂直各向异性作用。采用周期边界条件，以模拟无穷大二维系统。模拟中发现，当着垂直外加磁场强度满足3.9Tesla≤B≤4.1Tesla时，在T＜1.8K的低温区内诱发出反铁磁斯格明子阵列；略增大或减弱外加磁场，反铁磁斯格明子晶体便消失，被如简单的反铁磁结构取代。

图3、4示出在B＝4Tesla，T＝1K时的反铁磁斯格明子晶格。它具有完美的对称性，空间周期λ＝7。每个斯格明子中心区域内自旋的z分量平均值最小；而在平行于二对角线的方向上，二相邻斯格明子中心连线中点处的自旋的z分量值最大；相邻自旋磁矩的xy投影取向相反。这些都赋予二维斯格明子晶体的“反铁磁”特征。

实例一中B₁＝0.12Tesla，实例二中B₂＝4Tesla，D与J的比值D/J相近，但是B₂/B₁≈33.333。因此，要观测到反铁磁斯格明子阵列，需要施加33倍的强磁场。此外，反铁磁磁斯格明子晶格仅存在于很窄的B区间，因此稳定性差。这些都给实验上获得和观测反铁磁斯格明子晶体带来了极大的困难。

实施例三：Néel型的二维铁磁斯格明子晶格(Ferromagnetic Skymionic Crystalof Néel-Type)

选用大小为30×30的正方格子，每个格点上的自旋5＝1，用周期边界条件以模拟无穷大的二维平面系统。本例与实例一使用相同的参数，即J＝1K，D＝1.02733K，T＝1K，垂直外加磁场B＝0.12Tesla，模拟结果如图5所示。有趣的是，尽管二种斯格明子分别为Bloch型和Néel型，但二晶格的空间周期都为λ＝10，30×30的方格上也产生18个斯格明子，它们也形成正六角密排结构。在每个斯格明子的中心区域内，各自旋的z分量都沿着-z方向，与外磁场方向相反；其周边区域的各自旋z分量与外磁场方向相同，所以都是典型的斯格明子。二种磁结构的不同之处在于，这里的Néel型斯格明子内的自旋的xy分量都指向其中心。

实施例四：Néel型的二维反铁磁斯格明子晶格(Antiferromagnetic SkymionicCrystal of Néel-Type)

考虑大小为28×28的正方格子，每个格点上的自旋S＝1，采用周期性边界条件以模拟无穷大的二维平面系统。使用与实例二相同的参数，即J＝-1K，D＝1K，B＝4Tesla，模拟结果，如图6所示。与图3、4相比，尽管二实例分别为Bloch型和Néel型的二维反铁磁斯格明子晶体，但二晶体的空间周期都为λ＝7，斯格明子也形成正方结构。同样，图6中每个斯格明子中心处自旋z分量的平均值最小；而在平行于二对角线的方向上，二相邻斯格明子中心连线中点处的自旋的z分量值最大；且相邻自旋的xy投影取向相反。这些都使得斯格明子晶体具有“反铁磁”的特征。除边缘区域外，各斯格明子内的自旋的xy投影都沿着径向反向排列，所以属于Néel型。

为了模拟正方格子上的反铁磁斯格明子，R.Keesman等人令J＝-1，D/|J|＝1，H/|J|＝4，在8×8的有限大小的方格子上模拟出单个反铁磁斯格明子[3]。然而，Tretiakov等人的理论研究表明，此种斯格明子存在于非零温度，在外加磁场中会更加稳定[24]。所以，实例二、四充分说明，当着斯格明子的空间尺度很小时，经典理论的描述不再精确和可靠，必须采用量子理论和方法。

此外，图2与图5具有对偶性，图4也与图6具有对偶性，表明自然规律的完美对称性，进一步说明OQMC方法及其计算结果的正确性。

随着温度和外加磁场强度的变化，斯格明子及其晶格的形状、周期也随之变化。因篇幅所限，这里仅列出数张典型的代表图。

以上所述仅为本发明的数个实施例而已，并不用以限制本发明，本发明不仅可用于模拟有限大小的磁性纳米系统，而且可用于模拟更一般的二维、三维、几乎所有磁性材料。因此，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种优化的量子蒙特-卡洛方法，其特征在于，所述优化的量子蒙特卡-洛方法包括以下步骤：

S1、对磁系统哈密顿量进行量子化处理，即系统哈密顿量中的自旋或磁矩为量子力学算符；

S2、模拟始于高于磁临界点某一温度T＝T₀，故令系统中所有自旋的空间取向随机分布；

S3、模拟的每一步开始，都随机选取一个自旋，令其旋转一个随机的空间立体角度；

S4、根据Metropolis算法，判断此自旋的新取向是否被接受，如被接受，则更新局域中其他自旋的能态，并执行下一步，如不被接受，则直接执行下一步；

S5、判断本次循环是否结束，即模拟的步数是否等于系统中的自旋总数目，如是，则执行下一步，如否，则返回步骤S3；

S6、判断本次循环是否满足收敛条件或循环数是否大于某个给定的整数，如是，则执行下一步，如否，则返回步骤S3；

S7、根据量子理论计算和输出微观磁结构及磁系统的宏观物理量。

2.根据权利要求1所述的优化的量子蒙特-卡洛方法，其特征在于，该优化的量子蒙特-卡洛方法还包括以下步骤：

S8、令T＝T-△T，△T为温度循环步长,如T低于设定的最低温度值T_f，则计算结束，否则返回S3。

3.根据权利要求2所述的优化的量子蒙特-卡洛方法，其特征在于，所述步骤S4中还包括以下步骤：

S4-1、如涉及边界自旋，则需运用周期性边界条件；

S4-2、运用量子理论计算自旋的平均值及能量。

4.根据权利要求3所述的优化的量子蒙特-卡洛方法，其特征在于，通过每一步中更新局域中其他自旋的能态，使得模拟更贴近实际的物理过程，故能迅速收敛于系统的平衡态。

5.根据权利要求4所述的优化的量子蒙特-卡洛方法，其特征在于，在任意温度下，可以计算出系统中各处磁矩的不同幅值及方向。

6.根据权利要求1-5任一项所述的优化的量子蒙特-卡洛方法，其特征在于，本方法能很容易成功地模拟出二维Bloch型和Néel型的铁磁和反铁磁斯格明子晶体。

Images