First we need to load in Python the appropriate bits from our magres.atoms
module, namely the MagresAtoms
class. This will contain our calculation result.
from __future__ import print_function
from magres.atoms import MagresAtoms
from numpy import mean
Next, we'll load our .magres
calculation output file. This could be from any code, but this one in particular is from a CASTEP calculation.
atoms = MagresAtoms.load_magres('../samples/ethanol-all.magres')
We now have our atoms loaded into the atoms
variable. This behaves like a normal list
with a bit extra. Let's check how many we have.
print("We have", len(atoms), "atoms")
We have 9 atoms
That's the correct number for ethanol and our structure, which is a single ethanol molecule in a large box, emulating a vacuum.
Let's iterate through the atoms
structure and print out the unreferenced isotropic magnetic shielding for each one. We do this by accessing the ms.iso
property on each one.
for atom in atoms:
print(atom, atom.ms.iso)
1H1 29.5926190591 1H2 30.2560510084 1H3 30.1027534841 1H4 26.9800272217 1H5 27.3904129776 1H6 31.9849757497 13C1 156.467218182 13C2 109.857140107 17O1 268.028520177
Let's set an arbitrary reference on the hydrogens and print the chemical shieldings for just those.
atoms.species('H').set_reference(10.0)
for atom in atoms.species('H'):
print(atom, atom.ms.cs)
1H1 -19.5926190591 1H2 -20.2560510084 1H3 -20.1027534841 1H4 -16.9800272217 1H5 -17.3904129776 1H6 -21.9849757497
We've used a selector over the atoms
structure, in this case the species
selector. This returns a sub list of atoms that match the given species.
There are other selectors, such as the within
selector. This returns all atoms within a given radius of a particular point and respects periodicity, i.e. you'll get multiple images of the same atom in a crystal. Here we'll find all atoms within 2 Å of the C1 atom and print out the isotropic and anisotropic magnetic shielding.
for atom in atoms.within(atoms.C1, 2.0):
print(atom, atom.ms.iso, atom.ms.aniso)
1H1 29.5926190591 8.94151483127 1H2 30.2560510084 8.18850823039 1H3 30.1027534841 7.28759653709 13C1 156.467218182 33.7970367459 13C2 109.857140107 70.2510167561
There are a number of other conventions available in addition to the isotropic and anisotropic shieldings, such as span
and skew
. Read more in the module documentation.
You can access attributes directly on collections of atoms, returning a list of that property. E.g., to get all the hydrogen atom's isotropic magnetic shieldings:
atoms.species('H').ms.iso
ListPropertyView([29.592619059099999, 30.256051008433332, 30.102753484133331, 26.980027221666671, 27.390412977566669, 31.98497574966667])
This allows us to concisely take the mean of the magnetic shieldings of all the hydrogens bonded to each of the C1, C2, and O1 atoms by using the bonded
property on each atom. This gives, if present, the collection of atoms bonded to that atom.
print("C1 H mean ms iso = ", mean(atoms.C1.bonded.species('H').ms.iso))
print("C2 H mean ms iso = ", mean(atoms.C2.bonded.species('H').ms.iso))
print("O1 H mean ms iso = ", mean(atoms.O1.bonded.species('H').ms.iso))
C1 H mean ms iso = 29.9838078506 C2 H mean ms iso = 28.7851386496 O1 H mean ms iso = 31.9849757497
You can also directly access the magnetic shielding tensor sigma
and its eigenvalues and eigenvectors.
print("Magnetic shielding tensor, sigma")
print(atoms.C1.ms.sigma)
print()
print("Eigenvectors of sigma")
print(atoms.C1.ms.evecs)
print()
print("Eigenvalues of sigma")
print(atoms.C1.ms.evals)
Magnetic shielding tensor, sigma [[ 150.36339808 -10.15515726 -0.2283921 ] [ 0.26262154 160.89621961 23.09329641] [ 7.60589797 15.67968987 158.14203685]] Eigenvectors of sigma [array([-0.42122012, -0.90608125, -0.03987944]), array([-0.62630201, 0.25878997, 0.73537307]), array([ 0.65598735, -0.33473051, 0.67648805])] Eigenvalues of sigma [137.26423106282505, 153.13884747124231, 178.99857601293263]