Projection Method for Saddle Points of Energy Functional in H −1 Metric
aa r X i v : . [ m a t h . NA ] N ov Projection Method for Saddle Points of Energy Functional in H − Metric
Shuting Gu , Ling Lin , Xiang Zhou Abstract
Saddle points play important roles as the transition states of activated process ingradient system driven by energy functional. However, for the same energy functional,the saddle points, as well as other stationary points, are different in different metricssuch as the L metric and the H − metric. The saddle point calculation in H − metricis more challenging with much higher computational cost since it involves higher orderderivative in space and the inner product calculation needs to solve another Possionequation to get the ∆ − operator. In this paper, we introduce the projection ideato the existing saddle point search methods, gentlest ascent dynamics (GAD) anditerative minimization formulation (IMF), to overcome this numerical challenge dueto H − metric. Our new method in the L metric only by carefully incorporates asimple linear projection step. We show that our projection method maintains the sameconvergence speed of the original GAD and IMF, but the new algorithm is much fasterthan the direct method for H − problem. The numerical results of saddle points inthe one dimensional Ginzburg-Landau free energy and the two dimensional Landau-Brazovskii free energy in H − metric are presented to demonstrate the efficiency ofthis new method. Keywords : saddle point, transition state, projection method, gentlest ascent dy-namics
Mathematics Subject Classification (2010)
Primary 65K05, Secondary 82B051.
Introduction
Saddle points have important physical meaning and have been of broad interest inchemistry, physics, biology and material sciences. In computational chemistry [22],one of the most important objects on the potential energy surface is the transitionstate, a special type of the saddle point with index-1, which is defined as the criticalpoint with only one unstable direction. Such transition states are the bottlenecks onthe most probable transition paths between different local wells. In recent years, alarge number of numerical methods have been proposed and developed to efficientlycompute these saddle points. Generally speaking, there are two classes: path-findingmethods and surface-walking methods. The former includes the string method [20, 9]and the nudged elastic band method [15]. These methods are to search the so-calledminimum energy path (MEP). The points along the MEP with locally maximum energyvalue are then the index-1 saddle points. The later methods include the eigenvectorfollowing method [5], the dimer method [14], the activation-relaxation techniques [19], the gentlest ascent dynamics(GAD) [10] and the iterative minimization formulation(IMF) [11, 12]. They evolve a single state on the potential energy surface along theunstable direction, for example, the min-mode direction.There are different fixed points on different potential energy surfaces. Here we willaddress that even for the same energy functional, different stationary points (metastablestates) and saddle points can be obtained in different metrics such as the L metricand the H − metric. We take the Ginzburg-Landau free energy on a bounded domainΩ for example F ( φ ) = Z Ω h κ |∇ φ ( x ) | + f ( φ ) i dx, f ( φ ) = ( φ − / , (1)and the following two gradient flows are commonly used in physics models, dependingon which metric is used for the gradient.(1) In L metric: the (non-conserved) Allen-Cahn (AC) equation [1] ∂φ∂t = − δFδφ ( φ ) = κ ∆ φ − ( φ − φ ); (2)and(2) In H − metric: the (conserved) Cahn-Hilliard (CH) equation [4] ∂φ∂t = ∆ δFδφ = − κ ∆ φ + ∆( φ − φ ) . (3)Here δFδφ is the first order variation of F in the L sense. Nowadays, the Allen-Cahn andCahn-Hilliard equations have been widely used in many complicated moving interfaceproblems in materials science and fluid dynamics through a phase-field approach, forinstance, [21, 6, 2, 3].The inner product and the norm in H − metric can be rewritten in terms of the L product as follows: k φ k H − = (cid:10) ( − ∆) − φ, φ (cid:11) L , h φ, ψ i H − = (cid:10) ( − ∆) − φ, ψ (cid:11) L , (4)where ( − ∆) − , a bounded positive self-adjoint linear operator, is the inverse of − ∆subject to certain boundary condition[7]. The dynamics (2) and (3) are the gradientflows of the same energy functional (1) in L metric and H − metric, respectively. It isclear that these two gradient flows have distinctive dynamics and properties. The Cahn-Hilliard equation (3) preserves the mass R Ω φ dx while the Allen-Cahn does not. We areintertested in the stationary states of these two dynamics. With the same boundarycondition, the stationary states of dynamics (2) (with the sufficient regularity such asin the Sobolev H (Ω) space) are the stationary states of the dynamics (3), but not viceversa.It takes more computational cost to calculate the stationary points in CH equationthan that in AC equation, because the dynamics in the H − metric (3) is a fourthorder derivative equation in space, two order higher than that in the L metric (2).What is worse is that any computation involving the inner product calculation in H − metric needs to calculate the ∆ − operator (see (4)) by solving a Poisson equation. Soif one only wants to locate the fixed points(stationary points or saddle points) insteadof capturing the time evolution in H − metric, it is much less efficient to use dynamicsin H − metric such as the CH equation.Since the main difference between the dynamics in L metric and H − metric iswhether the mass is conservation, our idea to handle the above challenges is to add mass conservation constrains into the L metric dynamics. This conservation canbe enforced by a projection operator. In the work of [18], the projected Allen-Cahnequation ∂u∂t = P ( − δFδu ) (5)was proposed as a counterpart of the Cahn-Hillard equation to search different phasesin diblock copolymers. P in (5) is the orthogonal projection operator onto the confinedsubspace satisfying the mass conservation. The Cahn-Hilliard equation (3) and theprojected Allen-Cahn equation (5) then both preserve the mass, although the gradient-descent trajectories and the transition paths are different[25]. One important fact isthat (5) and (3) share the same stationary points (metastable states) and the saddlepoints if they have the same mass. [25] further compared the stochastic models arisingfrom these two dynamics ((3) and (5)) for the noise-induced transitions. They showedthe subtle difference in transition rates and minimum energy paths in the two stochasticmodels. For our purpose of locating the saddle point in this article, we utilize theequivalence of saddle points of (3) and (5) and solve the saddle points of the Cahn-Hilliard equation (3) by solving the projected Allen-Cahn equation (5) (in L metric).Compared with the stationary points, people are more concerned about the saddlepoints for rare event study. In [13], the IMF has been applied to locate the saddle pointof an energy functional in H − metric directly. However, as mentioned before, this di-rectly is quite expensive in computation. Considering the equivalence of the fixed pointsfor CH equation and projected AC equation, we propose to locate the saddle points ofthe AC equation with the mass conservation constrain. Recently, several methods havebeen developed to locate the saddle point with constrains. [8] developed a constrainedstring method for finding the saddle points subject to constraints. [24] studied the con-strained shrinking dimer dynamics to locate saddle points associated with an energyfunctional defined on a constrained manifold. [16] considered noise-induced transitionpaths in randomly perturbed dynamical systems on a smooth manifold. Besides, thepapers of iterative minimization formulation (IMF) [11, 12] have included the discus-sions on the projection idea for saddle point on manifold. But in these works, theconstraints are externally imposed and thus require higher computational cost thanthe unconstrained problems. Our motivation here is totally different. The questionwe considered here is essentially an unconstrained problem since the mass is conservedautomatically in H − metric. We transform a difficult unconstrained problem into aless challenging constraint problem and we only work on the orthogonal projection formass conservation. This method can reduce the computational cost efficiently since itcan not only avoid a higher order equation solving but also escape from the ∆ − oper-ator calculation. Furthermore, we verify that the projected IMF can ensure the sameconvergence rate as the original IMF. Finally, we remind the readers that if one reallylooks for the noise-induced transition paths in H − sense, then the true dynamics likethe CH equation (3) is still necessary, although our method for saddle points can assistthis path-finding task; see details in [25].The paper is organized as follows. Section 2 is a short review of two main methodsfor saddle points: the IMF and the GAD. In section 3, we first present the applicationof the IMF in the H − metric, and then propose the mathematical formulation ofthe projected IMF and the convergence result of the projected IMF. The projectedGAD is also presented here. In section 4, in order to validate the efficiency of our newmethod, we test two numerical examples: the saddle points of the one dimensional Ginzburg-Landau free energy and the two dimensional Landau-Bravoskii free energyin H − metric. Finally we make the conclusion.2. Review
In this section, we will review two main methods for saddle points: the IMF and theGAD, from which the projected IMF and the projected GAD in the next Section willbe proposed.2.1.
Iterative minimization formulation(IMF).
We first review the iteration min-imization formulation (IMF) in [11]. Suppose M is a function space equipped with thenorm k·k and the inner product h· , ·i . The IMF to locate the saddle point of an energyfunctional F ( φ ) is the following iteration v ( k +1) = argmin k v k =1 D φ, H ( φ ( k ) ) φ E , (6) φ ( k +1) = argmin φ L ( φ ; φ ( k ) , v ( k +1) ) , (7)where H = δ φ F is the second order variational operator of F , and L ( φ ; φ ( k ) , v ( k +1) ) = (1 − α ) F ( φ ) + αF (cid:16) φ − D v ( k +1) , φ − φ ( k ) E v ( k +1) (cid:17) − βF (cid:16) φ ( k ) + D v ( k +1) , φ − φ ( k ) E v ( k +1) (cid:17) . (8) α and β are two parameters, and α + β >
1. Two special choices for α and β are:(i) ( α, β ) = (2 , , then L ( φ ; φ ( k ) , v ) = − F ( φ ) + 2 F ( φ − (cid:10) v, φ − φ ( k ) (cid:11) v ); (ii) ( α, β ) =(0 , , then L ( φ ; φ ( k ) , v ) = F ( φ ) − F ( φ ( k ) + (cid:10) v, φ − φ ( k ) (cid:11) v ). (6) is called the “rotationstep” and (7) is the “translation step”. The main properties of the auxiliary objectivefunctional L ( φ ; φ ( k ) , v ) when α + β > Theorem 1 ([11]) . Suppose that φ ∗ is a (non-degenerate) index-1 saddle point of thefunctional F ( φ ) , and the auxiliary functional L is defined by (8) with α + β > , then (1) a neighbourhood U of φ ∗ exists such that for any φ ∈ U , L ( φ ; φ ( k ) , v ) is strictlyconvex in φ ∈ U and thus has a unique minimum in U ; (2) define the mapping Φ : φ ∈ U → Φ( φ ) ∈ U to be the unique minimizer of L in U for any φ ∈ U . Further assume that U contains no other stationary points of F exceptfor φ ∗ . Then the mapping Φ has only one fixed point φ ∗ ; (3) the mapping φ → Φ( φ ) has a quadratic convergence rate. Gentlest ascent dynamics(GAD).
The GAD for a gradient system ˙ φ = − δ φ F ( φ )is ˙ φ = − δ φ F ( φ ) + 2 h δ φ F ( φ ) , v ih v, v i v, (9a) γ ˙ v ( t ) = − δ φ F ( φ ) v + (cid:10) v, δ φ F ( φ ) v (cid:11) v, (9b)where δ φ F ( φ ) is the second order variational derivative of the energy functional F ( φ ). γ > γ means a fast dynamics for the directionvariable v ( t ) towards the steady state. For a frozen φ , this steady state is the minmode of δ φ F ( φ ): the eigenvector corresponds to the smallest eigenvalue of δ φ F ( φ ). Theorem 2 ([10]) . The (linearly) stable critical point of the GAD (9) corresponds tothe index-1 saddle point of the original dynamics ˙ φ = − δ φ F ( φ ) , i.e., (1) If ( φ ∗ , v ∗ ) is a stable critical point of the GAD, then φ ∗ is a saddle point of F ( φ ) ; (2) If φ ∗ is an index-1 saddle point of F ( φ ) with the eigenvector v ∗ , then ( φ ∗ , v ∗ ) isa stable critical point of the GAD. Remark 1.
In the IMF, there are two levels of iterations: the rotation step and thetranslation step. In general, it requires many iteration steps to get φ ( k +1) for thetranslation step, but it is not necessary to do so in practice. If the two subproblems ofthe IMF moves forward only one iteration step, the IMF becomes exactly the GAD incontinuous time limit minimizations. Main methods
We present the main methods of projection here by starting with the formulation inthe H − space where the inner product h· , ·i becomes h· , ·i H − .3.1. The IMF in H − metric. Formally, the IMF in the spatially extended systemto locate the saddle point of F ( φ ) , φ ∈ C (Ω) in H − metric is: v ( k +1) = argmin k v k H − =1 D v, e H ( φ ( k ) ) v E H − , (10) φ ( k +1) = argmin L ( φ ; φ ( k ) , v ( k +1) ) , (11)where e H = δ φ F ( φ ) | H − = − ∆ δ φ F ( φ ) = − ∆ H , (12)and L ( φ ; φ ( k ) , v ( k +1) ) = (1 − α ) F ( φ ) + αF (cid:16) φ − D v ( k +1) , φ − φ ( k ) E H − v ( k +1) (cid:17) − βF (cid:16) φ ( k ) + D v ( k +1) , φ − φ ( k ) E H − v ( k +1) (cid:17) . (13)Recall H = δ φ F ( φ ) is the second order variational operator of F w.r.t. φ in L metric.For convenience, in this paper, we take α = 0, β = 2, then L ( φ ) = F ( φ ) − F ( ˆ φ ) , with ˆ φ = φ ( k ) + D v ( k +1) , ( φ − φ ( k ) ) E H − v ( k +1) , (14)where the inner product in H − metric is defined by the L product: h u, v i H − = (cid:10) ( − ∆) − u, v (cid:11) L .In H − metric, φ is mass conserved, R Ω φ ( x ) dx = m . So any eigenvectors of e H satisfies R Ω ψ ( x ) dx = 0. Thus the eigenvalue problem (10) can be rewritten as ( e H ( φ ) ψ = λψ, R Ω ψ ( x ) dx = 0 , (15)subject to some boundary condition. Define the Rayleigh quotient e R ( ψ ) = D ψ, e H ψ E H − k ψ k H − , and thus the min-mode is the minimizer of the problemargmin ψ (cid:26) e R ( ψ ) : Z Ω ψ dx = 0 , k ψ k H − = 1 (cid:27) . (16)After the min-mode is obtained, the subproblem of minimizing the auxiliary functional(11) is then solved by evolving the gradient flow ∂φ∂t = ∆ δLδφ ( φ ) = ∆ (cid:18) δFδφ ( φ ) (cid:19) + 2 (cid:28) δFδ ˆ φ ( ˆ φ ) , v (cid:29) L v, (17)where ˆ φ is defined in (14). By solving (16) and (17), one can get the saddle point of F ( φ ) in H − metric. The readers can refer to [13] for details.For the IMF in the H − metric, we can see in (14) that the H − inner productcalculation requires to get the − ∆ − operator first. This can be transformed to aPoisson equation ∆ w = − u and it takes large computational cost. This is why weconsider the projected method in this note to locate the saddle point in H − metric.In the next two subsections, we first present the projected IMF and then propose theprojected GAD.3.2. The Projected IMF.
In this part, we propose the projected iterative minimiza-tion formulation to locate the saddle point of an energy functional F ( φ ) in the H − metric. Since the mass is preserved in H − metric, we introduce the projection PP u := u − | Ω | Z Ω u ( x ) dx (18)onto the linear subspace H = (cid:8) u ∈ L : R Ω u ( x ) dx = 0 (cid:9) . One can show that P hasthe following properties:(1) P = P ;(2) P u ∈ H , ∀ u ∈ L ;(3) P v = v, ∀ v ∈ H ;(4) h v, P w i L = h P v, w i L , ∀ v ∈ L and ∀ w ∈ H .In fact, P u = P ( u − | Ω | Z Ω u ( x ) dx )= u − | Ω | Z Ω u ( x ) dx − | Ω | Z Ω ( u ( x ) − | Ω | Z Ω u ( x ) dx ) dx = P u, ∀ u. Besides, one can easily show P u ∈ H , since Z Ω P u dx = Z Ω ( u − | Ω | Z Ω u ( x ) dx ) dx = 0 , ∀ u. (19)The third property is obvious. For the last one, when w ∈ H , v ∈ L , we have h v, P w i L − h P v, w i L = h v, w i L − h P v, w i L = h v − P v, w i L . Since P v ∈ H , v ∈ L and P is the projection from L to H , we obtain v − P v ∈ H ⊥ ⇒ h v − P v, w i L = 0 , that is, h v, P w i L = h P v, w i L . For the rotation step (10) in the IMF, the following equivalence can be get D v, e H v E H − = h v, − ∆ H v i H − = (cid:10) v, − ( − ∆) − ∆ H v (cid:11) L = h v, H v i L . When v ∈ H , by the last two properties of P , we have h v, H v i L = h P v, HP v i L = h v, PHP v i L , so the eigenvector problem (10) can be equivalent transformed to v ( k +1) = argmin h v, PHP v i L without regarding to the “length” (norm) of v since only the direction of v matters inthe translation step.Besides, due to the equivalence of saddle points of F ( φ ) in H − metric (Cahn-Hilliardequation) and in L metric with projection (projected Allen-Cahn equation), all theterms in equation (10) and (11) including the (first-order and second-order) variation,the inner-product and the norm in H − metric can be changed to the L metric ontothe confined subspace H . It can be verified that the relation of the variations in L space and its subspace H is µ = P µ , b H = PHP , where µ and µ are the first-order variations of F ( φ ) in H and L , respectively. b H and H are the second-order variations of F ( φ ) in H and L , respectively.So the projected IMF written in terms of h· , ·i L is v ( k +1) = argmin k v k L =1 D v, b H ( φ ( k ) ) v E L , (20) φ ( k +1) = argmin R Ω φ ( x ) dx = m L ( φ ; φ ( k ) , v ( k +1) ) , (21)where b H = PHP , H = δ φ F ( φ ). L ( φ ) = F ( φ ) − F ( ˆ φ ), withˆ φ = φ ( k ) + D v ( k +1) , ( φ − φ ( k ) ) E L v ( k +1) . (22)The eigenvector problem (20) is equivalent to ( PHP ψ = λψ, R Ω ψ ( x ) dx = 0 , (23)subject to some boundary condition. In this paper, we consider the periodic boundarycondition only. The Rayleigh quotient in this case is b R ( ψ ) = h ψ, PHP ψ i L k ψ k L , and thus the min-mode of b H is the minimizer of the problemargmin ψ (cid:26) b R ( ψ ) : Z Ω ψ dx = 0 , k ψ k L = 1 (cid:27) . (24)After the min-mode is obtained, the subproblem of minimizing the auxiliary functional(21) is then solved by evolving the gradient flow ∂φ∂t = − P δLδφ ( φ ) = − P δFδφ ( φ ) + 2 (cid:28) v, δFδ ˆ φ ( ˆ φ ) (cid:29) L P v, (25) where ˆ φ = φ ( k ) + D v ( k +1) , ( φ − φ ( k ) ) E L v ( k +1) . We can show that (25) ensures mass conservation automatically ∂∂t Z Ω φ ( x ) dx = − Z Ω P δLδφ dx = 0 , thanks to (19). Thus, by solving (24) and (25), we can get the saddle point of F ( φ ) in H − metric. It is easy to find that (25) is two order lower in spatial derivative than(17), which is the gradient flow in H − metric directly. Besides, the inner product in(25) is in L metric which avoids the ∆ − operator calculation. We can also apply theconvex splitting method to (25) to construct a large time step size scheme as in [13].Denote this mapping for the iteration as Φ( φ ), we shall show that the Jacobianmatrix of Φ( φ ) in the projection sense vanishes at the index-1 saddle point. Thisimplies that the projected IMF is of quadratic convergence rate.3.3. Convergence Results.Theorem 3.
Suppose that φ ∗ is a (non-degenerate) index-1 saddle point of the func-tional F ( φ ) , which satisfies that the second order variational derivative δ φ F ( φ ) is con-tinuous. For each φ , v ( φ ) is the normalized eigenvector corresponding to the smallesteigenvalue of the matrix b H = PHP , H = δ φ F , i.e., v ( φ ) = argmin k u k =1 u T b H ( φ ) u. Take α = 0 , β = 2 , and the auxiliary functional is L ( φ ) = F ( φ ) − F ( ˆ φ ) , then (1) φ ∗ is local minimizer of L ( φ ; φ ∗ , v ) ; (2) a neighbourhood U of φ ∗ exists such that for any φ ∈ U , L ( φ ; φ ( k ) , v ) is strictlyconvex in φ ∈ U and thus has a unique minimum in U ; (3) define the mapping Φ : φ ∈ U → Φ( φ ) ∈ U to be the unique minimizer of L in U for any φ ∈ U . Further assume that U contains no other stationary points of F exceptfor φ ∗ . Then the mapping Φ has only one fixed point φ ∗ ; (4) Φ( φ ) is differentiable in U and P Φ ′ ( φ ∗ ) P = 0 . Thus the mapping φ → Φ( φ ) has alocal quadratic convergence rate.Proof. The proof of the first three conclusions can be generalized from the finite spaceto the infinite space based on the proof of Theorem 3.1 in [11] without difficuty. Themain difference is the proof of the quadratic convergence rate. Here, we only give thedetails of the final conclusion.In fact, the first order variational derivative of L ( φ ; φ ( k ) , v ( φ )) can be calculated as δ φ L ( φ ; φ ( k ) , v ( φ )) = δ φ F ( φ ) − D v, δ ˆ φ F ( ˆ φ ) E L v. At each φ ( k ) ∈ U , the mapping Φ( φ ( k ) ) satisfies the first order equation P δ φ L (Φ( φ ( k ) ) , φ ( k ) , v ( φ ( k ) )) =0, that is, P δ φ F (Φ( φ ( k ) )) − D v ( φ ( k ) ) , δ ˆ φ F ( ˆ φ ) E L P v ( φ ( k ) ) = 0 . (26) Take derivative w.r.t. φ ( k ) on both sides of (26), we get P H (Φ( φ ( k ) )) P Φ ′ ( φ ( k ) ) P − D v ( φ ( k ) ) , δ ˆ φ F ( ˆ φ ) E L P J ( φ ( k ) ) P − D P J ( φ ( k ) ) P , δ ˆ φ F ( ˆ φ ) E L P v ( φ ( k ) ) − D v ( φ ( k ) ) , PH ( ˆ φ ) P ˆ φ ′ E L P v ( φ ( k ) ) = 0 , (27)where H = δ φ F, J ( φ ( k ) ) = ∂v ( φ ( k ) ) ∂φ ( k ) andˆ φ ′ ( φ ( k ) ) = P + D v ( φ ( k ) ) , Φ( φ ( k ) ) − φ ( k ) E L P J ( φ ( k ) ) P + D P J ( φ ( k ) ) P , Φ( φ ( k ) ) − φ ( k ) E L v ( φ ( k ) )+ D v ( φ ( k ) ) , P Φ ′ ( φ ( k ) ) P − P E L v ( φ ( k ) ) . Let φ ( k ) = φ ∗ be the saddle point, we have Φ( φ ∗ ) = φ ∗ , ˆ φ = φ ∗ , δ ˆ φ F ( φ ∗ ) = 0 andˆ φ ′ = P + (cid:10) v ( φ ( k ) ) , P Φ ′ ( φ ∗ ) P − P (cid:11) L v ( φ ∗ ) , thus (27) becomes PH ( φ ∗ ) P Φ ′ ( φ ∗ ) P = 2 h v ( φ ∗ ) , PH ( φ ∗ ) P [ P + h v ( φ ∗ ) , P Φ ′ ( φ ∗ ) P − P i L v ( φ ∗ )] i L P v ( φ ∗ ) , which can be simplified as( PH ( φ ∗ ) P − λvv T ) P Φ ′ ( φ ∗ ) P = 0 , (28)by denoting u T v = h u, v i , applying PHP v = λv and PHP ( I − vv T ) P = 0. (28) impliesthat P Φ ′ ( φ ∗ ) P = 0 . One can carry out the second order derivative of Φ( φ ) at φ ∗ further and observe that P Φ ′′ ( φ ∗ ) P = 0 does not trivially hold. Thus the iteration φ → Φ( φ ) locally convergesto φ ∗ with the quadratic rate. (cid:3) Remark 2.
Theorem 3 is also applicable for any auxiliary functional L only if α + β > . Here we take α = 0 , β = 2 just for convenience. α and β are defined in Section 2.1. Remark 3.
We are dealing with a linear constraint here, so the projection P is thestandard orthogonal projection. For general nonlinear constraints giving rise to a sub-manifold, the projection should follow the geodesic distance on the submanifold exactlyto ensure the quadratic convergence rate in IMF [12] . Projected GAD.
We now present the projected GAD to calculate the saddlepoint of the energy functional F ( φ ) in H − metric. For comparison, we put the originalGAD (9) in Section 2.2 here: ∂φ∂t = − δ φ F ( φ ) + 2 h δ φ F ( φ ) , v i L h v, v i L v, (29) γ ∂v∂t = − δ φ F ( φ ) v + (cid:10) v, δ φ F ( φ ) v (cid:11) L v. (30)By using the projection P in (18), the projected GAD is given as follows: ∂φ∂t = − P δ φ F ( φ ) + 2 h δ φ F ( φ ) , v i L h v, v i L P v, (31) γ ∂v∂t = − P δ φ F ( φ ) P v + (cid:10) v, P δ φ F ( φ ) P v (cid:11) L v. (32)By integrating w.r.t. x on both sides of (32), we get γ ∂∂t Z Ω v dx = − Z Ω P δ φ F ( φ ) P v dx + (cid:10) v, P δ φ F ( φ ) P v (cid:11) L Z Ω v dx = (cid:10) v, P δ φ F ( φ ) P v (cid:11) L Z Ω v dx, this is an ordinary differential equation of R Ω v dx . Considering the initial condition, R Ω v dx = 0, one can easily get Z Ω v ( x ) dx = 0 , ∀ v, and thus P v = v − Z Ω v ( x ) dx = v, ∀ v. Furthermore, by using the last property of P , h v, P w i L = h P v, w i L , ∀ v ∈ L , ∀ w ∈ H . the projected GAD can be rewritten as ∂φ∂t = − P δ φ F ( φ ) + 2 h δ φ F ( φ ) , v i L h v, v i L P v, (33) γ ∂v∂t = − P δ φ F ( φ ) P v + (cid:10) v, δ φ F ( φ ) v (cid:11) L v. (34)By solving the equation (33) and (34), we can calculate the saddle points of F ( φ ) in H − metric. 4. Numerical example
In this section, we will illustrate the above projection method by locating the transi-tion state of the one dimensional Ginzburg-Landau free energy and the two dimensionalLandau-Brazovskii free energy in the H − metric.4.1.
1D example: Ginzburg-Landau free energy.
Consider the one dimensionalGinzburg-Landau free energy on [0 , F ( φ ) = Z h κ ∂φ∂x ) + f ( φ ) i dx, (35)where φ ( x ) is an order parameter and κ > f ( φ ) = ( φ − /
4. The first and thesecond order variation of F ( φ ) can be calculated as δFδφ ( φ ) = − κ ∆ φ + φ − φ,δ Fδφ ( φ ) = − κ ∆ + 3 φ − H . So the projected IMF is v ( k +1) = argmin k v k =1 D v, PHP ( φ ( k ) ) v E L , (36) φ ( k +1) = argmin R Ω φ ( x ) dx = m L ( φ ; φ ( k ) , v ( k +1) ) , (37)with L ( φ ) = F ( φ ) − F ( ˆ φ ), ˆ φ is defined in (22). The second minimization sub-problem(37) is solved by evolving the gradient flow: ∂φ∂t = − P δ φ L ( φ ) , where − P δ φ L ( φ ) = − P (cid:2) − κ ∆ φ + ( φ − φ ) (cid:3) +2 D v, − κ ∆ ˆ φ + ( ˆ φ − ˆ φ ) E L P v, here φ = φ ( k +1) . And the projected GAD is ∂φ∂t = − P ( − κ ∆ φ + ( φ − φ )) + 2 (cid:10) v, − κ ∆ φ + ( φ − φ ) (cid:11) L h v, v i L P v, (38) γ ∂v∂t = − P ( − κ ∆+(3 φ − P v + (cid:10) v, − κ ∆ v +(3 φ − v (cid:11) L v. (39)We apply the finite difference scheme to achieve the numerical example. For the pro-jected IMF, we further construct the following convex splitting scheme (40) to discrete(37) in time. φ n +1 − φ n ∆ t = P (cid:2) κ ∆ φ − φ − h v, φ i v (cid:3) n +1 + P h − φ + 3 φ + 2 D v, − κ ∆ ˆ φ + ˆ φ E v i n . (40)In the numerical test, we take κ = 0 .
04, the initial mass m = 0 .
6, and the meshgrid is { x i = ih, i = 0 , , , . . . , N } . h = 1 /N. N = 100, ∆ t = 0 .
1. We use the periodicboundary condition in this example. We find that the saddle point of F ( φ ) in the H − metric calculated by the projected IMF or the projected GAD is exactly the same asthe result in [13] which applies the IMF in the H − metric directly, see Figure 1a.Besides, the quadratic convergence rate can also be observed when using the projectedIMF; see Figure 1b for the convergence result. In order to illustrate the advantageof this method, we make comparison of the CPU time required for the same iterationnumber between the projected IMF in L metric and the original IMF in H − metric.Table 1 shows the results for various initial states φ , φ and φ . One can find thatthe projected IMF in L metric can save almost half computational cost compared withthe IMF in H − metric, especially for the large inner iteration number.4.2.
2D example: Landau-Brazovskii free energy.
In this section, we study thenucleation problem of phase transition in diblock copolymers [17, 23], which have at-tracted a lot attention because of their various and abundant microstructures. Themodel is described by the two-dimensional Landau-Brazovskii energy functional of theorder parameter φ , F ( φ ) = Z Ω ξ φ ( r )] + Φ( φ ) d r , (41) x φ ( x ) initial stateProj IMF/GADIMF in H −1 (a) transition states −15 −10 −5 outer cycle k k P δ φ F ( φ ( k ) ) k L (b) convergence rate Figure 1. (A): initial state (dashed line); transition state by pro-jected IMF or projected GAD (green line); transition state by IMF in H − metric (red line). (B): The decay of the error k P δ φ F ( φ ( k ) ) k L measured by the L norm of the projected force at each cycle k .iterN 1 e e e e e φ IMF 16.11 32.06 80.52 158.75 320.99Projected IMF 9.14 18.14 45.94 91.07 182.68 φ IMF 15.73 32.42 80.55 158.87 316.28Projected IMF 9.30 18.26 45.27 90.75 179.99 φ IMF 16.02 33.21 80.24 160.02 325.60Projected IMF 9.17 18.43 45.71 91.13 183.00
Table 1.
CPU time (seconds) comparison. “IMF” means the originalIMF in H − metric; “Projected IMF” is in L metric.defined on Ω = [0 , π √ ] × [0 , π ], where Φ( φ ) = τ φ − γ φ + φ . The parameters are τ = − . , ξ = 1 . , γ = 0 .
25. We hope to calculate the transition state of F ( φ ) in the H − metric. First we calculate the first and the second order variations as follows δ φ F ( φ ) = ξ (∆ + 1) φ ( r ) + Φ ′ ( φ ) ,δ φ F ( φ ) = ξ (∆ + 1) + Φ ′′ ( φ ) := H , where Φ ′ ( φ ) = τ φ − γ φ + φ , Φ ′′ ( φ ) = τ − γφ + φ . So the projected IMF is v ( k +1) = argmin k v k =1 D v, PHP ( φ ( k ) ) v E L , (42) φ ( k +1) = argmin R Ω φ ( x ) dx = m L ( φ ; φ ( k ) , v ( k +1) ) , (43)with L ( φ ) = F ( φ ) − F ( ˆ φ ), ˆ φ is defined in (22). The second minimization sub-problemis solved by evolving the gradient flow: ∂φ∂t = − P δ φ L ( φ ) , where − P δ φ L ( φ ) = − P (cid:2) ξ (∆ + 1) φ + Φ ′ ( φ ) (cid:3) + 2 D v, ξ (∆ + 1) ˆ φ + Φ ′ ( ˆ φ ) E L P v, here φ = φ ( k +1) . And the projected GAD is ∂φ∂t = − P (cid:2) ξ (∆ + 1) φ + Φ ′ ( φ ) (cid:3) + 2 (cid:10) ξ (∆ + 1) φ + Φ ′ ( φ ) , v (cid:11) L h v, v i L P v, (44) γ ∂v∂t = − P (cid:2) ξ (∆ + 1) + Φ ′′ ( φ ) (cid:3) P v + (cid:10) v, ξ (∆ + 1) v + Φ ′′ ( φ ) v (cid:11) L v. (45)For this two-dimensional numerical example, we consider the periodic boundarycondition again. And for the convenience of saving computational cost, we apply thefast Fourier transform (FFT) for the two-dimensional case. We take the mesh points N x = N y = 64 and the time step size ∆ t = 0 .
1. The transition state can be obtainedby the projected IMF (or the projected GAD) in Figure 2a. The quadratic convergencerate can also be obtained for the projected IMF shown in Figure 2b. Similarly to theone-dimensional case, in order to illustrate the effect of the projected method, wemake comparison with the original IMF in H − metric. We fix various inner iterationnumber for both cases and compare the required CPU time. Table 2 shows the CPUtime comparison between the projected IMF and the original IMF in H − metric withvarious initial states. Results show that the projected method for this example cansave almost one-third computational cost. -1-0.500.511.5 (a) transition state -8 -6 -4 -2 (b) convergence rate Figure 2. (A): Transition state of the Landau Brazovskii free energyin the H − metric by the projected IMF. (B): The decay of the error k P δ φ F ( φ ( k ) ) k L measured by the L norm of the projected force ateach cycle k . 5. Conclusion
In this work, we present the projected method for the IMF and the GAD to calculatethe transition states of some energy functional in the H − metric. By introducing anorthogonal projection operator onto the confined subspace satisfying the mass conser-vation, the saddle points in H − metric calculation can be transformed equivalently iterN 5 e e e e e e φ IMF 12.62 15.27 17.67 20.35 22.71 25.84Projected IMF 9.03 11.05 12.71 14.07 16.06 17.62 φ IMF 12.71 15.56 17.93 20.67 22.97 25.32Projected IMF 9.05 10.71 12.67 14.47 16.34 17.76 φ IMF 12.77 15. 82 17.39 19.90 22.01 25.51Projected IMF 9.07 11.15 12.95 14.56 16.19 18.07
Table 2.
CPU time (seconds) comparison. “IMF” means the originalIMF in H − metric; “Projected IMF” is in L metric.to the saddle points in L metric with projection. This method can reduce muchcomputational cost. Since it leads to a lower order spatial derivative equation for thetranslation step in the IMF compared with that in H − metric directly; more impor-tantly, it avoids the ∆ − operator calculation. The same phenomenon can be obtainedfor the projected GAD. This projected method maintains the same convergence speedof the original GAD and IMF, but the new algorithm is much faster than the directmethod for H − problem. Acknowledgement
SG acknowledges the support of NSFC 11901211, the youth innovative talent projectof Guangdong province 2018KQNCX055 and the young teacher scientific research culti-vation fund of South China Normal University 18KJ17. LL acknowledges the supportof NSFC 11871486. XZ acknowledges the support of Hong Kong RGC GRF grants11337216 and 11305318. References [1]
S. M. Allen and J. W. Cahn , A microscopic theory for antiphase boundary motion and itsapplication to antiphase domain coarsening , Acta Metallurgica, 27 (1979), pp. 1085 – 1095.[2]
P. W. Bates and F. Chen , Spectral analysis and multidimensional stability of traveling wavesfor nonlocal allen–cahn equation , Journal of Mathematical Analysis and Applications, 273 (2002),pp. 45–57.[3]
M. Brachet and J. Chehab , Fast and stable schemes for phase fields models , (2020).[4]
J. W. Cahn and J. E. Hilliard , Free energy of a nonuniform system. i. interfacial free energy ,The Journal of Chemical Physics, 28 (1958), pp. 258–267.[5]
G. M. Crippen and H. A. Scheraga , Minimization of polypeptide energy : XI. the method ofgentlest ascent , Arch. Biochem. Biophys., 144 (1971), pp. 462–466.[6]
F. P. C. Dang H and P. L. A , Saddle solutions of the bistable diffusion equation , Ztschrift F¨urAngewandte Mathematik Und Physik Zamp, 43 (1992), pp. 984–998.[7]
D. A. Dawson and J. G¨artner , Large deviations, free energy functional and quasi-potential fora mean field model of interacting diffusions , vol. 78, Memoirs of American Mathematical Society,1989.[8]
Q. Du and L. Zhang , A constrained string method and its numerical analysis , Commun. Math.Sci., 7 (2009), pp. 1039–1051.[9]
W. E, W. Ren, and E. Vanden-Eijnden , String method for the study of rare events , Phys. Rev.B, 66 (2002), p. 052301.[10]
W. E and X. Zhou , The gentlest ascent dynamics , Nonlinearity, 24 (2011), p. 1831.[11]
W. Gao, J. Leng, and X. Zhou , An iterative minimization formulation for saddle point search ,SIAM J. Numer. Anal., 53 (2015), pp. 1786–1805.[12] ,
Iterative minimization algorithm for efficient calculations of transition states , Journal ofComputational Physics, 309 (2016), pp. 69 – 87.[13]
S. Gu and X. Zhou , Convex splitting method for the calculation of transition states of energyfunctional , J. Comput. Phys., 353 (2018), pp. 417–434.[14]
G. Henkelman and H. J´onsson , A dimer method for finding saddle points on high dimensionalpotential surfaces using only first derivatives , J. Chem. Phys., 111 (1999), pp. 7010–7022.[15]
H. J`onsson, G. Mills, and K. W. Jacobsen , Nudged elasic band method for finding min-imum energy paths of transitions , in Classical and Quantum Dynamics in Condensed PhaseSimulations, B. J. Berne, G. Ciccotti, and D. F. Coker, eds., New Jersey, 1998, LERICI, VillaMarigola,Proceedings of the International School of Physics, World Scientific, p. 385.[16]
T. Li, X. Li, and X. Zhou , Finding transition pathways on manifolds , Multiscale Model Simul.,14 (2016), pp. 173–206.[17]
T. Li, P. Zhang, and W. Zhang , Nucleation rate calculation for the phase transition of diblockcopolymers undr stochastic cahn-hilliard dynamics , Multiscale Model. Simul., 11 (2013), pp. 385–409.[18]
L. Lin, X. Cheng, W. E, A. Shi, and P. Zhang , A numerical method for the study of nucleationof ordered phases , J. Comput. Phys., 229 (2010), pp. 1797–1809.[19]
N. Mousseau and G. Barkema , Traveling through potential energy surfaces of disordered ma-terials: the activation-relaxation technique , Phys. Rev. E, 57 (1998), p. 2419.[20]
W. Ren and E. Vanden-Eijnden , A climbing string method for saddle point search , J. Chem.Phys., 138 (2013), p. 134105.[21]
J. Shen and X. Yang , Numerical approximations of allen-cahn and cahn-hilliard equations ,Discrete and Continuous Dynamical Systems, 28 (2010), pp. 1669–1691.[22]
D. J. Wales , Energy Landscapes with Application to Clusters, Biomolecules and Glasses , Cam-bridge University Press, 2003.[23]
S. M. Wise, C. Wang, and J. Lowengrub , An energy-stable and convergent finite-differencescheme for the phase field crystal equation , SIAM J. Numer. Anal., 47 (2009), pp. 2269–2288.[24]
J. Zhang and Q. Du , Constrained shrinking dimer dynamics for saddle point search with con-straints , J. Comput. Phys., 231 (2012), pp. 4745–4758.[25]