Geometry of a water molecule

	Tutorials
Tutorials	Calculating molecular properties

Introduction

In this example, you will be guided through some of the most basic features of ATK. In summary, you will learn how to:

define the geometric configuration for a water molecule
visualize the molecule in 3D
obtain physical information about the defined molecule

Prerequisites

Before you start going through this tutorial, we recommend that the following requirements are met:

You should be familiar with using a text editor to open, edit and save plain text files.
A working installation of ATK already is available on your system. If this is not the case, please contact your local system administrator or consult the installation guide.

We emphasize that this tutorial requires no prior knowledge about either ATK, NanoLanguage, or Python. As we go along, we will introduce some simple ATK/NanoLanguage concepts, discuss their specific usage, and then gradually piece these together to some real applications.

In order to do this, we will also need to introduce a few basic Python concepts, but we will explain their functionality and meaning along the way.

Whenever we mention elements from either NanoLanguage or Python, links pointing to further information is given directly in the text: If you click these, you will be directed to a page which provides further details, explanations, and usage examples for a given concept.

If you are a complete Python novice, we suggest that you follow the links whenever you encounter a new Python concept, and then return to the tutorial text once you feel updated.

Should you need to learn a little more Python before you move on, we have written a very short Python primer, Python basics, which provides you with a brief description of the most common concepts used in Python.

Preliminaries

Before we can start to use ATK, we need a precise definition of the physical system that we want to work with. For simple molecules, this boils down to specifying:

the atoms from which the molecule is composed
the positions of the individual atoms in the molecule

In our case, the system is a water molecule, which we have sketched in the small diagram shown Figure 9.

$A water molecule. The geometry of the molecule is determined by the respective coordinates 16 6 4cffd93a2a92f0bd0fbd840c5a79bdec (x,y,z)_{\mathrm{H}_1} , 16 5 cca541f9ab5d146fa8fca482a7bbecef (x,y,z)_{\mathrm{H}_2} , and 16 5 f243ff9c683c45593e4b16fe87eabdc8 (x,y,z)_{\mathrm{O}} of the two hydrogen and the oxygen atoms.$

Figure 9: A water molecule. The geometry of the molecule is determined by the respective coordinates $(x,y,z)_{\mathrm{H}_1}$ , $(x,y,z)_{\mathrm{H}_2}$ , and $(x,y,z)_{\mathrm{O}}$ of the two hydrogen and the oxygen atoms.

As illustrated in Figure 9, a specification of our system could be:

The molecule consists of two hydrogen atoms and one oxygen atom (indicated as $\mathrm{H}_1$ , $\mathrm{H}_2$ , and $\mathrm{O}$ on the figure).
The position in three dimensional space of the individual atoms is given by the three coordinate sets $(x,y,z)_{\mathrm{H}_1}$ , $(x,y,z)_{\mathrm{H}_2}$ , and $(x,y,z)_{\mathrm{O}}$ .

We now have a concise definition of both the molecule constituents and its overall geometry, and should be ready to define and build the molecule in our first ATK script.

Building the molecule

Now, start your favorite text editor and open an empty file named water.py. Add the following line to the top of the file and save it to the disk

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

This single line is a Python import statement that will instruct ATK to load the ATK.KohnSham module. The contents of the ATK.KohnSham module will give us access to a large set of ATK objects and functions for doing both simple and advanced DFT calculations. Here, we will focus on the ones needed for building the water molecule.

You may run the script by issuing the following command from the command line

atk water.py

Clearly, not much happens so far. The script instructs ATK to load the ATK.KohnSham module, exits from ATK and returns to the command line.

Let us provide some real content to the script related to our molecule by adding the additional lines

# Set up elements and positions
elm = [ Oxygen, Hydrogen, Hydrogen ]         
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]

These two new statements define two variables: elm and pos. Both of them are of the Python type list. The variable elm contains a list of the atoms of the water molecule, whereas the variable pos holds the three dimensional Cartesian positions of the corresponding atoms.

