伊辛模型 [Python]
Ising model [Python]
我正在尝试模拟 Barabasi-Albert 网络中的 Ising 相变,并尝试复制一些可观测值的结果,例如在 Ising 网格模拟中观察到的磁化强度和能量。但是,我在解释我的结果时遇到了麻烦:不确定物理学是否错误或实施中是否存在错误。这是一个最小的工作示例:
import numpy as np
import networkx as nx
import random
import math
## sim params
# coupling constant
J = 1.0 # ferromagnetic
# temperature range, in units of J/kT
t0 = 1.0
tn = 10.0
nt = 10.
T = np.linspace(t0, tn, nt)
# mc steps
steps = 1000
# generate BA network, 200 nodes with preferential attachment to 3rd node
G = nx.barabasi_albert_graph(200, 3)
# convert csr matrix to adjacency matrix, a_{ij}
adj_matrix = nx.adjacency_matrix(G)
top = adj_matrix.todense()
N = len(top)
# initialize spins in the network, ferromagnetic
def init(N):
return np.ones(N)
# calculate net magnetization
def netmag(state):
return np.sum(state)
# calculate net energy, E = \sum J *a_{ij} *s_i *s_j
def netenergy(N, state):
en = 0.
for i in range(N):
for j in range(N):
en += (-J)* top[i,j]*state[i]*state[j]
return en
# random sampling, metropolis local update
def montecarlo(state, N, beta, top):
# initialize difference in energy between E_{old} and E_{new}
delE = []
# pick a random source node
rsnode = np.random.randint(0,N)
# get the spin of this node
s2 = state[rsnode]
# calculate energy by summing up its interaction and append to delE
for tnode in range(N):
s1 = state[tnode]
delE.append(J * top[tnode, rsnode] *state[tnode]* state[rsnode])
# calculate probability of a flip
prob = math.exp(-np.sum(delE)*beta)
# if this probability is greater than rand[0,1] drawn from an uniform distribution, accept it
# else retain current state
if prob> random.random():
s2 *= -1
state[rsnode] = s2
return state
def simulate(N, top):
# initialize arrays for observables
magnetization = []
energy = []
specificheat = []
susceptibility = []
for count, t in enumerate(T):
# some temporary variables
e0 = m0 = e1 = m1 = 0.
print 't=', t
# initialize spin vector
state = init(N)
for i in range(steps):
montecarlo(state, N, 1/t, top)
mag = netmag(state)
ene = netenergy(N, state)
e0 = e0 + ene
m0 = m0 + mag
e1 = e0 + ene * ene
m1 = m0 + mag * mag
# calculate thermodynamic variables and append to initialized arrays
energy.append(e0/( steps * N))
magnetization.append( m0 / ( steps * N))
specificheat.append( e1/steps - e0*e0/(steps*steps) /(N* t * t))
susceptibility.append( m1/steps - m0*m0/(steps*steps) /(N* t *t))
print energy, magnetization, specificheat, susceptibility
plt.figure(1)
plt.plot(T, np.abs(magnetization), '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Average Magnetization per spin')
plt.figure(2)
plt.plot(T, energy, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Average energy')
plt.figure(3)
plt.plot(T, specificheat, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Specific Heat')
plt.figure(4)
plt.plot(T, susceptibility, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Susceptibility')
simulate(N, top)
观测结果:
- Magnetization trend as a function of temperature
- Specific Heat
我已经尝试对代码进行大量注释,以防万一我忽略了什么,请询问。
问题:
- 磁化趋势是否正确?磁化随温度升高而降低,但无法确定相变的临界温度。
- 随着温度升高,能量接近于零,这似乎与在伊辛网格中观察到的一致。为什么我得到负比热值?
- 如何选择 monte carlo 步数?这只是基于网络节点数的命中和试用吗?
编辑:02.06::反铁磁配置的模拟失效
首先,由于这是一个编程站点,我们来分析一下程序。您的计算效率非常低,这使得探索更大的图变得不切实际。您的情况下的邻接矩阵是 200x200 (40000) 个元素,只有大约 3% 的非零元素。将其转换为密集矩阵意味着在评估 montecarlo
例程中的能量差和 netenergy
中的净能量时需要更多的计算。以下代码在我的系统上的执行速度提高了 5 倍,并且预计更大的图形会带来更好的加速:
# keep the topology as a sparse matrix
top = nx.adjacency_matrix(G)
def netenergy(N, state):
en = 0.
for i in range(N):
ss = np.sum(state[top[i].nonzero()[1]])
en += state[i] * ss
return -0.5 * J * en
注意因子中的 0.5 - 因为邻接矩阵是对称的,所以每对自旋被计算两次!
def montecarlo(state, N, beta, top):
# pick a random source node
rsnode = np.random.randint(0, N)
# get the spin of this node
s = state[rsnode]
# sum of all neighbouring spins
ss = np.sum(state[top[rsnode].nonzero()[1]])
# transition energy
delE = 2.0 * J * ss * s
# calculate transition probability
prob = math.exp(-delE * beta)
# conditionally accept the transition
if prob > random.random():
s = -s
state[rsnode] = s
return state
注意跃迁能量中的因子 2.0 - 您的代码中缺少它!
这里有一些 numpy
索引魔法。 top[i]
是节点 i 的稀疏邻接行向量,top[i].nonzero()[1]
是非零元素的列索引(top[i].nonzero()[0]
是行索引,它们都等于 0,因为它是一个行向量)。 state[top[i].nonzero()[1]]
因此是节点 i.
的相邻节点的值
现在说物理。热力学性质是错误的,因为:
e1 = e0 + ene * ene
m1 = m0 + mag * mag
真的应该是:
e1 = e1 + ene * ene
m1 = m1 + mag * mag
并且:
specificheat.append( e1/steps - e0*e0/(steps*steps) /(N* t * t))
susceptibility.append( m1/steps - m0*m0/(steps*steps) /(N* t *t))
应该是:
specificheat.append((e1/steps/N - e0*e0/(steps*steps*N*N)) / (t * t))
susceptibility.append((m1/steps/N - m0*m0/(steps*steps*N*N)) / t)
(你最好早点把能量和磁化平均)
这使热容量和磁化率达到正值。请注意敏感性分母中的单个 t
。
现在程序(希望)是正确的,让我们谈谈物理学。对于每个温度,您都从一个完全向上的自旋状态开始,然后让它一次演化一个自旋。显然,除非温度为零,否则此初始状态远离热平衡,因此系统将开始向对应于给定温度的状态 space 的部分漂移。此过程称为热化,在此期间收集静态信息毫无意义。您必须始终将给定温度下的模拟分为两部分 - 热化和实际显着 运行。需要多少次迭代才能达到平衡?很难说 - 采用能量的移动平均值并在它变得(相对)稳定时进行监控。
其次,更新算法每次迭代更改一次自旋,这意味着程序将探索状态 space 非常缓慢,您将需要大量迭代才能获得良好的近似值分区函数。旋转 200 次,1000 次迭代可能就足够了。
其余问题真的不属于这里。
我正在尝试模拟 Barabasi-Albert 网络中的 Ising 相变,并尝试复制一些可观测值的结果,例如在 Ising 网格模拟中观察到的磁化强度和能量。但是,我在解释我的结果时遇到了麻烦:不确定物理学是否错误或实施中是否存在错误。这是一个最小的工作示例:
import numpy as np
import networkx as nx
import random
import math
## sim params
# coupling constant
J = 1.0 # ferromagnetic
# temperature range, in units of J/kT
t0 = 1.0
tn = 10.0
nt = 10.
T = np.linspace(t0, tn, nt)
# mc steps
steps = 1000
# generate BA network, 200 nodes with preferential attachment to 3rd node
G = nx.barabasi_albert_graph(200, 3)
# convert csr matrix to adjacency matrix, a_{ij}
adj_matrix = nx.adjacency_matrix(G)
top = adj_matrix.todense()
N = len(top)
# initialize spins in the network, ferromagnetic
def init(N):
return np.ones(N)
# calculate net magnetization
def netmag(state):
return np.sum(state)
# calculate net energy, E = \sum J *a_{ij} *s_i *s_j
def netenergy(N, state):
en = 0.
for i in range(N):
for j in range(N):
en += (-J)* top[i,j]*state[i]*state[j]
return en
# random sampling, metropolis local update
def montecarlo(state, N, beta, top):
# initialize difference in energy between E_{old} and E_{new}
delE = []
# pick a random source node
rsnode = np.random.randint(0,N)
# get the spin of this node
s2 = state[rsnode]
# calculate energy by summing up its interaction and append to delE
for tnode in range(N):
s1 = state[tnode]
delE.append(J * top[tnode, rsnode] *state[tnode]* state[rsnode])
# calculate probability of a flip
prob = math.exp(-np.sum(delE)*beta)
# if this probability is greater than rand[0,1] drawn from an uniform distribution, accept it
# else retain current state
if prob> random.random():
s2 *= -1
state[rsnode] = s2
return state
def simulate(N, top):
# initialize arrays for observables
magnetization = []
energy = []
specificheat = []
susceptibility = []
for count, t in enumerate(T):
# some temporary variables
e0 = m0 = e1 = m1 = 0.
print 't=', t
# initialize spin vector
state = init(N)
for i in range(steps):
montecarlo(state, N, 1/t, top)
mag = netmag(state)
ene = netenergy(N, state)
e0 = e0 + ene
m0 = m0 + mag
e1 = e0 + ene * ene
m1 = m0 + mag * mag
# calculate thermodynamic variables and append to initialized arrays
energy.append(e0/( steps * N))
magnetization.append( m0 / ( steps * N))
specificheat.append( e1/steps - e0*e0/(steps*steps) /(N* t * t))
susceptibility.append( m1/steps - m0*m0/(steps*steps) /(N* t *t))
print energy, magnetization, specificheat, susceptibility
plt.figure(1)
plt.plot(T, np.abs(magnetization), '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Average Magnetization per spin')
plt.figure(2)
plt.plot(T, energy, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Average energy')
plt.figure(3)
plt.plot(T, specificheat, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Specific Heat')
plt.figure(4)
plt.plot(T, susceptibility, '-ko' )
plt.xlabel('Temperature (kT)')
plt.ylabel('Susceptibility')
simulate(N, top)
观测结果:
- Magnetization trend as a function of temperature
- Specific Heat
我已经尝试对代码进行大量注释,以防万一我忽略了什么,请询问。
问题:
- 磁化趋势是否正确?磁化随温度升高而降低,但无法确定相变的临界温度。
- 随着温度升高,能量接近于零,这似乎与在伊辛网格中观察到的一致。为什么我得到负比热值?
- 如何选择 monte carlo 步数?这只是基于网络节点数的命中和试用吗?
编辑:02.06:
首先,由于这是一个编程站点,我们来分析一下程序。您的计算效率非常低,这使得探索更大的图变得不切实际。您的情况下的邻接矩阵是 200x200 (40000) 个元素,只有大约 3% 的非零元素。将其转换为密集矩阵意味着在评估 montecarlo
例程中的能量差和 netenergy
中的净能量时需要更多的计算。以下代码在我的系统上的执行速度提高了 5 倍,并且预计更大的图形会带来更好的加速:
# keep the topology as a sparse matrix
top = nx.adjacency_matrix(G)
def netenergy(N, state):
en = 0.
for i in range(N):
ss = np.sum(state[top[i].nonzero()[1]])
en += state[i] * ss
return -0.5 * J * en
注意因子中的 0.5 - 因为邻接矩阵是对称的,所以每对自旋被计算两次!
def montecarlo(state, N, beta, top):
# pick a random source node
rsnode = np.random.randint(0, N)
# get the spin of this node
s = state[rsnode]
# sum of all neighbouring spins
ss = np.sum(state[top[rsnode].nonzero()[1]])
# transition energy
delE = 2.0 * J * ss * s
# calculate transition probability
prob = math.exp(-delE * beta)
# conditionally accept the transition
if prob > random.random():
s = -s
state[rsnode] = s
return state
注意跃迁能量中的因子 2.0 - 您的代码中缺少它!
这里有一些 numpy
索引魔法。 top[i]
是节点 i 的稀疏邻接行向量,top[i].nonzero()[1]
是非零元素的列索引(top[i].nonzero()[0]
是行索引,它们都等于 0,因为它是一个行向量)。 state[top[i].nonzero()[1]]
因此是节点 i.
现在说物理。热力学性质是错误的,因为:
e1 = e0 + ene * ene
m1 = m0 + mag * mag
真的应该是:
e1 = e1 + ene * ene
m1 = m1 + mag * mag
并且:
specificheat.append( e1/steps - e0*e0/(steps*steps) /(N* t * t))
susceptibility.append( m1/steps - m0*m0/(steps*steps) /(N* t *t))
应该是:
specificheat.append((e1/steps/N - e0*e0/(steps*steps*N*N)) / (t * t))
susceptibility.append((m1/steps/N - m0*m0/(steps*steps*N*N)) / t)
(你最好早点把能量和磁化平均)
这使热容量和磁化率达到正值。请注意敏感性分母中的单个 t
。
现在程序(希望)是正确的,让我们谈谈物理学。对于每个温度,您都从一个完全向上的自旋状态开始,然后让它一次演化一个自旋。显然,除非温度为零,否则此初始状态远离热平衡,因此系统将开始向对应于给定温度的状态 space 的部分漂移。此过程称为热化,在此期间收集静态信息毫无意义。您必须始终将给定温度下的模拟分为两部分 - 热化和实际显着 运行。需要多少次迭代才能达到平衡?很难说 - 采用能量的移动平均值并在它变得(相对)稳定时进行监控。
其次,更新算法每次迭代更改一次自旋,这意味着程序将探索状态 space 非常缓慢,您将需要大量迭代才能获得良好的近似值分区函数。旋转 200 次,1000 次迭代可能就足够了。
其余问题真的不属于这里。