In this example, we will discuss the ATK set-up for a particular two-probe system, namely a dithiol-benzene (DTB) ring coupled to two gold (Au) fcc (111) electrodes. This system was studied in detail in [13].

Here, we will pay specific attention on how the basic elements of the two-probe system, i.e. the two electrodes and the scattering region, are constructed.

The main purpose of this tutorial is to provide you with a detailed understanding about the building blocks used for defining a two-probe set-up. Setting up a two-probe system is often a non-trivial business, and the easiest way is often to use the graphical user interface of Virtual NanoLab for this task.

However, if you wish to either use ATK directly for defining two-probe systems, or want to understand the steps taken by VNL, the relevant information is provided in this tutorial.

After you have finished the tutorial, you should be able to understand

This tutorial requires no prior knowledge about two-probe systems. You should, however, be familiar with basic Python and ATK usage corresponding to the concepts and ideas discussed in the tutorials

We will also make heavy use of the array slicing features available in the NumPy module for setting up and manipulating data structures. If these concepts are new to you, please consult this page before you proceed any further.

The main focus here will be on the ATK techniques used for defining the atomic constituents and the geometry of a two-probe system. For details regarding advanced two-probe calculations, we suggest that you read the tutorial Li-H₂-Li two-probe system.

A two-probe system consists of three units: two electrodes and a scattering region positioned in between the two electrodes:

The two-probe system that we will discuss here consists of two identical Au crystals with a dithiol-benzene molecule positioned in between the surfaces of the two Au crystals. The conformation of the system is shown in Figure 47.

Before we discuss how the entire two-probe configuration is defined, we will start by describing how its individual building blocks, namely the Au electrodes, the Au (111) layers, and the dithiol-benzene molecule, are constructed in ATK. We start out by describing the Au crystal structure in the next section.

The atoms in an Au crystal are ordered according to a face-centered cubic (fcc) Bravais lattice, i.e. a structure with eight Au atoms positioned at the lattice corners and six atoms placed on each of the six lattice faces. For gold, the lattice constant a, which corresponds to the length of any edge on the fcc lattice, is equal to a = 4.08 Å. A typical fcc lattice is shown in Figure 48. The process, where we translate the fcc lattice along the x, y, and z axis by integer multiples of the lattice constant a, will generate the entire fcc crystal. The fcc crystal lattice obtained by this process is illustrated in Figure 49.

Even though it may look as if the corner and face lattice points are not equivalent, they in fact are. This may be seen by cutting the Bravais cell along the plane q spanned by the three points (a,0,0), (0,a,0), and (0,0,a). Figure 50 shows an fcc lattice cell sliced by two such planes. The plane that we cut along is called a (111) plane, where the term (h k l) are the so-called Miller indices.

The structure of the crystal surface that appears when a large crystal (generated by translations of the fcc lattice) is being cut along a (111) plane is shown in Figure 51.

As you can see, the lattice points in a (111) plane are ordered in a two-dimensional tridiagonal lattice, where each lattice point is surrounded by a hexagonal arrangement of six nearest neighbors separated by a distance d equal to

Try for yourself to identify d as being equal to one half the length of a face diagonal of the fcc lattice. We may, in other words, also regard an fcc crystal as a stacking of (111) layers. To understand the exact stacking order of the layers, we can try to look at the geometrical conformation of a sequence of such layers. This is illustrated in Figure 52, where the arrangement of a sequence of three successive (111) layers is shown. If we label the layers from top to bottom as A, B, and C, you can see that the B layer is positioned underneath the “holes” that appear in the center of the tridiagonal patches of the A layer. However, no matter how we position the B layer, every second hole of the A layer will be empty. As you can see on the top view of the three-layer structure shown in Figure 52, the holes in the A layer that are left empty by the B layer are exactly covered by the third layer, the C layer. The layer underneath the C layer will be an A layer resulting in a ...ABCABC... packing of the fcc crystal. Crystals with layers packed according to this ABC pattern are also referred to as cubic closed-packed (ccp) structures.

The last piece of information we need in order to finish our description of the (111) crystal, is to determine the layer separation d_sep between the successive (111) layers. As shown in Figure 53, we see that d_sep can be identified as the height of the tetrahedron spanned by

Since the length of the vertices of the tetrahedron are equal to the Au-Au bond length d, the height of the tetrahedron and thereby the layer separation d_sep may be determined as

Now that we understand the finer details of a (111) Au crystal, we are ready to use ATK for defining the Au surfaces that constitutes the layers of the scattering region of the two-probe system. The first step is to load the ATK.TwoProbe and the NumPy modules into ATK and define some useful constants related to the Au crystal

from ATK.TwoProbe import *
from numpy import *

# Define constants for the fcc lattice
a     = 4.078250
d     = a/sqrt(2)
d_sep = a/sqrt(3)

