In [67]:
import onsager.OnsagerCalc as onsager
import numpy as np
from scipy.constants import physical_constants
kB = physical_constants['Boltzmann constant in eV/K'][0]
In [2]:
def pos2site(uvec, uindex, N, Nsites):
"""
Takes site (indexed by uindex) and supercell vector uvec and returns site
index, assuming PBC
"""
return uindex + Nsites*(uvec[0]%N[0] + N[0]*(uvec[1]%N[1] + N[1]*(uvec[2]%N[2])))
In [3]:
def map2sites(N, Nsites, sitelist, f):
"""
Takes an array f defined for unique sites (given by sitelist) and maps it onto
an array corresponding to the number of sites in the supercell given by N.
:param N: int[3] specifying the supercell
:param Nsites: number of sites in the unit cell
:param sitelist: list of lists of sites; grouped by symmetry equivalence
:param f: f, indexed by symmetry equivalence
:return fmap: f, mapped out to all the site indices in the supercell
"""
Nsuper = Nsites*N[0]*N[1]*N[2]
fmap = np.zeros(Nsuper)
invmap = []
for i in range(Nsites):
for lis,fval in zip(sitelist, f):
if i in lis: break
for nsuper in [np.array([n0, n1, n2])
for n2 in range(N[2])
for n1 in range(N[1])
for n0 in range(N[0])]:
fmap[pos2site(nsuper, i, N, Nsites)] = fval
return fmap
In [4]:
def breakdownjumpnetwork(crys, chem, jumpnetwork):
"""
Takes a crystal and jumpnetwork for a particular interstitial (specified by chem)
and returns a "simplified" version of the jumpnetwork, suitable for a supercell.
:param crys: Crystal object
:param chem: integer index corresponding to chemistry of jump atom
:param jumpnetwork: list of lists of ((i,j), dx) tuples: i->j and displacement
:return trans: list, indexed by site; in each, indexed by transition; in that
a tuple of (endsite, duvec, dx, Etransindex) -- endsite index,
delta uvec, displacement and index of transition energy
:return Etrans: list--for each transition (corresponding to Etransindex), which
unique transition from jumpnetwork it was
"""
trans = [[] for n in crys.basis[chem]]
Etrans = []
for Etransunique, jumplist in enumerate(jumpnetwork):
for ((i,j), dx) in jumplist:
# analyze! determine the uvec
uvec = (np.round(np.dot(crys.invlatt, dx) -
crys.basis[chem][j] + crys.basis[chem][i])).astype(int)
# either we've seen this jump before or we haven't
Etransindex = -1
for (i0, u0, dx0, E0) in trans[j]:
if (i==i0) and np.all(np.equal(uvec, -u0)):
Etransindex = E0
break
if Etransindex < 0:
Etransindex = len(Etrans)
Etrans.append(Etransunique)
trans[i].append((j, uvec, dx, Etransindex))
return trans, Etrans
In [5]:
def makesupercellKMC(N, trans, Etrans):
"""
Takes an analyzed set of jumps (in terms of a trans and Etrans list) and for
an N[0] x N[1] x N[2] supercell with Nsites per unit cell, constructs a
nearly complete KMC transition matrix.
:param N: vector of supercell size
:param trans: list, indexed by site; in each, indexed by transition; in that
a tuple of (endsite, duvec, dx, Etransindex) -- endsite index,
delta uvec, displacement and index of transition energy
:param Etrans: list--for each transition (corresponding to Etransindex), which
unique transition from jumpnetwork it was
:return transsuper: list, indexed by supercell site; in each, indexed by transition;
in that a tuple of (endsite, dx, Etransindex)
:return Etranssuper: list--for each unique transition in the supercell of the
corresponding index in jumpnetwork
"""
Nsites = len(trans)
Ntrans = len(Etrans)
transsuper = [[] for n in range(N[0]*N[1]*N[2]*Nsites)]
Etranssuper = [0 for n in range(N[0]*N[1]*N[2]*Ntrans)]
for i, translist in enumerate(trans):
for (j, du, dx, Etransindex) in translist:
for nsuper in [np.array([n0, n1, n2])
for n2 in range(N[2])
for n1 in range(N[1])
for n0 in range(N[0])]:
ni = pos2site(nsuper, i, N, Nsites)
nj = pos2site(nsuper + du, j, N, Nsites)
Eindex = min(pos2site(nsuper, Etransindex, N, Ntrans),
pos2site(nsuper + du, Etransindex, N, Ntrans))
transsuper[ni].append((nj, dx, Eindex))
Etranssuper[Eindex] = Etrans[Etransindex]
return transsuper, Etranssuper
In [6]:
def makeKMCmatrices(transsuper, pre, betaene, preT, betaeneT, computebias=False):
"""
Takes in the list, indexed by supercell site, of transitions and an indexing
of transition states, along with site and transition thermodynamics
(prefactors and energies) and return the KMC matrices
:param transsuper: list, indexed by supercell site; in each, indexed by transition;
in that a tuple of (endsite, dx, Etransindex)
:param pre: site prefactors
:param betaene: site energies / kB T
:param preT: transition state prefactors
:param betaeneT: transition state energy / kB T
:param computebias: whether to return the "bias" vector--for testing
:return transarray: array of transitions indices [site][t]
:return transcumprob: array of transition cumulative probabilities [site][t]
:return transdelta: array of displacement vector for transitions [site][t][3]
:return escapetime: array of escape times [site]
:return sitecumprob: array of cumulative probability for sites [site]
:return biasvect: rate-bias vector (or velocity vector); *should be zero*
"""
Nsite = len(transsuper)
Ntrans = max(len(t) for t in transsuper) # maximum number of transitions
transarray = np.zeros((Nsite,Ntrans), dtype=int)
transcumprob = np.zeros((Nsite,Ntrans))
transdelta = np.zeros((Nsite,Ntrans,3))
escapetime = np.zeros(Nsite)
biasvect = np.zeros(3)
# first up, site probabilities:
minbetaene = min(betaene) # avoiding underflow/overflow just in case
siteprob = np.array([p*np.exp(minbetaene-bE) for p,bE in zip(pre,betaene)])
sitecumprob = np.cumsum(siteprob)
siteprob *= 1./sitecumprob[-1] # normalize site probabilities
sitecumprob *= 1./sitecumprob[-1]
# now, all the transition information
for i,translist in enumerate(transsuper):
bE = betaene[i]
p = pre[i]
for t, (j, dx, Eindex) in enumerate(translist):
transarray[i][t] = j
transdelta[i][t] = dx
transcumprob[i][t] = preT[Eindex]*np.exp(bE-betaeneT[Eindex])/p
biasvect += siteprob[i]*transcumprob[i][t]*dx
for t in range(len(translist), Ntrans):
transarray[i][t] = j # safety -- should not be accessed, but...
transdelta[i][t] = dx # safety -- should not be accessed, but...
transcumprob[i] = np.cumsum(transcumprob[i])
escapetime[i] = 1./transcumprob[i][-1]
transcumprob[i] *= escapetime[i]
if not computebias:
return transarray, transcumprob, transdelta, escapetime, sitecumprob
else:
return transarray, transcumprob, transdelta, escapetime, sitecumprob, biasvect
In [74]:
def runKMC(transarray, transcumprob, transdelta, escapetime, sitecumprob,
Nstep=4, Nrun=256, seed=None):
"""
Take all the output from makeKMCmatrices, and actually run a KMC simulation.
:param transarray: array of transitions indices [site][t]
:param transcumprob: array of transition cumulative probabilities [site][t]
:param transdelta: array of displacement vector for transitions [site][t][3]
:param escapetime: array of escape times [site]
:param sitecumprob: array of cumulative probability for sites [site]
:param Nstep: number of steps to run in a given KMC trajectory (1 KMC step = # sites)
:param Nrun: number of independent runs (needed to get good statistics)
:param seed: seed for RNG; if None, don't reseed
:return D: stochastic estimate of diffusivity from runs
:return dD: stochastic estimate of error on diffusivity from runs
"""
if seed is not None: np.random.seed(seed)
Nsite = transcumprob.shape[0]
Ntrans = transcumprob.shape[1]
D = np.zeros((3,3))
dD = np.zeros((3,3))
for nrun in range(Nrun):
dr = np.zeros(3) # displacement
T = 0 # time
# select an initial state
u = np.random.random()
for i, P in enumerate(sitecumprob):
if u < P: break
# get our random numbers
for u in np.random.random(Nsite*Nstep):
for t, P in enumerate(transcumprob[i]):
if u < P: break
dr += transdelta[i][t]
T += escapetime[i]
i = transarray[i][t]
D0 = np.outer(dr,dr)*(0.5/T) # <RR>/2t
D += D0
dD += D0*D0
invN = 1./Nrun
D *= invN
dD = np.sqrt(dD)*invN
return D, dD
In [8]:
def findneigh(crys, solute, chem, cutoff):
"""
Construct a list of neighboring sites of a specific chemistry to a
"solute" (identified by an index) in a crys within a cutoff. The returned
set is a list of lists of tuples, grouped by symmetry equivalence.
:param crys: Crystal, specifying the structure
:param solute: double index of the atom to map neighbors
:param chem: chemistry index to search
:param cutoff: maximum distance for neighbor search
:return neighlist: list of lists of tuples (index, np.array([u1, u2, u3]))
where inde is the atom index, and u1,u2,u3 is the unit cell
"""
r2 = cutoff*cutoff
nmax = [int(np.round(np.sqrt(crys.metric[i,i])))+1
for i in range(3)]
supervect = [ np.array([n0, n1, n2])
for n2 in range(-nmax[2],nmax[2]+1)
for n1 in range(-nmax[1],nmax[1]+1)
for n0 in range(-nmax[0],nmax[0]+1) ]
neighlist = []
u0 = crys.basis[solute[0]][solute[1]]
PG = crys.pointG[solute[0]][solute[1]]
for wyckset in crys.Wyckoff:
for (c,i) in wyckset:
if c == chem:
u1 = crys.basis[chem][i]
du = u1-u0
for n in supervect:
dx = crys.unit2cart(n, du)
if np.dot(dx, dx) > 0 and np.dot(dx,dx) < r2:
# see if there's already a set to which it belongs...
found = False
for lis in neighlist:
lrep = lis[0]
indpair = (chem,lrep[0])
if indpair in wyckset:
for g in PG:
unew, (c,j) = crys.g_pos(g, lrep[1], indpair)
if j == i and np.all(unew == n):
# belongs to this symmetry class
lis.append((i, n))
found = True
break
if found: break
if not found:
# new symmetry class
neighlist.append([(i, n)])
return neighlist
In [9]:
def insertsolute(neighlist, DEneigh, N, Nsites, Ntranssuper, transsuper, xDE = None):
"""
Takes in a neighborlist for a single solute, the interaction energies for each
type of neighbor, and the information about the supercell (N, Nsites, supercell
transition information), and constructs the change in site energy and transition
state energies as vectors that can be added to those vectors, using the
linear interpolation of migration barrier (LIMB) approximation. Allows for the
optional parameter x that specifies exactly how much DE gets added to the
transition state from the endpoints (default = 1/2)
:param neighlist: list of lists of tuples (ind, [u1,u2,u3]), grouped by symmetry
equivalent neighbors
:param DEneigh: list of length(neighlist), with the interaction energies
:param N: N[0] x N[1] x N[2] supercell
:param Nsites: number of sites in the unit cell
:param Ntranssuper: number of transition state energies in supercell
:param transsuper: list, indexed by supercell site; in each, indexed by transition;
in that a tuple of (endsite, dx, Etransindex)
:param xDE: [Nsites,Nsites] array; for site type i transitioning to j,
DE[i]*xDE[i][j] gets added to the transition state between i and j
Should obey xDE[i][j] + xDE[j][i] = 1; default = 1/2
:return Denesuper: vector of changes in site energies
:return DeneTsuper: vector of changes in transition energies
"""
Denesuper = np.zeros(len(transsuper))
DeneTsuper = np.zeros(Ntranssuper)
if xDE is None:
xDE = 0.5*np.ones((Nsites, Nsites))
for lis, DE in zip(neighlist, DEneigh):
for (i, u) in lis:
isite = pos2site(u, i, N, Nsites)
Denesuper[isite] += DE
for j, dx, Etransindex in transsuper[isite]:
DeneTsuper[Etransindex] += DE*xDE[i][j%Nsites]
return Denesuper, DeneTsuper
In [10]:
with open("bin/HCP-interstitial.yaml", "r") as in_f:
dict_def = onsager.crystal.yaml.load(in_f)
crys = onsager.crystal.Crystal.fromdict(dict_def) # pull out the crystal part of the YAML
# sites:
sitelist = dict_def['sitelist']
pre = dict_def['Prefactor']
ene = dict_def['Energy']
dipole = dict_def['Dipole']
# jumps
jumpnetwork = dict_def['jumpnetwork']
preT = dict_def['PrefactorT']
eneT = dict_def['EnergyT']
dipoleT = dict_def['DipoleT']
# we don't do any checking here... just dive on in
chem = dict_def['interstitial']
# create our calculator
interstitial = onsager.Interstitial(crys, chem, sitelist, jumpnetwork)
In [11]:
preT
Out[11]:
In [12]:
map2sites([2,2,2], 6, sitelist, ene)
Out[12]:
In [13]:
# this is a little python magic to directly pass breakdownjumpnetwork output...
makesupercellKMC([1,1,1], *(breakdownjumpnetwork(crys, chem, jumpnetwork)))
Out[13]:
In [14]:
N=[2,2,2]
Nsites=6
transsuper,Etranssuper = makesupercellKMC(N, *(breakdownjumpnetwork(crys, chem, jumpnetwork)))
presuper = map2sites(N, Nsites, sitelist, pre)
enesuper = map2sites(N, Nsites, sitelist, ene)
preTsuper = np.array([preT[i] for i in Etranssuper])
eneTsuper = np.array([eneT[i] for i in Etranssuper])
runKMC(*(makeKMCmatrices(transsuper, presuper, enesuper, preTsuper, eneTsuper)),
Nstep=1, Nrun=128*128)
Out[14]:
In [15]:
interstitial.diffusivity(pre, ene, preT, eneT)
Out[15]:
In [16]:
(Out[14][0]-Out[15])/Out[14][1]
Out[16]:
In [17]:
import math
1-math.erfc(1.1154)
Out[17]:
In [18]:
hcpTi = onsager.crystal.Crystal.HCP(2.933, 4.638/2.933)
print(hcpTi)
In [19]:
# add octahedral, hexahedral, and crowdion sites
hcpTiO = hcpTi.addbasis(hcpTi.Wyckoffpos(np.array([0.,0.,0.])) + \
hcpTi.Wyckoffpos(np.array([2./3.,1./3.,0.25])) + \
hcpTi.Wyckoffpos(np.array([0.5,0.,0.])))
print(hcpTiO)
In [20]:
hcpTiO.sitelist(1)
Out[20]:
In [21]:
hcpTiOrawjumps = hcpTiO.jumpnetwork(1, 2.9, 0.25)
len(hcpTiOrawjumps)
Out[21]:
In [22]:
for j, jump in enumerate(hcpTiOrawjumps):
print(j, np.sqrt(np.dot(jump[0][1], jump[0][1])), jump[0])
In [23]:
hcpTiOjumps = [hcpTiOrawjumps[n] for n in [0, 1, 3, 7]] # o-o o-h o-c h-c
for j, jump in enumerate(hcpTiOjumps):
print(j, np.sqrt(np.dot(jump[0][1], jump[0][1])), jump[0])
In [24]:
hcpTiOsitelist = hcpTiO.sitelist(1)
print(hcpTiOsitelist)
In [25]:
hcpTiOene = [1.88, 0., 1.19] # energies in eV
hcpTiOpre = [16.84/12.21, 1., 10.33/5.58] # prefactors for sites (o->* / *->o)
hcpTiOeneT = [3.25, 2.04, 2.16, 1.19 + 0.94] # transition energies in eV
hcpTiOpreT = [11.76e12, 10.33e12, 16.84e12, 10.27e12 * (10.33/5.58)] # prefactors in Hz
In [26]:
hcpTiOsoluteneigh = findneigh(hcpTiO, (0,0), 1, 2.8)
for lis in hcpTiOsoluteneigh:
ind, u = lis[0]
print(ind, len(lis), hcpTiO.pos2cart(u, (1,ind)) - hcpTiO.pos2cart(np.zeros(3,dtype=int), (0,0)))
In [27]:
# setup the LIMB approximation so that the transition state is closer to
# crowdion (3/4) than octahedral or hexahedral (1/4). All others remain halfway
xDE = 0.5*np.ones((10, 10))
for i in range(0,4):
for j in range(4,10):
xDE[i][j], xDE[j][i] = 0.25, 0.75
print(xDE)
In [28]:
N=[4,4,3]
Nsites=10
transsuper,Etranssuper = makesupercellKMC(N, *(breakdownjumpnetwork(hcpTiO, 1, hcpTiOjumps)))
presuper = map2sites(N, Nsites, hcpTiOsitelist, hcpTiOpre)
enesuper = map2sites(N, Nsites, hcpTiOsitelist, hcpTiOene)
preTsuper = np.array([hcpTiOpreT[i] for i in Etranssuper])
eneTsuper = np.array([hcpTiOeneT[i] for i in Etranssuper])
# crowd-far, crowd-near, oct, hex-c, hex-basal: Mn
Denesuper, DeneTsuper = insertsolute(hcpTiOsoluteneigh,
[0.38, -0.76, 0.44, 0.04, -0.22],
N, Nsites, len(Etranssuper), transsuper)
beta = 1.0/0.075 # (900K)^-1
D0,dD0 = runKMC(*(makeKMCmatrices(transsuper,
presuper,
beta*enesuper,
preTsuper,
beta*eneTsuper)),
Nstep=1, Nrun=128*128)
Dc,dDc = runKMC(*(makeKMCmatrices(transsuper,
presuper,
beta*(enesuper+Denesuper),
preTsuper,
beta*(eneTsuper+DeneTsuper))),
Nstep=4, Nrun=128*128)
print(D0)
print(Dc)
In [29]:
print(dD0)
print(dDc)
In [34]:
(Dc-D0)
Out[34]:
In [30]:
tl, Et = breakdownjumpnetwork(hcpTiO, 1, hcpTiOjumps)
[ [tr[3] for tr in lis] for lis in tl]
Out[30]:
In [36]:
print(Et)
In [76]:
KMCmatrices = makeKMCmatrices(transsuper, presuper, beta*(enesuper+Denesuper),
preTsuper, beta*(eneTsuper+DeneTsuper))
%timeit runKMC(*KMCmatrices, Nstep=4, Nrun=16*16)
In [80]:
# check that the total bias vector for the supercell == 0 (detailed balance)
makeKMCmatrices(transsuper, presuper, beta*(enesuper+Denesuper),
preTsuper, beta*(eneTsuper+DeneTsuper), True)[5]
Out[80]:
In [82]:
runKMC(*KMCmatrices, Nstep=4, Nrun=16*16)
Out[82]:
In [ ]: