Dear LAMMPS users,

I want to follow-up on a question that was discussed on the forum: http://lammps.sandia.gov/threads/msg34227.html

I would like to control the pressure on a non-periodic dimension using a moving harmonic wall, as suggested by Steve in the above discussion.

The “barostat” works for some time and then becomes chaotic and crash.

I made a toy system made of water molecules periodic in y, z but not in x where I put a harmonic wall at the position defined by the following equation (I consider only the xlo side here):

x = x0 + v*t

The velocity of the wall is a function of time and adjusted every cycle from the velocity at the previous cycle (with a loop):

v(i) = v(i-1) + gamma*dv

Gamma is given by:

gamma = (Pext - P)/Pext

with Pext the external pressure and P the actual wall pressure computed from the force on the wall via: P = 68568.415*abs(fwall)/area (in real units) and fwall is time-averaged over 10 steps corresponding to one cycle.

Everything work as expected at the beginning and the pressure equilibrates well however after ~500 cycles, the pressure starts to diverge and rapidly becomes chaotic and the simulation crash.

I would appreciate any advice to solve this.

I know there are other simulation options but I just would like to understand this behavior.

I attach a plot of the actual pressure P (Atm) as a function of cycles (10 MD steps) with Pext = 1000 Atm as well as the corresponding input, force field and data files.

I use REAX/C & LAMMPS-17Nov16.

Thank you and regards,

Nicolas

water-Pext.in (3 KB)

water.data (40 KB)

ffield (28 KB)

Dear LAMMPS users,

I want to follow-up on a question that was discussed on the forum:

http://lammps.sandia.gov/threads/msg34227.html

I would like to control the pressure on a non-periodic dimension using a

moving harmonic wall, as suggested by Steve in the above discussion.

The “barostat” works for some time and then becomes chaotic and crash.

I made a toy system made of water molecules periodic in y, z but not in x

where I put a harmonic wall at the position defined by the following

equation (I consider only the xlo side here):

x = x0 + v*t

The velocity of the wall is a function of time and adjusted every cycle

from the velocity at the previous cycle (with a loop):

v(i) = v(i-1) + gamma*dv

Gamma is given by:

gamma = (Pext - P)/Pext

with Pext the external pressure and P the actual wall pressure computed

from the force on the wall via: P = 68568.415*abs(fwall)/area (in real

units) and fwall is time-averaged over 10 steps corresponding to one

cycle.

Everything work as expected at the beginning and the pressure equilibrates

well however after ~500 cycles, the pressure starts to diverge and rapidly

becomes chaotic and the simulation crash.

I would appreciate any advice to solve this.

when you adjust the wall to the full extent in every step, you open

yourself to feedback and your plot looks like that happened to you. the

system you construct is not very different from an oscillator.

you want a damped or better and overdamped oscillator. check out the

berendsen thermostat algorithm. in short, you don't want to apply gamma,

but gamma times the damping factor, which should be of the order of 1000 to

100000 time steps.

axel.

Dear Axel,

I went over Berendsen’s paper and slightly modified my input file accordingly:

mu = 1-(timestep/tau)*(P-Pext)

v_i = v_(i-1)*mu

With tau = 10,000, I obtain reasonable equilibration of the pressure.

I attach P in Atm = f(cycle) on both walls with Pext = 1,000 Atm.

Thank you.

Nicolas