Class: Phase Field Crystal

In this class, we simulate a crystal using the phase-field crystal methodology. The phase-field crystal (PFC) model is a mesoscopic model that describes the dynamics of a crystal. The model is based on a free energy functional, which is minimized to find the equilibrium state of the system. The free energy functional is a function of the phase field \(\psi\), which is a real-valued field that describes the local density of the crystal. The PFC model is a powerful tool for studying the dynamics of crystals, and it has been used to study a wide range of phenomena, including the formation of defects, the dynamics of grain boundaries, and the growth of crystals.

file: comfit/models/phase_field_crystal.py 
class: PhaseFieldCrystal

Variables and parameters

The phase field \(\psi\).

pfc.psi

Basic model

The PFC methodology is based on postulating a free energy

The PFC free energy

\[ \mathcal F = \int d\boldsymbol{r} \tilde f(\psi, \nabla \psi, ...) \]

with the goal of the state that minimizes \(\mathcal F\) has a certain symmetry.
In this documentation, we will present six such models to represent different crystalline structures in 1, 2 and 3 dimensions, all of which are on the form

The PFC free energy density

\[ \tilde f( \psi, \nabla \psi, ...) = \frac{1}{2} (\mathcal L(\nabla) \psi)^2 + \frac{1}{2} \texttt{r} \psi^2 + \frac{1}{3} \texttt t \psi^3 + \frac{1}{4} \texttt v \psi^4, \]

where \(\mathcal L (\nabla)\) is a gradient operator dependent on the dimension and target symmetry, listed below

Model	Derivative operator \(\mathcal L\)
1D periodic	\(\mathcal L_1 = (1+\nabla^2)\)
2D triangular	\(\mathcal L_1 = (1+\nabla^2)\)
2D square	\(\mathcal L_1 \mathcal L_2 = (1+\nabla^2) (2+\nabla^2)\)
3D bcc	\(\mathcal L_1 = (1+\nabla^2)\)
3D fcc	\(\mathcal L_1 \mathcal L_{4/3} = (1+\nabla^2) (4/3+\nabla^2)\)
3D simple cubic	\(\mathcal L_1 \mathcal L_2 \mathcal L_3 = (1+\nabla^2) (2+\nabla^2) (3+\nabla^2)\)

Table: Phase-field crystal models. \(\mathcal L_X = (X+\nabla^2)\).

Of particular interest is the PFC chemical potential, given by

The PFC chemical potential

\[ \tilde \mu_c = \frac{\delta \mathcal F}{\delta \psi} = \left ( \mathcal L^2 \psi + \texttt r \psi + \texttt t \psi^2 + \texttt v \psi^3 \right ) \]

The PFC model was introduced inin Ref.¹. The primary field is \(\psi\), a real valued field

pfc.psi 

Crystal symmetry

The way in which this field can model a crystalline lattice is by the ground state having that particular symmetry. For instance, the figure below shows the (real) field \(\psi\) for a (a) 2D triangular, and a (b) 2D square lattice.

Figure: Ground state of the PFC model for (a) 2D triangular and (b) 2D square lattice.

The crystalline symmetry is described by a Bravais lattice, which consists of vectors that point to peaks on the lattice. The Fourier transform of a field that a Bravais lattice symmetry will itself have a lattice symmetry, which is called the reciprocal lattice, and is useful to express the field in terms of a Fourier series.

Figure: Bravais lattices \(\mathcal B\) and their reciprocal lattices \(\mathcal R\) for square and triangular symmetry. In each case, \(\{\mathbf a^{(n)}\}_{n=1}^2\) are primitive lattice vectors and \(\{ \mathbf q^{(n)}\}_{n=1}^2\) primitive reciprocal lattice vectors (RLVs) that satisfy \(\mathbf a^{(n)} \cdot \mathbf q^{(m)} = 2\pi \delta_{nm}\). Amended and reprinted from Ref. ² with permission.

Using the reciprocal lattice, one can express a field which has a Bravais lattice symmetry as a Fourier series in terms of the reciprocal lattice vectors

\[ \psi (\mathbf r) = \sum_{\mathbf q^{(n)} \in \mathcal R} \eta_n e^{\mathfrak i \mathbf q^{(n)} \cdot \mathbf r}, \]

the reciprocal lattice is \(d\) dimensions generated by \(d\) primitive RLVs \(\lbrace\boldsymbol{q}^{(n)}\rbrace_{n=1}^d\), satisfying

The primitive BLV-RLV orthogonality relation.

\[ \boldsymbol{q}^{(n)} \cdot \boldsymbol{a}^{(m)} = 2\pi \delta_{nm}, \quad n,m \leq d. \]

There are infinitely many BLVs and RLVs, but of particular interest are those with the smallest magnitude; which are sometimes named primary RLVs/BLVs. These values are encoded in the instance properties a,q, e.g.,

import comfit as cf

pfc = cf.PhaseFieldCrystal2DTriangular(1,1)
print(pfc.a)
print(pfc.q)
>> 
[[3.627, 6.283], [7.255, 0.0], [3.627, -6.283]]
[[0, 1], [0.866, -0.5], [-0.866, -0.5]]

---

Example: The 1D PFC - Crystal symmetry

The simplest example of a PFC model is the 1D periodic model, which simply looks like a sine wave.

The crystal symmetry in this case refers to \(\psi\) having a periodicity of \(a_0=2\pi\), and so the Bravais lattice is simply \(\{ a^{(n)} \} = \{...,-2\pi,-\pi,0,\pi,2\pi,... \}\). The primary BLVs are thus given by

\[ \mathcal B_{\textrm{per}}^{(1)} = \{ a^{(1)} = 2\pi,~ a^{(-1)} = -2\pi \} \]

The reciprocal lattice is in this case simply \(\{ q^{(n)} \} = \{...,q^{(-2)} =-2,~ q^{(-1)} = -1,~ q^{(0)} = 0,~ q^{(1)} = 1,~ q^{(2)} =2,... \}\), the primary set being

\[ \mathcal R_{\textrm{per}}^{(1)} = \{ q^{(1)} = 1,~ q^{(-1)} = -1 \} \]

and any field with that periodicity can be expressed as a Fourier series in terms of these RLVs by

\[ \psi(x) = \psi_0 + A (e^{\mathfrak iq^{(1)} x} + e^{\mathfrak i q^{(-1)} x}) + B (e^{\mathfrak iq^{(2)} x} + e^{\mathfrak i q^{(-2)} x}) + ..., \]

Cutting the Fourier series off at the primary RLVs mode is called the one-mode approximation

\[ \psi^{\textrm{eq}} \approx \psi_0 + A \sum_{q^{(n)} \in \mathcal R^{(1)}} e^{\mathfrak i q^{(n)} x}, \]

where \(\mathcal R_{\textrm{per}}^{(1)}\) as the primary RLVs.

Primitive vs primary LV

Note the subtle difference between a primitive LV and a primary LV. The primitive BLVs generate the lattice and satisfy the orthogonal relationship with the primitive RLVs which generate the reciprocal lattice. There are only \(d\) primitive LVs in \(d\) dimensions whereas there may be any number of primary lattice vectors. For instance, there are six primary BLVs for the triangular lattice and four primary lattice vectors for the 2D square lattice. In the following, the primary BLVs and RLVs have been chosen so that the first \(d\) primary LVs are primitive LVs that satisfy the orthogonality relation.

Below are all the lattice constants and primary RLVs and BLVs for the models we consider.

1D periodic2D Triangular2D square3D body-centered cubic3D face-centered cubic3D simple cubic

Lattice constant

\[ a_0 = 2\pi \]

Primary RLVs

\[ \mathcal R_{\textrm{per}}^{(1)} = \left \lbrace \begin{array}{l} q^{(1)} = 1 \\ q^{(-1)}= - q^{(1)} \end{array} \right \rbrace \]

Primary BLVs

\[ \mathcal B_{\textrm{per}}^{(1)} = \left \lbrace \begin{array}{l} a^{(1)} = a_0 \\ a^{(-1)} = -\mathbf a^{(1)} \end{array} \right \rbrace \]

Lattice constant

\[ a_0 = \frac{4\pi}{\sqrt 3} \]

