Our first step is to set up the Bravais lattice associated with the silicon crystal. The following script shows how this is done using NanoLanguage:

from ATK.KohnSham import *

# Create Bravais lattice 
lattice_constant = 5.43*Angstrom
fcc_lattice = FaceCenteredCubic(lattice_constant)

The first line (after the import statement) defines a keyword a, which is set to the lattice constant of silicon. After that, we define the lattice for the silicon crystal as a FaceCenteredCubic object, which we then assign to the variable fcc_lattice. NanoLanguage has native support for all 14 Bravais lattices, so we never have to to specify a Bravais lattice in terms of its unit cell vectors.

So far, we have just created an empty lattice without any atoms. The atoms of the crystal are specified using a BulkConfiguration object in the following way:

# Specify bulk configuration with 2 silicon atoms
# as basis in a FCC Bravais lattice
si_fcc = BulkConfiguration(
    bravais_lattice = fcc_lattice,
    elements = 2*[Silicon],
    fractional_coordinates = [
        [0.00, 0.00, 0.00], 
        [0.25, 0.25, 0.25]
        ]
    )

The arguments to BulkConfiguration are:

There are two NanoLanguage methods for specifying atomic positions, namely

We generally recommend the latter approach, as relative coordinates are independent of lattice constants, but cases exists where Cartesian coordinates will be easier to use. See the reference documentation for BulkConfiguration for more specific details.

Now that we have defined the bulk crystal, it is a good idea to make sure that it has been set up properly. For convenience, let us define a general Python function that takes one single argument, namely the bulk configuration, and prints atomic symbols for all atoms, the corresponding positions, as well as the unit cell vectors. Below is a Python function printBulkConfig() that takes care of the task

from ATK.KohnSham import *

TAB = "\t"

def printAtomsInBulk(bulkSample):
    for (atom,pos) in zip(bulkSample.elements(), bulkSample.cartesianCoordinates()): 
        print atom.symbol(), TAB,
        for i in range( len(pos) ):
            print pos[i].inUnitsOf(Ang), TAB,
        print

def printLatticeVectors(lattice):
    for vector in lattice.primitiveVectors():
        for i in range( len(vector) ):
            print vector[i].inUnitsOf(Ang), TAB,
        print

def printBulkConfig(bulkSample):
    print "\nPrimitive vectors (in Angstrom)"
    printLatticeVectors(bulkSample.bravaisLattice())
    print "\nAtom list (positions in Angstrom)"
    printAtomsInBulk(bulkSample)
    print

printBulkConfig.py

Download and save the above script under the name printBulkConfig.py. The function printBulkConfig() can then be imported and used subsequently by adding the two lines

# Display Cartesian coordinates for silicon crystal
from printBulkConfig import printBulkConfig
printBulkConfig(si_fcc)

to the main script.

Finally, we store the silicon crystal sample in a VNLFile, so we can reuse it later in a different script. This is accomplished using the following two lines:

# Store crystal structure in a VNL file for later use
vnl_file=VNLFile("si.vnl")
vnl_file.addToSample(si_fcc,"si_fcc")

The entire script for setting up the fcc silicon crystal now reads

from ATK.KohnSham import *

# Create Bravais lattice 
lattice_constant = 5.43*Angstrom
fcc_lattice = FaceCenteredCubic(lattice_constant)

# Specify bulk configuration with 2 silicon atoms
# as basis in a FCC Bravais lattice
si_fcc = BulkConfiguration(
    bravais_lattice = fcc_lattice,
    elements = 2*[Silicon],
    fractional_coordinates = [
        [0.00, 0.00, 0.00], 
        [0.25, 0.25, 0.25]
        ]
    )

# Display Cartesian coordinates for silicon crystal
from printBulkConfig import printBulkConfig
printBulkConfig(si_fcc)

# Store crystal structure in a VNL file for later use
vnl_file=VNLFile("si.vnl")
vnl_file.addToSample(si_fcc,"si_fcc")

