something wrong in NPT

[fanheng]

Sorry to bother. I don’t see anything wrong in the script. But this reminds me a concern when I run a similar system in DPD.

I almost the same thing. One solid slab in the center and liquids around.

My concern is the pressure coupling is based on the total pressure tensors in the entire simulation box including the contributions from the liquid-solid interactions (is this correct?). If this is the case, the pressure in the liquid phase will not equal to the reference pressure we set in the pressure coupling but a lower value especially for smaller systems (I suppose most of the simulations systems can be considered small…). I confirmed it with local pressure measurement based on per-atom stress. The pressure in the bulk liquid phase is lower than the reference pressure (I excluded the interaction between solid as well, if that is included, the pressure will be way off), though not much, but still big enough to influence some sensitive quantaties (e.g. interfacial tension).

And this situation is somehow get more serious if I only apply pressure coupling to z direction. The pressure tensor component Pxx,Pyy will shift to lower values making the pressure in the bulk even lower. (I don’t know if it is true in MD with real atoms). I used iso in pressure coupling instead of in z even though It is not so correct to apply iso.

So is there any way that the pressure coupling can be implemented without using the total pressure tensor? The only way I can think of now is making a caliberation curve.

Cheers,
Heng

2011-09-23

[Andrew_Jewett]

Dear Heng

Thank you for the many interesting questions. Your question made me
realise that I have not tested my system sufficiently rigorously to
know my method works well at high pressure. (See below.)

Sorry to bother. I don't see anything wrong in the script. But this reminds
me a concern when I run a similar system in DPD.

I almost the same thing. One solid slab in the center and liquids around.

My concern is the pressure coupling is based on the total pressure tensors
in the entire simulation box including the contributions from the
liquid-solid interactions (is this correct?).

I think I understand your first question.
The answer to this question depends on whether or not you want the
coordinates of the solid to move as the system expands or contracts at
constant pressure:

1) If you want the solid atoms to expand/contract, then -I-think- you
want to include ALL the forces when you calculate the virial (the
pressure tensor). This includes the contributions from the liquid
atoms pushing on the solid atoms. (In other words: calculate the
pressure tensor in the ordinary way.)

2) However if you do not want the atoms in the solid to move, then you
should not include any forces acting on forces acting on the solid
when you calculate the virial.

We must use "dilate partial" when we apply "fix npt" to the atoms in
the liquid to prevent rescaling of the solid atoms. Here is the
relevant part of my script:

fix 2 mobile npt temp 300.0 300.0 100.0 z 1.0 1.0 1000.0 dilate partial

For this to work 3 things must be true:

a) Forces between atoms in the solid should not be included in the virial.

In my script I use this command to prevent these forces from being calculated:

   neigh_modify exclude group immobile immobile

  ("immobile" is the name of the group of atoms in the solid.)

a) Forces from the liquid pushing on the solid also should not be
included in the virial.

In my script I use this command:

fix 1 immobile rigid single force * off off off torque * off off off

I hope this command has the effect of excluding these forces from the virial
(in addition to keeping the "immobile" group stationary).

Unfortunately I am not sure if this is correct. I am guessing.
(Trung or Tony please correct me if I am wrong. Exactly how does "fix
rigid" modify the virial?)

c) However forces from the solid pushing on the liquid should be
included in the virial.
My "proof" of this is located in equation 9 of the PDF I have attached.

If this is the case, the
pressure in the liquid phase will not equal to the reference pressure we set
in the pressure coupling but a lower value especially for smaller systems (I
suppose most of the simulations systems can be considered small...).

In my simulations, I am satisfied if the liquid has the correct pressure.

(By the way, I don't care what is the pressure inside the solid group
because I use "fix npt ... dilate partial" to prevent the solid atoms
from expanding or contracting. This is option "2)" above. Those atoms
don't move and I don't allow LAMMPS to calculate the forces between
them.)

