Getting started with Atoms

This example will walk you through getting started with using quippy by playing with a small molecule. It complements the (older) introductory tutorial which is longer and goes more in depth, though it is recommended that new users start here first. This example is especially recommended for those who want to get going and working with Atoms as quickly as possible.

First, make sure you have QUIP and quippy properly installed. If you can run the cell below, then congratulations! If not, then see the Installation instructions.

In [1]:
import numpy as np

import quippy
from quippy import Atoms

Remember, this is an interactive tutorial. You are encouraged to run the cells yourself and even try editing the code to see what else you can do!

Now let’s try getting our hands on an Atoms object; this object is required to do almost anything useful in quippy.

In [2]:
struct = quippy.Atoms('methane.xyz')

OK, now what do we do with it? Let’s look at the Atoms source file and see what sort of information is encoded in it:

In [3]:
!cat methane.xyz
5
cutoff=-1.00000000 nneightol=1.20000000 pbc="T T T" Lattice="10.00000000       0.00000000       0.00000000       0.00000000      10.00000000       0.00000000       0.00000000       0.00000000      10.00000000" Properties=species:S:1:pos:R:3:Z:I:1
C               0.00000000      0.00000000      0.00000000       6
H               0.00000000      0.00000000      1.07000000       1
H               0.62414000      0.79255000     -0.35666000       1
H              -0.99844000      0.14425000     -0.35667000       1
H               0.37430000     -0.93680000     -0.35667000       1

We have positions, atomic symbols, atomic numbers, and a unit cell. There’s also some other information in there that we don’t need right now. But let’s see how to access these properties from quippy:

In [4]:
print("Positions:")
print(struct.get_positions())
print("Unit cell:")
print(struct.get_cell())
print("Atomic numbers:")
print(struct.get_atomic_numbers())
Positions:
[[ 0.       0.       0.     ]
 [ 0.       0.       1.07   ]
 [ 0.62414  0.79255 -0.35666]
 [-0.99844  0.14425 -0.35667]
 [ 0.3743  -0.9368  -0.35667]]
Unit cell:
[[ 10.   0.   0.]
 [  0.  10.   0.]
 [  0.   0.  10.]]
Atomic numbers:
[6 1 1 1 1]