Primary RLVs

\[ \mathcal R_{\textrm{tri}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf q^{(1)} = (\sqrt 3/2,-1/2) \\ \mathbf q^{(2)} = (0,1) \\ \mathbf q^{(3)} = (-\sqrt 3/2,-1/2) \\ \mathbf q^{(-n)} = - \mathbf q^{(n)}|_{n=1,2,3} \\ \end{array} \right \rbrace } \]

Primary BLVs

\[ \mathcal B_{\textrm{tri}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf a^{(1)} = a_0(1,0) \\ \mathbf a^{(2)} = a_0(1/2,\sqrt 3/2) \\ \mathbf a^{(3)} = a_0(1/2,-\sqrt 3/2) \\ \mathbf a^{(-n)}=-\mathbf a^{(n)}|_{n=1,2,3} \\ \end{array} \right \rbrace } \]

Lattice constant $$ a_0 = 2\pi $$

Primary RLVs

\[ \mathcal B_{\textrm{sq}}^{(1)} = \left \lbrace \begin{array}{l} \mathbf a^{(1)} = a_0(1,0) \\ \mathbf a^{(2)} = a_0(0,1) \\ \mathbf a^{(-1)},\mathbf a^{(-2)} \end{array} \right \rbrace \]

Primary Bravais Lattice vectors

\[ \mathcal R_{\textrm{sq}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf q^{(1)} = (1,0) \\ \mathbf q^{(2)} = (0,1) \\ \mathbf q^{(-n)} = - \mathbf q^{(n)} |_{n=1,2} \\ \end{array} \right \rbrace} \]

\[ \mathcal R_{\textrm{sq}}^{(2)} = \left \lbrace \begin{array}{l} \mathbf q^{(3)} = (1,-1) \\ \mathbf q^{(4)} = (1,1) \\ \mathbf q^{(-n)} = - \mathbf q^{(n)} |_{n=3,4} \\ \end{array} \right \rbrace \]

Lattice constant

\[ a_0 = 2\pi \sqrt {2} \]

Primary BLVs

