虽然 LDA 和 GGA 计算，对于 Si（以及大部分其他半导体材料）得不到合适的带隙，但这些“简单”的泛函却能够很好的预测几何结构性质。在本文中，将演示如何使用 QuantumATK 研究所谓的Si（100）表面不对称二聚重构。

本文着重关注物理以及如何正确的设置模型，而非如何操作 VNL。本文假定读者已经有了建模的经验（例如如何切开表面），以及设置计算、查看结果，因此提到具体步骤，不再详尽的解释。

1，启动 VNL 并打开 Builder

2，从数据库添加 alpha-硅晶体结构（Add ‣ From Database）

3，打开Builders ‣ Surface (Cleave)

• 通过输入对应的米勒指数，沿 100 面切开晶体
• 因为二聚重构只发生在 2×1 超胞，所以需要创建较大的表面单胞，而不是用最小表面单胞。因此在第二页窗口中设置$V_2=2u_2$

在表面建模工具的最后一页选择 Slab 厚度为 4.5 层（默认真空厚度是合适的）

4，将结构沿所有方向居中（Coordinate Tools ‣ Center）

5，选择 z 坐标值最小的两个 Si 原子，点击几次 Rattle 工具 ，对这两个原子的位置产生随机的小扰动。这是因为完全对称的状态是超稳定的结构，最好不要从这种状态开始优化

6，选择 Z 坐标值最大的两个 Si 原子，使用 加氢钝化，从而使得所有的 Si 原子饱和

生成的结构如下图所示：

• 依次插入 New Calculator 和 OptimizeGeometry。
• Calculator 中除了 k-point sampling 之外所有参数使用默认值。设置 k-point sampling 为 9×9（当然 c 方向一个 k 点就够了，因为是层状结构）—— 精度非常重要，因为对称、不对称二聚状态的能量差很小¹⁾。
• 优化，限制住底部的两个 Si 原子以及 4 个 H 原子。如果可能的话，限制弛豫的自由度总是很重要的，这样能够投影为纯的平移和转动。
• 确保勾选了Save trajectory，并设置轨迹文件的名字

• 设置输出文件的名字，现在已经准备好运行了。保存脚本！

将脚本送入 Job Manager，或者将其上传到集群之后运行。该计算并行化非常好，根据 MPI 并行节点数差别，计算耗时约为 1 个小时或几个小时。

计算完成之后，选择NetCDFVNL中的文件，查看第二个Bulk configuration，ID为gID001——这是优化之后的结构。 将其拖到Builder，可以立即看到得到的结构确实是不对称二聚态。鼠标划框选取相关的原子，点击Coordinates－Z-matrix即可显示键长、键角。分别查看原始的结构和优化之后的结构，进行对比。

分析结构弛豫的轨迹也很有意义。选择LabFloor中，你之前保存的轨迹文件，点击Viewer大卡，可以看到优化的动画过程（点开大图）:

如果我们对弛豫过程对每一步中的总能量和力作图，可以看到总能量或多或少是在持续下降的，这符合我们采用的BFGS能量最小化方法的期望。BFGS的初衷是为了能量的最小化，而不是仅仅跟随力的方向变化。下面的小脚本能够读取轨迹数据，并为每一步的总能量和力最大值作图：

traj = nlread('si_100_traj.nc', Trajectory)[0]
energies = []
max_forces = []
for i in range(len(traj)):
    energies.append(traj.imageEnergy(i).inUnitsOf(eV))
    forces = traj.imageForces(i).inUnitsOf(eV/Angstrom)
    force_magnitude = (forces**2).sum(axis=1)**0.5
    max_forces.append(numpy.max(forces))
 
# Plot the energy and max. force curve using pylab.
import pylab
 
pylab.figure()
pylab.subplot(2, 1, 1)
pylab.plot(range(len(traj)), energies)
pylab.xlabel('Optimization step')
pylab.ylabel('Energy (eV)')
pylab.grid(True)
 
pylab.subplot(2, 1, 2)
pylab.plot(range(len(traj)), max_forces)
pylab.xlabel('Optimization step')
pylab.ylabel('Maximum Forces (eV/Ang)')
pylab.grid(True)
 
pylab.show()

作图结果大致如下图所示，精确的曲线和顶部的Si原子初始结构中，精确的位移量有关系。

结合上图以及轨迹动画，可以了解很多信息。首先，除非力的收敛标准小于0.2 eV/Å，不然的话，你会以为对称的二聚态是能量最小值，因为优化的过程中是经过了这个态的！在轨迹动画中也能看到这一点。虽然此处得到这样结果，但实际上，众所周知如果不是使用了非常精确的方法（例如足够多的k-points），以及足够大的真空的话，对称二聚态在能量上会更占优势。

特别地，大约在step 13（见上面的轨迹图所示）达到对称二聚态这个局域最低点，但能量最小化算法仍然试图进一步降低能量，这时候力再次增大，走向不对称的结果。Step 15之后，真正得到不对称结构，在大约Step 20的时候，最终形成不对称二聚体。

可以使用下面的脚本来对原子层的垂直弛豫进行可视化：

traj = nlread('si_100_traj.nc', Trajectory)[0]
number_of_steps = len(traj)
z_coordinates = []
 
for i in range(number_of_steps):
    coordinates = traj.image(i).cartesianCoordinates().inUnitsOf(Angstrom)
    # Append all z-coordinates, apart from the fixed bottom layer.
    z_coordinates.append(coordinates[:16, 2].tolist())
 
# Convert to numpy array.
z_coordinates = numpy.array(z_coordinates)
# Plot the z-coordinates using pylab.
import pylab
 
pylab.figure()
# Loop over all unconstrained atoms.
for i in range(16):
    pylab.plot(range(number_of_steps), z_coordinates[:, i], '--*')
pylab.xlabel('Optimization step')
pylab.ylabel('Z-coordinates (Ang)')
pylab.grid(True)
 
pylab.show()

有一点非常重要：不对称的出现，是由长程作用导致的，这也是短程方法，比如tight-binding或者经典力场只能弛豫到对称二聚体到原因。

本例演示了QuantumATK能够通过DFT结构优化，成功地预测相当复杂的表面结构重构。Si(100)面首先形成一个对称二聚体，它反过来诱导出电子层面的不对称，从而我们在再次优化结构之后，找到基态不对称结构。这种不对称是一种长程效应，来自表面以下的少量几层原子。二聚体Si晶格在二聚体下面被向下压缩。