Class: Phase Field Crystal
In this class, we simulate a crystal using the phasefield crystal methodology. The phasefield 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 realvalued 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.
Variables and parameters
The phase field \(\psi\).
Basic model
The PFC methodology is based on postulating a free energy
The PFC free energy
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
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: Phasefield crystal models. \(\mathcal L_X = (X+\nabla^2)\).
Historical context: model presented in Ref.^{1}.
The primary field is \(\psi\), a real valued field
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, \(\{\vec a^{(n)}\}_{n=1}^2\) are primitive lattice vectors and \(\{ \vec q^{(n)}\}_{n=1}^2\) primitive reciprocal lattice vectors (RLVs) that satisfy \(\vec a^{(n)} \cdot \vec q^{(m)} = 2\pi \delta_{nm}\). Amended and reprinted from Ref. [^Skogvoll2023SymmetryTopology] 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
the reciprocal lattice is \(d\) dimensions generated by \(d\) primitive RLVs \(\lbrace\boldsymbol{q}^{(n)}\rbrace_{n=1}^d\), satisfying
The primitive BLVRLV orthogonality relation.
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
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
and any field with that periodicity can be expressed as a Fourier series in terms of these RLVs by
Cutting the Fourier series off at the primary RLVs mode is called the onemode approximation
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.
Lattice constant
Primary RLVs
Primary BLVs
Lattice constant
Primary RLVs
Primary BLVs
Lattice constant $$ a_0 = 2\pi $$
Primary RLVs
Primary Bravais Lattice vectors
Lattice constant
Primary BLVs
Primary RLVs
Lattice constant
Primary BLVs
Primary RLVs
Lattice constant
Primary BLVs
Primary RLVs
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 (coarsegrain) 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 onemode approximation, i.e.,
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]\)
To insert this into the free energy density, we need to calculate some auxiliary quantities.
so
Now, to calculate this, we need to make use of a condition of resonance, which is that under the coarsegraining operation, terms with periodicity average to zero^{2}. So we get
Similarly, the other terms give
which gives
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.
which gives \(A=0\) or
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
The pfc
instance will be initialized with a list called eta0
, which consists of the equilibrium values of the amplitudes to begin with.
Default resolution:
Default model parameters \((r,\psi_0)\): \((0.3,0.3)\)
Free energy per unit cell (calculation document):
Equilibrium amplitude (conserved)
Equilibrium amplitude equations (unconserved)
Default resolution: \([7,12]^{1}a_0\)
Default model parameters \((r,\psi_0)\): \((0.3,0.3)\)
Free energy per unit cell (calculation document.):
Equilibrium amplitude equations (conserved)
Equilibrium amplitude equations (unconserved)
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):
Equilibrium amplitude equations
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):
Equilibrium amplitude (conserved)
Equilibrium amplitude equations (unconserved)
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):
Equilibrium amplitude equations
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):
Equilibrium amplitude equations:
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\).
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 fewmode 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
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
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. ^{2} with permission.
Elasticity
This quantity is the Poisson ratio, under some conditions. I think. To be fixed. \(\(\nu = \frac{\lambda}{(d1)\lambda + 2\mu + \gamma}\)\)
Elastic constants are saved in these quantities
Stress tensor
The equations for calculating the stress tensor are found in Ref. ^{6}, 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\).
Stress tensor
Elastic constants
Stress tensor
Elastic constants
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
Stress tensor
Elastic constants
Stress tensor
Elastic constants
Stress tensor
Elastic constants
Strains
In order to calculate the strain of a PFC configuration, we follow a generalized procedure as that outlined in Ref. ^{2}. The ground state of the PFC can be written as
where \(\mathcal R^{(n)}\) is the nth closest modes, i.e., the full reciprocal lattice can be written as
We define the following quantity
Assuming that the modeapproximations 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
The strain of the phasefield crystal can be calculated using the following formula The strain of the phasefield crystal can be calculated using the following formula^{2}
where \(\mathcal S = \langle \nabla \psi \nabla \psi \rangle\) is called the structure function ^{7}
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 ^{10}.
It is therefore often useful to subtract the mean value for visualization purposes.
Dislocations
The Burgers vector is defined by
Burgers vector definition
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 onebody 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 (\(\vec b \perp \vec t\)), and (c) a 3D simple cubic lattice with a screw dislocation (\(\vec b \parallel \vec t)\).
In all cases, a circulation (green) that is righthanded with respect to the tangent vector \(\vec t\), i.e., following a path around the dislocation, gives rise to a connection error: the Burgers vector \(\vec b\).
Reprinted from Ref. [^Skogvoll2023SymmetryTopology] 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
where \(s_n\) is the charge associated with the Burgers vector \(\boldsymbol{b}\)
Dislocation charge
Using the primary BLVs we get the charges summarized below
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\) 
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\) 
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\) 
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\) 
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 ^{8}
The dislocation density tensor
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
In the language of this
Equations of motion
Conserved evolution
The basic equation motion presented in Ref. ^{1} can be derived by postulating a simple mechanism for free energy minimzation under the constraint of mass conservation.
The PFC evolution equation (evolve_PFC
)
where
The PFC chemical potential
which gives $$ {{\omega }{\mathfrak f}}= \boldsymbol{k}^2 (\texttt r + {{\mathcal L}^2) $$}
Unconserved dynamics
The unconserved PFC evolution equation (evolve_PFC_unconserved
)
Evolution at mechanical equilibrium
The PFC evolution at mechanical equilibrium (evolve_PFC_mechanical_equilibrium
)
Step 1:
Step 2:
where \(\boldsymbol{u}^\delta\) is the solution to
where
and the elastic constants tensor \(\mathcal C\) is given by
This equation is solved in Fourier space by
where the Greens function \(G_{\mathfrak f ~ ij}\) is given in Ref. ^{9} as
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. ^{7}, a hydrodynamic approach was derived. A simplified twoparameter model was proposed
The hydrodynamic PFC model (evolve_PFC_hydrodynamic
)
We insert for \(\tilde \mu_c\) and write it in matrix form to emphasize the linear and nonlinear parts
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 preconfigured 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
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 (defaultself.Lmin/4
). 
position
which specifies the position of the center of the inclusion (defaultself.rmid
). 
rotation
a float specifying the rotation of the inclusion wrt. the zaxis (defaultnp.pi/6
).
If initialized in 3 dimensions, the inclusion will be a cylinder extended along the zaxis.
This creates a fourgrain 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 zaxis.
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 zaxis. For a given rotation, we can demodulate the PFC with respect to a triplet of reciprocal lattice vectors.
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
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 \(\matfrak e_{xx}\).
In this case, these variables are updated according to
In two dimensions, the strain tensor is given by three components \(\matfrak e_{xx}\), \(\matfrak e_{xy}\), and \(\matfrak e_{yy}\). In principles, these components cannot be set arbitrarily, since the components must satisfy the compatibilit equations
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
The coordinates will transform according to
which is is given by
We save the components of the strain matrix \(T = \begin{array}{ll}
1 + \matfrak e_{xx} & \matfrak e_{xy} \\
\matfrak e_{xy} & 1 + \matfrak e_{yy} \\
\end{array}
\begin{array}{ll}\) in the strain_matrix
variable.
The wave vectors are updated according to
The grid components
Normally, the grid components x
, y
, and z
are set as 1dimensional 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 multidimensional 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

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 ↩↩

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

Emdadi, A., Asle Z., Mohsen and Asadi, E. (2016). Revisiting Phase Diagrams of TwoMode PhaseField Crystal Models. Computational Materials Science. 123, 139147. https://doi.org/10.1016/j.commatsci.2016.06.018 ↩

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

Wu, KA. and Karma, A. (2007). PhaseField Crystal Modeling of Equilibrium BccLiquid Interfaces. Physical Review B. 76, 18, 184107. https://doi.org/10.1103/PhysRevB.76.184107 ↩

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

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

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 ↩

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 ↩

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