机器学习二维系统相变

Abstract

凝聚态物理是研究诸如电子、原子核、原子、磁偶极子、量子比特等组成的复杂系统。其复杂性体现在随粒子数指数增长的相空间，让人想到了机器学习中的纬度诅咒。尽管存在这样的诅咒，机器学习依然可以在数据集上体现出很好的分类、回归性质。这里，我们用现代的机器学习技术（全链接、卷积网络），用以分类众多不同的哈密顿体系相与相变。通过调用库，神经网络可以训练用于区分序参量。

这是2017年Nature phsics上的一篇文献，其证明了机器学习方法可以在凝聚态物理中处理相变问题。区别于高能物理，凝聚态物理是研究多体系统的物理，也是公认研究人数最多，范围最广的物理分支。相变理论又因其存在临界现象，又是凝聚态物理中最有趣的一个分支。二阶朗道相变理论给出了自发对称性，让人感到物理美的地方。文献里主要是在研究二维铁磁系统，也就是Ising model。对Ising model做了二分类，分为有磁性和无磁性。方法简单暴力，而且精度很高。文献里的模型，没啥特殊的物理意义（复现最大的意义就是让我复习一下统计物理+paddle......）。接下来介绍一些Ising model的基本理论，随后对其进行分类。（项目中不给出任何公式，大家是非科班生，不必执迷于物理意义的部分。这里仅仅给出定性的描述。要给出一个很好的物理图像，需要掌握很多么门物理专业的课程。作者自己都不敢说自己这些基础的物理图像做得多好。）

一 相变理论

• 系综理论
• Ising Model
• 蒙特卡洛模拟

系综理论

由概率理论为基础的统计力学，完全不同于经典力学方法。这种引入统计手段为描述力学性质的方法，为20世纪初量子力学的诞生埋下伏笔。正是这样的背景，统计力学和经典理论、量子理论不同，却又相同，量子物理结合电动力学发展出的量子电动力学，非常好得描述了单粒子系统的性质。量子物理结合统计力学，也就是我们的凝聚态物理。分析力学中大家寻找哈密顿量，又或者是构造拉格朗日函数，就得到了系统含时演化的所有信息。量子力学中大家解薛定谔方程，也是为了找到一个波函数，用于完全描述系统。统计物理里也一样，这次我们寻找的函数叫做配分函数。配分函数代表了每一个state在整个phase space里出现的概率。

举个栗子：抛硬币

一个硬币

0000000000000000000000000000000000000000000000000000000000000（手动分割线）

0000000000000000000000000000000000000000000000000000000000000（手动分割线）

假如你手里有一个硬币，显而易见，其state有俩种：正面朝上和反面向下。那么你咋能知道每一个state出现的概率呢？其一种思路是拿着这个硬币一直抛，最终得到俩个state出现的频率，最终趋于概率（大数定律）。还有一种思路，拿一箱子硬币，然后晃箱子。开箱后开始计数，一次性得到概率。我们假设俩种方法得到的结果一致，也就是含时平均等于系综平均，又称为各态遍历假说。这是统计理论的支柱，靠着这一条基本假设，解决了多粒子系统力学耦合无法描述的问题。

我们就假设，手上有一百个这样相同边界条件的系统。那么其中某一个特定state出现了几次，就是这个state在space里的概率密度。上文的硬币，明显就是1/2。

Codiga

可自定义的静态代码分析检测工具

下载

Ising Model

Ising Model是最简单又好用的统计模型。其存在量子物理中的自旋，可又是建立在经典系统理论上。那ta到底讲了个啥？来做一个小实验。拿起一个磁铁，然后用打火机烧一烧。随后发现，磁铁的磁性竟然消失了！！？这其中到底发生了什么？？

我们的Ising model正是在说这样的事情。我们说有磁性的磁铁，变成了没有磁性的废铁，这一过程发生了相变。引入物理量，磁化强度作为序参量，描述这一现象。当温度低于临界值，磁化强度不为0，当高于某一温度，磁化强度恒为0。我们称这一过程，发生了相变。Ising Model，把所有原子视为固定的。（固体活动范围小，忽略）。那么电子也是固定的。但是电子存在自旋，每一个电子都贡献一个磁场强度，最后就是一个大磁铁了。