\[ \mathcal B_{\textrm{bcc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf a^{(1)} = a_0(-1,1,1)/2 \\ \mathbf a^{(2)} = a_0(1,-1,1)/2 \\ \mathbf a^{(3)} = a_0(1,1,-1)/2\\ \mathbf a^{(4)} = a_0(1,1,1)/2\\ \mathbf a^{(-n)}=-\mathbf a^{(n)}|_{n=1,...,4} \\ \end{array} \right \rbrace } \]

Primary RLVs

\[ \mathcal R_{\textrm{bcc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf q^{(1)} = (0,1,1)/\sqrt 2 \\ \mathbf q^{(2)} = (1,0,1)/\sqrt 2 \\ \mathbf q^{(3)} = (1,1,0)/\sqrt 2\\ \mathbf q^{(4)} = (0,-1,1)/\sqrt 2\\ \mathbf q^{(5)} = (-1,0,1)/\sqrt 2\\ \mathbf q^{(6)} = (-1,1,0)/\sqrt 2\\ \mathbf q^{(-n)} = - \mathbf q^{(n)}|_{n=1,...6} \\ \end{array} \right \rbrace } \]

Lattice constant

\[ a_0 = 2\pi \sqrt 3 \]

Primary BLVs

\[ \mathcal B_{\textrm{fcc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf a^{(1)} = a_0(0,1,1)/2 \\ \mathbf a^{(2)} = a_0(1,0,1)/2 \\ \mathbf a^{(3)} = a_0(1,1,0)/2\\ \mathbf a^{(4)} = a_0(0,-1,1)/2\\ \mathbf a^{(5)} = a_0(-1,0,1)/2\\ \mathbf a^{(6)} = a_0(-1,1,0)/2\\ \mathbf a^{(-n)}=-\mathbf a^{(n)}|_{n=1,...,6} \\ \end{array} \right \rbrace } \]

Primary RLVs

\[ \mathcal R_{\textrm{fcc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf q^{(1)} = (-1,1,1)/\sqrt 3 \\ \mathbf q^{(2)} = (1,-1,1)/\sqrt 3 \\ \mathbf q^{(3)} = (1,1,-1)/\sqrt 3\\ \mathbf q^{(4)} = (1,1,1)/\sqrt 3\\ \mathbf q^{(-n)} = - \mathbf q^{(n)}|_{n=1,...,4}\\ \end{array} \right \rbrace } \]

\[ \mathcal R_{\textrm{fcc}}^{(4/3)} = { \left \lbrace \begin{array}{l} \mathbf q^{(5)} = (2,0,0)/\sqrt 3 \\ \mathbf q^{(6)} = (0,2,0)/\sqrt 3 \\ \mathbf q^{(7)} = (0,0,2)/\sqrt 3\\ \mathbf q^{(-n)} = - \mathbf q^{(n)} |_{n=5,6,7} \\ \end{array} \right \rbrace } \]

Lattice constant

\[ a_0 = 2\pi \]

Primary BLVs

\[ \mathcal B_{\textrm{sc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf a^{(1)} = a_0(1,0,0) \\ \mathbf a^{(2)} = a_0(0,1,0) \\ \mathbf a^{(3)} = a_0(0,0,1)\\ \mathbf a^{(-n)} = - \mathbf a^{(n)}|_{n=1,2,3} \\ \end{array} \right \rbrace } \]

Primary RLVs

\[ \mathcal R_{\textrm{sc}}^{(1)} = { \left \lbrace \begin{array}{l} \mathbf q^{(1)} = (1,0,0) \\ \mathbf q^{(2)} = (0,1,0) \\ \mathbf q^{(3)} = (0,0,1)\\ \mathbf q^{(-n)} = - \mathbf q^{(n)} |_{n=1,2,3}\\ \end{array} \right \rbrace } \]

\[ \mathcal R_{\textrm{sc}}^{(2)} = { \left \lbrace \begin{array}{l} \mathbf q^{(4)} = (0,1,1) \\ \mathbf q^{(5)} = (1,0,1) \\ \mathbf q^{(6)} = (1,1,0)\\ \mathbf q^{(7)} = (0,-1,1) \\ \mathbf q^{(8)} = (-1,0,1) \\ \mathbf q^{(9)} = (-1,1,0)\\ \mathbf q^{(-n)} = - \mathbf q^{(n)}|_{n=4,...,9} \\\\ \end{array} \right \rbrace } \]

\[\mathcal R_{\textrm{sc}}^{(3)} = { \left \lbrace \begin{array}{l} \mathbf q^{(10)} = (-1,1,1) \\ \mathbf q^{(11)} = (1,-1,1) \\ \mathbf q^{(12)} = (1,1,-1)\\ \mathbf q^{(13)} = (1,1,1)\\ \mathbf q^{(-n)} = - \mathbf q^{(n)}|_{n=10,...,13} \\ \end{array} \right \rbrace } \]

The ground state

To find the ground state of a PFC, one inserts a particular mode approximation into the free energy density, average over a unit cell (coarse-grain) and minimizes it with respect to the amplitudes.

Example: The 1D PFC - Ground state

To find the ground state of the 1D PFC, we assume that we can express the field in the one-mode approximation, i.e.,

\[ \psi^{\textrm{eq}}(x) = \psi_0 + A (e^{\mathfrak i q x} + e^{\mathfrak i q x}), \]

where \(q\) is, for now, an arbitrary wave vector which we will show to be the primary RLV of the lattice. We want to insert this into the free energy density \(\tilde f\) and integrate over a unit cell (UC), which in this case is simply \([0,2\pi]\)

\[ \frac{1}{2\pi} \mathcal F_{UC} = \frac{1}{2\pi} \int_0^{2\pi} \tilde f(\psi^{\textrm{eq}}) dx \equiv \langle \tilde f(\psi^{\textrm{eq}}) \rangle_{UC} = \frac{1}{2} \langle (\mathcal L {\psi^{\textrm{eq}}})^2 \rangle_{UC} + \frac{1}{2} \texttt{r} \langle {\psi^{\textrm{eq}}}^2 \rangle_{UC} + \frac{1}{4} \langle {\psi^{\textrm{eq}}}^4 \rangle_{UC} \]

To insert this into the free energy density, we need to calculate some auxiliary quantities.

\[ \mathcal L \psi^{\textrm{eq}} = (1+\nabla^2) \psi^{\textrm{eq}} = \psi_0 + (1-q^2) A (e^{\mathfrak i q x} + e^{\mathfrak i q x}) \]

so

\[ \langle (\mathcal L \psi^{\textrm{eq}})^2 \rangle_{UC} = \langle (\psi_0 + (1-q^2) A (e^{\mathfrak i q x} + e^{\mathfrak i q x})(\psi_0 + (1-q^2) A (e^{\mathfrak i q x} + e^{\mathfrak i q x})) \rangle_{UC} \]

Now, to calculate this, we need to make use of a condition of resonance, which is that under the coarse-graining operation, terms with periodicity average to zero². So we get

\[ \langle (\mathcal L \psi^{\textrm{eq}})^2 \rangle_{UC} = \psi_0^2 + 2 A^2 (1-q^2)^2 \]

Similarly, the other terms give

\[ \langle {\psi^{\textrm{eq}}}^2 \rangle_{UC} = \psi_0^2 + 2 A^2 \]

\[ \langle {\psi^{\textrm{eq}}}^3 \rangle_{UC} = \psi_0^3 + 6 \psi_0 A^2 \]

\[ \langle {\psi^{\textrm{eq}}}^4 \rangle_{UC} = \psi_0^4 + 12 \psi_0^2 A^2 + 6 A^4 \]

which gives

\[ \mathcal F_{UC} = 2\pi \left( \frac{1}{2} (\psi_0^2 + 2 A^2 (1-q^2)^2) + \frac{1}{2} \texttt{r} (\psi_0^2 + 2 A^2) + \frac{1}{3} \texttt t (\psi_0^3 + 6 \psi_0 A^2) + \frac{1}{4} \texttt v (\psi_0^4 + 12 \psi_0^2 A^2 + 6 A^4) \right ) \]

First of all, we see that the equilibrium value must have \(q=1\), since the free energy is minimized by this choice. Now, if the PFC is evolving according to conservative dynamics \(\partial_t \psi = \nabla^2 (\delta \mathcal F/\delta \psi)\) (which is the standard PFC evolution equation as explained further down in the document), the average value \(\psi_0\) will be conserved and therefore considered a simulation constant. In this case, we find the equilibrium value of \(A\), by differentiating wrt. \(A\) and setting to zero.

\[ \partial_A \mathcal F_{UC} = 2 \texttt{r} A + 4 \texttt t \psi_0 A + 6 \texttt v \psi_0^2 A + 6\texttt v A^3 = 0, \]

which gives \(A=0\) or

\[ \texttt{r} + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2 + 3 \texttt v A^2 = 0, \]

\[ A^2 = \frac{1}{3\texttt v} (-\texttt{r} - 2 \texttt t \psi_0 - 3\texttt v \psi_0^2) \]

\[ A = \pm \frac{1}{\sqrt{3\texttt v}} \sqrt{- \texttt{r} - 2 \texttt t \psi_0 - 3 \psi_0^2} \]

The sign does not matter in the case of the 1D PFC, since it will only affect the global sign of the amplitude, flipping the cosine wave. In higher dimensions, however, it requires more subtle treatment, and typically, one inserts both solutions into the free energy \(\mathcal F_{UC}\) to see which is lower

If instead, we are evolving the PFC under not conserved dynamics \(\partial_t \psi = - \delta \mathcal F/\delta \psi\), also the zero mode \(\psi_0\) will reach an equilibrium value. In this case, then, we need to solve the system of equations given by

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & \psi_0 + \texttt r \psi_0 + \texttt t \psi_0^2 + 2 \texttt t A^2 + \texttt v \psi_0^3 + 6\texttt v \psi_0 A^2 &=0 \\ \partial_{A} \mathcal F_{UC} = 0: & \texttt{r} + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2 + 3 \texttt v A^2 &= 0\\ \end{array} \right \rbrace \]

The pfc instance will be initialized with a list called eta0, which consists of the equilibrium values of the amplitudes to begin with.

1D periodic2D triangular2D square3D body-centered cubic3D face-centered cubic3D simple cubic

Default resolution:

Default model parameters \((r,\psi_0)\): \((-0.3,-0.3)\)

Free energy per unit cell (calculation document):

\[ \mathcal F_{UC} = \frac{1}{2} \psi_0^2 + \frac{1}{2} \texttt{r} (\psi_0^2 + 2 A^2) + \frac{1}{3} \texttt t (\psi_0^3 + 6 \psi_0 A^2) + \frac{1}{4} \texttt v (\psi_0^4 + 12 \psi_0^2 A^2 + 6 A^4) \]

Equilibrium amplitude (conserved)

\[ \partial_{A} \mathcal F_{UC} = 0: \quad A =\pm \frac{1}{\sqrt {3\texttt v}} \sqrt{ -\texttt r - 2 \texttt t \psi_0 - 3 \texttt v \psi_0^2} \]

Equilibrium amplitude equations (unconserved)

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 2 A^2 (\texttt t + 3 \texttt v \psi_0) + \psi_0 (1 + \texttt r + \texttt t \psi_0 + \texttt v \psi_0^2) &=0 \\ \partial_{A} \mathcal F_{UC} = 0: & \texttt r + 3 A^2 \texttt v + \psi_0 (2 \texttt t + 3 \texttt v \psi_0) &= 0\\ \end{array} \right \rbrace \]

Default resolution: \([7,12]^{-1}a_0\)

Default model parameters \((r,\psi_0)\): \((-0.3,-0.3)\)

Free energy per unit cell (calculation document.):

\[ \mathcal F_{UC} = \frac{\pi^2}{3\sqrt 3} \left (270 A^4 \texttt v + 48 A^3 (\texttt t + 3 \texttt v \psi_0) + \psi_0^2 (6 + 6 \texttt r + 4 \texttt t \psi_0 + 3 \texttt v \psi_0^2) + 36 A^2 (\texttt r + \psi_0 (2 \texttt t + 3 \texttt v \psi_0))\right) \]

Equilibrium amplitude equations (conserved)

\[ \left \lbrace \partial_{A} \mathcal F_{UC} = 0: A=0 ~ \textrm{or} ~ \quad A= \frac{1}{15\texttt v} \left ( -\texttt t - 3 \texttt v \psi_0 \pm \sqrt{ \texttt t^2 - 15 \texttt r \texttt v - 24 \texttt t \texttt v \psi_0 - 36 \texttt v^2 \psi_0^2} \right ) \right \rbrace \]

Equilibrium amplitude equations (unconserved)

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 12 A^3 \texttt v + 6 A^2 (\texttt t + 3 \texttt v \psi_0) + \psi_0 (1 + \texttt r + \texttt t \psi_0 + \texttt v \psi_0^2) &=0 \\ \partial_{A} \mathcal F_{UC} = 0: & A(\texttt r + 15 A^2 \texttt v + 2 A (\texttt t + 3 \texttt v \psi_0) + \psi_0 (2 \texttt t + 3 \texttt v \psi_0)) &= 0\\ \end{array} \right \rbrace \]

Default resolution: \([7,7]^{-1}a_0\)

Default model parameters \((\texttt r,\psi_0)\): \((-0.3,-0.3)\)

Free energy per unit cell (calculation document):

\[ \mathcal F_{UC} = \frac{\pi^2}{3} (108 A^4 \texttt v + 108 B^4 \texttt v + \psi_0^2 (24 + 6 \texttt r + 4 \texttt t \psi_0 + 3 \texttt v \psi_0^2) + 24 B^2 (\texttt r + \psi_0 (2 \texttt t + 3 \texttt v \psi_0)) + 24 A^2 (\texttt r + 18 B^2 \texttt v + 4 B (t + 3 \texttt v \psi_0) + \psi_0 (2 \texttt t + 3 \texttt v \psi_0))) \]

Equilibrium amplitude equations

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 4 B^2 (\texttt t + 3 \texttt v \psi_0) + \psi_0 (4 + \texttt r + \texttt t \psi_0 + \texttt v \psi_0^2) + 4 A^2 (\texttt t + 3 \texttt v (2 B + \psi_0)) &= 0 \\ \partial_{A} \mathcal F_{UC} = 0: & A (\texttt r + 9 A^2 \texttt v + 18 B^2 \texttt v + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2 + 4 B (\texttt t + 3 \texttt v \psi_0)) &= 0 \\ \partial_{B} \mathcal F_{UC} = 0: & 9 B^3 \texttt v + 2 A^2 (t + 3 \texttt v \psi_0) + B (\texttt r + 18 A^2 \texttt v + 2 \texttt \texttt t \psi_0 + 3 \texttt v \psi_0^2) &= 0\\ \end{array} \right \rbrace \]

Default resolution: \([7,7,7]^{-1}a_0\)

Default model parameters \((\texttt r,\psi_0)\): \((-0.3,-0.325)\)

Free energy per unit cell (calculation document):

\[ \mathcal F_{UC} = \frac{4 \sqrt 2 \pi^3}{3} \left (1620 A^4 \texttt v + 192 A^3 (t + 3 \texttt v \psi_0) + \psi_0^2 (6 + 6 \texttt r + 4 t \psi_0 + 3 \texttt v \psi_0^2) + 72 A^2 (\texttt r + \psi_0 (2 t + 3 \texttt v \psi_0)) \right ) \]

Equilibrium amplitude (conserved)

\[ \left \lbrace \partial_{A} \mathcal F_{UC} = 0: A=0 ~ \textrm{or} ~ \quad A = \frac{1}{45\texttt v} \left ( -2 \texttt t - 6 \texttt v \psi_0 \pm \sqrt{ 4 \texttt t^2 - 45 \texttt r \texttt v - 66 \texttt t \texttt v \psi_0 - 99 \texttt v^2 \psi_0^2} \right) \right \rbrace \]

Equilibrium amplitude equations (unconserved)

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 48 A^3 \texttt v + 12 A^2 (\texttt t + 3 \texttt v \psi_0) + \psi_0 (1 + \texttt r + \texttt t \psi_0 + \texttt v \psi_0^2) &=0 \\ \partial_{A} \mathcal F_{UC} = 0: & A (\texttt r + 45 A^2 \texttt v + 4 A (\texttt t + 3 \texttt v \psi_0) + \psi_0 (2 \texttt t + 3 \texttt v \psi_0)) &= 0\\ \end{array} \right \rbrace \]

Default resolution: \([11,11,11]^{-1}a_0\)

Default model parameters \((\texttt r,\psi_0)\): \((-0.3,-0.325)\)

Free energy per unit cell (calculation document):

\[ \mathcal F_{UC} = \frac{2 \pi^3}{\sqrt 3} (1944 A^4 \texttt v + 810 B^4 \texttt v + \psi_0^2 (32 + 18 \texttt r + 12 \texttt t \psi_0 + 9 \texttt v \psi_0^2) + 108 B^2 (r + \psi_0 (2 \texttt t + 3 \texttt v \psi_0)) + 144 A^2 (\texttt r + 36 B^2 \texttt v + 6 B (t + 3 \texttt v \psi_0) + \psi_0 (2 \texttt t + 3 \texttt v \psi_0))) \]

Equilibrium amplitude equations

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 54 B^2 (t + 3 \texttt v \psi_0) + \psi_0 (16 + 9 \texttt r + 9 t \psi_0 + 9 \texttt v \psi_0^2) + 72 A^2 (t + 3 \texttt v (3 B + \psi_0)) &= 0 \\ \partial_{A} \mathcal F_{UC} = 0: & A (\texttt r + 27 A^2 \texttt v + 36 B^2 \texttt v + 2 t \psi_0 + 3 \texttt v \psi_0^2 + 6 B (t + 3 \texttt v \psi_0)) &= 0 \\ \partial_{B} \mathcal F_{UC} = 0: & 15 B^3 \texttt v + 4 A^2 (t + 3 \texttt v \psi_0) + B (\texttt r + 48 A^2 \texttt v + 2 t \psi_0 + 3 \texttt v \psi_0^2) &= 0\\ \end{array} \right \rbrace \]

