二维铁磁单层材料CrI3以及Cr2Ge2Te6的发现，进一步激发了大家寻找室温的二维铁磁材料的兴趣。在第一性原理研究磁性材料的过程中，常常需要对材料的居里温度进行评估。通常我们用经典自旋的海森堡模型描述材料的电子自旋相互作用，有效哈密顿量写为（这里以CrI3为例，各向异性轴为z轴，ref[3]）：
求解这个模型我们常采用基于Metropolis算法的Monte Carlo模拟。
本文具体介绍了如何利用Metropolis算法基于上述哈密顿量求得居里温度。更重要的是，在这里介绍2019年由J D Alzate-Cardona等人发展的自适应算法（adaptive algorithm）(ref [1])，该方法比起常用的随机取样，可以有效控制取样的接受率在0.5左右，从而大大缩短收敛时间，减少达到热平衡所需要的Monte Carlo 步数 （MC steps）（大约仅仅需要原来的1/10甚至更少,详细见下文）。
工具：python3, matlab
（python代码请在Linux系统上运行，或者注销掉包含os.system()的命令行，任何问题请下方留言）
代码：链接：https://pan.baidu.com/s/1tQTYXTo_UiaH7-qZq_UiVw
提取码：m3at
1.1 Metropolis算法
Metropolis算法详细地介绍请查看ref [2], 这里只贴一个流程图。简单来说就是首先根据材料结构定义一套格子，每个格子上有一个自旋。每一个MC step就是首先在自旋空间中取样，取样后通过上述哈密顿量求得体系取样前后能量，若能量减小，用取样的自旋代替原来格子上的自旋。否则进行（6）的判断，反复进行这个过程，我们就可以说体系在温度 T 下达到热平衡，我们得到了体系在温度 T 下的系综（数学证明超出本文范畴）。这样通过计算体系的平均磁化（Magnetizaiton）和比热（
），我们可以得到体系的居里转变温度。
1.2 利用数组描述自旋
我们以CrI3为例，如图，CrI3每个晶胞有2个原子，我们标记为0和1，同样的我们需要标记每个原胞，用（ii，jj）表示。这样，我们可以定义第（ii，jj）原胞上的第kk个原子具有的自旋为
=
(
)，为了避免麻烦，我们统一设自旋矢量的长度为1 （注意此时的J的值相应改变）。为了描写这个自旋系统，我们就需要定义一个四维数组（4d array）data_spin = np.zeros([L,L,atom_num,3])。
我们只考虑最近邻的贡献，很容易发现在六角晶格中第0个自旋和1个自旋的最近邻原子描述方法不同。（ii,jj,0）原子的最近邻有三个，为（ii,jj,1）,(ii+1,jj,1),(ii,jj-1,1);（ii,jj,1）原子的最近邻有三个，为（ii,jj,0）,(ii-1,jj,0),(ii,jj+1,0)。注意考虑周期性边界条件，详细见附的python代码
这种方法可以很容易写出每个原子的最近邻，次近邻等等，而且对所有格子都适用。不要采用ref[2]中的搞成一维数组的考虑方法，非常麻烦，而且计算速度并不会快多少，推测为ref[2]成书时计算机内存有限的无奈之举。
2.1 自适应算法（adaptive algorithm）取样
现在来到本文的重点，我们注意到每个MC step的第一步就是取样，假设我们对体系所有N=atom_num*L*L个格点遍历了一遍，根据1.1中的流程，不是每次取样都能代替原来格点上的自旋，定义接受率（acceptance ratio，AR）为被代替格点的数目除以N。对于Ising模型，自旋仅有2个取向，很容易达到接受率为50%，所以Ising模型的MC模拟十分高效，能够很快达到热平衡。但是对于海森堡自旋模型，我们的S取样是连续的，这就造成接受率特别低（CrI3，T=1K下，一般为0.0几）。这样为了达到热平衡，我们往往需要很高的MC steps。
于是有人提出高斯取样方法，原理也很简单，就是在原有自旋方向周围取点，取点概率为正态分布。取点公式为：
为取样前第i个格点上自旋，
用来调节正态分布的宽度，它的值越大说明圆锥越宽，越接近随机取样。
为3*1的数组（矢量），每个分量都由标准正态分布随机数函数np.random.randn()产生。这样可以取点的范围由整个球变为一个圆锥，如下图。
显然
是随着体系变化的函数（与温度T，参数J，晶格种类都有关），合适的
可以提升MC的接受率，达到加快收敛的目的。
ref[1]的自适应算法的妙处就在于，它能随着MC step调节
的大小，使得接受率R总是在0.5附近。它定义了一个factor，
. 其中R为每次遍历一遍体系统计出的接受率。下一次MC step遍历开始时，改变
（详细看ref[1]和代码CrI3_
MC_impro.py）
我们设定L = 20,温度T=1，初始每个格点自旋方向随机选取。比较两种取样的收敛情况。如图，左边的维随机取样，遍历了体系20000次，从磁化强度看仍然没有收敛，看体系的能量，左边大约为-12左右，右边为-12.5能量更低。注意，同等条件下，右边的自适应算法取样只遍历了2000次，收敛效果却大大超出左边。我们看左边的接受率普遍在1e-2量级，右边却稳定在了0.5，这就是其高收敛速度的根本原因。
我们用自适应取样算法计算体系磁化强度和比热随温度变化曲线为下图，得到居里温度为 50 K左右 ，符合文献ref[3]的结果。其中L=25，初始状态采用0温下的基态（有序铁磁态），我们选择：预热2e4次遍历，总共2.2e4次遍历。
我们看到即使经过这么多次遍历，得到的随温度变化的曲线仍然很不光滑,似乎说明尽管通过上述讨论我们知道自适应算法取样可以快速达到热平衡，但是曲线光滑程度与体系还是与热扰动有关，也许增大L，加大平衡后求平均值得步数（这里是2000步）可以使得效果更好，读者可使用代码自己尝试，欢迎在下方留言。
另外，我们也应该看到自适应取样厉害的地方，下图磁化M vs. T图，考虑磁各向异性前后，随机取样的曲线形状居然没多大变化，这与自适应取样不同。形状的变化表示在低温下，由于z方向的磁各项异性，自旋更喜欢在z轴附近扰动，因此曲线变化应该更加平缓，右图反映了这个过程。
（下图中我都是从初始温度T=1e-4,np.linspace(1e-4,150,51)设定模拟温度，这似乎导致初始温度热容数据变得很奇怪，画图时我扣除了这个点。另外也许从高温降温到低温模拟效果会更好，有待尝试）
reference:
[1] J D Alzate-Cardona etc. Optimal phase space sampling for Monte Carlo simulations of Heisenberg spin systems J. Phys.: Condens. Matter 31 (2019) 095802 (10pp)
[2] M. E. J. NEWMAN Monte Carlo Methods in Statistical Physics
[3] Huang, C. et al. Toward Intrinsic Room-Temperature Ferromagnetism in Two-Dimensional Semiconductors. Journal of the American Chemical Society, doi:10.1021/jacs.8b07879 (2018).
