 # Efficient numerical methods for computing the stationary states of phase field crystal models

Finding the stationary states of a free energy functional is an important problem in phase field crystal (PFC) models. Many efforts have been devoted for designing numerical schemes with energy dissipation and mass conservation properties. However, most existing approaches are time-consuming due to the requirement of small effective time steps. In this paper, we discretize the energy functional and propose efficient numerical algorithms for solving the constrained non-convex minimization problem. A class of first order approaches, which is the so-called adaptive accelerated Bregman proximal gradient (AA-BPG) methods, is proposed and the convergence property is established without the global Lipschitz constant requirements. Moreover, we design a hybrid approach that applies an inexact Newton method to further accelerate the local convergence. One key feature of our algorithm is that the energy dissipation and mass conservation properties hold during the iteration process. Extensive numerical experiments, including two three dimensional periodic crystals in Landau-Brazovskii (LB) model and a two dimensional quasicrystal in Lifshitz-Petrich (LP) model, demonstrate that our approaches have adaptive time steps which lead to a significant acceleration over many existing methods when computing complex structures.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The phase field crystal (PFC) model is an important approach to describe many physical processes and material properties, such as the formation of ordered structures, nucleation process, crystal growth, elastic and plastic deformations of the lattice, dislocations [9, 23]. More concretely, let the order parameter function be , the PFC model can be expressed by a free energy functional

 E[ϕ(r);Θ]=G[ϕ(r);Θ]+F[ϕ(r);Θ], (1)

where are the physical parameters, is the bulk energy with polynomial type or log-type formulation and is the interaction energy that contains higher-order differential operators to form ordered structures [7, 18, 28]. A typical interaction potential function for a domain is

 G[ϕ]=1|Ω|∫Ω[m∏j=1(Δ+q2j)ϕ]2dr,  m∈N (2)

which can be used to describe the pattern formation of periodic crystals, quasicrystals and multi-polynary crystals [18, 24, 20]. In order to understand the theory of PFC models as well as predict and guide experiments, it requires to find stationary states and construct phase diagrams of the energy functional Eq. 1. Denote to be a feasible space, the phase diagram is obtained via solving the minimization problem

 minE[ϕ(r);Θ], s.t. ϕ∈V, (3)

with different physical parameters , which brings the tremendous computational burden. Therefore, within an appropriate spatial discretization, the goal of this paper is to develop efficient and robust numerical methods for solving Eq. 3 with guaranteed convergence while keeping the desired dissipation and conservation properties during the iterative process.

Most existing numerical methods for computing the stationary states of PFC models can be classified into two categories. One is to solve the steady nonlinear Euler-Lagrange equations of

Eq. 3 through different spatial discretization approaches. The other class aims at solving the nonlinear gradient flow equation by using the numerical PDE methods. In these approaches, great efforts have been made for keeping the energy dissipation which is crucial for convergence. Typical energy stable schemes to gradient flows include convex splitting and stabilized factor methods, and recently developed invariant energy quadrature, and scalar auxiliary variable approaches for a modified energy [35, 31, 26, 36, 25]. It is noted that the gradient flow approach is able to describe the quasi-equilibrium behavior of PFC systems. Numerically, the gradient flow is discretized in both space and time domain via different discretization techniques and the stationary state is obtained with a proper choice of initialization.

Under an appropriate spatial discretization scheme, the infinite dimensional problem Eq. 3 can be formulated as a minimization problem in a finite dimensional space. Thus, there may exist alternative numerical methods that can converge to the steady states quickly by using modern optimization techniques. For example, similar ideas have shown success in computing steady states of the Bose-Einstein condensate  and the calculation of density functional theory [30, 19]. In this paper, in order to keep the mass conservation property, an additional constraint is imposed in Eq. 3

and the detail will be given in the next section. Inspired by the recent advances of first order methods which have been successfully applied in image processing and machine learning, we propose an adaptive accelerated Bregman proximal gradient (AA-BPG) method for computing the stationary states of