Default resolution: \([5,5,5]^{-1}a_0\)

Default model parameters \((\texttt r,\psi_0)\): \((-0.3,-0.325)\)

Free energy per unit cell (calculation document):

\[ \mathcal F_{UC} = \frac{2\pi^3}{3} \left (48 C^2 \texttt r + 270 A^4 \texttt v + 1620 B^4 \texttt v + 576 A^3 C \texttt v + 648 C^4 \texttt v + 96 C^2 \texttt t \psi_0 + 6 (36 + \texttt r + 24 C^2 \texttt v) \psi_0^2 + 4 \texttt t \psi_0^3 + 3 \texttt v \psi_0^4 + 192 B^3 (\texttt t + 3 \texttt v \psi_0) + 576 A B C (t + 3 \texttt v (3 B + \psi_0)) + 36 A^2 (\texttt r + 96 B^2 \texttt v + 36 C^2 \texttt v + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2 + 8 B (\texttt t + 3 \texttt v \psi_0)) + 72 B^2 (r + 54 C^2 \texttt v + \psi_0 (2 \texttt t + 3 \texttt v \psi_0)) \right ) \]

Equilibrium amplitude equations:

\[ \left \lbrace \begin{array}{lrl} \partial_{\psi_0} \mathcal F_{UC} = 0: & 8 C^2 \texttt t + 48 B^3 \texttt v + 144 A B C \texttt v + 36 \psi_0 + \texttt r \psi_0 + 24 C^2 \texttt v \psi_0 + \texttt t \psi_0^2 + \texttt v \psi_0^3 + 12 B^2 (t + 3 \texttt v \psi_0) + 6 A^2 (t + 3 \texttt v (4 B + \psi_0)) &= 0 \\ \partial_{A} \mathcal F_{UC} = 0: & 15 A^3 \texttt v + 24 A^2 C \texttt v + 8 B C (t + 3 \texttt v (3 B + \psi_0)) + A (\texttt r + 96 B^2 \texttt v + 36 C^2 \texttt v + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2 + 8 B (t + 3 \texttt v \psi_0)) &= 0 \\ \partial_{B} \mathcal F_{UC} = 0: & 45 B^3 \texttt v + 4 B^2 (t + 3 \texttt v \psi_0) + 2 A (A + 2 C) (t + 3 \texttt v \psi_0) + B (\texttt r + 48 A^2 \texttt v + 72 A C \texttt v + 54 C^2 \texttt v + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2) &= 0\\ \partial_{C} \mathcal F_{UC} = 0: & 27 C^3 \texttt v + C (\texttt r + 27 A^2 \texttt v + 81 B^2 \texttt v + 2 \texttt t \psi_0 + 3 \texttt v \psi_0^2) + 6 A (A^2 \texttt v + 9 B^2 \texttt v + B (t + 3 \texttt v \psi_0)) &= 0\\ \end{array} \right \rbrace \]

We refer to these amplitudes \(A,B,C\) as proto amplitudes and they are calculated by the functions calc_proto_amplitudes_conserved (by default) and calc_proto_amplitudes_unconserved in the pfc class. The proto amplitudes are saved in the pfc instance as eta0. The function calc_free_energy_from_proto_amplitudes calculates the free energy \(\mathcal F_{UC}\) from the proto amplitudes, including \(\psi_0\).

Straining the ground state to equilibrium

While the q-vector given above and the proto amplitudes are good approximations to the ground state, they are not exact. In fact, if the PFC is initialized with these values, it will experience a slight residual stress ¹⁰. To account for this, when initializing a PFC, the code will automatically strain the field to find the optimal value of \(q_0\), which is typically a few percent off the primary RLV. This is done numerically by simulating a \(1\times 1\) unit cell and using the function self.conf_strain_to_equilibrium in the pfc class.

In this function, we apply the following distortion

\[ \mathcal u_{ij} = \epsilon \delta_{ij} \]

Given the equilibrium configuration, we can find more accurate estimates for the equilibrium amplitudes, but evolving the PFC for a few time steps and calculating the amplitudes from the demodulation of the field with the primary RLVs. This is done in the calc_strained_amplitudes function in the pfc class, which also calculates new values for the elastic constants (see below).

