Predictor - corrector method

Amongst the different versions of integrators, the predictor-corrector ( PC ) [54-56], method was chosen for our simulations. The goal is to solve the second order ordinary differential equations:

$$\ddot x_i=f_i(x_i,\dot x_i,t)$$ (8.1)

where

$$\dot x=\frac {dx}{dt}$$ (8.2)

and

$$f(t) \equiv f(x(t),\dot x(t),t)$$ (8.3)

This is Newton's equations of motion for the $x$ component of the $i$th atom, and the same equation could be written for $y$ and $z$. We omit the $i$ index now for simplicity. The PC method is composed of three steps: prediction, evaluation, and correction. In particular, from the initial position $x(t)$ and velocity $\dot x(t)$ at time t, the steps are as follows:

Prediction
1.) Predict the position $x(t+h)$, where $h$ is the time step of integration:

$$x(t+h)=x(t)+h\dot x(t)+h^2\sum_{i=1}^{k-1} \alpha _i f(t+h(1-i))$$ (8.4)

In our simulation we use k=4, and therefore we can write explicitly the formula, using the $\{ \alpha_i \}$ from [54] :

$$x(t+h)=x(t)+h \dot x(t)+h^2 \left \{ \frac{19}{24} f(t) - \frac{10}{24} f(t-h)+\frac{3}{24} f(t-2h) \right \}$$ (8.5)

2.) Predict the velocities:

$$\dot x(t+h)=\frac{x(t+h)-x(t)}{h}+h\sum_{i=1}^{k-1} \alpha' _i f(t+h(1-i))$$ (8.6)

and therefore:

$$\dot x(t+h)=\frac{x(t+h)-x(t)}{h}+h\left \{ \frac{27}{24} f(t) - \frac{22}{24} f(t-h)+\frac{7}{24} f(t-2h) \right \}$$ (8.7)

Evaluation
3.) Evaluate the force using the predicted values:

$$f(t+h)\equiv f(x(t+h))$$ (8.8)

Correction
4.) The final step is to correct the predictions by using some combinations of the predicted and previous values of position and velocity:

$$x(t+h)=x(t)+h\dot x(t)+h^2\sum_{i=1}^{k-1} \beta _i f(t+h(2-i))$$ (8.9)

or

$$x(t+h)=x(t)+h\dot x(t)+h^2 \left \{ \frac{3}{24} f(t+h) + \frac{10}{24} f(t)-\frac{1}{24} f(t-h) \right \}$$ (8.10)

and for velocity:

$$\dot x(t+h)=\frac{x(t+h)-x(t)}{h}+h\sum_{i=1}^{k-1} \beta' _i f(t+h(2-i))$$ (8.11)

and therefore:

$$\dot x(t+h)=\frac{x(t+h)-x(t)}{h}+h\left \{ \frac{7}{24} f(t+h) + \frac{6}{24} f(t)-\frac{1}{24} f(t-h) \right \}$$ (8.12)

The set of parameters $\{\alpha_i,\alpha'_i,\beta_i,\beta'_i \}~~i=1,2,3$ promotes numerical stability of the algorithm. Gear [56] determined their values by applying the predictor-corrector algorithm to linear differential equations and analyzing the stability of the method.

The values of $\{\alpha_i,\alpha'_i,\beta_i,\beta'_i \}$ were chosen specially to make the local truncation error of order of $O(h^4)$, and the global error for the second order differential equations is order of $O(h^3)$. It worth to mention, that interactions are evaluated using the results of the predictor step, but they are not re-evaluated again following the corrector step. It was shown that mean error, induced by the absence of force re-evaluation, is insignificant [55].

The main ingredient of the predictor - corrector method is the corrector step, which accounts for a feedback mechanism. The feedback can damp the instabilities that might be introduced by the predictor step. This method provides both the positions and velocities of the atoms at the same time and it can be used to calculate forces, which depend explicitly on the velocity, this is in turn usually needed in algorithms which control temperature and pressure.