创建无定形材料结构模型

引言

无定形固体材料在广泛的技术领域里的重要性越来越高。晶体结构通常都能够很好的定义和描述，而要对特定无定形材料进行原子尺度的模拟，建立结构模型是相当困难的一件事。本例示范如何使用 VNL 创建无定形结构，并使用 QuantumATK 在不同精度下进行分子动力学模拟。

无定形结构，顾名思义就是组分原子、分子没有任何周期性顺序的结构。无定形结构常见于各种材料中，例如玻璃、聚合物。其力学性质、电子性质，通常与晶体区别很大。一个典型的例子是无定形二氧化硅的热胀系数比对应的石英晶体要小一个数量级。

因此，在无定形材料方面投入了越来越多的人力，以根据需求来定制材料和器件性质。在这方面，原子尺度的模拟能够辅助理解纳米和微观尺度无定形结构行为。但有时候对给定的材料，不容易得到合适的无定形结构模型。和对应的晶体相比，无定形结构并没有特定的单胞或超胞来表示，因为其本质带着很大的随机性。因此对于无定形的情况，通常需要考虑一系列同类结构，但这导致后续的计算量就很大，而往往无法进行了。

本例中，将会学习如何生成代表性无定形结构，并将其作为后续 DFT 或经典力场模拟的输入结构。将会演示如何使用经典力场分子动力学模拟获得无定形材料的起始结构，并通过半经验或 DFT 分子动力学得到改进——如果关心的是电子结构的话。

由于没有一个普遍通用的创建无定形结构的诀窍，对某种给定材料，一个最适合的方案就需要考虑很多方面的因素，例如：有没有合适的经验力场，超胞的尺寸想要多大，可用的计算资源有多少，等等。因此本例只提供产生无定形结构的指南，而不是非常具体切实的方案。对用户手里的实际研究案例，需要酌情调整。

本例模拟无定形二氧化硅，也就是熔融的二氧化硅（石英）。这是应用最广泛的无定形材料，比如用于常规或特种玻璃，或者半导体工业。本例还会介绍如何创建晶体区域和无定形区域之间的界面。

如果你不熟悉分子动力学模拟，那么可以从【教程：分子动力学基础】开始，这样可以帮助你熟悉 QuantumATK 的分子动力学功能。

注意

本例演示在合理的时间范围内，如何通过少量的操作步骤产生无定形结构。对实际中的模拟，要得到结构缺陷较少的可靠的无定形样本，需要的模拟时间就长的多了。比较典型地，一个熔融-退火的循环至少需要几百 ps，甚至几纳秒。

使用经典力场分子动力学生成无定形结构

生成无定形结构最普遍的方法就是仿照实验方法，使用分子动力学模拟熔融晶体之后，再次降温¹⁾,²⁾。虽然冷却速率不可能像实验中的冷却速率那么低，大部分情况下，得到的结果还是令人满意的。产生的无定形结构的质量，当然在很大程度上依赖于力场能否很好的描述这种材料，尤其是偏离平衡态条件时。

设置几何结构

您应该从 SiO₂ 晶体出发，本例使用方石英晶体结构。之所以选用方石英，是因为方石英的晶体结构很容易转换成正交盒子，这样进行模拟会更方便，石英就没有这么方便。另外，方石英的密度与实验测量的无定形二氧化硅密度很接近，为 2.2g/cm³。

打开 Builder，点击 Add ‣ From Database，搜索 SiO2。选择 cristobalite，点击右下角“＋”将其添加到 Builder 窗口中：

点击窗口右边的 Bulk Tools ‣ Supercell，点击 Conventional 和 Transform，那么左边窗口的方石英晶体就变成了立方超胞。

这样的超胞对于无定形结构来说太小，因为这个盒子是满足周期性边界条件的，因此这样实际上还是具有一定晶体的属性。通过创建超胞增大盒子：点击 Bulk Tools ‣ Repeat，设置为 3 x 3 x 3 的重复。这样新的盒子包括 648 个原子，可以用于分子动力学模拟的初始结构。

设置模拟

