CSCI 653 Final Project

1. Introduction

    Conventional Molecular Dynamics (MD) based on empirically fitted potential has been established and widely used in all the fields of researches for a long time. It has been able to successfully describe many materials' mechanical properties in certain range. However, since it's fitted to the bulk solid properties, it often fails when encountered with changed situations. For example, systems with defects or propagating cracks will have partial charges of atoms vary according to local environment.

    To make the conventional MD simulation more transferable, and to make the description more realistic, such that it can describe certain chemical reaction properties, bond order and charge variables are recently included many researches. In our current projects, we have been using the Reactive Force Field (ReaxFF) to describe the explosion of RDX and also the water enhanced indentation on amorphous silica glass (aSiO2).

    However, in ReaxFF (or most of the potential including charge variable descriptions), the computation costs in determing the equalized electronegativity often becomes the bottle neck of total computation time for the MD simulation. The purpose of the current project is then to improve the efficiency of QEQ calculation in our current code.

2. Algorithms

    To obtain the dynamic charge distribution in ReaxFF, currently it's done by the Unconstrained Two-Vector Algorithm for Constrained Energy Minization[1]. This algorithm uses two kinds of fictitious charge under charge-neutrality constraint. Electrostatic energy is minized by conjugate gradient (CG) method. The number of iteration needed for conjugate gradient to converge greatly depends on the system's initial conditions. In the project, the CG procedure is tested on a Hydrogen covered pyramid shape diamond indenter (a smaller version is shown in Fig.1). To maintain the charge neutrality of the system, subroutine QEq is implemented at each MD step. Each call of the QEq subroutine needs about 70 iterations of CG steps (See Fig. 3). The time cost for the QEq subroutine is often even more than the calculation of 4-body interactions.

     In order to speed up the total calculation speed, one way is to reduce the iteration number needed in each MD step for conjugate gradient energy minimization. The method of Preconditioned Conjugate Gradient (PCG) has been proposed for different situations for conjugate gradient minimization. Here we would also like to use PCG for improving the efficiency of QEQ subroutine.

    There are many ways to choose the preconditioning matrix for PCG. Here we would like to use something similar to what we have done before as dividing the total matrix into long range interaction part and short range interaction part, and use the short range part as the preconditioning matrix. Results are discussed in the following sections.

    An interesting thing noticed in implementing the PCG algorithm is that, by updating the conjugate direction iteratively along, instead of calling an external function to calculate the gradient from the derivate of the energy, can readily improve the QEq subroutine's efficiency about 40%. This result is shown below in section 3.

    For details of the CG and PCG algorithms, please refer to Dr. Nakano's paper [2].

3. Results

    1. The following figure shows the comparison of QEq time cost between the schemes a, b, and c. From the graph it can be seen that scheme b and c both have significant improvement on the original scheme a. But difference between scheme b and c is not big. And the averaged time cost for the three schemes are shown below in Table. 1. 

    2. This figure shows the comparison of number of iterations in each CG call between the schemes a, b, and c. This plot shows that the number of iterations between the three schemes are almost the same, except for that scheme (b) and (c) have their local minimums at different time steps. scheme (a) and (b) have the same result because they are basically the same procedure, except for used different methods for obtaining the conjugate gradient direction. scheme (c) actually serves as a test scheme, because the preconditioning matrix used here is simply the diagonal matrix of the full Hessian matrix. So, the improvement on the iteration number is not obvious.

    3. his figure shows the comparison of total MD time for each MD step between the schemes a, b, and c. The plot again shows the similar result as given in figure 2. The averaged time cost for the three schemes are shown below in Table. 2.

    4. It should also be mentioned here that scheme (d) was not forgot in the project. However, it has been found out that the convergence of the preconditioning procedure never comes out good by simply choosing a shorter range of interaction as the preconditioning matrix. The reason is currently considered as that the preconditioning matrix may not be a positive definite matrix. Because from the scheme (c), it has been shown that for a simple diagonal matrix the preconditioning works, even though no significant improvement of iteration number realized.

4. Conclusion

References:

[1] Write-up for Unconstrained Two-Vector Algorithm for Constrained Energy Minimization, Ken-ichi Nomura, not published.

[2] Parallel multilevel preconditioned conjugate-gradient approach to variable-charge molecule dynamics, Aiichiro Nakano, Comp. Phys. Comm, 104 (1997) 59-69.