bulk-setup.py

Download and save the bulk-setup.py script to your disk, and process it with ATK by issuing the command

atk bulk-setup.py

This will print the following lines of output:

Primitive vectors (in Angstrom)
0.0     2.715   2.715   
2.715   0.0     2.715   
2.715   2.715   0.0     

Atom list (positions in Angstrom)
Si      0.0     0.0     0.0     
Si      1.3575  1.3575  1.3575

and generate the VNLFile si.vnl, which will now be located on your disk. The file si.vnl can be imported into Virtual NanoLab if you would like to study a 3D rendering of the silicon crystal. A 3D rendering of the silicon structure is in Figure 18.

The fact that the periodicity is an inherent feature of a crystal structure, has very important consequences for the numerical treatment of such structures; these must be considered when we calculate periodic systems in ATK.

The fact that a crystal is invariant in certain directions with respect to translations by a unit cell length implies that the wave vector is a conserved quantity (note that, you can have 1D and 2D periodic systems also). As a result, the eigenfunctions obey Bloch's theorem and are therefore labeled by the wave vector k (in 1D, k is just a scalar, the wave number k). This also means that the Schrödinger equation can be solved independently for each possible value of k.

When constructing the total electron density , it is, however, necessary to know the solution in all points of the first Brillouin zone, since

where is the occupation of the state k in the 'th band and is the corresponding Bloch function.

In practice, the integral is approximated by a sum, which is then evaluated over a finite number of k-point samples in the Brillouin zone. Clearly, as the number of k-points increases, the sum approximates the 'true' value of the above integral to a greater accuracy.

In the present version of ATK, the k-point sampling is specified in terms of how many points that should be used in the three directions of the unit cell vectors. These directions need not always coincide with the Cartesian axes.

The k-point sampling is set up by calling the function brillouinZoneIntegrationParameters() as

brillouinZoneIntegrationParameters([n1,n2,n3])

with a list argument containing three positive, non-zero integers, where each integer corresponds to the number of k-points for the respective axis directions.

This sets up a mesh with k-points uniformly spaced in the first Brillouin zone (i.e. in the reciprocal lattice) following the Monkhorst-Pack scheme [7].

Exactly how many k-points to choose depends to a large extent on the geometry and dimensionality of the system. For example, the bigger the unit cell is, the fewer points are generally needed. Please note, that it is only necessary to sample the directions in which the crystal is periodic. A carbon nanotube, for example, only requires a large k-point sampling along the tube axis (and a single point, of course, in the two transverse directions).

Let us use the script below to check how many k-points are needed in order to obtain an acceptable value of the total energy in a silicon crystal. By acceptable we mean within a certain degree of accuracy.

from ATK.KohnSham import *

# Restore previously created silicon structure
si_file = VNLFile("si.vnl")
si_fcc  = si_file.readAtomicConfigurations()["si_fcc"]

# Define basis set and mesh cut-off
basis = basisSetParameters(type=SingleZeta)
mesh  = electronDensityParameters(mesh_cutoff=75.0*Rydberg)

print '\nConverging bulk fcc Si in number of k-points'
print 'kpts\t\tTotal energy (eV)\t(rel. diff)'
print '---------------------------------------------'

convergence_criterion = 1.0e-5

kpts = 0
E = Eprev = 1.0*eV

# Loop determine number of k-points for convergence
while ( (kpts<1) or (abs( (Eprev-E)/E ) > convergence_criterion) ):
    # Increase number of k-points for each iteration
    kpts = kpts + 1
    Eprev = E
    kgrid = brillouinZoneIntegrationParameters(
        monkhorst_pack_parameters=3*[kpts]
        )
    dft_method = KohnShamMethod(
        basis_set_parameters=basis,
        electron_density_parameters=mesh,
        brillouin_zone_integration_parameters=kgrid
        )
    # Perform actual calculation
    if kpts > 1:
        # Use the previous scf as initial calculation
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method,
            initial_calculation=scf,
            )
    else:
        # For (1,1,1) there is no initial calculation to use
        scf = executeSelfConsistentCalculation(
            atomic_configuration=si_fcc,
            method=dft_method
            )
    E = calculateTotalEnergy(scf)
    if kpts > 1:
        print 3*(kpts,), '\t', E.inUnitsOf(eV), '\t\t', abs((Eprev-E)/E)
    else:
        print 3*(kpts,), '\t', E.inUnitsOf(eV)
        