将结构拖到 ScriptGenerator，添加一个 NewCalculator，以及 Optimize ‣ MolecularDynamics 到 Script 面板。双击 NewCalculator，选择 ATK-Classical 计算器，使用 Pedone_Fe2_2006 力场³⁾。去掉 Print 和 Save 的 ☑️，点击 OK 结束 NewCalculator 设置。

双击 MolecularDynamics，为了熔化晶体，需要将温度设置的非常高，比如 5000K。因为这个初始的模拟主要的目的是为了将原子的排布随机化，因此这一步不需要固定压强。因此 Type 选择 NVT 系综（关于 QuantumATK 分子动力学模拟的系综，更详细的文档参考【分子动力学基础】）。NVT 系综中，Langevin 是一个很好的选择，系统与热浴紧密耦合，因此不会出现某个原子过热而导致系统不稳定的情况。

设置 Steps 为 20000，Log interval 为 5000，保存轨迹的文件名为 SiO2_5000K_low_density_traj.nc。在 Initial Velocity 栏，选择Maxwell-Boltzmann，并设置温度为 5000K。

在 Langevin 栏里面，Reservoir temperature 同样设置为 5000K，其余保持默认设置即可。取消 Save 和 Print 的 ☑️（因为上面的设置里面，原子轨迹已经保存过了）。

使用按钮，将脚本送到 Editor 进行手动修改。

修改脚本

有的时候，尤其是当你试图熔融、随机化相对较小的超胞的晶体结构，即使温度远高于实验熔点也会比较困难。除力场的因素外，这种现象的原因是盒子的周期性边界条件实际上为晶体增加了额外的“稳定”性。一般而言，超胞越小，这种周期边界带来的人为“稳定”效果就越强。另外，要记得 NVT 系综的体积是固定的，因此在高温的时候，体系的压强非常大，这实际上也提升了熔点，对某些情况可能会导致倾向于其他的晶相而不是液相。

修改脚本中定义晶格的部分，如下所示：

# Set up lattice
# Define a scale factor
scale_cell = 1.1
# Scale all 3 cell vectors by this factor
vector_a = [21.498*scale_cell, 0.0, 0.0]*Angstrom
vector_b = [0.0, 21.498*scale_cell, 0.0]*Angstrom
vector_c = [0.0, 0.0, 21.498*scale_cell]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

比例因子 scale_cell 在这里不太重要。如果你发现 5000K 熔化之后还残余部分晶体特征，那么可以适当增大。

把修改后的脚本发送到 Job Manager 并运行。

大约需要几分钟时间，计算就能完成。在 LabFloor 中选择生成的轨迹文件，使用右侧的 Movie Tool 可以看到原子以混乱的方式到处运动，频繁地穿越盒子的边界。为了看的更清楚，可以把动画的最后一帧，通过 Movie Tool 窗口右下角的按钮送到 Builder 里面。

在 Builder 里面，点击 Bulk Tools ‣ Wrap ‣ Apply，可以将所有穿过周期性边界的原子都投射回到盒子中。

最后的结构看起来类似如下：

由于局部区域原子簇的密度发生了这种人为性的降低，所以分布不变得均匀。不过对之前的晶体结构进行完全的随机化这个目的完全达到了。

为了恢复到最初的密度，将该结构发送到 ScriptGenerator，完全一样的重复上面的步骤（准备 NewCalculator 和 MolecularDynamics 脚本，并送到 Editor 中），但比例因子 scale_cell 设置为 1.1，而不像之前那样，为了降低密度而设置成几倍。这个过程将会压缩体系的密度到原先的较为均匀的密度（总体密度实际上没有改变）。温度仍然设置为 5000K、20000 步运行计算，达到平衡，将 Movie Tool 中最后一帧结构送到 Builder 中，并 Wrap 结构。

冷却结构、调整密度

观察生成的结构的局部原子排布，你会看到大量的结构缺陷，例如配位缺失的原子（例如氧原子周围只有一个硅原子）。

这样生成的无定形体系实际上是材料在高温下的液态的一个快照。而最终目的往往是为了生成固态的无定形材料结构，也就是体系温度远低于玻璃化温度T_g。因此，你需要进行一系列的降温模拟，将体系冷却到想要的温度。

