User guide: crystals and lattices

Handling crystal models is the main feature of the crystals library. This is done through the Crystal class, a representation of a crystal including unit cell atoms, lattice parameters, and other goodies.

Since working with Crystal instances is so important, there are many ways to construct them.

Constructing a Crystal object

Creating a Crystal object can be done most easily from a Crystal Information File (CIF, .cif):

>>> from crystals import Crystal
>>> diamond = Crystal.from_cif('diamond.cif') 

crystals also has an internal database of CIF files. Valid names are stored in Crystal.builtins and can be constructed like so:

>>> Crystal.builtins 
frozenset({'Ac',
           'Ag',
           'Al',
           ...
           'alpha-Mn',
           'b-Pu',
           'diamond',
           'gamma-Pu'
           })
>>> 'Au' in Crystal.builtins
True
>>> Crystal.from_database('Au') 
< Crystal object with following unit cell:
    Atom Au @ (0.00, 0.00, 0.00)
    Atom Au @ (0.00, 0.50, 0.50)
    Atom Au @ (0.50, 0.00, 0.50)
    Atom Au @ (0.50, 0.50, 0.00)
Lattice parameters:
    a=4.078Å, b=4.078Å, c=4.078Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
    Au: 100.000% >

RCSB Protein DataBank files are even easier to handle; simply provide the 4-letter identification code and the structure file will be taken care of by crystals:

>>> hemoglobin = Crystal.from_pdb('1gzx')
>>> print(hemoglobin)
< Crystal object with following unit cell:
    Atom C  @ (-0.50, 0.17, -0.05)
    Atom C  @ (-0.49, 0.19, -0.16)
    Atom C  @ (-0.49, 0.13, -0.15)
    Atom C  @ (-0.48, 0.30, -0.09)
    Atom C  @ (-0.48, 0.17, -0.06)
    Atom C  @ (-0.48, 0.17, -0.15)
    Atom C  @ (-0.48, 0.17, -0.04)
    Atom C  @ (-0.48, 0.23, -0.15)
    Atom C  @ (-0.47, 0.13, -0.16)
    Atom C  @ (-0.47, 0.23, -0.12)
      ... omitting 4554 atoms ...
      ... use repr() to show the full cell ...
Lattice parameters:
    a=97.050Å, b=99.500Å, c=66.110Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
     C: 64.724%
    Fe: 0.088%
     N: 17.090%
     O: 17.835%
     S: 0.263% >

Another convenient way to construct a Crystal is through the Crystallography Open Database:

>>> # Default is the latest revision
>>> vo2 = Crystal.from_cod(1521124)
>>> # Revisions are accessible as well
>>> old_vo2 = Crystal.from_cod(1521124, revision = 140771)

the Materials Project provides another avenue where to get crystal structures. You will need an API key from your account dashboard:

>>> fe2o3 = Crystal.from_mp(api_key="xxxxxxxxxxxxxxxx", query = "Fe2O3") 
>>> print(fe2o3) 
< Crystal object with following unit cell:
    Atom Fe @ (0.89, 0.37, 0.67)
    Atom Fe @ (0.99, 0.75, 0.96)
    Atom Fe @ (0.49, 0.75, 0.79)
    Atom Fe @ (0.25, 0.51, 0.71)
    Atom Fe @ (0.13, 0.39, 0.92)
    Atom Fe @ (0.88, 0.12, 0.59)
    Atom Fe @ (0.38, 0.38, 0.84)
    Atom Fe @ (0.62, 0.62, 0.66)
    Atom Fe @ (0.61, 0.87, 0.58)
    Atom Fe @ (0.38, 0.38, 0.16)
    ... omitting 150 atoms ...
    ... use repr() to show the full cell ...
Lattice parameters:
    a=8.525Å, b=8.525Å, c=25.593Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
     O: 60.000%
    Fe: 40.000% >