Eq. 3

. In each iteration, the AA-BPG method updates the estimation of the order parameter function by solving linear equations which have closed form when using the pseudo-spectral discretization and chooses step sizes by using the line search algorithm initialized with the Barzilai-Borwein (BB) method

. Meanwhile, a restart scheme is proposed such that the iterations satisfy energy dissipation property and it is proved that the generated sequence converges to a stationary point of Eq. 3 without the assumption of the existence of global Lipschitz constant of the bulk energy . Moreover, an inexact regularized Newton method is applied for further accelerating the local convergence. More specifically, an inexact preconditioned conjugate gradient method is designed for solving the regularized Newton system efficiently. Extensive numerical experiments have demonstrated that our approach can quickly reach the vicinity of an optimal solution with moderately accuracy, even for very challenge cases.

The rest of this paper is organized as follows. In Section 2, we present the PFC models considered in this paper, and the projection method discretization. In Section 3, we present the AA-BPG method for solving the constrained non-convex optimization with proved convergence. In Section 4, two choices of Bregman distance are proposed and applied for the PFC problems. In Section 5, we design a hybrid method combining with AA-BPG methods and an inexactly regularized Newton method together to further accelerate the calculation. Numerical results are reported in Section 6 to illustrate the efficiency and accuracy of our algorithms. Finally, some concluding remarks are given in Section 7.

## 2 Problem formulation

### 2.1 Physical models

Two classes of PFC models are considered in the paper. The first one is the Landau-Brazovskii (LB) model which can characterize the phase and phase transitions of periodic crystals

. It has been discovered in many different scientific fields, e.g., polymeric materials . In particular, the energy functional of LB model is

 ELB[ϕ(r)]=1|Ω|∫Ω⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ξ22[(Δ+1)ϕ]2G(ϕ)+τ2!ϕ2−γ3!ϕ3+14!ϕ4F(ϕ)⎫⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪⎭dr, (4)

where is a real-valued function which measures the order of system in terms of order parameter. is the bounded domain of system, is the bare correlation length, is the dimensionless reduced temperature, is phenomenological coefficient. Compared with double-well bulk energy , the cubic term in the LB functional helps us study the first-order phase transition.

The second one is the Lifshitz-Petrich (LP) model that can simulate quasiperiodic structures, such as the bi-frequency excited Faraday wave , and explain the stability of soft-matter quasicrystals [17, 11]. Since quasiperiodic structures are space-filling without decay, it is necessary to define the average spacial integral over the whole space as where is the ball centred at origin with radii . The energy functional of LP model is given by

 ELP[ϕ(r)]=−∫⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩c2[(Δ+q21)(Δ+q22)ϕ]2G(ϕ)+ε2ϕ2−κ3ϕ3+14ϕ4F(ϕ)⎫⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪⎭dr, (5)

where is the energy penalty, and are phenomenological coefficients.

Furthermore, we impose the following mean zero condition of order parameter on LB and LP systems to ensure the mass conservation, respectively.

 1|Ω|∫Ωϕ(r)dr=0 or −∫ϕ(r)dr=0. (6)

The equality constraint condition is from the definition of the order parameter which is the deviation from average density. Besides, this equality implies that the system works in canonical ensemble such that the number of particles in the system is conserved.

### 2.2 Projection method discretization

In this section, we introduce the projection method , a high dimensional interpretation approach which can avoid the Diophantine approximation error in computing quasiperiodic systems, to discretize the LB and LP energy functional. It is noted that the stationary states in LB model is periodic, and thus it can be discretized by the Fourier pseudo-spectral method which is a special case of projection method. Therefore, we only consider the projection method discretization of the LP model Eq. 5. We immediately have the following orthonormal property in the average spacial integral sense

 −∫eik⋅re−ik′⋅rdr=δkk′,   ∀k,k′∈Rd. (7)

For a quasiperiodic function, we can define the Bohr-Fourier transformation as



 (8)

