<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Fe-C diffusion and elastodiffusivity
Taking data from R.G.A. Veiga, M. Perez, C. Becquart, E. Clouet and C. Domain, Acta mater. 59 (2011) p. 6963 doi
Step1:  Create BCC lattice (lattice constant in nm).
Step2:  Elastic constants converted from GPa ($10^9$ J/m$^3$) to eV/(atomic volume).
Step3:  Add carbon interstitial sites at octahedral sites in the lattice. This code (1) gets the set of symmetric Wyckoff positions corresponding to the single site $[00\frac12]$ (first translated into unit cell coordinates), and then adds that new basis to our Fe crystal to generate a new crystal structure, that we name "FeC".
Step4:  Next, we construct a diffuser based on our interstitial. 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.
Step5:  Next, we assemble our data
Step6:  We look at the diffusivity $D$, the elastodiffusivity $d$, and the activation volume tensor $V$ over a range of temperatures from 300K to 1200K.
First, we calculate all of the pieces, including the diffusivity prefactor and activation barrier. As we only have one barrier, we compute the barrier at $k_\text{B}T = 1$.
Step7:  Activation volume. We plot the isotropic value (change in diffusivity with respect to pressure), but also the $V_{xxxx}$, $V_{xxyy}$, and $V_{xyxy}$ terms. Interestingly, the $V_{xxxx}$ term is negative---which indicates that diffusivity along the $[100]$ direction increases with compressive stress in the $[100]$ direction.
<ASSISTANT_TASK:>
Python Code:
import sys
sys.path.extend(['../'])
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')
%matplotlib inline
import onsager.crystal as crystal
import onsager.OnsagerCalc as onsager
from scipy.constants import physical_constants
kB = physical_constants['Boltzmann constant in eV/K'][0]
Explanation: Fe-C diffusion and elastodiffusivity
Taking data from R.G.A. Veiga, M. Perez, C. Becquart, E. Clouet and C. Domain, Acta mater. 59 (2011) p. 6963 doi:10.1016/j.actamat.2011.07.048
Fe in the body-centered cubic phase, $a_0 = 0.28553\text{ nm}$; C sit at octahedral sites, where the transition states between octahedral sites are represented by tetrahedral sites. The data is obtained from an EAM potential, where $C_{11} = 243\text{ GPa}$, $C_{12}=145\text{ GPa}$, and $C_{44} = 116\text{ GPa}$. The tetrahedral transition state is 0.816 eV above the octahedral site, and the attempt frequency is taken as 10 THz ($10^{13}\text{ Hz}$).
The dipole tensors can be separated into parallel and perpendicular components; the parallel direction points towards the closest Fe atoms for the C, while the perpendicular components lie in the interstitial plane. For the octahedral, the parallel component is 8.03 eV, and the perpendicular is 3.40 eV; for the tetrahedral transition state, the parallel component is 4.87 eV, and the perpendicular is 6.66 eV.
End of explanation
a0 = 0.28553
Fe = crystal.Crystal.BCC(a0, "Fe")
print(Fe)
Explanation: Create BCC lattice (lattice constant in nm).
End of explanation
stressconv = 1e9*1e-27*Fe.volume/physical_constants['electron volt'][0]
c11, c12, c44 = 243*stressconv, 145*stressconv, 116*stressconv
s11, s12, s44 = (c11+c12)/((c11-c12)*(c11+2*c12)), -c12/((c11-c12)*(c11+2*c12)), 1/c44
print('S11 = {:.4f} S12 = {:.4f} S44 = {:.4f}'.format(s11, s12, s44))
stensor = np.zeros((3,3,3,3))
for a in range(3):
    for b in range(3):
        for c in range(3):
            for d in range(3):
                if a==b and b==c and c==d: stensor[a,b,c,d] = s11
                elif a==b and c==d: stensor[a,b,c,d] = s12
                elif (a==d and b==c) or (a==c and b==d): stensor[a,b,c,d] = s44/4
