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

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, noreduce=True)

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’

• noreduce – should we pass on lattice and basis as is? (default=True)

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(rot, trans, cartrot, indexmap)

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