尽管真实的无定形结构总是会包含一定数量的结构缺陷，但这种降温模拟的过程应当尽可能慢，以尽可能的减少体系缺陷的数量²⁾。但很遗憾的是，模拟时间的限制导致降温速率与实验上相比还是非常高的。下面的计算中所使用的设置，是为了快速的演示如何生成结构的过程，因此可能不是最优设置。要得到更可靠的无定形结构的话，你可能需要增加每个温度步骤跌模拟时间（也就是增大MD步数）。

另外一个重要的参数，是无定形结构的密度。基本上有两种可能的方法来处理。

一种办法是预先调整体系的体积，从而得到期望的、实验上的无定形密度，并在冷却的过程中固定体积¹⁾。这种方法更适用于单步计算非常耗时的 MD 模拟，例如 DFT-MD 模拟。这种方法会导致最终体系有比较大的残余应力。 

另外一种办法是冷却模拟的过程中，连续地调整密度，也就是将体系与恒压器耦合，恒压器压强作为为环境压强²⁾。这种方式得到的密度与选取的力场的好坏有关系，需要检查密度与实验值是不是一致。如果偏离很大的话，可以考虑固定体积后重新模拟。

本例中，示范的是第二种方法。

把前面模拟得到的最后一帧结构，使用按钮送入到 Builder 中，用于冷却模拟。和前面一样，使用右边窗口中的 Wrap 工具，把所有原子都映射到单胞中。将结构送入 ScriptGenerator 中。添加一个 NewCalculator，以及 Optimize ‣ MolecularDynamics。在 NewCalculator 中使用和前面的模拟相同的设置。

打开 MolecularDynamics 设置，type 选择选择 NPT Melchionna，steps 设置为 5000，log interval 设置为 500，并将轨迹文件命名为SiO2_5000K_NPT_traj.nc。initial velocities 设置为 Maxwell-Boltzmann，5000K。

Reservoir temperature 设置为 5000K，并设置 Thermostat time scale 为 20fs，与 Barostat time scale 一致。设置 Bulk modulus 为 37GPa，这个数值对应着石英的弹性模量。

将脚本送入 Editor。

为了实现自动地完成多次温度降低过程的模拟，可以对脚本中 MD 模拟的部分添加一个 for 循环，如下：

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------
# Define a loop over the temperatures
for t in [5000.0, 4000.0, 3000.0, 2000.0, 1500.0, 900.0, 300.0]:
    initial_velocity = MaxwellBoltzmannDistribution(
        # Change the temperature here:
        temperature=t*Kelvin,
        remove_center_of_mass_momentum=True
    )
    method = NPTMelchionna(
        time_step=1*femtoSecond,
        # Change the temperature here:
        reservoir_temperature=t*Kelvin,
        external_stress=1*bar,
        thermostat_timescale=20*femtoSecond,
        barostat_timescale=20*femtoSecond,
        bulk_modulus=370000*bar,
        initial_velocity=initial_velocity,
        mask=[[True, False, False], [False, True, False], [False, False, True]]
    )
    # Define a new trajectory filename for each temperature
    trajectory_filename = 'SiO2_%iK_NPT_traj.nc' % (t)
    md_trajectory = MolecularDynamics(
        bulk_configuration,
        constraints=[],
        trajectory_filename=trajectory_filename,
        steps=5000,
        log_interval=500,
        method=method
    )
    bulk_configuration = md_trajectory.lastImage()

注意，脚本中，需要循环执行的部分，每行都需要缩进，这是 Python 语言的要求。

将修改后的脚本送入 Job Manager，并运行之。这个任务需要一些时间。

对于这个体系，这个模拟不会有问题。在 NPT 高温模拟中，如果你发现原胞的体积剧烈地不停震荡，或者剧烈地形成真空孔结构，这个时候可以考虑在温度降低到大致为玻璃化温度范围之后，模拟高温、恒容的时候，打开恒压器耦合（也就是 Molecular Dynamics 窗口 stress coupling 的部分）。

分析结构

模拟结束之后，你可以通过 Movie Tool 或者 Viewer 看一看分子动力学轨迹。不过因为很多原子运动越过了盒子边界，所以直接通过这种方式来估计结构比较困难。