Explanation: Elastic constants converted from GPa ($10^9$ J/m$^3$) to eV/(atomic volume).
End of explanation
uoct = np.dot(Fe.invlatt, np.array([0, 0, 0.5*a0]))
FeC = Fe.addbasis(Fe.Wyckoffpos(uoct), ["C"])
print(FeC)
Explanation: Add carbon interstitial sites at octahedral sites in the lattice. This code (1) gets the set of symmetric Wyckoff positions corresponding to the single site $[00\frac12]$ (first translated into unit cell coordinates), and then adds that new basis to our Fe crystal to generate a new crystal structure, that we name "FeC".
End of explanation
chem = 1  # 1 is the index corresponding to our C atom in the crystal
sitelist = FeC.sitelist(chem)
jumpnetwork = FeC.jumpnetwork(chem, 0.6*a0)  # 0.6*a0 is the cutoff distance for finding jumps
FeCdiffuser = onsager.Interstitial(FeC, chem, sitelist, jumpnetwork)
print(FeCdiffuser)
Explanation: Next, we construct a diffuser based on our interstitial. 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.
End of explanation
Dconv = 1e-2
vu0 = 10*Dconv
Etrans = 0.816
dipoledict = {'Poctpara': 8.03, 'Poctperp': 3.40, 
              'Ptetpara': 4.87, 'Ptetperp': 6.66}
FeCthermodict = {'pre': np.ones(len(sitelist)), 'ene': np.zeros(len(sitelist)),
                 'preT': vu0*np.ones(len(jumpnetwork)), 
                 'eneT': Etrans*np.ones(len(jumpnetwork))}
# now to construct the site and transition dipole tensors; we use a "direction"--either
# the site position or the jump direction--to determine the parallel and perpendicular
# directions.
for dipname, Pname, direction in zip(('dipole', 'dipoleT'), ('Poct', 'Ptet'), 
                                     (np.dot(FeC.lattice, FeC.basis[chem][sitelist[0][0]]),
                                      jumpnetwork[0][0][1])):
    # identify the non-zero index in our direction:
    paraindex = [n for n in range(3) if not np.isclose(direction[n], 0)][0]
    Ppara, Pperp = dipoledict[Pname + 'para'], dipoledict[Pname + 'perp']
    FeCthermodict[dipname] = np.diag([Ppara if i==paraindex else Pperp 
                                      for i in range(3)])
for k,v in FeCthermodict.items():
    print('{}: {}'.format(k, v))
Explanation: Next, we assemble our data: the energies, prefactors, and dipoles for the C atom in Fe, matched to the representative states: these are the first states in the lists, which are also identified by the tags above.
A note about units: If $\nu_0$ is in THz, and $a_0$ is in nm, then $a_0^2\nu_0 = 10^{-2}\text{ cm}^2/\text{s}$. Thus, we multiply by Dconv = $10^{-2}$ so that our diffusivity is output in cm<sup>2</sup>/s.
End of explanation
Trange = np.linspace(300, 1200, 91)
Tlabels = Trange[0::30]
Dlist, dDlist, Vlist = [], [], []
for T in Trange:
    beta = 1./(kB*T)
    D, dD = FeCdiffuser.elastodiffusion(FeCthermodict['pre'], 
                                        beta*FeCthermodict['ene'], 
                                        [beta*FeCthermodict['dipole']],
                                        FeCthermodict['preT'], 
                                        beta*FeCthermodict['eneT'], 
                                        [beta*FeCthermodict['dipoleT']])
    Dlist.append(D[0,0])
    dDlist.append([dD[0,0,0,0], dD[0,0,1,1], dD[0,1,0,1]])
    Vtensor = (kB*T/(D[0,0]))*np.tensordot(dD, stensor, axes=((2,3),(0,1)))
    Vlist.append([np.trace(np.trace(Vtensor))/3, 
                  Vtensor[0,0,0,0], Vtensor[0,0,1,1], Vtensor[0,1,0,1]])
D0 = FeCdiffuser.diffusivity(FeCthermodict['pre'], 
                             np.zeros_like(FeCthermodict['ene']), 
                             FeCthermodict['preT'], 
                             np.zeros_like(FeCthermodict['eneT']))
D, dbeta = FeCdiffuser.diffusivity(FeCthermodict['pre'], 
                                   FeCthermodict['ene'], 
                                   FeCthermodict['preT'], 
                                   FeCthermodict['eneT'],
                                   CalcDeriv=True)
