磁性隧道结中的自旋输运

本实例将为您展示如何对磁性隧道结（MTJs）的电子输运性质进行模拟和分析（通常用于自旋电子学应用的研究）。你将考虑一个Fe|MgO|Fe磁性隧道结，它是一个比较复杂的自旋极化体系，MacLaren和他的合作者们首次研究了这个体系[BZSM01]。

您将使用QuantumATK来学习共线/非共线自旋相关输运性质，包括电子传输，隧穿磁阻和自旋转移矩。你也将变得熟悉AdaptiveGrid方法，它是对于k点取样的传统Monkhorst-Pack型格点的一种自适应选择。

使用在Builder插件中的Magnetic Tunnel Junction来创建一个Fe|MgO|Fe隧道结很容易。然而，器件构型需要进行结构优化，而这不是本实例的主要目的。我们为此提供了优化好的器件中心区的QuantumATK Python脚本，central_region.py，您可以使用它来构建MTJ器件。如果您想学习如何进行几何弛豫，您也可以跳到 Relaxing the device central region部分。

现在，打开VNL，创建一个新项目并命名，选中它然后点击Open。下载central_region.py 并将其保存在Project Folder中。然后将脚本拖入VNL Builder中，然后使用Device from Bulk工具栏为中心区域添加电极，从而创建Fe|MgO|Fe器件。

默认的电极长度5.732 Å 对于当前的目的基本够用。这对应了铁的晶格常数的两倍。

接下来您将对MTJ左右两边的铁部分的平行自旋极化执行一个自洽NEGF-DFT计算。您也将计算透射谱与平行于输运方向的k矢量的函数关系。

将器件构型送至 Script Generator（为此使用按钮），然后为脚本添加下列模块：

• New Calculator
• Initial State
• MullikenPopulation
• Forces
• TransmissionSpectrum

然后将默认输出文件名修改为 mgo_para.nc 。

接下来，调整一些脚本模块的设置：

• Calculator:

• 选择SGGA交换关联方法，以执行一个基于PBE方法的自旋极化GGA计算。
• 将电子温度增加到1200K，以缓和自洽收敛。
• 使用一个7x7x100的k点网格。
• 为铁元素选择SingleZetaPolarized 基组。

• Initial State:

• 选择User spin作为初始态的类型。
• 设置O原子和Mg原子的相对自旋为0。

• TransmissionSpectrum:

• 为了节省时间，在本实例中您将仅计算在费米能级的输运。通过设置E0 = E1 = 0.0 eV 和仅一个能量点来完成这个计算。
• 为k点选取151×151 的网格。

保存脚本为mgo_para.py并运行计算，例如使用Job Manager。对于单核这个计算会花费一个小时，但是通过并行计算可能会得到一个很好的加速，对于4 MPI进程来说可以将计算时间缩短到大约20分钟。如果需要，您也可以下载脚本：mgo_para.py。

输出数据文件mgo_para.nc应该会在Project Files栏中显示。确保它被选中，以使数据项可以在LabFloor中显示。然后选择Forces项并使用在Panel Bar右手边的Text Representation工具查看受力的一个报告。

属于电极外延区域的原子所受的力没有被算出，所以这些被设置为0。对于其它原子，受力很小。

您也可以使用Text Representation工具去查看计算的密立根布居：

最左边栏列出了每个原子每个自旋轨道的总密立根布居。为了确保每个原子被正确地自旋极化，检查这些结果是很重要的。注意所有铁原子被极化为多数自旋向上，而氧原子基本上未被极化。紧邻界面（两端都是）的镁原子与铁原子形成反平行耦合，而远离界面的镁原子未被极化。

这个报告也显示了每个壳层以及内部的独立轨道的电荷布居。

使用Transmission Analyzer工具栏来分析计算得出的透射谱。对绘图选项进行一些修改会使得透射谱更容易分析：

• 右击图右手边的k依赖透射系数，并选择color map名为Spectral。
• 在Curves菜单中，勾选Up和Down。
• 使用Active Spin菜单，您可以选择自旋向上电子或自旋向下电子的k依赖透射。

对于计算的结果，从上图可以很清楚地看到自旋向下（少数）透射在量级上被显著抑制了，在布里渊区某些点显现出一些峰值，对应共振隧穿，而自旋向上（多数）透射显示了通常的势垒隧穿。这些都与Fe|MgO|Fe隧道结公认的隧穿机制图像相吻合（比较[BZSM01]中的Figures 6和10）。

