The following tutorial demonstrates how to perform a simple spin-polarized calculation on a molecule. In summary, this involves

In this example, we assume you are familiar with the process of creating basic atomic configurations in NanoLanguage as described in Geometry of a water molecule. If you are not familiar with the basic concepts of NanoLanguage we encourage you to study that tutorial first. We will also use Python functions and for-loops to automate and enhance our NanoLanguage script.

We start out by creating a function createOxygenMolecule() which automatically generates our molecular configuration taking a single parameter as input namely the bond distance.

def createOxygenMolecule(atom_separation):
    '''
    Simple function for creating a custom
    MoleculeConfiguration
    @param atom_separation - bond-length between
           oxygen atoms
    '''
    atoms = 2*[Oxygen]
    coordinates =  [
        ( 0.0, 0.0, atom_separation/2) ,
        ( 0.0, 0.0,-atom_separation/2)
        ]*Ang
    oxygen_molecule = MoleculeConfiguration(
        elements = atoms,
        cartesian_coordinates = coordinates
        )
    return oxygen_molecule

Instead of explicitly specifying the position of the atoms using coordinates, we specify the inter-atomic distance and use this to construct our MoleculeConfiguration object.

To anticipate what we will do later, we will construct a method which generates and initializes our KohnShamMethod.

def createKohnShamMethod():
    '''
    @return Predefined KohnShamMethod
    '''
    dft_method=KohnShamMethod( 
        basis_set_parameters=basisSetParameters(),
        exchange_correlation_type=GGA.PBE,
        )
    return dft_method

The KohnShamMethod object returned from this method has an non-spin-polarized electron density and a default basis set (DoubleZetaPolarized). We will use this method to calculate the total energy of the oxygen molecule for different bond distances. This will tell us what the equilibrium distance between the two atoms is. The following lines of NanoLanguage show how this could be done

# Create a dictionary for storing bond lengths and
# associated bond lengths
from numpy import arange
energy_dict ={}

print 'Separation (Angstrom)\tTotal Energy (eV)' 
# Iterate over a list of increasing bond lengths
for distance in  arange(1.10,1.30,0.02):
    dft_method = createKohnShamMethod()
    scf_oxygen = dft_method.apply(
        createOxygenMolecule(distance)
        )
    total_energy = calculateTotalEnergy(scf_oxygen).inUnitsOf(eV)
    print "%4.2f %s %9.4f"%(distance,'\t\t\t',calculateTotalEnergy(scf_oxygen).inUnitsOf(eV))
    energy_dict[total_energy] = distance

# Find the bond length with the minimum energy
energy_list = energy_dict.keys()
energy_list.sort()
print "Minimum total energy: %9.4f eV" % (energy_list[0])
print 'O-O equilibrium distance:',energy_dict[energy_list[0]],'A'

First, we import the arange() function from the NumPy module. This function is similar to the range() function as it is used to create a list of non-integer numbers by specifying a start, an end and an interval value. In this case a list containing numbers ranging from 1.1 to 1.3 using 0.01 as the step size. For each iteration step, the for-loop will create an oxygen molecule with an increasing inter-atomic distance, perform a self-consistent field (SCF) calculation for this molecule and, finally, calculate the total energy. All of this information is printed to your terminal while the script is being processed by ATK. The complete script is given below

from ATK.KohnSham import *

def createOxygenMolecule(atom_separation):
    '''
    Simple function for creating a custom
    MoleculeConfiguration
    @param atom_separation - bond-length between
           oxygen atoms
    '''
    atoms = 2*[Oxygen]
    coordinates =  [
        ( 0.0, 0.0, atom_separation/2) ,
        ( 0.0, 0.0,-atom_separation/2)
        ]*Ang
    oxygen_molecule = MoleculeConfiguration(
        elements = atoms,
        cartesian_coordinates = coordinates
        )
    return oxygen_molecule

def createKohnShamMethod():
    '''
    @return Predefined KohnShamMethod
    '''
    dft_method=KohnShamMethod( 
        basis_set_parameters=basisSetParameters(),
        exchange_correlation_type=GGA.PBE,
        )
    return dft_method

# Create a dictionary for storing bond lengths and
# associated bond lengths
from numpy import arange
energy_dict ={}

