这是本文档旧的修订版!
虽然LDA和GGA计算,对于Si(以及大部分其他半导体材料)得不到合适的带隙,但这些“简单”的泛函却能够很好的预测几何结构性质。在本文中,将演示如何使用ATK研究所谓的Si(100)表面不对称二聚重构。
本文着重关注物理以及如何正确的设置模型,而非如何操作VNL。本文假定读者已经有了建模队经验(例如切开表面),以及设置计算、查看结果,因此提到具体步骤,不再非常详尽的解释。
1,启动VNL并打开Builder
2,从数据库添加alpha硅晶体结构(Add ‣ From Database)
3,打开Builders ‣ Surface (Cleave)
* 通过输入对应的米勒指数,沿100面切开晶体 * 因为二聚重构只发生在2x1超胞,所以需要创建较大的表面单胞,而不是用最小表面单胞。因此在第二页窗口中设置V<sub>2</sub>=2u<sub>2</sub>
在表面建模工具的最后一页选择Slab厚度为4.5层(默认真空厚度是合适的)
4,将结构沿所有方向居中(Coordinate Tools ‣ Center)
5,选择z坐标值最小的两个Si原子,点击几次Rattle工具,对这两个原子的位置产生随机的小扰动。这是因为完全对称的状态是超稳定的结构,最好不要从这种状态开始优化
6,选择Z坐标值最大的两个Si原子,使用加氢钝化,从而使得所有的Si原子饱和
生成的结构如下图所示:
1,依次插入New Calculator和OptimizeGeometry。
2,Calculator中除了k-point sampling之外所有参数使用默认值。设置k-point sampling为9×9(当然c方向一个k点就够了,因为是层状结构)——精度非常重要,因为对称、不对称二聚状态的能量差很小1)。
3,优化,限制住底部的两个Si原子以及4个H原子。如果可能的话,限制弛豫的自由度总是很重要的,这样能够投影为纯的平移和转动。
4,确保勾选了Save trajectory,并设置轨迹文件的名字
5,设置输出文件的名字,现在已经准备好运行了。保存脚本!
将脚本送入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或者经典力场只能弛豫到对称二聚体到原因。