confirmed it with local pressure measurement based on per-atom stress. The
pressure in the bulk liquid phase is lower than the reference pressure (I
excluded the interaction between solid as well, if that is included, the
pressure will be way off), though not much, but still big enough to
influence some sensitive quantaties (e.g. interfacial tension).
...
So is there any way that the pressure coupling can be implemented without
using the total pressure tensor? The only way I can think of now is making a
caliberation curve.

Perhaps there are two problems that you could be having:

i) LAMMPS is not reporting the pressure correctly (or you are
computing it the wrong way).

ii) LAMMPS is not choosing the correct volume. (There is a problem
with pressure regulation.)

I don't yet know the answer to the first problem (i). I saw several
posts about this topic recently, and when I find out, I will post a
followup message.

However I hope that the input script example file I posted this
morning solves the second problem (ii).

Unfortunately, I am not sure if the pressure in my liquid was equal to
my target pressure because in my tests the systems were initially
equilibrated at very low pressure (1atm). The target pressure could
have been 2atm, and I would not have noticed because the natural
fluctuations in pressure are so much larger (500atm).

In order to learn if "fix rigid" works at high pressures, tonight I
ran two additional simulations of water at high pressure=10000atm,
starting from the original volume (at 1 atm).

"solid-liquid" system. In this simulation, the atoms in the middle
layer were held stationary using "fix rigid" (as explained earlier.)

"pure-liquid" system. In this control system, all of the atoms were
allowed to move.

After pressure equilibration at the higher pressure, the pure-liquid
system compressed by 20%, and the liquid layer in the solid+liquid
system was compressed by a similar percentage ~20% of it's original
volume (+/- 5%). This makes me hopeful that "fix rigid" is working.

However, this is not a very accurate test. More testing is needed.

If I have time to setup more accurate high pressure tests, I will do
it and post the results here.

Cheers.

Andrew

Heng, if you have time, I'm curious to see what happens if you try
these 3 lines in your lammps script:
fix 1 immobile rigid single force * off off off torque * off off off
neigh_modify exclude group immobile immobile
fix 2 mobile npt temp 300.0 300.0 100.0 z 1.0 1.0 1000.0 dilate partial
(You may need to change the temperature and pressure.)

virial_immobile_atoms_sm.pdf (101 KB)

[Andrew_Jewett]

This is what I do in all of my input scripts and in my simulations.

