Comparing to correlation coefficients from William D. Carlson and Clark R. Wilson, Phys Chem Minerals 43, 363-369 (2016) doi:10.1007/s00269-016-0800-2
Garnet structure includes pyrope, which we use as our example structure, with space group 230 (Ia3d) with stoichiometry Mg3Al2Si3O12. The occupied Wyckoff positions for this are (lattice constant $a_0$=1.1459 nm):
| Wykcoff site | chemistry | position |
|---|---|---|
| 24c | Mg | 1/8 0 1/4 |
| 16a | Al | 0 0 0 |
| 24d | Si | 3/8 0 1/4 |
| 96h | O | .03284 .05014 .65330 |
Data from G. V. Gibbs and J. V. Smith, "Refinement of the crystal structure of synthetic pyrope." American Mineralogist 50 2023-2039 (1965), PDF.
In [1]:
import sys
sys.path.extend(['../'])
import numpy as np
import onsager.crystal as crystal
import onsager.OnsagerCalc as onsager
Create garnet crystal (lattice constant in nm). Wyckoff positions cut and pasted from Bilbao crystallographic server.
In [2]:
# a0 = 1.1459
# alatt = a0*np.eye(3)
a0 = 1.
alatt = a0*np.array([[-0.5,0.5,0.5],[0.5,-0.5,0.5],[0.5,0.5,-0.5]])
invlatt = np.array([[0,1,1],[1,0,1],[1,1,0]])
x,y,z = (.03284,.05014,.65330)
uMg = ((1/8,0,1/4),(3/8,0,3/4),(1/4,1/8,0),(3/4,3/8,0),
(0,1/4,1/8),(0,3/4,3/8),(7/8,0,3/4),(5/8,0,1/4),
(3/4,7/8,0),(1/4,5/8,0),(0,3/4,7/8),(0,1/4,5/8))
uAl = ((0,0,0),(1/2,0,1/2),(0,1/2,1/2),(1/2,1/2,0),
(3/4,1/4,1/4),(3/4,3/4,3/4),(1/4,1/4,3/4),(1/4,3/4,1/4))
uSi = ((3/8,0,1/4),(1/8,0,3/4),(1/4,3/8,0),(3/4,1/8,0),
(0,1/4,3/8),(0,3/4,1/8),(3/4,5/8,0),(3/4,3/8,1/2),
(1/8,1/2,1/4),(7/8,0,1/4),(0,1/4,7/8),(1/2,1/4,1/8))
uO = ((x,y,z),(-x+1/2,-y,z+1/2),(-x,y+1/2,-z+1/2),(x+1/2,-y+1/2,-z),
(z,x,y),(z+1/2,-x+1/2,-y),(-z+1/2,-x,y+1/2),(-z,x+1/2,-y+1/2),
(y,z,x),(-y,z+1/2,-x+1/2),(y+1/2,-z+1/2,-x),(-y+1/2,-z,x+1/2),
(y+3/4,x+1/4,-z+1/4),(-y+3/4,-x+3/4,-z+3/4),(y+1/4,-x+1/4,z+3/4),(-y+1/4,x+3/4,z+1/4),
(x+3/4,z+1/4,-y+1/4),(-x+1/4,z+3/4,y+1/4),(-x+3/4,-z+3/4,-y+3/4),(x+1/4,-z+1/4,y+3/4),
(z+3/4,y+1/4,-x+1/4),(z+1/4,-y+1/4,x+3/4),(-z+1/4,y+3/4,x+1/4),(-z+3/4,-y+3/4,-x+3/4),
(-x,-y,-z),(x+1/2,y,-z+1/2),(x,-y+1/2,z+1/2),(-x+1/2,y+1/2,z),
(-z,-x,-y),(-z+1/2,x+1/2,y),(z+1/2,x,-y+1/2),(z,-x+1/2,y+1/2),
(-y,-z,-x),(y,-z+1/2,x+1/2),(-y+1/2,z+1/2,x),(y+1/2,z,-x+1/2),
(-y+1/4,-x+3/4,z+3/4),(y+1/4,x+1/4,z+1/4),(-y+3/4,x+3/4,-z+1/4),(y+3/4,-x+1/4,-z+3/4),
(-x+1/4,-z+3/4,y+3/4),(x+3/4,-z+1/4,-y+3/4),(x+1/4,z+1/4,y+1/4),(-x+3/4,z+3/4,-y+1/4),
(-z+1/4,-y+3/4,x+3/4),(-z+3/4,y+3/4,-x+1/4),(z+3/4,-y+1/4,-x+3/4),(z+1/4,y+1/4,x+1/4))
# tovec = lambda x: np.array(x)
# tovec2 = lambda x: np.array((x[0]+1/2,x[1]+1/2,x[2]+1/2))
tovec = lambda x: np.dot(invlatt, x)
pyrope = crystal.Crystal(alatt, [[vec(w) for w in ulist for vec in (tovec,)]
for ulist in (uMg, uAl, uSi, uO)],
['Mg','Al','Si','O'])
# print(pyrope)
Next, we construct a diffuser based on vacancies for our Mg ion. We need to create a sitelist (which will be the Wyckoff positions) and a jumpnetwork for the transitions between the sites. There are tags that correspond to the unique states and transitions in the diffuser. The first cutoff is $\sim 0.31a_0$, but that connects half of the Mg cation sites to each other; increasing the cutoff to $\sim 0.51a_0$ introduces a second network that completes the connections.
In [3]:
chem = 0 # 0 is the index corresponding to our Mg atom in the crystal
cutoff = 0.31*a0 # had been 0.51*a0
sitelist = pyrope.sitelist(chem)
jumpnetwork = pyrope.jumpnetwork(chem, cutoff)
Mgdiffuser = onsager.VacancyMediated(pyrope, chem, sitelist, jumpnetwork, 1)
print(Mgdiffuser)
Quick analysis on our jump network:
In [4]:
for jlist in jumpnetwork:
Z = 0
dx2 = np.zeros((3,3))
for (i,j), dx in jlist:
if i==0:
Z += 1
dx2 += np.outer(dx,dx)
print("coordination number:", Z)
print(dx2)
print("1/3 Tr dx dx:", dx2.trace()/3)
print("dx^2:", np.dot(dx,dx))
Next, we assemble our data: the energies and prefactors, for a VMg in pyrope for our representative states and transitions: these are the first states in the lists, which are also identified by the tags above. As we are computing a tracer, we make the choice to set $\nu_0 = 1/Z$ where $Z=4$ is the coordination number.
In [5]:
nu0 = 0.25
Etrans = 0.
# we don't need to use the tags, since there's only one site and jump type, and
# we want to build a tracer.
Mgthermodict = {'preV': np.ones(len(sitelist)),
'eneV': np.zeros(len(sitelist)),
'preT0': nu0*np.ones(len(jumpnetwork)),
'eneT0': Etrans*np.ones(len(jumpnetwork))}
Mgthermodict.update(Mgdiffuser.maketracerpreene(**Mgthermodict))
for k,v in Mgthermodict.items():
print('{}: {}'.format(k, v))
We compute the Onsager matrices, and look at $-L_\text{ss}/L_\text{sv}$ to get our correlation coefficient.
Note: we can define $f$ (for our tracer) as the ratio of $L_\text{ss}$ to $Z (\delta x)^2 w_2 c_\text{v}c_\text{s}/6 = \frac{1}{16}\nu_0 a_0^2$ in this case, the same as what we get for $L_\text{vv}$ and $-L_\text{sv}$.
In [6]:
Lvv, Lss, Lsv, L1vv = Mgdiffuser.Lij(*Mgdiffuser.preene2betafree(1, **Mgthermodict))
print(Lvv)
print(Lss)
print(Lsv)
print(L1vv)
print("Correlation coefficient:", -Lss[0,0]/Lsv[0,0])
Compare with tabulated GF data from Carlson and Wilson paper. They use the notation $(l,m,n)$ for a $\mathbf{\delta x}$ vector that is $a_0(l\hat x+m\hat y+n\hat z)/8$. We will need to find a corresponding site that lands at that displacement from our origin site.
Unfortunately, it looks like in two cases ((800), (444)) there are two distinct sites that are mapped in that displacement vector, which have different GF values; the CW reported values appear to be the averaged values. In two other cases, ((640), (420)) the reported values are half of what the computed values are here.
As Carlson and Wilson used a stochastic approach to compute their GF values, all of their other data has errors $\sim 10^{-4}$.
In [7]:
# tabulated data from paper
CarlsonWilsonGFdata = \
{(0,0,0): 2.30796022, (2,1,1): 1.30807261, (3,3,2): 0.80669536,
(4,2,0): 0.40469085, (4,4,4): 0.50242046, (5,3,2): 0.56195744,
(6,1,1): 0.56071092, (6,4,0): 0.22460654, (6,5,3): 0.42028488,
(6,5,5): 0.40137897, (7,2,1): 0.44437878, (8,0,0): 0.41938675}
In [8]:
print('CW index\tdx match\tGF (FT eval)\tGF(CW stoch.)\terror')
GF = Mgdiffuser.GFcalc # get our GF calculator; should already have rates set
basis = pyrope.basis[chem]
x0 = np.dot(alatt, basis[0])
for vec,gCW in CarlsonWilsonGFdata.items():
dx0 = np.array(vec,dtype=float)/8
nmatch, Gave, Gmatch = 0, 0, {}
for g in pyrope.G:
dx = np.dot(g.cartrot, dx0)
j = pyrope.cart2pos(x0+dx)[1]
if j is not None and j[0]==chem and j[1]<6:
G = GF(0, j[1], dx)
Gmatch[tuple((8*dx).astype(int))] = G
nmatch += 1
Gave += G
Gave /= nmatch
for t,G in Gmatch.items():
print('{}\t{}\t{:.12f}\t{:.8f}\t{:.4e}'.format(vec, t, -G, gCW, abs(G+gCW)))
print('{}\taverage value\t{:.12f}\t{:.8f}\t{:.4e}'.format(vec, -Gave, gCW, abs(Gave+gCW)))