\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
\errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage[]{algorithm2e}
\begin{document}
\title{Global energy minimization of multi-component phases with internal degrees of freedom}
\author{Richard Otis}
\affil{Affiliation not available}
\date{\today}
\maketitle
\selectlanguage{english}
\begin{abstract}
Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enimad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.%
\end{abstract}%
\section{Introduction}\label{introduction}
Miscibility gap detection is a crucial feature in thermodynamic
calculation software to accurately calculate the energy of phases
containing regions of compositional instability, i.e., spinodals, and is
commonly handled through global minimization (GM) of the Gibbs energy.
Inside the spinodal region, all points on the energy surface have
negative curvature and will be thermodynamically driven towards
demixing. The cause of miscibility gaps in non-ideal solutions is the
presence of unfavorable interactions between components that overwhelm
the entropically-driven ideal mixing contribution to the Gibbs energy.
From a computational perspective miscibility gaps pose a challenge
because they mean that the same phase may appear on multiple points on
the equilibrium tie hyperplane but at different compositions. In these
cases the software must increase the total degrees of freedom by
creating multiple composition sets of the same phase, potentially up to
the limit specified by the Gibbs phase rule. For multi-component systems
the topology of the energy surfaces can become quite complex, since
high-dimensional tangent hyperplanes can interact with the energy
surface. Moreover, when handling phases with internal degrees of
freedom, i.e., sublattices, it is possible for points on the global
energy surface (the composition) to be close together while being far
apart in their internal coordinates (the constitution). This is referred
to as an internal miscibility gap and is the cause of order-disorder
phase transitions.
While few authors have discussed fully generalized GM schemes, there has
been work in efficient sampling in low dimensions \cite{Emelianenko_2006} and
tie hyperplane calculation for the multi-component case
\cite{Perevoshchikova_2012}. However, there has not previously been detailed
discussion of methods for solving the multi-component case with multiple
sublattices.
Reliable convergence to the global energy minimum also requires a good
choice of starting points, i.e., phase compositions, for the
minimization procedure.
\section{Efficient sampling of composition
space}\label{efficient-sampling-of-composition-space}
Determinability is a desirable property from a reproducibility
perspective because it ensures that, in principle, users with the same
version of the software will always converge to exactly the same
solution. Quasi-random sampling is a technique that attempts to combine
the statistical advantages of pseudo-random sampling with a fully
deterministic algorithm. A quasi-random sequence is a deterministic
sequence of ``random-like'' values that can be generated by an
algorithm. Unlike psuedo-random number generators, quasi-random
sequences are explicitly designed for What makes quasi-random sequences
superior to pseudo-random numbers for sampling is that one can construct
quasi-random sequences which are biased toward ``spreading out,'' i.e.,
covering more of the domain. The Halton sequence is a quasi-random
sequence of numbers generated over \((0,1)\) using a prime
number as a base \cite{Halton_1964}. For \(N \lt 5\), where
\(N\) is the number of degrees of freedom, the
\(N\)-dimensional Halton sequence is low-discrepancy,
meaning it is very similar to the uniform distribution in terms of
domain coverage over \((0,1)^N\). For larger \(N\),
a scheme for permuting the coefficients should be used to prevent linear
correlations between dimensions \cite{hess2003performance, Chi2005}. Because values in
the Halton sequence can be computed out of order, mesh refinement
routines can continually draw from the sequence by keeping track of the
index of the last computed member.
The Halton sequence alone does not satisfy our purposes for sampling
composition space because it samples the \(N\)-hypercube,
implying that all \(N\) dimensions are independent, but
each sublattice's degrees of freedom must sum to unity in order to obey
the site fraction balance constraint. Instead the
\((N-1)\)-simplex must be sampled within each sublattice. The
most straightforward approach is to normalize each generated point so
that each sublattice's coordinates sum to unity, but this will tend to
oversample the middle of the simplices (Figure \ref{fig:normalized}).
Woronow \cite{Woronow1993} showed that uniformly distributed variates on
the \((N-1)\)-simplex follow a \(N\)-dimensional
symmetric Dirichlet distribution, which for the case of the full
composition range simplifies to \(N\) samples from a
standard exponential distribution divided by their sum. Because the
Halton sequence has low discrepancy, it can be substituted for the
uniform distribution. Exponentially distributed variates can be
generated with the inverse transformation \(-\ln(U)\), where
\(U\) is a uniform random or quasi-random variate. For
phases with multiple sublattices, all degrees of freedom are sampled
independently and then each sublattice's degrees of freedom grouped and
normalized to unity. This is the approach used in the present work and
is summarized as follows.
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
Draw the first \(K\) points from a
\(N\)-dimensional scrambled Halton sequence using the
first \(N\) primes. Store the result in a
\(K \times N\) matrix, \(M\).
\item
Compute \(-\ln(M)\) elementwise and store as the new
\(M\).
\item
Group together columns corresponding to each sublattice. Compute the
rowwise sum of each group and divide each element by its corresponding
group's sum.
\end{enumerate}
All end-members for each phase are also appended to the result of the
sampling algorithm, with the exception of pure-vacancy end-members.
Another approach called rejection sampling, which involves sampling the
\(N\)-hypercube and throwing out all points that fall
outside the \((N-1)\)-simplex, is not considered in this work.
The problem is that the fraction of rejected points quickly increases
with \(N\), making it very inefficient for multi-component
systems.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/normalized5/normalized5}
\caption{{\label{fig:normalized} Normalized exponential distributions produce better coverage of composition space versus the uniform distribution.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/halton4/halton4}
\caption{{\label{fig:halton} The standard Halton sequence has problems with
correlation between some dimensions. Deterministic scrambling alleviates
this.%
}}
\end{center}
\end{figure}
\section{Numerical Approach}\label{numerical-approach}
The general Gibbs energy minimization problem can be stated in the
following way, making use of the compound energy formalism to decompose
phases into independent sublattices. (A similar construction for
Helmholtz energy requires only a change of variables.)
\[
\textrm{min } G_\textrm{m}(T, P, f^i, y^i_{s, j}) \textrm{ where } c_k (T, P, f^i, y^i_{s, j}) = 0
\]
\(G_\textrm{m}\) is the molar Gibbs energy of the system,
\(T\) is the temperature, \(P\) is the
pressure, \(f^i\) is the fraction of phase \(i\)
and \(y^i_{s, j}\) is the site occupation fraction of component
\(j\) in sublattice \(s\) of phase
\(i\). \(c_k\) represents all equality
constraints. For example, fixing the temperature generates a constraint
\(T - \tilde{T} = 0\), where \(\tilde{T}\) is the user-specified
temperature. Site fraction balance constraints, \(\sum_s y^i_{s, j} - 1 = 0\), are
always present for each sublattice. Other constraints such as charge
balance can be constructed in the same way.
This formulation only considers equality constraints explicitly. The
inequality constraints, \(f^i \gt 0 \) and \(y^i_{s, j} > 0 \), are
left implicit and instead handled by a separate step discussed later.
Consider the multi-component, multi-phase case when \(T\)
and \(P\) are fixed along with some combination fixed of
\(\mu_j\) and \(x_j\), the chemical potential and
mole fraction of components, respectively. The following general
approach to energy minimization can be applied.
\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
Sample a set of points in the internal degrees of freedom of each
phase and compute the molar Gibbs energy. Include refined points from
previous iterations, if any.
\item
We are converged if the largest change in chemical potentials from the
previous iteration is less than some value. We cannot use the change
in energy because of dilute cases.
\item
Map the internal degrees of freedom of each point to global
composition space using the relation \(x_j = \sum_s b_s \frac{y_{s, j}}{1-y_{s, Va}}\), where
\(y_{s, Va}\) is the site fraction of vacancies in sublattice
\(s\).
\item
Construct an initial \(J\)-simplex, where
\(J\) is the number of components. Each vertex is
located at a pure composition for each component, i.e., the first
\(J\) simplex coordinates comprise a \(J \times J\)
identity matrix and therefore bound composition space. Each vertex
also has an energy coordinate. If any \(\mu_j\) are
specified the corresponding energy coordinate for that vertex is fixed
at that value. Otherwise, it can be set to any value, e.g., twice the
maximum energy (half if negative), as long as it causes the simplex to
lie strictly above the energy surface.
\item
Compute the driving force, \(\Delta G = G^z - \sum_j \mu_j x^z_j\) by computing the
distance of every point in the system to the current simplex. Find the
point \(z\) for which \(\Delta G\) is maximized.
The \(\mu_j\) here are the hyperplane coefficients of the
candidate simplex.
\item
Replace each vertex of the candidate simplex by the point
\(z\), one at a time, until a new simplex satisfying the
mass balance constraints is found. This is the new candidate simplex.
\item
Recompute the chemical potentials of the candidate simplex.
\item
Refine the internal degrees of freedom of the candidate simplex
vertices subject to the new potentials by solving a Newton-Raphson
sub-problem.
\item
Return to Step 1, including all our refined points in the next
iteration.
\end{enumerate}
Table of constraint types and formulae
\subsection{Alternative Formulation}\label{alternative-formulation}
The refinement sub-problem can also be solved by constructing an
augmented Hessian matrix including all phases' degrees of freedom and
their constraints. Iterative steps toward the solution are then computed
using the Newton-Raphson method. This is sometimes called the
``Newton-Lagrange'' method and is detailed in standard texts, e.g.,
section 18.1 of \cite{Nocedal2006}. A similar approach for phase diagram
calculation was reported by Lukas \cite{Lukas1982}.
\[ F(T, P, f^{i}, y^i_{s, j}) = \sum_i f^{i}G^{i}_{m}(T, P, y^i_{s, j})\] \[ L(T, P, f^{i}, y^i_{s, j}, \lambda_k) = F(T, P, f^{i}, y^i_{s, j}) - \sum_k \lambda_k c_k(T, P, f^{i}, y^i_{s, j})\] \[ \left[ \begin{array}{ccc}
W_k & -A^T_k \\\
A_k & 0 \end{array} \right]\left[\begin{array}{ccc}
p_k \\\
p_\lambda \end{array} \right] =
\left[\begin{array}{ccc}
-\nabla F_k + A^T_k \lambda_k \\\
-c_k \end{array} \right]\]
\(p\) is the Newton step direction for the next iteration.
\(W_k\) is the Hessian matrix of \(L\) with
respect to all degrees of freedom except \(\lambda_k\).
\(A_k\) is the Jacobian matrix of the constraints,
\([\nabla c_1, \nabla c_2, ..., \nabla c_m]^T\). The key advantage of this approach is you can
optimize the degrees of freedom of all phases, including phase fractions
and chemical potentials, simultaneously. The disadvantages are most
apparent in large multi-component systems and systems with miscibility
gaps; in principle \(c+2\) copies of \emph{each phase} in
the system need to be entered or else some procedure for adding and
removing composition sets is required, in which case the system of
equations must be rebuilt every time the stable set of phases changes.
However, this approach can be useful if the stable set of phases is
already known.
\subsection{Algorithms}
\begin{algorithm}[H]
\KwData{this text}
\KwResult{how to write algorithm with \LaTeX2e }
initialization\;
\While{not at end of this document}{
read current\;
\eIf{understand}{
go to next section\;
current section becomes this one\;
}{
go back to the beginning of current section\;
}
}
\caption{{How to write algorithms}}
\end{algorithm}
\selectlanguage{english}
\FloatBarrier
\nocite{*}
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}
\end{document}