CrystalStars

CrystalStars:

The crystalStars module defines the classes corresponding to stars (in this case, for solute-vacancy complexes that are equivalent by space group symmetry), and vector stars (the inclusion of a vector basis on the stars). These modules are primarily responsible for all the symmetry analysis, and converting that into matrix forms for rapid numerical evaluation as needed.

Stars module, modified to work with crystal class

Classes to generate star sets, double star sets, and vector star sets; a lot of indexing functionality.

NOTE: The naming follows that of stars; the functionality is extremely similar, and this code was modified as little as possible to translate that functionality to crystals which possess a basis. In the case of a single atom basis, this should reduce to the stars object functionality.

The big changes are:

• Replacing NNvect star (which represents the jumps) with the jumpnetwork type found in crystal

• Using the jumpnetwork_latt representation from crystal

• Representing a “point” as a solute + vacancy. In this case, it is a tuple (s,v) of unit cell indices and a vector dx or dR (dx = Cartesian vector pointing from solute to vacancy; dR = lattice vector pointing from unit cell of solute to unit cell of vacancy). This is equivalent to our old representation if the tuple (s,v) = (0,0) for all sites. Due to translational invariance, the solute always stays inside the unit cell

• Using indices into the point list rather than just making lists of the vectors themselves. This is because the “points” now have a more complex representation (see above).

crystalStars.PSlist2array(PSlist)

Take in a list of pair states; return arrays that can be stored in HDF5 format

Parameters

PSlist – list of pair states

Return ij

int_array[N][2] = (i,j)

Return R

int[N][3]

Return dx

float[N][3]

class crystalStars.PairState

A class corresponding to a “pair” state; in this case, a solute-vacancy pair, but can also be a transition state pair. The solute (or initial state) is in unit cell 0, in position indexed i; the vacancy (or final state) is in unit cell R, in position indexed j. The cartesian vector dx connects them. We can add and subtract, negate, and “endpoint” subtract (useful for determining what Green function entry to use)

Parameters

• i – index of the first member of the pair (solute)

• j – index of the second member of the pair (vacancy)

• R – lattice vector pointing from unit cell of i to unit cell of j

• dx – Cartesian vector pointing from first to second member of pair

static PairState_constructor(loader, node)

Construct a GroupOp from YAML

static PairState_representer(dumper, data)

Output a PairState

__add__(other)

Add two states: works if and only if self.j == other.i (i,j) R + (j,k) R’ = (i,k) R+R’ : works for thinking about transitions… Note: a + b != b + a, and may be that only one of those is even defined

__eq__(other)

Test for equality–we don’t bother checking dx

__hash__()

Hash, so that we can make sets of states

__ne__(other)

Inequality == not __eq__

__neg__()

Negation of state (swap members of pair) - (i,j) R = (j,i) -R Note: a + (-a) == (-a) + a == 0 because we define what “zero” is.

__sane__(crys, chem)

Determine if the dx value makes sense given everything else…

__str__()

Human readable version

__sub__(other)

Add a negative: a-b points from initial of a to initial of b if same final state (i,j) R - (k,j) R’ = (i,k) R-R’ Note: this means that (a-b) + b = a, but b + (a-b) is an error. (b-a) + a = b

__xor__(other)

Subtraction on the endpoints (sort of the “opposite” of a-b): a^b points from final of b to final of a if same initial state (i,j) R ^ (i,k) R’ = (k,j) R-R’ Note: b + (a^b) = a but (a^b) + b is an error. a + (b^a) = b

classmethod fromcrys(crys, chem, ij, dx)

Convert (i,j), dx into PairState

classmethod fromcrys_latt(crys, chem, ij, R)

Convert (i,j), R into PairState

g(crys, chem, g)

Apply group operation.

Parameters

• crys – crystal

• chem – chemical index

• g – group operation (from crys)

Return g*PairState

corresponding to group operation applied to self

iszero()

Quicker than self == PairState.zero()

classmethod zero(n=0, dim=3)

Return a “zero” state

class crystalStars.StarSet(jumpnetwork, crys, chem, Nshells=0, originstates=False, lattice=False)

A class to construct crystal stars, and be able to efficiently index.