我们现在考虑Fe|MgO|Fe隧道结左右铁部分自旋反平行排列的情况。之前计算的平行排列构型的自洽态将被用于当作反平行态的一个初始猜测。相比于完全以从头做起来开始反平行计算，这将加速SCF的收敛。回到Script Generator。如果您在设置平行排列计算之后没有关掉它，那么只需要进行几个修改。不然的话，将脚本mgo_para.py拖入Scripter然后开始对脚本进行设置，同上面设置平行排列计算完全一样。

打开Initial State模块，做如下来设置反平行自旋构型和使用平行计算作为器件密度矩阵的初始猜想：

• 再次设置所有铁原子的相对自旋为1，设置镁原子和氧原子为0。
• 将隧道结右手侧所有铁原子的初始自旋设置为-1。
• 点击Use old calculation勾选框并键入mgo_para.nc作为所载入旧的计算的文件名。

最后，改变默认输出文件名mgo_anti.nc并确保TransmissionSpectrum分析模块同自旋平行计算所使用的相同。 保存ATK Python脚本为mgo_anti.py并运行之。

注意到每个原子的受力仍然不可忽略。检查密立根布居来确保每个原子具有预期的自旋极化。

再次使用Transmission Analyzer来分析电子透射。在本例中，由于系统的镜面对称，上下（多数和少数）自旋的透射是相同的。我们再一次发现k依赖透射系数存在被零透射区域围绕的高度局域的峰。

为了理解为什么自旋向上和自旋向下的透射是相同的，首先要知道由于时间反演对称，从左到右的传播总是等同于从右到左的传播。透射谱‘上’分量对应从左电极传播到右电极的‘上’电子。透射谱‘下’分量对应从右电极传播到左电极的‘下’电子。

由于自旋反对称，左电极中的上电子和右电极中的下电子都是多数通道。由于器件的镜面对称，左多数通道到右少数通道的传播（第一个透射系数自旋分量）和右多数通道到左少数通道的传播（第二个透射系数自旋分量）是相同的。

因此，对于对称器件，两个自旋通道的相等是零偏压下反对称计算的一个重要的检查。

隧道磁电阻（TMR）的定义是

$\mathrm{TMR} = \frac{G_P-G_{AP}}{G_{AP}}$,

其中$G_P$是平行自旋排列隧道结的电导，$G_{AP}$是反平行自旋排列的电导。

电导可以通过它们各自的透射谱来计算。下面的脚本实现了这个计算。您可以在这里下载：tmr.py。

# Calculate conductance for parallel spin
transmission_para = nlread('mgo_para.nc', TransmissionSpectrum)[0]
conductance_para_uu = transmission_para.conductance(spin=Spin.Up)
conductance_para_dd = transmission_para.conductance(spin=Spin.Down)
conductance_para = conductance_para_uu + conductance_para_dd
 
# Calculate conductance for anti-parallel spin
transmission_anti = nlread('mgo_anti.nc', TransmissionSpectrum)[0]
conductance_anti_uu = transmission_anti.conductance(spin=Spin.Up)
conductance_anti_dd = transmission_anti.conductance(spin=Spin.Down)
conductance_anti = conductance_anti_uu + conductance_anti_dd
 
print 'Conductance Parallel Spin (Siemens)'
print 'Up=%8.2e, Down=%8.2e' % (conductance_para_uu.inUnitsOf(Siemens),
                                conductance_para_dd.inUnitsOf(Siemens))
print 'Total = %8.2e' % (conductance_para.inUnitsOf(Siemens))
print
 
print 'Conductance Anti-Parallel Spin (Siemens)'
print 'Up=%8.2e, Down=%8.2e' % (conductance_anti_uu.inUnitsOf(Siemens),
                                conductance_anti_dd.inUnitsOf(Siemens))
print 'Total = %8.2e' % (conductance_anti.inUnitsOf(Siemens))
print
 
print 'TMR (optimistic)  = %8.2f percent' % \
      (100.*(conductance_para-conductance_anti)/conductance_anti)
print 'TMR (pessimistic) = %8.2f percent' % \
      (100.*(conductance_para-conductance_anti)/(conductance_para+conductance_anti))

执行脚本，找到TMR：

Conductance Parallel Spin (Siemens)
Up=4.11e-09, Down=2.79e-10
Total = 4.39e-09
 
Conductance Anti-Parallel Spin (Siemens)
Up=4.22e-11, Down=4.00e-11
Total = 8.22e-11
 
TMR (optimistic)  =  5240.12 percent
TMR (pessimistic) =    96.32 percent