Other constructors are supported. See the reference for the Crystal class for more details.

Constructing a Crystal object by hand

If you don’t have a file on hand, or want to create an idealized crystal, consider building a Crystal object by hand.

To do this, you need:

1. iterable of Atom objects, with coordinates. These atoms must be the full unit cell;

2. three lattice vectors;

As an example, let’s create the simplest crystal structure known: alpha-Polonium (simple cubic):

>>> from crystals import Crystal, Atom
>>> import numpy as np
>>>
>>> lattice_vectors = 3.35 * np.eye(3)
>>> unitcell = [Atom('Po', coords = [0,0,0])]
>>>
>>> polonium = Crystal(unitcell, lattice_vectors)
>>> polonium
< Crystal object with following unit cell:
    Atom Po @ (0.00, 0.00, 0.00)
Lattice parameters:
    a=3.350Å, b=3.350Å, c=3.350Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
    Po: 100.000% >

In the case where atoms are given as an asymmetric unit cell and a set of symmetry operators, you can use the symmetry_expansion() function to generate a set of unique atoms (even if some symmetry operators might be redundant). The generated set of atoms can be passed to the constructor of Crystal.

Converting a Crystal to other formats

You can use the crystals package to convert crystal structures from one format to another. Currently, you can write a structure either to an .xyz (Crystal.to_xyz()) file, a Crystallography Information Framework .cif (Crystal.to_cif()), or to an ase.Atoms structure (Crystal.to_ase()). Here is an example:

>>> from crystals import Crystal
>>> import numpy as np
>>>
>>> # Create a crystal structure by hand
>>> lattice_vectors = 3.35 * np.eye(3)
>>> unitcell = [Atom('Po', coords = [0,0,0])]
>>> polonium = Crystal(unitcell, lattice_vectors)
>>>
>>> # Convert to CIF
>>> polonium.to_cif('polonium.cif') 

Would you like to convert to another format that is not supported yet? Please raise an issue!

Crystal attributes

The Crystal object provides some interfaces for easy structure manipulation. First, a Crystal is an iterable:

>>> from crystals import Crystal
>>> graphite = Crystal.from_database('C')
>>>
>>> for atm in graphite: 
...     print(repr(atm)) 
...
< Atom C  @ (0.33, 0.67, 0.25) >
< Atom C  @ (0.00, 0.00, 0.75) >
< Atom C  @ (0.00, 0.00, 0.25) >
< Atom C  @ (0.67, 0.33, 0.75) >

The len() of a Crystal is the unit cell size (in number of atoms):

>>> hemoglobin = Crystal.from_pdb('1gzx')
>>> len(hemoglobin)
4564

The Crystal class is a set-like container; checking containership (with the builtin in statement) is very fast:

>>> graphite = Crystal.from_database('C')
>>> carbon = next(iter(graphite))                   # Pick any atom from the crystal
>>>
>>> carbon in graphite 
True

Crystal instances can be equated to each other:

>>> gold = Crystal.from_database('Au')
>>> silver = Crystal.from_database('Ag')
>>>
>>> gold == silver
False

Just like lists and other container types, a Crystal instance is False if empty, and True otherwise:

>>> mycrystal = Crystal.from_database('Cu')
>>> if mycrystal: # equivalent to: if len(mycrystal) > 0:
...     pass

Structures can be extracted from a Crystal instance by making use of its superclass, AtomicStructure. For example, all atoms satisfying a certain condition can be found using Crystal.satisfying:

>>> vo2 = Crystal.from_database('vo2-m1')
>>>
>>> vo2.satisfying( lambda atom: atom.element == 'V' )
< AtomicStructure object with following orphan atoms:
    Atom V  @ (0.24, 0.53, 0.53)
    Atom V  @ (0.24, 0.97, 0.03)
    Atom V  @ (0.76, 0.03, 0.97)
    Atom V  @ (0.76, 0.48, 0.47) >

