This section looks at using Python to automatically process .magres output files, specifically looking at the example of processing disordered systems.

Full documentation for the Python module is available here.

For a alternative introduction see http://tfgg.me/magres-python/


To set up the python environment you will need to unload the intel compilers and load python2.7

module unload intel-compilers
module load module load python/2.7__gcc-4.8

Note that you will have to re-load the intel-compiler module before re-running CASTEP module load intel-compilers.

You can either use git to clone the latest version of the repository or download the package as a .zip. If you're using cygwin, you probably want to do the latter unless you've installed git.

git clone git://github.com/tfgg/magres-format.git
cd magres-format


curl -L "https://github.com/tfgg/magres-format/archive/master.zip" > master.zip
unzip master.zip
cd magres-format-master

Now run setup.py to install the module and associated scripts. If you want to install it locally (e.g. if you're on a cluster and donít have sudo), run

python setup.py install --user

and it should now be available.

If you want to install it system wide, run

sudo python setup.py install

Python should already be installed on your system. This Python module also requires the numpy module.


The module should now be installed. In the install folder, cd into the "samples" folder. There should be a number of .magres files. Start an interactive Python prompt by running the "python" command with no arguments.

Load the "ethanol-nmr.magres" samples file into a MagresAtoms structure with the following code.

>>> from magres.atoms import MagresAtoms
>>> atoms = MagresAtoms.load_magres('ethanol-nmr.magres')

Now try printing out all the atoms and their isotropic magnetic shieldings (note the 'tab' indentation in the second line)

>>> for atom in atoms:
...   print atom, atom.ms.iso

Try printing out all the proton spans and skews

>>> for atom in atoms.species('H'):
...   print atom, atom.ms.span, atom.ms.skew

Here we've used the species selector to just loop over the hydrogen atoms in the system.

Try loading .magres files from some of the calculations that you've performed in the previous sections. These are some of the available properties on atoms if the appropriate calculation has been performed:

  • atom
    • atom.position - Cartesian position of the atom in angstrom
    • atom.species - Elemental symbol of the atom
    • atom.label - Arbitrary label given to the atom in the calculation
    • atom.index - Index of the atom within the label
    • atom.isotope - The atom's isotope. You can change this by assigning a new isotope, e.g. atom.isotope = 63
    • atom.gamma - The gyromagnetic ratio for this particular isotope
    • atom.Q - The quadrupole moment for this particular isotope
  • atom.efg
    • atom.efg.V - The full electric field gradient tensor
    • atom.efg.Cq - The Cq of the V tensor
    • atom.efg.eta - The asymmetry
  • atom.ms
    • atom.ms.sigma - The full magnetic shielding tensor
    • atom.ms.iso - The isotropic magnetic shielding
    • atom.ms.aniso - The anisotropic magnetic shielding
    • atom.ms.eta - The shielding asymmetry
  • atom.isc
    • atom1.isc[atom2].K - The reduced coupling tensor
    • atom1.isc[atom2].J - The full coupling tensor
    • atom1.isc[atom2].K_iso - The reduced isotropic coupling.

And so on. You can read about all the available properties and their definitions in the documentation.

Example 5 - Gehlenite electric field gradients (EFGs)

Gehlenite is a disordered aluminosilicate consisting of two types of aluminium site and one type of silicon site. We want to look at aluminium EFGs to quantify the amount of disorder in the material, in particular the degree of violation of the "Loewenstein" rule (avoidance of Al-O-Al linkages) and its structural effects. This section is based on analysis performed in Elucidation of the Al/Si Ordering in Gehlenite Ca2Al2SiO7 by Combined 29Si and 27Al NMR Spectroscopy/Quantum Chemical Calculations.

Two sets of EFG calculations have been performed. Both have 5*5 relaxed structures, representing five different central aluminium sites with five random occupations of other sites. The first set, "low", has no Al-O-Al or Si-O-Si bonds present, the other "random" has some present. We want to aggregate aluminium EFGs according to their site type (indexed by number of silicon second-neighbours) and compute statistics.

Download the files here and unzip them somewhere. Download the example script and put it in the same directory.

curl -L "http://tfgg.github.io/magres-format/workshop/tutorial_disordered/gehlenite_efgs.zip" > gehlenite_efgs.zip
unzip gehlenite_efgs.zip
curl -L "http://tfgg.github.io/magres-format/workshop/tutorial_disordered/Cq_dist.py" > Cq_dist.py

Inspect Cq_dist.py with your favourite text editor and try running it on the two sets of structures:

python Cq_dist.py files/orig
python Cq_dist.py files/low

Note the means and standard deviations for the sites. What does this tell you about the presence of Loewenstein rule violating sites and disorder?

Some key techniques used in this script:

  • We used magres.utils.find_all to quickly locate all the .magres files in the given directory and make a MagresAtoms object for each.
  • We loop over each atoms object and loop over every aluminium atom in each by using the atoms.species('Al') selector.
  • In order to classify the aluminium site as T1 or T2 (see paper), we count the number of aluminium and silicon neighbours within 3.5 Angstroms (chosen for this specific system) using the atoms.species(['Al', 'Si']).within(Al_atom, 3.5) selector.
  • We further refine this list of neighbours to just silicon atoms in order to count them for n_si using the neighbours.species('Si') selector.
  • We access the C_q value for the aluminium atom by access the efg.Cq property on it. This is available on all atoms with EFG data.


If the matplotlib Python module is installed on your computer, you could try plotting a histogram of the EFGs with the following code at the end of the script:

import matplotlib.pyplot as plt
plt.hist(Al_Cqs['T1'][0], bins=20)
plt.hist(Al_Cqs['T1'][1], bins=20)
plt.hist(Al_Cqs['T1'][2], bins=20)
plt.hist(Al_Cqs['T1'][3], bins=20)
plt.hist(Al_Cqs['T1'][4], bins=20)

Extension: Calculating structural parameters

Modify the script to calculate Ghose's tetrahedral distortion parameter (Ghose, S.; Tsang, T. Am. Mineral. 1973, 58, 748): the mean tangent of the absolute deviation of the O-Al-O bond angles from the tetrahedral angle (109.47į). Does this correlate with the aluminium Cqs?

Hint 1: You can get the positions of atoms by accessing their ".position" property. Combine this with the dot product (numpy.dot) to calculate angles.

Hint 2: Get the four neighbouring oxygen atoms with atoms.species('O').within(atom_Al, 2.5)

A working script which does this is available at http://tfgg.github.io/magres-format/workshop/tutorial_disordered/Cq_ghose.py . Run it like

python Cq_ghose files/low > ghose.dat

and use gnuplot to plot a scatterplot of tetrahedral distortion vs. aluminium Cq.