"""
This file was adapted from ASE.
Module to read and write atoms in cif file format.
See http://www.iucr.org/resources/cif/spec/version1.1/cifsyntax for a
description of the file format. STAR extensions as save frames,
global blocks, nested loops and multi-data values are not supported.
"""
import re, os
import shlex
import warnings
import numpy as np
from muesr.core.magmodel import MM
from muesr.core.spg import Spacegroup
from muesr.core.nprint import nprint, nprintmsg
from muesr.core.isstr import isstr
from muesr.i_o.cif.crystal import crystal
from muesr.core.spg import spacegroup_from_data
from muesr.i_o.cif.cell import cellpar_to_cell
# Old conventions:
old_spacegroup_names = {'Abm2': 'Aem2',
'Aba2': 'Aea2',
'Cmca': 'Cmce',
'Cmma': 'Cmme',
'Ccca': 'Ccc1'}
[docs]def load_cif(sample, filename, reset_muon=True, reset_sym=True):
"""
Loads the structural information from a Cif file.
:param muesr.core.sample.Sample sample: the sample object.
:param str filename: the cif file path (path + filename).
:param bool reset_muon: if true the muon positions is reinitialized.
:param bool reset_sym: if true the symmetry is reinitialized. ALWAYS TRUE
:returns: True if succesfull, false otherwise
:rtype: bool
"""
filename = str(filename)
f = open(os.path.expanduser(filename),'r')
#nprint ("Problems with file name...file not found!", 'warn')
#return False
#print(read_cif(f,0))
atoms, sym = read_cif(f,0) # selectd index 0
f.close()
if atoms:
sample._reset(cell=True,sym=True,magdefs=True,muon=reset_muon)
sample.cell = atoms
sample.sym = sym
return True
else:
nprint ("Atoms not loaded!", 'warn')
return False
[docs]def load_mcif(sample, filename, reset_muon=True, reset_sym=True):
"""
Loads both the crystalline structure and the magnetic order from a
mcif file.
N.B.: This function is EXPERIMENTAL.
.. note::
Only the first lattice and magnetic structure is loaded.
mCif files hosting multiple structures are not supported.
:param muesr.core.sample.Sample sample: the sample object.
:param str filename: the mcif file path (path + filename).
:param bool reset_muon: if true the muon positions is reinitialized.
:param bool reset_sym: if true the symmetry is reinitialized.
:returns: True if succesfull, false otherwise
:rtype: bool
"""
# DEFINITION OF UNITS AND SETTINGS: http://cmswiki.byu.edu/wiki/Magnetic_Coordinates
# new link http://magcryst.org/resources/magnetic-coordinates/
f = open(os.path.expanduser(filename),'r')
data = parse_cif(f)
f.close()
#print('Parsing data from: ' + data[0][0])
tags = data[0][1]
# Load symmetry
#sym = tags2spacegroup(tags)
# load cell info
a = tags['_cell_length_a']
b = tags['_cell_length_b']
c = tags['_cell_length_c']
alpha = tags['_cell_angle_alpha']
beta = tags['_cell_angle_beta']
gamma = tags['_cell_angle_gamma']
# Find magnetic atoms
mag_atoms_labels = tags['_atom_site_moment_label']
# load mag moments
mag_atoms_moments = np.zeros([len(mag_atoms_labels),3],dtype=np.complex)
scaled_positions = np.array([tags['_atom_site_fract_x'],
tags['_atom_site_fract_y'],
tags['_atom_site_fract_z']]).T
scaled_positions = np.mod(scaled_positions, 1.)
all_scaled_pos = np.copy(scaled_positions)
symbols = []
if '_atom_site_type_symbol' in tags:
labels = tags['_atom_site_type_symbol']
else:
labels = tags['_atom_site_label']
for s in labels:
# Strip off additional labeling on chemical symbols
m = re.search(r'([A-Z][a-z]?)', s)
symbol = m.group(0)
symbols.append(symbol)
fcs = np.zeros_like(all_scaled_pos,dtype=np.complex)
for i, al in enumerate(tags['_atom_site_label']):
if al in tags['_atom_site_moment_label']:
mi = tags['_atom_site_moment_label'].index(al)
fcs[i] = [complex(tags['_atom_site_moment_crystalaxis_x'][mi]),
complex(tags['_atom_site_moment_crystalaxis_y'][mi]),
complex(tags['_atom_site_moment_crystalaxis_z'][mi])
]
# THESE ARE IN CRYSTAL AXIS COORDINATE SYSTEM!!!
# bohr magneton units are used
# the magnetic metric tensor is M = L.G.L^(-1), which is unitless.
# NOW GO TO REDUCED LATTICE COORDINATE SYSTEM TO DO THE SYMMETRY
L = np.diag([1./a,1./b,1./c])
fcs = np.dot(fcs,L)
# we copy the fourier components that were present in the mcif.
# the others will be obtained from symmetry operations
all_fcs = np.copy(fcs)
for j, m_a_p in enumerate(scaled_positions):
for cent in tags['_space_group_symop.magn_centering_xyz']: # magn_centering_xyz
rc,tc,trc = parse_magn_operation_xyz_string(cent)
# apply centering and go back to the unit cell
cm_a_p = (np.dot(rc,m_a_p)+tc)%1.
#print ('cmap is :' + str(cm_a_p))
for s in tags['_space_group_symop.magn_operation_xyz']: # magn_operation_xyz
r,t,tr = parse_magn_operation_xyz_string(s)
symp = (np.dot(r,cm_a_p)+t)%1.
#print('Symp is: '+ str(symp))
# check if this position was already present
for l, pp in enumerate(all_scaled_pos):
if np.allclose(symp,pp,rtol=1e-3):
break
else:
# Append position just found with symmetry
symbols.append(symbols[j])
all_scaled_pos = np.append(all_scaled_pos,[symp],axis=0)
# go to crystal units
crysfc = trc*np.linalg.det(rc)*tr*np.linalg.det(r)*np.dot(r, np.dot(rc,fcs[j]))
all_fcs = np.append(all_fcs ,[crysfc],axis=0)
#
mag_crys2car = cellpar_to_cell([a, b, c, alpha, beta, gamma], (0,0,1), None)
# go to cartesian coordinates
cartfc = np.dot(all_fcs,mag_crys2car)
sample._reset(muon=reset_muon,sym=reset_sym)
# symmetry is already introduced when parsing magnetism
sample.cell, _ = crystal(symbols=symbols, basis=all_scaled_pos, cellpar=[a, b, c, alpha, beta, gamma])
# initialization needs the number of atoms in the unit cell
nmm=MM(len(all_scaled_pos),sample._cell.get_cell())
# magnetic moments are specified in cartesian coordinates.
# position in crystal coordinates...not nice but simple!
nmm.fc_set(cartfc)
# propagation vector is 0 since the cell "contains" the magnetic
# structure, no incommensurate structures, sorry!
nmm.k=np.array([0.,0.,0.])
sample.mm=nmm
return True
[docs]def convert_value(value):
"""Convert CIF value string to corresponding python type."""
value = value.strip()
if re.match('(".*")|(\'.*\')$', value):
return value[1:-1]
elif re.match(r'[+-]?\d+$', value):
return int(value)
elif re.match(r'[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?$', value):
return float(value)
elif re.match(r'[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?\(\d+\)$',
value):
return float(value[:value.index('(')]) # strip off uncertainties
elif re.match(r'[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?\(\d+$',
value):
warnings.warn('Badly formed number: "{0}"'.format(value))
return float(value[:value.index('(')]) # strip off uncertainties
else:
return value
[docs]def parse_multiline_string(lines, line):
"""Parse semicolon-enclosed multiline string and return it."""
assert line[0] == ';'
strings = [line[1:].lstrip()]
while True:
line = lines.pop().strip()
if line[:1] == ';':
break
strings.append(line)
return '\n'.join(strings).strip()
[docs]def parse_singletag(lines, line):
"""Parse a CIF tag (entries starting with underscore). Returns
a key-value pair."""
kv = line.split(None, 1)
if len(kv) == 1:
key = line
line = lines.pop().strip()
while not line or line[0] == '#':
line = lines.pop().strip()
if line[0] == ';':
value = parse_multiline_string(lines, line)
else:
value = line
else:
key, value = kv
return key, convert_value(value)
[docs]def parse_loop(lines):
"""Parse a CIF loop. Returns a dict with column tag names as keys
and a lists of the column content as values."""
header = []
line = lines.pop().strip()
while line.startswith('_'):
tokens = line.split()
header.append(tokens[0].lower())
if len(tokens) == 1:
line = lines.pop().strip()
else:
line = ' '.join(tokens[1:])
break
columns = dict([(h, []) for h in header])
if len(columns) != len(header):
seen = set()
dublicates = [h for h in header if h in seen or seen.add(h)]
warnings.warn('Duplicated loop tags: {0}'.format(dublicates))
tokens = []
while True:
lowerline = line.lower()
if (not line or
line.startswith('_') or
lowerline.startswith('data_') or
lowerline.startswith('loop_')):
break
if line.startswith('#'):
line = lines.pop().strip()
continue
if line.startswith(';'):
t = [parse_multiline_string(lines, line)]
else:
if len(header) == 1:
t = [line]
else:
t = shlex.split(line, posix=False)
line = lines.pop().strip()
tokens.extend(t)
if len(tokens) < len(columns):
continue
if len(tokens) == len(header):
for h, t in zip(header, tokens):
columns[h].append(convert_value(t))
else:
warnings.warn('Wrong number of tokens: {0}'.format(tokens))
tokens = []
if line:
lines.append(line)
return columns
[docs]def parse_items(lines, line):
"""Parse a CIF data items and return a dict with all tags."""
tags = {}
while True:
if not lines:
break
line = lines.pop()
if not line:
break
line = line.strip()
lowerline = line.lower()
if not line or line.startswith('#'):
continue
elif line.startswith('_'):
key, value = parse_singletag(lines, line)
tags[key.lower()] = value
elif lowerline.startswith('loop_'):
tags.update(parse_loop(lines))
elif lowerline.startswith('data_'):
if line:
lines.append(line)
break
elif line.startswith(';'):
parse_multiline_string(lines, line)
else:
raise ValueError('Unexpected CIF file entry: "{0}"'.format(line))
return tags
[docs]def parse_block(lines, line):
"""Parse a CIF data block and return a tuple with the block name
and a dict with all tags."""
assert line.lower().startswith('data_')
blockname = line.split('_', 1)[1].rstrip()
tags = parse_items(lines, line)
return blockname, tags
[docs]def parse_cif(fileobj):
"""Parse a CIF file. Returns a list of blockname and tag
pairs. All tag names are converted to lower case."""
if isstr(fileobj):
fileobj = open(fileobj)
lines = [''] + fileobj.readlines()[::-1] # all lines (reversed)
blocks = []
while True:
if not lines:
break
line = lines.pop()
line = line.strip()
if not line or line.startswith('#'):
continue
blocks.append(parse_block(lines, line))
return blocks
[docs]def read_cif(fileobj, index, store_tags=False, primitive_cell=False,
subtrans_included=True):
"""Read Atoms object from CIF file. *index* specifies the data
block number or name (if string) to return.
If *index* is None or a slice object, a list of atoms objects will
be returned. In the case of *index* is *None* or *slice(None)*,
only blocks with valid crystal data will be included.
If *store_tags* is true, the *info* attribute of the returned
Atoms object will be populated with all tags in the corresponding
cif data block.
If *primitive_cell* is true, the primitive cell will be built instead
of the conventional cell.
If *subtrans_included* is true, sublattice translations are
assumed to be included among the symmetry operations listed in the
CIF file (seems to be the common behaviour of CIF files).
Otherwise the sublattice translations are determined from setting
1 of the extracted space group. A result of setting this flag to
true, is that it will not be possible to determine the primitive
cell.
"""
blocks = parse_cif(fileobj)
# Find all CIF blocks with valid crystal data
images = []
for name, tags in blocks:
try:
atoms, spg = tags2atoms(tags, store_tags, primitive_cell,
subtrans_included)
images.append([atoms,spg])
except KeyError:
pass
for data in images[index]:
yield data
[docs]def write_cif(fileobj, images, format='default'):
"""Write *images* to CIF file."""
if isstr(fileobj):
fileobj = paropen(fileobj, 'w')
if hasattr(images, 'get_positions'):
images = [images]
for i, atoms in enumerate(images):
fileobj.write('data_image%d\n' % i)
a, b, c, alpha, beta, gamma = atoms.get_cell_lengths_and_angles()
if format == 'mp':
comp_name = atoms.get_chemical_formula(mode='reduce')
sf = split_chem_form(comp_name)
formula_sum = ''
ii = 0
while ii < len(sf):
formula_sum = formula_sum + ' ' + sf[ii] + sf[ii + 1]
ii = ii + 2
formula_sum = str(formula_sum)
fileobj.write('_chemical_formula_structural %s\n' %
atoms.get_chemical_formula(mode='reduce'))
fileobj.write('_chemical_formula_sum "%s"\n' % formula_sum)
fileobj.write('_cell_length_a %g\n' % a)
fileobj.write('_cell_length_b %g\n' % b)
fileobj.write('_cell_length_c %g\n' % c)
fileobj.write('_cell_angle_alpha %g\n' % alpha)
fileobj.write('_cell_angle_beta %g\n' % beta)
fileobj.write('_cell_angle_gamma %g\n' % gamma)
fileobj.write('\n')
if atoms.pbc.all():
fileobj.write('_symmetry_space_group_name_H-M %s\n' % '"P 1"')
fileobj.write('_symmetry_int_tables_number %d\n' % 1)
fileobj.write('\n')
fileobj.write('loop_\n')
fileobj.write(' _symmetry_equiv_pos_as_xyz\n')
fileobj.write(" 'x, y, z'\n")
fileobj.write('\n')
fileobj.write('loop_\n')
if format == 'mp':
fileobj.write(' _atom_site_type_symbol\n')
fileobj.write(' _atom_site_label\n')
fileobj.write(' _atom_site_symmetry_multiplicity\n')
fileobj.write(' _atom_site_fract_x\n')
fileobj.write(' _atom_site_fract_y\n')
fileobj.write(' _atom_site_fract_z\n')
fileobj.write(' _atom_site_occupancy\n')
else:
fileobj.write(' _atom_site_label\n')
fileobj.write(' _atom_site_occupancy\n')
fileobj.write(' _atom_site_fract_x\n')
fileobj.write(' _atom_site_fract_y\n')
fileobj.write(' _atom_site_fract_z\n')
fileobj.write(' _atom_site_thermal_displace_type\n')
fileobj.write(' _atom_site_B_iso_or_equiv\n')
fileobj.write(' _atom_site_type_symbol\n')
scaled = atoms.get_scaled_positions()
no = {}
for i, atom in enumerate(atoms):
symbol = atom.symbol
if symbol in no:
no[symbol] += 1
else:
no[symbol] = 1
if format == 'mp':
fileobj.write(
' %-2s %4s %4s %7.5f %7.5f %7.5f %6.1f\n' %
(symbol, symbol + str(no[symbol]), 1,
scaled[i][0], scaled[i][1], scaled[i][2], 1.0))
else:
fileobj.write(
' %-8s %6.4f %7.5f %7.5f %7.5f %4s %6.3f %s\n' %
('%s%d' % (symbol, no[symbol]),
1.0,
scaled[i][0],
scaled[i][1],
scaled[i][2],
'Biso',
1.0,
symbol))
#### Addition for magnetic CIF files ####
[docs]def parse_magn_operation_xyz_string(tag):
""" Parse the symmetry operations of the magnetic part of mcif """
# this could probably be used also for cif
r1,r2,r3,p=tag.split(',')
r = np.zeros([3,3]) # rotations
t = np.zeros([3]) # translations
# compile traslation matrix
t[0] = convert_to_float(r1[-4:]) if '/' in r1 else 0.
t[1] = convert_to_float(r2[-4:]) if '/' in r2 else 0.
t[2] = convert_to_float(r3[-4:]) if '/' in r3 else 0.
# compile rotation matrix
for i, line in enumerate([r1,r2,r3]):
x=y=z=0.
m = re.search(r'-?\d?x',line)
if m:
xs=m.group()
x = 1. if ('x' in xs and not '-' in xs) else -1.
m = re.search(r'\d',xs)
if m:
x *= float(m.group())
m = re.search(r'-?\d?y',line)
if m:
ys=m.group()
y = 1. if ('y' in ys and not '-' in ys) else -1.
m = re.search(r'\d',ys)
if m:
y *= float(m.group())
m = re.search(r'-?\d?z',line)
if m:
zs=m.group()
z = 1. if ('z' in zs and not '-' in zs) else -1.
m = re.search(r'\d',zs)
if m:
z *= float(m.group())
r[i,0] = x
r[i,1] = y
r[i,2] = z
return (r,t,float(p))
[docs]def convert_to_float(frac_str):
"This function converts fractions to float, ex. -1/3 = -0.33333..."
try:
return float(frac_str)
except ValueError:
try:
num, denom = frac_str.split('/')
except ValueError:
return None
try:
leading, num = num.split(' ')
except ValueError:
return float(num) / float(denom)
if float(leading) < 0:
sign_mult = -1
else:
sign_mult = 1
return float(leading) + sign_mult * (float(num) / float(denom))