To make it easier, take a look at the is_element() function:

>>> from crystals import is_element
>>>
>>> vo2 = Crystal.from_database('vo2-m1')
>>> vo2.satisfying( is_element('O') )
< AtomicStructure object with following orphan atoms:
    Atom O  @ (0.10, 0.21, 0.20)
    Atom O  @ (0.10, 0.29, 0.70)
    Atom O  @ (0.39, 0.69, 0.29)
    Atom O  @ (0.39, 0.81, 0.79)
    Atom O  @ (0.61, 0.19, 0.21)
    Atom O  @ (0.61, 0.31, 0.71)
    Atom O  @ (0.90, 0.71, 0.30)
    Atom O  @ (0.90, 0.79, 0.80) >

If a Crystal was generated from a file, the path to its file can be retrieved from the source attribute:

>>> hemoglobin = Crystal.from_pdb('1gzx')
>>> print(hemoglobin.source) 
'(...omitted...)\pdb1gzx.ent'

Crystal instances have a nice string representation (with str()), and a complete string representation (with repr()):

>>> vo2 = Crystal.from_database('vo2-m1')
>>> print(vo2)         # Short string representation
< Crystal object with following unit cell:
    Atom O  @ (0.10, 0.21, 0.20)
    Atom O  @ (0.10, 0.29, 0.70)
    Atom O  @ (0.39, 0.69, 0.29)
    Atom O  @ (0.39, 0.81, 0.79)
    Atom O  @ (0.61, 0.19, 0.21)
    Atom O  @ (0.61, 0.31, 0.71)
    Atom O  @ (0.90, 0.71, 0.30)
    Atom O  @ (0.90, 0.79, 0.80)
    Atom V  @ (0.24, 0.53, 0.53)
    Atom V  @ (0.24, 0.97, 0.03)
      ... omitting 2 atoms ...
      ... use repr() to show the full cell ...
Lattice parameters:
    a=5.743Å, b=4.517Å, c=5.375Å
    α=90.000°, β=122.600°, γ=90.000°
Chemical composition:
     O: 66.667%
     V: 33.333% >

Crystal instances can be converted to NumPy arrays as well:

>>> import numpy
>>> numpy.array(Crystal.from_database('Si'))
array([[14.  ,  0.  ,  0.  ,  0.  ],
       [14.  ,  0.  ,  0.5 ,  0.5 ],
       [14.  ,  0.25,  0.25,  0.25],
       [14.  ,  0.25,  0.75,  0.75],
       [14.  ,  0.5 ,  0.  ,  0.5 ],
       [14.  ,  0.5 ,  0.5 ,  0.  ],
       [14.  ,  0.75,  0.25,  0.75],
       [14.  ,  0.75,  0.75,  0.25]])

arr will contain one row per unit cell atom:

You can calculate what the asymmetric cell of a Class is with the Crystal.asymmetric_cell() method:

>>> Crystal.from_database('C').asymmetric_cell() 
{< Atom C  @ (0.00, 0.00, 0.25) >,
 < Atom C  @ (0.67, 0.33, 0.75) >}

Atomic Number	Fractional coordinates
Z1	x1	y1	z1
Z2	x2	y2	z2
Z3	x3	y3	z3
…

Lattice vectors and reciprocal space

Once a Crystal object is ready, you can manipulate the lattice parameters via the Lattice super-class. Let’s use the built-in example of graphite:

>>> from crystals import Crystal
>>> import numpy as np
>>>
>>> graphite = Crystal.from_database('C')
>>>
>>> a1, a2, a3 = graphite.lattice_vectors
>>> a1
array([ 2.13388659e+00, -1.23200000e+00,  1.50876486e-16])
>>> b1, b2, b3 = graphite.reciprocal_vectors
>>> b1
array([2.94447949, 0.        , 0.        ])
>>>
>>> np.dot(a1, b1)  # 2 pi
6.283185307179586