In this paper, we carry out the above computation in a higher dimension using the projection method which is based on the fact that a -dimensional quasicrystal can be embedded into an -dimensional periodic structure () . The dimension is the number of linearly independent numbers over the rational number field. Using the projection method, the order parameter can be expressed as

 ϕ(r)=∑h∈Zn^ϕ(h)ei[(P⋅Bh)⊤⋅r],  r∈Rd, (9)

where is invertible, related to the -dimensional primitive reciprocal lattice. The corresponding computational domain in physical space is , . The projection matrix depends on the property of quasicrystals, such as rotational symmetry . If consider periodic crystals, the projection matrix becomes the

-order identity matrix, then the projection reduces to the common Fourier spectral method. The Fourier coefficient

satisfies

 X:={(^ϕ(h))h∈Zn:^ϕ(h)∈C, ∑h∈Zn|^ϕ(h)|<∞}. (10)

In practice, let , and

 XN:={^ϕ(h)∈X:^ϕ(h)=0, for all |hj|>Nj/2, j=1,2,…,n}. (11)

The number of elements in the set is . Together with Eq. 7 and Eq. 9, the discretized energy function Eq. 5 is

 Eh(^Φ)=Gh(^Φ)+Fh(^Φ), (12)

where and are the discretized interaction and bulk energies:

 Gh(^Φ) =c2∑h1+h2=0[q21−(PBh)⊤(PBh)]2[q22−(PBh)⊤(PBh)]2^ϕ(h1)^ϕ(h2) (13) Fh(^Φ) =ε2∑h1+h2=0^ϕ(h1)^ϕ(h2)−κ3∑h1+h2+h3=0^ϕ(h1)^ϕ(h2)^ϕ(h3) +14∑h1+h2+h3+h4=0^ϕ(h1)^ϕ(h2)^ϕ(h3)^ϕ(h4),

and , , , , . It is clear that the nonlinear terms in are -dimensional convolutions in the reciprocal space. A direct evaluation of these convolution terms is extremely expensive. Instead, these terms are simple multiplication in the -dimensional physical space. Similar to the pseudo-spectral approach, these convolutions can be efficiently calculated through FFT. Moreover, the mass conservation constraint Eq. 6 is discretized as

 e⊤1^Φ=0, (14)

where . Therefore, we obtain the following finite dimensional minimization problem

 min^Φ∈CNEh(^Φ)=Gh(^Φ)+Fh(^Φ), s.t. e⊤1^Φ=0. (15)

For simplicity, we omit the subscription in and in the following context. According to Eq. 13, denoting as the discretized Fourier transformation matrix, we have

 ∇G(^Φ)=D^Φ, ∇F(^Φ)=F−1NΛFN^Φ (16) ∇2G(^Φ)=D, ∇2F(^Φ)=F−1NΛ(′)FN, (17)

where is a diagonal matrix with nonnegative entries and are also diagonal matrices but related to . In the next section, we propose the adaptive accelerated Bregman proximal gradient (AA-BPG) method for solving the constrained minimization problem Eq. 15.

## 3 The AA-BPG method

Consider the minimization problem that has the form

 minxE(x)=f(x)+g(x), (18)

where is proper but non-convex and is proper, lower semi-continuous and convex. Let the domain of to be , we make the following assumptions.

###### Assumption

is bounded below and for any , the sub-level set is compact.

Let be a strongly convex function such that and . Then, it induces the Bregman divergence defined as

 Dh(x,y)=h(x)−h(y)−⟨∇h(y),x−y⟩, ∀(x,y)∈domh×intdomh. (19)

It is noted that and if and only if due to the strongly convexity of . Furthermore, as . In recent years, Bregman distance based proximal methods [3, 6] have been proposed and applied for solving the Eq. 18 in a general non-convex setting . Basically, given the current estimation and step size , it updates via

 xk+1=argminx{g(x)+⟨x−xk,∇f(xk)⟩+1αkDh(x,xk)}. (20)

