# 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, mic=True)
```

```
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, mic=True) for i in range(1, 5)]
```

```
Out[6]:
```

```
[1.0700000000000001,
1.0699965409757173,
1.0700018621479124,
1.0700038406005841]
```

*(Oh, and in case you’re wondering about the ``mic=True`` part, that
just means to use the minimum image convention. It’s only really
relevant when we use periodic systems, but it’s included here for
completeness.)*

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. Fortran, unlike Python and C, starts counting from one
instead of zero (cue the holy wars), so we just gave QUIP an atom index
that doens’t even exist – and that can have unpredictable results. Now,
the QUIP-derived functions are often useful, but for now it’ll just be
too confusing to keep this indexing distinction in our heads - let’s
stick with ASE until we have good reason to do otherwise.

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
`FortranArray`

s). 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.

```
In [ ]:
```

```
```