这是本文档旧的修订版!
在这篇教程中,你会学习到如何利用“刚性”(rigid body)结构限制和简单的优化步骤,来快速、可靠的弛豫一个器件结构的内坐标。
器件结构内坐标的弛豫,是一个冗长费时的过程。主要有两个原因:
双电极体系由三个区域组成:左右两个电极,以及一个由电极扩展层和中间散射区组成的中心区域。在第一性原理计算中,每一部分的几何结构需要单独优化,但是同时中心区域的两电极扩展层必须和相应的电极匹配。
在第二步中,我们必须在确保电极扩展层内坐标不变的前提下最小化器件总能量。在ATK里至少有两种方法实现。
首先是完全弛豫电极结构,并在此基础上构建双电极器件结构。之后,我们需要取出并弛豫中心区域的晶胞,但是必须对电极扩展层的原子进行“刚性”结构限制。之后可以恢复添加两端电极,重新得到器件(或是界面)结构。这个方法计算效率比较高,通常可以得到大部分的弛豫结构。我们称这种方法为“Bulk Rigid Relaxation (BRR)”。 BRR方法不计算整个器件的总能量,因此重新得到的器件结构可能不是能量最低的几何构型。部分电子态性质可能会对这种差异敏感。
更有力的方法则必须使器件处于能量最小化状态(相对于内坐标和中心区域长度),即在改变中心区域长度的同时进行整个双电极器件的计算。这是一个一维原子链最小化的问题,在每一步都针对原子位置作弛豫。因此,我们称之为1DMIN。为此,我们发展了用ATK进行这种最小化的方案,使之尽可能有效地处理这个问题,但是这种方法的计算量远大于BRR方法。因此,我们建议在对中心区域使用BRR弛豫之后,再进行1DMIN计算。
这篇教程介绍了利用BRR和1DMIN方法对电极和中心区域的弛豫。教程中的实例是用经典力场方法对银-金接触表面进行弛豫。利用经典势可以使计算更快捷,同样的方法可以适用于DFT计算。通常,我们会建议利用DFT计算来检验ATK-Classical计算的结果。
提示 我们对上面的步骤作了一个图形化的解释,可以在此处下载相关的文献 Infographics
生成一个“device_relaxation”的新项目。我们从一个未弛豫的双电极结构开始。其结构可以参照此处教程 Building an interface between Ag(100) and Au(111),或者从此处下载结构脚本如何弛豫器件体系的几何结构。 在项目中保存py文件,并将其拖拽至Builder中。 两电极结构如下图所示。它由左电极Ag(100)、中心区域、右电极Au(111)构成。利用插件,将器件分割成左电极、右电极和中心区域。
电极的几何结构是由Ag(100)和Au(111)构成(参考Building an interface between Ag(100) and Au(111))。其左右晶胞的结构采用了实验上的几何数据。我们假设实验结构也是理论上能量最低的结构。(实际理论研究中需要弛豫晶胞以验证这一点)
在接触器件构建中,我们选择Strain second surface(对第二个表面施加应力)处理表面晶胞。这意味着在形成界面的过程中,Ag(100)结构将会保持实验上的结构,而Au(111)表面将在x、y方向上将被拉伸。当晶胞在x、y方向上拉伸之后,它将按照材料的泊松比在z方向上弛豫。
接下来,我们将把这种弛豫应用于Au(111)晶胞,方法是在保持x、y方向固定的情况下对z方向上的应力进行弛豫。获得的Au(111)结构将会应用到中心区域的金原子中。
打开Builder并重命名右电极结构为Au111,双击激活结构,并发送至Scripter,添加:
打开New Calculator,选择ATK-Classical计算方法,并选择“EAM_Zhou_2004”势。点击OK。 下一步,点开 OptimizeGeometry :
发送脚本至 Job Manager 并进行计算,计算耗时不到 1分钟,当执行完之后,检查log文件。 我们将会在 LabFloor看到一个“Au111.nc”的nc文件,它包含了弛豫前后的两个结构的结构数据。选中两个结构,并拖拽至Viewer。在Viewer中,鼠标选中晶胞,我们会发现,弛豫前的结构,C方向矢量为7.064 Å ,弛豫后为7.011 Å。这个结果同样可以在“Au111.log”文件中获得。 弛豫后的结果压缩了0.75%。为了应用到中心区域,打开Builder,双击打开中心区域结构。
现在,我们仅需要进行中心器件弛豫。当进行几何优化时,我们希望电极扩展层可以自由移动,因此我们需要在两端添加真空层: 鼠标点击绘制窗口的背景区域,确保没有原子被选中。点开右侧的Bulk Tools ‣ Stretch Cell ,利用C-vector slider,延长晶胞30%的长度,或者简单将因子设置为1.3。然后点击Apply应用,晶胞就会相应的增大。 点开Coordinate Tools ‣ Center,取消勾选A、B,点击Apply,可将中间区域制动设置到中心。
我们要固定电极扩展层的原子,但为了在后面便于使用“Device from bulk”工具添加电极时确定电极的长度,我们需要比常规的电极扩展层多固定一层。为此,选中最外层的5层银原子,点开Selection Tools ‣ Tags。在“Click here to write a new tag…”输入框中,输入Left。
接下来,同样选择最外层的4个金原子, 标记为Right.
发送结构至Script Generator,设置如下参数:
注意 可以对一个标记的区域选择“Fixed”,另一个选择“Rigid”,这样可以减少BFGS弛豫步骤,从而节省计算时间。
发送脚本至Job Manager,并运行计算。计算非常快,一旦计算结束,我们将会看到central.nc文件出现在LabFloor上。
可以从Viewer中看到优化的轨迹动画。
优化后的中心区域结构保存在central.nc文件中的gID0002中。这是我们对中心区域结构的最初猜测。将其拖拽至Builder中。
打开Device Tools ‣ Device From Bulk,将bulk结构转换为器件。保留默认参数即可。
结构如下图所示,右键或F2重命名为“device”。
我们现在可以搜索器件全局最小能量结构了。通过Save as,将器件结构保存为NetCDF文件,并命名为“device_in.nc”。
现在我们将用SciPy minimization程序搜索可以使能量最小化的中心区域长度。我们将从device_in.nc的器件结构开始。 我们需要两个脚本: device_relaxation/device.py和 device_relaxation/optimize.py。将两个脚本和device_in.nc放在同一个目录下。第一个脚本是读取结构并设置计算参数,第二个脚本包含了计算所需的类和函数的定义。
注意 利用device.py 计算其他体系,需要修改脚本,设定相应的势以便能够处理新体系。
两个脚本都可以在VNL的Project Files中看到。拖拽device.py至Jobs,并运行脚本。
一旦计算完成,将会出现一个数据图窗口。蓝色点是弛豫结构的总能量。晶胞的最初长度是28.51 Å,其他的蓝点是算法模拟出来的结果。其中位于28.58 Å的红色点为中心区域能量最低的点。优化程序只简单地绘制出所示的势能曲线,并找到它的最小值。
手动关闭弹出的窗口可以结束脚本运行。
弛豫器件结构保存在输出文件“device_out.nc”,并包含总能量。红色点处的结构(最终找到的能连个最低点)是最后一个器件结构(id最大的那一个)。
在Viewer中,能量最小点的结构如图所示:
这个例子里,1DMIN弛豫后结构与BRR弛豫的结构相比较,变化十分细微。这可以看出,通常情况BRR方法的结果已经足够精确,并不需要再进行1DMIN的步骤。1DMIN仅用于极其精确的计算。
下载pdf, Infographics。点击顶部按钮查看图片。
本文翻译:闫强