蒙特卡洛模拟的前提（复式格子CrI3的居里温度蒙特卡洛模拟）
蒙特卡洛模拟的前提（复式格子CrI3的居里温度蒙特卡洛模拟）86-98行是计算平均磁矩，这里注意最后结果取平均后还要取绝对值，不然会出现整体翻转的现象。MC函数中是主要的蒙特卡洛方法，这里随机选一个点，先判断是A还是B，然后对他进行自旋翻转，然后再进行能量判断，是否翻转这个点。代码import mathimport random######init set###nx=30ny=30eKB=11610mca=100000mcb=100000J1= 0.000675TBEG=30 TEND=160step=5spinA=[[3]*nx for i in range(ny)]spinB=[[3]*nx for i in range(ny)]#initdef MC(T): global spinA global spinB x = random.randint(0 nx-1) y = random.randint(0 ny-1) xa=x-1 xb=x 1
preface
之前介绍过蒙特卡洛方法的思想，其实就是随机过程，利用统计概率的规律采样取样。但是之前只是介绍了简单晶格的Ising model，这一节1以CrI3为例，介绍一下复式格子的Ising model如何写代码。
代码以及解释
对于复格子，我们要构建多个数组（每一个格点构建一个数组）。如下图，我们只考虑Cr原子，这里两个不等价的格点分别用紫色跟蓝色表示。只考虑最近邻效应。
代码
import mathimport random######init set###nx=30ny=30eKB=11610mca=100000mcb=100000J1= 0.000675TBEG=30 TEND=160step=5spinA=[[3]*nx for i in range(ny)]spinB=[[3]*nx for i in range(ny)]#initdef MC(T): global spinA global spinB x = random.randint(0 nx-1) y = random.randint(0 ny-1) xa=x-1 xb=x 1 ya=y-1 yb=y 1 if x==0: xa=nx-1 if x==nx-1: xb=0 if y==0: ya=ny-1 if y==ny-1: yb=0 lat_choose = random.random if lat_choose <= 0.5: Ej = spinB[xa][y] spinB[xb][y] spinB[x][yb] dE = 2 * Ej * spinA[x][y] * J1 if lat_choose > 0.5: Ej = spinA[xa][y] spinA[xb][y] spinA[x][ya] dE = 2 * Ej * spinB[x][y] * J1 if dE<0 and lat_choose <= 0.5: spinA[x][y]=-spinA[x][y] if dE<0 and lat_choose > 0.5: spinB[x][y]=-spinB[x][y] if dE>0: r=random.random if r < math.exp(-dE * eKB/T): if lat_choose <= 0.5: spinA[x][y]=-spinA[x][y] if lat_choose > 0.5: spinB[x][y]=-spinB[x][y]def calE: global E E=0 for i in range(nx): for j in range(ny): xa=i-1 xb=i 1 ya=j-1 yb=j 1 if i==0: xa=nx-1 if i==nx-1: xb=0 if j==0: ya=ny-1 if j==ny-1: yb=0 Eja= (spinB[xa][j] spinB[xb][j] spinB[i][yb]) * spinA[i][j] Ejb= (spinA[xa][j] spinA[xb][j] spinB[i][ya]) * spinB[i][j] H=-(Eja Ejb)*J1 E=E H E=E/2.0def cals: global Mi Mi=0 for i in range(nx): for j in range(ny): Mi=Mi spinA[i][j] spinB[i][j] Mi=Mi/(nx*ny*2)########### t-M-calculation ##################################for t in range(TBEG TEND step): M=0 Mi=0 for j in range(mca): MC(t) for k in range(mcb): MC(t) cals M=Mi M M=M/mcb M=abs(M) print(t M)###########t-Cv-calculation###################################for t in range(TBEG TEND step):# E=0# E_total=0# E_2total=0# for j in range(mca):# MC(t)# for j in range(mcb):# MC(t)# calE# E_total=E_total E # E_2total=E*E E_2total# E_total=E_total/(mcb 0.0)# E_2total=E_2total/(mcb 0.0)# Cv=(E_2total-E_total*E_total)/(t*t 0.0)# print(t Cv)
逻辑思路
第15-16行定义两个数组，就是给两个格点A、B的自旋。
MC函数中是主要的蒙特卡洛方法，这里随机选一个点，先判断是A还是B，然后对他进行自旋翻转，然后再进行能量判断，是否翻转这个点。
86-98行是计算平均磁矩，这里注意最后结果取平均后还要取绝对值，不然会出现整体翻转的现象。
100-114是用来统计比热的，这里注释掉了，如果读者还要看比热曲线，去掉注释就可以。
输出结果
30 3.035 2.99998996666659940 2.99988863333355445 2.999881633333343350 2.998778733335073355 2.99712720000251960 2.99341996666845165 2.988127566668197270 2.970977300000532475 2.961977000000336780 2.934361966667259685 2.901108866666506490 2.8536209000003195 2.8057225333332316100 2.614236533332927105 2.4220778333332262110 1.9734956666666663115 1.4144269999999945120 0.8720862666666452125 0.6744532333333494130 0.029764500000001286135 0.06705653333332263140 0.19025083333333828145 0.16574436666666764150 0.11107360000000503155 0.031159566666671294
分析
对于”磁矩-温度“曲线中，磁矩下降最快的点就是居里温度点。当然也可通过看比热的拐点来判断（比热更准确）。对于反铁磁，奈尔温度只能通过比热（或磁阻）来确定。
这里可以看到比实验值大很多，因为Ising model的相空间只有up down，少了很多自由度，而且从Ising model的形式来看，默认其MAE为无穷大了。（对于这一点有问题的读者可以私信计算凝聚态物理公众号的小编，点赞后会回复^_^）
