Implementation of Nose-Hoover thermostat using Predictor Corrector in Molecular Dynamics

In molecular dynamics (MD) simulations temperature can be controlled via numerous methods. Out of those methods the most popular method is Nose-Hoover (NH) thermostat. The steps that are to be followed to implement NH equations in MD using Predictor Corrector algorithm is a bit challenging. So, I decided to lay an outline for the same. The steps are:

1. Prior to entering the main MD loop, evaluate forces [through a function called force( ) ]once and initialize the damping parameters (call it ${\eta}$) at 0.

2. In the main MD loop (let us say you want N simulations) call a function [let us say integrate( ) ] that evaluates the updated position and velocities

3. Within integrate( ):
    3.1. Predictor Steps: 
          3.1.a) Evaluate force( ) and let us say is stored in a vector called ${F}$

          3.1.b) Predict the updated positions through:

 ${x\left(t+h\right) = x\left(t\right) + v\left(t\right)h + \dfrac{1}{2}a\left(t\right)h^2 + \dfrac{1}{6}b\left(t\right)h^3}$

          3.1.c) Predict the updated velocities through:

 ${v\left(t+h\right) = v\left(t\right) + a\left(t\right)h + \dfrac{1}{2}b\left(t\right)h^2}$

          3.1.d) Predict the updated forces through:

 ${F\left(t+h\right) = F\left(t\right) + b\left(t\right)h}$

          3.1.e) Evaluate whether the predicted positions satisfy the periodic boundary conditions

   3.2.  Corrector Steps:
          3.1.a) Evaluate force( ) let us say the force output is stored in vector ${F_{new}\left(t+h\right)}$

          3.1.b) Update the new force to reflect the effect from bath:

${F_{new}\left(t+h\right)=F_{new}\left(t+h\right)-\eta v\left( t+h \right)}$ 

          3.1.c) Calculate the difference between the ${F_{new}}$ and ${F}$:

  ${factor = F_{new}\left(t+h\right) - F\left(t+h\right)}$

          3.1.d) Correct the updated positions through:

${x_c\left(t+h\right) = x\left(t+h\right) + p \times factor h^2}$

          3.1.e) Correct the updated velocities through:

${v_c\left(t+h\right) = v\left(t+h\right) + q \times factor h}$     

          3.1.f) Correct the updated forces through:

${F_c\left(t+h\right) = F\left(t+h\right) + r \times factor}$

          3.1.g) Correct the updated b through:

 ${b_c\left(t+h\right) = b\left(t+h\right) + \dfrac{s \times factor}{h}}$ 

          3.1.h) Check for periodic boundary conditions with the updated corrected positions

   3.3.  Update eta with the expression using NH equations of motion

4. Repeat step 3 till N simulations are complete.

NOTE: In Gear's algorithm: ${p = 1/12, q = 5/12, r = 1}$ and ${s = 1}$
PS: Please feel free to contact if you need any help or if you find any mistakes