Introductory Tutorial

This tutorial assumes you’ve successfully installed quippy and that all the tests pass. If not, see Installation of QUIP and quippy. We make use matplotlib, scipy and the quippy atomeye extension module in this tutorial, but you can skip other these sections if you don’t have these components installed.

Other tutorials that may be useful:

  1. If you’re new to Python, you might like to start with the offical Python tutorial.
  2. quippy depends on and is compatible with the Atomic Simulation Environnment (ASE), so their tutorials are a useful resource. In particular, there’s a good introductory Python tutorial as part of the ASE documentation.
  3. The numpy project also provides an (unfinished) tutorial. If you follow this tutorial, be aware that quippy’s FortranArray class uses one- rather than zero-based indexing, in order to fit in better with Fortran. numpy also has an online reference manual
  4. Numpy for Matlab Users may be useful if you’re familiar with MatLab.

There’s a couple of different ways you can use quippy: either by typing commands interactively in Python (or, better, ipython), or by writing a script. It’s easier to play around with things interactively, which is what we’ll do in this tutorial, but when you want to do something more than once it’s worth saving the commands in a script file. Just create a text file with the extension .py.

Getting Started

Let’s start by firing up Python and importing quippy:

$ python
>>> from qlab import *

The >>> is the Python prompt, at which you type commands. Here we’ve asked Python to import all the names defined in qlab, which is a convenience module that gives us access to all the quippy functionality.

Alternatively, if you have the ipython improved interactive shell installed, you can start quippy and load all the required modules with a helper script:

$ quippy

Now, let’s create an Atoms object. This is the fundamental class in quippy, and it represents a collection of atoms which together make up a crystal or molecule. We’ll make an 8-atom silicon unit cell:

>>> dia = diamond(5.44, 14)

The first argument to diamond() is the lattice constant in Angstrom, and the second is the atomic number. Let’s poke around a little inside our new Atoms object. Here are two equivalent ways to determine the number of atoms:

>>> print len(dia)
8
>>> print dia.n
8

We can save our new atoms to an XYZ file, and then re-read it by passing the filename to the Atoms constructor. Let’s do that, and check that what we get back is the same as the original:

>>> dia.write('diamond.xyz')
>>> dia2 = Atoms('diamond.xyz')
>>> print dia == dia2
True
>>> print dia is dia2
False

Note that dia == dia2 evaluates to True since dia2 is an identical copy of dia, but dia is dia2 is False since there are two separate Atoms objects, residing at different memory locations.

Warning

Be aware that the assignment operator (i.e. the = sign) in Python does not behave as you may expect from other programming languages: assignment always creates references not copies.

In Python, variables are names for values, and assignment binds a name to a value. For example the sequence of commands

>>> a = 1
>>> b = a
>>> print a
1
>>> print b
1
>>> print a == b
True
>>> print a is b
True

creates one value (1) with two names (a and b). The comparison operator == checks whether two names corresond to the same value, wheras the is operator checks if they correspond to the same address. The expression a is b is equivalent to id(a) == id(b), where id() is a Python built in function which returns a unique number identifying a the memory address to which a reference points.

If you forget about this you can get into difficulties when dealing with mutable objects such as lists or Atoms objects:

>>> a = []
>>> b = a
>>> print a is b
True
>>> a.append(1)
>>> print a
[1]
>>> print b
[1]

Here a and b are two names for the same list, so appending a value to a modifies b as well. If this is not what you want, create a copy: for list objects this is done with the slice operator:

>>> a = []
>>> b = a[:]
>>> print a is b
False
>>> a.append(1)
>>> print a
[1]
>>> print b
[]

To make a copy of a more complex object like an array or an Atoms instance you should use the copy() method, e.g.

>>> a = diamond(5.44, 14)
>>> b = a.copy()
>>> print a == b
True
>>> print a is b
False

Instead, use the copy() method, like in the second statement above.

For more information, consult the offical Python tutorial.

Various different file formats are supported, see Supported File Formats.

If you have the optional atomeye module installed, you can visualise your new Atoms object by creating an AtomsViewer:

