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

by   Kai Jiang, et al.
Tsinghua University
Xiangtan University

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.


page 1

page 2

page 3

page 4


An efficient method for computing stationary states of phase field crystal models

Computing stationary states is an important topic for phase field crysta...

An accelerated method for computing stationary states of multicomponent phase field crystal model

Finding the stationary states of a free energy functional is an essentia...

Projection Method for Saddle Points of Energy Functional in H^-1 Metric

Saddle points play important roles as the transition states of activated...

A "parallel universe" scheme for crack nucleation in the phase field approach to fracture

Crack nucleation is crucial in many industrial applications. The phase f...

Ground states and their characterization of spin-F Bose-Einstein condensates

The computation of the ground states of spin-F Bose-Einstein condensates...

Field model for complex ionic fluids: analytical properties and numerical investigation

In this paper, we consider the field model for complex ionic fluids with...

Field model for complex ionic fluids

In this paper, we consider the field model for complex ionic fluids with...

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


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


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


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 [32] 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

[2]. 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 

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


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 [28], 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 [18], 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


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.


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 [13], 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


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



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 () [10]. 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


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 [10]. 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



In practice, let , and


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


where and are the discretized interaction and bulk energies:


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


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


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


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


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.


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[8] defined as


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 [16]. Basically, given the current estimation and step size , it updates via


Under suitable assumptions, it is proved in [16] that the iterates has similar convergence property as the traditional proximal gradient method [4] 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


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 [22]. Therefore, we propose a restart algorithm that leads to a convergent algorithm for solving Eq. 18 with energy dissipation property. Given , define


we reset if the following does not hold


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 [2] estimation, i.e.


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


The detailed estimation method is presented in Algorithm 2.

0:  , , , , and
1:  while the stop criterion is not satisfied do
2:     Update
3:     Update
4:     Estimate by Algorithm 2
5:     Calculate via Eq. 22
6:     if Eq. 23 holds then
7:        .
8:     else
9:        Reset .
10:     end if
11:     .
12:  end while
Algorithm 1 AA-BPG Algorithm

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

Definition 2

is called a stationary point of if

Definition 3

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


0:  , , and and
1:  Initialize by BB step Eq. 24.
2:  for  do
3:     Calculate via Eq. 22
4:     if Eq. 25 holds or  then
5:        break
6:     else
8:     end if
9:  end for
10:  Output .
Algorithm 2 Estimation of at

3.3 Convergence property

Throughout this section, we impose the next assumption on .


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 [3].

Lemma 1 ([3])

If is -relative smooth with respect to , then


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


then there exists some such that



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


From the optimal condition of Eq. 28, we have

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


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


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


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


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


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


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

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

where the last inequality is from and .

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


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


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


and implies as . As a result,

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


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


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

Thus, by the convexity of , we have


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 [5].


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



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



The proof is in Appendix A.

It is known from [5] 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


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


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.


The subproblem Eq. 22 is reduced to


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


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


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


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

Using Eq. 47 and Eq. 16, we know

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

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


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


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


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