The standard three lengths and angles description of a lattice is also accessible:

>>> graphite = Crystal.from_database('C')
>>> a, b, c, alpha, beta, gamma = graphite.lattice_parameters

The unit cell volume (and by extensions, density) is also accessible:

>>> graphite = Crystal.from_database('C')
>>> vol = graphite.volume   # Angstroms cubed
>>> density = vol/len(graphite)

As a lattice can be fully described by its basis vectors, a NumPy array can be created from a Lattice instance.

Space-group Information

The lattice system of a Lattice or Crystal instance is also available via the lattice_system attribute:

>>> vo2 = Crystal.from_database('vo2-m1') # Monoclinic M1 VO2
>>> vo2.lattice_system
<LatticeSystem.monoclinic: 2>

Better control on length tolerances is available via the lattice_system() function.

Thanks to spglib, we can get space-group information from a Crystal instance:

>>> from pprint import pprint # pretty printing
>>>
>>> gold = Crystal.from_database('Au')
>>> pprint(gold.symmetry())
{'centering': <CenteringType.face_centered: 'F'>,
 'hall_number': 523,
 'hall_symbol': '-F 4 2 3',
 'hm_symbol': 'Fm-3m',
 'international_full': 'F 4/m -3 2/m',
 'international_number': 225,
 'international_symbol': 'Fm-3m',
 'pointgroup': 'm-3m'}

In the above example, spg_info is a dictionary with the following keys:

1. 'international_symbol': International Tables of Crystallography space-group symbol (short);

2. 'international_full': International Tables of Crystallography space-group full symbol;

3. 'hall_symbol' : Hall symbol;

4. 'hm_symbol' : Hermann-Mauguin symbol;

5. 'pointgroup' : International Tables of Crystallography point-group;

6. 'international_number' : International Tables of Crystallography space-group number (between 1 and 230);

7. 'hall_number' : Hall number (between 1 and 531).

Each of those items are also available directly from the Crystal instance. The Hall number of a crystal structure is located in the Crystal.hall_number attribute, the short international symbol is located in the Crystal.international_symbol. attribute, and so on.

Symmetry operations

You can get the matrix symmetry operations directly from the Crystal class:

>>> cryst = Crystal.from_database('C')
>>> first_symop = cryst.symmetry_operations()[0]
>>>
>>> print(first_symop)
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

Symmetry operations are described using 4x4 affine matrices, where the rotation is the top 3x3 block, and the translation is the right-most column. Example with iteration:

>>> for m in cryst.symmetry_operations():
...     rotation = m[:3, :3]
...     translation = m[:3, -1]
...
>>>

Symmetry operations in reciprocal space are also made available via Crystal.reciprocal_symmetry_operations().

Grouping atoms by site-symmetry

For various reasons, it might be useful to know which of the atoms in a Crystal are related by some site symmetry, for example Wyckoff letters or crystallographic orbits. This can be done with the Crystal.groupby() method, which groups atoms according to their site symmetry. For example:

>>> graphite = Crystal.from_database('C')
>>> groups = graphite.groupby(by="crystallographic_orbits")
>>> for orbit, atoms in groups.items():                                         
...     print(f"Orbit {orbit}: {atoms}")                                        
...                                                                             
Orbit 0: [< Atom C  @ (0.00, 0.00, 0.25) >, < Atom C  @ (0.00, 0.00, 0.75) >]   
Orbit 2: [< Atom C  @ (0.33, 0.67, 0.25) >, < Atom C  @ (0.67, 0.33, 0.75) >]   

Supported site-symmetry measures are currently “crystallographic_orbits”, “wyckoffs”, and “equivalent_atoms”. See the spglib documentation for a description of these measures of symmetry.

Cell refinements and manipulations

Again, through spglib, we can create different versions of Crystal instances. For example, a primitive Crystal can be created using the Crystal.primitive() method:

