目录

双电极体系中如何正确设置电极

简介

双电极体系的电极建模看似简单，但实际上也可能对计算有重要的影响。清楚的了解与电极有关的基本概念，有利于计算得到更可靠的结果。电极建模的最关键的一点是电极的长度。下面将会介绍，电极长度太短会可能会导致结果不准确或计算无法收敛。

电极单胞

电极单胞需要符合以下规则：

ATK通常能够自动判断单胞是否满足这两个条件。

周期性与真空

电极单胞还须满足以下规则：

例如，创建fcc格子的[111]方向晶胞时，我们会得到Z方向有三层A-B-C重复单元。

整个电极体系在Z方向可以是非周期性的。这有两点含义，一方面，体系在有限偏压下，即使完整周期结构也变成了不对称的非周期结构；另一方面，左、右电极可以是不同种材料。

电极体系在X、Y方向则必须是周期性的，对使用FFT2D方法求解泊松方程是如此，对MultiGrid方法求解泊松方程时仍然如此。

研究一维、二维体系时，用户应该在X、Y方向留下足够的真空填充，以避免重复单元之间的基组重叠和静电相互作用。

研究三维表面的分子结时，同样也需要考虑同样的问题，电极表面XY方向的周期应足够大，以保证重复单元中的分子间距大于应有的相互作用距离，这样我们研究的体系才是孤立分子的导电性，而没有引入分子间相互作用。

XY方向上的周期性在模拟层状结构、或层状堆叠的界面结构时非常有用，因为我们可以用一个比较小的单胞研究这类体系。

Z方向电极长度的影响

上文提到对于fcc[111]方向的电极单胞为ABC三层原子，多数文章中用到的电极都是如此。

但是如果考虑更简单的体系，比如实例教程中的Li-H2-Li体系，似乎可以采用“单个原子”作为电极周期：

然而，用这样的模型进行计算，可能会遇到不收敛的问题，或者即使“运气好”收敛了，也很难保证结果的可靠性。这种在建模时就埋下的隐患极有可能会导致计算失败，而且很难察觉。

导致这个模型失败的原因是：由于电极太短，电极的费米能级与中间区域的费米能级无法匹配。为了避免处理无限大体系，计算时先分别处理两个有限的周期电极体系，然后再将电极与同样是非周期有限大小的中间区域耦合。

因此，算法要求电极单胞中的原子仅与Z方向最近邻晶胞中的原子相互作用。这样的假设忽略了除此以外的相互作用。这使得双电极体系与周期电极体系的哈密顿量无法匹配，因为后者在计算时考虑了所有的相互作用。结果导致电极的费米能级与双电极体系的平衡态费米能级无法匹配。

由于采用SIESTA有限基组（作用范围有限），以上要求在计算有限大小晶胞时可以严格遵守。用高斯基组或平面波基组是无法做到这一点的。

上面提到的相互作用有两类：

直接相互作用很重要，必须包括，否则结果会受到明显的影响。二阶的相互作用项则在很多时候对结果的影响微乎其微，可以忽略。唯独对于态密度在费米面处剧烈变化的体系，可能需要对计算准确的费米能级有比较大的影响。

足以包含所有相互作用的标准可以定义为：

如果只考虑直接重叠的影响，基本上就是基组半径的二倍。（通常这个标准也足以包括二阶相互作用。）

如何确保合适的电极长度

实际操作时的问题时：如何判断电极是否长度合适？答案很明显，找出体系中最长的相互作用半径，将电极长度与之比较。

ATK-DFT中Li原子的基组是4.04Å，projector最大半径是2.01Å，相互作用半径为6.05Å。Li原子间距为2.9Å，因此，电极应该包含三个原子（周期8.7Å）就保证了大于二倍的基组作用半径。当电极包含四个原子时，尽管长度仍然不足最大相互作用半径的二倍，但相互作用也大体上都被包括进来。

快速判断电极长度的脚本

