Table of Contents
In this example, you will learn how to
set up periodic structures, exemplified by a silicon crystal,
calculate eigenstates and band structure of silicon
using NanoLanguage. We will also go into some of the numerical details specific for periodic systems, since the understanding of these is important in order to calculate the systems correctly with ATK.
It will be assumed that you are familiar with the most basic operations in NanoLanguage, such as
specifying chemical elements,
setting up atom positions,
saving samples to VNL files,
using some of the functionality in the ATK.KohnSham module.
If these concepts are new to you, we recommend that you consult the tutorial Geometry of a water molecule before you continue any further.
A crystal is characterized by being periodic. If we take the unit cell associated with the crystal, and repeat this in all three spatial dimensions, this operation will generate the actual crystal. Put differently, we say that the crystal is invariant to translational symmetry operations. Our first task will be to see how such a periodic structure is constructed using ATK.
We will do this step-by-step by creating a NanoLanguage script which
Sets up a bulk silicon crystal
Prints out the crystal configuration
Saves the crystal sample to a VNLFile for later use.
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
absolute Cartesian coordinates
relative coordinates expressed in the three unit vectors
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
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")
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.
Figure 18: The structure of the silicon fcc crystal, plotted using Virtual NanoLab. The unit cell of the crystal is indicated with a wire box.
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'
Before we actually run the script, a few remarks are in place regarding its structure
First, we load the silicon sample from our previous run and store this
configuration in the variable si_fcc
.
The basis set and the mesh cut-off are stored in the variables
basis
and mesh
using the respective
functions basisSetParameters() and electronDensityParameters(). Note, to speed up
the calculation a bit, we have chosen both a small basis set and mesh cut-off.
A Python while
loop is constructed. In each iteration step,
the number of k-points is incremented by 1 in all directions, after which the
total system energy E
is calculated. When the relative
difference between the current energy E
and the energy
Eprev
of the previous step is lower than the variable
convergence_criterion
, the while-loop is terminated.
In each iteration step, the method for the self-consistent calculation is set
using a KohnShamMethod object, which is
initialized with the chosen basis set, mesh cut-off, and the number of
k-points as specified by the respective variables basis
,
mesh
, and kgrid
. The returned method
object for solving the Kohn-Sham equations is assigned to the variable
method
.
A self-consistent calculation is now performed for the silicon system. This is done in two steps:
For the first set of k-points (1,1,1), we calculate the electron
density from scratch and change the Boolean value of
use_old
from False
to
True
.
In every subsequent calculation we reuse the density calculated in the
previous calculation by passing the variable old_scf
to executeSelfConsistentCalculation() using the optional keyword
initial_calculation
. This reduces the number of steps
required for the calculation to converge.
In every step we pass the system si_fcc
as the argument to
executeSelfConsistentCalculation() along
with the parameters previously specified in our dft_method
.
The total energy is evaluated using the function calculateTotalEnergy(), which is called with
the object calculation
returned by
the executeSelfConsistentCalculation()
function.
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.
In this example, we will demonstrate some general steps involved in ATK based band structure calculations. To illustrate the basic principles, it is sufficient to calculate and plot the band structure for a small portion of the Brillouin zone. The power of NanoLanguage, however, easily enables you to create a script that generates and plots the band structure for the entire zone, for example in order to produce a plot similar to the one shown for silicon in Figure 19.
Figure 19: The band structure of silicon, plotted in the entire Brillouin zone along the traditional symmetry directions for the fcc lattice. The data and plot script for the 2D plotting program gnuplot were generated entirely with NanoLanguage.
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)
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")
The quantum numbers given to the VNLFile function are given as a tuple sequence of
the band index (a non-negative integer, less than the total number of bands), and
k-points, in units of the reciprocal lattice vectors.
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 first two arguments are the two points that should be connected.
The third argument is an optional integer N
, that specifies
the number of points to include in the path between the two points
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
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.
Figure 21: The GGA band structure of silicon between the Γ and X point, generated from the data produced by the presented NanoLanguage script.
The reason for this discrepancy is, in this case, that the SingleZeta basis set (which for silicon contains the sp^{3} 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.