模拟气相沉积薄膜生长过程

在本实例中你将学到如何使用Virtual Nanolab和Atomistic Toolkit中的分子动力学（MD）模拟来模拟在基底上生长气相沉积薄膜。

气相沉积是一种广泛应用的技术，用以在基底材料上生长晶体薄膜或者无定型固态薄膜。模拟这些过程可以帮助理解其生长微观机制，和过程参数（压力，温度等）对原子结构和动力学的影响。

在本实例中你将模拟在晶体碳化硅（SiC）基底上生长SiC，与文献[1]中的研究类似。你将使用半经验势，相比从头算方法（譬如DFT方法），它能进行更大系统和更大时间规模的模拟。本例中的模拟主要是物理气相沉积（PVD）方法。当然，如果包含了描述化学反应的合适的力场，或者有充足的计算资源来使用从头算方法，化学气相沉积（CVD）过程原则上也是可以进行模拟的。

对于本实例来说，你需要熟悉在分子动力学基础中描述的分子动力学基本功能。你将学到如何使用Python脚本语言和QuantumATK分子动力学常规中的挂钩功能（hook functionality）来运行高级的模拟。

运行一个沉积模拟需要一些不同的技术。与通常的平衡态分子动力学模拟最主要的不同是，随着蒸汽原子或分子的引入，整个系统的粒子数随之增加。在QuantumATK中有两种策略来实现：

1.为每个新引进的原子或者分子运行一次新的模拟。通过在已沉积原子/分子上连续添加新的原子/分子可以实现整个沉积的模拟。

2.将所有需要沉积的原子或分子放在模拟晶胞的库（reservoir）中。对于每个新的沉积过程，从库中取出一个原子把它放在基底上方。

第一个方法的好处是只有真正需要的原子才出现在模拟晶胞中，这提高了模拟效率。然而，由于原子数在变化，在QuantumATK中不可能将整个模拟保存在一个MD Trajectory中用以后续的可视化和分析。所以，本实例选取了第二种方法。当准备模拟时，我们必须注意在库中的原子不会与系统活跃部分有显著的相互作用，尤其是吸附发生的表面。在本实例中，库将会呈现为放在紧挨着基底底部的晶体。

碳化硅晶体有很多构型。在本实例中，基底将会以2H-SiC晶体结构来构造，可以在VNL的database里找到。你当然也可以用别的SiC结构来构建你的基底晶体。

打开Builder，点击Add‣From Database，搜索SiC，添加Wurtzite结构。

对此模拟我们考虑SiC晶体的（0001）表面。使用面板栏（panel bar）右手边的Builders‣Surface（Cleave）插件来构建这个表面。选择（0001）密勒指数和2*2的表面晶格。最后，选择一个Periodic（bulk-like）层外晶胞矢量和包含一层的厚度。点击Finish完成分裂过程。

所得的晶胞在x-y平面上仍具有六角形状，这在分子动力学模拟中会有所不便。 由于六角对称性，表面晶格可以很容易的转化为正交晶格。为此，打开Bulk Tools‣Lattice Parameters，选择Keep‣Cartesian坐标，将B矢量x分量调为0。

使用Bulk Tools‣Wrap插件将原子包入新的模拟晶胞中，你将会看到晶体结构会保持不变。

经过分裂（cleaving）和包裹（wrapping），顶层的终端未必是理想的。在本实例中终端会出现悬挂键，即一重配位的碳原子，这在模拟当中极有可能是不稳定的。为了得到更好的终端，使用Coordinate Tools‣Translate将整个构型平移，使用-0.2Å的z-translation，最后通过Bulk Tools ‣ Wrap将构型包入晶胞里。最终的终端看起来应该如下图，在底部和顶部形成一个更稳定的终端。 接下来通过Bulk Tools ‣ Repeat产生尺寸合适的超胞。在A-B平面使用3*3的重复。C方向矢量必须进行一定的重复，以产生足够晶体层的基底，同时含有足够原子数的库来产生所需厚度的沉积薄膜。对于本例来讲，基底选择3层。为了使模拟时间适度的短，你将使用仅一层的库，从而使C轴总的重复数为4。如果你还想进行实际的生长模拟，你可以增加库的层数（和沉积模拟中的步数）。