>>> view(dia)
_images/si8.png

An AtomEye window should pop up with a 3D representation of the silicon crystal. Right clicking on an atom prints information about it:

Atom 3
------

pos             =  [ 2.72  2.72  0.  ]
species         =  Si
z               =  14

The positions and atomic numbers of each atom are available in the pos and z properties:

>>> print dia.z
[14 14 14 14 14 14 14 14]
>>> print dia.pos.T
[[ 0.    0.    0.  ]
 [ 1.36  1.36  1.36]
 [ 2.72  2.72  0.  ]
 [ 4.08  4.08  1.36]
 [ 2.72  0.    2.72]
 [ 4.08  1.36  4.08]
 [ 0.    2.72  2.72]
 [ 1.36  4.08  4.08]]

Note that we print the transpose() of the position array using the T attribute, since the actual dia.pos array has 3 rows and len(dia) == 8 columns, which isn’t very convenient for printing. Atomic properties are stored in arrays (actually, in a special type of array called a FortranArray), so it’s easy to access parts of the data. Array indexing and slicing works in the same way as in Fortran, except that square brackets are used:

>>> print dia.pos[1,1]   # x coordinate of atom 1
0.0
>>> print dia.pos[:,1]   # position of atom 1
[ 0.  0.  0.]
>>> print dia.pos[1]     # alternative syntax for position of atom 1
[ 0.  0.  0.]
>>> print dia.pos[1,:]   # all the x coordinates
[ 0.    1.36  2.72  4.08  2.72  4.08  0.    1.36]
>>> print dia.z[1:3]        # atomic numbers of atoms 1 to 3 inclusive
[14 14 14]

You can also do fancier indexing as we’ll see below.

Manipulating Atoms

Let’s make a more complex structure to play with, a supercell() of 3 x 3 x 3 alpha_quartz() unit cells:

>>> unit = alpha_quartz(4.92, 5.40, 0.4697, 0.4135, 0.2669, 0.1191)
>>> aq = supercell(unit, 3, 3, 3)
>>> view(aq) # only do this if AtomEye plugin is installed
_images/alphaquartz.png

This cell contains 243 atoms. Let’s look at the lattice vectors, which are stored in a 3 x 3 matrix with a = aq.lattice[:,1], b = aq.lattice[:,2] and c = aq.lattice[:,3]:

>>> print aq.lattice
[[  7.38         7.38         0.        ]
 [-12.78253496  12.78253496   0.        ]
 [  0.           0.          16.2       ]]

We can convert from a cartesian representation to cell lengths and angles to confirm that the cell has trigonal symmetry, using a function called get_lattice_params(). To get online help on any anything from within Python, you type help(name). In ipython, it’s even easier, just postfix the name with a question mark (or two question marks for even more information). Here’s the help for /get_lattice_params():

>>> help(get_lattice_params)
get_lattice_params(lattice)
    Wrapper around Fortran :func:`get_lattice_params_`

    Returns parameters of `lattice` as 6-tuple (a,b,c,alpha,beta,gamma).

From this we can see that this function takes a 3x3 matrix lattice and returns six numbers: the three cell lengths and three cell angles. Let’s call it with the lattice associated with our new alpha-quartz Atoms object:

>>> a, b, c, alpha, beta, gamma =  get_lattice_params(aq.lattice)
>>> print a, b, c
14.76 14.76 16.2
>>> print degrees((alpha, beta, gamma))
90.0 90.0 120.0

So we have \(a = b \ne c\), and \(\alpha = \beta \ne \gamma\), i.e. trigonal symmetry.

Let’s try out some of the more advanced indexing operations. Firstly, logical operations on arrays return arrays, which can be used for indexing:

>>> print (aq.z == 8).sum()      # number of Oxygen atoms
162
>>> print (aq.z == 14).sum()     # number of Silicon atoms
81
>>> aq.pos[aq.z == 8]      # positions of the Oxygen atoms
...

Secondly, it’s also possible to specify a list of indices, and to use negative numbers to count backwards from the end of an array:

>>> print aq.pos[:, [1,-1]].T   # positions of first and last atoms
[[ 1.155462   -2.00131889  3.6       ]
 [ 9.544062   -1.7618594   8.35686   ]]

You can also assign to sections of an array using fancy indexing. For example, the following codes adds a charge property to our alpha quartz object and sets the charge of the Oxygen atoms to -1.4 and that of the Silicon atoms to +2.8

aq.add_property('charge', 0.0)
aq.charge[aq.z == 8]  = -1.4
aq.charge[aq.z == 14] =  2.8

If you have atomeye installed, you can colour the atoms by the new charge property using:

aux_property_coloring('charge')

Here’s an example of making a new Atoms object containing only the atoms with positive y coordinates. We first ensure all atoms are within the first unit cell using map_into_cell():

>>> aq.map_into_cell()
>>> aq2 = aq.select(aq.pos[2,:] > 0)
>>> view(aq2) # only if AtomEye is available
_images/alphaquartz2.png

Note

AtomEye has a different idea of the first unit-cell than quippy. I’ve shifted the cell by lattice coordinates of (0.5, 0.5, 0.5) to make this image. With the extended version of AtomEye included in quippy, you can do this by pressing Shift+z.

Using an AtomsList to analyse results

One common use of quippy is for post-processing the results of simulations conducted with QUIP and libAtoms. It can read and write Atoms objects to and from a variety of file formats, with extended XYZ and NetCDF being the natively supported formats. Let’s see how we might construct a simple script to load and analyse some simulation data.

The first example is si-1000.xyz, which contains snapshots from a molecular dynamics simulation of silicon with the Stillinger-Weber potential, starting from a randomised velocity distribution at a temperature of 1000K (see Potentials and DynamicalSystems example).

The AtomsList class is used to deal with a series of configurations. We can open the XYZ file and see how many frames it contains like this:

>>> al = AtomsList('si-1000.xyz')
>>> print len(al)
100

If you have the AtomEye extension, you can open a viewer with the view() function at this point, to allow visualisation of the trajectory: use Insert and Delete to move forwards and backwards:

>>> view(al)

Note

The previous two steps can be accompished in one go using the quippy script from the shell command prompt:

$ quippy si-1000.xyz

In this case the trajectory is loaded directly into an AtomsListViewer class, with a name taken from the basname of the input filename, e.g. si_1000 in this case. You can check the name of the current viewer with the gcv() function:

>>> print gcv().name
'si_1000'
>>> print len(si_1000)

For the purposes of following along with this tutorial, it’s simplest to create the al AtomsList as described above.

An AtomsList is essentially nothing more than a list of Atoms objects. You can access the component objects by indexing, e.g. al[i] is the ith Atoms object within al.

Let’s check the conservation of the total energy in this simulation. The first frame in the data we’ve loaded has the following properties and parameters:

>>> print al[0].properties.keys()
['species', 'pos', 'Z', 'travel', 'mass', 'move_mask', 'damp_mask',
 'thermostat_region', 'avg_ke', 'velo', 'acc', 'avgpos', 'oldpos', 'force']

>>> print al[0].params.keys()
['energy', 'time']

The energy is the total potential energy of each configuration, and velo and mass the atomic velocities and masses respectively. We can get a list of all the energies in a couple of different ways:

>>> energies = [ at.energy for at in al ]       # Method 1
>>> print energies
...
>>> print al.energy                             # Method 2
...

Method 1 uses a Python feature called list comprehension to build a list with one entry for each frame in our AtomsList. The second method works in the same way but saves typing by automatically constructing the loop for us, returning a list of energies.

This trick also works for the Atoms property arrays., e.g. al.velo gives us an array with the first axis corresponding to the frame within the file (length 100 here), the second to the spatial dimension (length 3), and the third to the atom number (length 216 here)

>>> print al.velo.shape
(100, 3, 216)
>>> print al.velo[0, :, :] # velocities of first configuration
...

To plot the potential energy as a function of time (assuming matplotlib is available), run the commands:

>>> from pylab import *
>>> plot(al.time, al.energy)
>>> xlabel('Time / fs')
>>> ylabel('Energy / eV')
>>> savefig('potentialenergy.png')

The last command is only necessary if you want to save the plot. See the matplotlib documentation for a lot more information about plotting.

_images/potentialenergy.png

We can calculate the kinetic energy of each frame with another list comprehension:

>>> ke = [.5*sum(at.mass*at.velo.norm2()) for at in al]

Alternatively, we could compute the kinetic energy using array arithmetic:

>>> ke = sum(.5*array(al.mass)*sum(array(al.velo)**2,axis=1), axis=1)

The array() is used here to combine the trajectory frames into a single array, and to convert from Fortran one-based indexing to the standard numpy zero-based indexing expected by matplotlib’s plot() function. If you want to retain the one-based indexing, using farray() instead, which returns a FortranArray.

Note the partial sums: the inner sum is over axis 1 of al.velo which is the spatial dimension, and the outer sum is over axis 1 of .5*array(al.mass)*sum(array(al.velo)**2,axis=1) which is the atom number (here of length 216). Let’s add plots of the kinetic and total energies to our graph:

>>> plot(al.time, ke)                              # Plot kinetic energy
>>> plot(al.time, ke + al.energy)                  # Plot total energy

Here’s the graph we get after adding a legend with.

>>> legend(['Potential energy', 'Kinetic energy', 'Total energy'],loc='center right')
_images/totalenergy.png

We can see that the total energy is fairly well conserved, indicating that the time-step used for this simulation was appropriate. Let’s compare the distribution of atomic speeds in the second half of this simulation (to allow some time for equilibriation) with what we’d expect from the Maxwell-Boltzmann distribution of speeds, given by

\[f(v)\,\mathrm{d}v = 4 \pi \left( \frac{m}{2 \pi k_B T} \right)^{3/2} v^2 \exp \left[ -\frac{mv^2}{2 k_B T} \right] \mathrm{d}v\]

The following function implements this distribution:

def max_bolt(m,T,v):
   "Maxwell-Boltmann distribution of speeds at temperature T for particles of mass m"
   return 4*pi*(m/(2*pi*BOLTZMANN_K*T))**(3.0/2.0)*(v**2)*exp(-m*v**2/(2*BOLTZMANN_K*T))

In this function v can be either a scalar or an array of velocities. The first line of the function after the declaration is a documentation string (or docstring), used for providing online help. According to the principle of equipartition of energy, we expect half of the initial kinetic energy to be converted to potential energy, so that after thermalisation the temperature of our system should be close to 500K. We can plot the expected distribution at a temperature of 500K like this:

>>> vs = linspace(0.0,0.02,100)  # min, max, number of points
>>> plot(vs, max_bolt(at.mass[1], 500.0, vs))

(If you want to clear the previous plot first, use the clf() matplotlib function). To add a histogram of the speeds in the last frame of our AtomsList to the plot, we make use of matplotlibs hist() function:

>>> hist(al[-1].velo.norm(), normed=True, bins=20)

The norm() function returns the norm of each element in a 3 x N matrix as an array of length N. The agreement is not that good as we don’t have enough data. Let’s average over the last 50 snapshots to improve the sampling:

>>> hist(hstack([at.velo.norm() for at in al[-50:]]), normed=True, bins=20, alpha=0.5)

hstack() is a numpy function which stacks a sequence of arrays horizontally to form a composite array. The alpha parameter to hist() causes the bars to be partially transparent. After setting the axis labels and legends, we get the figure shown below.

>>> xlabel(r'Velocity / $\AA$ fs$^{-1}$')
>>> ylabel('Frequency')
>>> legend(['Maxwell-Boltzmann Distribution', 'Simulation Results'])

Note

LaTeX mathematical symbols can be used in matplotlib plots by enclosing them between dollar signs: $...$. Note how we put an r before the string containing the x-axis label. This causes Python to treat it as a raw string, meaning that backslash characters do not have to be escaped.