print '\nConverged with',kpts,'k-points in each direction\n'

bulk-kpoint.py

Before we actually run the script, a few remarks are in place regarding its structure

Observe, that the above script explicitly uses a single zeta basis set. In this case, this simple basis set is sufficient for determining the optimal k-point sampling.

Now, download and save the script, and process it through ATK. You will get the following output

Converging bulk fcc Si in number of k-points
kpts            Total energy (eV)       (rel. diff)
---------------------------------------------
(1, 1, 1)       -337.095933988
(2, 2, 2)       -356.154925831          0.0535132058008
(3, 3, 3)       -355.723973423          0.00121147979962
(4, 4, 4)       -356.327501856          0.00169374642619
(5, 5, 5)       -356.273593339          0.000151312132773
(6, 6, 6)       -356.329983582          0.000158252871966
(7, 7, 7)       -356.32241598           2.1238075445e-05
(8, 8, 8)       -356.330060267          2.14528268082e-05
(9, 9, 9)       -356.328781041          3.59001521208e-06

Converged with 9 k-points in each direction

From the results, we see, that a (9,9,9)-grid gives excellent convergence for the silicon crystal. Inspecting the numbers a bit closer reveals that the accuracy is acceptable even for a (4,4,4)-grid. It is noteworthy how small the change in the total energy is when going from a (4,4,4) to a (5,5,5) grid, or more generally from an even number in each direction to an odd number. Clearly the (5,5,5) calculation takes longer time, but does not improve the accuracy much. This is a typical feature when using a small number of k-points in all directions.

Reusing the old density can be used to save considerable time. Say, you have an old calculation which was converged using a given set of k-points but now you want to use an even larger mesh to obtain a better accuracy. Using the initial_calculation keyword this can be done in a few iterations instead of several steps if ATK had to perform the calculation from scratch.

The first step is to perform the self-consistent calculation of the silicon crystal. We know that a k-point mesh of (4,4,4) is sufficient to obtain converged total energies; in the tutorial Calculating molecular properties, we examined the effect of choosing a different basis set and exchange-correlation functional. For improved accuracy, we will here use the GGA-PBE- approximation for the exchange-correlation functional, as well as the DoubleZetaPolarized basis set. The following script starts the self-consistent field calculation and stores the results:

from ATK.KohnSham import *
import ATK

# Restore previously defined Silicon structure
si_fcc = VNLFile("si.vnl").readAtomicConfigurations()["si_fcc"]

#Specify NetCDF file for storing results
ATK.setCheckpointFilename('si_gga_dzp.nc')

# Define new KohnShamMethod using a non-default
# exchange-correlation potential (GGA)
method = KohnShamMethod(
    exchange_correlation_type = GGA.PBE,
    basis_set_parameters = basisSetParameters(type=AtomicOrbitals.DZP),
    brillouin_zone_integration_parameters = \
        brillouinZoneIntegrationParameters(monkhorst_pack_parameters = [4,4,4])
    )

# Perform the actual calculation
gga_calculation = method.apply(atomic_configuration = si_fcc)

bulk-gga_calculation.py

The results of the calculation are stored in the NetCDF file si_gga_dzp.nc. Note, that the checkpoint file is overwritten without warning if it already exists! More information can be found in the reference manual under setCheckpointFilename().

After the calculation has finished, we can use the file si_gga_dzp.nc as the starting point for calculating the band structure.

One of the analysis functions for studying bulk samples is the calculateEigenstates() function. By specifying the quantum number for specific states at a desired k-vector, it is possible to visualize the eigenstates of silicon using VNL.