现在你需要使用标签（tags）标记系统的不同部分。首先选择系统的上两层作为热化基底，如库（reservoir）方法原理图所示，并将他们标记为“基底”。然后选择接下来的一层晶体层为基底固定的底层，并标记为“bottom”。最后，将下面这些晶体层作为库并标记为“reservoir”。

在构建平板构型之前，块体系统需要使用在实际沉积模拟中所使用的势来进行优化，以得到优化的横向晶格常数。为此将所准备的块体晶体送到ScriptGenerator，然后添加一个New Calculator和Optimize ‣ OptimizeGeometry模块。在New Calculator设置中，选择ATK-Classical并选择Parameter set ‣ Tersoff_SiC_2005[2]。点掉Print和Save并关闭窗口设置。在OptimizeGeometry设置中作如下更改。

运行模拟并将最终优化的构型从LabFloor拖至Builder。

最后，我们需要在块体晶体上产生一个表面。通过增加优化的块体系统晶胞的高度，这将很容易实现。为此，打开Bulk Tools ‣ Lattice Parameters。优化之后，晶胞的一些非对角元可能会不为零。为了方便，将他们重设为零，并保持分数坐标不变。为了调整晶胞高度，确保你选择了Keep ‣ Cartesian坐标。增加C矢量的z分量为100 Å来得到反应表面和库之间足够的真空层。

在运行沉积模拟之前，表面必须被热化到想要的温度。在本例中，你将控制基底温度为2400K，如文献[1]中所建议。将所准备的构型从Builder送到Script Generator。添加一个New Calculator和Optimize‣MolecularDynamics模块。在New Calculator设置中，做与之前优化模拟相同的设置。在MolecularDynamics设置栏中，选择MolecularDynamics‣Type‣NVT Nose Hoover Chain，steps选择100000，log interval 选择5000步，并输入一个合适的轨迹文件名。选择Initial Velocity‣Type‣Maxwell-Boltzmann初始速度为2400K并设置虚拟热浴温度（reservoir temperature）为2400K。（注意在这里reservoir temperature不是指模拟晶胞中的粒子库，而是指虚拟热浴温度，参见分子动力学基础）。点掉Save和Print复选框。

通过点击Add constraints按钮打开约束窗口栏。基底底层和库原子都应该在平衡过程中保持不动。为此选择Constraint‣Fixed约束这两个标记组并关闭约束窗口栏。 通过Job Manager或者在使用atkpython的终端中来运行模拟，使用Movie Tool将最终的快照送回Script Generator。

沉积模拟包含了一个标准的，长的分子动力学模拟，采用了附加功能来对每个分子动力学步骤之后进行自定义操作。（参见参考手册中的分子动力学）。本实例中，一个原子定期的从库中取出并以一个朝向反应表面的随机的热速度放置在表面上方。这个自定义的功能必须在Python程序语言的脚本中定义。

在Script Generator中设置分子动力学模拟，添加一个New Calculator和一个Optimize‣Molecular Dynamics模块，调整计算设置为前面所示。在Molecular Dynamics窗口栏中，选择NVT Nose Hoover Chain类型，步数为800000，log interval为2500步，并将轨迹文件命名。设置虚拟热浴温度为2400K并选择Configuration Velocities为初始速度。点掉Save和Print。将模拟脚本送到Editor以进行必要的修改。

首先你必须引入wrap（）函数从而把坐标包入模拟晶胞。

from NL.CommonConcepts.Configurations.Utilities import wrap

然后你需要定义实施挂钩函数的类型（总体过程类似于例子碳纳米管的拉伸）。首先，定义类型名和构造器方法。