_images/velocitydistribution.png

Potentials and DynamicalSystems

Now let’s look at how to run simulations entirely using quippy. This will only be practical for small systems which run quickly, but can be very useful to get a feel for what’s going on. For this example we’ll build up the tools necessary to generate the si-1000.xyz file used above. The code for this example is contained in the script file:examples/md.py in the quippy source tree.

The first step is to define the initial atomic configuration. We’ll use a 3 x 3 x 3 supercell of the 8-atom silicon unit cell:

>>> s = supercell(diamond(5.44, 14), 3, 3, 3)

Next we need to create a Potential object to represent the interatomic force model to be used for calculating energies and forces. The Potential object is created by with two string arguments, called args_str and param_str. The first of these identifies the type of potential, and the second is an XML string containing the parameters. QUIP contains a database of parameter files in the share/Parameters. We’ll use a Stillinger-Weber interatomic potential (args_str=”IP SW”) with the original parameters published in Phys. Rev. B 31, 5262 (1984). Python triple-quoted strings can be split over multiple lines, so we can include the XML string easily:

>>> pot = Potential('IP SW', param_str="""<SW_params n_types="1">
 <comment> Stillinger and Weber, Phys. Rev. B  31 p 5262 (1984)</comment>
 <per_type_data type="1" atomic_num="14" />

 <per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
       p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />

 <per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
 lambda="21.0" gamma="1.20" eps="2.1675" />
 </SW_params>
 """)

Alternatively, you can read the XML string from a file:

>>> pot = Potential('IP SW', param_filename='sw.xml')

If no XML parameters are given, the QUIP parameter database (located in the directory /share/Parameters) is searched for matching parameters:

>>> pot = Potential('IP SW')

Our new potential object has several methods, the most important being calc() which is used whenever we want to calculate energies, forces or virial stresses. calc() operates on an Atoms object. We ask for the quantities we are interested in, which can either be returned in arrays or stored within the Atoms object. Let’s try to calculate the energy of our silicon cell:

>>> energy = farray(0.0)   # Create a rank-0 array to store energy
>>> pot.calc(s, energy=energy)  # Do the calculation
>>> print energy
-936.325908705

Note

In earlier versions of quippy, the following error may be raised:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jk2/lib/python2.6/site-packages/quippy/oo_fortran.py", line 438, in <lambda>
    wrapmethod = lambda name: lambda self, *args, **kwargs: self._runroutine(name, *args, **kwargs)
  File "/home/jk2/lib/python2.6/site-packages/quippy/oo_fortran.py", line 316, in _runroutine
   res = fobj(*newargs, **newkwargs)
RuntimeError: Atoms_N_Neighbours: Atoms structure has no connectivity data. Call calc_connect first.

The RuntimeError exception gets raised whenever something goes wrong within a Fortran routine. In this case we need to calculate the atomic connectivity before we can evaluate the energy. To do this we need to choose a cutoff distance, which we can do in sensible fashion by adding a “crust” distance to the cutoff of the potential itself (the crust should be chosen so that no atom moves further than this distance inbetween recalculations of the connectivity).

>>> s.set_cutoff(pot.cutoff() + 2.0)    # crust of 2.0 Angstrom
>>> s.calc_connect()
>>> pot.calc(s, energy=energy)
>>> print energy
-936.325908705

If we now calculate forces for our crystal, we find their all almost zero by symmetry in the unperturbed bulk configuration:

>>> pot.calc(s, force=True)
>>> print s.force.max()
1.74703475792e-14

As an example, let’s make a copy of the Atoms object and randomise the positions a little.

>>> c = s.copy()
>>> c.pos += numpy.random.uniform(-0.1, 0.1, 3*c.n).reshape(3,c.n)
>>> c.calc_connect()
>>> pot.calc(c, force=True)

The extended AtomEye plugin can draw arrows to represent these forces using the draw_arrows() function:

>>> view(c)
>>> draw_arrows('force')
_images/forces.png

