Crystal

Crystal:

The crystal module defines the crystal class, and GroupOp for group operations.

Crystal class

Class to store definition of a crystal, along with some analysis 1. geometric analysis (nearest neighbor displacements) 2. space group operations 3. point group operations for each basis position 4. Wyckoff position generation (for interstitials)

crystal.CombineTensorBasis(b1, b2, symmetric=True)

Combines (intersects) two tensor spaces into one; uses SVD to compute null space.

Parameters

• b1 – list of tensors

• b2 – list of tensors

Return tensorbasis

list of 2nd-rank symmetric tensors making up the basis

crystal.CombineVectorBasis(b1, b2)

Combines (intersects) two vector spaces into one.

Parameters

• b1 – (dim, vect) – dimensionality (0..3), vector defining line direction (1) or plane normal (2)

• b2 – (dim, vect)

Return dim

dimensionality, 0..3

Return vect

vector defining line direction (1) or plane normal (2)

class crystal.Crystal(lattice, basis, chemistry=None, spins=None, NOSYM=False, noreduce=False, threshold=1e-08)

A class that defines a crystal, as well as the symmetry analysis that goes along with it. Now includes optional spins. These can be vectors or “scalar” spins, for which we need to consider a phase factor. In general, they can be complex. Ideally, they should have magnitude either 0 or 1.

Specified by a lattice (3 vectors), a basis (list of lists of positions in direct coordinates). Can also name the elements (chemistry), and specify spin degrees of freedom.

classmethod BCC(a0, chemistry=None)

Create a body-centered cubic crystal with lattice constant a0

Parameters

a0 – lattice constant

Return BCC crystal

classmethod FCC(a0, chemistry=None)

Create a face-centered cubic crystal with lattice constant a0

Parameters

a0 – lattice constant

Return FCC crystal

FullVectorBasis(chem=None)

Generate our full vector basis, using the information from our crystal

Parameters

chem – (optional) chemical index to consider; otherwise return a list of such

Return VBfunctions

(list) of our unique vector basis lattice functions, normalized; each is an array (NVbasis x Nsites x 3)

Return VVouter

(list) of ouf VV “outer” expansion (NVbasis x NVbasis for each chemistry)

classmethod HCP(a0, c_a=1.632993161855452, chemistry=None)

Create a hexagonal closed packed crystal with lattice constant a0, c/a ratio c_a

Parameters

• a0 – lattice constant

• c_a – (optional) c/a ratio, default=ideal \(\sqrt{8/3}\)

Return HCP crystal

SymmTensorBasis(ind)

Generates the symmetric tensor basis corresponding to an atomic site

Parameters

ind – tuple index for atom

Return tensorbasis

list of 2nd-rank symmetric tensors making up the basis

VectorBasis(ind)

Generates the vector basis corresponding to an atomic site

Parameters

ind – tuple index for atom

Return dim

dimensionality, 0..3

Return vect

vector defining line direction (1) or plane normal (2)

Wyckoffpos(uvec)

Generates all the equivalent Wyckoff positions for a unit cell vector.

Parameters

uvec – 3-vector (float) vector in direct coordinates

Return Wyckofflist

list of equivalent Wyckoff positions

__init__(lattice, basis, chemistry=None, spins=None, NOSYM=False, noreduce=False, threshold=1e-08)

Initialization; starts off with the lattice vector definition and the basis vectors. While it does not explicitly store the specific chemical elements involved, it does store that there are different elements.

Parameters

• lattice – array[3,3] or list of array[3] (or 2 if 2-dimensional) lattice vectors; if [3,3] array, then the vectors need to be in column format so that the first lattice vector is lattice[:,0]

• basis – list of array[3] or list of list of array[3] (or 2 if 2-dimensional) crystalline basis vectors, in unit cell coordinates. If a list of lists, then there are multiple chemical elements, with each list corresponding to a unique element

• chemistry – (optional) list of names of chemical elements

• spins – (optional) list of numbers (complex) / vectors or list of list of same spins for individual atoms; if not None, needs to match the basis. Can either be scalars or vectors, corresponding to collinear or non-collinear magnetism

• NOSYM – turn off all symmetry finding (except identity)

• noreduce – do not attempt to reduce the atomic basis

• threshold – threshold for symmetry equivalence (stored in the class)

__repr__()

String representation of crystal (lattice + basis)

__str__()

Human-readable version of crystal (lattice + basis)

__weakref__

list of weak references to the object (if defined)

addbasis(basis, chemistry=None, spins=None)

Returns a new Crystal object that contains additional sites (assumed to be new chemistry). This is intended to “add in” interstitial sites. Note: if the symmetry is to be maintained, should be the output from Wyckoffpos().

Parameters

• basis – list (or list of lists) of new sites

• chemistry – (optional) list of chemistry names

