Calculating molecular properties

	Tutorials
Geometry of a water molecule	Spin-polarized oxygen

Introduction

In this tutorial, we go into detail with some of the more advanced features of ATK. If you have studied the last sections of the tutorial Geometry of a water molecule, you already would have met some of these, namely

In the tutorial Geometry of a water molecule, we did not provide you with a detailed account of the above functions and classes - but here, we will explain their usage in detail. We will stick with the water molecule as our physical system, and use this as our test case. In particular, we will take a closer look at the following challenges

calculating total energy and forces
finding Kohn-Sham eigenvalue spectra
visualizing the Kohn-Sham eigenstates

Prerequisites

If ATK and NanoLanguage are still new to you, we suggest that you have a closer look at the Geometry of a water molecule tutorial before you continue. Once you have been updated on the basic concepts, you should be ready to follow and understand the ideas discussed here.

Getting started

Let's start with the basic script that was constructed in Geometry of a water molecule for setting up a water molecule, namely

# Import the KohnSham module from ATK
from ATK.KohnSham import *
   
# Define function for handling molecule setup   
def waterConfiguration(angle,bondLength):
    from math import sin,cos
       
    theta = angle.inUnitsOf(radians)
    positions = [(0.0, 0.0, 0.0)*Angstrom,
                 (1.0, 0.0, 0.0)*bondLength,
                 (cos(theta), sin(theta), 0.0)*bondLength]
    elements = [Oxygen] + [Hydrogen]*2
       
    return MoleculeConfiguration(elements,positions)
    
# Setup water molecule
h2o = waterConfiguration(104.5*degrees, 0.958*Angstrom)

# Open a VNL file and add the molecule to it
vnl_file = VNLFile("h2o.vnl")
vnl_file.addToSample(h2o, "h2o")

water02.py

where we have used the function waterConfiguration() for setting up a water molecule as a function of the O-H bond length and the H-O-H bending angle. The function waterConfiguration() returns a MoleculeConfiguration object, which we store in the variable h2o. Note, that we have used the equilibrium values for both the bond length and angular argument. The h2o object contains the complete definition of the water molecule. Finally, we save the molecule using a VNLFile object. We may then reload the molecule later from other scripts using the NanoLanguage construction

h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

The method readAtomicConfigurations() of the VNLFile object returns a dictionary containing the atomic configurations as values and the name they were stored under as the key. In the above example we restore the atomic configuration with the key h2o from the dictionary.

We will use the generated VNLFile h2o.vnl as a starting point for all subsequent calculations. Therefore, run the above script with the command

atk water02.py

in order for the VNLFile h2o.vnl to be created on your disk if it has not already been created. If the script fails because the file was already present and contained an identical atomic configuration, you may safely move on to the next step using that file.

Solving the electronic structure problem

The first task we will tackle for the water molecule, is to determine the electron structure of the molecule in its equilibrium state. In ATK, this is done by performing a self-consistent field (SCF) calculation for the molecule by solving the Kohn-Sham equations. We do that by adding the following two lines of code

method = KohnShamMethod()
scf = method.apply(h2o)

Even though these two lines of code may appear 'innocent' at first, they in fact instruct ATK to perform a series of relatively complicated steps. In the following, we will try to expose some of the gory details that actually takes place.

First, to solve the Kohn-Sham equations that describe the electronic structure of the water molecule, we create a KohnShamMethod object, which we assign to the variable method.

After we have created the KohnShamMethod object, we can use its functionality to determine the electronic structure of the atomic configuration (in this case, a water molecule). To do this, we use the method apply(), which is a method of the KohnShamMethod object. When KohnShamMethod is called, we must provide information about the specific atomic configuration that we wish to study. In this case, we would like to study an object of the type MoleculeConfiguration, namely, the water molecule. Since the water molecule object was assigned to the variable h2o, we perform the self-consistent field calculation by passing h2o as an argument to KohnShamMethod - in other words

scf = method.apply(h2o)

When we use the apply() method of the KohnShamMethod object, apply() returns a SelfConsistentCalculation object. In our case, we have assigned this object to the variable called scf. As we will see, ATK contains several analysis functions for querying a wide variety of properties related to the solution of the Kohn-Sham equations. One of these is calculateTotalEnergy(), which returns the total energy of the system at hand.

