Page 1 of 1
MD simulations with constraints
Posted: Mon Feb 06, 2023 4:03 am
by ruoyuwang
Hi,
I noticed on the VASP page that you use the leap-frog and SHAKE algorithms when performing the constrained molecular dynamics. However, when considering the thermostat, are you still using the leap-frog algorithms? Or do we need to use the velocity verlet algorithm? Thank you!
Best,
Ruoyu
Re: MD simulations with constraints
Posted: Mon Feb 06, 2023 2:58 pm
by andreas.singraber
Hello Ruoyu,
welcome to the VASP forum! We are not entirely sure where your question is heading because the leapfrog and velocity-Verlet algorithm are
closely related. In VASP there is only the Verlet integrator implemented, hence there is also no
INCAR tag to change it. In order to use constrained molecular dynamics you only need to define your constraints in the
ICONST file and use any of the
MDALGO=1,2,3 methods to thermostat and/or barostat your system. Then, the SHAKE algorithm is automatically applied in your simulation.
Best,
Jonathan Lahnsteiner
Andreas Singraber
Re: MD simulations with constraints
Posted: Tue Feb 07, 2023 7:11 am
by ruoyuwang
Hi,
Thank you for your reply! I found that you use leap-frog and SHAKE algorithms with the Nose-Hoover thermostat. However, you use velocity verlet and RATTLE algorithms with Langevin dynamics. How do you consider the choice between SHAKE and RATTLE algorithms? Thank you!
Best,
Ruoyu
Re: MD simulations with constraints
Posted: Wed Feb 08, 2023 1:51 pm
by toms_bucko1
Hello Ruoyu,
RATTLE is just a variant of SHAKE adapted for the use with the velocity Verlet algorithm. Hence, the thermostat algorithms implemented with the leapfrog integrator (MDALGO=1|2) use SHAKE while those that come with velocity Verlet use RATTLE. Note that both integrators are equivalent provided a fixed time step is used (this can be easily seen by suitably rearranging time propagation sequences, see e.g.
https://en.wikipedia.org/wiki/Leapfrog_integration), which is always the case in MD with VASP. You can convince yourself about the equivalence of both integrators by running short MD simulations (<50 steps) in NVE ensemble using POSCAR in which velocities are explicitly defined (e.g., POSCAR from a previous MD run). For that, you can use MDLAGO=1;ANDERSEN_PROB=0 (leapfrog) and MDALGO=3;LANGEVIN_GAMMA=0 (velocity Verlet). By checking potential energies (grep F OSZICAR|awk '{print $9}') you'll see that both runs will generate identical sequences and the same holds true whether you use constraints or not.
All the best,
tb