[lammps-users] Fix aveforce piston barostat

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

First of all, my LAMMPS version is 3-Mar20

I would like to simulate a system of water confined between
two walls. My intention is to allow the walls to move in the “z” direction
whereas “x, y” dimensions remain during the simulations (Np_zAT ensemble).

The idea which I thought should work was as follows:

  1. define the groups

group water type 1 2
group wall1 type 3
group wall2 type 4

  1. integrate all of the atoms and thermostat only the water molecules
    fix verlet all nve
    fix trescale water temp/csvr 262 262 1000 123321

  2. Set the force of wall molecules to 0.0 in x and y direction

whereas not change in the z direction.

fix rigiddown wall1 setforce 0.0 0.0 NULL
fix rigidup wall2 setforce 0.0 0.0 NULL

  1. Apply the additional force to maintain the desired pressure.

fix pressdown wall1 aveforce NULL NULL -0.034727274
fix pressup wall2 aveforce NULL NULL 0.034727274

For the system which has the area of 1131 A^2

I assumed the pressure of 200 atm should be approx.

-0.034727274 J/(A-mol) for walls built of 96 molecules each.

However, the issue is that depending on the initial configurations
I obtain completely different values for the ZZ component of the pressure
tensor in this system, which are approx. 0 +/- 50 atm.

Does anyone have any suggestion on how to maintain
the desired pressure with this kind of “piston barostat”?
Or am I missing something?
Alternatively, is there another suitable way to do these kinds of simulations?

Best wishes,

2

I am sorry for the typo, but the force value here:

"or the system which has the area of 1131 A^2

I assumed the pressure of 200 atm should be approx.

-0.034727274 J/(A-mol) for walls built of 96 molecules each."

has the correct value but the units should be Kcal/(A-mol) as specified
in the documentation.

3

A couple of thoughts.

(a) What happens if you don’t use fix aveforce on the Z walls, and just freeze them, like the XY walls?
Does the pressure also vary depending on initial condition? I.e. the use of fix aveforce is not the issue.

(b) What happens if you run with fix aveforce for a while, then remove the fix, freeze the walls, and run some more.
Does the equilibrated pressure depend on the initial condition? Does the final position of the frozen Z wall depend
on the initial condition? I’m thinking the different pressures you are getting are simply a function of the wall
position, not fix aveforce?

(c) Depending on how you define pressure, there is probably a contribution to Pzz from fix aveforce itself.
With fix addforce there is a fix_modify virial yes option to add that extra force to the virial. See the bottom of
its doc page. That’s a relatively new option in LAMMPS. Not sure why we don’t provide that option with fix aveforce as well.

I’m CCing Aidan to comment on (c) whether that makes sense to do for fix aveforce.

Steve

4

Dear Steve,

After multiple further tests, I realized that the issue is different, i.e. very short runs and small mass of the walls
on which the force was exerted. It has led to large fluctuations. With large piston atoms mass, the fluctuations

are of the order of few atmospheres, so similar as for the regular global barostat.

For the pressure calculations, I did not add virial components from the additional force.

If one would want to do so, it can be easily done with one additional fix, i.e.

  1. setforce 0.0 0.0 null

  2. addforce 0.0 0.0 +/-{force}

  3. aveforce 0.0 0.0 null

and then it does not matter whether aveforce has this option or not.

Another major issue is how the LAMMPS is calculating the pressure, i.e. if we specify the
force that has to be exerted on the upper and lower wall to maintain the desired pressure,
then what volume should be taken into account, (i) entire system’s or (ii) only the fluid

confined between those two slabs? If it is the latter, then it is not trivial to calculate

that volume, as fluid can be contracted/expanded so this is a dynamical variable.

However, I assumed that this contraction/expansion are of a negligible order
so I assumed the vol_fluid to be initial volume and calculated the pressure as follows:

compute peratom all stress/atom temp
compute ptens all reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
variable pz equal -(c_ptens[3])/vol_fluid

and I output that with thermo command.

This approach leads to definitively closer results compared to the previous one, however,
it is not perfect and the lower the pressure I want to simulate the larger errors I receive.

The piston configuration is as follows:
N_p=192 of each upper and lower wall and those walls are not able to interact with one another,
but only with oxygen atoms of water molecules.

To determine the desired pressure, I assumed the following routine to estimate the forces:

f_z=(press*A_p)/N_p
where A_p is area of the piston, N_p is number of piston molecules and f_z is the force in the z-direction.
Therefore, for the pressure of 400 atmospheres, we need f_z=+/-0.034636126133563525 kcal/(mol-A).
In the above calculations with stress/atom pressure definition, the average however is approx 385 atmospheres
with the standard deviation of ~500 atmospheres…

If anyone could suggest me what am I doing incorrectly I would be very grateful.
This idea is for sure possible to be performed with LAMMPS (there are few papers
which have this approach being used)

Best wishes,

5

I think for this system you want the volume to be only for the fluid between the walls. The inner surface of the wall is acting
like the simulation box boundary in a periodic system.

Likewise for the kinetic component of the pressure you want only the fluid, not the wall atoms.
I.e. they shouldn’t contribute to N or the temperature calculation. You can define your own
temperature compute to work only on the fluid atoms and assign that to your own pressure compute
and output both with thermo if that helps.

Actually you could define multiple pressure computes, one with just the kinetic term, the other with just
the virial. The latter would allow you to scale the volume to subtract the thickness of the z-walls, e.g. by
multiplying the virial term by a scale factor in a formula you define.

So it’s a bit of work, but you should be able to output the value you are looking for. Whether it will be “right”
or not is another Q … :slight_smile:

Steve