Usage examples

Evenly spaced grid of k-points:

bz_int_param = brillouinZoneIntegrationParameters((5,5,5))

Typical parameters for a nanotube:

bz_int_param = brillouinZoneIntegrationParameters((1,1,100))

Typical parameters for a metallic electrode:

bz_int_param = brillouinZoneIntegrationParameters((3,3,50))
method = KohnShamMethod(
    brillouin_zone_integration_parameters=bz_int_param
    )
scf = method.apply(...)

Calculate the transmission spectrum:

from ATK.TwoProbe import *

scf = executeSelfConsistentCalculation(...)
bz_int_param = brillouinZoneIntegrationParameters((3,3))

trans_spectrum = calculateTransmissionSpectrum(
       self_consistent_calculation = scf,
       energies = [0.0]*electronVolt,
       brillouin_zone_integration_parameters=bz_int_param
    )

Notes

There are no default values for the k-point sampling used in bulk or electrode systems. The user must specify the Brillouin zone sampling for every such calculation. For transmission calculations, the default is (1,1), i.e. only the Γ point is included.

The number of k-points necessary for an adequate description of the system varies. For a two-probe system, a large number of k-points should generally be used along the transport direction to obtain an accurate description of the electronic structure of the infinite electrodes. See the tutorial for further details and examples.

ATK employs time-reversal symmetry to reduce the number of unique k-points that need to be included in the Brillouin zone integration.

For bulk calculations in which the Brillouin zone is sampled using only the Γ k-point, ATK makes use of a real eigenvalue solver. A complex eigenvalue solver is instead used in the general case, in which more k-points are used to sample the Brillouin zone. Depending on the system size and basis set, a real eigenvalue solver can provide a significant speed-up for bulk calculations.

Obtaining Monkhorst-Pack grid

Having performed a bulk calculation suppose you would like to know the actual k-points used for comparison reasons or to obtain the k-points weights. Below, we provide two scripts (along with code comments) which combined can be used to obtain the grid generated by ATK for a set of Monkhorst-Pack[7] parameters and a bulk atomic configuration. The first script (MonkhorstPack.py) is the actual file in which you specify the Brillouin zone and how many k-points you want to find, whereas the second script (bulk_utilities.py) contains a set of library functions that MonkhorstPack.py needs access to. These two files should be in the same directory.

def generateMonkhorstPackGrid (n_abc, R, time_reversal_symmetry):
    import math
    from numpy import array, arange, ones

    '''
    Returns the k-points of a (na x nb x nc) Monkhorst-Pack grid
    for (reciprocal) lattice vectors R = [A,B,C]
    where na = n_abc[0] etc
    Optionally, the grid can be reduced by time-reversal symmetry,
    which removes -k where k is present, except for the Gamma point
    The general formulae for the grid are Eqs. 3-4 in
    H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976)
    which state that
        (3)  ur = (2r-N-1)/2N   where   r=1,2,...,N
    where N is the number of states in each direction.
    These integers are passed as na, nb, nc to this function
    The k-points are then
        (4)  k = up * A + uq * B + ur * C
    where A, B, C are the reciprocal lattice vectors
    '''
    
    '''
    In the case of time-reversal symmetry, there are half as many
    k-points plus the Gamma point, if included. The total number
    of points na*nb*nc is odd if the Gamma point is included,
    thus (na*nb*nc+1)/2 using integer division will give the
    correct number of points.
    '''
    n_product = n_abc[0]*n_abc[1]*n_abc[2]
    total_number_of_kpoints = n_product
    if (time_reversal_symmetry):
        total_number_of_kpoints = (n_product+1)/2

    # Create the numbers ur according to Eq. (3)
    ua = (2*arange(1,n_abc[0]+1)-n_abc[0]-1)/(2.*n_abc[0])
    ub = (2*arange(1,n_abc[1]+1)-n_abc[1]-1)/(2.*n_abc[1])
    uc = (2*arange(1,n_abc[2]+1)-n_abc[2]-1)/(2.*n_abc[2])

    if (time_reversal_symmetry):
        '''
        Impose time-reversal symmetry by removing one half of ua
        (or half of ub, if ua has only a single point, 
        or half of uc if both ua and ub have a single point, 
        or remove none if there is just a single point)
        '''
        if   len(ua)>1:
            ua = ua[0:(len(ua)+1)/2]
        elif len(ub)>1:
            ub = ub[0:(len(ub)+1)/2]
        elif len(uc)>1:
            uc = uc[0:(len(uc)+1)/2]

    # Weights
    if (time_reversal_symmetry):
        wG = 1./n_product
        w0 = 2.*wG
        # If na*nb*nc is odd, it means that the Gamma
        # point (0,0,0) is included.
        # It has half of the weight of the other points
        # if time-reversal symmetry is used
        if (n_product % 2):
            W = [w0,]*(total_number_of_kpoints-1)+[wG]
        else:
            W = [w0,]*total_number_of_kpoints
    else:
        w0 = 1./(total_number_of_kpoints)
        W = [w0,]*total_number_of_kpoints
    
    # K-points
    K = []
    for i in range(len(ua)):
        for j in range(len(ub)):
            for k in range(len(uc)):
                # Calculate k-points according to Eq. (4)
                K.append(ua[i]*R[0]+ub[j]*R[1]+uc[k]*R[2])
                '''
                When time-reversal is used, the ua vectors are ordered
                such that when we hit the Gamma point, the remaining
                points are time-reversed copies of the already
                included points, and should not be included
                '''
                if (time_reversal_symmetry) & \
                       (ua[i]==0) & (ub[j]==0) & (uc[k]==0):
                    return K,W
    return K,W

