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

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

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

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

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

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

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

第一个方法的好处是只有真正需要的原子才出现在模拟晶胞中，这提高了模拟效率。然而，由于原子数在变化，在VNL-ATK中不可能将整个模拟保存在一个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

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

        # 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
# ------------------------------------------------------------