氦泡在氚化钛中的建模

背景介绍

钛是一种良好的储氢元素,也可以用来储存氚(氢的同位素),在核物理领域应用十分广泛,但是氚会衰变生成氦元素,这些氦原子在氚化钛晶体中会逐渐聚集形成氦泡,并最终影响基体的力学性质，也就是材料的老化，为了研究其微观性质，我们建立了氚化钛中含氦泡的模型。

建模过程

基础晶体建模，这里在每一个钛晶胞中添加一个氚原子。

import numpy as np

def bulid_structure(Lattice_Constant, Structure_Type, x, y, z):
    """
    该函数用于创建常见金属晶体结构！
    :param Lattice_Constant: 晶格常数
    :param Structure_Type: 金属晶体结构：FCC, BCC, HCP
    :param x: x方向拓展个数
    :param y: y方向拓展个数
    :param z: z方向拓展个数
    :return: 原子坐标 N-3 array
    """

    if Structure_Type == 'FCC':
        # Cubic FCC basis
        basis = np.array([[1.0, 0.0, 0.0],
                        [0.0, 1.0, 0.0],
                        [0.0, 0.0, 1.0]]) * Lattice_Constant

        base_atoms = np.array([[0.0, 0.0, 0.0],
                            [0.5, 0.5, 0.0],
                            [0.5, 0.0, 0.5],
                            [0.0, 0.5, 0.5]]) * Lattice_Constant

    elif Structure_Type == 'HTi':
        basis = np.array([[1.0, 0.0, 0.0],
                        [0.0, np.sqrt(3), 0.0],
                        [0.0, 0.0, np.sqrt(8 / 3)]]) * Lattice_Constant

        base_atoms = np.array([[0.0, 0.0, 0.0],
                            [0.5, 0.5 * np.sqrt(3), 0.0],
                            [0.5, np.sqrt(3) * 5 / 6, 0.5 * np.sqrt(8 / 3)],
                            [0.0, 1 / 3 * np.sqrt(3), 0.5 * np.sqrt(8 / 3)],
                            [0.25, np.sqrt(3) * 5 / 12, 1 / 4 * np.sqrt(8 / 3)]]) * Lattice_Constant
        
        atype = [1, 1, 1, 1, 2]
        positions = []
        atypelist = []
        for i in range(x):
            for j in range(y):
                for k in range(z):
                    base_position = np.array([i, j, k], dtype=float)
                    cart_position = np.dot(basis, base_position)
                    atypelist += atype
                    for atom in base_atoms:
                        positions.append(cart_position + atom)

        positions = np.array(positions)
        atypelist = np.array(atypelist)
        return positions, atypelist

    positions = []
    for i in range(x):
        for j in range(y):
            for k in range(z):
                base_position = np.array([i, j, k])
                cart_position = np.inner(basis.T, base_position)
                for atom in base_atoms:
                    positions.append(cart_position + atom)
    return np.array(positions)

加入单个氦泡