print('D0: {:.4e} cm^2/s\nEact: {:.3f} eV'.format(D0[0,0], dbeta[0,0]/D[0,0]))
D, dD = np.array(Dlist), np.array(dDlist)
d11_T = np.vstack((Trange, dD[:,0])).T
d11pos = np.array([(T,d) for T,d in d11_T if d>=0])
d11neg = np.array([(T,d) for T,d in d11_T if d<0])
fig, ax1 = plt.subplots()
ax1.plot(1./Trange, D, 'k', label='$D$')
# ax1.plot(1./Trange, dD[:,0], 'b', label='$d_{11}$')
ax1.plot(1./d11pos[:,0], d11pos[:,1], 'b', label='$d_{11}$')
ax1.plot(1./d11neg[:,0], -d11neg[:,1], 'b--')
ax1.plot(1./Trange, dD[:,1], 'r', label='$d_{12}$')
ax1.plot(1./Trange, dD[:,2], 'g-.', label='$d_{44} = D$')
ax1.set_yscale('log')
ax1.set_ylabel('$D$ [cm$^2$/s]', fontsize='x-large')
ax1.set_xlabel('$T^{-1}$ [K$^{-1}$]', fontsize='x-large')
ax1.legend(bbox_to_anchor=(0.15,0.15,0.2,0.4), ncol=1, 
           shadow=True, frameon=True, fontsize='x-large')
ax2 = ax1.twiny()
ax2.set_xlim(ax1.get_xlim())
ax2.set_xticks([1./t for t in Tlabels])
ax2.set_xticklabels(["{:.0f}K".format(t) for t in Tlabels])
ax2.set_xlabel('$T$ [K]', fontsize='x-large')
ax2.grid(False)
ax2.tick_params(axis='x', top='on', direction='in', length=6)
plt.show()
# plt.savefig('Fe-C-diffusivity.pdf', transparent=True, format='pdf')
d11pos[0,0], d11neg[-1,0]
Explanation: We look at the diffusivity $D$, the elastodiffusivity $d$, and the activation volume tensor $V$ over a range of temperatures from 300K to 1200K.
First, we calculate all of the pieces, including the diffusivity prefactor and activation barrier. As we only have one barrier, we compute the barrier at $k_\text{B}T = 1$.
End of explanation
V = np.array(Vlist)
fig, ax1 = plt.subplots()
ax1.plot(1./Trange, V[:,0], 'k', label='$V_{\\rm{total}}$')
ax1.plot(1./Trange, V[:,1], 'b', label='$V_{11}$')
ax1.plot(1./Trange, V[:,2], 'r', label='$V_{12}$')
ax1.plot(1./Trange, 2*V[:,3], 'g', label='$V_{44}$')
ax1.set_yscale('linear')
ax1.set_ylabel('$V$ [atomic volume]', fontsize='x-large')
ax1.set_xlabel('$T^{-1}$ [K$^{-1}$]', fontsize='x-large')
ax1.legend(bbox_to_anchor=(0.3,0.3,0.5,0.2), ncol=2, 
           shadow=True, frameon=True, fontsize='x-large')
ax2 = ax1.twiny()
ax2.set_xlim(ax1.get_xlim())
ax2.set_xticks([1./t for t in Tlabels])
ax2.set_xticklabels(["{:.0f}K".format(t) for t in Tlabels])
ax2.set_xlabel('$T$ [K]', fontsize='x-large')
ax2.grid(False)
ax2.tick_params(axis='x', top='on', direction='in', length=6)
plt.show()
# plt.savefig('Fe-C-activation-volume.pdf', transparent=True, format='pdf')
print('Total volume: {v[0]:.4f}, {V[0]:.4f}A^3\nV_xxxx: {v[1]:.4f}, {V[1]:.4f}A^3\nV_xxyy: {v[2]:.4f}, {V[2]:.4f}A^3\nV_xyxy: {v[3]:.4f}, {V[3]:.4f}A^3'.format(v=V[-1,:], V=V[-1,:]*1e3*Fe.volume))
Vsph = 0.2*(3*V[-1,1] + 2*V[-1,2] + 4*V[-1,3])  # (3V11 + 2V12 + 2V44)/5
print('Spherical average uniaxial activation volume: {:.4f} {:.4f}A^3'.format(Vsph, Vsph*1e3*Fe.volume))
Explanation: Activation volume. We plot the isotropic value (change in diffusivity with respect to pressure), but also the $V_{xxxx}$, $V_{xxyy}$, and $V_{xyxy}$ terms. Interestingly, the $V_{xxxx}$ term is negative---which indicates that diffusivity along the $[100]$ direction increases with compressive stress in the $[100]$ direction.
End of explanation
<END_TASK>