AdaptiveGrid类为对于k点取样的传统的Monkhorst-Pack（MP）网格实现了可供替代的选择。正如QuantumATK参考手册条目AdaptiveGrid中所详述，自适应k点网格可被用于自动放大电子透射谱的重要特征。当k依赖透射以局部峰值为主导时（比如在Fe|MgO|Fe MTJ中），总的计算的透射会严重依赖这些峰值的解析，对此使用自适应k点网格是尤为有用的。

自适应算法本质上将布里渊区（BZ）按三角形划分并对所有三角形求积分。每个三角形以迭代的方式显著变小直到达到收敛，但当细化级别数达到 maximum_number_of_levels（默认值为20）时停止。有两种误差测定可供使用：

Absolute:当对所有三角形的BZ积分的绝对值低于所选的tolerance时收敛。 Relative:当对所有三角形的BZ积分的相对值低于所选的tolerance时收敛。

图显示了Fe|MgO|Fe MTJ中的总透射如何随平行/反平行的k依赖透射函数所使用的k点分辨率变化。将使用MP方法k点网格的结果与使用自适应k点网格的结果进行了比较。

这里考虑的MTJ中总透射的计算在不少于10³ Monkhorst–Pack k点时收敛，大约对应一个31×31网格。对于自适应k点网格，一个10_-2的relative tolerance似乎就已经很充足了。

脚本adaptive_grid.py同时使用一个稠密的MP网格和自适应网格方法计算了自旋方向平行的Fe|MgO|Fe器件的透射：

# -------------------------------------------------------------
# Load device configuration
# -------------------------------------------------------------
device_configuration = nlread('mgo_para.nc', DeviceConfiguration)[0]
 
# -------------------------------------------------------------
# Transmission Spectrum using a dense Monkhorst-Pack grid.
# -------------------------------------------------------------
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(0,0,1)*eV,
    kpoints=MonkhorstPackGrid(500,500),
    )
nlsave('adaptive_grid.nc', transmission_spectrum)
 
# -------------------------------------------------------------
# AdaptiveGrid.
# -------------------------------------------------------------
adaptive_grid = AdaptiveGrid(
    tolerance=0.01,
    error_measure=Relative)
 
# -------------------------------------------------------------
# Transmission Spectrum using the AdaptiveGrid.
# -------------------------------------------------------------
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(0,0,1)*eV,
    kpoints=adaptive_grid,
    )
nlsave('adaptive_grid.nc', transmission_spectrum)

Monkhorst-Pack网格有500×500的k点，而自适应网格被设置为10^-2的relative tolerance。下载并运行脚本。如果在一个具有4 MPI进程的现代笔记本电脑上进行并行执行则需要大约25分钟完成。

使用稠密Monkhorst–Pack网格得到的自旋向上透射谱与在Elecron transmission部分使用较稀疏网格得到的透射谱有所不同：我们可以看到透射谱中间的特征实际上是环形的，并在直接围绕Г点处为零透射。显然，需要一个很稠密的网格来解析这个事实。然而，注意到改进的k点取样仅将总k点平均自旋向上透射增加了4倍。

使用自适应网格评估的自旋向上透射谱如上图所示，黑点表明了k点取样。自适应算法明显地放大了k依赖透射的环形特征，导致环形透射特征周围的k点密度最大。

在密集取样的自旋向下透射中的4个峰值在对其中之一放大之前很难被观测。将使用MP网格（上图）得到的透射与使用自适应网格（下图）得到的透射进行比较，显然，自适应算法放大了透射峰值。然而，总k点平均自旋向上透射仅改变了5%。

自适应网格类也可被用来解析k依赖透射的一个特定的部分。下面所示的脚本（zoom.py）给出了一个例子，其中自适应网格被限制在特定的$k_A,k_B$范围内，并对自旋向上和自旋向下k依赖透射的特定特征分别进行放大。使用10^-3的relative tolerance来完全解析透射特征。

得出的自旋向下和自旋向下透射阐明如下。注意图片坐标轴的不同刻度。

# -------------------------------------------------------------
# Load device configuration
# -------------------------------------------------------------
device_configuration = nlread('mgo_para.nc', DeviceConfiguration)[0]
 
# -------------------------------------------------------------
# Transmission Spectrum zoom #1.
# -------------------------------------------------------------
adaptive_grid = AdaptiveGrid(
    kA_range=[-0.04, 0.04],
    kB_range=[-0.04, 0.04],
    tolerance=0.001,
    error_measure=Relative)
 
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(0,0,1)*eV,
    kpoints=adaptive_grid,
    )
