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/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/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.

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.