If we restore our DFT calculation from above (stored in the si_gga_dzp.nc file), the eigenstates can be calculated using a minimum of effort.

from ATK.KohnSham import *

# Restore old NetCDF results
gga_calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')

# Calculate eigenstates
eigenstates = calculateEigenstates(
    self_consistent_calculation = gga_calculation,
    quantum_numbers = ( (0,1), (0.375,0.375,0.75) )
)

# Store eigenstates in a VNL file
vnl_file=VNLFile("si.vnl")
for eigenstate in eigenstates:
    vnl_file.addToSample(eigenstate, "si_fcc")

bulk-eigenstates.py

The quantum numbers given to the VNLFile function are given as a tuple sequence of

As we see in the example, the band index can be a sequence in itself, to get the Bloch functions of several bands at the same k-point. Comparing with the list of symmetry points below, we see that the selected point is the 'K' point. The eigenstates are stored in the same VNLFile using the same sample name si_fcc as used above to associate the eigenstates with the silicon structure. If there are several band indices, the returned object eigenstates is a list of Bloch functions, which is why we loop over eigenstate when saving the plots to the VNL file.

The generated VNLFile can now be imported into Virtual NanoLab, which then can be used to generate a 3D rendering of the eigenstates similar to the plot shown in Figure 20.

Our next task is to calculate the band structure of silicon along a specific direction in the Brillouin zone.

Because the complete script for calculating and plotting the band structure of the silicon crystal is somewhat lengthy, we will go through the key elements of the script step-by-step. We start by creating a dictionary containing common symmetry points in the Brillouin zone for the fcc lattice

# Define common symmetry points in the Brillouin zone of silicon
fcc_symmetry_points = {
    'G' : (0.000, 0.000, 0.000),
    'X' : (0.500, 0.000, 0.500), 
    'W' : (0.500, 0.250, 0.750), 
    'L' : (0.500, 0.500, 0.500), 
    'K' : (0.375, 0.375, 0.750), 
    'U' : (0.625, 0.250, 0.625) 
}

The points are given in lattice coordinates independent of the lattice constant. ATK will automatically convert these points into Cartesian coordinates, using the reciprocal lattice vectors of the fcc lattice. We have used the label 'G' to denote the Γ-point of the Brillouin zone.

Now, we want to create a "path" of points in the Brillouin zone connecting two points in the above dictionary. For simplicity, we only consider the path between the Γ-point and and the point 'X'.

#Create line in Brillouin zone along which the
#band structure is calculated
start_point = fcc_symmetry_points['G']
end_point = fcc_symmetry_points['X']
band_k_points = lineOfVectorPoints(start_point,end_point,25)

Here, we have used a function lineOfVectorPoints() (defined below), which takes three arguments

The function lineOfVectorPoints() returns a list of N points positioned on a straight line between the two given points. The definition of lineOfVectorPoints() is simple, namely

def lineOfVectorPoints(start_point,end_point,N=20):
    '''
    @return: A list containing (default: 20) k-points between two
    specified points in the Brillouin zone.
    '''
    from numpy import array
    points = []
    for i in range(N):
        new_point = (float(i)/(N-1),)*3*(array(end_point) - array(start_point))+start_point
        points.append(new_point)
    return points

The statement

    from numpy import array

is necessary in order to convert the tuples of points in the fcc_symmetry_points dictionary to numeric arrays; otherwise, we would not be able to calculate the difference between the two points. The returned list includes the end-point; to exclude it (which is sometimes desirable), simply use N instead of N-1 in the denominator.

Here are the lines of the script is the part where the actual band structure gets calculated

#Restore old SCF calculation
calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')
band_structure = calculateEnergyBands(
    self_consistent_calculation = calculation,
    k_point_list = band_k_points)

The input to the function calculateEnergyBands() should be prepared in the form of a list of points, given as tuples ; this input corresponds exactly to the list returned by the user defined function lineOfVectorPoints() described above, which sets up the sequence of points between the symmetry points.

