Python和CUDA(C++)量子退火和伊辛二次算法模型

🎯要点

1. 简化量子退火或离散优化算法处理，使用张量网络模拟和动态系统方法及神经网络逼近。
2. 实现并行退火算法和CUDA支持下穷举搜索法。
3. 使用大都会算法模拟二维自旋玻璃伊辛模型并测量磁化率、比热容和能量。
4. 对比其他组合优化解方法，使用英伟达A100 GPU测试。

🍁自旋玻璃退火算法

🍪语言内容分比

🍇Python伊辛模型

伊辛模型（或伦茨-伊辛模型）以物理学家恩斯特·伊辛和威廉·伦茨的名字命名，是统计力学中铁磁性的数学模型。该模型由离散变量组成，这些变量表示原子“自旋”的磁偶极矩，可以处于两种状态之一（1 或 −1）。自旋排列成一个图形，通常是晶格（其中局部结构在所有方向上周期性重复），允许每个自旋与其邻居相互作用。相邻的一致自旋的能量低于不一致的自旋；系统趋向于最低能量，但热量会干扰这种趋势，从而产生不同结构相的可能性。该模型允许将相变识别为现实的简化模型。二维方晶格伊辛模型是显示相变的最简单的统计模型之一。

考虑一组 $\Lambda$ 个晶格点，每个点都有一组相邻点（例如一个图）形成一个 $d$ 维晶格。对于每个晶格点 $k\in \Lambda$，都有一个离散变量 ${\sigma }_k$，使得 σ k ∈ { − 1 , + 1 } \sigma_k \in\{-1,+1\} σk​∈{−1,+1}，表示点的自旋。自旋配置 σ = { σ k } k ∈ Λ \sigma=\left\{\sigma_k\right\}_{k \in \Lambda} σ={σk​}k∈Λ​ 是每个晶格点的自旋值分配。对于任何两个相邻位置 $i,j\in \Lambda$，存在相互作用 $J_{ij}$。此外，位置 $j\in \Lambda$ 有一个外部磁场 $h_j$ 与其相互作用。配置 $\sigma$ 的能量由哈密顿函数给出：
$$H(\sigma )=-\sum _{⟨ij⟩}{J}_{ij}{\sigma }_i{\sigma }_j-\mu \sum _j{h}_j{\sigma }_j$$
其中第一个和是对相邻自旋对的和（每对都计算一次）。符号 $⟨ij⟩$ 表示位置 $i$ 和 $j$ 是最近的邻居。磁矩由 $\mu$ 给出。请注意，上面哈密顿量的第二项中的符号实际上应该是正的，因为电子的磁矩与其自旋反向平行，但通常使用负项。配置概率由温度倒数为 $\beta \ge 0$ 的玻尔兹曼分布给出：
$$P_{\beta }(\sigma )=\frac{e^{-\beta H(\sigma )}}{Z_{\beta }}$$
其中 $\beta =1/\left(k_{B}T\right)$，归一化常数
$$Z_{\beta }=\sum _{\sigma }{e}^{-\beta H(\sigma )}$$
是配分函数。对于自旋函数 $f$（“可观测”），我们表示为
$$⟨f{⟩}_{\beta }=\sum _{\sigma }f(\sigma )P_{\beta }(\sigma )$$
$f$ 的期望（均值）。

配置概率 $P_{\beta }(\sigma )$ 表示系统（在平衡状态下）处于配置 $\sigma$ 状态的概率。

二维模型模拟

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

def initialstate(N):   
    state = 2*np.random.randint(2, size=(N,N))-1
    return state


def mcmove(config, beta):
    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 config


def calcEnergy(config):
    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):
    mag = np.sum(config)
    return mag

更改这些参数以实现更小（更快）的模拟

nt      = 88         
N       = 16         
eqSteps = 1024       
mcSteps = 1024       

T       = np.linspace(1.53, 3.28, nt); 
E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
n1, n2  = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N) 

for tt in range(nt):
    E1 = M1 = E2 = M2 = 0
    config = initialstate(N)
    iT=1.0/T[tt]; iT2=iT*iT;
    
    for i in range(eqSteps):         
        mcmove(config, iT)          

    for i in range(mcSteps):
        mcmove(config, iT)           
        Ene = calcEnergy(config)     
        Mag = calcMag(config)        

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

    E[tt] = n1*E1
    M[tt] = n1*M1
    C[tt] = (n1*E2 - n2*E1*E1)*iT2
    X[tt] = (n1*M2 - n2*M1*M1)*iT

f = plt.figure(figsize=(18, 10));   

sp =  f.add_subplot(2, 2, 1 );
plt.scatter(T, E, s=50, marker='o', color='IndRed')
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Energy ", fontsize=20);         plt.axis('tight');

sp =  f.add_subplot(2, 2, 2 );
plt.scatter(T, abs(M), s=50, marker='o', color='RoyalBlue')
plt.xlabel("Temperature (T)", fontsize=20); 
plt.ylabel("Magnetization ", fontsize=20);   plt.axis('tight');

sp =  f.add_subplot(2, 2, 3 );
plt.scatter(T, C, s=50, marker='o', color='IndianRed')
plt.xlabel("Temperature (T)", fontsize=20);  
plt.ylabel("Specific Heat ", fontsize=20);   plt.axis('tight');   

sp =  f.add_subplot(2, 2, 4 );
plt.scatter(T, X, s=50, marker='o', color='RoyalBlue')
plt.xlabel("Temperature (T)", fontsize=20); 
plt.ylabel("Susceptibility", fontsize=20);   plt.axis('tight');

👉更新：亚图跨际