• spins – (optional) list of spins

Return Crystal

new crystal object, with additional sites

calcmetric()

Computes the volume of the cell and the metric tensor

Return volume

cell volume

Return metric tensor

3x3

cart2pos(v)

Return the lattvec and index corresponding to an atomic position in cartesian coord.

Parameters

v – 3-vector (float) position in Cartesian coordinates

Return lattvec

3-vector (integer) lattice vector in direct coordinates,

Return (c,i)

tuple of matching basis atom; None if no match

cart2unit(v)

Return the lattvec and unit cell coord. corresponding to a position in cartesian coord.

Parameters

v – 3-vector (float) position in Cartesian coordinates

Return lattvec

3-vector (integer) lattice vector in direct coordinates,

Return uvec

3-vector (float) inside unit cell, in direct coordinates

center()

Center the atoms in the cell if there is an inversion operation present.

chemindex(chemistry)

Return index corresponding to chemistry; None if not present.

Parameters

chemistry – value to check

Return index

corresponding to chemistry

classmethod fromdict(yamldict)

Creates a Crystal object from a very simple YAML-created dictionary

Parameters

yamldict – dictionary; must contain ‘lattice’ (using row vectors!) and ‘basis’; can contain optional ‘lattice_constant’

Return Crystal(lattice.T, basis)

new crystal object

fullkptmesh(Nmesh)

Creates a k-point mesh of density given by Nmesh; does not symmetrize but does put the k-points inside the BZ. Does not return any weights as every point is equally weighted.

Parameters

Nmesh – mesh divisions Nmesh[0] x Nmesh[1] x Nmesh[2]

Return kpt

array[Nkpt][3] of kpoints

g_cart(g, x)

Apply a space group operation to a (Cartesian) vector position

Parameters

• g – group operation (GroupOp)

• x – 3-vector position in space

Return gx

3-vector position in space (Cartesian coordinates)

static g_direc(g, direc)

Apply a space group operation to a direction

Parameters

• g – group operation (GroupOp)

• direc – 3-vector direction

Return gdirec

3-vector direction

g_direc_equivalent(d1, d2, threshold=1e-08)

Tells us if two directions are equivalent by according to the space group

Parameters

• d1 – direction one (array[3])

• d2 – direction two (array[3])

• threshold – threshold for equality

Return equivalent

True if equivalent by a point group operation

g_pos(g, lattvec, ind)

Apply a space group operation to an atom position specified by its lattice and index

Parameters

• g – group operation (GroupOp)

• lattvec – 3-vector (integer) lattice vector in direct coordinates

• ind – two-tuple index specifying the atom: (atomtype, atomindex)

Return glatt

3-vector (integer) lattice vector in direct coordinates

Return gindex

tuple of new basis atom

static g_tensor(g, tensor)

Apply a space group operation to a 2nd-rank tensor

Parameters

• g – group operation (GroupOp)

• tensor – 2nd-rank tensor

Return gtensor

2nd-rank tensor

static g_vect(g, lattvec, uvec)

Apply a space group operation to a vector position specified by its lattice and a location in the unit cell in direct coordinates

Parameters

• g – group operation (GroupOp)

• lattvec – 3-vector (integer) lattice vector in direct coordinates

• uvec – 3-vector (float) vector in direct coordinates

• guvec – 3-vector (float) vector in direct coordinates

Return glatt

3-vector (integer) lattice vector in direct coordinates

genBZG()

Generates the reciprocal lattice G points that define the Brillouin zone.

Return Garray

array of G vectors that define the BZ, in Cartesian coordinates

genWyckoffsets()

Generate our Wykcoff sets.

Return Wyckoffsets

set of sets of tuples of positions that correspond to identical Wyckoff positions

gengroup()

Generate all of the space group operations. Now handles spins! Doesn’t store spin phase factors for each group operation, though.

Return Gset

frozenset of group operations

genpoint()

Generate our point group indices. Done with crazy list comprehension due to the structure of our basis.

Return Gpointlists

list of lists of frozensets of point group operations that leave a site unchanged

inBZ(vec, BZG=None, threshold=1e-05)

Tells us if vec is inside our set of defining points.

Parameters

• vec – array [3], vector to be tested

• BGZ – array [:,3], optional (default = self.BZG), array of vectors that define the BZ

• threshold – double, optional, threshold to use for “equality”

Return inBZ

False if outside the BZ, True otherwise

jumpnetwork(chem, cutoff, closestdistance=0)

Generate the full jump network for a specific chemical index, out to a cutoff. Organized by symmetry-unique transitions. Note that i->j and j->i are always related to one-another, but by equivalence of transition state, not symmetry. Now updated with closest-distance parameter.

Parameters

• chem – index corresponding to the chemistry to consider

• cutoff – distance cutoff