Here, the respective variables a, d, and d_sep correspond to the Au lattice constant, the inter-atomic distances between each Au-Au pair, and the layer separation.

We start by setting up a 3x3 layer identical to the one shown in Figure 54. This layer corresponds to the A layer of the crystal. The two vectors u and v defined by the following coordinates

and

can be used to generate the entire 3x3 grid by using the following nested Python for loop

# Define vector generators for 3x3 layers
u  = array([1.0,             0.0, 0.0])*d
v  = array([0.5,     sqrt(3)/2.0, 0.0])*d
du = array([0.5, 1.0/2.0/sqrt(3), 0.0])*d

# Setup the A layer
n = 3
A = reshape(zeros(3*n*n), (n*n,3))*1.0
for j in range(n):
    for i in range(n):
        A[i+n*j] = i*u + j*v

Here we have used a NumPy array for storing the location of the 9 lattice points of the 3x3 lattice. The vector du defined in the above script will be used for the following purpose: By making translations of the A layer, we may easily generate the remaining B and C layers simply by translating the A layer with the displacement vector du as in Figure 55. The coordinates of du are

The layer displacements can then be used in practice for generating the B and C layers by applying the advanced slicing capabilities of NumPy:

# Setup the B and C layers
B1 = A + 1.0*du
B2 = A + 1.0*du
C  = A + 2.0*du

Here we have defined to B layers, B1 and B2, since we later on will need to position a B layer both on the left and the right side of the DTB molecule. Notice also, that the layer arrays that we have generated so far all have their z-coordinate set equal to zero. Later, when the A, B, and C layers must placed at their respective positions in the scattering region, their z-coordinates can be set using the following NumPy slicing operation

B1[:,2] += d_sep

In this case, we shift the B layer B1 one layer separation to the right along the z axis.

The DTB (dithiol-benzene) molecule is a benzene molecule where two sulfur atoms have been substituted for hydrogen at carbon atom 1 and 4. The inter-atomic distances in the molecule are

In order to set up the DTB ring as shown in Figure 56, we will use a Python for loop and iterate through the angles 0, π/3, 2π/3,...,5π/3. The position of the C, H, and S atoms can then be specified as follows

# Setup DTB molecule
dCC  = 1.397 
dCH  = 1.084 
dCS  = 1.754
dDTB = 2*(dCC + dCS)

n = 6
DTB = reshape(zeros(2*3*n), (2*n,3))*1.0
for i in range(n):
    t = i*pi/3.0
    DTB[i] = array([0.0, sin(t), cos(t)])*dCC
    if (i%3 == 0):
        DTB[i+n] = array([0.0, sin(t), cos(t)])*(dCC + dCS)
    else:
        DTB[i+n] = array([0.0, sin(t), cos(t)])*(dCC + dCH)

# Translate left S atom to z = 0
DTB[:,2] += (dCC + dCS)

The entire molecule has been positioned in the z-y-plane of the three-dimensional space with the left-most S atom (located at θ = π) positioned at (0,0,0). Notice also how we use Pythons built-in integer modulo operator % to detect the angles where the two S atoms should be inserted.

You should also observe that the positions of the C, H, and S are stored into the DTB array according to the sequence

C-C-C-C-C-C-S-H-H-S-H-H

This insertion order becomes important later on, when the entire geometry of the scattering region has to be defined.

We now have all the ingredients, the Au layers and the DTB molecule, prepared for setting up the scattering region. Our "geometric" strategy will be to order the elements in the scattering region according to the sequence AB-DTB-BC. Here is the detailed recipe that we will follow:

Let us see how we can accomplish this set up using ATK . First we align all the elements in the sequence AB-DTB-BC at their proper positions along the z-axis by using the slicing features of the NumPy module

# Define AB-DTB-BC central region
dAuDTB = 1.7

B1[:,2]  += d_sep
DTB[:,2] += d_sep + dAuDTB
B2[:,2]  += d_sep + dAuDTB + dDTB + dAuDTB
C[:,2]   += d_sep + dAuDTB + dDTB + dAuDTB + d_sep
DTB      += 6.0*du

The last line makes a translation of the DTB molecule in the x,y-plane positioning the molecule at the second hole in the 3x3 C layer.

Now that all elements have their correct positions set, we store all the information into two list variables scattering_elements and scattering_positions; the first list is used for storing information about the atom types present in the scattering region, whereas the latter contains the positions of the respective atoms. Since the scattering region is composed of the sequence AB-DTB-BC and remembering that the atoms of DTB molecule was stored in the order

C-C-C-C-C-C-S-H-H-S-H-H

we can use ATK to define the scattering_elements and scattering_positions as

# Central region elements
scattering_elements  = 18*[Gold]
scattering_elements += 6*[Carbon] + [Sulfur] + 2*[Hydrogen] + [Sulfur] + 2*[Hydrogen]
scattering_elements += 18*[Gold]

