"""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.
The "latin-1" encoding is required by the IUCR specification.
"""
import re
import shlex
import warnings
import numpy as np
from ase import Atoms
from ase.parallel import paropen
from ase.spacegroup import crystal
from ase.spacegroup.spacegroup import spacegroup_from_data, Spacegroup
from ase.utils import basestring
from ase.data import atomic_numbers, atomic_masses
from ase.io.cif_unicode import format_unicode, handle_subscripts
# Old conventions:
old_spacegroup_names = {'Abm2': 'Aem2',
'Aba2': 'Aea2',
'Cmca': 'Cmce',
'Cmma': 'Cmme',
'Ccca': 'Ccc1'}
def convert_value(value):
"""Convert CIF value string to corresponding python type."""
value = value.strip()
if re.match('(".*")|(\'.*\')$', value):
return handle_subscripts(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 handle_subscripts(value)
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()
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)
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
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
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
def parse_cif(fileobj, reader='ase'):
"""Parse a CIF file. Returns a list of blockname and tag
pairs. All tag names are converted to lower case."""
if reader == 'ase':
return parse_cif_ase(fileobj)
elif reader == 'pycodcif':
return parse_cif_pycodcif(fileobj)
def parse_cif_ase(fileobj):
"""Parse a CIF file using ase CIF parser"""
blocks = []
if isinstance(fileobj, basestring):
fileobj = open(fileobj, 'rb')
data = fileobj.read()
if isinstance(data, bytes):
data = data.decode('latin1')
data = format_unicode(data)
data = [e for e in data.split('\n') if len(e) > 0]
if len(data) > 0 and data[0].rstrip() == '#\\#CIF_2.0':
warnings.warn('CIF v2.0 file format detected; `ase` CIF reader might '
'incorrectly interpret some syntax constructions, use '
'`pycodcif` reader instead')
lines = [''] + data[::-1] # all lines (reversed)
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
def parse_cif_pycodcif(fileobj):
"""Parse a CIF file using pycodcif CIF parser"""
blocks = []
if not isinstance(fileobj, basestring):
fileobj = fileobj.name
try:
from pycodcif import parse
except ImportError:
raise ImportError(
'parse_cif_pycodcif requires pycodcif ' +
'(http://wiki.crystallography.net/cod-tools/pycodcif/)')
data,_,_ = parse(fileobj)
for datablock in data:
tags = datablock['values']
for tag in tags.keys():
values = [convert_value(x) for x in tags[tag]]
if len(values) == 1:
tags[tag] = values[0]
else:
tags[tag] = values
blocks.append((datablock['name'], tags))
return blocks
def tags2atoms(tags, store_tags=False, primitive_cell=False,
subtrans_included=True, fractional_occupancies=True):
"""Returns an Atoms object from a cif tags dictionary. See read_cif()
for a description of the arguments."""
if primitive_cell and subtrans_included:
raise RuntimeError(
'Primitive cell cannot be determined when sublattice translations '
'are included in the symmetry operations listed in the CIF file, '
'i.e. when `subtrans_included` is True.')
cell_tags = ['_cell_length_a', '_cell_length_b', '_cell_length_c',
'_cell_angle_alpha', '_cell_angle_beta', '_cell_angle_gamma']
# If any value is missing, ditch periodic boundary conditions
has_pbc = True
try:
cell_values = [tags[ct] for ct in cell_tags]
a, b, c, alpha, beta, gamma = cell_values
except KeyError:
has_pbc = False
# Now get positions
try:
scaled_positions = np.array([tags['_atom_site_fract_x'],
tags['_atom_site_fract_y'],
tags['_atom_site_fract_z']]).T
except KeyError:
scaled_positions = None
try:
positions = np.array([tags['_atom_site_cartn_x'],
tags['_atom_site_cartn_y'],
tags['_atom_site_cartn_z']]).T
except KeyError:
positions = None
if (positions is None) and (scaled_positions is None):
raise RuntimeError('No positions found in structure')
elif scaled_positions is not None and not has_pbc:
raise RuntimeError('Structure has fractional coordinates but not '
'lattice parameters')
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)
# Symmetry specification, see
# http://www.iucr.org/resources/cif/dictionaries/cif_sym for a
# complete list of official keys. In addition we also try to
# support some commonly used depricated notations
no = None
if '_space_group.it_number' in tags:
no = tags['_space_group.it_number']
elif '_space_group_it_number' in tags:
no = tags['_space_group_it_number']
elif '_symmetry_int_tables_number' in tags:
no = tags['_symmetry_int_tables_number']
symbolHM = None
if '_space_group.Patterson_name_h-m' in tags:
symbolHM = tags['_space_group.patterson_name_h-m']
elif '_symmetry_space_group_name_h-m' in tags:
symbolHM = tags['_symmetry_space_group_name_h-m']
elif '_space_group_name_h-m_alt' in tags:
symbolHM = tags['_space_group_name_h-m_alt']
if symbolHM is not None:
symbolHM = old_spacegroup_names.get(symbolHM.strip(), symbolHM)
for name in ['_space_group_symop_operation_xyz',
'_space_group_symop.operation_xyz',
'_symmetry_equiv_pos_as_xyz']:
if name in tags:
sitesym = tags[name]
break
else:
sitesym = None
# The setting needs to be passed as either 1 or two, not None (default)
setting = 1
spacegroup = 1
if sitesym is not None:
if isinstance(sitesym, str):
sitesym = [sitesym]
subtrans = [(0.0, 0.0, 0.0)] if subtrans_included else None
spacegroup = spacegroup_from_data(
no=no, symbol=symbolHM, sitesym=sitesym, subtrans=subtrans,
setting=setting)
elif no is not None:
spacegroup = no
elif symbolHM is not None:
spacegroup = symbolHM
else:
spacegroup = 1
kwargs = {}
if store_tags:
kwargs['info'] = tags.copy()
if 'D' in symbols:
deuterium = [symbol == 'D' for symbol in symbols]
symbols = [symbol if symbol != 'D' else 'H' for symbol in symbols]
else:
deuterium = False
setting_name = None
if '_space_group_crystal_system' in tags:
setting_name = tags['_space_group_crystal_system']
elif '_symmetry_cell_setting' in tags:
setting_name = tags['_symmetry_cell_setting']
if setting_name:
no = Spacegroup(spacegroup).no
# rhombohedral systems
if no in (146, 148, 155, 160, 161, 166, 167):
if setting_name == 'hexagonal':
setting = 1
elif setting_name in ('trigonal', 'rhombohedral'):
setting = 2
else:
warnings.warn(
'unexpected crystal system %r for space group %r' % (
setting_name, spacegroup))
# FIXME - check for more crystal systems...
else:
warnings.warn(
'crystal system %r is not interpreted for space group %r. '
'This may result in wrong setting!' % (
setting_name, spacegroup))
occupancies = None
if fractional_occupancies:
try:
occupancies = tags['_atom_site_occupancy']
# no warnings in this case
kwargs['onduplicates'] = 'keep'
except KeyError:
pass
else:
try:
if not np.allclose(tags['_atom_site_occupancy'], 1.):
warnings.warn(
'Cif file containes mixed/fractional occupancies. '
'Consider using `fractional_occupancies=True`')
kwargs['onduplicates'] = 'keep'
except KeyError:
pass
if has_pbc:
if scaled_positions is None:
_ = Atoms(symbols, positions=positions,
cell=[a, b, c, alpha, beta, gamma])
scaled_positions = _.get_scaled_positions()
if deuterium:
numbers = np.array([atomic_numbers[s] for s in symbols])
masses = atomic_masses[numbers]
masses[deuterium] = 2.01355
kwargs['masses'] = masses
atoms = crystal(symbols, basis=scaled_positions,
cellpar=[a, b, c, alpha, beta, gamma],
spacegroup=spacegroup,
occupancies=occupancies,
setting=setting,
primitive_cell=primitive_cell,
**kwargs)
else:
atoms = Atoms(symbols, positions=positions,
info=kwargs.get('info', None))
if occupancies is not None:
# Compile an occupancies dictionary
occ_dict = {}
for i, sym in enumerate(symbols):
occ_dict[i] = {sym: occupancies[i]}
atoms.info['occupancy'] = occ_dict
if deuterium:
masses = atoms.get_masses()
masses[atoms.numbers == 1] = 1.00783
masses[deuterium] = 2.01355
atoms.set_masses(masses)
return atoms
[docs]def read_cif(fileobj, index, store_tags=False, primitive_cell=False,
subtrans_included=True, fractional_occupancies=True,
reader='ase'):
"""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.
If *fractional_occupancies* is true, the resulting atoms object will be
tagged equipped with an array `occupancy`. Also, in case of mixed
occupancies, the atom's chemical symbol will be that of the most dominant
species.
String *reader* is used to select CIF reader. Value `ase` selects
built-in CIF reader (default), while `pycodcif` selects CIF reader based
on `pycodcif` package.
"""
blocks = parse_cif(fileobj, reader)
# Find all CIF blocks with valid crystal data
images = []
for name, tags in blocks:
try:
atoms = tags2atoms(tags, store_tags, primitive_cell,
subtrans_included,
fractional_occupancies=fractional_occupancies)
images.append(atoms)
except KeyError:
pass
for atoms in images[index]:
yield atoms
def split_chem_form(comp_name):
"""Returns e.g. AB2 as ['A', '1', 'B', '2']"""
split_form = re.findall(r'[A-Z][a-z]*|\d+',
re.sub(r'[A-Z][a-z]*(?![\da-z])',
r'\g<0>1', comp_name))
return split_form
def write_enc(fileobj, s):
"""Write string in latin-1 encoding."""
fileobj.write(s.encode("latin-1"))
[docs]def write_cif(fileobj, images, format='default'):
"""Write *images* to CIF file."""
if isinstance(fileobj, basestring):
fileobj = paropen(fileobj, 'wb')
if hasattr(images, 'get_positions'):
images = [images]
for i, atoms in enumerate(images):
write_enc(fileobj, '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)
write_enc(fileobj, '_chemical_formula_structural %s\n' %
atoms.get_chemical_formula(mode='reduce'))
write_enc(fileobj, '_chemical_formula_sum "%s"\n' %
formula_sum)
# Do this only if there's three non-zero lattice vectors
if atoms.number_of_lattice_vectors == 3:
write_enc(fileobj, '_cell_length_a %g\n' % a)
write_enc(fileobj, '_cell_length_b %g\n' % b)
write_enc(fileobj, '_cell_length_c %g\n' % c)
write_enc(fileobj, '_cell_angle_alpha %g\n' % alpha)
write_enc(fileobj, '_cell_angle_beta %g\n' % beta)
write_enc(fileobj, '_cell_angle_gamma %g\n' % gamma)
write_enc(fileobj, '\n')
write_enc(fileobj, '_symmetry_space_group_name_H-M %s\n' %
'"P 1"')
write_enc(fileobj, '_symmetry_int_tables_number %d\n' % 1)
write_enc(fileobj, '\n')
write_enc(fileobj, 'loop_\n')
write_enc(fileobj, ' _symmetry_equiv_pos_as_xyz\n')
write_enc(fileobj, " 'x, y, z'\n")
write_enc(fileobj, '\n')
write_enc(fileobj, 'loop_\n')
# Is it a periodic system?
coord_type = 'fract' if atoms.pbc.all() else 'Cartn'
if format == 'mp':
write_enc(fileobj, ' _atom_site_type_symbol\n')
write_enc(fileobj, ' _atom_site_label\n')
write_enc(fileobj, ' _atom_site_symmetry_multiplicity\n')
write_enc(fileobj, ' _atom_site_{0}_x\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_{0}_y\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_{0}_z\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_occupancy\n')
else:
write_enc(fileobj, ' _atom_site_label\n')
write_enc(fileobj, ' _atom_site_occupancy\n')
write_enc(fileobj, ' _atom_site_{0}_x\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_{0}_y\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_{0}_z\n'.format(coord_type))
write_enc(fileobj, ' _atom_site_thermal_displace_type\n')
write_enc(fileobj, ' _atom_site_B_iso_or_equiv\n')
write_enc(fileobj, ' _atom_site_type_symbol\n')
if coord_type == 'fract':
coords = atoms.get_scaled_positions().tolist()
else:
coords = atoms.get_positions().tolist()
symbols = atoms.get_chemical_symbols()
occupancies = [1 for i in range(len(symbols))]
# try to fetch occupancies // rely on the tag - occupancy mapping
try:
occ_info = atoms.info['occupancy']
for i, tag in enumerate(atoms.get_tags()):
occupancies[i] = occ_info[tag][symbols[i]]
# extend the positions array in case of mixed occupancy
for sym, occ in occ_info[tag].items():
if sym != symbols[i]:
symbols.append(sym)
coords.append(coords[i])
occupancies.append(occ)
except KeyError:
pass
no = {}
for symbol, pos, occ in zip(symbols, coords, occupancies):
if symbol in no:
no[symbol] += 1
else:
no[symbol] = 1
if format == 'mp':
write_enc(fileobj,
' %-2s %4s %4s %7.5f %7.5f %7.5f %6.1f\n' %
(symbol, symbol + str(no[symbol]), 1,
pos[0], pos[1], pos[2], occ))
else:
write_enc(fileobj,
' %-8s %6.4f %7.5f %7.5f %7.5f %4s %6.3f %s\n'
% ('%s%d' % (symbol, no[symbol]),
occ,
pos[0],
pos[1],
pos[2],
'Biso',
1.0,
symbol))