大致看看局部的原子排布，你会发现几乎整个结构是由共享顶点的、扭曲很小的四面体 SiO₄ 单元组成的。另外可以看到少量的配位缺陷。

为了定量估计，在 LabFloor 栏选择 SiO2_300K_NPT_traj.nc，并使用 MD Analyzer 工具打开。这样你可以使用各种方法分析 MD 轨迹。对于无定形体系，一般来说会特别关心径向分布函数（radial distribution function，RDF）：点击 Load，增加截断半径 R_max 到 8.0Å，Element 1 和Element 2 分别选择 Silicon 和 Oxygen，点击 Plot 显示这两种元素之间的RDF。

和晶体结构尖锐的单线型峰有所不同，你会发现 RDF 峰被展宽了，这是无定形结构的特征。你可以为任意两种不同的元素对之间的 RDF 作图。

选择 Coordination number 分析可以分析结构缺陷的数量。Element 1 和Element 2 分别选择 Silicon 和 Oxygen，计算硅的氧原子配位数。设置截断半径为 2.5Å，这是 RDF 曲线的第一个峰（半径最小的峰）的位置。点击 Plot 创建配位数分布函数。这样你会看到大多数硅原子的配位数是正常的 4，极少量的硅原子配位数较低。然后切换元素的选择，计算氧原子的配位数。点击 Add Line，在现有的曲线基础上，增加新的曲线。类似地，可以看到氧原子的主要配位数是2，有少量的低配位氧原子。MD Analyzer 提供了大量的各种分析数据，这些分析通过简单的几次鼠标点击就可以得到。大家可以大胆地尝试这些功能。例如，可以看看 Silicon-Oxygen-Silicon 键角的 Angular distribution function，或者中子散射函数。

生成的无定形二氧化硅的密度，可以使用盒子的体积计算出来，盒子体积可以在Builder 里面的 Bulk Tools ‣ Lattice Parameters 查到，并能查到元素的原子个数以及各自的质量。对于本例，得到的密度约为 2.6 g/cm³。这个值比实验值（2.2 g/cm³）略高一点，这应该被认为是力场的缺点导致的。要得到与实验值完全一致的密度值，你可以考虑调整盒子尺寸，并固定体积之后重新冷却——使用 NVT Nose-Hoover Chain 系综，温度线性下降，这样是比较合适的。

最后的结果可以作为输入结构用于后面的静态的或 MD 模拟，以计算其他感兴趣的性质。不过尤其是对于静态计算，一定不要忘记，生成的这个结构，仅仅是无定形结构体系的一个可能的快照而已。要得到更有意义的结果，应该考虑几个不同样本的平均值。

改进无定形结构

对经典力场 MD 模拟得到的结构，可以通过进一步的精细计算改进来克服力场可能存在的不足。如果你接下来是打算使用半经验或量子力学计算，比如 Tight-Bonding(TB)，或者DFT 计算，那么这样做就非常有用。

最简单快捷的方法，是使用 TB 或 DFT 进行结构优化，将结构优化到最近邻的能量最低点。这能改进局部的结构变形，但对可能存在的连接缺陷就很困难。

更高级的改进，可以通过较短的 TB/DFT MD 退火、冷却过程实现¹⁾。因为这种模拟非常耗时，所以只是对于比较小的体系，比如 100 到 200 原子的体系比较合适。

接下来，演示一个简单的使用 ATK-DFT 计算器模拟小的二氧化硅体系的例子。由于DFT-MD 模拟时间所限，要得到可靠的恒压平衡结果很困难，而且计算费事费力，所以本例使用恒容的经典力场和 DFT MD 模拟。密度保持为方石英的值，该值与无定形二氧化硅的实验值非常接近。

设置经典力场模拟

在进行 DFT 模拟之前，首先必须使用经典力场 MD 模拟得到一个无定形的单胞。过程与前面一节所述基本一致。

首先，和前一节一样，准备好方石英晶体盒子，但这里 Builder 里面的 Repeat Tools 中，只使用 2x2x1 的超胞。这样的到一个盒子其中包含 96 个原子。盒子的尺寸是 14.33 x 14.33 x 7.17 Å³，也就是 z 方向比其他方向薄——这有可能会加强这个方向的有序度。为了避免这种缺点，我们保持体积不变，将盒子转换为立方形的。最终盒子的尺寸为 11.375 x 11.375 x 11.375 Å³。在高温模拟中，体系是流体状态，你需要通过两步来实现盒子形状的转换。

