在本实例中你将学到如何使用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

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

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

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