The following three lines of code appends the energy bands of our silicon crystal to a list. This may then be used for exporting the band structure to other applications.

#Extract calculated bands into a list
energy_bands = []
for band_index in range(band_structure.numberOfBands()):
    energy_bands.append(band_structure.band(band_index))

The number of bands in the band structure depends on the crystal lattice basis (i.e. the atoms in the unit cell) as well as on the size of the basis set (in our case DoubleZetaDoublePolarized). The latter quantity is difficult to know beforehand had we not written the script ourselves, and we therefore use the method numberOfBands() to query for the number of bands contained in the band structure. Once we know the number of bands, we can easily iterate through them one by one, and extract the band structure for each band. To accomplish this, we use the method band(), which takes the band index as argument. The band index ranges from 0 to n-1, where n is the total number of bands.

We now combine the small pieces of code described above into the following complete script for calculating the energy bands of silicon.

from ATK.KohnSham import *

def lineOfVectorPoints(start_point,end_point,N=20):
    '''
    @return: A list containing (default: 20) k-points between two
    specified points in the Brillouin zone.
    '''
    from numpy import array
    points = []
    for i in range(N):
        new_point = (float(i)/(N-1),)*3*(array(end_point) - array(start_point))+start_point
        points.append(new_point)
    return points

# Define common symmetry points in the Brillouin zone of silicon
fcc_symmetry_points = {
    'G' : (0.000, 0.000, 0.000),
    'X' : (0.500, 0.000, 0.500), 
    'W' : (0.500, 0.250, 0.750), 
    'L' : (0.500, 0.500, 0.500), 
    'K' : (0.375, 0.375, 0.750), 
    'U' : (0.625, 0.250, 0.625) 
}

#Create line in Brillouin zone along which the
#band structure is calculated
start_point = fcc_symmetry_points['G']
end_point = fcc_symmetry_points['X']
band_k_points = lineOfVectorPoints(start_point,end_point,25)

#Restore old SCF calculation
calculation = restoreSelfConsistentCalculation('si_gga_dzp.nc')
band_structure = calculateEnergyBands(
    self_consistent_calculation = calculation,
    k_point_list = band_k_points)

#Extract calculated bands into a list
energy_bands = []
for band_index in range(band_structure.numberOfBands()):
    energy_bands.append(band_structure.band(band_index))

#print band energies in units of eV
for i in range(len(band_k_points)):
    print band_k_points[i][0],
    for band in range(band_structure.numberOfBands()):
        print energy_bands[band][i].inUnitsOf(eV),
    print

bulk-calculateBands.py

The script will print the energy bands of silicon to your terminal. If you want to store the information in a file, you should redirect the script output to a file in the following manner:

atk bulk-calculateBands.py > si_bands.dat

The first data column contains the value of (band_k_points[i][0]), followed by the energy eigenvalues for each band associated with the chosen k-point in each column.

The obtained band structure of silicon is shown in Figure 21. The plot indicates that silicon would have an indirect band gap of approximately 2 eV. Experimentally, however, silicon is known to have an indirect band gap of 1.13 eV.

The reason for this discrepancy is, in this case, that the SingleZeta basis set (which for silicon contains the sp³ valence electrons) is not versatile enough to describe a silicon crystal. In particular, it does not have enough angular degrees of freedom to correctly form the strong tetragonal covalent bonds present in the diamond structure of silicon. Going to DoubleZeta, one would only add radial degrees of freedom, and thus it is necessary as a minimum to have a SingleZetaPolarized (but ideally DoubleZetaPolarized) to get accurate results for silicon.

If you re-run the above calculation with the DoubleZetaPolarized basis set (remember to change the name of the checkpoint file!), you should find a band gap below 0.6 eV, i.e. the band gap is severely underestimated. This is a well known artifact of DFT calculations using either the LDA or GGA-PBE exchange correlation functional.