To perform a molecular dynamics calculation using our potential, we need to create a DynamicalSystem object from our Atoms object. This stores all the dynamical variables and has methods to advance the simulation using the Velocity-Verlet algorithm. We’ll construct a DynamicalSystem, and randomise the initial velocities to 1000K:

>>> ds = DynamicalSystem(s)
>>> ds.rescale_velo(1000.0)

The molecular dynamics loop would then look something like this:

n_steps          = 100           # Number of simulation steps
dt               = 1.0           # Time-step, in fs
connect_interval = 10            # How frequently to update connectivity
ds.atoms.calc_connect()
for n in range(n_steps):
    ds.advance_verlet1(dt)
    pot.calc(ds.atoms, force=True, energy=True)
    ds.advance_verlet2(dt, ds.atoms.force)
    ds.print_status(epot=ds.atoms.energy)
    if n % connect_interval == 0:
        ds.atoms.calc_connect()

advance_verlet1() advances the velocities by half the time-step \(\delta t\) and the positions by the full time-step. We then calculate at new forces from our potential at \(\mathbf{r}(t+\delta t)\). Finally advance_verlet2() advances the velocities by the second half time-step. print_status() prints a status line with the current time, temperature, total energy, etc.

Note

DynamicalSystem has a run() method which implements an MD loop like the example and returns a list of snapshots:

traj    = ds.run(pot, dt=1.0, n_steps=1000, save_interval=10)
energies = [ at.energy for at in traj ]

The run() method also has an trajectory argument which can be used in conjunction with write_interval to write frames to a file as the simulation runs, for example to write every 100th frame to traj.xyz:

outfile = CInOutput('traj.xyz', OUTPUT)
traj = ds.run(pot, dt=1.0, n_steps=1000, trajectory=outfile, write_interval=100)
outfile.close()

See also the Dynamics wrapper around the DynamicalSystem class, which makes it compatible with the ASE MolecularDynamics interface.

Structural optimisation

If you want to optimise the positions and/or the lattice of an atomic configuration, you’ll need to use either the potential.Potential.minim() method or the potential.Minim helper class.

In this tutorial we’ll consider the Stillinger-Weber potential we defined above. Let’s use this to relax the positions and lattice of a cubic Si cell using conjugate-gradients (cg), to a tolerance of \(|\mathbf{f}|^2 < 10^{-7}\), using a maximum of 100 steps.

>>> pot = Potential('IP SW')
>>> at0 = diamond(5.44, 14)
>>> at0.set_cutoff(pot.cutoff() + 2.0)   #crust of 2.0 Angstrom
>>> at0.calc_connect()
>>> pot.minim(at0, 'cg', 1e-7, 100, do_pos=True, do_lat=True)
Welcome to minim()
space is 33 dimensional
Method: Conjugate Gradients
cg iter = 1 df^2 = .0119620309597580  f = -34.6787373594354591 max(abs(df)) = .0631454167240399
 dcos = .0000000000000000 q = .0000000000000000
resetting eps to .0012500000000000
*** Resetting conjugacy
cg iter = 3 df^2 = .0000000350553860  f = -34.6799999956617953 max(abs(df)) = .0001080977737070
 dcos = -1.0000000000000000 q = 1.0000000000003677
  Linmin: Ea = -34.6799999956617953 a = .0000000000000000
  Linmin done 9 bracket and 5 steps Ex= -34.6799999993383778 x= .2097614050995762
Converged after step 4
Goodbye from minim()

minim() overwrites the original atoms object with the new positions and lattice coordinates, so we can check the lattice properties of the minimised configuration:

>>> a, b, c, alpha, beta, gamma = get_lattice_params(at0.lattice)
>>> print (a, b, c)
(5.4309497787368581, 5.4309497787368581 5.4309497787368581)
>>> print degrees((alpha,beta,gamma))
(90.0, 90.0, 90.0)

i.e. we get a cubic cell with lattice constant of 5.4309498 A, which is correct for this potential.

To use the potential.Minim class instead, which implements the ASE optimizer inferface, do the following:

>>> pot = Potential('IP SW')
>>> at0 = diamond(5.44, 14)
>>> at0.set_calculator(pot)
>>> minimiser = Minim(at0)
>>> minimiser.run(fmax=1e-3)
>>> print np.diag(at0.cell)
[ 5.43094978  5.43094978  5.43094978]

Note the Atoms.set_calculator() call, and that it is no longer necessary to set the Atoms.cutoff, since Minim takes care of this.

Bulk Modulus Calculations

In this example, we will calculate the bulk modulus of silicon using the Stillinger-Weber potential. We’ll do this by fitting the Birch-Murnaghan equation of state, which assumes the bulk modulus is of the form

\[B = B_0 + B'_0 P\]

The energy as a function of volume is

\[E(V) = E_0 + \frac{9V_0B_0}{16} \left\{ \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^3B_0^\prime + \left[\left(\frac{V_0}{V}\right)^\frac{2}{3}-1\right]^2 \left[6-4\left(\frac{V_0}{V}\right)^\frac{2}{3}\right]\right\}.\]

Here is a Python function which implements this relation:

def birch_energy(vo,eo,bo,bop,v):
   t = (vo/v)**.6666666666666667 - 1.0
   e = eo + 1.125*bo*vo*t*t*(1.0+0.5*(bop-4.0)*t)
   return e

We want to calculate the energy for a variety of cell volumes. The routine below takes an initial configuration at0 and first compresses and then expands it in increments of eps. To achieve this we exploit the Atoms.set_lattice() method making sure that scale_positions is true so that all atomic positions are scaled appropriately. Finally, we set the calculator and calculate the energy using Atoms.get_potential_energy() and cell volume using Atoms.get_volume().

def compress_expand(at0, pot, eps=1e-3, relax=False):
    volumes = []
    energies = []
    for i in frange(-5,5):
        at = at0.copy()
        at.set_lattice((1.0 + eps*i) * at.lattice, scale_positions=True)
        at.set_calculator(pot)
        volumes.append(at.get_volume())
        energies.append(at.get_potential_energy())

    return (volumes, energies)

We can use this function to plot an energy-volume curve:

>>> pot = Potential('IP SW')
>>> at0 = diamond(5.43, 14)
>>> volumes, energies = compress_expand(at0, pot)
>>> plot(volumes, energies, 'bo', label='Unrelaxed')
_images/energyvolume1.png

Fitting the equation of state requires the scipy.optimize module, so you need scipy installed to run this example. We need to define an error function, which will make use of the birch_energy() routine defined above.

from scipy.optimize import leastsq

def errfunc(params, volume, energy):
   return birch_energy(*params, v=volume) - energy

The *params entry in the argument list tells Python to use the sequence params for the first four arguments to birch_energy(); we override the volume specifically.

Here’s a function to estimate initial values for the four parameters eo, vo, bo and bop and carry out the least-squares fit:

def birch_fit(energies, volumes):
   energies = farray(energies)
   volumes  = farray(volumes)

   eo = energies.min()
   vo = volumes[energies.argmin()]
   bo = 30.0
   bop = 1.0

   (vo,eo,bo,bop), success = leastsq(errfunc, (vo, eo, bo, bop), args=(volumes, energies))
   print 'Volume vo =', vo, 'A^3'
   print 'Energy eo =', eo, 'eV'
   print 'Bulk modulus bo = ', bo, 'eV/A^3 =', bo*GPA, 'GPa'
   print 'dB/dP (T=0) bop = ', bop

   return (vo, eo, bo, bop)

We can carry out the fit by calling this function:

>>> vo, eo, bo, bop = birch_fit(energies, volumes)
Volume vo = 160.187030023 A^3
Energy eo = -34.6800000521 eV
Bulk modulus bo =  0.632823864562 eV/A^3 = 101.39103958 GPa
dB/dP (T=0) bop =  2.9275321357

We can add a plot of the fitted equations to our graph:

>>> vs = linspace(157, 163, 100)
>>> plot(vs, birch_energy(vo , eo , bo , bop , vs))
_images/energyvolume2.png