第一步包括设置模拟，与前面的设置方式一致，也就是添加 ATK-Classical calculator，力场设为 Pedone_Fe3_2006，MolecularDyncamics 中，type 设置为 Langevin，step number 设为 20000，温度设为 5000K。将产生的脚本送入 Editor。为了达到更低的密度，这样能够促进熔化，增大盒子 C 方向长度到最终的 11.375Å：

# Set up lattice
vector_a = [14.332, 0.0, 0.0]*Angstrom
vector_b = [0.0, 14.332, 0.0]*Angstrom
# Previous C-vector:
# vector_c = [0.0, 0.0, 7.166]*Angstrom
# New C-vector:
vector_c = [0.0, 0.0, 11.375]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

运行任务，并将 Movie Tools 中的最后一帧送到 Builder 中，Wrap 结构，设置与上一节中所示相同的模拟脚本，并送入 Editor。

盒子体积必须减小，以恢复原先的方石英密度。不重设盒子的 C 方向长度，现在改为设置 A、B 方向的长度也设置为 11.375Å，这样的到一个立方的盒子。运行任务。

将最后一帧结构送入 ScriptGenerator，添加与上一节相同的 ATK-Classical 的 calculator，以及 Optimize ‣ MolecularDynamics。在 MolecularDynamics 中，type 选为 NVT Nose-Hoover Chain。这种系综能够允许温度线性地下降。设置如下图所示：

这些热力学设置使得热浴温度从 5000K 在 50ps 中下降到 300K。提交计算任务。将最后一帧，通过 Movie Tool ‣ Send To ‣ Builder ‣ Wrap ‣ Send To，送到 ScriptGenerator。

设置DFT模拟

接下来使用 ATK-DFT 计算器进行 MD 模拟来退火，并将其冷却到室温。相对于能量最小化而言，MD 模拟在计算时间上相当昂贵。为了生成无定形结构，并不一定需要得到精确到收敛能量（能量极小值点）。因此，本例中建议到参数设置，在精度上略低于默认参数设置，以提高计算效率。开始 DFT 模拟之前，应该清楚这些模拟是非常耗时的，可能长达数天，因此推荐在集群上来完成。

设置模拟脚本：添加一个 NewCalculator、两个 two Optimize > MolecularDynamics。在 NewCalculator 中，选择 ATK-DFT calculator。打开 calculator 的详细设置，调整一些参数设置：在 Iteration control parameters 中设置 Tolerance 为 0.001 Hatree；在 Basis set/exchange correlation 中将交换相关改为 GGA；将 basis sets 从 DoubleZetaPolarized 改为 SingleZetaPolarized。

考虑拥有的硬件资源，尤其是要进行生产模拟的时候，可以选择更精确的参数设置。

点击 OK，关闭 NewCalculator 设置窗口。

双击打开第一个 MolecularDynamics，设置退火模拟温度为 3000K，选择 type 为 Langevin，设置 number of steps 为500，log interval 为 20，选择合适的轨迹文件名字，例如 traj-sio2-DFT-MD-3000K.nc，initial velocities 选择 Maxwell-Boltzmann，温度设为 3000K，reservoir temperature 也设为 3000K。去掉 Save 和 Print 的勾，关闭设置窗口。

双击打开第二个 MolecularDynamics，在该过程里面，需要将温度迅速地从 3000K 降低到 300K。同样地，Type 使用 Langevin。number of steps 设为1000，log interval 设为 20。initial velocities 同样设为 Maxwell-Boltzmann，温度为 3000K。在该过程中要进行冷却，因此设置 reservoir temperature 为 300K，这样系统从一个非常热的热浴变为室温热浴。考虑到要节省模拟时间，这种方法比 NVT Nose-Hoover chain 系综线性降温的效率要高。

模拟结束后，第二个过程的动画类似如下：

温度从初始 3000K 降低到室温 300K。这个降温过程仅为 1ps，降温速度非常快，因此和实验过程相比，包含一定的近似。