Under suitable assumptions, it is proved in  that the iterates has similar convergence property as the traditional proximal gradient method  while iteration Eq. 20 does not require the Lipschitz condition on . Motivated by the Nesterov acceleration technique [29, 4], we add an extrapolation step before Eq. 20 and thus the iterate becomes

 yk =xk+wk(xk−xk−1), (21) xk+1 =argminx{g(x)+⟨x−yk,∇f(yk)⟩+1αkDh(x,yk)},

where . It is noted that the minimization problems in Eq. 20 and Eq. 21 are well defined and single valued as is convex and is strongly convex. Although the extrapolation step accelerates the convergence in some cases, it may generate the oscillating phenomenon of the objective value that slows down the convergence . Therefore, we propose a restart algorithm that leads to a convergent algorithm for solving Eq. 18 with energy dissipation property. Given , define

 zk=argminx{g(x)+⟨x−yk,∇f(yk)⟩+1αkDh(x,yk)}, (22)

we reset if the following does not hold

 E(xk)−E(zk)≥c∥xk−xk+1∥2 (23)

for some constant . In the next section, we will show that Eq. 23 holds when . Overall, the AA-BPG algorithm is presented in Algorithm 1.

Step size estimation. In each step, is chosen adaptively by backtracking linear search method which is initialized by the BB step  estimation, i.e.

 αk=⟨sk,sk⟩⟨sk,vk⟩ or ⟨vk,sk⟩⟨vk,vk⟩, (24)

where and . Let be a small constant and be obtained from Eq. 22, we adopt the step size whenever the following inequality holds

 E(yk)−E(zk)≥η∥yk−zk∥2. (25)

The detailed estimation method is presented in Algorithm 2.

### 3.1 Convergence analysis

In this section, we focus on the convergence analysis of the proposed AA-BGP method. Before proceeding, we introduce some definitions used in analysis.

### 3.2 Definitions

###### Definition 1

Let be a proper and lower semi-continuous function and convex function. The subgradient of at is defined as

 ∂g(x)={u:f(y)−f(x)−⟨u,y−x⟩≥0,∀y∈domg}.

###### Definition 2

is called a stationary point of if

 0∈∂E(x)=∇f(x)+∂g(x).

###### Definition 3

A function is -relative smooth if there exists a strongly convex function such that

 Lf∇2h(x)−∇2f(x)⪰0,∀x∈intdomh. (26)

### 3.3 Convergence property

Throughout this section, we impose the next assumption on .

###### Assumption

There exists such that is -relative smooth with respect to a strongly convex function .

Under the Section 3.3, we have the following useful lemma as stated in .

###### Lemma 1 ()

If is -relative smooth with respect to , then

 f(x)−f(y)−⟨∇f(y),x−y⟩≤LfDh(x,y),∀x,y∈intdomh. (27)

Based on the above Lemma, the descent property of the iteration generated by Bregman proximal operator Eq. 22 is established as follows.

###### Lemma 2

Let and suppose the Section 3.3 holds. If

 z=argminx{g(x)+⟨x−y,∇f(y)⟩+1αDh(x,y)}, (28)

then there exists some such that

 E(y)−E(z)≥(1α−Lf)σ2∥z−y∥2. (29)

###### Proof

Since is strongly convex, there exists some constant such that is convex. Then, and we have

 Dh(z,x)=h(z)−h(y)−⟨∇h(y),z−y⟩≥σ2∥z−y∥2. (30)

From the optimal condition of Eq. 28, we have

 E(y) =f(y)+g(y)=[f(y)+⟨∇f(y),x−y⟩+1αDh(x,y)+g(x)]x=y ≥f(y)+⟨∇f(y),z−y⟩+1αDh(z,y)+g(z) ≥f(z)−LfDh(z,y)+1αDh(z,y)+g(z) =E(z)+(1α−Lf)Dh(z,y)≥E(z)+(1α−Lf)σ2∥z−y∥2,

