微观相结构的自动识别

微观相结构的自动识别

上次我们提到利用浓度数据的极值点进行相结构的识别,但是这样会遇到几个问题:首先,对于较为复杂的相极值点的数目相当多,同时不同嵌段比例以及支链的不同都会使各个相区的形状和大小发生变化,从而使极值点的位置和数目发生变化,给识别带来极大的困难;其次,如果浓度数据并不是按我们初始化的浓度进行排列的,而是进行了一定的平移或旋转变换,即相位的变化,这样极值点的位置也将发生变化,并且平移或旋转矩阵很难确定,因此又让这种识别方法变得几乎不可能。

考虑到微观相结构的周期性,利用傅里叶变换可以剔除掉相位的变化,从而将不同空间相位的浓度数据简化为傅里叶空间内单一的周期数据。实际上,对浓度的傅里叶变换等价于散射光的振幅。

因此类似于实验中的X射线衍射的方法可以利用傅里叶变换得到不同的衍射花纹,来确定微观相结构。但在傅里叶变换过程中,三维的傅里叶变化计算量很大,为了加快计算,因此我考虑对三维的晶胞进行三个方向的切片处理,得到三个二维的浓度矩阵,再做二维的离散傅里叶变换得到衍射图案。

在衍射图案中,在(0,0)主峰的周围,对于不同的相结构会出现不同的衍射图案,并且衍射图案是与微观相结构的对称性有关的,因此利用主峰周围其他峰的位置可以很清晰地判断微观相结构。以下为几个微观相结构A组分的衍射图案,为了简便只截取衍射图案的1/4:

σ相

BCC

Gyroid

FCC

A15

csHelix

在这里为了识别的简便化,因此降低傅里叶空间矩阵的大小,即在傅里叶变换时选取较小的u和v,来使花纹粗粒化,以Gyroid相为例,粗粒化之后三个面上A组分的衍射花纹分别为

除去原点位置的衍射峰,三个面上都有很明显的特征峰,可以根据这三个特征峰,即最大最小值的位置来作为这个相的特征进行判断,同时还有B和C两个组分的数据可以作为备用判断,这样就可以很方便的完成微观相行为的特征识别。最后附上一个python的简单的识别代码。

import numpy as np
import matplotlib.pyplot as plt

def setText(text):
    print(text)

#Defining size of pha
nx = 128
ny = 128
nz = 64

#Read pha file
phafile = open('pha.dat', 'r')
phalines = phafile.readlines()
pha = [[], [], []]
sf = [[], [], []]
for item in phalines:
    if not item == '\n':
        list = item.split(' ')
        pha[0].append(float(list[0]))
        pha[1].append(float(list[1]))
        pha[2].append(float(list[2]))
phafile.close()

pha = [np.array(item).reshape(nx, ny, nz) for item in pha]

#The size of Fourier Space divided by size of pha
d=1

#DFT
sf[0] = [np.fft.fft2(item[int(nx/2),::,::],[64*d,64*d]).real[0:int(4*d),0:int(4*d)] for item in pha]
sf[1] = [np.fft.fft2(item[::,int(ny/2),::],[64*d,64*d]).real[0:int(4*d),0:int(4*d)] for item in pha]
sf[2] = [np.fft.fft2(item[::,::,int(nz/2)],[64*d,64*d]).real[0:int(4*d),0:int(4*d)] for item in pha]
#Coordinate of sf: three slices, three components, 2D Fourier Square

#Eliminate the main peak at (0,0)
for item in sf:
    for subitem in item:
        subitem[0][0]=0

#Judging
if abs(sf[0][0][0][1] - np.min(sf[0][0])) < 0.1 and abs(sf[0][0][1][0] - np.min(sf[0][0])) < 0.1 and abs(sf[0][0][1][1] - np.max(sf[0][0])) < 0.1:
    setText('BCC')
elif abs(sf[0][0][0][2] - np.max(sf[0][0])) < 0.1 and abs(sf[0][0][2][0] - np.max(sf[0][0])) < 0.1 and abs(sf[0][0][1][1] - np.min(sf[0][0])) < 0.1:
    setText('FCC')
elif abs(sf[0][0][1][1] - np.min(sf[0][0])) < 0.1 and abs(sf[1][0][0][1] - np.max(sf[1][0])) < 0.1 and abs(sf[2][0][0][1] - np.max(sf[2][0])) < 0.1:
    setText('G')
elif abs(sf[0][0][2][0] - np.max(sf[0][0])) < 0.1 and abs(sf[0][0][1][2] - np.min(sf[0][0])) < 0.1 and abs(sf[1][0][0][2] - np.max(sf[1][0])) < 0.1:
    setText('A15')
elif abs(sf[0][0][0][1] - np.min(sf[0][0])) < 0.1 and abs(sf[0][0][1][1] - np.max(sf[0][0])) < 0.1 and abs(sf[1][0][2][0] - np.min(sf[1][0])) < 0.1:
    setText('csHelix')
elif abs(sf[0][0][1][0] - np.max(sf[0][0])) < 0.1 and abs(sf[0][0][2][1] - np.min(sf[0][0])) < 0.1 and abs(sf[1][0][2][1] - np.min(sf[1][0])) < 0.1:
    setText('csHelix2')
elif abs(sf[0][0][0][2] - np.max(sf[0][0])) < 0.1 and abs(sf[0][0][2][1] - np.min(sf[0][0])) < 0.1 and abs(sf[1][0][0][2] - np.max(sf[1][0])) < 0.1:
    setText('csσ')

#Plot
plt.imshow(sf[1][0])
plt.show()

下一步希望能够找到分析衍射图案和相结构对称性的关系,以从理论上解释这种判断方法的合理性。

 

One thought on “微观相结构的自动识别

发表评论