Python 方格子Ising模型模拟

使用最简单的逐点翻转法逐点模拟

import time
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt

def init_state(N):   ''' generates a random spin configuration for initial condition'''# np.random.randint(2, size=(N,N)) 表示产生N*N随机数,返回为列表state = 2*np.random.randint(2, size=(N,N))-1return state

def flipping(grid, beta):'''Monte Carlo move using Metropolis algorithm '''N = len(grid)for i in range(N):for j in range(N):# 产生0到N-1的随机数a = np.random.randint(0, N)b = np.random.randint(0, N)s = grid[a][b]E = grid[(a+1)%N][b] + grid[a][(b+1)%N] + grid[(a-1)%N][b] + grid[a][(b-1)%N]cost = 2*s*E# 如果能量降低接受翻转if cost < 0:s *= -1# 在0到1产生随机数，如果概率小于exp(-E/(kT))翻转elif rand() < np.exp(-cost*beta):s *= -1grid[a][b] = sreturn grid

def calculate_energy(grid):'''Energy of a given configuration'''energy = 0N = len(grid)for i in range(N):for j in range(N):S = grid[i][j]E = grid[(i+1)%N][j] + grid[i][(j+1)%N] + grid[(i-1)%N][j] + grid[i][(j-1)%N]# 负号来自哈密顿量energy += -E*S# 最近邻4个格点return energy/4def calculate_magnetic(grid):'''Magnetization of a given configuration'''mag = np.sum(grid)return mag

参数设置

nt      = 2**8        # 温度取点数量
N       = 2**4        # 点阵尺寸, N x N
eqSteps = 2**10       # MC方法平衡步数
mcSteps = 2**10       # MC方法计算步数Energy = [] # 内能
Magnetization  = [] # 磁矩
SpecificHeat = []  # 比热容/温度的涨落
Susceptibility = [] # 磁化率/磁矩的涨落T= np.linspace(1.2,3.8,nt)
T =list(T)n1 = 1/mcSteps*N*N
n2 = 1/mcSteps**2*N*N

开始模拟

time_start=time.time()
j = 0
for t in T:# 初始构型E = 0 # 每一温度下的能量M = 0 # 每一温度下的磁矩cv = 0 # 每一温度下的热容k = 0 # 每一温度下的磁化率config = init_state(N)# 热浴，达到平衡态for i in range(eqSteps):flipping(config, 1/t)# 抽样计算for i in range(mcSteps):flipping(config, 1/t)e = calculate_energy(config)m= calculate_magnetic(config)E += eM += mcv += e*ek += m*mCv = (n1*cv - n2*E*E)/(t*t)K = (n1*k - n2*M*M)/tEnergy.append(E*n1)Magnetization.append(M*n1)SpecificHeat.append(Cv)Susceptibility.append(K/(mcSteps**2*N*N))j +=1if j%10==0:print("已完成第%d步模拟" %j)
time_end=time.time()
print('totally cost',time_end-time_start)

处理画图

Magnetization = np.array(Magnetization)f = plt.figure(figsize=(18, 10)); # plot the calculated values    sp =  f.add_subplot(2, 2, 1 );
plt.plot(T, Energy, 'o', color="red");
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="blue");
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Magnetization ", fontsize=20);sp =  f.add_subplot(2, 2, 3 );
plt.plot(T, SpecificHeat, 'o', color="red");
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="blue");
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Susceptibility", fontsize=20);

Python 方格子Ising模型模拟相关推荐

二维方格子Ising模型代码
#include <iostream> #include <cstdlib> #include <cmath> #include <fstream> # ...
只显示小方格_不妨谈谈二维方格子吧
(想借该模型讲清直积态以及TB Model哈密顿量) 入手复杂的事物之前,从手算几个简单的例子开始是最好的. 所以,不妨谈谈二维方格子吧. 一个复杂的晶体,其内部的具有复杂的元胞(最小重复单元),每 ...
给定一个8*8的方格子，A点到B点的最短路径有多少条？
题目:给定一个8*8的方格子,如下图所示,求A点到B点的最短路径有多少条?用算法实现.(回溯法) 广度优先搜索只能找出一条最短路径答:从图中可以看出,A点到B点的最短路径为16,即A点横走8小格,纵 ...
有一个方格子，A点在左下角，B点在右上角，求A点到B点的最短路径数量
博主的第十篇博客,特此标记纪念. 概述有一个方格子,A点在左下角,B点在右上角,求A点到B点的最短路径有多少条. 思路我们这里假设方格子是8*8的方格.题目有时会混淆,A和B有没有占用左下角和右上 ...
一个8*8的方格子，A点在左下角，B点在右上角，求A点到B点的最短路径有多少条
题目: 一个8*8的方格子,A点在左下角,B点在右上角,求A点到B点的最短路径有多少条偶然看到一个求最短路径的题目, 感觉有趣所以分析了下, 此类问题不外乎递归排列.或者有规律的数列. 若只关 ...
Excel 方格子插件、DIY工具箱
小编通过使用这两款插件,极大的方便了Excel的操作,特意在此做一分享. 插件图示: DIY工具箱: 方格子插件: 下载地址:https://download.csdn.net/download/qq ...
深入理解python.md_跟黄哥学python序列文章之python方法链(method chaining)
跟黄哥学python序列文章之python方法链(method chaining) 写这篇文章来由,有朋友说下面这样的代码看不懂. choice = raw_input("please in ...
Python数字图像处理---1.1图像的像素格式与图像读写
目录前言图像像素格式图像读写前言本专栏面向所有希望或有兴趣从事数字图像处理工作.学习或研究的朋友,编程语言采用了当下最火的Python语言. Python是一种跨平台的计算机设计语言,也是一 ...
python可视化的图表汉字显示成框框_数据可视化——Matplotlib输出中文显示问题...
写在前面在学习可视化过程中,Matplotlib是其余Python可视化工具包的基础,是它们的老祖宗. Matplotlib是一个用于绘制高质量图形的Python第三方包,一般将其简写成mpl(ma ...

Python 方格子Ising模型模拟

Python 方格子Ising模型模拟相关推荐

最新文章

热门文章