用于脉冲激光沉积的等离子体求解器

发表于 更新于 分类于

模型描述

计算初始速度分布

分布函数

该模型有三种分布函数类型：normal、log-normal以及Maxwell-Boltzmann。

求解晶胞中原子获取的动能

\[\begin{equation} E_{laser}=A\left(N_{uc}\left(E_{b}+\sum_{i=1}^{m}\left(E_{k}\left(i\right)+\Delta E_{exc}\left(i\right)\right)\right)+Q\right) \end{equation}\]

单组分二元氧化物靶材计算示例: TiO2

所有类和函数均放在本文最后。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
import sys

# 指定您想要添加为模块导入路径的目录
pathClasses = 'path/to/classes'

# 添加指定路径到sys.path
if pathClasses not in sys.path:
    sys.path.append(pathClasses)

# 现在您可以导入该路径下的模块了

from UnitCells import UnitCells
# import PhysicalConstants as C
# import PeriodicTable as PT

UC = UnitCells()

uc = UC.TiO2

# 激光参量

# laserFluence：激光能量密度
laserFluence    = 2
# spotWidth：光斑宽，m
spotWidth       = 1E-3
# spotHeight：光斑长, m
spotHeight      = 2.1E-3

# energyLaser： 激光能量
energyLaser = (laserFluence * 1E4) * (spotWidth * spotHeight)

# 首先是最简单的情况，TiO2靶材

# 材料参量

# absorbedRatio： 靶材对248 nm等离子体的吸收率。对于TiO2而言，为0.66
absorbedRatio = 0.66

# heatTarget： 靶材被激光加热的温度，对于陶瓷靶而言非常小
heatTarget = 0

# 每脉冲烧蚀深度，m
ablationDepth   = 100E-9

# Ablation volume [m^3]
ablationVolume  = spotWidth * spotHeight * ablationDepth

# 靶材中晶胞占比，为1代表只有一种晶胞
ucRatio = 1

# 单晶靶材，densityRatio为1； 
densityRatio = 1

# nUCAblated(iUC)，iuc靶材中的晶胞数，
# nUCAblated靶材中每个晶胞被烧蚀的原子数
# 对于单组份靶材，烧蚀晶胞数 = （烧蚀体积/晶胞体积）
nUCAblated = (ablationVolume / uc.VOLUME)

# energyExcitation：
# Excitation energy per atom in unit cell [J]
#   * Default: 0
energyExcitation = 0

energyKinetic = ((absorbedRatio * energyLaser - heatTarget)/(nUCAblated*sum(uc.AMOUNT))) - energyExcitation - uc.ENERGY_FORMATION

print(energyKinetic)

得到的结果与Matlab代码一致。

多组分靶材计算示例: SrO + TiO2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
import sys
import numpy as np

# 指定您想要添加为模块导入路径的目录
pathClasses = 'path/to/classes'

# 添加指定路径到sys.path
if pathClasses not in sys.path:
    sys.path.append(pathClasses)

# 现在您可以导入该路径下的模块了

from UnitCells import UnitCells
# import PhysicalConstants as C
# import PeriodicTable as PT

UC = UnitCells()

#-------------------------------------------------------------------------------
# Material parameters
#-------------------------------------------------------------------------------

# Unit cell
#   * Preface with UC.
#   * Any mixture of crystals is possible, eg. [UC.TiO2 UC.Li2O]
#   * To add new materials or to view material propertie, see UnitCells.m
uc = [UC.TiO2, UC.SrO]

# Mixture ratio of species in target
#   * Ignore if target contains only one component
#   * Should be of same length as unit cell array
ucRatio = np.array([1, 1])

# 激光参量

# laserFluence：激光能量密度
laserFluence    = 2
# spotWidth：光斑宽，m
spotWidth       = 1E-3
# spotHeight：光斑长, m
spotHeight      = 2.1E-3

# energyLaser： 激光能量
energyLaser = (laserFluence * 1E4) * (spotWidth * spotHeight)

# absorbedRatio： 靶材对248 nm等离子体的吸收率。对于TiO2而言，为0.66
absorbedRatio = 0.66

# heatTarget： 靶材被激光加热的温度，对于陶瓷靶而言非常小。默认为0
heatTarget = 0

# 每脉冲烧蚀深度，m
ablationDepth   = 100E-9

# Ablation volume [m^3]
ablationVolume  = spotWidth * spotHeight * ablationDepth

# 单晶靶材，densityRatio为1； 
densityRatio = 1

# energyExcitation：
# Excitation energy per atom in unit cell [J]
#   * Default: 0
energyExcitation = 0