elm contains entries of the type PeriodicTable.Element. Note that you have to specify the full name of a given atom, i.e. you cannot use “H” or “O” as abbreviations for “Hydrogen” or “Oxygen”, respectively. All the elements from the periodic table are available in ATK. See ATK.PeriodicTable for a complete list.
pos contains a list of the positions of the atoms. In ATK, the position of each atom consists of three objects of the type PhysicalQuantity. The first object corresponds to the position along the x-axis, the second is the position along the y-axis and so on. A physical quantity is created by taking a number and multiplying by a unit. In this case, we have a good guess of the positions of the atoms in terms of Angstrom, so we use this guess and multiply by Angstrom. If we had a good idea of what the position should be in Bohr, we would use that number and multiply by Bohr.
Note that the order of the elements in the two lists are important. The first element in pos should correspond to the oxygen atom, since this is the first element in elm.

In ATK, all variables and parameters that correspond to a physical quantity must be specified as such. That is, it should be given in terms of a number multiplied by a mandatory unit. A broad range of relevant units are available from the ATK.Units module. However, most of them are also included in the ATK.KohnSham and ATK.TwoProbe modules.

A more detailed discussion on units and PhysicalQuantity objects can be found in the reference manual

Now, just as before, you may run the script by issuing the command

atk water.py

in which case the following set of actions takes place

atk loads the ATK.KohnSham module
the variables elm and pos are defined
atk exits, and returns to the command line

Now that we have set up definitions for the elements of the water molecule as well as the molecule geometry, we are ready to pass this information on to ATK in order to 'build' the actual molecule.

To do this, add the following line to the end of script

# Add them to a configuration
h2o = MoleculeConfiguration(elm,pos)

The last line in the script defines a new variable h2o, which becomes an object of the type MoleculeConfiguration. Since both the elements and their positions are passed on to the object at its creation, all necessary information about the water molecule is now stored in the h2o object. In this sense, the water molecule has been built and is now represented on the computer as an ATK object.

The entire script you have typed should now look like this

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

# Set up elements and positions
elm = [ Oxygen, Hydrogen, Hydrogen ]         
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]      

# Add them to a configuration
h2o = MoleculeConfiguration(elm,pos)

water.py

If you wish to download a copy of the script, click the link shown directly underneath the above script listing.

Observe also, how we have added comments to the various parts of the script. You do this by adding the symbol “#” in front of a line, which instructs ATK to ignore the entire line. Always add comments to your scripts: It greatly improves both readability and structure, and most important, it becomes much easier for other people to see what you have been doing.

As you will see in the following sections, the h2o object contains several features that will enable us to do both simple and complex manipulations with the water molecule.

Visualizing the molecule

Looking at a water molecule in terms of the two lists

elm = [ Oxygen, Hydrogen, Hydrogen ]
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]

provides a somewhat non-intuitive view of the molecule's geometry. This trend becomes increasingly dominant as the size of the molecule increases. Furthermore, when when you start using ATK for building larger molecules, there will be a much higher risk of making a mistake or typing error in one of the above lists. In addition, these errors will typically be very difficult to hunt down.

For this reason, it is always good practice to inspect your molecule with an external 3D viewer before you proceed any further. ATK provides this functionality via its graphical counterpart Virtual NanoLab (VNL). Exporting a molecule to VNL is simple and consists of two steps:

exporting the molecule to a VNLFile
importing the VNLFile into VNL

Let us see how these two steps are carried out in practice:

First, we have to export the molecule to a VNLFile. You do this by adding the following lines to your script

vnl_file = VNLFile('h2o.vnl')
vnl_file.addToSample(h2o, 'h2o')

If you run the script, the following takes place when the two last lines are processed by ATK

A variable called vnl_file is defined which is an object of the type VNLFile. The argument h2o.vnl to VNLFile is a string containing the name of the file in which we store information about our water molecule.
To store the information in the VNLFile, we call the method addToSample(), which takes two arguments: The MoleculeConfiguration object assigned to the variable h2o and a text string h2o associated with the object we store. Since you may store several configurations in a single VNLFile by using the method addToSample(), the text string can later be used to retrieve or load a specific configuration from the VNLFile.

The entire script for exporting the VNLFile should look like this

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

# Set up elements and positions
elm = [ Oxygen, Hydrogen, Hydrogen ]
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]

# Add them to a configuration
h2o = MoleculeConfiguration(elm,pos)

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

water-viz.py

Now, save the script and process it through ATK. Once the script has been processed, verify that a VNLFile named h2o.vnl is available on your disk.

In order to view the molecule in VNL, you should first start VNL and then open the VNLFile file h2o.vnl using the result browser.. A picture showing a 3D plot of the water molecule is shown in Figure 10

	Warning
	If you rerun the script `water-viz.py`, ATK will overwrite the existing sample called `h2o`, without issuing an error.

So in general, if you want to run scripts without overwriting samples, you either have to rename your sample, for example

vnl_file.addToSample(h2o, 'h2o-new')

or rename or copy your VNLFile before the next run. You can find a possible solution (one of many) for automatically back up an existing VNL file in the Tips and Tricks appendix of this manual.

3D plot of a water molecule rendered in VNL. The plot was generated by exporting the molecule data from an ATK script and then importing the data into VNL.

Figure 10: 3D plot of a water molecule rendered in VNL. The plot was generated by exporting the molecule data from an ATK script and then importing the data into VNL.

Querying the molecule

Once we have created a physical system in ATK, there are numerous things that we can do with the structure. In particular, we may

ask questions about the current state of the system
perform operations to alter the state of the system

Let us have a look at some simple operations of the first type using our water molecule as an example. Suppose that we did not know the elements present in the molecule, this information could be obtained using the following construction

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

# Set up elements and positions
elm = [ Oxygen, Hydrogen, Hydrogen ]
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]

# Add them to a configuration
h2o = MoleculeConfiguration(elm,pos)

# Print the elements in the molecule
for e in h2o.elements():
    print e.symbol()

water-elm.py

Running this script produces the following output

O
H
H

telling us that the molecule consists of the three atoms O, H, and H (in this case, we knew that already, of course). In summary, the following operations are performed by the script

First, we call the method elements() of the MoleculeConfiguration object which was assigned to the variable h2o. This returns a list of the elements contained in the molecule.
The variable e is used in a for-loop construction to address each member of the list.
The variable e is an Element object. For each element in the list, the method symbol() of the Element object is used to print out the atomic symbol for each element in the molecule.

Another relevant piece of information would be to obtain the physical positions of the oxygen and the two hydrogen atoms in the molecule. This can be achieved by using a construction very similar to the previous example, namely

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

# Set up elements and positions
elm = [ Oxygen, Hydrogen, Hydrogen ]
pos = [ ( 0.000, 0.000, 0.0)*Angstrom,
        ( 0.757, 0.586, 0.0)*Angstrom,
        (-0.757, 0.586, 0.0)*Angstrom ]

# Add them to a configuration
h2o = MoleculeConfiguration(elm,pos)

# Print the coordinates of the individual atoms
for coord in h2o.cartesianCoordinates():
    for c in coord:
        print c.inUnitsOf(Angstrom),
    print

water-cor.py

which generates the output

0.0 0.0 0.0
0.757 0.586 0.0
-0.757 0.586 0.0

in accordance with our original specification of the positions for the atoms in the water molecule. In this case, the following actions take place

The method cartesianCoordinates() of the MoleculeConfiguration object is called. This returns a list containing tuples of coordinates for each element in the molecule.
We access each tuple using a for-loop. In each iteration of the for loop, the tuple is assigned to the variable coord.
Each of these tuples contain the x-, y-, and z-coordinates of the atom position. Each coordinate is represented in the tuple as a PhysicalQuantity object.
We iterate through the tuple elements using an additional for loop; In each loop iteration, the tuple elements are assigned to the variable c. Finally, we print out the coordinates using the method inUnitsOf() of the PhysicalQuantity object.

Finding an energetic minimum

Now that we know how to set up a molecular structure, we are ready to do some actual calculations on our water molecule. You have probably already noticed that the way we specify the geometry of the water molecule contains some redundant elements. Providing coordinates for all three atoms is not necessary, as is shown in Figure 11.