>>> gold = Crystal.from_database('Au')
>>> print(gold)
< Crystal object with following unit cell:
    Atom Au @ (0.00, 0.00, 0.00)
    Atom Au @ (0.00, 0.50, 0.50)
    Atom Au @ (0.50, 0.00, 0.50)
    Atom Au @ (0.50, 0.50, 0.00)
Lattice parameters:
    a=4.078Å, b=4.078Å, c=4.078Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
    Au: 100.000% >
>>>
>>> primitive_gold = gold.primitive() # this is a whole new Crystal instance
>>> print(primitive_gold)
< Crystal object with following unit cell:
    Atom Au @ (0.00, 0.00, -0.00)
Lattice parameters:
    a=2.884Å, b=2.884Å, c=2.884Å
    α=60.000°, β=60.000°, γ=60.000°
Chemical composition:
    Au: 100.000% >

Notice how the primitive structure is much simpler.

Idealized versions of Crystal objects are also made available. Let’s take the example of iron arsenide:

>>> iron_arsenide = Crystal.from_database('FeAs')
>>> print(iron_arsenide)
< Crystal object with following unit cell:
    Atom As @ (0.20, 0.58, 0.25)
    Atom As @ (0.30, 0.08, 0.75)
    Atom As @ (0.70, 0.92, 0.25)
    Atom As @ (0.80, 0.42, 0.75)
    Atom Fe @ (0.00, 0.20, 0.25)
    Atom Fe @ (0.50, 0.70, 0.75)
    Atom Fe @ (0.50, 0.30, 0.25)
    Atom Fe @ (1.00, 0.80, 0.75)
Lattice parameters:
    a=5.440Å, b=6.026Å, c=3.371Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
    As: 50.000%
    Fe: 50.000% >
>>> idealized = iron_arsenide.ideal() # this is a whole new Crystal instance
>>> print(idealized)
< Crystal object with following unit cell:
    Atom As @ (0.20, 0.25, 0.42)
    Atom As @ (0.30, 0.75, 0.92)
    Atom As @ (0.70, 0.25, 0.08)
    Atom As @ (0.80, 0.75, 0.58)
    Atom Fe @ (0.00, 0.25, 0.80)
    Atom Fe @ (0.50, 0.75, 0.30)
    Atom Fe @ (0.50, 0.25, 0.70)
    Atom Fe @ (1.00, 0.75, 0.20)
Lattice parameters:
    a=5.440Å, b=3.371Å, c=6.026Å
    α=90.000°, β=90.000°, γ=90.000°
Chemical composition:
    As: 50.000%
    Fe: 50.000% >

The atomic coordinates have been swapped to follow conventions of the appropriate space-group.

Scattering utilities

Lattice objects have a few methods that make life easier when dealing with scattering data and modeling.

The conversion between Miller indices and scattering vectors is available:

>>> from crystals import Crystal
>>> graphite = Crystal.from_database('C')
>>>
>>> # Behavior inherited from Lattice superclass
>>> e1 = (1, 0, 0)
>>> G = graphite.scattering_vector(e1)
>>> graphite.miller_indices(G)
array([1., 0., 0.])

Since version 1.3, the calculation of Miller indices or scattering vectors is much faster for tables of indices, where every row corresponds to a reflection/scattering vector:

>>> from crystals import Crystal
>>> import numpy as np
>>>
>>> graphite = Crystal.from_database('C')
>>> indices = np.array([
...    [1,0,0],
...    [1,1,0],
...    [1,2,1],
...    [2,0,0],
... ])
>>> graphite.scattering_vector(indices)
array([[2.94447949, 0.        , 0.        ],
       [4.41671923, 2.54999404, 0.        ],
       [5.88895897, 5.09998807, 0.93625172],
       [5.88895897, 0.        , 0.        ]])

Supercells