The amplitude approximation - deviations from the ground state

So far, we have only looked at the equilibrium state of the PFC, which has a specific symmetry, and the few-mode approximations that can be used to describe the PFC in this state. For small deviations from the equilibrium state, like in the presence of few dislocations and small strains, we can use the amplitude approximation, in which we assume that the field can be written as

The amplitude approximation

\[ \psi \approx \bar \psi + \sum_{\mathbf q^{(n)} \in \mathcal R^{(1)}} \eta_n(\mathbf r) e^{\mathfrak i \mathbf q^{(n)} \cdot \mathbf r}, \]

where \(\eta_n\) are slowly varying complex fields. These fields can be found by demodulating the field \(\psi\) with the primary RLVs, i.e.,

Demodulation

\[ \eta_n = \langle \psi e^{- \mathfrak i \boldsymbol{q}^{(n)} \cdot \boldsymbol{r}} \rangle \]

The figure below shows an example of the evolution of the PFC for a triangular lattice and the corresponding amplitude fields.

Figure: Snapshots of an example evolution of the 2D triangular PFC model at (top row) \(t=0\) and (bottom row) \(t=600\). Parameters used were \((r,\psi_0) = (-0.3,-0.3)\). The columns show the demodulated fields \(\bar \psi\) and \(\{ \eta_n \}_{n=1}^{3}\), where the complex fields are shown by their phase \(\theta_n\) and brightness corresponding to the magnitude \(|\eta_n|\). Taken from Ref. ² with permission.

Elasticity

This quantity is the Poisson ratio, under some conditions. I think. To be fixed. \(\(\nu = \frac{\lambda}{(d-1)\lambda + 2\mu + \gamma}\)\)

Elastic constants are saved in these quantities

pfc.el_mu
pfc.el_lambda
pfc.el_gamma
pfc.el_nu

Stress tensor

The equations for calculating the stress tensor are found in Ref. ⁶, but we list them below, together wtih the elastic constants for each of the models. The stress tensor \(h_{ij}\) and its associated elastic constants interms of amplitudes (\(A,B,C\)) of the mode expansion for different PFCmodels. Here, \(\mathcal L_X = X+\nabla^2\). The elastic constants canbe expressed in Voigt notation by \(C_{11} = \lambda+2\mu+\gamma\),\(C_{12} = \lambda\), \(C_{44}=\mu\).

Example: The 1D PFC - Stress tensor

In order to derive the stress tensor, suppose we have a state of the PFC given by \(\psi\). We may deform this state in a mass-conserving fashion by a displacement field \(u\) (which is a scalar field in one dimension) as follows