def printMonkhorstPackGrid(K,W):
    for i in range(len(K)):
        print "%.5f  \t%.5f  \t%.5f  \t%.5f" %\
              (K[i][0],K[i][1],K[i][2],W[i])

def MonkhorstPackGrid (n_abc,
                       bulk_configuration,
                       time_reversal_symmetry=True,
                       verbosity_level=0):
    from bulk_utilities import reciprocalVectors
    R = reciprocalVectors(bulk_configuration)
    (K,W) = generateMonkhorstPackGrid (n_abc, R, time_reversal_symmetry)
    if verbosity_level>=10:
        printMonkhorstPackGrid(K,W)
    return (K,W)

MonkhorstPack.py

from ATK.KohnSham import *
from numpy import *

def cross (a,b):
    return array([a[1]*b[2]-a[2]*b[1],
                  a[2]*b[0]-a[0]*b[2],
                  a[0]*b[1]-a[1]*b[0]])

def invertVectors (v):
    v1 = v[0]
    v2 = v[1]
    v3 = v[2]
    norm = dot(v1,cross(v2,v3))
    r1 = cross(v2,v3)
    r2 = cross(v3,v1)
    r3 = cross(v1,v2)
    return 2.*pi*array([r1,r2,r3])/norm

def primitiveVectors (lattice):
    V = identity(3)*1.
    for i in range(3):
        for j in range(3):
            V[i,j] = lattice.primitiveVectors()[i][j].inUnitsOf(Ang)
    return V
    
def reciprocalVectors (bulk_config):
    return invertVectors(primitiveVectors(bulk_config.bravaisLattice()))
    
def convertCartesianCoordsToLattice (bulk_config,cartesian_coordinates):
    import numpy
    from numpy.linear_algebra import inverse as matrixinverse
    V = numpy.transpose(reciprocalVectors(bulk_config))
    return numpy.matrixmultiply(
        matrixinverse(V),
        numpy.transpose(
            cartesian_coordinates
            )
        )

bulk_utilities.py

How to use these two scripts should be apparent from the example shown below. Here, we specify a bulk bcc iron crystal and use it to find 27 k-points:

import ATK
from ATK.KohnSham import *
from MonkhorstPack import *

ATK.setVerbosityLevel(10)

atom = [Iron]

position = [(0.0,0.0,0.0)]*Angstrom

unitcell = BodyCenteredCubic(2.87*Angstrom)

iron_crystal = BulkConfiguration(
    bravais_lattice=unitcell,
    elements=atom,
    cartesian_coordinates=position
    )

(k_points,weights) = MonkhorstPackGrid(
    n_abc=[3,3,3],
    bulk_configuration=iron_crystal,
    time_reversal_symmetry=True,
    verbosity_level=ATK.verbosityLevel()
    )

MonkhorstPack-example.py

This script will store the actual coordinates of the k-points in the variable k_point whereas the weight of each k-point is stored in the variable weights. In the example above we have used the names of the optional and required parameters. Required arguments to the MonkhorstPackGrid function includes a valid BulkConfiguration and a three-dimensional set of Monkhorst-Pack parameters which specifies how dense a grid to return.

Try running the example above and see for yourself how all k-points are generated. The generated k-point-grid along with weights is printed if the verbosity level is higher than or equal to 10. Here, we have used the built-in verbosity functionality of ATK to specify it but you could just as well have typed an integer number. The following output should display on your screen if you run the script above:

-1.45951        -1.45951        -1.45951        0.07407
-0.72975        -0.72975        -1.45951        0.07407
 0.00000         0.00000        -1.45951        0.07407
-0.72975        -1.45951        -0.72975        0.07407
 0.00000        -0.72975        -0.72975        0.07407
 0.72975         0.00000        -0.72975        0.07407
 0.00000        -1.45951         0.00000        0.07407
 0.72975        -0.72975         0.00000        0.07407
 1.45951         0.00000         0.00000        0.07407
-1.45951        -0.72975        -0.72975        0.07407
-0.72975         0.00000        -0.72975        0.07407
 0.00000         0.72975        -0.72975        0.07407
-0.72975        -0.72975         0.00000        0.07407
 0.00000         0.00000         0.00000        0.03704

Additional information can be found in the code comments and in the original article [7].

Further examples