#-------------------------------------------------------------------------------
# Material calculations
#-------------------------------------------------------------------------------

if len(uc) == 1:  # If target consists of only one unit cell
    # Set mixture ratio to 1
    ucRatio = 1 
    # Single crystal density
    scDensity = uc.DENSITY
    # Unit cell volume
    ucVolume = uc.VOLUME
else:
    # Normalize unit cell ratio 
    ucRatio = ucRatio / sum(ucRatio)
    # Weighted average single crystal density
    scDensity = sum(uc[i].DENSITY * ucRatio[i] for i in range(len(uc)))
    # Weighted average unit cell volume
    ucVolume = sum(uc[i].VOLUME * ucRatio[i] for i in range(len(uc)))

if len(uc) != len(ucRatio):
    raise ValueError("Number of elements of mixture ratio array must be the same as the unit cell array.")

# nUCAblated(iUC)，iuc靶材中的晶胞数，
# nUCAblated靶材中每个晶胞被烧蚀的原子数
nUCAblated = (ablationVolume / ucVolume)*ucRatio*densityRatio

# 以下内容后续应当写在initialVelocityDistribution.py中

for iUC in range(len(uc)):  # Python中的循环从0开始
    # 计算每个单元格的动能
    energyKinetic = ((absorbedRatio * energyLaser - heatTarget)
                    / (nUCAblated[iUC] * sum(uc[iUC].AMOUNT))
                    ) - energyExcitation - uc[iUC].ENERGY_FORMATION
    print("Energy Kinetic for UC", iUC+1, ":", energyKinetic)

得到的结果与Matlab代码一致.

原子的平均动能

\[\begin{equation} E_k\left(i\right)=\frac12m_i\nu_i^2 \end{equation}\]

1
2
3
4
5
6
7
# Loop through the elements in the unit cell
for iElement in range(len(uc[iUC].AMOUNT)):
    # Find index of element
    elementIndex = np.where(uc[iUC].ELEMENTS[iElement].NUMBER == np.array([atomX.NUMBER for atomX in atomUC]))[0]
    
    # Calculate the average velocity of each element in the unit cell
    veloRMS = np.sqrt(2 * energyKinetic / uc[iUC].ELEMENTS[iElement].MASS)

得到的结果与Matlab代码一致.

麦克斯韦-玻尔兹曼分布

\[\begin{equation} f_x(v)=\frac1{A_v}v^2\exp\left(-\frac{m_x v^2}{2k_BT}\right) \end{equation}\] 其中，\(A_v\) 是归一化常数，\(m_x\) 是原子质量，\(k_B\) 是玻尔兹曼常数，\(T\) 是热力学温度。

温度可以与每个粒子的平均动能相关联： \[\begin{equation} \bar{E}_{kin}=\frac{3}{2}k_{B}T\rightarrow k_{B}T=\frac{2}{3}\bar{E}_{kin} \end{equation}\] 其中，\(\bar{E}_{kin}\) 是平均动能。

因此最终的表达式： \[\begin{equation} f_x(v)=\frac1{A_v}v^2\exp\left(-\frac{3m_x v^2}{4\bar{E}_{kin,x}}\right) \end{equation}\]

具体代码见下面函数章节的

函数

getUniqueElements

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
import numpy as np
import sys

# 指定您想要添加为模块导入路径的目录
pathClasses = '/Users/zhouguangdi/Documents/data_for_software/Python/等离子体/PLD/classes'

# 添加指定路径到sys.path
if pathClasses not in sys.path:
    sys.path.append(pathClasses)


