"""Modules for handling crystallographic Spacegroup related operations."""
from functools import reduce
from jarvis.core.lattice import Lattice
from jarvis.core.atoms import Atoms
import spglib
from jarvis.core.specie import Specie
import numpy as np
from numpy import sin, cos
import itertools
# from fractions import gcd
from jarvis.core.utils import gcd
# from numpy import gcd
# from math import gcd
import os
[docs]def unique_rows_2(a):
"""Remove duplicate rows."""
order = np.lexsort(a.T)
a = a[order]
diff = np.diff(a, axis=0)
ui = np.ones(len(a), "bool")
ui[1:] = (diff != 0).any(axis=1)
return a[ui]
[docs]def symmetrically_distinct_miller_indices(max_index=3, cvn_atoms=None):
"""Get unique miller indices for max_index."""
# Need to work on this
r1 = list(range(1, max_index + 1))
r2 = list(range(-max_index, 1))
r2.reverse()
r3 = r1 + r2
# r.reverse()
r = r3
# print ('sorted',r,sorted(r))
conv_hkl_list = [
miller
for miller in itertools.product(r, r, r)
if any([i != 0 for i in miller])
]
spg = Spacegroup3D(cvn_atoms)._dataset
rot = spg["rotations"]
done = []
unique_millers = []
# print (conv_hkl_list)
# print (sorted(conv_hkl_list, key=lambda x: x[0], reverse=True))
for i in conv_hkl_list:
d = abs(reduce(gcd, i))
miller = tuple([int(k / d) for k in i])
for j in rot:
prod = list(np.dot(miller, j))
if prod not in done:
unique_millers.append(i)
done.append(prod)
uniq = unique_rows_2(np.array(unique_millers))
return uniq
wyckoff_file = str(os.path.join(os.path.dirname(__file__), "Wyckoff.csv"))
[docs]def parse_wyckoff_csv(wyckoff_file):
"""Parse Wyckoff.csv from spglib.
There are 530 data sets. For one example:
9:C 1 2 1:::::::
::4:c:1:(x,y,z):(-x,y,-z)::
::2:b:2:(0,y,1/2):::
::2:a:2:(0,y,0):::
"""
rowdata = []
points = []
hP_nums = [433, 436, 444, 450, 452, 458, 460]
for i, line in enumerate(wyckoff_file):
if line.strip() == "end of data":
break
rowdata.append(line.strip().split(":"))
# 2:P -1 ::::::: <-- store line number if first element is number
if rowdata[-1][0].isdigit():
points.append(i)
points.append(i)
wyckoff = []
for i in range(len(points) - 1): # 0 to 529
symbol = rowdata[points[i]][1] # e.g. "C 1 2 1"
if i + 1 in hP_nums:
symbol = symbol.replace("R", "H", 1)
wyckoff.append({"symbol": symbol.strip()})
# When the number of positions is larger than 4,
# the positions are written in the next line.
# So those positions are connected.
for i in range(len(points) - 1):
count = 0
wyckoff[i]["wyckoff"] = []
for j in range(points[i] + 1, points[i + 1]):
# Hook if the third element is a number (multiplicity), e.g.,
#
# 232:P 2/b 2/m 2/b::::::: <- ignored
# ::8:r:1:(x,y,z):(-x,y,-z):(x,-y+1/2,-z):(-x,-y+1/2,z)
# :::::(-x,-y,-z):(x,-y,z):(-x,y+1/2,z):(x,y+1/2,-z) <- ignored
# ::4:q:..m:(x,0,z):(-x,0,-z):(x,1/2,-z):(-x,1/2,z)
# ::4:p:..2:(0,y,1/2):(0,-y+1/2,1/2):(0,-y,1/2):(0,y+1/2,1/2)
# ::4:o:..2:(1/2,y,0):(1/2,-y+1/2,0):(1/2,-y,0):(1/2,y+1/2,0)
# ...
if rowdata[j][2].isdigit():
pos = []
w = {
"letter": rowdata[j][3].strip(),
"multiplicity": int(rowdata[j][2]),
"site_symmetry": rowdata[j][4].strip(),
"positions": pos,
}
wyckoff[i]["wyckoff"].append(w)
for k in range(4):
if rowdata[j][k + 5]: # check if '(x,y,z)' or ''
count += 1
pos.append(rowdata[j][k + 5])
else:
for k in range(4):
if rowdata[j][k + 5]:
count += 1
pos.append(rowdata[j][k + 5])
# assertion
# for w in wyckoff[i]['wyckoff']:
# n_pos = len(w['positions'])
# n_pos *= len(lattice_symbols[wyckoff[i]['symbol'][0]])
# assert n_pos == w['multiplicity']
return wyckoff
[docs]def read_wyckoff_csv(filename):
"""Read wyckoff_csv file."""
with open(filename) as wyckoff_file:
return parse_wyckoff_csv(wyckoff_file)
[docs]def get_wyckoff_position_operators(hall_number):
"""Get all Wyckoff operations for Hall number."""
wyckoff = read_wyckoff_csv(wyckoff_file)
operations = wyckoff[hall_number - 1]
return operations
[docs]class Spacegroup3D(object):
"""
Provide spacegroup related data for Atoms object.
Currently uses spglib to derive spacegroup
related information for 3D materials mainly
"""
def __init__(self, atoms=[], dataset={}, symprec=1e-2, angle_tolerance=5):
"""
Following information are needed for Spacegroup3D.
If dataset is not provided, the default dataset is used.
Args:
atoms: jarvis.core.Atoms
dataset: spacegroup dataset
symprec: symmetry precision
angle_tolerance: angle tolerance
"""
self._dataset = dataset
self._atoms = atoms
self._symprec = symprec
self._angle_tolerance = angle_tolerance
if self._dataset == {}:
spg = self.spacegroup_data()
self._dataset = spg._dataset
[docs] def spacegroup_data(self):
"""Provide spacegroup data from spglib."""
phonopy_atoms = (
self._atoms.lattice_mat,
self._atoms.frac_coords,
self._atoms.atomic_numbers,
)
dataset = spglib.get_symmetry_dataset(
phonopy_atoms,
symprec=self._symprec,
angle_tolerance=self._angle_tolerance,
)
"""
keys = ('number',
'hall_number',
'international',
'hall',
'choice',
'transformation_matrix',
'origin_shift',
'rotations',
'translations',
'wyckoffs',
'site_symmetry_symbols',
'equivalent_atoms',
'mapping_to_primitive',
'std_lattice',
'std_types',
'std_positions',
'std_rotation_matrix',
'std_mapping_to_primitive',
'pointgroup')
"""
return Spacegroup3D(
atoms=self._atoms,
symprec=self._symprec,
angle_tolerance=self._angle_tolerance,
dataset=dataset,
)
@property
def space_group_symbol(self):
"""Get spacegroup symbol."""
# spg = self.spacegroup_data()
return self._dataset["international"]
@property
def space_group_number(self):
"""Get spacegroup number."""
# spg = self.spacegroup_data()
return self._dataset["number"]
@property
def primitive_atoms(self):
"""Get primitive atoms."""
phonopy_atoms = (
self._atoms.lattice_mat,
self._atoms.frac_coords,
self._atoms.atomic_numbers,
)
lattice, scaled_positions, numbers = spglib.find_primitive(
phonopy_atoms, symprec=self._symprec
)
elements = self._atoms.elements
el_dict = {}
for i in elements:
el_dict.setdefault(Specie(i).Z, i)
prim_elements = [el_dict[i] for i in numbers]
prim_atoms = Atoms(
lattice_mat=lattice,
elements=prim_elements,
coords=scaled_positions,
cartesian=False,
)
return prim_atoms
@property
def refined_atoms(self):
"""Refine atoms based on spacegroup data."""
phonopy_atoms = (
self._atoms.lattice_mat,
self._atoms.frac_coords,
self._atoms.atomic_numbers,
)
lattice, scaled_positions, numbers = spglib.refine_cell(
phonopy_atoms, self._symprec, self._angle_tolerance
)
elements = self._atoms.elements
el_dict = {}
for i in elements:
el_dict.setdefault(Specie(i).Z, i)
ref_elements = [el_dict[i] for i in numbers]
ref_atoms = Atoms(
lattice_mat=lattice,
elements=ref_elements,
coords=scaled_positions,
cartesian=False,
)
return ref_atoms
@property
def crystal_system(self):
"""Get crystal system."""
n = self._dataset["number"]
def f(i, j):
return i <= n <= j
cs = {
"triclinic": (1, 2),
"monoclinic": (3, 15),
"orthorhombic": (16, 74),
"tetragonal": (75, 142),
"trigonal": (143, 167),
"hexagonal": (168, 194),
"cubic": (195, 230),
}
crystal_sytem = None
for k, v in cs.items():
if f(*v):
crystal_sytem = k
break
return crystal_sytem
@property
def lattice_system(self):
"""Get lattice system."""
n = self._dataset["number"]
system = self.crystal_system
if n in [146, 148, 155, 160, 161, 166, 167]:
return "rhombohedral"
elif system == "trigonal":
return "hexagonal"
else:
return system
@property
def point_group_symbol(self):
"""Get pointgroup."""
return self._dataset["pointgroup"]
@property
def conventional_standard_structure(
self, tol=1e-5, international_monoclinic=True
):
"""
Give a conventional cell according to certain conventions.
The conventionss are defined in Setyawan, W., & Curtarolo,
S. (2010). High-throughput electronic band structure calculations:
Challenges and tools. Computational Materials Science,
49(2), 299-312. doi:10.1016/j.commatsci.2010.05.010
They basically enforce as much as possible
norm(a1)<norm(a2)<norm(a3)
Returns:
The structure in a conventional standardized cell
"""
struct = self.refined_atoms
latt = struct.lattice
latt_type = self.lattice_system
sorted_lengths = sorted(latt.abc)
sorted_dic = sorted(
[
{"vec": latt.matrix[i], "length": latt.abc[i], "orig_index": i}
for i in [0, 1, 2]
],
key=lambda k: k["length"],
)
if latt_type in ("orthorhombic", "cubic"):
# you want to keep the c axis where it is
# to keep the C- settings
transf = np.zeros(shape=(3, 3))
if self.space_group_symbol.startswith("C"):
transf[2] = [0, 0, 1]
a, b = sorted(latt.abc[:2])
sorted_dic = sorted(
[
{
"vec": latt.matrix[i],
"length": latt.abc[i],
"orig_index": i,
}
for i in [0, 1]
],
key=lambda k: k["length"],
)
for i in range(2):
transf[i][sorted_dic[i]["orig_index"]] = 1
c = latt.abc[2]
elif self.space_group_symbol.startswith(
"A"
): # change to C-centering to match Setyawan/Curtarolo convention
transf[2] = [1, 0, 0]
a, b = sorted(latt.abc[1:])
sorted_dic = sorted(
[
{
"vec": latt.matrix[i],
"length": latt.abc[i],
"orig_index": i,
}
for i in [1, 2]
],
key=lambda k: k["length"],
)
for i in range(2):
transf[i][sorted_dic[i]["orig_index"]] = 1
c = latt.abc[0]
else:
for i in range(len(sorted_dic)):
transf[i][sorted_dic[i]["orig_index"]] = 1
a, b, c = sorted_lengths
latt = Lattice.orthorhombic(a, b, c)
elif latt_type == "tetragonal":
# find the "a" vectors
# it is basically the vector repeated two times
transf = np.zeros(shape=(3, 3))
a, b, c = sorted_lengths
for d in range(len(sorted_dic)):
transf[d][sorted_dic[d]["orig_index"]] = 1
if abs(b - c) < tol and abs(a - c) > tol:
a, c = c, a
transf = np.dot([[0, 0, 1], [0, 1, 0], [1, 0, 0]], transf)
latt = Lattice.tetragonal(a, c)
elif latt_type in ("hexagonal", "rhombohedral"):
# for the conventional cell representation,
# we allways show the rhombohedral lattices as hexagonal
# check first if we have the refined structure shows a rhombohedral
# cell
# if so, make a supercell
a, b, c = latt.abc
if np.all(np.abs([a - b, c - b, a - c]) < 0.001):
struct.make_supercell(((1, -1, 0), (0, 1, -1), (1, 1, 1)))
a, b, c = sorted(struct.lattice.abc)
if abs(b - c) < 0.001:
a, c = c, a
new_matrix = [
[a / 2, -a * np.sqrt(3) / 2, 0],
[a / 2, a * np.sqrt(3) / 2, 0],
[0, 0, c],
]
latt = Lattice(new_matrix)
transf = np.eye(3, 3)
elif latt_type == "monoclinic":
# You want to keep the c axis where it is to keep the C- settings
if self.space_group_symbol.startswith("C"):
transf = np.zeros(shape=(3, 3))
transf[2] = [0, 0, 1]
sorted_dic = sorted(
[
{
"vec": latt.matrix[i],
"length": latt.abc[i],
"orig_index": i,
}
for i in [0, 1]
],
key=lambda k: k["length"],
)
a = sorted_dic[0]["length"]
b = sorted_dic[1]["length"]
c = latt.abc[2]
new_matrix = None
for t in itertools.permutations(list(range(2)), 2):
m = latt.matrix
latt2 = Lattice([m[t[0]], m[t[1]], m[2]])
lengths = latt2.abc
angles = latt2.angles
if angles[0] > 90:
# if the angle is > 90 we invert a and b to get
# an angle < 90
a, b, c, alpha, beta, gamma = Lattice(
[-m[t[0]], -m[t[1]], m[2]]
).parameters
transf = np.zeros(shape=(3, 3))
transf[0][t[0]] = -1
transf[1][t[1]] = -1
transf[2][2] = 1
alpha = np.pi * alpha / 180
new_matrix = [
[a, 0, 0],
[0, b, 0],
[0, c * cos(alpha), c * sin(alpha)],
]
continue
elif angles[0] < 90:
transf = np.zeros(shape=(3, 3))
# print ('464-470')
transf[0][t[0]] = 1
transf[1][t[1]] = 1
transf[2][2] = 1
a, b, c = lengths
alpha = np.pi * angles[0] / 180
new_matrix = [
[a, 0, 0],
[0, b, 0],
[0, c * cos(alpha), c * sin(alpha)],
]
if new_matrix is None:
# print ('479-482')
# this if is to treat the case
# where alpha==90 (but we still have a monoclinic sg
new_matrix = [[a, 0, 0], [0, b, 0], [0, 0, c]]
transf = np.zeros(shape=(3, 3))
for c in range(len(sorted_dic)):
transf[c][sorted_dic[c]["orig_index"]] = 1
# if not C-setting
else:
# try all permutations of the axis
# keep the ones with the non-90 angle=alpha
# and b<c
new_matrix = None
for t in itertools.permutations(list(range(3)), 3):
m = latt.matrix
a, b, c, alpha, beta, gamma = Lattice(
[m[t[0]], m[t[1]], m[t[2]]]
).parameters
if alpha > 90 and b < c:
a, b, c, alpha, beta, gamma = Lattice(
[-m[t[0]], -m[t[1]], m[t[2]]]
).parameters
transf = np.zeros(shape=(3, 3))
transf[0][t[0]] = -1
transf[1][t[1]] = -1
transf[2][t[2]] = 1
alpha = np.pi * alpha / 180
new_matrix = [
[a, 0, 0],
[0, b, 0],
[0, c * cos(alpha), c * sin(alpha)],
]
continue
elif alpha < 90 and b < c:
# print ('510-515')
transf = np.zeros(shape=(3, 3))
transf[0][t[0]] = 1
transf[1][t[1]] = 1
transf[2][t[2]] = 1
alpha = np.pi * alpha / 180
new_matrix = [
[a, 0, 0],
[0, b, 0],
[0, c * cos(alpha), c * sin(alpha)],
]
if new_matrix is None:
# print ('523-530')
# this if is to treat the case
# where alpha==90 (but we still have a monoclinic sg
new_matrix = [
[sorted_lengths[0], 0, 0],
[0, sorted_lengths[1], 0],
[0, 0, sorted_lengths[2]],
]
transf = np.zeros(shape=(3, 3))
for c in range(len(sorted_dic)):
transf[c][sorted_dic[c]["orig_index"]] = 1
if international_monoclinic:
# The above code makes alpha the non-right angle.
# The following will convert to proper international convention
# that beta is the non-right angle.
op = [[0, 1, 0], [1, 0, 0], [0, 0, -1]]
transf = np.dot(op, transf)
new_matrix = np.dot(op, new_matrix)
beta = Lattice(new_matrix).beta
if beta < 90:
op = [[-1, 0, 0], [0, -1, 0], [0, 0, 1]]
transf = np.dot(op, transf)
new_matrix = np.dot(op, new_matrix)
latt = Lattice(new_matrix)
elif latt_type == "triclinic":
# we use a LLL Minkowski-like reduction for the triclinic cells
struct = struct.get_lll_reduced_structure()
a, b, c = latt.abc # lengths
alpha, beta, gamma = [np.pi * i / 180 for i in latt.angles]
new_matrix = None
test_matrix = [
[a, 0, 0],
[b * cos(gamma), b * sin(gamma), 0.0],
[
c * cos(beta),
c * (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma),
c
* np.sqrt(
sin(gamma) ** 2
- cos(alpha) ** 2
- cos(beta) ** 2
+ 2 * cos(alpha) * cos(beta) * cos(gamma)
)
/ sin(gamma),
],
]
def is_all_acute_or_obtuse(m):
recp_angles = np.array(Lattice(m).reciprocal_lattice().angles)
return np.all(recp_angles <= 90) or np.all(recp_angles > 90)
if is_all_acute_or_obtuse(test_matrix):
transf = np.eye(3)
new_matrix = test_matrix
test_matrix = [
[-a, 0, 0],
[b * cos(gamma), b * sin(gamma), 0.0],
[
-c * cos(beta),
-c * (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma),
-c
* np.sqrt(
sin(gamma) ** 2
- cos(alpha) ** 2
- cos(beta) ** 2
+ 2 * cos(alpha) * cos(beta) * cos(gamma)
)
/ sin(gamma),
],
]
if is_all_acute_or_obtuse(test_matrix):
transf = [[-1, 0, 0], [0, 1, 0], [0, 0, -1]]
new_matrix = test_matrix
test_matrix = [
[-a, 0, 0],
[-b * cos(gamma), -b * sin(gamma), 0.0],
[
c * cos(beta),
c * (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma),
c
* np.sqrt(
sin(gamma) ** 2
- cos(alpha) ** 2
- cos(beta) ** 2
+ 2 * cos(alpha) * cos(beta) * cos(gamma)
)
/ sin(gamma),
],
]
if is_all_acute_or_obtuse(test_matrix):
transf = [[-1, 0, 0], [0, -1, 0], [0, 0, 1]]
new_matrix = test_matrix
test_matrix = [
[a, 0, 0],
[-b * cos(gamma), -b * sin(gamma), 0.0],
[
-c * cos(beta),
-c * (cos(alpha) - cos(beta) * cos(gamma)) / sin(gamma),
-c
* np.sqrt(
sin(gamma) ** 2
- cos(alpha) ** 2
- cos(beta) ** 2
+ 2 * cos(alpha) * cos(beta) * cos(gamma)
)
/ sin(gamma),
],
]
if is_all_acute_or_obtuse(test_matrix):
transf = [[1, 0, 0], [0, -1, 0], [0, 0, -1]]
new_matrix = test_matrix
latt = Lattice(new_matrix)
new_coords = np.dot(transf, np.transpose(struct.frac_coords)).T
new_struct = Atoms(
lattice_mat=latt.matrix,
elements=struct.elements,
coords=new_coords,
cartesian=False,
)
return new_struct
[docs]def parse_xyz_string(xyz_string):
"""
Convert xyz info to translation and rotation vectors.
Adapted from pymatgen.
Args:
xyz_string: string of the form 'x, y, z', '-x, -y, z',
'-2y+1/2, 3x+1/2, z-y+1/2', etc.
Returns:
translation and rotation vectors.
"""
from jarvis.core.utils import parse_xyz_string
return parse_xyz_string(xyz_string)
[docs]def operate_affine(cart_coord=[], affine_matrix=[]):
"""Operate affine method."""
affine_point = np.array([cart_coord[0], cart_coord[1], cart_coord[2], 1])
return np.dot(np.array(affine_matrix), affine_point)[0:3]
[docs]def get_new_coord_for_xyz_sym(frac_coord=[], xyz_string=""):
"""Obtain new coord from xyz string."""
from jarvis.core.utils import get_new_coord_for_xyz_sym
return get_new_coord_for_xyz_sym(
frac_coord=frac_coord, xyz_string=xyz_string
)
[docs]def check_duplicate_coords(coords=[], coord=[]):
"""Check if a coordinate exists."""
from jarvis.core.utils import check_duplicate_coords
return check_duplicate_coords(coords=coords, coord=coord)
"""
if __name__ == "__main__":
x = get_wyckoff_position_operators(488)
print (x)
import sys
#sys.exit()
box = [[2.715, 2.715, 0], [0, 2.715, 2.715], [2.715, 0, 2.715]]
coords = [[0, 0, 0], [0.25, 0.25, 0.25]]
elements = ["Si", "Si"]
Si = Atoms(lattice_mat=box, coords=coords, elements=elements)
spg = Spacegroup3D(atoms=Si) # .spacegroup_data()
#print(spg.space_group_symbol)
#print(spg.space_group_number)
#primt = spg.primitive_atoms
#print("primt", primt)
#print("cryst_sys", spg.crystal_system)
#print("point group", spg.point_group_symbol)
#print("cvn", spg.conventional_standard_structure)
cvn = spg.conventional_standard_structure
ml=symmetrically_distinct_miller_indices(max_index=3, cvn_atoms=cvn)
print ('miller=',ml,len(ml))
# ml_pmg=get_symmetrically_distinct_miller_indices(cvn.pymatgen_converter(),max_index=3)
# print ('millerpmg=',ml_pmg,len(ml_pmg))
# pmg=Si.pymatgen_converter()
# from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
# spg=SpacegroupAnalyzer(pmg)
# print (spg.get_space_group_symbol(),spg.get_space_group_number())
# print (pmg.get_primitive_structure())
"""