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