Note, that ATK introduces a clear-cut separation between the atomic configuration and the specific method that we use to solve the electronic structure problem. This gives you several advantages. For example, you may

create several atomic configurations, and use the same method to solve the electronic structure problem for each configuration
create several methods to solve the electronic structure problem for one particular atomic configuration

Here is a complete script summarizing the key ideas discussed so far using calculateTotalEnergy() to obtain the total system energy

# Import the KohnSham module from ATK
from ATK.KohnSham import *
   
# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

# Perform a selc-consistent calculation
method = KohnShamMethod()
scf = method.apply(h2o)

print 'Total energy = ', calculateTotalEnergy(scf)

water02-etot.py

To run the script, type the following at the command line

$ atk water02-etot.py

which results in the following output once the calculation finishes (this will take some time so be patient):

Total energy =  -465.269986233 eV

In short, the method variable which we assigned to the object KohnShamMethod takes care of all details related to the actual calculation. In practical terms, this involves

creating basis functions to represent the solution of the Kohn-Sham equations
loading adequate pseudo-potentials shipped with ATK for treating the interaction of the Kohn-Sham states with the screened nuclei
creating real-space grids used to represent the effective Kohn-Sham potential and electron density

So far, none of these elements have been exposed - we completely ignored these finer details, and left everything to their internal default values. Often, however, using pure default settings for a calculation will not suffice to achieve acceptable results. In the next section, we will show how basis functions, pseudo-potentials, and grid settings can be fully controlled by fine-tuning various method parameters.

Increasing verbosity

Calculating the total energy of water the way we did above, takes a lot of time where nothing appears to happen. It is possible to supervise the convergence of the self-consistent calculation by increasing the verbosity level of ATK during the calculation. By inserting the following two lines

import ATK
ATK.setVerbosityLevel(1)

at the top of your script, ATK will be more verbose about how the SCF-calculation is proceeding.

If you insert the two lines above, your script should look like this:

# Import the KohnSham module from ATK
from ATK.KohnSham import *
import ATK
ATK.setVerbosityLevel(1)
   
# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

# Perform a selc-consistent calculation
method = KohnShamMethod()
scf = method.apply(h2o)

print 'Total energy = ', calculateTotalEnergy(scf)

water02-verb.py

and running it gives rise to the following output during the calculation:

# sc  0 : Fermi Energy =    0.00000 Ry
# sc  1 : Fermi Energy =   -0.21613 Ry  dRho =  1.0317E+00
# sc  2 : Fermi Energy =   -0.17129 Ry  dRho =  6.1782E-02
# sc  3 : Fermi Energy =   -0.01558 Ry  dRho =  2.1349E-01
# sc  4 : Fermi Energy =   -0.03431 Ry  dRho =  1.0472E-01
# sc  5 : Fermi Energy =    0.00167 Ry  dRho =  5.4534E-02
# sc  6 : Fermi Energy =    0.00428 Ry  dRho =  4.9847E-03
# sc  7 : Fermi Energy =    0.00590 Ry  Etot =  -34.19669 Ry  dRho =  2.5855E-03
# sc  8 : Fermi Energy =    0.00561 Ry  Etot =  -34.19669 Ry  dRho =  6.2618E-04  dEtot = -2.3161E-07 Ry
Total energy =  -465.269986229 eV
"""

Investigating parameter effects

When we study a system in ATK, there are several parameters that control the accuracy of the numerical calculations. So far, we have relied on the default settings of these, but often you have to fine tune several parameters to obtain results with sufficient accuracy. In other situations, the default settings are far too strict for the problem at hand resulting in lengthy calculations, that possibly could have been carried out much faster by altering a few critical parameters.

In this section, we will take a closer look at the following important quantities

Mesh cut-off
Basis set size
Exchange-correlation function

and describe a few general ATK constructions for changing parameter setups for your calculations. We start by looking at how we control the mesh cut-off.

Optimizing the real-space grid density

In the solution of the Kohn-Sham equation, the electron density and the Kohn-Sham potential are calculated on a real space grid. The accuracy of the calculation can be increased by increasing the resolution of the real-space grid. The grid resolution is defined by the so-called mesh cut-off, which is an energy that determines the spacing ( $\Delta x$ ) between real space grid points according to this formula

$\displaystyle \Delta x = \frac{h}{\sqrt{2m_eE}}$

The default mesh cut-off used in ATK is 150 Rydberg.

It should be noted that the size of the mesh cut-off generally is determined by the chemical species of the system and the associated pseudo-potentials, implying that the mesh cut-off generally can be determined using a small basis set and any one of the multiple exchange-correlation functionals.

Suppose that we would like to find the mesh cut-off most suitable for our water molecule system. A strategy for solving this task could for example be

generate a list of mesh cut-offs
iterate through them using a for loop
initialize the Kohn-Sham method for each mesh cut-off
calculate the total energy for each specific mesh cut-off
locate the mesh cut-off where the total system energy does not change anymore

The first thing we should do, is to set up the geometry of our water molecule. You have seen this before: We use the user-defined function waterConfiguration() and store the returned MoleculeConfiguration in the variable h2o. The next task, is then to use this object, and calculate its total energy as a function of the mesh cut-off.

When we construct the KohnShamMethod object with no arguments, like

method = KohnShamMethod()

all parameter values are just initialized to their default setting. To override this, we must supply our preferred parameter settings when the KohnShamMethod is constructed. In ATK, there is a broad list of parameter arguments that can passed on to KohnShamMethod. Note that all of these are defined as optional variables. The one which is relevant for setting the mesh cut-off is the variable electron_density_parameters. The is done by first using the function electronDensityParameters(), which is invoked like this

mcut = 200*Rydberg
el = electronDensityParameters(mcut)

and then passing the variable el to the KohnShamMethod object upon its construction. In other words, the code snippet

mcut = 200*Rydberg
el = electronDensityParameters(mesh_cutoff=mcut)
method = KohnShamMethod(electron_density_parameters=el)

prepares a run, where the Kohn-Sham equations will be solved on a grid using a mesh cut-off equal to 200 Rydberg. Note that electronDensityParameters() also accepts variable specified by their name.

All we need to do to complete our task, is to nest this code structure inside a for loop by iterating through a list of mesh cut-off values. So, putting all the above pieces together, here is a script for locating the minimum cut-off

# Import the KohnSham module from ATK
from ATK.KohnSham import *
   
# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

print '    MeshCutoff (Rydberg)     TotalEnergy (eV)'
print '------------------------------------------------'
for mcut in range(30,200,10)*Rydberg:
    el = electronDensityParameters(mesh_cutoff=mcut)
    method = KohnShamMethod(electron_density_parameters=el)
    print '\t', mcut.inUnitsOf(Rydberg), '\t\t\t', \
          calculateTotalEnergy(method.apply(h2o)).inUnitsOf(eV)

water02-mesh.py

Running the script generates the output

    MeshCutoff (Rydberg)     TotalEnergy (eV)
------------------------------------------------
        30.0                    -465.732523778
        40.0                    -465.480505002
        50.0                    -465.207166349
        60.0                    -465.245531699
        70.0                    -465.290437537
        80.0                    -465.273032958
        90.0                    -465.270366088
        100.0                   -465.269134142
        110.0                   -465.268292704
        120.0                   -465.269941954
        130.0                   -465.27032458
        140.0                   -465.270140252
        150.0                   -465.269986233
        160.0                   -465.270026007
        170.0                   -465.269922622
        180.0                   -465.269886349
        190.0                   -465.269939292

Apparently, choosing a mesh cut-off larger than 90 Rydberg will have no significant effect on the calculated energy. In other words, the default setting uses a mesh cut-off equal to 150 Rydberg, but in our case, this is not necessary. We may safely set the mesh cut-off to 90 Rydberg and thereby gain a significant boost in the CPU time required to complete the calculation.

When you, like in the example above, incrementally converge towards a more and more accurate value it can save time, if the results from the previous calculation are reused i.e. the calculation using a mesh cut-off of 200 Rydberg took its starting point in the calculation using 190 Rydberg. This can be done by slightly altering the above example:

# Import the KohnSham module from ATK
from ATK.KohnSham import *
import datetime
# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

print '    MeshCutoff (Rydberg)     TotalEnergy (eV)'
print '------------------------------------------------'
use_old=False
for mcut in range(30,200,10)*Rydberg:
    el = electronDensityParameters(mesh_cutoff=mcut)
    method = KohnShamMethod(electron_density_parameters=el)
    if use_old:
        density = executeSelfConsistentCalculation(
            method=method,
            atomic_configuration=h2o,
            initial_calculation=old_scf
            )
    else:
        density = executeSelfConsistentCalculation(
            method=method,
            atomic_configuration=h2o,
            )
        use_old=True
    old_scf=density
    print '\t', mcut.inUnitsOf(Rydberg), '\t\t\t', \
          calculateTotalEnergy(density).inUnitsOf(eV)

water02-mesh-restart.py

Here, we have replaced the apply() method with the executeSelfConsistentCalculation() method and placed it inside a conditional code structure called an if-statement. What basically happens is that we check to see if a previous calculation has been performed. If results from a previous calculation are available (use_old=True) we use those as the starting point. Note the named argument initial_calculation which is assigned the value of old_scf. Otherwise, we calculate the electron density using the default initial starting point. By using previous calculations as the initial guess in a similar subsequent calculation you can, depending on the situation, save considerable calculation time.

Selecting a suitable basis set

In ATK, you can choose among five different basis set types, namely

where the basis sets are ordered with respect to increasing completeness, which is a term used for how many functions that they use to describe the electron density at the atoms. The default basis set in ATK is the DoubleZetaPolarized and we recommend that you use this basis set or the more complete DoubleZetaDoublePolarized unless you know what you are doing or run into memory or time problems when performing the calculation. You can find a more detailed description of the basis sets in the AtomicOrbitals section of the reference manual.

When we defined the mesh cut-off, we used the optional variable electron_density_parameters and the helper function electronDensityParameters(). ATK offers an equivalent construction for choosing the basis set type for a given molecule, namely basis_set_parameters and the helper function electronDensityParameters(). Here is how we use them, if you want to to override the default setting and choose the basis set type SingleZeta for all atoms in the molecule

bsp = basisSetParameters(type=SingleZeta)
method = KohnShamMethod(basis_set_parameters=bsp)

Let us calculate the total energy of the water molecule for the different choices of basis. To do this, we set up a list containing a short-hand notation string for each basis set as well as the basis set itself.By iterating over all list elements, we may calculate the total system energy for each of the basis set types available in ATK. Here is the script

# Import the KohnSham module from ATK
from ATK.KohnSham import *
   
# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

# List of basis set and their names 
BasisSetList = [['SZ'  , SingleZeta],              
                ['SZP' , SingleZetaPolarized],
                ['DZ'  , DoubleZeta],              
                ['DZP' , DoubleZetaPolarized],     
                ['DZDP', DoubleZetaDoublePolarized]]

# Variable for holding the previous energy
E0 = None

print '    Basis set    TotalEnergy (eV)    Difference    '
print '---------------------------------------------------'

for basis in BasisSetList:
    name,basis_type = basis
    bsp = basisSetParameters(type=basis_type)
    method = KohnShamMethod(basis_set_parameters=bsp)
    E = calculateTotalEnergy(method.apply(h2o)).inUnitsOf(eV)

    if E0:
        print "     %-9s %15.6f %15.6f" % (name, E, (E-E0))
    else:
        print "     %-9s %15.6f"        % (name, E)
    E0 = E

water02-basis.py

which gives the output (where the difference between the energies of the current and the previous basis also is given)

    Basis set    TotalEnergy (eV)    Difference    
---------------------------------------------------
     SZ            -461.465517
     SZP           -463.712300       -2.246783
     DZ            -463.731880       -0.019580
     DZP           -465.269986       -1.538106
     DZDP          -465.639658       -0.369672

As you can see, the energy differences keep decreasing as the accuracy and completeness of the basis set increases in consistency with the variational principle.

The exchange correlation functional

In ATK, there are three exchange-correlation functionals available, namely

LDA.PZ : Perdew-Zunger parametrization of the local density approximation.
GGA.PBE: Perdew-Burke-Ernzerhof parametrization of the generalized gradient approximation.
GGA.revPBE: Revised Perdew-Burke-Ernzerhof parametrization of the generalized gradient approximation.

The default exchange correlation functional is set to LDA.PZ. In the previous examples, we defined parameter values using an optional variable and a helper function. In this case, things are even simpler, since we choose the exchange correlation functional when the KohnShamMethod object is constructed. Here is what to do

method = KohnShamMethod(exchange_correlation_type=GGA.PBE)

From the mesh cut-off example, we already know that a mesh cut-off of the order 100 is sufficient for the water molecule. In addition, the previous example demonstrated that the basis set DoubleZetaDoublePolarized is the superior choice in ATK. So, let us finish things off with a somewhat more complicated example, where the exchange correlation functional is involved. Here is our strategy

Set up a water molecule
Define the Kohn-Sham method object KohnShamMethod with the basis set DoubleZetaDoublePolarized and the optimized value for the mesh cut-off.
For each exchange correlation functional, calculate the equilibrium values of the O-H bond length and H-O-H bending angle using a self-consistent field calculation.

Solving the first two points are simple, and we have already seen several examples of setting this up. The third point, performing a self-consistent field calculation, is accomplished by the ATK function calculateOptimizedAtomicGeometry(). This function takes several arguments, and in our case, we will have to pass

a water molecule object (we generate this using the user defined function waterConfiguration())
a method object for solving the Kohn-Sham equations (this will be the object returned by KohnShamMethod)

Additional parameters for calculateOptimizedAtomicGeometry() can be set using the helper function geometricOptimizationParameters() using constructions completely similar to the ones used for setting the mesh cut-off and choosing the basis set type. Here we will just settle for the default values.

A script that solves the task at hand is displayed here

# Import the KohnSham module from ATK
from ATK.KohnSham import *

# Function for handling molecule setup   
def waterConfiguration(angle,bondLength):
    from math import sin,cos

    theta = angle.inUnitsOf(radians)
    positions = [(0.0, 0.0, 0.0)*Angstrom,
                 (1.0, 0.0, 0.0)*bondLength,
                 (cos(theta), sin(theta), 0.0)*bondLength]
    elements = [Oxygen] + [Hydrogen]*2

    return MoleculeConfiguration(elements,positions)


# Function for calculating H-O-H bond angle
def waterBendingAngle(water):
    from math import acos

    coords = water.cartesianCoordinates()
    oh1 = coords[1] - coords[0]
    oh2 = coords[2] - coords[0]
    prod = dot(oh1,oh2)/( vLength(oh1)*vLength(oh2) )

    return acos(prod)*radians


# Function for calculating O-H bond length
def waterBondLength(water):
    coords = water.cartesianCoordinates()
    oh1 = coords[1] - coords[0]

    return vLength(oh1)


# Dot product between two vectors
def dot(u,v):
    return (u*v).sum()
    

# Length of a vector containing elements of type PhysicalQuantity
def vLength(u):
    from math import sqrt

    return sqrt(((u*u).sum()).inUnitsOf(Angstrom**2))*Angstrom

# We are going to loop through this dictionary
ExCorDict = {
    'LDA.PZ'     : LDA.PZ,
    'GGA.PBE'    : GGA.PBE,
    'GGA.revPBE' : GGA.revPBE
    }

# Setup water molecule
h2o = waterConfiguration(100.0*degrees, 1.1*Angstrom)

# Set basis set and mesh cutoff
el  = electronDensityParameters(mesh_cutoff=100*Rydberg)
bsp = basisSetParameters(type=DoubleZetaDoublePolarized)

# Choose default geometric optimization parameters
geometric_optimization_parameters = geometricOptimizationParameters()

print '    Functional     H-O-H angle/deg      O-H Bondlength/A    '
print '------------------------------------------------------------'

for (name,xcf) in ExCorDict.iteritems():
    # Setup Kohn-Sham method
    dft_method = KohnShamMethod(
        exchange_correlation_type=xcf,
        basis_set_parameters=bsp,
        electron_density_parameters=el
        )

    h2o_opt = calculateOptimizedAtomicGeometry(
        atomic_configuration=h2o,
        method=dft_method,
        optimization_parameters=geometric_optimization_parameters
        )

    angle  = waterBendingAngle(h2o_opt).inUnitsOf(degrees)
    length = waterBondLength(h2o_opt).inUnitsOf(Angstrom)

    print "     %-10s %15.3f %19.3f" % (name, angle, length)

water02-opt.py

Some of the structures and helper functions used in the script were introduced in this section in the tutorial Geometry of a water molecule. Notice also, how the Python dict method itername() is used for accessing both the key and value of ExCorDict.

If you run the script, the following equilibrium values are obtained for the three different exchange correlations functional

    Functional     H-O-H angle/deg      O-H Bondlength/A
------------------------------------------------------------
     LDA.PZ             104.741               0.973
     GGA.revPBE         104.063               0.968
     GGA.PBE            104.008               0.968

The Kohn-Sham eigenvalue spectrum

Another relevant property that often is important to study is the eigenvalue spectrum of a given molecule or structure. In ATK, we can do this by using the analysis function calculateMolecularEnergySpectrum(). The steps you have to take are simple:

Set up the KohnShamMethod object.
Apply the method object to a given structure using the method apply() from the KohnShamMethod object.
Use the object returned by apply() as input to the function calculateMolecularEnergySpectrum().

A simple test script, where we compare the energy spectra obtained using two different exchange correlation functionals, could for example be constructed like this

# Import the KohnSham module from ATK
from ATK.KohnSham import *

# Load molecule from VNLFile
h2o = VNLFile("h2o.vnl").readAtomicConfigurations()["h2o"]

# Set up Kohn-Sham for two exchange correlation functionals
gga_calculation = KohnShamMethod(exchange_correlation_type=GGA.PBE).apply(h2o)
lda_calculation = KohnShamMethod(exchange_correlation_type=LDA.PZ).apply(h2o)

# Perform the actual calculations
gga_spectrum = calculateMolecularEnergySpectrum(gga_calculation).energies()
lda_spectrum = calculateMolecularEnergySpectrum(lda_calculation).energies()

# eV = electronVolt

print '   Energy eigenvalues for water (in eV)\n'
print '-----------------------------------------'
print '  LDA\t\t  GGA\t\t  diff\t'
print '-----------------------------------------'

for (E1,E2) in zip(lda_spectrum,gga_spectrum):
    e1 = E1.inUnitsOf(eV)
    e2 = E2.inUnitsOf(eV)
    print '%7.3f\t\t%7.3f\t\t%7.3f' % (e1,e2,(e1-e2))

water02-egv.py

In the above for loop, we have used the built-in Python function zip(): zip() takes two list arguments, whose elements are returned in a new list ordered as pairs of Python tuples. If you run the script, ATK outputs the following eigenvalue spectrum for the water molecule

   Energy eigenvalues for water (in eV)

-----------------------------------------
  LDA             GGA             diff
-----------------------------------------
-24.243         -24.376           0.133
-12.470         -12.442          -0.028
 -8.187          -8.129          -0.058
 -6.140          -5.988          -0.152
  1.811           1.977          -0.166
  4.195           4.342          -0.147
 11.994          12.330          -0.336
 13.019          13.530          -0.511
 17.469          18.449          -0.980
 19.089          20.234          -1.145
 19.532          20.539          -1.007
 25.733          26.771          -1.038
 26.259          26.906          -0.647
 28.061          28.973          -0.912
 33.451          34.662          -1.211
 35.165          36.330          -1.165
 47.420          48.609          -1.189
 50.064          50.683          -0.619
 65.164          65.309          -0.145
 66.721          66.836          -0.115
 71.480          72.006          -0.526
 89.808          89.971          -0.163
108.047         108.790          -0.743

Visualizing Kohn-Sham eigenstates

The next issue we will touch upon is how to use ATK to calculate molecular orbitals. The orbitals are the eigenstates (or eigenfunctions) whose respective eigenvalues correspond to the energy levels that we calculated in the previous section. In its ground state, the eight electrons of the water molecule will occupy the four orbitals with the lowest energy levels. Examining the output and using the data from the calculation employing the LDA exchange correlation functional, we see that these orbitals correspond to the energy levels -24.243 eV, -12.470 eV, -8.187 eV, and -6.140 eV. The orbital with the highest energy (-6.140 eV) is the so-called HOMO (Highest Occupied Molecular Orbital), whereas the so-called LUMO (Lowest Unoccupied Molecular Orbital) corresponds to the energy level of 1.811 eV. The properties of these two orbitals are often especially important, because electron transitions typically takes place between these to states.

Often, you also want to visualize the geometric properties of certain orbitals. In ATK, this is done using the following procedure:

# calculate HOMO and LUMO eigenstates
egvs = calculateEigenstates(
    self_consistent_calculation=scf,
    quantum_numbers=( [3,4] )
    )

The above procedure contains the following steps:

First we perform a self-consistent field calculation using the water molecule stored in our VNLFile and the apply() method of the KohnShamMethod object. The result of this calculation is stored in the variable scf.
The value of the variable scf is then passed on to the function calculateEigenstates(). The second argument to calculateEigenstates() is a list of quantum numbers specifying the orbitals we want to be calculated; since Python list items are numbered starting from 0, we provide a list contain the numbers 3 (HOMO) and 4 (LUMO). calculateEigenstates() returns a list, where each element contains the data needed for plotting the respective orbitals. We store this information in the variable egvs.

Finally, the content of egvs can then be added to a VNLFile using the method addToSample().

results = VNLFile('results.vnl')
results.addToSample(h2o, 'h2o')
results.addToSample(egvs[0], 'h2o','0')
results.addToSample(egvs[1], 'h2o','1')

The generated VNLFile can then be loaded into the Virtual NanoLab (VNL), and used to generate a 3D rendering of the HOMO and LUMO states. The entire script reads

# Import the KohnSham module from ATK
from ATK.KohnSham import *

# Load molecule from VNLFile
h2o = VNLFile('h2o.vnl').readAtomicConfigurations()['h2o']

method = KohnShamMethod()
scf    = method.apply(h2o)

# calculate HOMO and LUMO eigenstates
egvs = calculateEigenstates(
    self_consistent_calculation=scf,
    quantum_numbers=( [3,4] )
    )

results = VNLFile('results.vnl')
results.addToSample(h2o, 'h2o')
results.addToSample(egvs[0], 'h2o','0')
results.addToSample(egvs[1], 'h2o','1')

water02-orb.py

Figure 13: The HOMO state of the water molecule rendered using VNL.

Figure 14: The LUMO state of the water molecule rendered using VNL.

Mulliken populations

Our final example in this tutorial illustrates how to calculate the Mulliken population for a given molecule or structure. The outcome of a Mulliken population analysis is a description of the electron charge distribution for each individual atom in the molecule. In ATK, this analysis is performed using the function calculateMullikenPopulation(), which takes the object returned from a self-consistent field calculation as argument. Here is the principal usage

# Define and perform DFT
dft = KohnShamMethod()
scf = executeSelfConsistentCalculation()

# Calculate the Mulliken population per atom:
mulliken = calculateMullikenPopulation(scf).toArray()

First we perform a self-consistent field calculation, whose return object is stored in the variable scf. The variable scf is then passed on to calculateMullikenPopulation(), which returns a MullikenPopulation object. To access the data, we use the the method toArray() of the MullikenPopulation object, which returns a list containing the Mulliken population. This can then be passed on for further analysis.

Here is a complete script which prints the values of the Mulliken population to the screen

# Import the KohnSham module from ATK
from ATK.KohnSham import *

# Load molecule from VNLFile
h2o = VNLFile('h2o.vnl').readAtomicConfigurations()['h2o']

# Specify high verbosity 
runtime_params = runtimeParameters(
    verbosity_level = 10
    )

# Define and perform DFT
dft = KohnShamMethod()
scf = executeSelfConsistentCalculation(
    h2o,
    dft,
    runtime_params
    )

# Calculate the Mulliken population per atom:
mulliken = calculateMullikenPopulation(scf).toArray()

# Print Mulliken population to a text file:
for (elm,pop) in zip( h2o.elements(), mulliken ):
    print '%s: %s' % (elm.symbol(),pop)

water02-mul.py

Running it through ATK produces the following output

O: 5.88901735645
H: 1.05566628484
H: 1.05531635871

where the numbers are ordered according to the atom ordering O, H, and H. Clearly, these numbers deviate very little from the original electron charge of the separate atoms, illustrating that the water molecule only exhibits a (very) weak electron charge polarization.