howking的个人博客分享 http://blog.sciencenet.cn/u/howking

博文

2维伊辛模型 数值模拟

已有 577 次阅读 2019-10-3 09:41 |个人分类:计算物理|系统分类:科研笔记| 计算物理、凝聚态物理


2维伊辛模型 数值模拟


2维伊辛模型 数值模拟


 等 70 人赞同了该文章


算法
♦ 选择一个格点i, 其自旋将考虑作翻转 si =-si.

♦ 计算与此翻转相联系的能量变化 dH.

♦ 计算这一翻转的转移几率 w.

♦ 产生一在 [0,1] 之间均匀分布的随机数 x.

♦ 如果x<w, 则翻转该自旋, 否则, 保持不变. 不论何种情况, 其结果都作为一新的状态.

♦ 分析该状态, 为计算平均值收集数据.


以下代码在2D中模拟Ising模型。

能量,磁化,比热和磁化率已被计算和绘制。


from __future__ import divisionimport numpy as npfrom numpy.random import randimport matplotlib.pyplot as plt#----------------------------------------------------------------------def initialstate(N):   
    ''' generates a random spin configuration for initial condition'''
    state = 2*np.random.randint(2, size=(N,N))-1
    return statedef mcmove(config, beta):
    '''Monte Carlo move using Metropolis algorithm '''
    for i in range(N):
        for j in range(N):
                a = np.random.randint(0, N)
                b = np.random.randint(0, N)
                s =  config[a, b]
                nb = config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N]
                cost = 2*s*nb
                if cost < 0:
                    s *= -1
                elif rand() < np.exp(-cost*beta):
                    s *= -1
                config[a, b] = s
    return configdef calcEnergy(config):
    '''Energy of a given configuration'''
    energy = 0
    for i in range(len(config)):
        for j in range(len(config)):
            S = config[i,j]
            nb = config[(i+1)%N, j] + config[i,(j+1)%N] + config[(i-1)%N, j] + config[i,(j-1)%N]
            energy += -nb*S
    return energy/4.def calcMag(config):
    '''Magnetization of a given configuration'''
    mag = np.sum(config)
    return magnt      = 2**8        # 温度点数量N       = 2**4        # 点阵尺寸, N x NeqSteps = 2**10       # MC方法平衡步数mcSteps = 2**10       # MC方法计算步数n1, n2  = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)tm = 2.269;    T=np.random.normal(tm, .64, nt)T  = T[(T>1.2) & (T<3.8)];    nt = np.size(T)Energy       = np.zeros(nt);   Magnetization  = np.zeros(nt)SpecificHeat = np.zeros(nt);   Susceptibility = np.zeros(nt)---------------------------------------------------------------------for m in range(len(T)):
    E1 = M1 = E2 = M2 = 0
    config = initialstate(N)
    iT=1.0/T[m]; iT2=iT*iT;
    
    for i in range(eqSteps):         # equilibrate
        mcmove(config, iT)           # Monte Carlo moves

    for i in range(mcSteps):
        mcmove(config, iT)           
        Ene = calcEnergy(config)     # calculate the energy
        Mag = calcMag(config)        # calculate the magnetisation

        E1 = E1 + Ene
        M1 = M1 + Mag
        M2 = M2 + Mag*Mag 
        E2 = E2 + Ene*Ene

        Energy[m]         = n1*E1
        Magnetization[m]  = n1*M1
        SpecificHeat[m]   = (n1*E2 - n2*E1*E1)*iT2
        Susceptibility[m] = (n1*M2 - n2*M1*M1)*iTf = plt.figure(figsize=(18, 10)); # plot the calculated values    sp =  f.add_subplot(2, 2, 1 );plt.plot(T, Energy, 'o', color="#A60628");plt.xlabel("Temperature (T)", fontsize=20);plt.ylabel("Energy ", fontsize=20);sp =  f.add_subplot(2, 2, 2 );plt.plot(T, abs(Magnetization), 'o', color="#348ABD");plt.xlabel("Temperature (T)", fontsize=20);plt.ylabel("Magnetization ", fontsize=20);sp =  f.add_subplot(2, 2, 3 );plt.plot(T, SpecificHeat, 'o', color="#A60628");plt.xlabel("Temperature (T)", fontsize=20);plt.ylabel("Specific Heat ", fontsize=20);sp =  f.add_subplot(2, 2, 4 );plt.plot(T, Susceptibility, 'o', color="#348ABD");plt.xlabel("Temperature (T)", fontsize=20);plt.ylabel("Susceptibility", fontsize=20);

计算能量差是最费时的工作,对于Ising 模型, 由于能量差只能取很少几个数值, 我们可以予先算好存起来以节省计算量. 这一技巧不仅适用于Ising 模型, 也适用于其它分立变量的模型如

Potts 模型等.

例如, 对于正方格子上的Ising模型, 每个自旋(设为s0)有四个近邻, 设为s1, s2, s3, s4, 则s0与近邻的相互作用能量为:

翻转一个自旋可能的能量差为

于是

求出dE0 和exp(-\beta d E0)的几个数值并保存下来,可以节省很多计算时间.



yang元祐:Ising 模型zhuanlan.zhihu.com图标





http://blog.sciencenet.cn/blog-566204-1200447.html

上一篇:爱因斯坦求和约定
下一篇:[转载]从密度泛函微扰理论(DFPT)分析声子及其晶体性质 -Stefano Baroni (1)

1 王安良

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-12-6 18:49

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部