\[ \psi' = \psi - \partial_x (\psi u) \]

We insert this into the expression for the free energy giving a new free energy \(\mathcal F'\).

\[ \mathcal F' = \int dx \left ( \frac{1}{2} ((1+\partial_x^2) \psi'))^2 + \frac{1}{2} \texttt r \psi'^2 + \frac{1}{3} \texttt t \psi'^3 + \frac{1}{4} \texttt v \psi'^4 \right ) \]

If we take \(u\) to be an infinitesimal exploratory field, we can expand in \(u\), which gives

\[ \mathcal F' = \mathcal F[\psi] - \int dx (1+\partial_x^2)\psi (1+\partial_x^2) \partial_x (\psi u_x) - \texttt r \psi \partial_x (\psi u_x) - \texttt t \psi^2 \partial_x (\psi u_x) - \texttt v \psi^3 \partial_x (\psi u_x) \]

\[ = \mathcal F[\psi] - \int dx (1+\partial_x^2)\psi (1+\partial_x^2) \partial_x (\psi u_x) + (\tilde \mu_c - \mathcal L^2) \partial_x (\psi u_x) \]

\[ \nabla \cdot h = \left \langle \tilde \mu_c \nabla \psi - \nabla \tilde f \right \rangle \]

1D periodic2D triangular2D square3D body-centered cubic3D face-centered cubic3D simple cubic

Stress tensor

Elastic constants

\[ \lambda = \]

\[ \mu = \]

\[ \gamma = \]

Stress tensor

\[ h_{ij} = -2\left \langle (\mathcal L_1 \psi) \partial_{ij} \psi \right \rangle \]

Elastic constants

\[\lambda = 3A^2\]

\[\mu = 3A^2\]

\[\gamma = 0\]

Stress tensor $$ h_{ij} = -2\left \langle (\mathcal L_1 \mathcal L_2 \psi)(\mathcal L_1 + \mathcal L_2) \partial_{ij} \psi \right \rangle$$

Elastic constants

\[\lambda = 16B^2\]

\[\mu = 16B^2\]

\[\gamma = 8A^2-32B^2\]

Stress tensor

\[h_{ij} = -2\left \langle (\mathcal L_1 \psi) \partial_{ij} \psi \right \rangle\]

Elastic constants

\[\lambda = 4 A^2\]

\[\mu = 4 A^2\]

\[\gamma = -4A^2\]

Stress tensor

\[h_{ij} = -2\left \langle (\mathcal L_1 \mathcal L_{\frac 4 3} \psi)(\mathcal L_1 + \mathcal L_{\frac 4 3}) \partial_{ij} \psi \right \rangle\]

Elastic constants

\[\lambda = \frac{32}{81} A^2\]

\[\mu = \frac{32}{81} A^2\]

\[\gamma = \frac{32}{81} (2B^2 - A^2)\]

Stress tensor

\[h_{ij} = -2\left \langle (\mathcal L_1 \mathcal L_2 \mathcal L_3 \psi)(\mathcal L_2 \mathcal L_3 + \mathcal L_1 \mathcal L_3 + \mathcal L_1 \mathcal L_2) \partial_{ij} \psi \right \rangle\]

Elastic constants

\[\lambda = 16 B^2 + 128 C^2\]

\[\mu = 16 B^2 + 128 C^2\]

\[\gamma = 32 A^2 - 16 B^2 - 256 C^2\]

As one sees, the expressions for the stress tensor all on the format

\[ h_{ij} = -2\left \langle (\mathcal L \psi)(\mathcal L_{\sum} \partial_{ij} \psi \right \rangle, \]

where \(\mathcal L\) is as defined previously (as the product \(\mathcal L_1 ...\)), and \(\mathcal L_{\sum}\) is a particular sum of \(\mathcal L_i\) operators. \(\mathcal L\) is calculated (in Fourier space) by the function calc_L_f, and the \(\mathcal L_{\sum}\) is calculated by the function calc_L_sum_f in the pfc class. The stress tensor is then calculated by the function calc_microscopic_stress_tensor in the pfc class.

Equilibrium elastic constants

As mentioned above, the equilbrium state of the PFC is not exactly given by \(q_0=1\), and the proto amplitudes are not exact. However, in the process of deriving the equilibrium amplitudes after straining, we can also use the \(1\times 1\) unit cell to calculate the elastic constants. This is done numerically by deforming the unit cell and calculating the increase in the elastic energy. To illustrate this, consider a distortion of pure shear, i.e.,

\[ \mathfrak u = \begin{pmatrix} 0 & \epsilon & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \]

The corresponding strain is given by

\[ \varepsilon = \begin{pmatrix} 0 & \frac 1 2 \epsilon & 0 \\ \frac 1 2 \epsilon & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \]

Under this deformation, the elastic energy

\[ F_{el} = \frac{1}{2} \mathcal C_{ijkl} \varepsilon_{ij} \varepsilon_{kl} \]

is given by

\[ F_{el} = \mu e_{ij} e_{ij} = \frac{1}{2} \mu \epsilon^2. \]

We can apply the distortion \(\mathfrak u\) with the function conf_apply_distortion in the pfc class. The result of this procedure for a \(4\times 4\) square PFC is shown in the figure below

Figure: The elastic constants of a 2D square PFC model as a function of the strain \(\epsilon\) in the \(x\)-direction. The fit shows the numerical fit of the elastic constant \(\mu\) to the applied strain, showing excellent agreement.

While this works for the elastic constant \(\mu\), the full elastic free energy is determined by three elastic constants, \(\lambda\), \(\mu\), and \(\gamma\)

\[ F_{el} = F_{el} = \frac{1}{2} \lambda (\varepsilon_{kk})^2 + \mu \varepsilon_{ij} \varepsilon_{ij} + \frac{1}{2} \gamma (\varepsilon_{11}^2 + \varepsilon_{22}^2 + \varepsilon_{33}^2). \]

To find these, we apply a series of deformations to the unit cell and use scipy.optimize.curve_fit to fit the elastic energy to the applied strain. This is done in the calc_strained_amplitudes and set to the pfc instance as el_lambda, el_mu, and el_gamma upon initialization. The values are also printed to the terminal at initialization.

The expression for the elastic constants in terms of the proto amplitudes were derived by inserting the amplitude approximation into the expression for the stress tensor. It is interesting to see that elastic constants obtained by using the equilibrium values of amplitudes match the elastic constants obtained by straining the PFC to equilibrium and performing the fit.

Strains

In order to calculate the strain of a PFC configuration, we follow a generalized procedure as that outlined in Ref. ². The ground state of the PFC can be written as

\[ \psi = \psi_0 + \eta_1 \sum_{\mathbf q^{(n)} \in \mathcal R^{(1)}} e^{i \mathbf q^{(n)} \cdot \mathbf r} + \eta_2 \sum_{\mathbf q^{(n)} \in \mathcal R^{(2)}} e^{i \mathbf q^{(n)} \cdot \mathbf r} +..., \]

where \(\mathcal R^{(n)}\) is the nth closest modes, i.e., the full reciprocal lattice can be written as

\[ \mathcal R = \{\mathbf 0 \} \cup \left (\cup_{n=1}^\infty \mathcal R^{(n)} \right) \]

We define the following quantity

\[ \Phi = N_1q_1^2\eta_1^2 + N_2q_2^2\eta_2^2 + ... = \sum_{\mathbf q^{(n)} \in \mathcal R} \eta_{n,0}^2|\mathbf q^{(n)}|^2 \]

Assuming that the mode-approximations hold, we can calculate \(\Phi\) from the equilibrium amplitudes which are stored in the PhaseFieldCrystal instance, which is done routinely in the initialization of the class and stored in pfc.Phi. Consider now the structure function of the equilibrium state

\[ \mathcal S_{ij} = \langle \partial_i \psi \partial_j \psi\rangle \]

\[ = \eta_1^2 \sum_{\mathbf q^{(n)} \in \mathcal R^{(1)}} q^{(n)}_iq^{(n)}_j + \eta_2^2 \sum_{\mathbf q^{(n)} \in \mathcal R^{(2)}} q_i^{(n)}q_j^{(n)}+... \]

\[ = \eta_1^2 \frac{N_1 q_1^2}{d} \delta_{ij} + \eta_2^2 \frac{N_2 q_2^2}{d} \delta_{ij}+... \]

\[ = \frac{\Phi}{d}\delta_{ij} \]

The strain of the phase-field crystal can be calculated using the following formula The strain of the phase-field crystal can be calculated using the following formula²

\[ \varepsilon_{ij} = \frac{1}{2}\delta_{ij} - \frac{d}{2\Phi} \mathcal S_{ij}, \]

where \(\mathcal S = \langle \nabla \psi \nabla \psi \rangle\) is called the structure function ⁷ and is implemented in the calc_strain_tensor method of the PhaseFieldCrystal class. It should be noted that there is some residual strain due to the equilibrium periodicity of the lattice not matching exactly the periodicity of the simulation domain ¹⁰. It is therefore often useful to subtract the mean value for visualization purposes.

strain = pfc.calc_strain_tensor()
strain = strain - np.mean(strain, axis=0)

Dislocations

The Burgers vector is defined by

Burgers vector definition

\[ \oint_{\partial \mathcal M} d \boldsymbol{u} = -\boldsymbol{b}. \]

The minus sign in this convention reflects the fact that we consider the Burgers vector to be the disconnection error from the ending point to the starting point, when going an oriented path around the dislocation,

Burgers vector definition: (a) The one-body density of a crystalline solid containing an edge dislocation in a 2D square lattice (superimposed), (b) a 3D simple cubic lattice with an edge dislocation (\(\mathbf b \perp \mathbf t\)), and (c) a 3D simple cubic lattice with a screw dislocation (\(\mathbf b \parallel \mathbf t)\). In all cases, a circulation (green) that is right-handed with respect to the tangent vector \(\mathbf t\), i.e., following a path around the dislocation, gives rise to a connection error: the Burgers vector \(\mathbf b\).
Reprinted from Ref. ² with permission.

In three dimensions, the path defining the dislocation is given by the direction of the dislocation line tangent \(\boldsymbol{t}\). Multiplying this equation by the reciprocal lattice vector \(-\boldsymbol{q}^{(n)}\), we get

\[\oint_{\partial \mathcal M} d (-\boldsymbol{q}^{(n)} \cdot \boldsymbol{u}) = (\boldsymbol{b} \cdot \boldsymbol{q}^{(n)}) \equiv 2\pi s_n,\]

where \(s_n\) is the charge associated with the Burgers vector \(\boldsymbol{b}\)

Dislocation charge

\[ s_n = \frac{1}{2\pi} \boldsymbol{b} \cdot \boldsymbol{q}^{(n)}. \]

Using the primary BLVs we get the charges summarized below

2D triangular

Burgers vector \(\mathbf b\)	\(s_1\)	\(s_2\)	\(s_3\)
\(\mathbf a^{(1)} = a_0 (1,0)\)	\(\color{red} 1\)	\(0\)	\(\color{blue}-1\)
\(\mathbf a^{(2)} = a_0 (1/2,\sqrt 3/2)\)	\(0\)	\(\color{red} 1\)	\(\color{blue}-1\)
\(\mathbf a^{(3)} = a_0 (1/2,-\sqrt 3/2)\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(0\)

2D square

Burgers vector \(\mathbf b\)	\(s_1\)	\(s_2\)		\(s_3\)	\(s_4\)
\(\mathbf a^{(1)} = a_0 (1,0)\)	\(\color{red} 1\)	\(0\)		\(\color{red} 1\)	\(\color{red} 1\)
\(\mathbf a^{(2)} = a_0(0,1)\)	\(0\)	\(\color{red} 1\)		\(\color{blue}-1\)	\(\color{red} 1\)

3D body-centered cubic

Burgers vector \(\mathbf b\)	\(s_1\)	\(s_2\)	\(s_3\)	\(s_4\)	\(s_5\)	\(s_6\)
\(\mathbf a^{(1)} = a_0/2 (-1,1,1)\)	\(\color{red} 1\)	\(0\)	\(0\)	\(0\)	\(\color{red} 1\)	\(\color{red} 1\)
\(\mathbf a^{(2)} = a_0/2 (1,-1,1)\)	\(0\)	\(\color{red} 1\)	\(0\)	\(\color{red} 1\)	\(0\)	\(\color{blue}-1\)
\(\mathbf a^{(3)} = a_0/2 (1,1,-1)\)	\(0\)	\(0\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(\color{blue}-1\)	\(0\)
\(\mathbf a^{(4)} = a_0/2 (1,1,1)\)	\(\color{red} 1\)	\(\color{red} 1\)	\(\color{red} 1\)	\(0\)	\(0\)	\(0\)

3D face-centered cubic

Burgers vector \(\mathbf b\)	\(s_1\)	\(s_2\)	\(s_3\)	\(s_4\)		\(s_5\)	\(s_6\)	\(s_7\)
\(\mathbf a^{(1)} = a_0/2(0,1,1)\)	\(\color{red} 1\)	\(0\)	\(0\)	\(\color{red} 1\)		\(0\)	\(\color{red} 1\)	\(\color{red} 1\)
\(\mathbf a^{(2)} = a_0/2(1,0,1)\)	\(0\)	\(\color{red} 1\)	\(0\)	\(\color{red} 1\)		\(\color{red} 1\)	\(0\)	\(\color{red} 1\)
\(\mathbf a^{(3)} = a_0/2 (1,1,0)\)	\(0\)	\(0\)	\(\color{red} 1\)	\(\color{red} 1\)		\(\color{red} 1\)	\(\color{red} 1\)	\(0\)
\(\mathbf a^{(4)} = a_0/2(0,-1,1)\)	\(0\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(0\)		\(0\)	\(\color{blue}-1\)	\(\color{red} 1\)
\(\mathbf a^{(5)} = a_0/2 (-1,0,1)\)	\(\color{red} 1\)	\(0\)	\(\color{blue}-1\)	\(0\)		\(\color{blue}-1\)	\(0\)	\(\color{red} 1\)
\(\mathbf a^{(6)} = a_0/2(-1,1,0)\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(0\)	\(0\)		\(\color{blue}-1\)	\(\color{red} 1\)	\(0\)

3D simple cubic

Burgers vector \(\mathbf b\)	\(s_1\)	\(s_2\)	\(s_3\)		\(s_4\)	\(s_5\)	\(s_6\)	\(s_7\)	\(s_8\)	\(s_9\)		\(s_{10}\)	\(s_{11}\)	\(s_{12}\)	\(s_{13}\)
\(\mathbf a^{(1)} = a_0(1,0,0)\)	\(\color{red} 1\)	\(0\)	\(0\)		\(0\)	\(\color{red} 1\)	\(\color{red} 1\)	\(0\)	\(\color{blue}-1\)	\(\color{blue}-1\)		\(\color{blue}-1\)	\(\color{red} 1\)	\(\color{red} 1\)	\(\color{red} 1\)
\(\mathbf a^{(2)} = a_0 (0,1,0)\)	\(0\)	\(\color{red} 1\)	\(0\)		\(\color{red} 1\)	\(0\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(0\)	\(\color{red} 1\)		\(\color{red} 1\)	\(\color{blue}-1\)	\(\color{red} 1\)	\(\color{red} 1\)
\(\mathbf a^{(3)} = a_0(0,0,1)\)	\(0\)	\(0\)	\(\color{red} 1\)		\(\color{red} 1\)	\(\color{red} 1\)	\(0\)	\(\color{red} 1\)	\(\color{red} 1\)	\(0\)		\(\color{red} 1\)	\(\color{red} 1\)	\(\color{blue}-1\)	\(\color{red} 1\)

The dislocation density tensor

Given a PFC configuration, the dislocation density tensor may be calculated as ⁸

The dislocation density tensor

\[ \alpha_{ij} = \frac{2d}{N\eta_0^2} \sum_{n=1}^N D_i^{(n)} q_j^{(n)} = \frac{2 \pi d}{N} \sum_{n=1}^N \rho_i^{(n)} q_j^{(n)} \]

where \(\boldsymbol{D}^{(n)}\) is calculated from \(\boldsymbol{\psi }= (\Re(\eta_n), \Im(\eta_n))\) from the \(N\) primary reciprocal lattice vectors \(\boldsymbol{q}^{(n)}\) and

\[ \rho^{(n)}_i = \frac{1}{\pi \eta_0^2} D^{(n)}_i \]

In the language of this

Equations of motion

Conserved evolution

The basic equation motion presented in Ref. ¹ can be derived by postulating a simple mechanism for free energy minimzation under the constraint of mass conservation.

The PFC evolution equation (evolve_PFC)

\[ \partial_t \psi = \nabla^2 \tilde \mu_c, \]

The linear part \(\mathcal L^2 + \texttt r\) of the chemical potential (not multiplied by psi) is calculated by

pfc.calc_chemical_potential_linear_part_f()

The PFC chemical potential

which gives

\[ {{\omega }_{\mathfrak f}}= -\boldsymbol{k}^2 (\texttt r + {{\mathcal L}_{\mathfrak f}}^2) \]

\[ N = \nabla^2 (\psi^3) \]

Unconserved dynamics

The unconserved PFC evolution equation (evolve_PFC_unconserved)

\[ \partial_t \psi = - \tilde \mu_c \]

Evolution at mechanical equilibrium

The PFC evolution at mechanical equilibrium (evolve_PFC_mechanical_equilibrium)

Step 1:

\[\psi(t+\Delta t) = \textrm{Integrate($\Delta t$):} (\partial_t \psi)\]

Step 2:

\[\psi(t + \Delta t, \boldsymbol{r}) \leftarrow \psi(t + \Delta t, \boldsymbol{r} - \boldsymbol{u}^\delta),\]

where \(\boldsymbol{u}^\delta\) is the solution to

\[ \partial_{j} \mathcal C_{ijkl} \partial_l u_k = - g^\psi_i, \]

where

\[ g^\psi_i = -\partial_j \mathfrak h_{ij}^\psi, \]

and the elastic constants tensor \(\mathcal C\) is given by

\[ \mathcal C_{ijkl} = \lambda \delta_{ij}\delta_{kl} + 2\mu \delta_{k(i} \delta_{j)l} + \gamma \delta_{ijkl}. \]

This equation is solved in Fourier space by

\[ {u}_{\mathfrak f ~ i}^\delta = G_{\mathfrak f ~ ij} g_{\mathfrak f ~ j}^\psi, \]

where the Greens function \(G_{\mathfrak f ~ ij}\) is given in Ref. ⁹ as

\[ {G_{\mathfrak f ~ij}} (\mathbf k) =\frac{1}{\mathbf k^2}\left ( \frac{\delta_{ij}}{\mu + \gamma \kappa_{(i)}^2} - \frac{\kappa_i \kappa_j}{(\mu + \gamma \kappa_{(i)}^2 )(\mu + \gamma \kappa_{(j)}^2 )} \frac{\mu+\lambda}{1+ \sum_{l=1}^3 \frac{\mu + \lambda}{\mu+\gamma \kappa_l^2} \kappa_l^2 }\right ), \]

with \(\boldsymbol \kappa = \mathbf k/|\mathbf k|\) there is no implicit summation over indices \((i),(j)\).
By defining \(k_3=0\), this equation is also valid for the triangular and square symmetry in two dimensions.

Note that due to the asymmetry of the elastic constants, this method can only be used for small deviations of the lattice orientation.

Hydrodynamic PFC evolution

In Ref. ⁷, a hydrodynamic approach was derived. A simplified two-parameter model was proposed

The hydrodynamic PFC model (evolve_PFC_hydrodynamic)

\[ \partial_t \psi = \nabla^2 \tilde \mu_c - \boldsymbol{v} \cdot \nabla \psi \]

\[ \partial_t \boldsymbol{v} = \frac{1}{\rho_0} (\nabla \cdot h + \Gamma_S \nabla^2 \boldsymbol{v} + \boldsymbol{f}^{\textrm(ext)}) \]

We insert for \(\tilde \mu_c\) and write it in matrix form to emphasize the linear and non-linear parts

\[ \partial_t \begin{pmatrix} \psi \\ v_1 \\ v_2 \\ v_3 \end{pmatrix} = \begin{pmatrix} \nabla^2 (r + \mathcal L^2) \psi \\ \frac{1}{\rho_0}\Gamma_S \nabla^2 v_1 \\ \frac{1}{\rho_0}\Gamma_S \nabla^2 v_2 \\ \frac{1}{\rho_0}\Gamma_S \nabla^2 v_3 \end{pmatrix} + \begin{pmatrix} \nabla^2 ( \psi^3) - \boldsymbol{v} \cdot \nabla \psi \\ \frac{1}{\rho_0} \left ( \left \langle \tilde \mu_c \partial_x \psi - \partial_x \tilde f \right \rangle + f_x^{(ext)}\right )\\ \frac{1}{\rho_0} \left ( \left \langle \tilde \mu_c \partial_y \psi - \partial_y \tilde f \right \rangle + f_y^{(ext)}\right ) \\ \frac{1}{\rho_0} \left ( \left \langle \tilde \mu_c \partial_z \psi - \partial_z \tilde f \right \rangle + f_z^{(ext)} \right ) \end{pmatrix} \]

Configurations

The PFC class has a number of methods for generating initial conditions.

Dislocation dipoles and loops

The simplest setup for the PFC is to insert a dislocation dipole or dislocation loop.

Polycrystal configurations

To set a polycrystal configuration, one typically calculates the PFC from a set of rotated reciprocal lattices. This is done in the calc_PFC_from_amplitudes method, using the rotation keyword argument. Then, one constructs a field with the desired orientation by setting the field in the regions corresponding to the rotated reciprocal lattice to the rotated field. In order to avoid numerical artifacts with the interface, it is recommended to smooth the interface by evolving the PFC for a few time steps according to classical PFC dynamics.

Example: Creating a circular inclusion

To create a circular inclusion in 2 dimensions, you can run the following code.

import comfit as cf
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp

pfc = cf.PhaseFieldCrystal2DSquare(30,30)
pfc.dt=0.05

# This creates a standard orientation of the crystal
pfc.conf_PFC_from_amplitudes()

# Create the rotated field
psi_rotated = pfc.calc_PFC_from_amplitudes(rotation=np.pi/6)

# Specify the region centered at the mid position with radius 6 a0.
inclusion_region = pfc.calc_region_disk(pfc.rmid, 6*pfc.a0)

# Set the rotated field in the inclusion region
pfc.psi[inclusion_region] = psi_rotated[inclusion_region]
pfc.psi_f = sp.fft.fftn(pfc.psi)

#Smooth the interface
tau = 10
pfc.evolve_PFC(round(tau/pfc.dt)) 

pfc.plot_field(pfc.psi)
plt.show()

This will produce the following image.

The class comes with a variety of pre-configured setups for polycrystalline configurations, which are created using the conf_create_polycrystal method. This function takes as argument a keyword argument type, which is a string specifying the type of polycrystal to create. The available types are

circularfour_grain

This creates a circular inclusion valid in 2 and 3 dimensions. The type opens for two additional keywords:

• radius which specifies the radius of the inclusion (default self.size_min/4).

• position which specifies the position of the center of the inclusion (default self.rmid).

• rotation a float specifying the rotation of the inclusion wrt. the z-axis (default np.pi/6).

If initialized in 3 dimensions, the inclusion will be a cylinder extended along the z-axis.

This creates a four-grain configuration valid in 2 and 3 dimensions. The type allows for no additional keywords.

If initialized in 3 dimensions, the four grains will be extended along the z-axis.

Calculating the orientation field

In many simulations with polycrystals, it is useful to calculate the orientation field. This is a vector field that points in the direction of the crystal orientation. To calculate the orientation field, we rely on the amplitude approximation of the PFC field, demodulate with respect to different orientations and estimate the orientation based on the magnitude of these amplitudes.

Consider the 2D case for a triangular PFC, in which the rotation angle \(\theta\) is given with respect to the z-axis. For a given rotation, we can demodulate the PFC with respect to a triplet of reciprocal lattice vectors.

\[ \{ R(\theta) \mathbf q^{(n)} \}, \]

where \(R(\theta)\) is the rotation matrix and \(\mathbf q^{(n)}\) is the primary reciprocal lattice vector of the unrotated PFC. This will give three amplitudes \(\eta_n\), the combined magnitude

\[ \Phi^2 = \sum_{n=1}^3 |\eta_n|^2 \]

of which will indicate to which degree the orientation of the PFC aligns with that direction.

Straining the PFC

In some simulations, it is useful to prescribe a strain to the PFC. This is done by updating the grid on which the PFC is defined.

In one dimension, the only strain that can be applied is a compression or expansion of the grid, given by \(\mathfrak e_{xx}\).

In this case, these variables are updated according to

self.k[0] = self.k[0]/(1+strain)
self.dif[0] = self.dif[0]/(1+strain)
self.x = self.x*(1+strain)

In two dimensions, the strain tensor is given by three components \(\mathfrak e_{xx}\), \(\mathfrak e_{xy}\), and \(\mathfrak e_{yy}\). In principles, these components cannot be set arbitrarily, since the components must satisfy the compatibilit equations

\[ \partial_y^2 \mathfrak e_{xx} + \partial_x^2 \mathfrak e_{yy} = 2 \partial_x \partial_y \mathfrak e_{xy}. \]

In this case, however, we are limiting ourselves to a constant strain, in which case the compatibility equation is automatically satisfied and the deformation is given by

\[ \begin{array}{ll} u_x = \mathfrak e_{xx} x + \mathfrak e_{xy} y, \\ u_y = \mathfrak e_{xy} x + \mathfrak e_{yy} y. \end{array} \]

The coordinates will transform according to

\[ \begin{array}{ll} x \rightarrow x + u_x, \\ y \rightarrow y + u_y. \end{array} \]

which is is given by

\[ \begin{array}{ll} x' \\ y' \end{array} = \begin{array}{ll} 1 + \mathfrak e_{xx} & \mathfrak e_{xy} \\ \mathfrak e_{xy} & 1 + \mathfrak e_{yy} \\ \end{array} \begin{array}{ll} x \\ y \end{array} \]

We save the components of the strain matrix \(T = \begin{array}{ll} 1 + \mathfrak e_{xx} & \mathfrak e_{xy} \\ \mathfrak e_{xy} & 1 + \mathfrak e_{yy} \\ \end{array} \begin{array}{ll}\) in the distortion_matrix variable. The wave vectors are updated according to

\[ \begin{array}{ll} k_x' \\ k_y' \end{array} = T^{-1} \begin{array}{ll} k_x \\ k_y \end{array} \]

The grid components

Normally, the grid components x, y, and z are set as 1-dimensional arrays, and the same with self.k[0], self.k[1], and self.k[2]. When applying an external strain, however, the grid components will vary in the different directions and will become multi-dimensional arrays. This increases the memory usage of the simulation.

The strain thus contains three components strain[0], strain[1], and strain[2].

, and the grid is updated according to

---

1. Elder, K. R., & Grant, M. (2004). Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals. Physical Review E, 70(5), 051605. https://doi.org/10.1103/PhysRevE.70.051605 ↩↩

2. Skogvoll, V. (2023). Symmetry, topology, and crystal deformations: a phase-field crystal approach. Doctoral Thesis. University of Oslo https://www.duo.uio.no/handle/10852/102731 ↩↩↩↩↩↩

3. Emdadi, A., Asle Z., Mohsen and Asadi, E. (2016). Revisiting Phase Diagrams of Two-Mode Phase-Field Crystal Models. Computational Materials Science. 123, 139-147. https://doi.org/10.1016/j.commatsci.2016.06.018 ↩

4. Wu, K. A., Adland, A. and Karma, A. (2010). Phase-Field-Crystal Model for Fcc Ordering. Physical Review E. 81, 6, 06101. https://doi.org/10.1103/PhysRevE.81.061601 ↩

5. Wu, K-A. and Karma, A. (2007). Phase-Field Crystal Modeling of Equilibrium Bcc-Liquid Interfaces. Physical Review B. 76, 18, 184107. https://doi.org/10.1103/PhysRevB.76.184107 ↩

6. Skogvoll, V., Skaugen, A. and Angheluta, L. (2021). Stress in Ordered Systems: Ginzburg-Landau-type Density Field Theory. Physical Review B. 103, 22, 224107. https://doi.org/10.1103/PhysRevB.103.224107 ↩

7. Skogvoll, V., Salvalaglio, M. and Angheluta, L. (2022). Hydrodynamic Phase Field Crystal Approach to Interfaces, Dislocations and Multi-Grain Networks. Modelling and Simulation in Materials Science and Engineering. https://doi.org/10.1088/1361-651X/ac9493 ↩↩

8. Skogvoll, V., Angheluta, L., Skaugen, A., Salvalaglio, M. and Viñals, J. (2022). A Phase Field Crystal Theory of the Kinematics of Dislocation Lines. Journal of the Mechanics and Physics of Solids. 166, 104932. https://doi.org/10.1016/j.jmps.2022.104932 ↩

9. Dederichs, P. H. and Leibfried, G. (1969). Elastic Green's function for anisotropic cubic crystals. Physical Review. 188, 3, 1175. https://doi.org/10.1103/PhysRev.188.1175 ↩

10. Punke, M., Skogvoll, V., & Salvalaglio, M. (2023). Evaluation of the elastic field in phase-field crystal simulations. PAMM, 23(3), e202300213. https://doi.org/10.1002/pamm.202300213 ↩↩