$A simpler representation of the molecule geometry. The geometry is determined by the bond length 14 3 aca67a1a450a55ae646afd1a32844f23 d_\mathrm{OH} of the O-H bond and the bending angle 14 0 5b11e05f05dbbee9ab42e5ffa80feff6 \theta between the three atoms H-O-H.$

Figure 11: A simpler representation of the molecule geometry. The geometry is determined by the bond length $d_\mathrm{OH}$ of the O-H bond and the bending angle $\theta$ between the three atoms H-O-H.

As you can see, the geometry is completely defined, if we just specify the length $d_\mathrm{OH}$ of the O-H bond and the bending angle $\theta$ between the three atoms H-O-H. A simple Python function to automate this process of 'converting' the two numbers $d_\mathrm{OH}$ and $\theta$ to the respective Cartesian coordinates of the individual atoms could look like this

# Define function for handling molecule setup   
def waterConfiguration(angle,bondLength):
    from math import sin,cos
    
    theta = angle.inUnitsOf(radians)
    pos = [(0.0, 0.0, 0.0)*Angstrom,
           (1.0, 0.0, 0.0)*bondLength,
           (cos(theta), sin(theta), 0.0)*bondLength]
    elm = [Oxygen] + [Hydrogen]*2
    
    return MoleculeConfiguration(elm,pos)

Here, we have used one of the advantages of ATK, namely that we always have the power of Python behind us: The geometric layout of the molecule is now taken care of by the user defined Python function waterConfiguration() - all we have to do now is to pass the bond length $d_\mathrm{OH}$ and the bending angle $\theta$ to waterConfiguration(), which then returns a MoleculeConfiguration with everything set up correctly. This approach greatly simplifies matters for several reasons

We reduce the probability of errors, since the function only has to be entered once.
The Python function can be reused to generate as many different MoleculeConfiguration objects as you like.
Trivial code segments are removed from the main part of the script, increasing both its readability and focus.

Calculating the energetic minimum

The experimental values for the equilibrium bond length $d^\mathrm{eq}_\mathrm{OH}$ and the bending angle $\theta^\mathrm{eq}_\mathrm{ }$ are reported in the NIST Chemistry WebBook as

$\displaystyle d^\mathrm{eq}_\mathrm{OH} = 0.958\; \text{Å} \quad\text{and}\quad \theta^\mathrm{eq} = 104.5^\circ$

Let us use the function waterConfiguration() to scan a range of bending angles, and then let ATK check that the reported equilibrium angle $\theta = 104.5^\circ$ in fact corresponds to an energetic minimum of the molecule. A simple script to carry out this task would be

# 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)
    pos = [(0.0, 0.0, 0.0)*Angstrom,
           (1.0, 0.0, 0.0)*bondLength,
           (cos(theta), sin(theta), 0.0)*bondLength]
    elm = [Oxygen] + [Hydrogen]*2
    
    return MoleculeConfiguration(elm,pos)

   
# Choose DFT method
method = KohnShamMethod()

# Scan different bending angles and calculate energy
for i in range(30,181,5):
    theta = i*degrees
    h2o = waterConfiguration(theta, 0.958*Angstrom)
    scf = method.apply( h2o )
    print 'Angle = ', theta, ' Total Energy = ',calculateTotalEnergy(scf)

water-min.py

In this case, we

set the bond length to the experimental value $d^\mathrm{eq}_\mathrm{OH} = 0.958\; \text{Å}$ ,
scan the range of angles $\theta \in [30:180]$ incrementing the angle by 5 after each loop iteration,
for each angle $\theta$ , we set up a molecule with this particular bending angle,
perform a DFT calculation,
and finally calculate the total energy of the molecule.

As you can see in Figure 12, the calculation predicts an energetic minimum which corresponds nicely with the reference angle $\theta^\mathrm{eq} = 104.5$ .

$The variation of the total energy of the water molecule as a function of the bending angle 14 0 5b11e05f05dbbee9ab42e5ffa80feff6 \theta . The energy variation shows a minimum close to the experimental minimum 14 1 9da0732c8858a141be7d8b7dca09e39f \theta^\mathrm{eq} = 104.5^\circ .$