def getUniqueElements(uc, ucRatio):
    # GETUNIQUEELEMENTS Get unique elements present in target and their
    # normalized amount
    """
    获取目标中存在的独特元素及其规范化数量。
    """
    # 获取unit cell中的元素
    atomUCTemp = [ucx.ELEMENTS for ucx in uc]
    atomUCTemp = np.concatenate(atomUCTemp)

    # 根据unit cell之间的比率获取原子数量
    nAtomUCTemp = [amount * ratio for amount, ratio in zip([ucx.AMOUNT for ucx in uc], ucRatio)]
    nAtomUCTemp = np.concatenate(nAtomUCTemp)
    # 规范化原子数量数组
    total_uc_ratio = sum(ucRatio)
    nAtomUCTemp = [x / total_uc_ratio for x in nAtomUCTemp]
    min_amount = min(nAtomUCTemp)
    nAtomUCTemp = [x / min_amount for x in nAtomUCTemp]

    
    # 获取排序后的独特原子数组
    element_numbers = np.array([atomX.NUMBER for atomX in atomUCTemp])
    element_indices = np.argsort(element_numbers)
    unique_element_numbers, element_indices = np.unique(element_numbers, return_index=True)
    sorted_element_indices = element_indices #element_indices[np.argsort(element_indices)]

    # 检查是否有氧 (Z = 8)
    oxygen_indices = np.where(unique_element_numbers == 8)[0]
    if oxygen_indices.size > 0:
        # 如果氧不在首位，我们才将其移动到起始位置
        if oxygen_indices[0] != 0:
            sorted_element_indices = np.concatenate(([sorted_element_indices[oxygen_indices[0]]], np.delete(sorted_element_indices, oxygen_indices[0])))

    # 将独特原子插入数组
    atom_uc = [atomUCTemp[i] for i in element_indices]
    
    # 初始化原子数量数组
    n_atom_uc = np.zeros(len(unique_element_numbers))
    
    # Look through unique atoms
    for i, atom in enumerate(atom_uc):
        # 将每个原子的计算比例插入原子数量数组
        n_atom_uc[i] = sum(n for at, n in zip(atomUCTemp, nAtomUCTemp) if at.NUMBER == atom.NUMBER)
    

    #print(atom_uc[0].SYMBOL,atom_uc[1].SYMBOL,atom_uc[2].SYMBOL, n_atom_uc)

    return atom_uc, n_atom_uc

angularDistribution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import numpy as np

def angularDistribution(angle, angleDelta, radiusDelta, cosPowerFit, nAtomAblated):
    """
    Calculate the angular plasma particle distribution
    """
    
    # Number of angles
    nAngle = len(angle)
    
    # Pre-allocate memory
    nParticleAngle = np.zeros(nAngle)
    
    # Compute initial angular plasma particle distribution (eq. 5)
    for iAngle in range(nAngle - 1):
        nParticleAngle[iAngle] = ((4/3)*np.pi * radiusDelta**3) * \
            abs(-np.cos(np.radians(angle[iAngle + 1])) + np.cos(np.radians(angle[iAngle]))) * \
            np.cos(np.radians(angle[iAngle]))**cosPowerFit
    
    # Normalize and multiply by the total number of ablated atoms
    nParticleAngle = (nParticleAngle / np.sum(nParticleAngle)) * nAtomAblated
    # print(nParticleAngle)
    return nParticleAngle

initialVelocityDistribution

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
# %INITIALVELOCITYDISTRIBUTION Generate initial particle velocity ditribution
# %   uc                  : Array of type UnitCells, containing unit cells defined
# %                           in the class.
# %   ucRatio             : Array of numbers same length as uc, containing the
# %                           ratio between the unit cells in the target.
# %   atomUC              : Array containing the unique elements present in the target
# %   nUCAblated          : Array containing the number of unit cells ablated per
# %                           unit cell present in the target.
# %   velo                : Velocity array
# %   initVeloDistWidth   : Velocity distribution width, only valid for
# %                           normal and log-normal distributions.
# %   absorbedRatio       : Ratio of laser energy absorbed by the target.
# %   energyLaser         : Laser energy [J]
# %   heatTarget          : Heat dissipation into the target [J]
# %   energyExcitation    : Excitation energy per atom in the target [J]
# %   plotBool            : Boolean that determined is results get plotted or not.
# %   nMin                : Number of plasma particle threshold.
# %   nParticleTotal      : Total number of ablated particles.
# %   distributionType    : Initial velocity distribution type:
# %                             'normal', 'log-normal', 'Maxwell-Boltzmann'

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