用户可以使用最后的结构进行 DFT 结构优化计算，并计算用户感兴趣的性质。

创建晶体／无定形界面

晶体／无定形界面结构是一个研究热点，界面附近的局域性质尚不清楚。使用 VNL 的 Builder 可以很容易地创建这种结构。

如果晶体与无定形区域是同一种材料，你可以首先准备好晶体结构。之后通过Builder 里面的 Surface（Cleave）来创建构成界面的表面，并创建一个拷贝。第二个结构可以用来产生无定形结构——因为这样，它们的表面已经是精确匹配了。

下面的例子考虑石英和无定形二氧化硅之间的界面。初始结构通过 Surface（Cleave）来创建，选择的是（110）面，并将表面的基矢分别设为原值的4倍和3倍，得到新的原包。创建 Periodic (bulk-like) 超胞，层数设为 12 层。对结构以及晶格矢量进行优化，calculator 仍然使用 ATK-Classical Pedone_2006_Fe3，并将最后得到的结构保存为 NetCDF 文件。将该文件拖到 Builder 中，并使用该结构作为无定形结构的初始结构。

使用前面的流程制备无定形结构盒子，不过只设置 z 轴 pressure coupling，这样盒子的横向尺寸会保持不变。单向的 pressure coupling 是在MolecularDynamics 设置中选定：使用 NPT Melchionna，在 Stress coupling 的 3×3 勾选栏中，只勾选 zz。

创建好无定形超胞之后，使用 Builder 中 Builders ‣ Interface tool，将优化好的晶体和无定形结构结合到一起：简单地从 Stash 里面将两个结构分别拖到 Interface 工具里面的两个窗口里面即可。Interface 工具会自动检测两个结构的合适的结合方式。

点击 Create 按钮，完成创建。将该界面结构送入 ScriptGenerator。

放大界面结构，可以看到大量的缺陷，比如悬挂键。大多数这样的缺陷会被经典力场的退火模拟消除掉，即使力场并不适合用于化学反应。

添加一个 NewCalculator，一个 Optimize ‣ OptimizeGeometry，一个Optimize ‣ MolecularDynamics。设置想要的 calculator（本例中为Pedone_2006_Fe3-potential），双击打开 OptimizeGeometry 设置。

创建好界面之后，有可能有很多人为原因导致的距离过近的原子，这会导致后面的 MD 模拟不稳定。要去掉界面上这些巨大的力，需要先进行结构优化。Force tolerance 不需要很低，可以设置为 0.5eV/Å，Maximum number of steps 设为 1000，去掉 Print 的勾选，output file 命名为 interface_optimized.nc，并关闭 OptimizeGeometry 设置窗口。

双击打开第一个 MolecularDynamics 设置窗口。在 MD 退火中，需要固定晶体的中央区域的原子，这样保证晶体区域不至于被熔化。因此需要设置如下：点击 Add constraints 打开 Constraints Editor。在弹出窗口内通过鼠标画框选中晶体的中间区域的原子，点击 Add tag from selection，constraints 选择 Fixed：

关闭 Constraints Editor。

在 MolecularDynamics 设置窗口，选择 NVT Nose-Hoover Chain，设置 step number 为 50000，设置 log interval 为 2000。首先应该对结构在一个足够高的温度下退火，这样能够使得原子在局部进行重排，促进界面附近的配位缺陷消除，但温度也不能太高，不能把全局结构熔化了。对本例，1500K 是合适的。在initial velocity 和 reservoir temperature 中设置成该温度。

接下来将 Optimize ‣ MolecularDynamics，设置为冷却模拟。打开Constraints Editor，前面选中的原子团 “Selection 0” 应该已经出现在列表中了。为这一组原子选择 Fixed，关闭窗口。MolecularDynamics 窗口中，选择NVT Nose- Hoover Chain 系综，50000steps，log interval 为 2000，Configuration velocities 设置为 initial velocities。Reservoir temperature 设为 1500K，Final temperature 设为 300K，这样通过模拟降温达到室温。

把脚本送入 Job Manager，运行脚本。