For various reasons, creating a supercell from a crystal might be desirable. The process is very easy. Let’s imagine we want to create a 2x2x2 supercell of graphite:

>>> from crystals import Crystal
>>> graphite = Crystal.from_database('C')
>>>
>>> for atm in sorted(graphite): 
...     print(repr(atm))         
...                              
< Atom C  @ (0.00, 0.00, 0.25) >
< Atom C  @ (0.00, 0.00, 0.75) >
< Atom C  @ (0.33, 0.67, 0.25) >
< Atom C  @ (0.67, 0.33, 0.75) >
>>>
>>> for atm in sorted(graphite.supercell(2,2,2)): 
...    print(repr(atm))                           
...                                               
< Atom C  @ (0.00, 0.00, 0.25) >
< Atom C  @ (0.00, 0.00, 0.75) >
< Atom C  @ (0.00, 0.00, 1.25) >
< Atom C  @ (0.00, 0.00, 1.75) >
< Atom C  @ (0.00, 1.00, 0.25) >
< Atom C  @ (0.00, 1.00, 0.75) >
< Atom C  @ (0.00, 1.00, 1.25) >
< Atom C  @ (0.00, 1.00, 1.75) >
< Atom C  @ (0.33, 0.67, 0.25) >
< Atom C  @ (0.33, 0.67, 1.25) >
< Atom C  @ (0.33, 1.67, 0.25) >
< Atom C  @ (0.33, 1.67, 1.25) >
< Atom C  @ (0.67, 0.33, 0.75) >
< Atom C  @ (0.67, 0.33, 1.75) >
< Atom C  @ (0.67, 1.33, 0.75) >
< Atom C  @ (0.67, 1.33, 1.75) >
< Atom C  @ (1.00, 0.00, 0.25) >
< Atom C  @ (1.00, 0.00, 0.75) >
< Atom C  @ (1.00, 0.00, 1.25) >
< Atom C  @ (1.00, 0.00, 1.75) >
< Atom C  @ (1.00, 1.00, 0.25) >
< Atom C  @ (1.00, 1.00, 0.75) >
< Atom C  @ (1.00, 1.00, 1.25) >
< Atom C  @ (1.00, 1.00, 1.75) >
< Atom C  @ (1.33, 0.67, 0.25) >
< Atom C  @ (1.33, 0.67, 1.25) >
< Atom C  @ (1.33, 1.67, 0.25) >
< Atom C  @ (1.33, 1.67, 1.25) >
< Atom C  @ (1.67, 0.33, 0.75) >
< Atom C  @ (1.67, 0.33, 1.75) >
< Atom C  @ (1.67, 1.33, 0.75) >
< Atom C  @ (1.67, 1.33, 1.75) >

Supercells are generated by copying unit cells along crystallographic axes. For example, a 2x3x5 supercell will be extended by 2x along the a1 lattice vector, extended by 3x along the a2 lattice vector, and extended by 5x along the a3 lattice vector.

Supercell objects can be used anywhere a Crystal object can be used. With the addition of all methods and attributes from the Crystal class, a Supercell object also has two extra attributes:

1. Supercell.dimensions: Integer supercell dimensions. For example, for a 3x2x1 supercell, the dimensions are (3,2,1).

2. Supercell.scaled_lattice_vectors: Lattice vectors scaled by the Supercell.dimensions.

Compatibility with ASE

The Atomic Simulation Environment is a powerful tool. You can harness its power and convert between ase.Atoms and crystals.Crystal at will.

To create an ase.Atoms object from a Crystal, use the Crystal.ase_atoms() method:

>>> from ase.calculators.abinit import Abinit 
>>> from crystals import Crystal
>>>
>>> gold = Crystal.from_database('Au')
>>> ase_gold = gold.to_ase(calculator = Abinit(...)) 

All keywords of the ase.Atoms constructor are supported. To get back to a Crystal instance:

>>> gold2 = Crystal.from_ase(ase_gold)