Takes in a jumpnetwork, which is used to construct the corresponding stars, a crystal object with which to operate, a specification of the chemical index for the atom moving (needs to be consistent with jumpnetwork and crys), and then the number of shells.

In this case, shells = number of successive “jumps” from a state. As an example, in FCC, 1 shell = 1st neighbor, 2 shell = 1-4th neighbors.

__add__(other)

Add two StarSets together; done by making a copy of one, and iadding

__contains__(PS)

Return true if PS is in the star

__iadd__(other)

Add another StarSet to this one; very similar to generate()

__init__(jumpnetwork, crys, chem, Nshells=0, originstates=False, lattice=False)

Initiates a star set generator for a given jumpnetwork, crystal, and specified chemical index. Does not include “origin states” by default; these are PairStates that iszero() is True; they are only needed if crystal has a nonzero VectorBasis.

Parameters

• jumpnetwork – list of symmetry unique jumps, as a list of list of tuples; either ((i,j), dx) for jump from i to j with displacement dx, or ((i,j), R) for jump from i in unit cell 0 -> j in unit cell R

• crys – crystal where jumps take place

• chem – chemical index of atom to consider jumps

• Nshells – number of shells to generate

• originstates – include origin states in generate?

• lattice – which form does the jumpnetwork take?

__str__()

Human readable version

__weakref__

list of weak references to the object (if defined)

addhdf5(HDF5group)

Adds an HDF5 representation of object into an HDF5group (needs to already exist).

Example: if f is an open HDF5, then StarSet.addhdf5(f.create_group(‘StarSet’)) will (1) create the group named ‘StarSet’, and then (2) put the StarSet representation in that group.

Parameters

HDF5group – HDF5 group

copy(empty=False)

Return a copy of the StarSet; done as efficiently as possible; empty means skip the shells, etc.

diffgenerate(S1, S2, threshold=1e-08)

Construct a starSet using endpoint subtraction from starset S1 to starset S2. Will include zero. Points from vacancy states of S1 to vacancy states of S2.

Parameters

• S1 – starSet for start

• S2 – starSet for final

• threshold – threshold for sorting magnitudes (can influence symmetry efficiency)

generate(Nshells, threshold=1e-08, originstates=False)

Construct the points and the stars in the set. Does not include “origin states” by default; these are PairStates that iszero() is True; they are only needed if crystal has a nonzero VectorBasis.

Parameters

• Nshells – number of shells to generate; this is interpreted as subsequent “sums” of jumplist (as we need the solute to be connected to the vacancy by at least one jump)

• threshold – threshold for determining equality with symmetry

• originstates – include origin states in generate?

jumpnetwork_omega1()

Generate a jumpnetwork corresponding to vacancy jumping while the solute remains fixed.

Return jumpnetwork

list of symmetry unique jumps; list of list of tuples (i,f), dx where i,f index into states for the initial and final states, and dx = displacement of vacancy in Cartesian coordinates. Note: if (i,f), dx is present, so if (f,i), -dx

Return jumptype

list of indices corresponding to the (original) jump type for each symmetry unique jump; useful for constructing a LIMB approximation, and needed to construct delta_omega

Return starpair

list of tuples of the star indices of the i and f states for each symmetry unique jump

jumpnetwork_omega2()

Generate a jumpnetwork corresponding to vacancy exchanging with a solute.

Return jumpnetwork

list of symmetry unique jumps; list of list of tuples (i,f), dx where i,f index into states for the initial and final states, and dx = displacement of vacancy in Cartesian coordinates. Note: if (i,f), dx is present, so if (f,i), -dx

Return jumptype

list of indices corresponding to the (original) jump type for each symmetry unique jump; useful for constructing a LIMB approximation, and needed to construct delta_omega

Return starpair

list of tuples of the star indices of the i and f states for each symmetry unique jump

classmethod loadhdf5(crys, HDF5group)

Creates a new StarSet from an HDF5 group.

Parameters

• crys – crystal object–MUST BE PASSED IN as it is not stored with the StarSet

• HDFgroup – HDF5 group

Return StarSet

new StarSet object

starindex(PS)

Return the index for the star to which pair state PS belongs; None if not found

