Visualisation

Scientific visualisation is the representation of data graphically so that you can gain an understanding and insight into the data. It is useful in virtually every domain of science. Among the first courses where 3D visualisation helps is in gaining an insight into bonds in chemistry.

The simplest molecule would be a hydrogen molecule ion. You can approximate its orbital functions as a linear combination of atomic orbitals. The intention in this article is to take a very simple item and explore the various ways of visualisation.

2d Plots

Let the two atoms be at +/-a on the x-axis. The simplest two orbitals will be the bonding and anti-bonding 1s-orbitals. Ignoring the normalisation factors, you can plot the electron density on the x axis as follows.

import scipy as np

a = 0.7 # half the distance between atoms

a0 = 0.53 # Bohr radius in Angstrom

x = np.mgrid[-5:5:100j]

sa = np.exp(-np.sqrt((x - a)**2/a0))

sb = np.exp(-np.sqrt((x + a)**2/a0))

import matplotlib.pyplot as plt

plt.subplot(2,1,1)

plt.plot(x,sa + sb, '-', label='Bond')

plt.plot(x,sa , '-', label='A')

plt.plot(x,sb , '-', label='B')

plt.legend()

plt.subplot(2,1,2)

plt.plot(x, np.absolute(sa - sb), '-', label='Anti-bond')

plt.plot(x,sa , '-', label='A')

plt.plot(x,sb , '-', label='B')

plt.legend()

plt.show()

You have used mgrid instead of arange to create an array. The syntax has a hack that if the 3rd parameter is complex, it is interpreted as the number of points rather than the step size.

You have used the Matplotlib and it gives a static 2-d image as you can see in Figure 1. If you wanted to explore the function in other directions, you would have to create a separate plot for each them.



3 Dimensional Visualisation and Mayavi

What if you wish to explore the electron density in a plane? Mayavi comes to your aid. It is a project initiated by Prabhu Ramachandran and currently enhanced and supported by http://www.enthought.com/. Mayavi offers very easy to use Python and Scipy bindings to the vistualisation toolkit, VTK, http://www.vtk.org/. You will need to install Mayavi package. It is available in the Fedora and Ubuntu repositories.

You can plot the function values of the bonding orbital in the xy plane as follows:

import scipy as np

x,y = np.mgrid[-3:3:100j, -3:3:100j]

a = 0.7

a0 = 0.53

r1=np.sqrt((x-a)**2 + y**2)

r2= np.sqrt((x+a)**2 + y**2)

values = np.exp(-r1/a0) + np.exp(-r2/a0)

from enthought.mayavi import mlab

mlab.surf(values, warp_scale='auto')

mlab.outline()

mlab.axes()

mlab.show()

The mlab module in mayavi interfaces to the scipy arrays. You create a two dimensional array using mgrid. The surf method plots the bonding function values against the x,y grid. The outline and axes are shown. Finally, the show method ensures that the image remains visible till you close the window.

You can interact with the scene as illustrated by two images in Figure 2.

  


A molecule is a 3 dimensional object. So, what if you wanted to explore the orbital function in x, y and z axes as well? One solution is to observe the contours.

import scipy as np

x,y,z = np.mgrid[-3:3:100j, -3:3:100j, -3:3:100j]

a = 0.7

a0 = 0.53

r1=np.sqrt((x-a)**2 + y**2 + z**2)

r2= np.sqrt((x+a)**2 + y**2 + z**2)

values = np.exp(-r1/a0) + np.exp(-r2/a0)

from enthought.mayavi import mlab

mlab.contour3d(values,opacity=0.5)

mlab.outline()

mlab.axes()

mlab.show()

You have now created a 3-dimensional grid. The values of the electron density at each point is given in values. The method contour3d plots the contours and you can interact with the image and view it from various perspectives. The opacity parameter is very important as it allows you to view the contours surrounded by other contours. One perspective is shown in Figure 3.


While contours are useful, you may wish to examine the values at each location. This is possible in Mayavi as in the following example.

import numpy as np

x,y,z = np.mgrid[-3:3:100j, -3:3:100j, -3:3:100j]

a = 0.7

a0 = 0.53

r1=np.sqrt((x-a)**2 + y**2 + z**2)

r2= np.sqrt((x+a)**2 + y**2 + z**2)

values = np.exp(-r1/a0) + np.exp(-r2/a0)

from enthought.mayavi import mlab

source = mlab.pipeline.scalar_field(values)

min = values.min()

max=values.max()

mlab.pipeline.volume(source,vmin=min, vmax=min + .5*(max-min))

mlab.show()

You create a source for a scalar field and the volume method maps the source values to colours. The range of values to be mapped is given by the vmin and vmax parameters. A result can be seen in Figure 4.



The power of Mayavi/VTK is that with just a few lines of code, you can create beautiful imagery which can help get an insight into the scientific data. Scipy/Mayavi combination allows you to focus on the data and the observations.

You can explore more examples of what is possible at http://code.enthought.com/projects/mayavi/docs/development/html/mayavi/auto/examples.html.

Comments