def initialVelocityDistribution(uc, ucRatio, atomUC, nUCAblated, velo, initVeloDistWidth,
                                 absorbedRatio, energyLaser, heatTarget, energyExcitation,
                                 plotBool,nmin, nParticleTotal, distributionType):

    # 检查分布类型是否无效，如果是，则抛出错误
    valid_distributions = ['normal', 'log-normal', 'Maxwell-Boltzmann']
    if distributionType not in valid_distributions:
        raise ValueError(f'Invalid velocity distribution, choose: {", ".join(valid_distributions)}')

    # 初始化粒子速度分布数组
    nParticleVeloDist = np.zeros((len(atomUC), len(velo)))


    for iUC in range(len(uc)):  # Python中的循环从0开始
        # 计算每个单元格的动能
        energyKinetic = ((absorbedRatio * energyLaser - heatTarget)
                        / (nUCAblated[iUC] * sum(uc[iUC].AMOUNT))
                        ) - energyExcitation - uc[iUC].ENERGY_FORMATION
        #print("Energy Kinetic for UC", iUC+1, ":", energyKinetic)

        # Loop through the elements in the unit cell
        for iElement in range(len(uc[iUC].AMOUNT)):
            # Find index of element
            elementIndex = np.where(uc[iUC].ELEMENTS[iElement].NUMBER == np.array([atomX.NUMBER for atomX in atomUC]))[0][0]
            
            # Calculate the average velocity of each element in the unit cell
            veloRMS = np.sqrt(2 * energyKinetic / uc[iUC].ELEMENTS[iElement].MASS)
            
            # 正态分布
            if distributionType == 'normal':
                veloDistTemp = np.exp(- (velo - veloRMS)**2 / (2 * initVeloDistWidth**2)) \
                                / (np.sqrt(2*np.pi) * initVeloDistWidth)

            # 对数正态分布
            elif distributionType == 'log-normal':
                # 对数正态分布的均值
                mu = np.log(veloRMS**2 / np.sqrt(veloRMS**2 + initVeloDistWidth**2))
                
                # 对数正态分布的方差
                sigma = np.sqrt(np.log(1 + (initVeloDistWidth**2 / veloRMS**2)))
                
                veloDistTemp = np.concatenate(([0], np.exp(- (np.log(velo[1:]) - mu)**2 / (2 * sigma**2)) \
                                                        / (velo[1:] * (np.sqrt(2*np.pi) * sigma))))

            # Maxwell-Boltzmann 分布
            elif distributionType == 'Maxwell-Boltzmann':
                veloDistTemp = velo**2 * np.exp(-(3*uc[iUC].ELEMENTS[iElement].MASS * (velo**2)) \
                                                / (4 * energyKinetic))

            # 根据被消融粒子数量归一化
            veloDistTemp = (veloDistTemp / np.sum(veloDistTemp)) * (nUCAblated[iUC] * uc[iUC].AMOUNT[iElement])
            # 为指定元素添加分布
            nParticleVeloDist[elementIndex, :] += veloDistTemp
            

    nParticleVeloDist = (nParticleVeloDist / np.sum(nParticleVeloDist)) * nParticleTotal
    
    if plotBool:
        # 字符串表示物质公式
        formulaString = ""
        for i, u in enumerate(uc):
            if i == len(uc) - 1:
                formulaString += u.FORMULA
            else:
                formulaString += u.FORMULA + ' - '

        # 单元格比例字符串
        if len(ucRatio) > 1:
            ratioString = ""
            uc_ratio_temp = ucRatio / np.min(ucRatio)
            for i, ratio in enumerate(uc_ratio_temp):
                if i == len(uc) - 1:
                    ratioString += str(ratio)
                else:
                    ratioString += str(ratio) + ' : '

        # 初始化图形
        fig, ax = plt.subplots()
        plt.get_current_fig_manager().set_window_title('Initial velocity distribution')
        fig.set_size_inches(7, 5)
        plt.xlabel('Velocity [m/s]')
        plt.ylabel('Number of particles')
        if len(ucRatio) > 1:
            plt.title('Initial velocity distribution of ' + formulaString + ' (' + ratioString + ')')
        else:
            plt.title('Initial velocity distribution of ' + formulaString)

        # 更高分辨率的速度轴
        velo_new = np.linspace(np.min(velo), np.max(velo), 100 * len(velo))

        # 遍历唯一元素
        for iElement in range(len(atomUC)):
            # 数据插值以匹配原始x轴更高分辨率的绘图
            interpolator = interp1d(velo, nParticleVeloDist[iElement, :], kind='cubic')
            
            # 生成插值点
            velo_new = np.linspace(min(velo), max(velo), 1000)
            n_particle_velo_dist_new = interpolator(velo_new)
            
            # 绘制结果
            plt.plot(velo_new, n_particle_velo_dist_new, linewidth=3, label=atomUC[iElement].SYMBOL)


        # 查找绘图的最大速度
        i_velo_lim = np.argwhere(nParticleVeloDist >= 0.001 * np.max(nParticleVeloDist)).max()

        # 启用图例
        plt.legend()
        plt.xlim([0, velo[i_velo_lim]])
        plt.ylim([0, 1.1 * np.max(nParticleVeloDist)])
        plt.show()

    return nParticleVeloDist
