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 ofEq. 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 ofEq. 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
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
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 , 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 () . 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 . If consider periodic crystals, the projection matrix becomes the
-order identity matrix, then the projection reduces to the common Fourier spectral method. The Fourier coefficientsatisfies
In practice, let , and
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  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 . Basically, given the current estimation and step size , it updates via
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
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
we reset if the following does not hold
Step size estimation. In each step, is chosen adaptively by backtracking linear search method which is initialized by the BB step  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.
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.
Let be a proper and lower semi-continuous function and convex function. The subgradient of at is defined as
is called a stationary point of if
A function is -relative smooth if there exists a strongly convex function such that
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 .
Lemma 1 ()
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.
Let and suppose the Section 3.3 holds. If
then there exists some such that
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)
Now, we are ready to establish the sub-convergence property of Algorithm 1.
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 .
Furthermore, the sub-sequence convergence can be strengthen by imposing the next assumption on which is known as the Kurdyka-Lojasiewicz (KL) property .
is the function, i.e. for all , there exist , a neighborhood of and such that for all , the following inequality holds,
The proof is in Appendix A.
4 AA-BPG method for solving PFC models
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.
Case (P4). In this case, the subproblem Eq. 22 is reduced to
where . The next lemma shows the optimal condition of minimizing Eq. 48.
The KKT conditions of Eq. 48 imply that there exists a Lagrange multiplier such that