# ------------------------------------------------------------
# Define a hook class to deposit atoms from the reservoir
# onto the active surface of the substrate.
# ------------------------------------------------------------
 
class VaporDepositionHook:
    def __init__(self,
                 deposition_interval,
                 configuration,
                 substrate_atoms,
                 vapor_temperature=300.0*Kelvin):
        """
        Constructor
        """
        # Store the parameters.
 
        # The interval after which a new atoms is deposited.
        self._deposition_interval = deposition_interval
 
        # Store the original coordinates. They will be used to constrain the reservoir atoms.
        self._original_coordinates = configuration.cartesianCoordinates().inUnitsOf(Angstrom)
 
        # The indices of the atoms that should not be deposited
        self._substrate_atoms = substrate_atoms
 
        # The vapor temperature.
        self._vapor_temperature = vapor_temperature
 
        # Deposition count.
        self._deposited_index = 0
 
        # Get the atoms in the reservoir, i.e. all with z-coordinates lower
        # than the fixed bottom of the substrate.
        self._lowest_substrate = self._original_coordinates[self._substrate_atoms, 2].min()
        self._reservoir_indices = numpy.where(self._original_coordinates[:, 2] < self._lowest_substrate - 0.1)[0].tolist()
 
        # Get the lattice vectors.
        cell = configuration.bravaisLattice().primitiveVectors().inUnitsOf(Angstrom)
        self._lx = cell[0, 0]
        self._ly = cell[1, 1]
        self._lz = cell[2 ,2]

在构造器中，你本质上为沉积存储有用的参数，比如从库中取出原子的间隔，和将原子约束在库中的原始坐标。 最后，你需要定义一个类函数，当挂钩函数被调用时访问该类函数。由于我们的挂钩是以Python类来实现的，call() 类函数必须被定义得使它表现的像一个函数。所有挂钩函数必须有 step, time, 和configuration参数。这个类函数应该以如下形式实现。

def __call__(self, step, time, configuration, forces, stress):
        """ Call the hook during MD. """
        # Get the coordinates.
        coordinates = configuration.cartesianCoordinates().inUnitsOf(Angstrom)
 
        # The reservoir is composed of all atoms that are below the bottom
        # layer of the substrate
        # or above the ceiling of the cell.
        self._reservoir_indices = numpy.where((coordinates[:, 2] < self._lowest_substrate - 0.1) |
                                              (coordinates[:, 2] > self._lz))[0]
 
        # Freeze the reservoir atoms and reset their positions at every step.
        velocities = configuration.velocities()
        velocities[self._reservoir_indices, :] = 0.0*Angstrom/fs
        configuration._changeAtoms(indices=self._reservoir_indices,
                                   positions=self._original_coordinates[self._reservoir_indices]*Angstrom)

首先，从当前构型中提取笛卡尔坐标。其次，库中原子是根据它们当前的位置是在基底之下或者在上晶胞边界之外来决定。对于所有的库原子速度设置为零（我们不想让它们移动），尤其是当库被移出原子不断的侵蚀。另外，在沉积模拟的开始，它们的位置被重设为原始位置。