print 'Separation (Angstrom)\tTotal Energy (eV)' 
# Iterate over a list of increasing bond lengths
for distance in  arange(1.10,1.30,0.02):
    dft_method = createKohnShamMethod()
    scf_oxygen = dft_method.apply(
        createOxygenMolecule(distance)
        )
    total_energy = calculateTotalEnergy(scf_oxygen).inUnitsOf(eV)
    print "%4.2f %s %9.4f"%(distance,'\t\t\t',calculateTotalEnergy(scf_oxygen).inUnitsOf(eV))
    energy_dict[total_energy] = distance

# Find the bond length with the minimum energy
energy_list = energy_dict.keys()
energy_list.sort()
print "Minimum total energy: %9.4f eV" % (energy_list[0])
print 'O-O equilibrium distance:',energy_dict[energy_list[0]],'A'

oxygen.py

The output printed to your screen should resemble the following:

Separation (Angstrom)   Total Energy (eV)
1.10                     -863.3836
1.12                     -863.6485
1.14                     -863.8520
1.16                     -864.0021
1.18                     -864.1057
1.20                     -864.1728
1.22                     -864.2015
1.24                     -864.1990
1.26                     -864.1696
1.28                     -864.1169
Minimum total energy: -864.2015 eV
O-O equilibrium distance: 1.22 A

The equilibrium bond length is the separation at minimum total energy and is given directly on the screen. Comparing with the experimentally determined equilibrium bond distance reported in the Computational Chemistry Comparison and Benchmark Database, we find good agreement between our calculated value and experiment (1.208 Angstrom).

We will now explore, how the introduction of a spin-polarized electron density changes the bond distance and total energy of the oxygen molecule.

Before we can use the spin-polarized electron density, we need to modify our createKohnShamMethod() function. We want to define the degree of spin-polarization on initializing our KohnShamMethod so we add the optional argument spin to the createKohnShamMethod() function. The modified function now consists of the following lines:

def createKohnShamMethod(spin):
    '''
    @param spin - Total spin for the electron density.
           Supplied with no unit.
    '''
    fixed_electron_density = electronDensityParameters(
        fixed_spin = spin*hbar
        )
    
    dft_method=KohnShamMethod(
        basis_set_parameters=basisSetParameters(),
        exchange_correlation_type=GGA.PBE,
        electron_density_parameters = fixed_electron_density
        )
    
    return dft_method

The spin argument mentioned above is used to set the electronDensityParameters(), which, when combined with the KohnShamMethod, will initialize our system with a spin-polarized electron density.

In ATK, you can specify the degree of spin-polarization using three different methods:

Your choice of method is specified by setting the electronDensityParameters() similar to what is shown above. We will use a fixed spin for all of our calculations involving oxygen. This will show how different spin values alters the state of the oxygen molecule. Had we used the initial_spin keyword to initialize the electronDensityParameters(), the self-consistent field calculation would change the spin-polarization of, instead of keeping it fixed, to a value optimized in accordance with our convergence criterion. We want to keep the spin-polarization fixed and see how this changes the molecular configuration and total energy. But what does this value represent? How is it kept fixed?

To illustrate the effects of using a spin-polarized electron density, we will calculate the total energy and selected eigenstates of the oxygen molecule for three different states: singlet, triplet, and quintet spin-configuration corresponding to zero, two, and four unpaired electrons in our oxygen molecule. The following NanoLanguage shows how this can be done.

# Create a list of corresponding bond-lengths and
# values of total spin that we want to compare.
spin_types = [
    ('singlet',0.0,1.22),
    ('triplet',1.0,1.22),
    ('quintet',2.0,1.93)] 

vnl_file = VNLFile('oxygen.vnl')
for (spinType,spinValue,ooDistance) in spin_types:
    # Create spin-polarized KohnShamMethod
    dft_oxygen = createKohnShamMethod(spinValue)
    # Create oxygen molecule
    oxygen_molecule = createOxygenMolecule(ooDistance) 
    # Perform actual calculation
    scf_oxygen = dft_oxygen.apply(oxygen_molecule)
    print 'Spin type:',spinType,
    print "Total energy: %8.4f eV" %(calculateTotalEnergy(scf_oxygen).inUnitsOf(Units.eV)) 
    # Calculate eigenstates and export to VNL file.
    # ID reveals the value of the fixed spin.
    homoLumoEigenstates = calculateEigenstates(
        self_consistent_calculation = scf_oxygen, 
        quantum_numbers=([6,7], Spin.Up)
        )
    vnl_file.addToSample(oxygen_molecule, 'oxygen'+spinType)
    for eigenstate in homoLumoEigenstates:
        vnl_file.addToSample(eigenstate, 'oxygen'+spinType)