check_electrode_length.py
def ranges(device_configuration,method):
    """
    Returns the max basis set and projector radii, in Angstrom
 
    @param device_configuration The device configuration to check.
    """
    # Import required modules.
    from NL.Calculators.BulkCalculatorInterface import minimalBasis
 
    calculator = device_configuration.calculator()
 
    # If no calculator is attached, use specified method with default parameters
    if calculator is None:
        # In order not to destroy the input, we should make a local copy in that case
        print "No calculator on configuration, using %s with default parameters" % method
        if method=='DFT':
            calculator = DeviceLCAOCalculator()
        else:
            calculator = DeviceHuckelCalculator()
 
    # If it is a device configuration we will need to take into account projectors.
    if isinstance(calculator, DeviceLCAOCalculator):
        # Get the exchange-correlation calculator.
        xc_calculator = calculator.exchangeCorrelation()._radialBackEngine()
        # Get the minimal basis.
        minimal_basis = minimalBasis(device_configuration, calculator.basisSet())
 
        # Find the basis range.
        basis_range = 0.0
        for basis in minimal_basis:
            b = basis._backEngine(xc_calculator)
            N = b.size()
            for i in range(N):
                basis_range = max(basis_range, b[i].range())
 
        # Find the projector range.
        projector_range = 0.0            
        for basis in minimal_basis:
            pseudopotential = basis.pseudopotential()
            pseudopotential.load()
            projectors = pseudopotential._projectors()
            N = projectors.size()
            for i in range(N):
                projector_range = max(projector_range, projectors[i].range())
 
        return ((basis_range*Bohr).convertTo(Angstrom),(projector_range*Bohr).convertTo(Angstrom))
 
    elif isinstance(calculator, DeviceHuckelCalculator):
        # Get the builder.
        builder = calculator._builder()(device_configuration, calculator)
        matrix_calculator = builder.createMatrixCalculator()
 
        # Find the basis ranges.
        basis_range = 0.0
        for e in numpy.unique(device_configuration.elements()):
            ranges = matrix_calculator.basisSet().ranges(e.symbol())
            for i in range(ranges.size()):
                basis_range = max(basis_range, ranges[i])
 
        return ((basis_range*Bohr).convertTo(Angstrom),0.)
 
 
def minimal_nextnext_distance(electrode):
    '''
    Computes the smallest distance between the atoms in the electrode
    and the atoms in the corresponding next-next neighbor electrode cell
 
    @param electrode   A BulkConfiguration
    @return            Minimal distance, in Angstrom
    '''
    # Length of electrode C vector in Z
    electrode_Cz = electrode.bravaisLattice().primitiveVectors()[2,2].inUnitsOf(Angstrom)
    # Original electrode coordinates
    elec_coords  = electrode.cartesianCoordinates().inUnitsOf(Angstrom)
    # Coordinates of atoms in next-next cell
    test_coords  = elec_coords+numpy.array([0,0,2.*electrode_Cz])
    # Compute the distance from each atom in next-next cell to all atoms in original cell (order N)
    distances = []
    for test in test_coords:
        distances.append(numpy.sqrt(((test-elec_coords)**2).sum(axis=1)))
    # Return the smallest such distance
    return numpy.min(numpy.array(distances))
 
 
def print_status(lr,dist,basis_range,projector_range):
    if dist < 2.*basis_range:
        print 'WARNING: Truncation of direct basis pair overlaps will occur in the %s electrode in this system. Your results will be incorrect!' % lr
    elif dist < 2.0*(basis_range + projector_range):
        print 'NOTE: Truncation of some second-order interactions will occur in the %s electrode in this system. This may influence convergence and the results.' % lr
    else:
        print 'The %s electrode is sufficiently long to include all interations.' % lr
 
 
def validateElectrodes(device_configuration, method):
    """
    Method for checking that the electrodes is long enough
 
    @param device_configuration The device configuration to check
    """
 
    if device_configuration==None:
        return None
 
    # Compute the maximum basis set and projector (for DFT only) radius
    basis_range,projector_range = ranges(device_configuration,method)
 
    left_mini_dist  = minimal_nextnext_distance(device_configuration.electrodes()[0])
    right_mini_dist = minimal_nextnext_distance(device_configuration.electrodes()[1])
 
    print 'Maximum basis set radius             = %8.4f Ang' % basis_range
    print 'Maximum basis + projector radius     = %8.4f Ang' % (basis_range + projector_range)
    print 'Basis set overlap range              = %8.4f Ang' % (2.0*basis_range)
    print 'Maximum interaction range            = %8.4f Ang' % (2.0*(basis_range + projector_range))
    print 'Left electrode, truncation distance  = %8.4f Ang' % left_mini_dist
    print 'Right electrode, truncation distance = %8.4f Ang' % right_mini_dist
    print '(Truncation occurs at the smallest atom-atom distance between electrode and its next-nearest neighbor cell)\n'
 
    print_status("LEFT",left_mini_dist,basis_range,projector_range)
    print
    print_status("RIGHT",right_mini_dist,basis_range,projector_range)
 
    return device_configuration
 
 
# Initialize the Analyzer
builder = Builder()
builder.title('Electrode length check')
 
# Set the Analysis function
builder.setConfigurationGenerator(validateElectrodes)
 
# Set up a Analysis widget interface
builder.configuration('device_configuration', label='Device configuration', text='Drop system here')
builder.choice('method', [
    ('DFT', 'DFT'),
    ('Extended Huckel', 'Extended Huckel'),
    ], 'DFT', 'Default method')
 

使用方法

  1. 下载脚本
  2. 在VNL中打开Custom Builder，将脚本拖入；
  3. 将双电极体系结构模型拖放到Drop system here；
  4. 程序根据模型脚本中设置的参数进行计算，或手动设置DFT还是Hückel方法；
  5. log窗口中会显示所有重要的长度，并给出电极长度是否合适的结论。

原始文章