(Remember, units in quippy are always Ångstrøms and electronvolts: http://libatoms.github.io/QUIP/units.html)

You might want get an idea of what this structure actually looks like. The native xyz format is supported by many open-source molecular viewers, like VMD and Avogadro. If you’ve gotten the atomeye module to work you can use that as well.

We’re still looking for a simple way to view molecular structures in the notebook environment; once we (the developers) get that working with quippy we’ll update this tutorial.

Let’s try another way to get information about the structure: This is a molecule, so we can check the bonds and angles.

In [5]:
struct.get_distance(0, 1)
Out[5]:
1.0700000000000001

OK, that’s a pretty reasonable C-H bond length in a methane molecule. What about the other ones?

In [6]:
[struct.get_distance(0, i) for i in range(1, 5)]
Out[6]:
[1.0700000000000001,
 1.0699965409757173,
 1.0700018621479124,
 1.0700038406005841]

And how about an H-C-H angle?

In [7]:
struct.get_angle(1, 0, 2)
Out[7]:
109.47090748192832

which is the correct “tetrahedral angle” between the hydrogens in a methane molecule (it’s equal to \(\arccos\left(-\frac{1}{3}\right)\)).

Note that the atom indices in these functions are zero_based, meaning the first atom (here, the carbon) is referred to by the index 0.

ASE and quippy

If you’ve read some of the documentation, you may be aware of the following alternative method for getting the distance between the two atoms:

In [8]:
struct.distance_min_image(0, 1)
Out[8]:
0.0

DANGER WILL ROBINSON! As you may have noticed, this isn’t the correct bond distance. In some circumstances the code may even just crash. This is because the above function is derived from the underlying Fortran code, QUIP, rather than the get_distance function from before, which came from ASE. The QUIP-derived functions can be very useful, but for now it’ll just be too confusing to use them - let’s stick with ASE.

So how do we tell whether a function is derived from QUIP or from ASE? Well, let’s take a look at this function’s help message:

In [9]:
help(struct.distance_min_image)
Help on method <lambda> in module quippy._atoms:

<lambda> lambda self, *args, **kwargs method of quippy.atoms.Atoms instance
    This interface calculates the distance between the nearest periodic images of two points (or atoms).
     Return minimum image distance between two atoms or positions.
     End points can be specified by any combination of atoms indices
     'i' and 'j' and absolute coordinates 'u' and 'w'. If 'shift' is
     present the periodic shift between the two atoms or points will
     be returned in it.

    Wrapper around Fortran interface ``distance_min_image`` containing multiple routines:

       .. function :: distance_min_image(v,w,[shift])

          :param v:
          :type v:  input rank-1 array('d') with bounds (3)
          :param w:
          :type w:  input rank-1 array('d') with bounds (3)
          :param shift:
          :type shift:  in/output rank-1 array('i') with bounds (3), optional

          :returns: **ret_distance8_vec_vec** --  float


          Routine is wrapper around Fortran routine ``distance8_vec_vec`` defined in file :git:`src/libAtoms/Atoms_types.f95`.

       .. function :: distance_min_image(i,j,[shift])

          :param i:
          :type i:  input int
          :param j:
          :type j:  input int
          :param shift:
          :type shift:  in/output rank-1 array('i') with bounds (3), optional

          :returns: **ret_distance8_atom_atom** --  float


          Routine is wrapper around Fortran routine ``distance8_atom_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.

       .. function :: distance_min_image(i,v,[shift])

          :param i:
          :type i:  input int
          :param v:
          :type v:  input rank-1 array('d') with bounds (3)
          :param shift:
          :type shift:  in/output rank-1 array('i') with bounds (3), optional

          :returns: **ret_distance8_atom_vec** --  float


          Routine is wrapper around Fortran routine ``distance8_atom_vec`` defined in file :git:`src/libAtoms/Atoms_types.f95`.

       .. function :: distance_min_image(v,j,[shift])

          :param v:
          :type v:  input rank-1 array('d') with bounds (3)
          :param j:
          :type j:  input int
          :param shift:
          :type shift:  in/output rank-1 array('i') with bounds (3), optional

          :returns: **ret_distance8_vec_atom** --  float


          Routine is wrapper around Fortran routine ``distance8_vec_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.

A bit confusing, yes, but take a look at the final line:

Routine is wrapper around Fortran routine ``distance8_vec_atom`` defined in file :git:`src/libAtoms/Atoms_types.f95`.

It says this function is a wrapper around something in the Fortran code, which means it’s from QUIP and only works with one-based indices (i.e. to this function, the first atom has the index 1).

In [10]:
struct.distance_min_image(1, 2)
Out[10]:
1.07

Whereas the ASE function has a different help string:

In [11]:
help(struct.get_distance)
Help on method get_distance in module ase.atoms:

get_distance(self, a0, a1, mic=False, vector=False) method of quippy.atoms.Atoms instance
    Return distance between two atoms.

    Use mic=True to use the Minimum Image Convention.
    vector=True gives the distance vector (from a0 to a1).

See that part, ‘in module ase.atoms’? That means this function comes from ASE and works with zero-based indices.

You should always do this check before using any Atoms function you’re not yet familiar with.

Useful QUIP equivalents

Just in case you ever want to work with the Fortran QUIP routines, here are some equivalents of the queries we’ve seen above:

In [12]:
[struct.distance_min_image(1, i) for i in range(2, 6)]
Out[12]:
[1.07, 1.0699965409757173, 1.0700018621479124, 1.070003840600584]
In [13]:
print(struct.cosine(1, 2, 3))
print(np.arccos(struct.cosine(1, 2, 3)) * 180 / np.pi)
-0.333328180365
109.470907482

Note that the order of the arguments in Atoms.cosine is different as well: The central atom index is now the first argument.

The properties can also be accessed as one-based arrays (called FortranArrays). Note that these are also the transpose of the ASE arrays, so the row and column indices are swapped.

In [14]:
struct.pos
Out[14]:
FortranArray([[ 0.     ,  0.     ,  0.62414, -0.99844,  0.3743 ],
              [ 0.     ,  0.     ,  0.79255,  0.14425, -0.9368 ],
              [ 0.     ,  1.07   , -0.35666, -0.35667, -0.35667]])
In [15]:
struct.lattice
Out[15]:
FortranArray([[ 10.,   0.,   0.],
              [  0.,  10.,   0.],
              [  0.,   0.,  10.]])
In [16]:
struct.Z
Out[16]:
FortranArray([6, 1, 1, 1, 1], dtype=int32)

Changing Atoms

There’s a more complete tutorial on manipulating Atoms objects in ASE at the ASE website; this is a shorter example just to get you started.

First, let’s try reorienting our methane molecule in space. The first hydrogen atom is already oriented along the Z axis; let’s try to rotate the molecule so that the second one points along the X axis.

In [17]:
r_h2 = struct.get_distance(0, 2, vector=True)
r_h2[2] = 0.0
np.arccos(r_h2[0] / np.linalg.norm(r_h2))
Out[17]:
0.90371859179284519

So we need to rotate it by about -0.9 radians about the Z axis.

But first, we want to keep the old methane molecule for comparison. Let’s copy it - note that this must be done explicitly; the Python assignment operator (=) only assigns a new reference to the underlying object! This can be one of the hardest things to understand for programmers coming from a different language, so make sure to read up on Python references and do some of their tutorials if you find this distinction confusing.

In [18]:
struct_old = struct.copy()
In [19]:
struct.rotate([0.0, 0.0, 1.0], -1.0 * np.arccos(r_h2[0] / np.linalg.norm(r_h2)))
In [20]:
struct.get_positions()
Out[20]:
array([[  0.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   1.07000000e+00],
       [  1.00880436e+00,  -1.11022302e-16,  -3.56660000e-01],
       [ -5.04400083e-01,   8.73653852e-01,  -3.56670000e-01],
       [ -5.04404280e-01,  -8.73653852e-01,  -3.56670000e-01]])

Exercise: Now compare the positions of the rotated structure with those of the new structure; you may want to compute the RMS difference.

Now change the cell and make a supercell

Generating structures from scratch

We could have also just used ASE functionality to pull a methane structure from its database:

In [21]:
import ase
from ase.build import molecule
In [22]:
me = molecule('CH4')
In [23]:
me.get_positions()
Out[23]:
array([[ 0.      ,  0.      ,  0.      ],
       [ 0.629118,  0.629118,  0.629118],
       [-0.629118, -0.629118,  0.629118],
       [ 0.629118, -0.629118, -0.629118],
       [-0.629118,  0.629118, -0.629118]])

Like the introductory tutorial at http://libatoms.github.io/QUIP/tutorial.html - use the diamond generator, view a supercell, do fun things with the structure…

Computing the energy

Introduce Potential. Use something like IP SW to get the energy of the silicon supercell. Direct user to the list of potentials and the next tutorial.