# Central region positions
scattering_positions  = concatenate((A,B1,DTB,B2,C)).tolist()
scattering_positions *= Angstrom

The elements of scattering_elements are objects of the type ATK.PeriodicTable. Notice how we can use the unary Python operator +=> to concatenate new elements to the list. This avoids having the entire definition of scattering_elements on one single line, making the script much more readable and structured. The variable scattering_positions is defined by concatenating all molecule positions from the sequence AB-DTB-BC into one NumPy array. Using the method tolist()> from the NumPy array class, this structure is converted to a regular Python list. Finally, we multiply the list by the unit ATK.Units, which converts the entire list to a list containing ATK ATK.Units objects.

The scattering region is now completely defined and is ready to be passed on to the constructor of a TwoProbeConfiguration. Before we can do that, however, we must define the two-probe electrode. We will do that in the next section.

The Au electrodes form a 3 × 3 structure with 9 atoms in each layer. Both the left and the right electrode contains 3 layers with an ABC layer stacking. This choice of sequence will ensure that the two electrodes are chosen to match the scattering DTB region which have given the structure (AB-DTB-BC). The complete two-probe structure therefore has the geometry (ABC)AB-DTB-BC(ABC).

In this example, we will construct a super cell with 3 × 3 = 9 atoms per layer. Instead of typing in the positions of the 27 Au atoms in the electrode, we will here benefit from the fact that the super cell calculation can be performed with 3 atoms and 9 k-points in the A/B-plane, instead of 27 atoms and 1 k-point.

The results will be the same, as you can convince yourself of by first applying Bloch's theorem and then folding the super cell Brillouin zone accordingly. The 27-atom case, however, requires much more CPU time, since the calculation time scales with the number of atoms to the third power, but only linearly with the number of k-points. Moreover, when running in a parallel computing environment, ATK is designed to parallelize the k-point sampling; in this case, the 27-atom calculation would, with access to 9 parallel nodes, take the same time as a 3-atom calculation.

In addition, ATK will automatically adjust the k-point grid to account for the zone folding, so the k-point sampling specified for the electrodes will correspond to the correct sampling for the complete system, i.e ultimately the composite two-probe system. To be specific, in the present example, the final two-probe calculation will be performed with a single k-point in the A/B plane (note that the sampling in the C-direction is a completely different matter), while the actual electrode calculation will be performed with 9 points.

Let us see how we put all these things into play using ATK . First, we have to define the three super cell elements as well as their cell positions. The three Au atoms in the super cell is just the first atoms in the respective A, B, and C layers corresponding to the bottom left atoms shown in Figure 55. Since the respective B and C layers are lifted by one and two layer separations relative to the A layer, we store the three positions in a list variable electrode_positions as follows

# Electrode elements and positions
electrode_elements  = 3*[Gold]
electrode_positions = [
    [0.0,      0.0,        0.0],
    [0.5*d, d*0.5/sqrt(3), 1*d_sep],
    [1.0*d, d*1.0/sqrt(3), 2*d_sep]
    ]*Angstrom

Information about the constituents of the electrode, the three Au atoms, are stored in the list variable electrode_elements. Observe also, that we must inform ATK about the unit used for "measuring" positions; we do that by multiplying electrode_positions by ATK.Units.

After this operation, we should specify the geometry of the super cell. The information we need to provide to ATK is the three primitive vectors that span the super cell. The vectors that span the x,y plane are simply the two vectors u and v shown in Figure 54. In addition, since the super cell contains three stacked layers (ABC), it should be repeated after three layer separations. The vector specifying this direction (3 × the layer separation along the z-separation), is the third primitive vector of the super cell. The whole information is stored in the list variable electrode_unitcell using the following construction

# Electrode unit cell
electrode_unitcell = [
    [d,               0.0, 0.0],
    [0.5*d, sqrt(3)/2.0*d, 0.0],
    [0.0,             0.0, 3.0*d_sep]
    ]*Angstrom

The electrode is almost finished. The last step needed is to pass on the variables to the constructor of the ATK object PeriodicAtomConfiguration, which is used for representing our semi-infinite periodic electrode. Here is how to do it

# Electrode
electrode = PeriodicAtomConfiguration(
    electrode_unitcell,
    electrode_elements,
    electrode_positions
    )

The electrode is now stored as a PeriodicAtomConfiguration in the variable electrode. In the next section, we will use this for setting up the actual two-probe system.

