First, let's import the magres.atoms
library and use it to load a CASTEP calculation output. This calculation was performed on an ethanol molecule between every nucleus with the following input files (jcoupling_site
set to C 1
as an example):
%block POSITIONS_ABS
ang
H -0.869341 -0.175289 -0.649057
H 0.183454 -0.923201 0.560623
H 0.869130 0.198626 -0.628783
H -1.129705 0.975874 1.565773
H -0.437769 2.079941 0.372865
H 1.061671 0.678653 2.298066
C -0.003000 -0.003000 -0.003000
C -0.246915 1.165107 0.938396
O 0.896314 1.459074 1.742964
%endblock POSITIONS_ABS
%block SPECIES_POT
H H_POT.ASC.DAT
C C_POT.ASC.DAT
O O_POT.ASC.DAT
%endblock SPECIES_POT
%block LATTICE_CART
ang
6.000000 0.000000 0.000000
0.000000 6.000000 0.000000
0.000000 0.000000 6.000000
%endblock LATTICE_CART
jcoupling_site: C 1
kpoint_mp_grid 1 1 1
kpoint_mp_offset 0.25 0.25 0.25
opt_strategy : speed
task : magres
magres_task : jcoupling
cut_off_energy : 80 ry
xc_functional : PBE
popn_calculate : true
spin_polarized : false
max_scf_cycles : 60
nextra_bands : 20
from __future__ import print_function
from magres.atoms import MagresAtoms, MagresAtomsView
ethanol_atoms = MagresAtoms.load_magres('../samples/ethanol-all.magres')
IPython's rich representations means that if you just put a set of atoms to the IPython notebook console, the magres.atoms
library will automatically produce an image of the structure. This is generated using pydot
, which is a Python interface to the graph (as in graph theory) plotting program graphviz
.
If present, it will also plot the (unreferenced) magnetic shielding below the atom label, and plot the largest indirect nuclear spin-spin interactions (J-coupling) as lines between nuclei. Red lines mean a positive sign on the reduced isotropic coupling and blue lines a negative one.
Let's also load in bonding information from a .castep file. This is a bit of a hack right now.
castep_file = open('../samples/ethanol.castep').read()
ethanol_atoms.calculate_bonds()
This is done by the magres.atoms.MagresAtomsView
class implementing a _repr_png_
method, returning raw png data.
You can also show just a subset of atoms, e.g. all atoms bonded to C1
plus C1
itself.
(ethanol_atoms.C1 + ethanol_atoms.C1.bonded)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-3-7364da3b33a7> in <module>() ----> 1 (ethanol_atoms.C1 + ethanol_atoms.C1.bonded) /Users/christian/Dropbox/PycharmProjects/magres-format/magres/atoms.py in __radd__(self, b) 327 def __radd__(self, b): 328 if type(b) is MagresAtom: --> 329 new_atoms = set(self.atoms + [b]) 330 331 return MagresAtomsView(list(new_atoms), self.lattice) TypeError: unhashable type: 'MagresAtom'
You might use this to investigate a sub-section of a system.
Here's the J-coupling network just between hydrogen atoms.
Let's try another, larger system. From the Wikipedia page on Alanine:
Alanine (abbreviated as Ala or A) is an α-amino acid with the chemical formula CH3CH(NH2)COOH. The L-isomer is one of the 20 amino acids encoded by the genetic code. Its codons are GCU, GCC, GCA, and GCG. It is classified as a nonpolar amino acid. L-Alanine is second only to leucine in rate of occurrence, accounting for 7.8% of the primary structure in a sample of 1,150 proteins. D-Alanine occurs in bacterial cell walls and in some peptide antibiotics.
Here we look at its solid-state crystalline form.
Let's load our calculation. This only has J-coupling calculations between some of the nuclei in the system.
alanine_atoms = MagresAtoms.load_magres('../samples/alanine-jc-all.magres')
#castep_file = open('../samples/alanine.castep').read()
alanine_atoms.calculate_bonds()
Again, all we need to do print the bonding network is output the atoms object to the IPython console.
alanine_atoms
<magres.atom.MagresAtoms - 52 atoms>
You can already see a few interesting features, such as the H-bonding mediated J-coupling between nitrogen and oxygen molecules, e.g. N2 to O1.
Let's look at that more closely. We can select all atoms within 2 Angstroms of N1, N2, O1 and O3 (the two H-bridges) by doing a .within
query and adding together the sets.
alanine_atoms.N1.bonded + alanine_atoms.N2.bonded + \
alanine_atoms.O1.bonded + alanine_atoms.O3.bonded + \
alanine_atoms.N1 + alanine_atoms.N2 + alanine_atoms.O1 + alanine_atoms.O3
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-6-e6894e7d24fa> in <module>() ----> 1 alanine_atoms.N1.bonded + alanine_atoms.N2.bonded + alanine_atoms.O1.bonded + alanine_atoms.O3.bonded + alanine_atoms.N1 + alanine_atoms.N2 + alanine_atoms.O1 + alanine_atoms.O3 /Users/christian/Dropbox/PycharmProjects/magres-format/magres/atoms.py in __add__(self, b) 321 322 elif type(b) is MagresAtom: --> 323 new_atoms = set(self.atoms + [b]) 324 325 return MagresAtomsView(list(new_atoms), self.lattice) TypeError: unhashable type: 'MagresAtom'
Let's recurse along the the bonding network to pluck out an individual alanine molecule.
def get_molecule(atom):
atoms = set()
def _get_molecule(atom1):
for atom2 in atom1.bonded:
if atom2 not in atoms:
atoms.add(atom2)
_get_molecule(atom2)
_get_molecule(atom)
return MagresAtomsView(list(atoms), atom.bonded.lattice)
get_molecule(alanine_atoms.N1)
<magres.atom.MagresAtomsView - 164 atoms>