Problems with popular crystal structure software

In order to successfully apply molecular dynamics simulations to real world problems, the systems under investigation must be represented by atomistic (or sometimes coarse-grained) models. Creating these models requires the knowledge of the initial positions of the atoms (or coarse-grained particles).

Many nanofluidic systems often require the modelling of a substrate in contact with a fluid. To model fluids, the initial positions of the atoms/particles can be randomly selected. In many cases, the substrate is usually some type of crystal structure. For example, I have been using a crystal substrate of alumina (aluminium oxide) in my systems.

Most of the time, one can use crystal structure software in order to quickly create these types of substrates. Most software usually uses .cif (Crystallographic Information File) files as input (the American Mineralogist Crystal Structure Database is a good source). The .cif files contain the structural information needed to create the unit cell of the crystal, such as the initial coordinates for the assymetric unit cell, the unit cell dimensions (lengths and angles) and the symmetry group operations for the particular type of crystal. Once the unit cell has been created, a substrate can be quickly built by replicating the unit cell any number of times in any dimension. In the past, I have used a piece of open source software called Avagadro to check the crystal structure of substrates I have built with my own codes, but recently, I found another piece of software called Mercury (free license version).

So, it was to my surprise that after using the same CIF input to both codes, two very different unit cells were calculated. To highlight the major differences between the output from both codes, the CIF for iron oxide has been used as input:

 

The image on the left is from Avagadro and the image on the right is from Mercury. This is worrying, as both should, in theory, calculate the same unit cell. From a simple visual inspection of both structures, you can tell they are very different…so the question is, which software has produced the correct unit cell?

The answer, is neither! While Avagadro preserves the stoic ratio of the crystal (this is very important, as it maintains charge neutrality), it calculates many overlapping atoms (mainly oxygen atoms), which is clearly incorrect. While the structure calculated by Mercury appears to be visually okay, i.e. no overlapping atoms, it seems to calculate a somewhat random number of atoms in the unit cell, not preserving the stoic ratio… I still have not been able to understand how or why Mercury calculates either the total number of atoms or the number of each atomic type.

So the take home message from this post should probably be not to simply trust the output of software. I normally write my own codes to calculate unit cells and then replicate to create substrates, so luckily, these "bugs" have not troubled me. Also, it should be noted, that not every type of crystal has been tested as input to both codes, so both may well produce indentical outputs for other crystals. It was just that I happened to stumble upon the differences for iron oxide!