The variable spin_types is a list, which contains three tuples corresponding to an identifier, a total spin and a bond-distance for the oxygen molecule, respectively. ATK uses atomic units and each unpaired electron therefore has a spin of 0.5 . When we fix the total spin of the molecule to 1.0 we instruct ATK to have two unpaired electrons with spin-up in the system. The output, resulting from the script listed below, shows that the ground state of the oxygen molecule is the triplet state (two spin-up electrons) since both the quintet (four spin-up electrons) and singlet (zero spin-up electrons) state have a higher total energy.

Spin type: singlet Total energy: -864.2015 eV
Spin type: triplet Total energy: -865.3799 eV
Spin type: quintet Total energy: -859.7940 eV

If you look at the tuples describing the three states above, you will notice that we have used a longer bond distance for the quintet state than for the triplet state. We leave it up to you to verify that the bond distance corresponding to the minimum total energy for the quintet state corresponds to 1.93 Angstrom.

The figures to the right show the HOMO of the triplet and the quintet state corresponding to an anti-bonding π*- and σ*-orbital, respectively, with the oxygen atoms shown as wire frames. A closer study of the eigenvalues and corresponding molecular orbitals would reveal that a spin-down electron in a degenerated and bonding π-orbital in the triplet state has been transferred to the anti-bonding σ*-orbital (HOMO) of the quintet state.

This explains the increase in bond length when we fix the spin at 2.0 . By populating the anti-bonding σ*-orbital, with an electron from a bonding π-orbital, we weaken the oxygen bond which in turn increases the bond distance from 1.23 Angstrom to 1.93 Angstrom.

Below you will find the script generating the VNL file used to plot the eigenstates of the oxygen molecule shown in the figures.

from ATK.KohnSham import *

def createOxygenMolecule(atom_separation):
    '''
    @param atom_seperation - bond-length for
           a oxygen molecule
    @return MoleculeConfiguration object for an Oxygen
            molecule.
    '''
    atoms = 2*[Oxygen]

    coordinates =  [
        ( 0.0, 0.0, atom_separation/2) ,
        ( 0.0, 0.0,-atom_separation/2)
        ]*Ang

    oxygen_molecule = MoleculeConfiguration(
        elements = atoms,
        cartesian_coordinates = coordinates
        )
    return oxygen_molecule

def createKohnShamMethod(spin):
    '''
    @param spin - Total spin for the electron density.
           Supplied with no unit.
    '''
    fixed_electron_density = electronDensityParameters(
        fixed_spin = spin*hbar
        )
    
    dft_method=KohnShamMethod(
        basis_set_parameters=basisSetParameters(),
        exchange_correlation_type=GGA.PBE,
        electron_density_parameters = fixed_electron_density
        )
    
    return dft_method

# Create a list of corresponding bond-lengths and
# values of total spin that we want to compare.
spin_types = [
    ('singlet',0.0,1.22),
    ('triplet',1.0,1.22),
    ('quintet',2.0,1.93)] 

vnl_file = VNLFile('oxygen.vnl')
for (spinType,spinValue,ooDistance) in spin_types:
    # Create spin-polarized KohnShamMethod
    dft_oxygen = createKohnShamMethod(spinValue)
    # Create oxygen molecule
    oxygen_molecule = createOxygenMolecule(ooDistance) 
    # Perform actual calculation
    scf_oxygen = dft_oxygen.apply(oxygen_molecule)
    print 'Spin type:',spinType,
    print "Total energy: %8.4f eV" %(calculateTotalEnergy(scf_oxygen).inUnitsOf(Units.eV)) 
    # Calculate eigenstates and export to VNL file.
    # ID reveals the value of the fixed spin.
    homoLumoEigenstates = calculateEigenstates(
        self_consistent_calculation = scf_oxygen, 
        quantum_numbers=([6,7], Spin.Up)
        )
    vnl_file.addToSample(oxygen_molecule, 'oxygen'+spinType)
    for eigenstate in homoLumoEigenstates:
        vnl_file.addToSample(eigenstate, 'oxygen'+spinType)

oxygen-spincomp.py

A different approach in determining the equilibrium state could have been to use the calculateOptimizedAtomicGeometry() method which we also used in the tutorials Geometry of a water molecule and Calculating molecular properties. The strategy will now be to:

We will reuse most of the configurations from the previous NanoLanguage script shown above, however, we will need to change the for-loop over the spin-polarization.

    coord = oxygen.cartesianCoordinates()
    oo_bond_length = coord[1][2] - coord[0][2]
    return math.fabs(oo_bond_length.inUnitsOf(Angstrom))

from numpy import arange
energy_dict ={}
geometry_opt_param = geometricOptimizationParameters()

print 'Separation (Angstrom)\tTotal Energy (eV)\tSpin (hbar)' 
for polarization in  arange(0.0,2.1,0.1):
#for polarization in  arange(0.0,0.1,0.1):
    dft_method = createKohnShamMethod(polarization)

    opt_oxygen = calculateOptimizedAtomicGeometry(
        atomic_configuration=createOxygenMolecule(),
        method=dft_method,
        optimization_parameters=geometry_opt_param
        )
    
    distance = determineBondLength(opt_oxygen)

The above script will calculate an optimized atomic geometry for a range of different degrees of polarization. We then calculate the total energy of this optimized geometry and output the result. In addition, we have created a function for finding the bond distance, determineBondLength() which is specified in the following way:

        )

    return dft_method

def determineBondLength(oxygen):

Here, we return the absolute value of the bond difference, using the fabs() function from the math module of Python. From the figure to the right, we see that the equilibrium bond distance increases with increasing spin-polarization. This is also what we saw before and therefore in agreement with our previous results.

The drawback of the above method using geometry optimization is that it performs two SCF-calculations per atomic configuration of the oxygen molecule: once, when optimizing and once because we need access to a SelfConsistentCalculation object to determine the total energy. On the other hand, this method requires no prior knowledge of approximate bond distances or spin states as opposed to what was used in the first example.

The output from the optimization script and the complete script for performing this alternative approach is given below.

from ATK.KohnSham import *

def createOxygenMolecule( atom_separation = 1.23):

    atoms = 2*[Oxygen]

    coordinates =  [
        ( 0.0, 0.0, atom_separation/2) ,
        ( 0.0, 0.0,-atom_separation/2)
        ]*Ang
    
    oxygen_molecule = MoleculeConfiguration(
        elements = atoms,
        cartesian_coordinates = coordinates
        )

    return oxygen_molecule

def createKohnShamMethod( spin = 0.0):

    total_electron_density = electronDensityParameters(
        fixed_spin = spin*hbar
        )

    iteration_control_parameters = iterationControlParameters(tolerance = 1e-6)
    
    dft_method=KohnShamMethod(
        basis_set_parameters=basisSetParameters(),
        exchange_correlation_type=GGA.PBE,
        electron_density_parameters = total_electron_density,
    iteration_control_parameters = iteration_control_parameters
        )

    return dft_method

def determineBondLength(oxygen):
    import math
    coord = oxygen.cartesianCoordinates()
    oo_bond_length = coord[1][2] - coord[0][2]
    return math.fabs(oo_bond_length.inUnitsOf(Angstrom))

from numpy import arange
energy_dict ={}
geometry_opt_param = geometricOptimizationParameters()

print 'Separation (Angstrom)\tTotal Energy (eV)\tSpin (hbar)' 
for polarization in  arange(0.0,2.1,0.1):
#for polarization in  arange(0.0,0.1,0.1):
    dft_method = createKohnShamMethod(polarization)

    opt_oxygen = calculateOptimizedAtomicGeometry(
        atomic_configuration=createOxygenMolecule(),
        method=dft_method,
        optimization_parameters=geometry_opt_param
        )
    
    distance = determineBondLength(opt_oxygen)
    scf_oxygen = dft_method.apply(createOxygenMolecule(distance))
    total_energy = calculateTotalEnergy(scf_oxygen).inUnitsOf(eV)
    print "%8.4f %s %8.4f %s %2.1f" %(distance,'\t\t',
                                      calculateTotalEnergy(scf_oxygen).inUnitsOf(eV),\
                                      '\t\t',polarization)
    energy_dict[total_energy] = (distance,polarization)

energy_list = energy_dict.keys()
energy_list.sort()
print "Minimum total energy: %8.4f eV" % (energy_list[0])
print "O-O equilibrium distance: %8.4f A" % (energy_dict[energy_list[0]][0])

oxygen-spinopt.py

In this example we have

If you want to see more examples of configuring a spin-polarized we configure a spin-polarized bulk system.