Hello,
I am trying to write a fix to implement a langevin piston barostat in lammps. I was wondering if anyone could help me understand what I need to do so that it interacts correctly with the shake command. I am looking at the unconstrained update
xshake[i][0] = x[i][0] + dtvv[i][0] + dtfmsqf[i][0]
Now, I am trying to understand if this equation needs to be modified to account for the change in the equations of motion due a barostat. Specifically, do I need to add a term to f[i][0] here to account for the acceleration due to volume expansion? In the langevin piston, the equations for the position and velocity are
dr/dt = v + r (1/3)d(volume)/dt
dv/dt = f/m -1/3v(d(volume)/dt)/volume
(Scott E. feller, constant pressure molecular dynamics simulation: the langevin piston method, 1995)
So i am wondering if I need to add -1/3v(d(volume)/dt)/volume and make the shake update
xshake[i][0] = x[i][0] + dtvv[i][0] + dtfmsqf[i][0] - dtvdtv0.51/3v*(d(volume)/dt)/volume
I see the comment in the code that the unconstrained update seems to work well enough for NPT, but I didn’t want to interpret this as no modification to the force is taking place without verifying with the lammps community.
Sorry for all the equations.
Conor
Hello,
I am trying to write a fix to implement a langevin piston barostat in
lammps. I was wondering if anyone could help me understand what I need to do
so that it interacts correctly with the shake command. I am looking at the
unconstrained update
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0]
Now, I am trying to understand if this equation needs to be modified to
account for the change in the equations of motion due a barostat.
Specifically, do I need to add a term to f[i][0] here to account for the
acceleration due to volume expansion? In the langevin piston, the equations
for the position and velocity are
dr/dt = v + r *(1/3)*d(volume)/dt
dv/dt = f/m -1/3*v*(d(volume)/dt)/volume
(Scott E. feller, constant pressure molecular dynamics simulation: the
langevin piston method, 1995)
So i am wondering if I need to add -1/3*v*(d(volume)/dt)/volume and make the
shake update
xshake[i][0] = x[i][0] + dtv*v[i][0] + dtfmsq*f[i][0] -
dtv*dtv*0.5*1/3*v*(d(volume)/dt)/volume
I see the comment in the code that the unconstrained update seems to work
well enough for NPT, but I didn't want to interpret this as no modification
to the force is taking place without verifying with the lammps community.
if you worry about such details, you probably shouldn't be using fix
shake, but rather fix rattle and take it from there. i am not an
expert in this, but my understanding from the detailed discussion on
the mailing list with the author of fix rattle suggest that shake does
at best an approximation of RATTLE, but that for proper time
integration explicit implementation of SHAKE and RATTLE are needed.
while the former is sufficient for a regular velocity time
integration, it is not for velocity verlet as implemented in LAMMPS.
RATTLE contains the additional steps required for velocity updates.
axel.