stateindex(PS)

Return the index of pair state PS; None if not found

symmatch(PS1, PS2)

True if there exists a group operation that makes PS1 == PS2.

symmequivjumplist(i, f, dx)

Returns a list of tuples of symmetry equivalent jumps

Parameters

• i – index of initial state

• f – index of final state

• dx – displacement vector

Return symmjumplist

list of tuples of ((gi, gf), gdx) for every group op

class crystalStars.VectorStarSet(starset=None)

A class to construct vector star sets, and be able to efficiently index.

All based on a StarSet

GFexpansion()

Construct the GF matrix expansion in terms of the star vectors, and indexed to GFstarset.

Return GFexpansion

array[Nsv, Nsv, NGFstars] the GF matrix[i, j] = sum(GFexpansion[i, j, k] * GF(starGF[k]))

Return GFstarset

starSet corresponding to the GF

__init__(starset=None)

Initiates a vector-star generator; work with a given star.

Parameters

starset – StarSet, from which we pull nearly all of the info that we need

__weakref__

list of weak references to the object (if defined)

addhdf5(HDF5group)

Adds an HDF5 representation of object into an HDF5group (needs to already exist).

Example: if f is an open HDF5, then StarSet.addhdf5(f.create_group(‘VectorStarSet’)) will

(1) create the group named ‘VectorStarSet’, and then (2) put the VectorStarSet representation in that group.

Parameters

HDF5group – HDF5 group

bareexpansions(jumpnetwork, jumptype)

Construct the bare diffusivity expansion in terms of the jumpnetwork. We return the reference (0) contribution so that the change can be determined; this is useful for the vacancy contributions. This saves us from having to deal with issues with our outer shell where we only have a fraction of the escapes, but as long as the kinetic shell is one more than the thermodynamics (so that the interaction energy is 0, hence no change in probability), this will work. The PS (pair stars) is useful for including the probability factor for the endpoint of the jump; we just call it the ‘probfactor’ below.

Note also: this currently assumes that the displacement vector does not change between omega0 and omega(1/2).

Parameters

• jumpnetwork – jumpnetwork of symmetry unique omega1-type jumps, corresponding to our starset. List of lists of (IS, FS), dx tuples, where IS and FS are indices corresponding to states in our starset.

• jumptype – specific omega0 jump type that the jump corresponds to

Return D0expansion