In order to "glue" the scattering region together with the electrodes, we need to specify the equivalent atoms of the electrodes and the scattering region. We do that using the argument equivalent_atoms which we, when defined, will pass to the constructor of a TwoProbeConfiguration object. The variable equivalent_atoms is defined as a list of the type [(x1,y1),(x2,y2)], where x1 and y1 are integer variables describing the equivalent atoms of the left electrode and the central regions. The integers x1 and y1 refer to the atomic list positions of the equivalent left electrode and central regions atoms. Similarly, the integers x2 and y2 have the same meaning for the equivalent right electrode and scattering region atoms.

In our case, we will use the convention

Since the center atom corresponds to the 5th atom in a 3x3 layer, we should hook up the electrodes to the 5th first and the 5th last atom specified in the list variable scattering_elements. Since scattering_elements contains 48 elements, we therefore have the following

Here is how we add all this information into a TwoProbeConfiguration object

# Two-probe configuration
two_probe = TwoProbeConfiguration(
    (electrode,electrode),
    scattering_elements,
    scattering_positions,
    electrode_repetitions=[[3,3],[3,3]],
    equivalent_atoms = ([0,4],[2,43])
    )

The arguments that we supply to TwoProbeConfiguration are the following:

In addition, by adding the following lines

# Export the two-probe system to an VNL file.
vnl_file=VNLFile('au-dtb-au.vnl')
vnl_file.addToSample(two_probe, 'au-dtb-au')

the two-probe structure gets written to a VNLFile, which can be used for visualizing the two-probe setup in the Virtual NanoLab.

That finishes our discussion of how to set up and assemble the building blocks of a two-probe system. The entire script for setting up the Au-DTB system is shown below. A VNL visualization of the system is shown in Figure 57.

from ATK.TwoProbe import *
from numpy import *

# Define constants for the fcc lattice
a     = 4.078250
d     = a/sqrt(2)
d_sep = a/sqrt(3)

# Define vector generators for 3x3 layers
u  = array([1.0,             0.0, 0.0])*d
v  = array([0.5,     sqrt(3)/2.0, 0.0])*d
du = array([0.5, 1.0/2.0/sqrt(3), 0.0])*d

# Setup the A layer
n = 3
A = reshape(zeros(3*n*n), (n*n,3))*1.0
for j in range(n):
    for i in range(n):
        A[i+n*j] = i*u + j*v

# Setup the B and C layers
B1 = A + 1.0*du
B2 = A + 1.0*du
C  = A + 2.0*du

# Setup DTB molecule
dCC  = 1.397 
dCH  = 1.084 
dCS  = 1.754
dDTB = 2*(dCC + dCS)

n = 6
DTB = reshape(zeros(2*3*n), (2*n,3))*1.0
for i in range(n):
    t = i*pi/3.0
    DTB[i] = array([0.0, sin(t), cos(t)])*dCC
    if (i%3 == 0):
        DTB[i+n] = array([0.0, sin(t), cos(t)])*(dCC + dCS)
    else:
        DTB[i+n] = array([0.0, sin(t), cos(t)])*(dCC + dCH)

# Translate left S atom to z = 0
DTB[:,2] += (dCC + dCS)

# Define AB-DTB-BC central region
dAuDTB = 1.7

B1[:,2]  += d_sep
DTB[:,2] += d_sep + dAuDTB
B2[:,2]  += d_sep + dAuDTB + dDTB + dAuDTB
C[:,2]   += d_sep + dAuDTB + dDTB + dAuDTB + d_sep
DTB      += 6.0*du

# Central region elements
scattering_elements  = 18*[Gold]
scattering_elements += 6*[Carbon] + [Sulfur] + 2*[Hydrogen] + [Sulfur] + 2*[Hydrogen]
scattering_elements += 18*[Gold]

# Central region positions
scattering_positions  = concatenate((A,B1,DTB,B2,C)).tolist()
scattering_positions *= Angstrom

# Electrode elements and positions
electrode_elements  = 3*[Gold]
electrode_positions = [
    [0.0,      0.0,        0.0],
    [0.5*d, d*0.5/sqrt(3), 1*d_sep],
    [1.0*d, d*1.0/sqrt(3), 2*d_sep]
    ]*Angstrom

# Electrode unit cell
electrode_unitcell = [
    [d,               0.0, 0.0],
    [0.5*d, sqrt(3)/2.0*d, 0.0],
    [0.0,             0.0, 3.0*d_sep]
    ]*Angstrom

# Electrode
electrode = PeriodicAtomConfiguration(
    electrode_unitcell,
    electrode_elements,
    electrode_positions
    )

# Two-probe configuration
two_probe = TwoProbeConfiguration(
    (electrode,electrode),
    scattering_elements,
    scattering_positions,
    electrode_repetitions=[[3,3],[3,3]],
    equivalent_atoms = ([0,4],[2,43])
    )

# Export the two-probe system to an VNL file.
vnl_file=VNLFile('au-dtb-au.vnl')
vnl_file.addToSample(two_probe, 'au-dtb-au')

au-dtb.py