In the tutorial Li-H₂-Li two-probe system, many of the basic aspects of a two-probe set-up are described using two Li-chains as electrodes and a hydrogen molecule positioned in between the electrodes. One of the important parameters in that system is the lattice constant of the lithium chain. In this example, we will show how ATK can be used to calculate this quantity.

To follow the discussion and the examples presented in this tutorial, you should be familiar with basis functions, the mesh cut-off, k-point sampling, and defining Kohn-Sham methods. If you are unfamiliar with these concepts, we recommend that you go through the tutorial Calculating molecular properties before continuing with the material presented here.

We start by setting up the various script elements that are required to determine the lattice constant of an infinite lithium chain. In summary, we will have to put the following elements into play

Once these elements are in place, we will define a for-loop that calculates the total system energy as a function of the lattice constant. The lattice constant associated with the lowest system energy should correspond to the equilibrium lattice constant of the lithium chain.

basis_set = basisSetParameters(type=SingleZeta)

electron_density_params = electronDensityParameters(
    mesh_cutoff=100*Rydberg
    )

bz_integration = brillouinZoneIntegrationParameters( (1,1,100) )

dft_method = KohnShamMethod(
    basis_set_parameters=basis_set,
    brillouin_zone_integration_parameters=bz_integration,
    electron_density_parameters=electron_density_params
    )

In order to find the lattice constant for the lithium chain, we will use the following chain of events

The unit cell is defined as a Python list, where each list element holds the coordinates of three unit cell vectors a, b, and c. To obtain an effectively one-dimensional lithium chain, the interatomic distances determined by the a and b axis must be chosen sufficiently large in order to neglect the Li-Li interactions along these directions. Here is how to do it

    unit_cell = [[9.0, 0.0, 0.0],
                [ 0.0, 9.0, 0.0],
                [ 0.0, 0.0,  a ]]*Angstrom

We use a variable a for defining the length of the unit cell in the chain direction. This will allow us to conveniently control the lattice constant of the lithium chain. The unit cell is stored in the variable unit_cell. To generate a periodic repetition of the unit cell, we pass on unit_cell to the constructor of a PeriodicAtomConfiguration object

    li_chain = PeriodicAtomConfiguration(unit_cell,[Lithium], 
                                         [ (0.0, 0.0, 0.0)*Angstrom ] )

which we then store in the variable li_chain. The lithium chain is now defined and its total energy can be calculated using

    # Execute the self-consistent calculation
    dft_calculation = dft_method.apply(li_chain)
    
    # Calculate the system energy
    energy = calculateTotalEnergy(dft_calculation)

All we need to do now is to wrap the above code snippets inside a for-loop, which then calculates the system energy over a range of lattice constant samples. So, assembling the initial settings and the for-loop leaves us with the following complete script

from ATK.KohnSham import *

# Set the basis set to be SZ
basis_set = basisSetParameters(type=SingleZeta)

# Define the density grid cut-off
electron_density_params = electronDensityParameters(
    mesh_cutoff=100*Rydberg
    )

# Define the k-point sampling
bz_integration = brillouinZoneIntegrationParameters( (1,1,100) )

# Define the DFT method
dft_method = KohnShamMethod(
    basis_set_parameters=basis_set,
    brillouin_zone_integration_parameters=bz_integration,
    electron_density_parameters=electron_density_params
    )

# Open file for storing the results
f = open('energy.dat', 'w')

# Loop over lattice constants
import numpy
for a in numpy.arange(2.9,3.1,0.01):

    # Create the unit cell for 1D chain
    unit_cell = [[9.0, 0.0, 0.0],
                [ 0.0, 9.0, 0.0],
                [ 0.0, 0.0,  a ]]*Angstrom

    # Create the atomic configuration
    li_chain = PeriodicAtomConfiguration(unit_cell,[Lithium], 
                                         [ (0.0, 0.0, 0.0)*Angstrom ] )
    
    # Execute the self-consistent calculation
    dft_calculation = dft_method.apply(li_chain)
    
    # Calculate the system energy
    energy = calculateTotalEnergy(dft_calculation)
    
    # Write lattice constant and energy to file and screen 
    s = '%.4f %.10f \n' % (a, energy.inUnitsOf(Rydberg))
    f.write(s)
    print s,

