Free Energy Curvature

Calculating the interdiffusivity from atomic mobility requires the chemical potential gradient which can be determined from the curvature of the free energy surface.

$$ \frac{\partial \mu_i}{\partial x_j} = \frac{\partial^2 G_m^\alpha}{\partial x_i \partial x_j} $$

For models representing solid solutions with a single sublattice, the curvature of the Gibbs free energy can be easily found. However, for most other models, the chemical potential gradient must be found numerically. This is done by first solving for equilibrium, then using chain rule differentiation to get the derivative with respect to any parameter.

Solving equilibrium for a single phase system is done by minimizing the Gibbs free energy subjected to mass balance and summation constraints ($\boldsymbol g$).

$$ \mathcal{L} = n^\alpha G_m^\alpha + \sum_{i}{\Lambda_i g_i} $$

For a single phase system, assuming no charged species, the Langragian becomes the following:

$$ \mathcal{L} = n^\alpha G_m^\alpha + \sum_{A}{\mu_A \left(x_A - n^\alpha \frac{M_A^\alpha}{M}\right)} + \sum_{s}{\lambda^s \left(1 - \sum_{i \epsilon s}{y_i^s} \right)} $$

The terms in the above equation can be divided into two groups: the fixed parameters (z) and the degrees of freedom ($\boldsymbol\theta$). The fixed parameters include temperature (T), pressure (P) and composition ($x_A$). The degrees of freedom include the site fractions ($y_i^s$), phase fraction ($n^\alpha$), chemical potential ($\mu_A$) and summation constraint multiplier ($\lambda^s$).

The $\boldsymbol\theta$ group will be further divided into $\boldsymbol{\mathcal{v}}$ (which includes $y_i^s$ and $n^\alpha$) and $\boldsymbol \Lambda$ (which includes the Langragian multipliers, $\mu_A$ and $\lambda^s$).

Solving for equilibrium involves finding the solution to the following equations.

$$ \frac{\partial \mathcal{L}}{\partial \mathcal{v}_i} = \frac{\partial}{\partial \mathcal{v}_i} (n^\alpha G_m^\alpha) + \sum_i{\Lambda_j \frac{\partial g_j}{\partial \mathcal{v}_i}} = 0 $$

$$ \frac{\partial \mathcal{L}}{\partial \Lambda_j} = g_j = 0 $$

Once equilibrium is found, a system of equations can be constructed using chain rule derivative to find the derivatives of $\boldsymbol\theta$ with respect to z.

$$ \sum_{j}{\frac{\partial^2 \mathcal{L}}{\partial \theta_i \partial \theta_j} \frac{\partial \theta_j}{\partial z_k}} = \frac{\partial^2 \mathcal{L}}{\partial \theta_i \partial z_k} $$

$$ \left[ \begin{array}{cc} \frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{\mathcal{v}}^2} & \frac{\partial \boldsymbol{g}}{\partial \boldsymbol{\mathcal{v}}} \newline \left(\frac{\partial \boldsymbol{g}}{\partial \boldsymbol{\mathcal{v}}} \right)^T & 0 \end{array} \right] \left[ \begin{array}{c} \frac{\partial \boldsymbol{\mathcal{v}}}{\partial z} \newline \frac{\partial \boldsymbol{\Lambda}}{\partial z} \end{array} \right] = \left[ \begin{array}{c} -\sum{\Lambda_i \frac{\partial^2 g_i}{\partial \boldsymbol{v} \partial z}} \newline \frac{\partial \boldsymbol{g}}{\partial z} \end{array} \right] $$

Single phase systems #

For a single phase system, the above system of equations become

$$ \left[ \begin{array}{cccc} \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y}^2} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right) & \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right) \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right)^T & 0 & 0 & \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right) \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right)^T & 0 & 0 & 0 \newline \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right)^T & \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right)^T & 0 & 0 \end{array} \right] \left[ \begin{array}{c} \frac{\partial \boldsymbol{y}}{\partial x_i} \newline \frac{\partial n^\alpha}{\partial x_i} \newline \frac{\partial \boldsymbol{\lambda}}{\partial x_i} \newline \frac{\partial \boldsymbol{\mu}}{\partial x_i} \end{array} \right] = \left[ \begin{array}{c} 0 \newline \vdots \newline 0 \newline \frac{\partial x_A}{\partial x_i} \newline \vdots \newline \frac{\partial x_N}{\partial x_i} \end{array} \right] $$

Where

$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y}^2} \right)_{ij} = n^\alpha \frac{\partial^2 G_m^\alpha}{\partial y_i \partial y_j} $$

$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial n^\alpha} \right)_i = \frac{\partial G_m^\alpha}{\partial y_i} - \sum_A{\mu_A \frac{\partial}{\partial y_i} \left(\frac{M_A^\alpha}{M} \right)} $$

$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\lambda}} \right)_{ij} = -1 \text{ if } y_i \epsilon \text{ sublattice j else } 0 $$

$$ \left(\frac{\partial^2 \mathcal{L}}{\partial \boldsymbol{y} \partial \boldsymbol{\mu}} \right)_{ij} = -n^\alpha \frac{\partial}{\partial y_i} \left(\frac{M_j^\alpha}{M} \right) $$

$$ \left(\frac{\partial^2 \mathcal{L}}{\partial n^\alpha \partial \boldsymbol{\mu}} \right)_j = -\frac{M_j^\alpha}{M} $$

$$ \frac{\partial x_A}{\partial x_i} = \delta_{Ai} $$

Defining the chemical potential gradient #

The FreeEnergyHessian module in kawin provides two methods: dMudX and partialdMudX. These methods give different definitions of the derivative of the chemical potential. Additionally, it will give different values than if the chemical potential gradient is determined using a finite difference approach using thermodynamic software. The partialdMudX method outputs a ($\text{n x n}$) matrix of the partial derivatives of the chemical potential with respect to the composition, treating all components as independent. Using a finite difference method to calculate the chemical potential gradient will calculate the total derivative (with reference element n) since most thermodynamic software only require n-1 components to be defined. This will give a ($\text{n x n-1}$) matrix. The dMudX method gives a ($\text{n-1 x n-1}$) matrix, where each index is the total derivative of the chemical potential with reference element n subtracted by the total derivative of the chemical potential of the reference element. This is more equivalent to the curvature (Laplacian) of the Gibbs free energy surface with the reference element. This method is required for finding $M^{-1}$ for the multicomponent growth rate.

$$ \text{partialdMudX} = \frac{\partial \mu_i}{\partial x_j} $$

$$ \text{dMudX} = \frac{\partial}{\partial x_j} (\mu_i - \mu_n) = (\frac{\partial \mu_i}{\partial x_j} - \frac{\partial \mu_i}{\partial x_n}) - (\frac{\partial \mu_n}{\partial x_j} - \frac{\partial \mu_n}{\partial x_n}) $$

$$ \text{finite difference} = \frac{\partial \mu_i}{\partial x_j} - \frac{\partial \mu_i}{\partial x_n} $$

References #

1. H. Larsson and M. Jansson, “Rate of change at equilibrium” Calphad: Computer Coupling of Phase Diagrams and Thermochemistry 51 (2015) p. 220