• closestdistance – closest distance allowed in transition (can be a list)

Return jumpnetwork

list of symmetry-unique transitions; each is a list of tuples: ((i,j), dx) corresponding to jump from \(i \to j\) with vector \(\mathbf{\delta x}\)

jumpnetwork2lattice(chem, jumpnetwork)

Convert a “standard” jumpnetwork (that specifies displacement vectors dx) into a lattice representation, where we replace dx with the lattice vector from i to j.

Parameters

• chem – index corresponding to the chemistry to consider

• jumpnetwork – list of symmetry-unique transitions; each is a list of tuples: ((i,j), dx) corresponding to jump from \(i \to j\) with vector \(\mathbf{\delta x}\)

Return jumplattice

list of symmetry-unique transitions; each is a list of tuples: ((i,j), R) corresponding to jump from i in unit cell 0 -> j in unit cell R

minlattice()

Try to find the optimal lattice vector definition for a crystal. Our definition of optimal is (a) length of each lattice vector is minimal; (b) the vectors are ordered from shortest to longest; (c) the vectors have minimal dot product; (d) the basis is right-handed.

Works recursively, and in-place.

nnlist(ind, cutoff)

Generate the nearest neighbor list for a given cutoff. Only consider neighbor vectors for atoms of the same type. Returns a list of cartesian vectors.

Parameters

• ind – tuple index for atom

• cutoff – distance cutoff

Return nnlist

list of nearest neighbor vectors

pos2cart(lattvec, ind)

Return the cartesian coordinates of an atom specified by its lattice and index

Parameters

• lattvec – 3-vector (integer) lattice vector in direct coordinates

• ind – two-tuple index specifying the atom: (atomtype, atomindex)

Return v

3-vector (float) in Cartesian coordinates

reduce(threshold=None)

Reduces the lattice and basis, if needed. Works (tail) recursively.

Parameters

threshold – threshold for symmetry comparison; default = self.threshold

Algorithm is slightly complicated: we attempt to identify if there is a internal translation symmetry in the crystal (called t) that applies to all sites. Once identified, we transform the lattice vectors and basis into the “reduced” form of the cell. We use tail recursion to continue until no further reduction is possible. Will usually require some “polishing” on the unit cell after the fact.

We try to do this efficiently: we check the GCD of the site counts (called M); if it’s 1, we kick out. We check translations against the smallest site set first.

We try to do this carefully: We make sure that our translation can be expressed rationally with M as the denominator; this helps protect against roundoff error. When we reduce the atomic basis, we average the values that match. Finally, as we reduce, we also change the self.threshold value accordingly so that recursion uses the same “effective” threshold.

reducekptmesh(kptfull, threshold=None)

Takes a fully expanded mesh, and reduces it by symmetry. Assumes every point is equally weighted. We would need a different (more complicated) algorithm if not true…

Parameters

• kptfull – array[Nkpt][3] of kpoints

• threshold – threshold for symmetry equality

Return kptsymm

array[Nsymm][3] of kpoints

Return weight

array[Nsymm] of weights (integrates to 1)

remapbasis(supercell)

Takes the basis definition, and using a supercell definition, returns a new basis

Parameters

supercell – integer array[3,3]

Return atomic basis

list of list of positions

simpleYAML(a0=1.0)

Creates a simplified YAML dump, in case we don’t want to output the full symmetry analysis

Return YAML

string dump

sitelist(chem)

Return a list of lists of Wyckoff-related sites for a given chemistry. Done with a single list comprehension–useful as input for diffusion calculation

Parameters

chem – index corresponding to chemistry to consider

Return symmequivsites

list of lists of indices that are equivalent by symmetry

strain(eps)

Returns a new Crystal object that is a strained version of the current.

Parameters

eps – strain tensor

Return Crystal

new crystal object, strained

unit2cart(lattvec, uvec)

Return the cartesian coordinates of a position specified by its lattice and unit cell coordinates

Parameters

• lattvec – 3-vector (integer) lattice vector in direct coordinates

• uvec – 3-vector (float) unit cell vector in direct coordinates

Return v

3-vector (float) in Cartesian coordinates

static vectlist(vb)

Returns a list of orthonormal vectors corresponding to our vector basis.

Parameters

vb – (dim, v)

Return vlist

list of vectors

crystal.FourthRankIsotropic(F)

Returns the average and shear values from orientational averaging of a symmetric fourth-rank tensor.

Parameters

F[a,b,c,d] – symmetric fourth-rank tensor (F[abcd]=F[abcd]=F[bacd]=F[cdab])

Return average

average value = (F11+2F12)/3, orientationally averaged

Return shear

shear value = F44, orientationally averaged

class crystal.GroupOp

A class corresponding to a group operation. Based on namedtuple, so it is immutable.