```            


## 类

### Atom

```python
class Atom:
    """
    ATOM Class holding atomic properties

    Attributes:
        NUMBER (int): Atomic number
        SYMBOL (str): Chemical symbol
        NAME (str): Name
        MASS (float): Atomic mass
        RADIUS (float): Van der Waals radius
        ENERGY_FI (float): First ionization energy [J]
    """

    def __init__(self, number, symbol, name, mass, radius, energy_fi):
        """
        Construct and instance of this class

        Args:
            number (int): Atomic number
            symbol (str): Chemical symbol
            name (str): Name
            mass (float): Atomic mass
            radius (float): Van der Waals radius
            energy_fi (float): First ionization energy [J]
        """
        self.NUMBER = number
        self.SYMBOL = symbol
        self.NAME = name
        self.MASS = mass
        self.RADIUS = radius
        self.ENERGY_FI = energy_fi

Material

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class Material:
    """
    Material Class holding material properties
    """
    
    def __init__(self, formula, elements, amount, volume, density, E_f, E_tot):
        """
        Initializes an instance of this class
        :param formula: Unit cell formula
        :param elements: Array of elements in unit cell
        :param amount: Array of amount of each element in unit cell
        :param volume: Unit cell volume [m^3]
        :param density: Single crystal density [kg / m^3]
        :param E_f: Formation energy per atom [J]
        :param E_tot: Total energy of unit cell [J]
        """
        self.FORMULA = formula
        self.ELEMENTS = elements
        self.AMOUNT = amount
        self.VOLUME = volume
        self.DENSITY = density
        self.ENERGY_FORMATION = E_f
        self.ENERGY_TOTAL = E_tot

PhysicalConstants

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
class PhysicalConstants:
    """
    Class containing physical constants

    Properties:
        BOLTZMANN: Boltzman constant [m2 kg / s2 K]
        PLANCK: Planck constant [m2 kg / s]
        SOL: Speed of light [m / s]
        AMU: Atomic mass unit [kg]
        EV: Electron volt [J]
    """

    BOLTZMANN = 1.38064852e-23
    PLANCK = 6.62607004e-34
    SOL = 299792458
    AMU = 1.66053906660e-27
    EV = 1.602176634e-19

PeriodicTable

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
from Atom import Atom
from PhysicalConstants import PhysicalConstants

class PeriodicTable:
    """
    PERIODICTABLE Class holding atomic properties
    Source: PTable.com
    """

    def __init__(self):
        self.properties = {
            "Li": Atom(3, "Li", "Lithium", 6.94 * PhysicalConstants.AMU, 167E-12, 5.391 * PhysicalConstants.EV),
            "O": Atom(8, "O", "Oxygen", 15.999 * PhysicalConstants.AMU, 48E-12, 13.618 * PhysicalConstants.EV),
            "Ar": Atom(18, "Ar", "Argon", 39.948 * PhysicalConstants.AMU, 71E-12, 15.760 * PhysicalConstants.EV),
            "Ti": Atom(22, "Ti", "Titanium", 47.867 * PhysicalConstants.AMU, 176E-12, 6.8280 * PhysicalConstants.EV),
            "Mn": Atom(25, "Mn", "Manganese", 54.938 * PhysicalConstants.AMU, 161E-12, 7.4343 * PhysicalConstants.EV),
            "Sr": Atom(38, "Sr", "Strontium", 87.62 * PhysicalConstants.AMU, 219E-12, 5.6952 * PhysicalConstants.EV),
            "La": Atom(57, "La", "Lanthanum", 138.91 * PhysicalConstants.AMU, 195E-12, 5.577 * PhysicalConstants.EV),
            "Ni": Atom(28, "Ni", "Nickel", 58.69 * PhysicalConstants.AMU, 135E-12, 7.646 * PhysicalConstants.EV),
        }

    def __getattr__(self, item):
        if item in self.properties:
            return self.properties[item]
        raise AttributeError(f"'{type(self).__name__}' object has no attribute '{item}'")

UnitCells

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
from PeriodicTable import PeriodicTable
from Material import Material
from PhysicalConstants import PhysicalConstants

PT = PeriodicTable()


class UnitCells:
    """
    UnitCells Class holding material unit cell data
    Source: materialsproject.org
    """
    
    def __init__(self):
        self.properties = {
            "TiO2": Material('TiO_2', 
                             [PT.Ti, PT.O],
                             [2,4],             
                             volume=66.292E-30, 
                             density=4130, 
                             E_f=3.464 * PhysicalConstants.EV, 
                             E_tot=56.3693 * PhysicalConstants.EV),
        



        }

    def __getattr__(self, item):
        if item in self.properties:
            return self.properties[item]
        raise AttributeError(f"'{type(self).__name__}' object has no attribute '{item}'")