nlsave('zoom.nc', transmission_spectrum)
 
# -------------------------------------------------------------
# Transmission Spectrum zoom #2.
# -------------------------------------------------------------
adaptive_grid = AdaptiveGrid(
    kA_range=[0.275, 0.285],
    kB_range=[-0.03, 0.03],
    tolerance=0.001,
    error_measure=Relative)
 
transmission_spectrum = TransmissionSpectrum(
    configuration=device_configuration,
    energies=numpy.linspace(0,0,1)*eV,
    kpoints=adaptive_grid,
    )
nlsave('zoom.nc', transmission_spectrum)

您将在这个部分计算自旋转移矩（STT），其中需要非共线自旋。实例Introduction to noncollinear spin给出关于如何执行STT计算和非共线自旋的更多信息。 接下来，您将：

1. 计算非线性器件构型；
2. 计算密立根布居和自旋转移矩。

脚本mgo_nonco_dc.py通过使用共线铁磁计算作为起点和一个密度混合方案（在混合之前对角化密度矩阵）计算了器件构型的非共线基态。请注意这个脚本将自旋极化交换关联势（SGGA.PBE）用非共线的（NCGGA.PBE）来替代。

# Read in the collinear calculation
device_configuration = nlread('mgo_para.nc', DeviceConfiguration)[0]
 
# Use the special noncollinear mixing scheme
iteration_control_parameters = IterationControlParameters(
    algorithm=PulayMixer(noncollinear_mixing=True)
    )
# Get the calculator 
calculator = device_configuration.calculator()
new_calculator = calculator(
       exchange_correlation=NCGGA.PBE,
       iteration_control_parameters = iteration_control_parameters
       )
# Define the spin rotation
theta = 120*Degrees
left_spins = [(i, 1, 0*Degrees, 0*Degrees) for i in range(3)]
center_spins = [(i+3, 1, theta*i/5, 0*Degrees) for i in range(6)]
right_spins = [(i+9, 1, theta, 0*Degrees) for i in range(3)]
spin_list = left_spins+center_spins+right_spins
initial_spin = InitialSpin(scaled_spins=spin_list)
 
# Setup the initial state 
device_configuration.setCalculator(
    calculator=new_calculator, 
    initial_spin=initial_spin,
    initial_state=device_configuration)
 
# Calculate and save
device_configuration.update()
nlsave("mgo_nonco_dc.nc", device_configuration)

对应左电极中原子自旋极化的自旋设置指向沿着输运Z轴方向，而在右电极中极化旋转了120度。在中心区域，角度是这两个值的内插值。注意到这只是初始自旋构型-实际的自旋极化矢量将被自洽地计算并因此可能会改变（您可以在下一个部分看到结果）。

当计算结束时，您会在LabFloor中找到非共线器件构型。为了计算密立根布居和自旋转移矩，去往 Script Generator并为脚本添加如下模块：

• Analysis from File
• MullikenPopulation
• SpinTransferTorque

然后，将默认输出文件名改为mgo_nonco_dc_analysis.nc。

调整一些模块的设置：

• Analysis From File:

• 选择mgo_nonco_dc.nc 文件。

• SpinTransferTorque:

• Contributions选择Left。这样，您将计算从左到右的线性响应电流。

使用发送按钮（ 图标）将脚本送往editor并按照如下所示将Monkhorst-Pack网格用Adaptive网格来替代。

# Setup adaptive grid object.
adaptive_grid = AdaptiveGrid(
        kA_range=[-0.5, 0.5],
        kB_range=[-0.5, 0.5],
        tolerance=1e-2,
        error_measure=Relative,
        )
# -------------------------------------------------------------
# Spin Transfer Torque
# -------------------------------------------------------------
spin_transfer_torque = SpinTransferTorque(
    configuration=device_configuration,
    energy=0*eV,
    kpoints=adaptive_grid,
    contributions=Left,
    energy_zero_parameter=AverageFermiLevel,
    infinitesimal=1e-06*eV,
    self_energy_calculator=RecursionSelfEnergy(),
    )
nlsave('mgo_nonco_dc_analysis.nc', spin_transfer_torque)

如果您愿意，您可以在这里下载脚本（mgo_nonco_dc_analysis.py）。 STT的计算应该使用足够稠密的k点取样以使器件右手侧的总扭矩得到一个收敛很好的值。在这个脚本中我们使用relative tolerance为10^-2的自适应网格。下图显示了不同relative tolerance下Fe|MgO|Fe MTJ中右侧总自旋转移矩的转移情况。