Intended to be used in combination with Crystal, we have a few operations that can be defined out-of-the-box.

Parameters

• rot – np.array(3,3) integer idempotent matrix

• trans – np.array(3) real vector

• cartrot – np.array(3,3) real unitary matrix

• indexmap – tuples of tuples, containing the atom mapping

static GroupOp_constructor(loader, node)

Construct a GroupOp from YAML

static GroupOp_representer(dumper, data)

Output a GroupOp

__add__(other)

Add a translation to our group operation

__eq__(other)

Test for equality–we use numpy.isclose for comparison, since that’s what we usually care about

__hash__()

Hash, so that we can make sets of group operations

__mul__(other)

Multiply two group operations to produce a new group operation

__ne__(other)

Inequality == not __eq__

__sane__()

Return true if the cartrot and rot are consistent and ‘sane’

__str__()

Human-readable version of groupop

__sub__(other)

Add a (negative) translation to our group operation

eigen()

Returns the type of group operation (single integer) and eigenvectors. 1 = identity 2, 3, 4, 6 = n- fold rotation around an axis negative = rotation + mirror operation, perpendicular to axis “special cases”: -1 = mirror, -2 = inversion

eigenvect[0] = axis of rotation / mirror eigenvect[1], eigenvect[2] = orthonormal vectors to define the plane giving a right-handed coordinate system and where rotation around [0] is positive, and the positive imaginary eigenvector for the complex eigenvalue is [1] + i [2].

Return type

integer

Return eigenvectors

list of [ev0, ev1, ev2]

classmethod ident(basis)

Return a group operation corresponding to identity for a given basis

incell()

Return a version of groupop where the translation is in the unit cell

inhalf()

Return a version of groupop where the translation is in the centered unit cell

inv()

Construct and return the inverse of the group operation

static optype(rot)

Returns the type of group operation (single integer): 1 = identity 2, 3, 4, 6 = n- fold rotation around an axis negative = rotation + mirror operation, perpendicular to axis “special cases”: -1 = mirror, -2 = inversion

Parameters

rot – rotation matrix (can be the integer rot)

Return type

integer

crystal.ProjectTensorBasis(tensor, basis)

Given a tensor, project it onto the basis.

Parameters

• tensor – tensor

• basis – list consisting of an orthonormal basis

Return tensor

tensor, projected

crystal.SymmTensorBasis(rottype, eigenvect)

Returns a symmetric second-rank tensor basis corresponding to the optype and eigenvectors for a GroupOp

Parameters

• rottype – output from eigen()

• eigenvect – eigenvectors

Return tensorbasis

list of 2nd-rank symmetric tensors making up the basis

crystal.VectorBasis(rottype, eigenvect)

Returns a vector basis corresponding to the optype and eigenvectors for a GroupOp

Parameters

• rottype – output from eigen()

• eigenvect – eigenvectors

Return dim

dimensionality, 0..3

Return vect

vector defining line direction (1) or plane normal (2)

crystal.Voigtstrain(e1, e2, e3, e4, e5, e6)

Returns a symmetric strain tensor from the Voigt reduced strain values.

Parameters

• e1 – xx

• e2 – yy

• e3 – zz

• e4 – yz + zx

• e5 – zx + xz

• e6 – xy + yx

Return strain

symmetric strain tensor

crystal.gcdlist(lis)

Returns the GCD of a list of integers

crystal.incell(vec)

Returns the vector inside the unit cell (in [0,1)**3)

Parameters

vec – 3-vector (unit coord)

Returns

3-vector

crystal.inhalf(vec)

Returns the vector inside the centered cell (in [-0.5,0.5)**3)

Parameters

vec – 3-vector (unit coord)

Returns

3-vector

crystal.isotropicFourthRank(average, shear)

Returns a symmetrized, isotropic fourth-rank tensor based on an average value and “shear” value

Parameters

• average – averaged value = (F11+2F12)/3

• shear – shear value = F44 = (F11-F12)/2

Return F[a,b,c,d]

isotropic fourth-rank tensor

crystal.maptranslation(oldpos, newpos, oldspins=None, newspins=None, threshold=1e-08)

Given a list of transformed positions, identify if there’s a translation vector that maps from the current positions to the new position.

The mapping specifies the index that the translated atom corresponds to in the original position set. If unable to construct a mapping, the mapping return is None; the translation vector will be meaningless.

If old/newspins are given then ONLY mappings that maintain spin are considered. This means that a loop is needed to consider possible spin phase factors.

Parameters

• oldpos – list of list of array[3]

• newpos – list of list of array[3], same layout as oldpos

• oldspins – (optional) list of list of numbers/arrays

• newspins – (optional) list of list of numbers/arrays

Return translation

array[3]

Return mapping

list of list of indices

crystal.ndarray_representer(dumper, data)

Output a numpy array