def onevoidHeTiH(R=250, Lbox=520, L_HTi=2.9511, L_He=2.98):
    x, y, z = np.round(Lbox / np.array([1, np.sqrt(3), np.sqrt(8 / 3)]) / L_HTi).astype(int)
    HTi, HTi_Type = bulid_structure(L_HTi, 'HTi', x, y, z)
    HTi -= np.mean(HTi, axis=0)
    index = np.where(HTi_Type==2)[0]
    N = len(index)
    np.random.seed(1000)
    index1 = np.random.choice(index, size = int(0.95*N), replace=False)
    select = np.ones(len(HTi), dtype=bool)
    select[index1] = False
    HTi = HTi[select]
    HTi_Type = HTi_Type[select]
    HeM = bulid_structure(L_He, 'FCC', 200, 200, 200)
    HeM -= np.mean(HeM, axis=0)
    HeM = HeM[np.sum(np.square(HeM), axis=1) < R**2]
    void = np.sum(np.square(HTi), axis=1) > R**2
    HTi = HTi[void]
    HTi_Type = HTi_Type[void]
    positions = np.r_[HTi, HeM]
    atypelist = HTi_Type.tolist() + [3]*len(HeM)
    ratio = np.round(len(HeM) / len(HTi_Type[HTi_Type==1]), 2)

    with open(f'H-Ti-{Lbox}-{L_He}-{R}-{ratio}.data', 'w') as op:

        # First line is a comment line
        op.write('HTi structure, written by HerrWu\n\n')
        op.write('{} atoms\n'.format(len(positions)))
        op.write('3 atom types\n')

        op.write('{} {} xlo xhi\n'.format(HTi[:,0].min(), HTi[:,0].min() + x * L_HTi * 1.0))
        op.write('{} {} ylo yhi\n'.format(HTi[:,1].min(), HTi[:,1].min() + y * L_HTi * np.sqrt(3)))
        op.write('{} {} zlo zhi\n\n'.format(HTi[:,2].min(), HTi[:,2].min() + z * L_HTi * np.sqrt(8 / 3)))

        # Atoms section
        op.write('Atoms\n\n')

        # Write each Position
        for idx, i in enumerate(positions):
            op.write(f'{idx+1} {atypelist[idx]} {i[0]} {i[1]} {i[2]}\n')
onevoidHeTiH(R=250, Lbox=550, L_HTi=2.9511, L_He=5.75)

结果如下：

加入多个氦泡

def numvoidHeTiH(nvoid, R=10, Lbox=100, L_HTi=2.9511, L_He=2.98):
    x, y, z = np.round(Lbox / np.array([1, np.sqrt(3), np.sqrt(8 / 3)]) / L_HTi).astype(int)
    HTi, HTi_Type = bulid_structure(L_HTi, 'HTi', x, y, z)
    HTi -= np.mean(HTi, axis=0)
    index = np.where(HTi_Type==2)[0]
    N = len(index)
    index1 = np.random.choice(index, size = int(0.95*N), replace=False)
    select = np.ones(len(HTi), dtype=bool)
    select[index1] = False
    HTi = HTi[select]
    HTi_Type = HTi_Type[select]
    
    HeM = bulid_structure(L_He, 'FCC', 80, 80, 80)
    HeM -= np.mean(HeM, axis=0)
    void = np.ones(len(HTi), dtype=bool)
    HeMpos = []
    print(np.random.randint(-30, 30, size=(nvoid, 3)))
    for coma in np.random.randint(-30, 30, size=(nvoid, 3)):
        void = np.logical_and(void, (np.sum(np.square(HTi - coma), axis=1) > R**2))
        hp = HeM.copy()
        hp = hp[np.sum(np.square(hp), axis=1) < R**2]
        hp += coma
        HeMpos.append(hp)
    HeMpos = np.vstack(HeMpos)

    HTi = HTi[void]
    HTi_Type = HTi_Type[void]

    positions = np.r_[HTi, HeMpos]
    atypelist = HTi_Type.tolist() + [3]*len(HeMpos)
    ratio = np.round(len(HeMpos) / len(HTi_Type[HTi_Type==1]), 2)

    with open(f'H-Ti-{Lbox}-{L_He}-{R}-{ratio}.data', 'w') as op:

        # First line is a comment line
        op.write('HTi structure, written by HerrWu\n\n')
        op.write('{} atoms\n'.format(len(positions)))
        op.write('3 atom types\n')

        op.write('{} {} xlo xhi\n'.format(HTi[:,0].min(), HTi[:,0].min() + x * L_HTi * 1.0))
        op.write('{} {} ylo yhi\n'.format(HTi[:,1].min(), HTi[:,1].min() + y * L_HTi * np.sqrt(3)))
        op.write('{} {} zlo zhi\n\n'.format(HTi[:,2].min(), HTi[:,2].min() + z * L_HTi * np.sqrt(8 / 3)))

        # Atoms section
        op.write('Atoms\n\n')

        # Write each Position
        for idx, i in enumerate(positions):
            op.write(f'{idx+1} {atypelist[idx]} {i[0]} {i[1]} {i[2]}\n')

numvoidHeTiH(nvoid=6)