M04A Conjugate Gradient Method in MD++

Manual 04A for MD++
Conjugate Gradient Method in MD++

Keonwook Kang and Wei Cai

Mar 19 , 2009

---

Conjugate Gradient Method

The conjugate gradient (CG) method was initially proposed as an iterative method to solve a large linear system, such as $Ax=b$. Later, Fletcher and Reeves adopted it to solve nonlinear optimization problems. In MD++, the CG method is used to find a (local) minimum structure by minimizing a given (usu. nonlinear) potential energy function.

Algorithmically, the structure of the CG method is composed of two nested loops. The outer loop tries different search directions, which are A-orthogoanl to each other for a linear system A. The inner loop aims to find the optimal step length which minimizes the object function along a given search direction. When certain conditions are met, the program exits out of the loops and stops.

$x_{k+1}=x_{k}+\alpha _{k}s_{k}$.

where $\alpha _{k}$ and $s_{k}$ are the step length and the search direction at k-th iteration, respectively, and k is the index of the outer loop.

Implementation of the CG method in MD++

For more information, see supplementary document: Vacancy formation energy by MD++ .

We can run this script file by typing

$ bin/fs gpp scripts/movacancy.script

First, MD++ creates a 5×5×5 perfect cubic crystal of Mo with edges along <100> directions. Then it fixes atom 0 by

   input = [ 1         # number of atoms to be fixed
             0 ]       # index of an atom to be fixed
   fixatoms_by_ID

and then removs this atom by the command removefixedatoms. The first number in the input array specifies the number of atoms to be removed. For example, if you would like to remove two atoms, say 3 and 8, you can do so by

input = [ 2            #  number of atoms to be fixed
          3 8 ]        #  index of atoms to be fixed
fixatoms_by_ID         #  fix a set of atoms by their index
removefixedatoms       # remove fixed atoms

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 214. The potential fuction is MEAM Siz with the cut-off radius r_c=6(Å) 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/Å), respectively. In the 1st column of table 1, 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_max (0.1, 0.01, and 0.001), 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.

Table 1. 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/Å). Compiler = intel compiler

	Eperf (eV)	Edefect (eV)	dmax (Å)	Erlx (eV)	Iter. No.	No. of Calls	CPU time
MD++ (PR+)	-1.00008000022397e+03	-9.00066105589930e+02	0.1	-9.04000029718829e+02	15	46	0.30002
			0.01	-9.04000167032814e+02	46	51	0.33202
			0.001	-9.03987772936424e+02	424	425	2.2201
LAMMPS (Serial Run)	-1.00008000022397e+03	-9.00066105589930e+02	0.1	-9.04000029718830e+02	15	60	0.75296
			0.01	-9.04000167032815e+02	46	96	1.2149
			0.001	-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 the scale coordinate (x^*) while LAMMPS does in the real coordinate (x). Because of this scaling effect, d_max must be devided by the box length L for MD++, or d^*_max = d_max/L. By the same reason, the force tolerance in MD++ needs to be multiplied by the box length, or ftol^* = L·ftol.

MD++ is equipped with another implementation of the CG method, called as "zxcgr". This "zxcgr" code was written by Dongyi Liao at MIT in 1999 based on IMSL's U2CGG.F, and used the Hestenes-Stiefel method and the quadratic spline interpolation. Now I will compare two implemetations; MD++(zxcgr) and MD++(PR+). 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. The numbers are same upto 1e-8. Nothe that the energy obtained by the PR+ method becomes lower when compard with the number from the same method with nonzero etol (=1e-6).

Table 2. 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/Å).

	dmax (Å)	Erlx (eV)	Iter. No.	No. of Calls	CPU time
MD++ (zxcgr)		-9.04000245866814e+02	31	45	0.31602
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 this scaling effect ever affect the relaxation of a structure in a non-cubic cell? Let's assume that we have a 6[100] x 3[010] x 3[001] Si diamond cubic structure stretched by 10% in x direction with two missing atoms (of atom ID 0 and 99).

in the Appendix.