Figure 12: The variation of the total energy of the water molecule as a function of the bending angle $\theta$ . The energy variation shows a minimum close to the experimental minimum $\theta^\mathrm{eq} = 104.5^\circ$ .

	Note
	The methods employed here for determining the equilibrium state of the water molecule were introduced for pedagogical reasons. In fact, the equilibrium state can be determined much more effectively by using the relaxation tools provided by ATK, as shown in the next section.

Relaxing the molecule

Suppose that we did not know the precise equilibrium value of either the O-H bond length $d^\mathrm{eq}_\mathrm{OH}$ or the bending angle $\theta^\mathrm{eq}$ . These two quantities can also be obtained by ATK by performing a geometric optimization of the water molecule.

To do this, we set up a water molecule with slightly wrong initial values of both $d_\mathrm{OH}$ and $\theta$ . We then let ATK perform a DFT calculation to find the equilibrium structure of the molecule. Here is a script that solves the task

# 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)
    pos = [(0.0, 0.0, 0.0)*Angstrom,
           (1.0, 0.0, 0.0)*bondLength,
           (cos(theta), sin(theta), 0.0)*bondLength]
    elm = [Oxygen] + [Hydrogen]*2

    return MoleculeConfiguration(elm,pos)


# 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)/( length(oh1)*length(oh2) )

    return acos(prod)*radians


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

    return length(oh1)


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

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

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


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

# Set DFT method with a low tolerance
low_tolerance = iterationControlParameters( tolerance = 1.0e-6 )
dft_method = KohnShamMethod(iteration_control_parameters = low_tolerance)

# Calculate the optimized geometry for H2O molecule.
opt_params = geometricOptimizationParameters(
    force_tolerance = 0.005*eV/Ang
    )
h2o_opt = calculateOptimizedAtomicGeometry(
    atomic_configuration = h2o,
    optimization_parameters=opt_params,
    method = dft_method
    )

# Calculate optimized bond length and bending angle
bondlength = waterBondLength(h2o_opt)
bendingangle = waterBendingAngle(h2o_opt)

# Print the results
print 'H-O-H angle    = ', bendingangle.inUnitsOf(degrees)
print 'O-H bondlength = ', bondlength

water-opt.py

We will not go into detail with the objects and functions:

- instead we refer to the tutorial Calculating molecular properties, where their use is explained in detail. As you may have noticed, we have also included some additional user defined functions waterBendingAngle() and waterBondLength() for calculating angles and bond lengths.

	Note
	In the script for the geometry optimization of a water molecule, the tolerance for the self-consistent calculations has been lowered to the value of $10^{-6}$ in order to achieve low values of the force acting on the single atoms.

Now, if you run the script, ATK will output the following results

$\displaystyle d_\mathrm{OH} = 0.975\; \text{Å} \quad\text{and}\quad \theta= 103.2^\circ$

in good correspondence with the reported experimental values.

Summing up

We hope that you have enjoyed your first encounter with ATK. At this point, you should be able to

set up simple molecular structures,
visualize your structure in a 3D plot,
obtain information about the geometry and nature of the atoms in your molecule,
and perform simple energy calculations and geometric optimizations on your structure.

In this simple walk-through, we have tried to emphasize some of the most basic and fundamental aspects needed when setting up an ATK script. ATK scripts can tackle much more challenging aspects than just setting up and manipulating the geometry of a water molecule. And combined with the power of Python, complex and involved problems can often be solved using very simple and elegant approaches.

Hopefully, we have managed to whet your appetite. In this case, we suggest that you move on to one of the more complicated tutorials:

Calculating molecular properties, where the calculation of more advanced molecular properties are explained.
The band structure of bulk silicon, where ATK is used to build and simulate a bulk silicon crystal.
Li-H₂-Li two-probe system, where we show how to construct a two-probe system.

Of course, you can also go ahead and define your own physical system and start solving and analyzing it using ATK.

Whatever your approach will be, we hope you will enjoy using ATK!