(Sorry. I don't think this was clear in my last post.)

Cheers
Andrew

[Andrew_Jewett]

One more correction. I meant to say:

(i) LAMMPS is not reporting the pressure of the liquid correctly.
(This pressure should equal the target pressure.)

Cheers
Andrew

[fanheng]

I did almost exact the same. I used dialte partial and I used exclude to rule out the interaction between solid atoms. There is only one thing that is different, I fixed those atoms by not include them in any integration (Acturally I tried to use fix rigid and turn all the force and torque flags to off as you did, but I don't know why it never worked out in my DPD simulation)

I think I agree with you in most part. Maybe I did not state my concern clearly.
For my case, (I think for most of the cases involving liquid), the purpose of using pressure coupling is to regulate the pressure in the liquid phase, in other words I want the pressure in the bulk liquid phase to reach the target reference pressure. However, when LAMMPS is doing the pressure coupling, it's not targeting at the the pressure in the bulk liquid phase, but the total pressure in the entire system. Then here is my concern. Although we turn off the interactions in the solid, we can't turn off the interaction between liquid and solid. Maybe I'm wrong in this part but the local pressure at the liquid-solid interface may be different with the pressure in the liquid, which will bias the total pressure
. I mean if the local pressure at the liquid-solid interface is larger than reference pressure then the pressure in the liquid has to be lower to make the total pressure equal to the reference pressure.

Your reply reminds me another concern. How did LAMMPS measure the volume when the pressure is calculated? I assume the volume of the box will be used. But I think in your case, the correct volume should be the volume that only contain mobile phase.

Thanks a lot

Heng

[akohlmey]

. I mean if the local pressure at the liquid-solid interface is larger than reference pressure then the pressure in the liquid has to be lower to make the total pressure equal to the reference pressure.

so what? remember stress is inferred from forces, so your system
_wants_ to change where the stress is large (due to large forces),
so you have to include those interactions.

Your reply reminds me another concern. How did LAMMPS measure the volume when the pressure is calculated? I assume the volume of the box will be used. But I think in your case, the correct volume should be the volume that only contain mobile phase.

i think this is exactly the source of so many misunderstandings
*the* reason why i keep repeating that npt is not really the ensemble
that you want to use and that is suitable for your purposes. it can
approximate it, but will never be correct, since you effectively have
an ill-defined problem. an npt MD ensemble is meant to represent
the *bulk* behavior of a system through a small sample.
by using npt mixed with a rigid integrator, you are breaking this
requirement. so it cannot work perfectly and it will get worse
the larger the deviation from the ideal scenario is.

ignoring the forces (and stress tensor contributions) from the
interactions inside the rigid body is a must (they are ignored
by the rigid integrator, too), however, the interaction between
the rigid objects and the rest *has* to be included and this is
where you have a source of error, that you cannot remove.
to avoid spurious motions and scaling of the rigid particles,
you have to use dilate partial, but that will also result in
spurious forces due to rescaling the positions of the npt
group particles while the rigid particles remain in place.
what you _really_ want is to have the rigid object particles
moved _with_ the rest, but only as a center of mass motions.
and, yet, even then you have a small error due to the
rigid object not expanding or shrinking at all while the
rest does. these effects become large the larger the
fraction of rigid particles is and the large the system becomes.

what you *actually* want to model is more like a piston
type setup, but with surface effects removed through pbc.
an npt ensemble is not that. there is a small but significant
difference there. the advantage of the piston type setup is
that you don't look at the average pressure of your total
system, but at the force per area on the movable part of
the piston and thus don't have to worry about all the details
that give you trouble with npt+rigid including not having
to expand or shrink your system.

now after all the bad news, the question is: what to do?
a) the simple way out: run a bigger system, do an initial
relaxation (with all the "tricks" discussed) and then switch
to nvt and stop worrying. the larger the system, the lower
the need to couple to size external fluctuation to properly
sample phase space, since larger fluctuations can happen
in larger system cells
b) pick the variant that has the least error and ignore them.
not what you generally want, but if you can quantify the
error inferred, it may not be too bad.
c) the hard way out: write a mixed npt atomistic/rigid integrator
that takes the intricacies of those systems into account
and removes the errors from dilate partial.
d) your ideas...

cheers,
     axel.

[Laurent_Joly]

Dear Axel and all,

I am very interested by this topic as I am just trying to investigate
a similar system (liquid reservoirs at a given pressure separated by
rigid or frozen membranes). At one point Andrew suggested to "add
[very stiff] harmonic bonds connecting the atoms you want to
immobilise together at the correct distances" to mimic the effect of
fix rigid, and apply fix npt to all atoms. Equivalently one could use
fix spring/self to mimic the effect of not integrating the equation of
motion for the frozen part of the system. This seems a very reasonable
approach to me. Do you see any flaw with this procedure?

Best,
Laurent

[akohlmey]

dear laurent,

Dear Axel and all,

I am very interested by this topic as I am just trying to investigate