这是数据集中的一张有序的图片。

这是无序的图片

黑色代表自旋向上，白色代表自旋向下，每一个像素点代表一个电子。 有序（序参量不为0）、无序给人的映像很直观，就是字面意思。有序代表秩序，无序则代表混沌。无序情况下熵值高，不确定性大。有序表现为聚集，自旋向下的聚成一块，自旋向上的聚为一块。对于论文复现，我们知道到这就够了。下面贴上Ising Model的蒙特卡洛模拟代码，用以生成这样的图片。值得一提的是，我们假设2.269度为相变临界点。高于2.269则无序（顺磁），低于则有序（铁磁）。

交代清楚任务，我们便要开始进行机器学习。首先的任务是进行蒙特卡洛仿真，模拟二维Ising model。运行这一段代码生成数据集，但时间过长，至少需要一个星期才能跑完。我把自己运行好的数据集上传到了data里，想要自己生成数据集的可以运行代码观察一下热容在临界点附近的突变。时间来不及可以直接略过。

#这一段代码是从网上直接copy的，原网址：https://rajeshrinet.github.io/blog/2014/ising-model/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 mag

nt      = 2**8        # 温度点数量N       = 4**2        # 点阵尺寸, N x N  eqSteps = 2**10       # MC方法平衡步数mcSteps = 2**4        # MC方法计算步数n1, n2  = 1.0/(mcSteps*N*N), 1.0/(mcSteps*mcSteps*N*N)
tm = 2.269;    T=np.random.normal(tm, .64, nt)          #设2.269为临界点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)*iT#计算四张图f = 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);

二 模型实现

• 数据集处理
• 模型组网
• 模型训练
• 模型验证

1.数据集处理

从非平衡态过度到平衡态需要走很多蒙特卡罗步，耗时极长。为了方便大家复现，我把自己生成的数据也打包上传到了项目里。

!unzip /home/aistudio/data/data110989/dataset.zip -d /home/aistudio/work/data

import osimport numpy as npfrom PIL import Imageimport matplotlib.pyplot as pltimport paddle# 生成图像列表################################################################拿出来保存数据的路径data_path = '/home/aistudio/work/data'character_folders = os.listdir(data_path)###############################################################################################################################当已经存在图像列表时，把图像列表删掉，以便于重新写一个列表。if(os.path.exists('./train_data.txt')):
    os.remove('./train_data.txt')if(os.path.exists('./test_data.txt')):
    os.remove('./test_data.txt')###############################################################################################################################character_folders是列表，是data_pathfor character_folder in character_folders:    with open('./train_data.txt', 'a') as f_train:        with open('./test_data.txt', 'a') as f_test:            if character_folder == '.DS_Store':                continue     #continue 跳出本次循环
            character_imgs = os.listdir(os.path.join(data_path,character_folder))
            count = 0 
            for img in character_imgs:                if img == '.DS_Store':                    continue
                if count%10 == 0:
                    f_test.write(os.path.join(data_path,character_folder,img) + '\t' + character_folder + '\n')                else:
                    f_train.write(os.path.join(data_path,character_folder,img) + '\t' + character_folder + '\n')
                count +=1print('列表已生成')###############################################################

 with open('train_data.txt',"r") as f:
            path_sample=f.readline()         #拿出来train_data里保存的第一个路径。
            sample=path_sample.split("\t")   #\t是空格，这里的文本默认有一个空格 以区分路径与标签 
            label=sample[1]  
            sample=sample[0]                 #刚才取出来的列表，第一个元素是路径。（python是从0开始）第二个是标签。我们拿前一个
            img=Image.open(sample)
            plt.imshow(img)
            img=np.array(img)            print(img)
            plt.show()            print(label)                     #打印标签