之后还应该进行一个 NPT 模拟来得到最终结构。这次不限制原子。使用 NPT Melchionna 系综，取消所有的 stress coupling 勾选，只保留 zz，这样只保证盒子的纵向垂直于界面。运行 300K、50000 步，最终得到室温的无定形／晶体界面结构。

其他例子

最后，这一节在技术上，提供其他产生无定形结构的例子，包括三氧化二铝、二氧化钛、二氧化铪。

Al2O3

创建三氧化二铝无定形结构，基本流程可以参考文献³ 和⁴⁾。使用相同的经典力场，可以选择 Matsui_CaMgAlSiO_1994，或者可以使用比较早的不包含该力场的 QuantumATK 版本，下载该操作向导附录的python文件。如果将其拷贝并粘贴到脚本的最前面，模拟需要的力场就有了。

为了设置这个体系，你需要导入 α-Al₂O₃，因此需要下载文件 al2o3_corundum_ams_data.zip。在 Builder 里面，点击Add ‣ From Files，选中该文件。为了更方便的控制 MD 模拟，三角晶系盒子需要改为立方盒子。有几种方式。可能的一种：Builders ‣ Surface(Cleave) tool。选择（10-10）米勒指数点击 Next 并选择 3×1 表面晶格。点击 Next 选择 Periodic (bulk-like) ，厚度为 4 层。点击 Finish 的到正交的超胞结构，包含 360 原子。

现在可以进行和上面相似的过程进行模拟，唯一的差别是，冷却过程中，密度固定到 3.175 g/cm³，这是文献⁴⁾建议的。整个流程如下：

进行 5000K 高温模拟 MD 模拟，并降低密度;
使用 NVT Nose-Hoover Chain 系综，冷却系统到 3000K;
压缩体系到希望的无定形密度;
退火到 3000K;
缓慢冷却体系到室温;
在室温进行最后的模拟，计算物理量。

将正交盒子送到 ScriptGenerator，设置基本的 MD 模拟脚本。添加一个 NewCalculator，两个 Optimize ‣ MolecularDynamics。

在 NewCalculator 中选择 ATK-Classical 以及 Matsui_CaMgAlSiO_1994 力场。

在第一个 MolecularDynamics 设置系综为 NVT Nose-Hoover Chain，150000 steps，log interval 为 5000，取消 Save trajectory 的勾选。initial velocities 设置为 Maxwell-Boltzmann-distribution，5000K。设置reservoir temperature 为 5000K，最终温度设为 3000K。其他参数保持默认即可。最后取消Save和Print的勾选。关闭窗口，将脚本送入 Editor。

第二个 MolecularDynamics 设置类似。但 steps 设置为 300000，reservoir temperature 设为 3000K，最终温度为 300K。

为了人为地降低密度，需要再次将三个晶格常数改为原先的 1.1 倍。

第一个 MolecularDynamics 产生 3000K 的液态结构，密度被降低。基于这个体系，密度需要被调整到目标密度 3.175 g/cm³。通常收工计算微观单胞的宏观密度，以及找到晶格矢量的精确比例因子比较麻烦。通过插入如下QuantumATK单元引擎代码到脚本中，第二个 MolecularDynamics 之前，可以很方便的实现（也就是bulk_configuration = md_trajectory.lastImage()这一行之后的内容）：

# Get the volume and the lattice vectors of the cell.
V = bulk_configuration.bravaisLattice().unitCellVolume()
lattice_vectors = bulk_configuration.bravaisLattice().primitiveVectors()
# Get the coordinates and the elements.
fractional_coordinates = bulk_configuration.fractionalCoordinates()
elements = bulk_configuration.elements()
# The total and species-resolved number of atoms.
N = len(elements)
N_Al = len(numpy.where(numpy.array(elements)==Aluminum)[0])
N_O  = len(numpy.where(numpy.array(elements)==Oxygen)[0])
# Calculate the current density
density = (N_Al*Aluminium.atomicMass() + N_O*Oxygen.atomicMass())/V
# Get the density in g/cm**3
density = density.inUnitsOf(kiloGram/Meter**3)/1000.0
target_density = 3.175
# Calculate the scaling factor for the lattice vectors.
scale_factor = (density/target_density)**(1.0/3)
# Scale the lattice vectors.
lattice_vectors = lattice_vectors*scale_factor
# Define a new cell.
new_lattice = UnitCell(lattice_vectors[0], lattice_vectors[1], lattice_vectors[2])
# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=new_lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
# The calculator is still defined and can readily be attached.
bulk_configuration.setCalculator(calculator)