it is a very interesting discussion (i wish there were more of this
kind on the list and less of the "my input doesn't work and i am
to lazy to debug it" questions), as it forces us to re-evaluate the
the very basics and intricacies of the methods that we often use
without questioning.

it also seems a very timely discussion, since there appears
to be a growing demand for running simulations of (hard) nanoparticles
dispersed in some kind of soft media.

a similar system (liquid reservoirs at a given pressure separated by
rigid or frozen membranes). At one point Andrew suggested to "add
[very stiff] harmonic bonds connecting the atoms you want to
immobilise together at the correct distances" to mimic the effect of
fix rigid, and apply fix npt to all atoms. Equivalently one could use
fix spring/self to mimic the effect of not integrating the equation of
motion for the frozen part of the system. This seems a very reasonable
approach to me. Do you see any flaw with this procedure?

as is using fix rigid. the trouble arises from mixing it with fix npt
and wanting to apply fix npt to only a part of the system.
both using stiff harmonic bonds or fix spring/self will suffer
from the same problems of creating spurious forces when
combined with fix npt. if you use dilate all, you get too large
"internal" forces (and thus stress) from stretching the bond
and you'll have to massively reduce the time step, if you use
fix spring/self and dilate partial, you still get the spurious
non-bonded forces between the fixed and unfixed particles
from the position rescaling.

i personally think that the cleanest way to handle this
kind of setup would be to go non-periodic, but that is
very expensive in terms of computer time. if there are
no long-range electrostatics involved, one could conceive
trying to use a "two-layer" movable wall, i.e. a wall with
one layer of fixed particles and a second that is attached
to the first with springs and that then would be strongly
thermalized and thus provide coupling to bulk. pressure
would then be computed as "macroscopic pressure"
from the total force on the wall per area. but i have not
thought that through entirely, so this approach may
have some skeletons hidden somewhere, too.

cheers,
     axel.

[Laurent_Joly]

Dear Axel,

i personally think that the cleanest way to handle this
kind of setup would be to go non-periodic, but that is
very expensive in terms of computer time. if there are
no long-range electrostatics involved, one could conceive
trying to use a "two-layer" movable wall, i.e. a wall with
one layer of fixed particles and a second that is attached
to the first with springs and that then would be strongly
thermalized and thus provide coupling to bulk. pressure
would then be computed as "macroscopic pressure"
from the total force on the wall per area. but i have not
thought that through entirely, so this approach may
have some skeletons hidden somewhere, too.

Alas in my case there is long-range electrostatics involved, as I'm
simulating aqueous electrolytes. However your suggestion of going
non-periodic is interesting, as I only wish to reach the required
pressure in the liquid, and then work at constant volume. So I could
could use frozen-liquid-water-pistons to put the system at the desired
pressure, then cut a box between the pistons, and go back to periodic
boundary conditions and work at constant volume.

Another possibility (with pbc) would be to measure the pressure
locally in the bulk region of the liquid and use the result as a
control parameter for a loop over the box size, using fix deform...

Best,
Laurent

[akohlmey]

dear laurent,

Dear Axel,

> i personally think that the cleanest way to handle this
> kind of setup would be to go non-periodic, but that is
> very expensive in terms of computer time. if there are
> no long-range electrostatics involved, one could conceive
> trying to use a "two-layer" movable wall, i.e. a wall with
> one layer of fixed particles and a second that is attached
> to the first with springs and that then would be strongly
> thermalized and thus provide coupling to bulk. pressure
> would then be computed as "macroscopic pressure"
> from the total force on the wall per area. but i have not
> thought that through entirely, so this approach may
> have some skeletons hidden somewhere, too.

Alas in my case there is long-range electrostatics involved, as I'm
simulating aqueous electrolytes. However your suggestion of going
non-periodic is interesting, as I only wish to reach the required
pressure in the liquid, and then work at constant volume. So I could
could use frozen-liquid-water-pistons to put the system at the desired
pressure, then cut a box between the pistons, and go back to periodic
boundary conditions and work at constant volume.

hmmmm... i suspect this will put you a bit off-target.
you'll be missing the long-range contributions,
particularly near the walls. since water is not very
compressible, your pressure will jump up by quite a
bit as soon as you go back to periodic.

also you'd have to use a cutoff that covers the entire
system to avoid artifacts and that may make it intractable.

Another possibility (with pbc) would be to measure the pressure
locally in the bulk region of the liquid and use the result as a
control parameter for a loop over the box size, using fix deform...

if you do have a sufficiently large bulk region,
that could work. if you are planning to do production
in fixed volume, then using any of the npt+rigid routes,
may not be as bad. if you can minimize the impact
of the unwanted forces you could pre-equilibrate
in the "automatic" way and then perhaps only need
some small "manual" refinements as you describe.

long-range electrostatics make this kind of system
really messy since you also have to worry about
charges in the rigid part of your system. in the
kspace part you cannot exclude them!

salut,
   axel.

[Andrew_Jewett]

You are right that lammps regulates the volume of the simulation by
attempting to set the pressure of the whole system to the target
pressure. (Not just the liquid.) I think I can address these points
in detail soon.

It is amazing how many things that can go wrong in a lammps input
script. I also had some difficulty using "fix rigid". In my case it
is because I made a mistake using the "group" command to define the
"mobile" and "immobile" groups. (Be careful using "group ... type i
j" command.)
I was able to find my own bug by simplifying my system, throwing away
things like "shake", and reducing the number of atom types and
molecule types. In your case, it looks like you have a barrier
penetrated by a nanotube, and tip3p water. There are many ways you can
simplify this.)

I'm setting up some new tests related to this question. When they are
finished, I'll try to give you a better reply.

Cheers
Andrew

[Andrew_Jewett]

This is interesting. I did not know about fix spring. From a casual
glance at the lammps source code, it appears that "fix spring" behaves
in a similar way to "dilate partial", because I don't think the
positions of the tether points are ever re-scaled during pressure
equilibration. (Please correct me if I'm wrong.)

(In contrast, in AMBER, harmonic anchor points are rescaled as the box
size changes. This is not always desirable.)

My only concern with using fix spring is that I'm not sure if the
spring's forces on the immobilized atoms are included in the virial.
(In the special case that you are using extermal forces like "fix
spring", instead of rigid constraints, like "fix rigid", then I THINK
we have to make sure that forces from those springs are also included
in the virial, along with all the other forces acting on the frozen
atoms. I could not tell if that his happening. If I had to guess, I
would guess no. Perhaps somebody else knows the answer? Perhaps I
can test this...)
Anyone know how does "fix spring" behave under NPT conditions?

Cheers
Andrew

[akohlmey]

This is interesting. I did not know about fix spring. From a casual
glance at the lammps source code, it appears that "fix spring" behaves
in a similar way to "dilate partial", because I don't think the
positions of the tether points are ever re-scaled during pressure
equilibration. (Please correct me if I'm wrong.)

you are correct. the tether positions are defined
when the fix is created and not changed.

(In contrast, in AMBER, harmonic anchor points are rescaled as the box
size changes. This is not always desirable.)

My only concern with using fix spring is that I'm not sure if the
spring's forces on the immobilized atoms are included in the virial.

they are not. but you seem to overlook the bigger problem
which results from scaling some positions but not the others.
with having fix npt applied to a rather incompressible
matter you generate a lot of spurious forces and stress
due to the (partial) rescaling of positions.

(In the special case that you are using extermal forces like "fix
spring", instead of rigid constraints, like "fix rigid", then I THINK
we have to make sure that forces from those springs are also included
in the virial, along with all the other forces acting on the frozen
atoms. I could not tell if that his happening. If I had to guess, I
would guess no. Perhaps somebody else knows the answer? Perhaps I
can test this...)

Anyone know how does "fix spring" behave under NPT conditions?

yes. it knows nothing about it.

axel.

[Laurent_Joly]

> i personally think that the cleanest way to handle this
> kind of setup would be to go non-periodic, but that is
> very expensive in terms of computer time. if there are
> no long-range electrostatics involved, one could conceive
> trying to use a "two-layer" movable wall, i.e. a wall with
> one layer of fixed particles and a second that is attached
> to the first with springs and that then would be strongly
> thermalized and thus provide coupling to bulk. pressure
> would then be computed as "macroscopic pressure"
> from the total force on the wall per area. but i have not
> thought that through entirely, so this approach may
> have some skeletons hidden somewhere, too.

Alas in my case there is long-range electrostatics involved, as I'm
simulating aqueous electrolytes. However your suggestion of going
non-periodic is interesting, as I only wish to reach the required
pressure in the liquid, and then work at constant volume. So I could
could use frozen-liquid-water-pistons to put the system at the desired
pressure, then cut a box between the pistons, and go back to periodic
boundary conditions and work at constant volume.

hmmmm... i suspect this will put you a bit off-target.
you'll be missing the long-range contributions,
particularly near the walls. since water is not very
compressible, your pressure will jump up by quite a
bit as soon as you go back to periodic.

also you'd have to use a cutoff that covers the entire
system to avoid artifacts and that may make it intractable.

Here is what I am planning to do: confine my mixed solid-liquid,
charge neutral system between two charge neutral frozen-liquid-water
walls. The system being periodic in the x and y directions, and
non-periodic in the z direction (with kspace modify slab 3.0), apply a
given force f=p*lx*ly to one wall with fix aveforce, and check that
the force acting on the opposite wall is p*lx*ly with the variable
function fcm(group,z). It seems to me that as long as each wall and
the confined system are charge neutral, and that there is enough
liquid so that there is a bulk liquid between the pistons and the
solid part of the confined system, this procedure should provide the
correct pressure inside the liquid.

But I agree that when turning back to pbc, I will slightly change the
density and that this may result in large pressure change. Hopefully
what I am trying to look at is not too sensitive to the system
pressure, and anyway the uncertainty on the pressure resulting from
this procedure can be evaluated.

I will also try to use pbc and fix spring/self, and see if I can find
a reasonable compromise between the stiffness of the springs and the
timestep.

Another possibility (with pbc) would be to measure the pressure
locally in the bulk region of the liquid and use the result as a
control parameter for a loop over the box size, using fix deform...

if you do have a sufficiently large bulk region,
that could work. if you are planning to do production
in fixed volume, then using any of the npt+rigid routes,
may not be as bad. if you can minimize the impact
of the unwanted forces you could pre-equilibrate
in the "automatic" way and then perhaps only need
some small "manual" refinements as you describe.

long-range electrostatics make this kind of system
really messy since you also have to worry about
charges in the rigid part of your system. in the
kspace part you cannot exclude them!

I just realized that this long-range part is indeed a pain because it
prevents me from computing the pressure locally, at least with the
current version of lammps (if I understood correctly the many
discussions on this topic that there has been recently...)

A bientôt,
Laurent

[akohlmey]

long-range electrostatics make this kind of system
really messy since you also have to worry about
charges in the rigid part of your system. in the
kspace part you cannot exclude them!

I just realized that this long-range part is indeed a pain because it
prevents me from computing the pressure locally, at least with the
current version of lammps (if I understood correctly the many
discussions on this topic that there has been recently...)

if you do a piston type setup, you don't have to compute
the local pressure (whatever that would be), but rather
only monitor the average force due to the liquid phase
on your piston wall. you'll basically switch from a
microscopic measuring of pressure via the stress tensor
to a macroscopic approach, which has none of the
technical problems, the other approach has.

axel.

[Laurent_Joly]

I just realized that this long-range part is indeed a pain because it
prevents me from computing the pressure locally, at least with the
current version of lammps (if I understood correctly the many
discussions on this topic that there has been recently...)

if you do a piston type setup, you don't have to compute
the local pressure (whatever that would be), but rather
only monitor the average force due to the liquid phase
on your piston wall. you'll basically switch from a
microscopic measuring of pressure via the stress tensor
to a macroscopic approach, which has none of the
technical problems, the other approach has.

I totally agree. In this paragraph I was referring to my suggestion of
measuring the pressure locally in the bulk liquid while keeping pbc.
This is not possible with long-range electrostatic, so I'll definitely
go for the piston type setup.

Best,
Laurent