array[3,3, Njump_omega0] the D0[a,b,jt] = sum(D0expansion[a,b, jt] * sqrt(probfactor0[PS[jt][0]]*probfactor0[PS[jt][1]) * omega0[jt])

Return D1expansion

array[3,3, Njump_omega1] the D1[a,b,k] = sum(D1expansion[a,b, k] * sqrt(probfactor[PS[k][0]]*probfactor[PS[k][1]) * omega[k])

biasexpansions(jumpnetwork, jumptype, omega2=False)

Construct the bias1 and bias0 vector expansion in terms of the jumpnetwork. We return the bias0 contribution so that the db = bias1 - bias0 can be determined. This saves us from having to deal with issues with our outer shell where we only have a fraction of the escapes, but as long as the kinetic shell is one more than the thermodynamics (so that the interaction energy is 0, hence no change in probability), this will work. The PS (pair stars) is useful for including the probability factor for the endpoint of the jump; we just call it the ‘probfactor’ below. Note: this used to be separated into bias1expansion, and bias2expansion,and had terms that are now in rateexpansions. Note also that if jumpnetwork_omega2 is passed, it also works for that. However, in that case we have a different approach for the calculation of bias1expansion: if there are origin states, they get the negative summed bias of the others.

Parameters

• jumpnetwork – jumpnetwork of symmetry unique omega1-type jumps, corresponding to our starset. List of lists of (IS, FS), dx tuples, where IS and FS are indices corresponding to states in our starset.

• jumptype – specific omega0 jump type that the jump corresponds to

• omega2 – (optional) are we dealing with the omega2 list, so we need to remove origin states? (default=False)

Return bias0expansion

array[Nsv, Njump_omega0] the gen0 vector[i] = sum(bias0expasion[i, k] * sqrt(probfactor0[PS[k]]) * omega0[k])

Return bias1expansion

array[Nsv, Njump_omega1] the gen1 vector[i] = sum(bias1expansion[i, k] * sqrt(probfactor[PS[k]] * omega1[k])

generate(starset, threshold=1e-08)

Construct the actual vectors stars

Parameters

starset – StarSet, from which we pull nearly all of the info that we need

generateouter()

Generate our outer products for our star-vectors.

Return outer

array [3, 3, Nvstars, Nvstars] outer[:, :, i, j] is the 3x3 tensor outer product for two vector-stars vs[i] and vs[j]

classmethod loadhdf5(SSet, HDF5group)

Creates a new VectorStarSet from an HDF5 group.

Parameters

• SSet – StarSet–MUST BE PASSED IN as it is not stored with the VectorStarSet

• HDFgroup – HDF5 group

Return VectorStarSet

new VectorStarSet object

originstateVectorBasisfolddown(elemtype='solute')

Construct the expansion to “fold down” from vector stars to origin states.

Parameters

elemtype – ‘solute’ of ‘vacancy’, depending on which site we need to reduce

Return OSindices

list of indices corresponding to origin states

Return folddown

[NOS, Nvstars] to map vector stars to origin states

Return OS_VB

[NOS, Nsites, 3] mapping of origin state to a vector basis

rateexpansions(jumpnetwork, jumptype, omega2=False)

Construct the omega0 and omega1 matrix expansions in terms of the jumpnetwork; includes the escape terms separately. The escape terms are tricky because they have probability factors that differ from the transitions; the PS (pair stars) is useful for finding this. We just call it the ‘probfactor’ below. Note: this used to be separated into rate0expansion, and rate1expansion, and partly in bias1expansion. Note also that if jumpnetwork_omega2 is passed, it also works for that. However, in that case we have a different approach for the calculation of rate0expansion: if there are origin states, then we need to “jump” to those; if there is a non-empty VectorBasis we will want to account for them there.

Parameters

• jumpnetwork – jumpnetwork of symmetry unique omega1-type jumps, corresponding to our starset. List of lists of (IS, FS), dx tuples, where IS and FS are indices corresponding to states in our starset.

• jumptype – specific omega0 jump type that the jump corresponds to

• omega2 – (optional) are we dealing with the omega2 list, so we need to remove origin states? (default=False)

Return rate0expansion

array[Nsv, Nsv, Njump_omega0] the omega0 matrix[i, j] = sum(rate0expansion[i, j, k] * omega0[k]); IF NVB>0 we “hijack” this and use it for [NVB, Nsv, Njump_omega0], as we’re doing an omega2 calc and rate0expansion won’t be used anyway.

Return rate0escape

array[Nsv, Njump_omega0] the escape contributions: omega0[i,i] += sum(rate0escape[i,k]*omega0[k]*probfactor0(PS[k]))

Return rate1expansion

array[Nsv, Nsv, Njump_omega1] the omega1 matrix[i, j] = sum(rate1expansion[i, j, k] * omega1[k])

Return rate1escape

array[Nsv, Njump_omega1] the escape contributions: omega1[i,i] += sum(rate1escape[i,k]*omega1[k]*probfactor(PS[k]))

crystalStars.array2PSlist(ij, R, dx)

Take in arrays of ij, R, dx (from HDF5), return a list of PairStates

Parameters

• ij – int_array[N][2] = (i,j)

• R – int[N][3]

• dx – float[N][3]

Return PSlist

list of pair states

crystalStars.doublelist2flatlistindex(listlist)

Takes a list of lists, returns a flattened list and an index array

Parameters

listlist – list of lists of objects

Return flatlist

flat list of objects (preserving order)

Return indexarray

array indexing which original list it came from

crystalStars.flatlistindex2doublelist(flatlist, indexarray)

Takes a flattened list and an index array, returns a list of lists

Parameters

• flatlist – flat list of objects (preserving order)

• indexarray – array indexing which original list it came from

Return listlist

list of lists of objects

crystalStars.zeroclean(x, threshold=1e-08)

Modify x in place, return 0 if x is below a threshold; useful for “symmetrizing” our expansions