这段代码会提取体积，计算密度，以及得到精确目标密度 3.175 g/cm³ 的比例因子。最后晶格矢量被按比例放大，产生了新的结构。第二个冷却过程 MD 模拟将会基于这个结构。

保存整个脚本为 python 文件（最终的脚本范例，见附件中），运行该脚本。大约需要几个小时才能完成。完成模拟之后，可以如常分析结构。

完整的脚本：script_al2o3_amorphous.zip。

TiO2

模拟无定形二氧化钛，需要使用 Matsui and Akaogi 力场⁵⁾。这个力场已经被大量地用于二氧化钛的晶体和无定形结构模拟，得到的性质与 DFT MD 模拟非常一致。这个力场目前VNL中没有，但可以在脚本里面非常简单的设置。直接下载脚本：matsui_tio.zip。拷贝并粘贴到后面将要提到的过程中产生脚本的顶部。

设置体系：打开 Builder，使用 Add ‣ From Database 添加 TiO2 的金红石晶体结构。使用Bulk tools ‣ Repeat 功能创建 4x4x6 重复，包含 576 原子。将结构送入 ScriptGenerator 并设置和三氧化二铝相同的基本的 MD 脚本。将脚本送入 Editor。

找到calculator的定义部分，如下：

potentialSet = Matsui_TiO_1991() # or any other potential class.
calculator = TremoloXCalculator(parameters=potentialSet)

将 potentialSet = …整行替换成上面的附件 Matsui_TiO.py 的内容。

当调整体积以便得到希望的密度的时候，确保你选用的脚本里面使用的是正确的化学计量法，正确的元素质量，以及正确的目标密度。无定形的二氧化钛的密度为3.80 g/cm³⁶⁾。

完整的脚本下载：script_tio2_amorphous.zip。

HfO2

对于二氧化铪，可以使用 Wang et al 力场⁷⁾，在 QuantumATK 中名为Wang_HfOZr_2012。

同样的，可以从 database 中添加 HfO2 结构。通过 Bulk tools ‣ Supercell plugin，点击 Conventional 和 Transform，将 FCC 晶格转为常规单胞。对体系进行 3x3x3 重复，得到合适的超胞尺寸。

和前面两个例子完全一样地设置脚本。确保选取的体系和元素参数一致，无定形结构密度使用 7.97 g/cm3 ⁷⁾。

完整的脚本下载：script_hfo2_amorphous.zip。

参考

¹⁾

E. A. Chagarov and A. C. Kummel: Generation of Realistic Amorphous Al2O3 And ZrO2 Samples By Hybrid Classical and First-Principle Molecular Dynamics Simulations. ECS Transactions 16, 773, 2008

²⁾

D. J. Cole, M. C. Payne, G. Csányi, S. M. Spearing, and L. Colombi Ciacchi: Development of a classical force field for the oxidized Si surface: Application to hydrophilic wafer bonding. J. Chem. Phys. 127, 204704, 2007

³⁾

A. Pedone, et al: A New Self-Consistent Empirical Interatomic Potential Model for Oxides, Silicates, and Silica-Based Glasses. J. Phys. Chem. B 110, 11780, 2006

⁴⁾

G. Gutierrez and B. Johansson: Molecular dynamics study of structural properties of amorphous Al2O3. Phys. Rev. B 65, 104202, 2002

⁵⁾

M Matsui and M. Akaogi: Molecular dynamics simulations of the structural and physical properties of the four polymorphs of TiO2. Mol. Sim. 6, 239, 1991

⁶⁾

D. Mergel et al.: Density and refractive index of TiO2 films prepared by reactive evaporation. Thin Solid Films 371, 218, 2000

⁷⁾

Y. Wang, F. Zahid, J. Wang, H. Guo: Structure and dielectric properties of amorphous high-k oxides: HfO2 , ZrO2 , and their alloys. Phys. Rev. B 85, 224110, 2012

目录