View source for M04A Conjugate Gradient Method in MD++

<H1 ALIGN="CENTER">
Conjugate Gradient Method in MD++</H1>
<DIV>

<P ALIGN="CENTER"><STRONG>Keonwook Kang and Wei Cai</STRONG></P>
</DIV>
<P ALIGN="CENTER">[[ Mar 19 ]], [[2009]]</P>
<P>
<BR><HR>

== Conjugate Gradient Method ==
The conjugate gradient (CG) method was initially proposed as an iterative method to solve a large linear system of equations,<ref>M. R. Hestenes and E. Stiefel, "Methods of conjugate gradients for solving linear systems", ''Journal of Research of the National Bureau of Standards'', '''49'''(1952) pp.409-436</ref> such as <math>\mathrm{A}\mathbf{x} = \mathbf{b}</math>, where A is an ''n'' by ''n'' symmetric (positive-definite) matrix, and '''b''' is an ''n'' by 1 vector. The basic idea is that the ''quadratic form'' of a function ''&Psi;''('''x'''),
{|border="0" align="center"
|<math> \Psi(\mathbf{x}) = \frac{1}{2}\mathbf{x}^T \mathrm{A}\mathbf{x} - \mathbf{b}^T\mathbf{x} + c</math> 
|}
is minimized by the solution of <math>\mathrm{A}\mathbf{x} = \mathbf{b}</math> if A is symmetric and positive-definite.
Later, Fletcher and Reeves<ref>R. Fletcher and C. M. Reeves, "Function minimization by conjugate gradients", ''Computer Journal'', ''7''(1964) pp.149-154</ref> adopted it to solve nonlinear optimization problems. In MD++, the CG method is used to find a meta-stable structure by minimizing a given (usu. nonlinear) potential energy function. 

The algorithm of the CG method is structured by two nested loops. The outer loop tries different search directions, and the inner loop aims to find the optimal step length which minimizes the object function along a given search direction. The inner loop may not be needed when the exact expression of the optimal step length can be given. If a certain condition is met, the program exits out of the loops and stops. The algorithm is shown in figure 1. 
[[Image:Logic_flow.jpg | frame | 100px | '''Figure 1''' Algorithm of the (nonlinear) CG method.]]
In figure 1, ''&Psi;''('''x''') is a function to be minimized, '''g''' is the gradient of ''&Psi;'', and '''s''' is the search direction. ''&alpha;''<sub>k</sub> is the optimal step length that minimizes the function ''&Psi;''('''x'''<sub>k</sub>+''&alpha;''<sub>k</sub>'''s'''<sub>k</sub>) at the ''k''-th iteration. ''&beta;''<sub>k+1</sub> is the parameter related with update of the search direction. 

Depending on how we choose ''&beta;'', the minimization methods can be different. For example, if ''&beta;''<sub>k</sub>=0, it is the steepest descent (SD) method. 

While there exists the analytical expression of ''&alpha;'' for a linear system, given as 