接下来的部分包含了库中的一个新的原子的实际沉积过程，它只在每个沉积间隔的开始执行一次。

 # Check if it is time for the deposition of a new atom from the reservoir.
        if (step % self._deposition_interval) == 0:
            # Get elements and velocities.
            elements = numpy.array(configuration.elements())
            velocities = configuration.velocities().inUnitsOf(Angstrom/fs)
 
            # If the reservoir is empty, do nothing.
            if len(self._reservoir_indices) == 0:
                return
 
            # Decide which element to deposit next.
            possible_elements = [Silicon, Carbon]
            element_index = self._deposited_index % 2
            next_element = possible_elements[element_index]
            reservoir_elements = elements[self._reservoir_indices]
 
            # Get the element candidates.
            possible_atoms = numpy.where(reservoir_elements == next_element)[0]
 
            # If not avaiblable, try the other element.
            if len(possible_atoms) == 0:
                next_element = possible_elements[element_index-1]
                possible_atoms = numpy.where(reservoir_elements == next_element)[0]
 
            # Return back the global atom indices.
            possible_atoms = self._reservoir_indices[possible_atoms]
 
            # Pick an atom at the bottom of the reservoir
            lowest_atom = numpy.argmin(coordinates[possible_atoms, 2])
            next_index = possible_atoms[lowest_atom]
 
            # Place the deposition atom above the surface at a random lateral position.
            new_coords = numpy.array([numpy.random.uniform()*self._lx,
                                      numpy.random.uniform()*self._ly,
                                      self._lz - 15.0])
 
            configuration._changeAtoms(indices=[next_index], positions=new_coords*Angstrom)
            wrap(configuration)
 
            # Draw random velocities for the atom according to
            # Maxwell-Boltzmann at the specified temperature.
            m = next_element.atomicMass()
            sig = self._vapor_temperature*boltzmann_constant/m
            sig = sig.inUnitsOf(Ang**2/fs**2)
            sig = numpy.sqrt(sig)
 
            new_velocity = numpy.random.normal(scale=sig, size=3)
            velocity_norm = numpy.linalg.norm(new_velocity)
 
            # Set the velocity so that it points towards the active surface.
            new_velocity = numpy.array([0.0, 0.0, -velocity_norm])
            velocities[next_index] = new_velocity
 
            # Set the new velocities on the configuration.
            configuration.setVelocities(velocities*Angstrom/fs)
 
            self._deposited_index += 1
 
# ------------------------------------------------------------
# End of hook class definition
# ------------------------------------------------------------

首先我们决定接下来沉积哪种元素。由于硅原子和碳原子交替地沉积，只需要查看当前沉积指数是奇数还是偶数就可以检测到。然后库原子中所有这个元素被挑选出来，选择在库最底的原子。如果出于某种原因库中没有所需的元素类型，我们会选择其他元素来代替。通过构型中的 _changeAtoms() 类函数产生新的随机横坐标和距晶胞上端15.0Å的高度。为了更好的可视化，使用 wrap()类函数将坐标包进晶胞中。为了建模沉积原子的蒸汽温度，按Maxwell-Boltzmann分布选取新的速度，并且为了保证原子驶向表面，速度矢量被重新对齐。为了简化，矢量总是垂直指向表面，但总的来说，不同的冲击角度也很容易实现。使用setVelocities() 类函数来给沉积原子分配新的速度和更新构型的速度。最后，运行沉积指数增加，完成类的定义。 完成VaporDepositionHook 类的定义后，新的实例开始，并指定了当前模拟的参数。这应当在脚本中的Molecular Dynamics模块紧前面完成。

# Initialize the hook class.
deposition_hook = VaporDepositionHook(deposition_interval=5000,
                                      configuration=bulk_configuration,
                                      substrate_atoms=substrate+bottom,
                                      vapor_temperature=2400.0*Kelvin)

为了运行模拟，在Script Generator中预定义的分子动力学模块基础上做两个小改动。在Nose-Hoover chain虚拟热浴类函数的定义中， reservoir_temperatures参数被修改，从而使只有被标记为“基底”的原子群才与虚拟热浴进行耦合。

# Define an NVT thermostat where only the substrate atoms are thermalized.
method = NVTNoseHoover(
    time_step=1*femtoSecond,
    reservoir_temperature=[('substrate', 2400*Kelvin)],
    thermostat_timescale=100*femtoSecond,
    heating_rate=0*Kelvin/picoSecond,
    chain_length=3,
    initial_velocity=initial_velocity
)

这是有必要的，确保气相中的沉积原子的动力学在被表面吸附前不被虚拟热浴所影响（也就是说，沉积原子应该处于NVE系综中）。最后，deposition_hook 项作为逻辑判据必须在 MolecularDynamics()函数中被通过。