where the second inequality is from Eq. 27 and the last inequality follows from Eq. 30.

###### Remark

Lemma 2 shows that the non-restart condition Eq. 23 and the linear search condition Eq. 25 are satisfied when

 0<α<¯α:=min(12c/σ+Lf,12η/σ+Lf) and 0<αmin≤¯α (31)

Therefore, the line search in Algorithm 2 stops in finite iterations, and thus the Algorithm 1 is well defined.

In the following analysis, we always assume that the parameter satisfies Eq. 31 for simplicity. Therefore, we can obtain the sufficient decrease property of the sequence generated by Algorithm 1.

###### Corollary 1

Suppose the Section 3 and Section 3.3 hold. Let be the sequence generated by the Algorithm 1. Then, and

 E(xk)−E(xk+1)≥c0∥xk−xk+1∥2, (32)

where .

The proof of Corollary 1 is a straightforward result as AA-BPG algorithm is well defined and the condition Eq. 23 or Eq. 25 holds at each iteration. Let be the closed ball that contains . Since , there exist such that

 ρh=supx∈B(x0)∥∇2h(x)∥,ρf=supx∈B(x0)∥∇2f(x)∥. (33)

Thus, we can show the subgradient of each step generated by Algorithm 1 is bounded by the movement of .

###### Lemma 3 (Bounded the subgradient)

Suppose Section 3 and Section 3.3 holds. Let be the sequence generated by Algorithm 1. Then, there exists such that

 dist(0,∂E(xk+1))≤c1(∥xk+1−xk∥+¯w∥xk−xk−1∥), (34)

where and , are defined as Eq. 33 and are constants defined in Algorithm 1 and Algorithm 2, respectively.

###### Proof

By the first order optimality condition of Eq. 21, we get

 0∈∇f(yk)+1αk(∇h(xk+1)−∇h(yk))+∂g(xk+1) ⟺ −∇f(yk)−1αk(∇h(xk+1)−∇h(yk))∈∂g(xk+1)

From Lemma 2, we have , then . It follows that

 dist(0,∂E(xk+1)) =infy∈∂g(xk+1)∥∇f(xk+1)+y∥ ≤∥∇f(xk+1)−∇f(yk)−1αk(∇h(xk+1)−∇h(yk))∥ ≤∥∇f(xk+1)−∇f(yk)∥+1αk∥∇h(xk+1)−∇h(yk)∥ ≤(ρf+ρhαk)∥xk+1−yk∥ ≤c1(∥xk+1−xk∥+¯w∥xk−xk−1∥),

where the last inequality is from and .

Now, we are ready to establish the sub-convergence property of Algorithm 1.

###### Theorem

Suppose Section 3 and Section 3.3 hold. Let be the sequence generated by Algorithm 1. Then, for any limit point of , we have .

###### Proof

From Corollary 1, we know and thus bounded. Then, the set of limit points of is nonempty. For any limit point , there exist a subsequence such that as . We know is a decreasing sequence. Together with the fact that is bounded below, there exists some such that as . Moverover, it has

 E(x0)−¯E=limK→∞K∑j=0(E(xj)−E(xj+1))≥c0limK→∞K∑j=0∥xj−xj+1∥2, (35)

and implies as . As a result,

 limk→∞∥xk−yk−1∥≤limk→∞(∥xk−xk−1∥+¯ω∥xk−1−xk−2∥)=0.

Together with Eq. 34, it implies that there exists such that

 limj→∞∥∇f(xkj)+ukj∥=0⇒limj→∞ukj=−∇f(x∗), (36)

as is continuous and when .

In the next, we prove . It is easy to know that for finite since . Thus, we have as . From Eq. 21, we know

 g(xkj)+⟨xkj−ykj−1,∇f(ykj−1)⟩+1αkD(xkj,ykj−1)≤g(x)+⟨x−ykj−1,∇f(ykj−1)⟩+1αkD(x,ykj−1),∀x. (37)