{|border="0" align="center"
| rowspan="2"| <math>\alpha_k=\Bigg\{</math>
|<math> \frac{\mathbf{r}^T_k\mathbf{r}_k}{\mathbf{r}^T_k\mathrm{A}\mathbf{r}_k}</math> 
| for SD,<ref>This expression can be obtained by the fact that &alpha; minimizes &Psi; when the directional derivative d&Psi;(''x''<sub>k+1</sub>)/d&alpha; is zero, which is effectively same with that two consecutive gradients (or residuals) are orthogonal, or '''g'''<sup>T</sup><sub>k+1</sub>'''g'''<sub>k</sub>=0. In this case, we don't need the inner loop.</ref>
|-
|<math>\frac{\mathbf{r}^T_k\mathbf{r}_k}{\mathbf{s}^T_k\mathrm{A}\mathbf{s}_k} </math>
| for CG,
|}

where '''r'''<sub>k</sub> is the residual and is defined as <math>\mathbf{r}_k \equiv \mathbf{b} - \mathrm{A}\mathbf{x}_k=-\mathbf{g}_k</math><ref>The residual is a transformed error, or ''A''(''x''-''x''<sub>k</sub>) = ''A'' ''e''<sub>k</sub>. At the same time the residual is the steepest descent direction, or ''r''<sub>k</sub> = -''g''(''x''<sub>k</sub>). As we approach the solution, the residual becomes zero.</ref>, for nonlinear problems ''&alpha;'' is generally obtained by a numerical method such as Newton's method, quasi-Newton method, secant method, interpolation method, or other similiar methods. These methods of finding an optimal ''&alpha;'' are called as the '''line search methods'''. For details of the line search method, refer Nocedal and Wright's '''Numerial Optimization'''<ref name="Nocedal">Jorge Nocedal and Stephen J. Wright, ''Numerical Optimization'' 2nd ed. Springer (2006)</ref>

In the CG method, the search directions are constructed from the residuals such that any two directions are to be ''A-orthogoanl'' or ''conjugate''.<ref>Two vectors '''d'''<sub>''i''</sub> and '''d'''<sub>''j''</sub> are A-orthogonal or conjugate if <math>\mathbf{d}_i^TA\mathbf{d}_j = 0</math> for different ''i'' and ''j''.</ref> For example, to construct a ''k''-th search direction '''s'''<sub>''k''</sub>, we take '''r'''<sub>''k''</sub> from which we subtract any non A-orthogonal components to all the previous search directions. In other words, let's express '''s'''<sub>''k''</sub> as
{|border="0" align="center"
|<math> \begin{array}{rcl} \mathbf{s}_k & = & \mathbf{r}_k + \sum_{i=0}^{k-1}\beta_{ki}\mathbf{s}_i \\  & = & -\mathbf{g}_k + \sum_{i=0}^{k-1}\beta_{ki}\mathbf{s}_i \end{array} </math>; 
|-
|&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <math>\mathbf{s}_0 = \mathbf{r}_0 = -\mathbf{g}_0</math>.
|}
and determine &beta;<sub>''ki''</sub> to satisfy the conjugacy, <math>\mathbf{s}_k^TA\mathbf{s}_i= 0</math> for different ''k'' and ''i''.
<!--{|border="0" align="center"
|<math>\mathbf{r}_i^T \mathbf{r}_j = 0</math>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; for &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<math>i \neq j</math>.
|}-->
From this condition, ''&beta;''<sub>k</sub> is obtained as

{|border="0" align="center"
|<math>\beta_{k+1} \equiv \frac{\mathbf{r}_{k+1}^T\mathbf{r}_{k+1}}{\mathbf{r}_{k}^T\mathbf{r}_{k}}</math>.
|}

For details of the derivation, refer J. R. Shewchuk's.<ref>J. R. Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (1994)</ref> The benefit of this choice is that a new search direction is A-orthogonal to all the previous search directions, or ''s''<sub>i</sub><sup>T</sup>''A'' ''s''<sub>j</sub> = 0 if i is not equal to j. In other words, we don't need to store all the previous search directions to produce a linearly independent new search direction. 

=== Non-linear CG method ===

For non-linear CG method, there exist different choices of ''&beta;'' as shown in table 1, in which the residual '''r'''<sub>k</sub> is now defined as the negative gradient or <math>\mathbf{r}_k = -\mathbf{g}_k</math>.

{| class="wikitable" style="text-align:center"
|+ Table 1. The update parameter ''&beta;'' in non-linear CG. For more different choices of ''&beta;'', refer Hager and Zhang's.<ref>W. W. Hager and H. Zhang, "A survey of nonlinear conjugate gradient methods", ''Pacific J. of Optimization'', '''2'''(2006) pp.35-58</ref>
|-
! Name !! <math>\beta</math> !! Reference !! Features
|-
! HS
| <math>\beta^{HS}_{k+1}= \frac{\mathbf{g}_{k+1}^T\mathbf{g}_{k+1} - \mathbf{g}_{k+1}^T\mathbf{g}_{k}} {(\mathbf{g}_{k+1} - \mathbf{g}_{k})^T\mathbf{s}_{k}}</math> <ref>For detailed derivation, refer Nocedal and Wright<ref name="Nocedal"/></ref> || Hestenes and Stiefel (1952)|| 
|-
! FR
| <math>\beta^{FR}_{k+1}= \frac{\mathbf{r}_{k+1}^T\mathbf{r}_{k+1}}{\mathbf{r}_{k}^T\mathbf{r}_{k}} = \frac{\mathbf{g}_{k+1}^T\mathbf{g}_{k+1}}{\mathbf{g}_{k}^T\mathbf{g}_{k}}</math>  || Fletcher and Reeves (1964)|| Considered as the 1st non-linear CG method. Same with ''&beta;'' in linear CG. <font color="red">Jamming problem</font><ref>While you minimize a general object function, you may have a bad search direction, which is almost orthogonal to the steepest descent direction, -'''g'''<sub>k</sub>. Once you have a bad direction, the FR method is likely to give you a next search direction close to the previous direction, or <math>\mathbf{s}_{k+1}\sim \mathbf{s}_{k}</math>, which is another bad direction. See Nocedal and Wright pp. 125 - 127 Behavior of the Fletcher-Reeves Method</ref>
|-
! PR
| <math>\beta^{PR}_{k+1}= \frac{\mathbf{r}_{k+1}^T(\mathbf{r}_{k+1}-\mathbf{r}_{k})}{\mathbf{r}_{k}^T\mathbf{r}_{k}} = \frac{\mathbf{g}_{k+1}^T(\mathbf{g}_{k+1}-\mathbf{g}_{k})}{\mathbf{g}_{k}^T\mathbf{g}_{k}}</math>  || Polak and Ribiere (1969), and Polyak (1969)|| If an object function &Psi; is strongly convex quadratic and the line search is exact, ''&beta;''<sup>PR</sup> = ''&beta;''<sup>FR</sup> because '''r'''<sub>''i''</sub><sup>T</sup>'''r'''<sub>''j''</sub>=0 for <math>i \neq j</math>. Can avoid jamming problem by effectively restarting CG with <math>\beta \simeq 0</math>. <ref>When a search direction is almost orthogonal to the steepest descent direction, ''&beta;''<sup>PR</sup> becomes almost zero (<math>\beta \sim 0</math>) and the search direction is close to the steepest descent direction.</ref>
|-
! PR+
| <math>\beta^{PR+}_{k+1}= \max\{\beta^{PR}_{k+1},0\}</math>  || Powell (1984)<ref>M. J. D. Powell, "Nonconvex minimization calculations and the conjugate gradient method", '''Numerical Ananlysis''' (Dundee, 1983), Lecture Notes in Mathematics, vol. 1066, Springer-Verlag, Berlin, 1984 pp.122-141</ref>|| rowspan=2 | Convergence guaranteed. <font color="blue">Restart CG if ''&beta;''<sub>k+1</sub><0.</font> Restarting CG means that all the previous search directions are to be forgotten and the new search direction is set to be the steepest descent direction.
|-
! HS+
| <math>\beta^{HS+}_{k+1}= \max\{\beta^{HS}_{k+1},0\}</math>  || 
|}

=== Hybrid method ===

For better performance, people may use a hybrid method, which can select different ''&beta;'' at each iteration.   One example of hybrid methods is PR-FR method, in which method ''&beta;'' is selected as

{|border="0" align="center"
|<math> \beta_k = \left\{ \begin{array}{rll} -\beta_k^{FR} & \mathrm{if} & \beta_k^{PR}<-\beta_k^{FR} \\ \beta_k^{PR} & \mathrm{if} & |\beta_k^{PR}| \leq \beta_k^{FR}\\ \beta_k^{FR} & \mathrm{if} &\beta_k^{PR}<\beta_k^{FR} \end{array}. \right. </math> 
|}

=== Stopping conditions ===

Below are listed the conditions to be used for stopping the CG program. You can terminate the CG routine if any one of the conditions is satisfied.

# Energy criterion
#;Exit the loop if the following inequality holds for a given energy tolerance, ''E''<sub>tol</sub>, which is in no unit. 
#;::::::<math>\frac{|E_{k+1} - E_k|}{|E_k|} < E_{\mathrm{tol}}</math> 
#Force criterion (''f''<sub>tol</sub>)
#;Exit the loop if the following inequality holds for a given force tolerance, ''f''<sub>tol</sub> (eV/&Aring;).
#;::::::<math> \mathbf{g}_{k+1}^T\mathbf{g}_{k+1} < f_{\mathrm{tol}}^2</math> 
#Max. evaluations
#;Exit the loop if the number of function-evaluation exceeds the maximum calling number, max<sub>call</sub>.
# Max. iterations
#;The loop will finish if the iteration number reaches the maximum value, max<sub>iter</sub>.

Conditions 1 to 3 may also be checked inside the inner loop for earlier exit and CG termination. Other condtions to be checked are
* the new search direction is uphill ('''s'''<sub>k+1</sub><sup>T</sup>'''g'''<sub>k+1</sub> >= 0) or downhill ('''s'''<sub>k+1</sub><sup>T</sup>'''g'''<sub>k+1</sub> < 0), and
* the step size &alpha; is zero or not.
If the new search direction is uphill, the search direction is set to be a steepest descent direction. If the step size &alpha; is zero, then stop the CG algorithm.

<!-- <math> x_{k+1} = x_k + \alpha_k s_k  </math>.
where <math>\alpha_k</math> and <math>s_k</math> are the step length and the search direction at k-th iteration, respectively, and ''k'' is the index of the outer loop.-->
<HR>

== Implementation of the CG method in MD++ == 

MD++ is equipped with two different CG methods. One is what we call as '''zxcgr''' and the other is '''PR+'''. Details of each method will be explained below. 

=== zxcgr ===
'''zxcgr''' code was a part of [http://www.vni.com/products/imsl/ IMSL(International Mathematics and Statistics Library)]'s [http://www.vni.com/products/imsl/documentation/fort06/math/NetHelp/default.htm?turl=umcgg.htm U2CGG.F]. It was reimplemented into C to be used in MD++.
<!--'''zxcgr''' code was written by Dongyi Liao at MIT in 1999 based on [http://www.vni.com/products/imsl/ IMSL(International Mathematics and Statistics Library)]'s [http://www.vni.com/products/imsl/documentation/fort06/math/NetHelp/default.htm?turl=umcgg.htm U2CGG.F]. -->
This code follows Powell's restart procedures for the CG method.<ref>M. J. D. Powell, "Restart Procedures for the Conjugate Gradient Method", ''Mathematical Programming'' '''12''' (1977) pp.241-254</ref> According to his paper, the search direction is updated as

{|border="0" align="center"
|<math>\mathbf{s}_k = -\mathbf{g}_k + \beta_k\mathbf{s}_{k-1} + \lambda_k\mathbf{s}_t \quad \mathrm{for} \quad k \geq 1 \quad \mathrm{and} \quad \mathbf{s}_0 = -\mathbf{g}_0</math>, 
|}
where <math>\beta_k = \beta_k^{HS}</math> and the extra term <math>\lambda_k\mathbf{s}_t</math> is related with how we restart CG. If we need to restart CG at ''k''-th iteration, we set ''t'' = ''k'' - 1 and &lambda;<sub>k</sub> = 0. Then the restarting direction simply becomes '''s'''<sub>''k''</sub> = -'''g'''<sub>''k''</sub> + &beta;<sub>''k''</sub>'''s'''<sub>''k''-1</sub>. After restarting CG, &lambda;<sub>''k''</sub> is defined as 

{|border="0" align="center"
|<math>\lambda_k = \frac{\mathbf{g}_k^T(\mathbf{g}_{t+1}-\mathbf{g}_t)} {\mathbf{s}_t^T(\mathbf{g}_{t+1}-\mathbf{g}_t)}, </math>
|}
where ''t'' is the iteration number that is one time earlier than the iteration number at which restarting happened. Thus, &lambda;<sub>k</sub> can be summarized as

{|border="0" align="center"
|<math>\lambda_k = \left\{ \begin{array}{ll} 0 & \mathrm{if} \quad k=t+1 \\ \frac{\mathbf{g}_k^T(\mathbf{g}_{t+1}-\mathbf{g}_t)} {\mathbf{s}_t^T(\mathbf{g}_{t+1}-\mathbf{g}_t)} & \mathrm{if} \quad k > t+1. \end{array} \right.</math>
|}
The initial value of ''t'' is zero and will not change until the 1st restart occurs. Naturally, you would think of '''what the restarting conditions are'''. One condition is that a new gradient '''g'''<sub>''k''+1</sub> loses enough orthogonality with the previous gradient '''g'''<sub>''k''</sub>. In other words, the quantities '''g'''<sub>''k''</sub><sup>T</sup>'''g'''<sub>''k''+1</sub> and '''g'''<sub>''k''+1</sub><sup>T</sup>'''g'''<sub>''k''+1</sub> become comparable. If the inequality
{|border="0" align="center"
|<math> \frac{|\mathbf{g}_k^T\mathbf{g}_{k+1}|}{|\mathbf{g}_{k+1}^T\mathbf{g}_{k+1}|} \geq 0.2 </math>
|}
holds, restart CG. The other condition is that a new search direction is not sufficiently downhill. As the search direction gets closer to the steepest descent direction or -'''g'''<sub>''k''+1</sub>, we can say the direction becomes more downhill. So, we can propose an inequality 
{|border="0" align="center"
|<math> -1.2\|\mathbf{g}_{k+1}\|^2 \leq \mathbf{s}_{k+1}^T\mathbf{g}_{k+1} \leq -0.8\|\mathbf{g}_{k+1}\|^2 </math>
|}
for the sufficient downhill condition. If this condition is not met, restart CG. The sufficient downhill condition can be rewritten as
{|border="0" align="center"
|<math> |\beta_{k+1}(\mathbf{s}_k^T\mathbf{g}_{k+1}) + \lambda_{k+1}(\mathbf{s}_t^T\mathbf{g}_{k+1})| \leq 0.2\|\mathbf{g}_{k+1}\|^2 </math>
|}
using the relation for the search direction update. There is one more restarting condition, <math>(k-t) \geq n</math>, where ''n'' is the number of degree of freedom of the system. This condition is rarely used in atomistic simulations because ''n'' is usually very huge.

Unlike the other CG methods, this variation of CG method provides a restarting direction that is not a steepest descent direction. This modification usually takes less iterations to get to the minimum.

'''The Line Search Method in zxcgr'''

Since the search direction is downhill when the line search starts, there are only three possible cases shown below.
<gallery caption="Figure 2" widths="220px" heights="200px" perrow="3">
Image:Linesearch1.jpg|(a) <math>f(\Delta\alpha)<f(0)</math> and <math>f'(\Delta\alpha)f'(0)>0</math>
Image:Linesearch2.jpg|(b) <math>f(\Delta\alpha)<f(0)</math> and <math>f'(\Delta\alpha)f'(0)<0</math>
Image:Linesearch3.jpg|(c) <math>f(\Delta\alpha)>f(0)</math> 
</gallery>
In case (a), keep searching the optimal ''&alpha;'' in the same direction. you may increase the step size for fast search.
In case (c), reduce the step size by half, and perform the next search. 
In case (b), interpolate ''f''(''&alpha;'') with the quadratic spline and set ''&alpha;''<sub>k</sub> that minimizes the spline.
Note that the slope of a function ''f''(''&alpha;''), d''f''/d''&alpha;'' corresponds to the dot product of the gradient vector at ''&alpha;'' with the search direction, or
{|border="0" align="center"
|<math>\frac{df(\alpha)}{d\alpha} = \mathbf{s}_k^T\mathbf{g}(\mathbf{x}_k + \alpha\mathbf{s}_k)</math>.
|}

In '''zxcgr''', the maximum iteration number for the line search is set as MAXLIN = 5 in the relax_zxcgr.cpp file. If the line search cannot find an ''&alpha;'' that reduces the function ''f'' within 5 consecutive iterations, the following error message will be printed: "CGRelax: Line search aborted"

=== PR+ ===
'''PR+''' routine was added by Keonwook Kang in 2009 based on LAMMPS's min_cg.cpp file. The update parameter is  ''&beta;''<sup>PR+</sup>. For the line search method, the backtracking method<ref>Nocedal and Wright, "Numerical Optimization" ch.3.1 pp.37 sufficient decrease and backtracking</ref> is used. The backtracking method does not guarantee that a function ''f'' minimizes at the obtained ''&alpha;''. This method allows any step size that can sufficiently decrease the function ''f''.
This decision can be justified if it is undesirable or impossible to exactly minimizes a function ''f'' or if we do not want to spend too much computation time to find an exact ''&alpha;''.

<!--<pre></pre>
For more information, see [[media:VFE.pdf‎ | supplementary document: Vacancy formation energy by MD++ ]].-->
<HR>

== Comparision of the CG methods in different implementation ==
In the first part of this section, I compare the computation results obtained from CG relxation between MD++ (ver.Mar202009) and LAMMPS (ver.Mar262009). The test system is a 3[100] x 3[010] x 3[001] Si diamond cubic structure stretched by 10% in all three directions with two missing atoms (of atom ID 0 and 99). The number of atoms in the test system is 8x3<sup>3</sup> - 2 = 214. The potential fuction is MEAM Siz with the cut-off radius ''r''<sub>c</sub>=6(&Aring;). The CG method is implemented by using the Polak-Ribiere method (PR+) and the backtracking method. The energy and the force tolerances are given as etol = 1e-6 (no unit) and ftol = 1e-8 (eV/&Aring;), respectively.
In the 1st column of table 2, the potential energies of the perfect structure from MD++ and LAMMPS are compared, and they are identical. In the 2nd column, the potential energies of the defect structure before relaxation are compared, and they are again same. In the 4th column, the energies of the relaxed structure are listed in different colors depending on the values of ''d''<sub>max</sub> (0.1, 0.01, and 0.001), which are related with the initial step length, and we can say that the both results are effectively identical. For this small test system, MD++ is faster than LAMMPS, because of less evaluation of the potential function. For a larger system, LAMMPS can perform better due to its ability of parallel run.

{| class="wikitable" style="text-align:center"
|+ Table 2. Relaxation of 3[100]x3[010]x3[001] Si structure with a defect. Potential file : MEAM Siz w/ r_c = 6 (A). CG method = (PR+) + (backtracking). Conditions: etol = 1e-6 (no unit), ftol = 1e-8 (eV/&Aring;). Compiler = intel compiler
|-
!   !! ''E''<sub>perf</sub> (eV) !! ''E''<sub>defect</sub> (eV) !! ''d''<sub>max</sub> (&Aring;) !! ''E''<sub>rlx</sub> (eV) !! Iter. No. !! No. of Calls !! CPU time
|-
! rowspan=3 valign="center" | MD++ (PR+)
| -1.00008000022397e+03 || -9.00066105589930e+02 || 0.1  || style="color:red"|-9.04000029718829e+02 || 15  || 46  || 0.30002
|- 
||   ||   || 0.01 || style="color:blue"|-9.04000167032814e+02 || 46  || 51  || 0.33202
|-
||   ||   || 0.001||style="color:green"|-9.03987772936424e+02  || 424 || 425 || 2.2201
|-
! rowspan=3 valign="center" | LAMMPS (Serial Run)
| -1.00008000022397e+03 || -9.00066105589930e+02 || 0.1 ||style="color:red"|-9.04000029718830e+02 || 15 || 60 || 0.75296
|-
||    ||   ||0.01|| style="color:blue"| -9.04000167032815e+02 || 46 || 96 || 1.2149
|-
||    ||   ||0.001||style="color:green"| -9.03987772936423e+02 || 424 || 848 || 10.281
|}

The distinct difference of the CG method between MD++ and LAMMPS lies in the fact that MD++ runs the CG method in scaled coordinate ('''x'''<sup>*</sup>) while LAMMPS does in real coordinate ('''x'''). Because of this scaling effect, ''d''<sub>max</sub> must be devided by the box length ''L'' for MD++, or ''d''<sup>*</sup><sub>max</sub> = ''d''<sub>max</sub>/''L''. By the same reason, the force tolerance in MD++ needs to be multiplied by the box length, or ftol<sup>*</sup> = ''L''&middot;ftol.

Now, I compare '''zxcgr''' and '''PR+''' in MD++. Since it is only the force tolerance that the "zxcgr" code considers to check the convergence, the energy tolerance will be set to be 0 in the comparison below. In table 3, the numbers are same upto 1e-8. Note that the energy obtained by the PR+ method becomes lower than the number (in table 2) from the same method with nonzero etol (=1e-6). This further relaxation can be understood, if we consider that in the case of nonzero energy tolerance a small step makes the energy change small enough to satisfy the energy tolerance and the system is prone to stop relaxing further even though it is able to do because of its nonzero gradient (or force).

{| class="wikitable" style="text-align:center"
|+ Table 3. Relaxation of 3[100]x3[010]x3[001] Si structure with a defect. Potential file : MEAM Siz w/ r_c = 6 (A). Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!   !! ''d''<sub>max</sub> (&Aring;) !! ''E''<sub>rlx</sub> (eV) !! Iter. No. !! No. of Calls !! CPU time
|-
! MD++ (zxcgr)
|     || |-9.04000245866814e+02 || 31  || 45  || 0.31602
|- 
! rowspan=3 valign="center" | MD++ (PR+)
| 0.1 ||-9.04000245866095e+02 || 35 || 124 || 0.70004
|-
||0.01||  -9.04000245869413e+02 || 63 || 117 || 0.70004
|-
||0.001|| -9.04000245869795e+02 || 460 || 483 || 2.6002
|}

At this point, we may raise a question. '''Does the scaling effect ever affect the relaxation of a structure in a non-cubic cell?''' Let's say that we have a 6[100] x 3[010] x 3[001] Si diamond cubic structure stretched by 20% in ''x'' direction with two missing atoms (of atom ID 0 and 99). The relaxation results are listed in table 4 below.

{| class="wikitable" style="text-align:center"
|+ Table 4. Relaxation of 6[100]x3[010]x3[001] Si structure with a defect. Potential file : MEAM Siz w/ r_c = 6 (A). Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). For the PR+ method, ''d''<sub>max</sub> = 0.1 (&Aring;). No. of atoms = 430
|-
!   !! ''E''<sub>perf</sub> (eV) !! ''E''<sub>defect</sub> (eV) !! ''E''<sub>rlx</sub> (eV) !! Iter. No. !! No. of Calls 
|-
! MD++ (zxcgr)
| -2.00016000044795e+03  || -1.86579478735442e+03 || -1.86764046397026e+03 || 129  || 175
|- 
! MD++ (PR+)
|  || || -1.86764046392565e+03 || 168 || 686 
|-
! LAMMPS
|  -2.00016000044795e+03 || -1.86579478735442e+03 || -1.86764046393067e+03 || 74 || 387
|}

The relaxed energy values are same upto 1e-6. 

<!-- in the [[#Equilibrium Vacancy Concentration|Appendix]]. -->
<HR>

== Issues with an object function of multiple local minima ==

'''If a structure has multiple local minima that are closely located, can we get the same relaxed structure with different CG methods?''' A Si slab can easily have multiple metastable structures because of the surface reconstruction and for that reason it is a good test system for our purpose. The simulation systems are 5[100]x5[010]x2[001] Si slabs: one with a surface atom of ID 142 intentionally removed and the other with not. (See figure 3.) Both systems are relaxed with zxcgr and PR+ CG methods, and the results are summarized in table 5.

The Si slab with no atom removed was able to reach the same metastable structure regardless of CG methods. However, in the other slab the removed atom triggered the surface atoms' bonding or de-bonding differently depending on the CG method. Apparently, we observe two different surface reconstruction patterns obtained from the same initial structure. 

<!--'''We may ask which structure is physically obtainable.''' The answer would be yes and no for both structures. Fow now, I think it is safe to say that none of them are trustable.-->

{| class="wikitable" style="text-align:center"
|+ '''Table 5''' Relaxation of a 5[100]x5[010]x2[001] Si slab. Potential file : MEAM Siz w/ r_c = 6 (A). Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!  rowspan=2 valign="center" | CG !! colspan="3" | System 1 (No atom removed) !! colspan="3" | System 2 (1 atom removed)
|-
!       Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) !! Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) 
|-
! zxcgr
|  13 ||  29 || -1.69159677459827e+03 || 189 || 328 || -1.72478726205782e+03
|- 
! PR+
|  44 || 183 || -1.69159677459797e+03  || 571 || 1255 || -1.72304961624534e+03
|}

{| border="0"
|+'''Figure 3'''
|-
!  !! rowspan="2" valign="center" align="center"| Initial structure of System 1 !! rowspan="2" valign="center" align="center"| Initial structure of System 2 !! colspan="2" align="center" | Relaxed structure of System 2
|-
!  !! align="center"| zxcgr !! align="center"|PR+ 
|-
| top view 
| align="center" width="150px" | [[Image:Thinslab-top.jpg|180px]]
| align="center" width="150px" | [[Image:Thinslab_rm1-top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-top.jpg|180px]]
|-
| side view
| align="center" width="150px" | [[Image:Thinslab-side.jpg|180px]]
| align="center" width="150px" | [[Image:Thinslab_rm1-side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-side.jpg|180px]]
|}

How do we interprete this behavior? If a given object function has multiple local minima that are very close, a slightly different search path can lead you different energy valleys. This idea can be confirmed by relaxing structures with another Si potential, Stillinger-Weber (SW), which is less bumpy than MEAM. The relaxation results are shown below in table 6 and figure 4. According to the results, both '''zxcgr''' and '''PR+''' methods give the same relaxed structure of the system 2 in terms of energy and the surface pattern. Note that the relaxed structures do not show any appreciable surface reconstruction. 

{| class="wikitable" style="text-align:center"
|+ '''Table 6''' Relaxation of a 5[100]x5[010]x2[001] Si slab with '''SW''' Potential. CG Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!  rowspan=2 valign="center" | CG !! colspan="3" | System 1 (No atom removed) !! colspan="3" | System 2 (1 atom removed)
|-
!       Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) !! Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) 
|-
! zxcgr
|  6 ||  10 || -1.62049999855975e+03 || 7 || 11 || -1.61586999856701e+03
|- 
! PR+
|  7 || 37 || -1.62049999855992e+03  || 7 || 37 || -1.61586999856631e+03
|}

{| border="0"
|+'''Figure 4''' Images of relaxed structures with '''SW''' potential
|-
!  !! colspan="2" align="center" | Relaxed structure of System 2
|-
!  !! align="center"| zxcgr !! align="center"|PR+ 
|-
| top view 
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-sw-top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-sw-top.jpg|180px]]
|-
| side view
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-sw-side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-sw-side.jpg|180px]]
|}

This comparison does not mean that SW potential performs better than MEAM or the other way around. We just want to point out that there is a feature of a potential function that may unintendedly affect the relaxation results, and we have to be very careful when taking them.

By the way, it is reported that the undesirable metastable states of MEAM potential can be removed by fine-tuning of parameters related with the multi-body screening function.<ref>Seunghwa Ryu, Christopher R. Weinberger, Michael I. Baskes, and Wei Cai, "Improved Modified Embedded-Atom Method Potentials for Gold and Silicon", submitted to ''Modelling and Simulation in Materials Science and Engineering'' (2009)</ref> Let's try C<sub>min</sub> = 0.8 and see what comes out of this change. The relaxation results are summarized below.

{| border="0"
|-
|
{| class="wikitable" style="text-align:center"
|+ '''Table 7''' The energy of the relaxed system 2 with '''MEAM Siz''' Potential. '''C'''<sub>'''min'''</sub>=0.8, CG Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!  rowspan=2 valign="center" | CG !! colspan="3" | System 2 (1 atom removed)
|-
!     Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) 
|-
! zxcgr
|  180 ||  304 || -1.72217359502167e+03
|- 
! PR+
|  577 || 1323 || -1.72217359501501e+03
|}
||
{| border="0"
|+'''Figure 5''' Images of the relaxed structure of the system 2 with '''MEAM Siz''' Potential. '''C'''<sub>'''min'''</sub>=0.8
|-
!  !! align="center"| zxcgr !! align="center"|PR+ 
|-
| top view 
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-cmin0.8-top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-cmin0.8-top.jpg|180px]]
|-
| side view
| align="center" width="150px" | [[Image:Rlx_rm1_zxcgr-cmin0.8-side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_rm1_prplus-cmin0.8-side.jpg|180px]]
|}
|}

Typically, people use the simulatetd annealing (SA) method to avoid the undesirable meta-stable states. In the simulated annealing, the system is heated up and slowly cools down to overcome energy barriers. Figure 6 shows exponential decreas of temperature in simulated annealing process. When temperature becomes lower than 100K, the CG relaxation starts. According to table 7 and figure 7, both the zxcgr and the PR+ relaxations find the same minimum structure, which may or may not be the global minimum structure.

{| border="0"
|-
|
{| border="0"
|+'''Figure 6''' 
|-
| align="center"  | [[Image:Sa_temperature.jpg|400px]]
|-
|}
{| class="wikitable" style="text-align:center"
|+ '''Table 7''' The energy of the relaxed system 2 after simulated annealing. '''MEAM Siz''' Potential. '''C'''<sub>'''min'''</sub>=<font color = 'red'>2.0</font>, CG Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!  rowspan=2 valign="center" | CG !! colspan="3" | System 2 (1 atom removed)
|-
!     Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) 
|-
! zxcgr
|  135 ||  169 || -1.73045698795218e+03
|- 
! PR+
|  173 || 818 ||  -1.73045698794833e+03
|}
||
{| border="0"
|+'''Figure 7''' Images of the relaxed structure of the system 2 after simulated annealing. '''MEAM Siz''' Potential with '''C'''<sub>'''min'''</sub>=<font color = 'red'>2.0</font>
|-
!  !! align="center"| zxcgr !! align="center"|PR+ 
|-
| top view 
| align="center" width="150px" | [[Image:Rlx_zxcgr_SA_top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_prplus_SA_top.jpg|180px]]
|-
| side view
| align="center" width="150px" | [[Image:Rlx_zxcgr_SA_side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_prplus_SA_side.jpg|180px]]
|}
|}

{| border="0"
|-
|
{| class="wikitable" style="text-align:center"
|+ '''Table 8''' The energy of the relaxed system 2 after simulated annealing. '''MEAM Siz''' Potential. '''C'''<sub>'''min'''</sub>=<font color = 'red'>0.8</font>, CG Conditions: etol = 0 (no unit), ftol = 1e-4 (eV/&Aring;). 
|-
!  rowspan=2 valign="center" | CG !! colspan="3" | System 2 (1 atom removed)
|-
!     Iter. No. !! No. of Calls !! ''E''<sub>rlx</sub> (eV) 
|-
! zxcgr
|  123 ||  168 || -1.72878420193421e+03
|- 
! PR+
|  180 || 854 ||  -1.72878420193527e+03
|}
||
{| border="0"
|+'''Figure 8''' Images of the relaxed structure of the system 2 after simulated annealing. '''MEAM Siz''' Potential with '''C'''<sub>'''min'''</sub>=<font color = 'red'>0.8</font>
|-
!  !! align="center"| zxcgr !! align="center"|PR+ 
|-
| top view 
| align="center" width="150px" | [[Image:Rlx_zxcgr_SA_cmin0p8_top.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_prplus_SA_cmin0p8_top.jpg|180px]]
|-
| side view
| align="center" width="150px" | [[Image:Rlx_zxcgr_SA_cmin0p8_side.jpg|180px]]
| align="center" width="150px" | [[Image:Rlx_prplus_SA_cmin0p8_side.jpg|180px]]
|}
|}

<!--<gallery caption="Fig.3" widths="200px" heights="200px" perrow="4">
Image:Linesearch1.jpg|(a) 
Image:Linesearch2.jpg|(b) 
Image:Linesearch3.jpg|(c) 
Image:Linesearch3.jpg|(d) 
Image:Linesearch1.jpg|(a) 
Image:Linesearch2.jpg|(b) 
Image:Linesearch3.jpg|(c) 
Image:Linesearch3.jpg|(d)
</gallery>-->
<HR>

== Notes ==

<references/>