md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=bottom,
    trajectory_filename='SiC_deposit_tutorial_md.nc',
    steps=800000,
    post_step_hook=deposition_hook,
    log_interval=2500,
    method=method
)

整个模拟完整的脚本可以在Script_deposit_SiC.py文件中找到。

通过Job Manager或着在一个终端中运行脚本。如此长的分子动力学模拟可能花费数小时来完成。你会注意到一些原子并没有被表面吸附，而是被散射回气相。不过没关系，这些原子向上运动将最终到达库，然后自动地被再一次视为正常的库原子。

下面这些图展示了沉积模拟不同阶段的快照。

最开始，只有少数原子被沉积在表面并被吸附形成单个吸附原子。随着越来越多的原子被沉积，在基底表面上开始形成第一层。第一层，至少部分地，像基底的晶体结构，而接着的几层展现的更像是无定型结构。这主要归结于所选的势（参见开始的小提示），而在某种程度上也归结于高的沉积频率（由 deposition_interval参数所规定）。一个大的 deposition_interval会降低沉积频率，从而使沉积原子的集成更可控。切记，如果你增加沉积间隔，你必须对沉积模拟总步数也有相同的增幅来确保库中所有的原子会被沉积。

注意在薄膜生长的同时，随着越来越多的原子被拿走，平板底部库区域是如何变小的。然而，脚本中特殊的约束保证了消失的结构不会使整个系统不稳定。

这个模拟技术原则上适用于各种材料。当然，你可能会需要在脚本中做一些修改，比如调整元素。除了单原子，你也可以沉积整个分子。对于生产模拟你或许需要增加库的层数来得到一个大的沉积薄膜。你还可以对其他参数比如温度或者沉积间隔进行调整。沉积频率，也就是沉积间隔的倒数，可以通过Hertz- Knudsen方程与蒸汽压建立联系（参见[3]）:

$\frac{dN_i}{dt}=\frac{p_i A}{\sqrt{2 \pi m_i k_B T}}$,

其中$p_i$表示原子/分子i的分压，A表示表面积，$m_i$表示原子/分子质量。

插入模拟中的参数，你将会发现相应的蒸汽压大约为66 bar，远远大于典型的实验的近真空状态。这方面显示了此类模拟的一个显著缺点，由于在多数情况下可达到的模拟时间不允许实际压力值的使用，致使时间分辨解析生长动力学很困难。本质上，这意味着粒子以过高的速率撞击，由于反应热从表面到受恒温器控制的块体基底传的不够快，将导致局部加热的增加。然而，所获的沉积薄膜的结构模型在多数情况下可被视为真实薄膜结构的一个很好的近似。你应当记住实行模拟的条件，并考虑这些条件是否会对沉积结构产生影响。 此外，要注意使用恒定的沉积间隔提供了另一个近似，因为实际过程遵循的是泊松分布而不是均匀分布。虽然泊松分布可以很容易在脚本中实现（使用比如python的 random.poisson函数），当沉积时间间隔选的足够高从而产生想要的实际蒸汽压，选择这种分布才会变得有意义。

[1] (1, 2, 3) Kang, K.-H; Eun, T.; Jun, M.-C.; Lee, B.-J.: Governing factors for the formation of 4H or 6H-SiC polytype during SiC crystal growth: An atomistic computational approach. J. Cryst. Growth, 389, 120 (2014)

[2] P. Erhart and K. Albe: Analytical potential for atomistic simulations of silicon, carbon, and silicon carbide. Physical Review B, 71, 035211, (2005)

[3] K. Reuter: First-Principles Kinetic Monte Carlo Simulations for Heterogeneous Catalysis: Concepts, Status, and Frontiers. in “Modeling Heterogeneous Catalytic Reactions: From the Molecular Process to the Technical System”, Wiley-VCH (2011)

本文翻译：王吉章