# Close file and exit
f.close()

li-chain.py

As you may have noticed, the script also includes a few lines which enables us to store the calculated date to a file energy.dat. You will see shortly how we can use this for post-processing the calculated data.

Now, if you run the script, the following output is produced

2.9000 -0.7588874823 
2.9100 -0.7589779279 
2.9200 -0.7590611298 
2.9300 -0.7591372387 
2.9400 -0.7592063806 
2.9500 -0.7592685775 
2.9600 -0.7593239155 
2.9700 -0.7593723696 
2.9800 -0.7594139393 
2.9900 -0.7594491366 
3.0000 -0.7594781002 
3.0100 -0.7595002064 
3.0200 -0.7595164557 
3.0300 -0.7595263519 
3.0400 -0.7595303386 
3.0500 -0.7595281172 
3.0600 -0.7595202662 
3.0700 -0.7595062710 
3.0800 -0.7594865796 
3.0900 -0.7594613974 
3.1000 -0.7594307260

You should also check that the exact same output also has been written to the file energy.dat.

As is shown in Figure 35, the lithium chain system clearly passes an energetic minimum as the lattice constant is varied. Simple inspection of the data reveals that the minimum is located inside the interval [2.85;2.95]. At this point, we could choose the center of the interval 2.90 Å and then use this as the value of the lattice constant in subsequent calculations.

Figure 35: Plot illustrating the variation of the total energy of a one-dimensional chain of lithium atoms as a function of the lattice constant. The equilibrium lattice constant is determined as the minimum of the energy curve.

On the other hand, we may use a relatively simple approach to refine the determined equilibrium value of the lattice constant by fitting our sampled data points to a second order polynomial of the form

We may then choose the minimum of

as the value of the lattice constant. We now set (the central point of the interval) and let = 0.01 Rydberg (the equidistant spacing between sample points). In addition, using a few algebraic manipulations leads to the following analytic expressions for the two constants and

and

A simple script that implements a user-defined function QuadraticMin() for locating and refining the minimum is shown below

# Locate minimum in X,Y data set, fit to quadratic
# polynomial P(x) = A*x**2 + B*x + C, and find the minimum
# xmin of P(x).
def QuadraticMin(X,Y):
    n = len(X)
    h = X[1] - X[0]

    x0 = None

    # Locate possible minimum
    for i in range(1,n-1):
        if ( (Y[i-1] > Y[i]) and (Y[i] < Y[i+1]) ):
            x0 = X[i]
            (ym,y0,yp) = (Y[i-1], Y[i], Y[i+1])
            break

    # No minimum
    if (x0 == None):
        print '%s' % ('No minimum found')
        return

    # Fit polynomial to minimum interval
    A = (yp + ym - 2.0*y0)/(2.0*h*h)
    B = (yp - ym)/(2.0*h)
    xmin = x0 - B/(2.0*A)
    print 'Quadratic fit: xmin = %g' % (xmin)

    return


# Read two-column data from file
f = open('energy.dat', 'r')
data = f.readlines()
f.close()

# Partition data into separate lists 
X = [float(i.split()[0]) for i in data]
Y = [float(i.split()[1]) for i in data]

QuadraticMin(X,Y)

quad-min.py

First we read in the the column data from the file energy.dat, which we then pass on as the arguments to QuadraticMin(). Running the script outputs the following value for the improved lattice constant

Quadratic fit: xmin = 3.04142

in units of Angstrom.

A plot comparing the fitted polynomial with the sampled data points is shown in Figure 36. As you can see, the sampled energy points follow the quadratic dependence to a very good approximation; certainly, this justifies using the determined value of as an even better approximation for the lattice constant of the one-dimensional lithium chain.

Figure 36: Plot showing a second order polynomial (solid line) fitted to the data shown in Figure 35 (black circles). The minimum of the polynomial (red circle) is used for the lattice constant of the one-dimensional lithium chain.