这个图片显示了使用relative tolerance为10^-2对于本实例中考虑的MTJ达到计算收敛看起来是足够的。我们也使用了一个更为稠密的Monkhorst-Pack网格计算了右侧的总自旋转移矩。

可以看出使用一个Monkhorst-Pack网格，总STT不会达到完全收敛。

当计算完成时，您将在LabFloor中得到MullikenPopulation和SpinTransferTorque项目。您可以使用text representation来检查密立根电荷和自旋极化矢量分量。

在非共线计算中，‘上’和‘下’布居的总和对应着通常的密立根电荷（电子数量），它们的差异-结合两个角度-形成了自旋极化矢量。这个矢量可以使用Viewer来可视化：

您也可以使用1D Projector来绘出STT的空间分量：

下图显示了MTJ构型和STT的“y”分量。

这部分为您显示如何创建和几何优化Fe|MgO|Fe隧道结的中心区域。实例Advanced device relaxation给出更多关于整体器件优化（包括优化的电极和最佳中心区域长度）的细节说明。为了简洁，我们只执行器件中心区的受力最小化。

首先，打开Builder。为了给Fe|MgO|Fe器件设置初始几何构型，您将使用MTJ Builder工具栏，它被特定地设计用以构建这类几何构型。这样找到它：Add ‣ From Plugin ‣ Magnetic Tunnel Junction (FeMgO-style)。

MTJ builder可以构建多种不同磁性隧道结几何构型。对于当前的目的，将barrier中number of layers改为6，将number of left surface layers和number of right surface layers改为2，保持其他参数为默认值不变。点击Preview按钮查看几何构型，然后点击Build将其加入Stash。

通过点击插件工具栏上的 Bulk工具将器件几何构型转换为一个块体构型。这只是移除电极，为您留下块体中心区域。

将器件中心区域送到 Script Generator并向脚本添加如下模块：

• New Calculator
• Initial State
• OptimizeGeometry
• Forces

同时，改变默认输出文件名为mgo_relax.nc。

Calculator和InitialState设置应与在Parallel spin部分所使用的设置很相似：

• Calculator:

• 选择SGGA交换关联方法。
• 增加电子温度到1200K。
• 设置k点网格为7x7x2。这个选择给出一个在x和y方向合理的k点取样。在z方向，您想要构建一个器件构型，所以在边缘的铁原子必须是块材类似的。在z方向选择2个k点使铁边缘原子的描述更像块体。
• 普通列表项目为铁元素选择SingleZetaPolarized 基组。

• Initial State:

• 选择User spin作为初始态的类型。
• 设置氧原子和镁原子的相对自旋为0。

接下来，打开 Optimization模块来指定中心区域（电极扩展处的原子）的前四个和最后四个原子应在优化过程中保持固定。这很重要-如果这些原子中的任何一个移动了，中心胞边缘将不再是周期结构，而在之后的产生和附着器件电极时需要中心胞边缘是周期结构的。 点击Add Constraints来打开Constraints Editor。之后进行如下操作：

• 通过按住Ctrl键并在原子周围点住鼠标画矩形来选择下图所示的原子。然后点击Add tag from Selection来为这些原子分配标签为“Selection 0”。
• 对右手边电极扩展的原子做同样的操作，但使得约束类型为“Fixed”。

将脚本保存为 mgo_relax.py 并使用Job Manager或者如下命令语句来运行：

atkpython mgo_relax.py > mgo_relax.log

如果需要，您也可以在此下载脚本：mgo_relax.py。这个计算在标准个人电脑上会耗时1-2个小时，但如果分布于4个并行MPI进程下会在20分钟内完成。

在LabFloor上找到计算数据项。选择Forces项（gID002），点击Text Representation插件。您将看到所有原子最终受力的报告。

注意所有x-y方向的力为了实用目的都为零。对于已经被优化的原子，它们的受力小于优化标准0.05 eV/Ang。对于电极延伸原子，受力略微大一些。而且，力矢量指向晶胞外，表明晶胞受压应力。增加更多表面层来增加晶胞z方向长度会降低应力，但在本例中对结果影响不大。

您现在可以使用优化好的器件中心区域来进行器件输运性质的计算，参考Getting started部分。

[BZSM01] (1, 2) W. H. Butler, X.-G. Zhang, T. C. Schulthess, and J. M. MacLaren. Spin-dependent tunneling conductance of Fe|MgO|Fe sandwiches. Physical Review B, 63(5):054416, 2001.doi:10.1103/PhysRevB.63.054416.

本文翻译：王吉章