Let and , we get . By the fact that is lower semi-continuous, it has .

Thus, by the convexity of , we have

 (38)

Let in (38) and using the facts , as and (36), we have and thus .

Furthermore, the sub-sequence convergence can be strengthen by imposing the next assumption on which is known as the Kurdyka-Lojasiewicz (KL) property .

###### Assumption

is the function, i.e. for all , there exist , a neighborhood of and such that for all , the following inequality holds,

 ψ′(E(x)−E(¯x))dist(0,∂E(x))≥1. (39)

###### Theorem

Suppose Section 3, Section 3.3 and Section 3.3 hold. Let be the sequence generated by Algorithm 1. Then, there exists a point such that

 limk→+∞xk=x∗,0∈∂E(x∗). (40)

###### Proof

The proof is in Appendix A.

It is known from  that many functions satisfy Section 3.3 including the energy function in PFC models. In the following context, we apply the AA-BPG method for solving the PFC models Eq. 15 by introducing two Bregman distances.

## 4 AA-BPG method for solving PFC models

The problem Eq. 15 can be reduced to Eq. 18 by setting

 (41)

where and if and otherwise. The main difficulty of applying Algorithm 1 is solving the subproblem Eq. 22 efficiently. In this section, two different strongly convex functions are chosen as

 h(x)=12∥x∥2  (P2)  and h(x)=a4∥x∥4+b2∥x∥2+1  (P4), (42)

where and (P2) and (P4) represent the highest order of the norm.

Case (P2). The Bregman distance of is reduced to the Euclidean distance, i.e.

 Dh(x,y)=12∥x−y∥2. (43)

The subproblem Eq. 22 is reduced to

 min^Φ G(^Φ)+⟨∇F(^Ψk),^Φ−^Ψk⟩+12αk∥^Φ−^Ψk∥2, s.t. e⊤1^Φ=0, (44)

where . Although the Eq. 44 is a constrained minimization problem, it has a closed form solution based on our discretization which leads to a fast computation.

###### Lemma 4

Given , if , the minimizer of Eq. 44, denoted by , is given by

 ^Φk+1=(I+αkD)−1(^Ψk−αkP1∇F(^Ψk)), (45)

where is defined in Eq. 16 and is the projection into the set .

###### Proof

The KKT conditions for this subproblem Eq. 22 can be written as

 ∇G(^Φk+1)+∇F(^Ψk)+1αk(^Φk+1−^Ψk)−ξke1=0, (46) e⊤1^Φk+1=0, (47)

where is the Lagrange multiplier. Taking the inner product with in Eq. 46, we obtain

 ξk=e⊤1(∇G(^Φk+1)+∇F(^Φk)−1αk^Ψk).

Using Eq. 47 and Eq. 16, we know

 e⊤1∇G(^Φk+1)=e⊤1(D^Φk+1)=0.

Together with , we have . Substituting it into Eq. 46, it follows that

 ^Φk+1 =(αkD+I)−1(^Ψk−αkP1∇F(^Ψk)).

It is noted that from the proof of Lemma 4, the feasibility assumption holds as long as which can be set in the initialization. The detailed algorithm is given in Algorithm 3 with .

Case (P4). In this case, the subproblem Eq. 22 is reduced to

 min^ΦG(^Φ)+⟨∇F(^Ψk),^Φ−^Ψk⟩+Dh(^Φ,^Ψk), s.t. e⊤1^Φ=0. (48)

where . The next lemma shows the optimal condition of minimizing Eq. 48.

###### Lemma 5

Given . If , the minimizer of Eq. 48, denoted by , is given by

 ^Φk+1=[αkD+(ap∗+b)I]−1(∇h(^Ψk)−αkP1∇F(^Ψk)), (49)

where is given in Eq. 16 and is a fixed point of .

###### Proof

The KKT conditions of Eq. 48 imply that there exists a Lagrange multiplier such that