import numpy as npfrom paddle.io import Datasetclass MyDataset(Dataset):
    """
    数据集类的定义
    """
    def __init__(self, mode='train_data',transform=None):
        """
        初始化函数
        """
        super(MyDataset,self).__init__()
        self.data = []
        self.transform=transform        with open(f'{mode}.txt') as f:            for line in f.readlines():
                info = line.strip().split('\t')     #\t指缩进   strip是去掉字符串前后的一些特殊符号
                if len(info) > 0:
                    self.data.append([info[0].strip(), info[1].strip()])                       
    def __getitem__(self, index):
        """
        根据索引获取单个样本
        """
        image_file, label = self.data[index]
        img = Image.open(image_file)  
        img = np.array(img)                          
        if self.transform is not None:   
            img=self.transform(img)        if label=="有序":
            label=1
        else:
            label=0
        
        return img, label    def __len__(self):
        """
        获取样本总数
        """
        return len(self.data)

import paddle.nn.functional as Ffrom paddle.vision.transforms import Compose, Resize, ToTensor, Normalize 
#compose的功能是，把要做的预处理以列表组合起来。#这里用了图像归一化处理和调整图像大小，还可以调用其他的。transform = Compose([ToTensor("CHW")]) #CHW可以理解为BGR格式。           #ToTensor是将输入的numpy数据转为paddle的数据类型，暂时可以不管这个。

train_dataset = MyDataset(mode='train_data',transform=transform)


test_dataset = MyDataset(mode='test_data',transform=transform)

2.模型组网

class MyModel(paddle.nn.Layer):
    def __init__(self):
        super(MyModel, self).__init__()        #(1)加一个卷积层。四个参数分别是(输入图像通道，卷积后输出图像通道数，卷积核大小，卷积步长)
        #一个filter经过一次抓取，生成一个feature map.不同的filter抓取不同特征。这里有2个filter
        #输出的尺寸：输入数据的大小-卷积核大小+2×边界填充0的数量）/移动步长+1
        #[（256-5+2*0）/1]+1=252 
        
        self.conv1 = paddle.nn.Conv2D(in_channels=1,   
                                      out_channels=6,
                                      kernel_size=5,
                                      stride=1)        #(2)加一个池化层。 
        #输出的尺寸：（输入数据的大小-过滤器大小）/移动步长+1
        #[(252-4)/2]+1=125
        self.pool1 = paddle.nn.MaxPool2D(kernel_size=4,
                                         stride=2)        #（3）
        #这一层的功能是把张量拉平成一维向量。输入是多少，输出就是多少，改变的仅仅是形状。
        self.flatten=paddle.nn.Flatten(start_axis=1,stop_axis=-1)        #（4）加一个线性层。
        #线性层的输入是6*125*125  6是因为我们的conv1有6个out_channels 
        self.fc1=paddle.nn.Linear(6*125*125,100)        #再叠一个relu
        self.relu1=paddle.nn.ReLU()
        self.fc2=paddle.nn.Linear(100,50)
        self.fc3=paddle.nn.Linear(50,10)
        self.fc4=paddle.nn.Linear(10,2)    #正向传播过程
    def forward(self, x):     #调用之前的网络
        x = self.conv1(x)
        x = self.pool1(x)
        x = self.flatten(x)     
        x = self.fc1(x)    
        x = self.relu1(x)
        x = self.fc2(x)
        x = self.fc3(x)
        x = self.fc4(x)        return x

model=paddle.Model(MyModel())
model.prepare(paddle.optimizer.Adam(parameters=model.parameters()),    #指定优化方法，这里选择了Adam
              paddle.nn.CrossEntropyLoss(),                            #选择损失函数，这里是交叉熵
              paddle.metric.Accuracy())                                #评估，计算正确率

model.summary((1,1,256,256))

3.模型训练

model.fit(train_dataset,    #使用的数据
          epochs=10,         #训练的轮数
          batch_size=50,   #单次送进多少个数据         
          verbose=1)        #用进度条的形式打印日志

4.模型验证

model.evaluate(test_dataset,verbose=1)