Combination of Setforce and NPT

Dear Developers and Users

I’ve simulated a system consists of a liquid placed on a mineral crystal and that structure was such that the upper side of the liquid contacting the lower bond of the crystal as a result of PBC. I tried to adjust the liquid pressure at the desired value by NPT ensemble only acting through z-axis, i.e., normal to the mineral surface and the liquids. However, the pressure does not converge to the set point and also, the density of the liquid is not equivalent to the experimental one. I already repeated the simulation at the same conditions without mineral and both pressures and liquid density achieved the correct values. Note that I used setforce to keep crystal structure rigid during the simulation.

Has anybody ever experienced such case? I guess the setforce and NPT could not act together. Is it right?

In this respect, I realized that mineral structure has slightly moved in contrary to the setforce.

Is it a bug?

Bests

What do you mean by “rigid”? If you apply setforce on the crystal, it will be frozen, and it will not take part to the dynamics of your system. It will behave more or less like a wall. The best way to check your system, is too output the force per atom in your dump. Then you can check if your atom are set to the force you want. I assume not. Julien

Dear Developers and Users

I’ve simulated a system consists of a liquid placed on a mineral crystal and that structure was such that the upper side of the liquid contacting the lower bond of the crystal as a result of PBC. I tried to adjust the liquid pressure at the desired value by NPT ensemble only acting through z-axis, i.e., normal to the mineral surface and the liquids. However, the pressure does not converge to the set point and also, the density of the liquid is not equivalent to the experimental one. I already repeated the simulation at the same conditions without mineral and both pressures and liquid density achieved the correct values. Note that I used setforce to keep crystal structure rigid during the simulation.

What do you mean by “rigid”? If you apply setforce on the crystal, it will be frozen, and it will not take part to the dynamics of your system. It will behave more or less like a wall.

when, oh when will people learn the difference between rigid and immobile?

Has anybody ever experienced such case? I guess the setforce and NPT could not act together. Is it right?

In this respect, I realized that mineral structure has slightly moved in contrary to the setforce.

The best way to check your system, is too output the force per atom in your dump. Then you can check if your atom are set to the force you want.

Is it a bug?

I assume not.

it is documented behavior. what has to be considered in this case as well, is that by default fix npt will apply the cell deformation due to pressure to all atoms (dilate option). in this case, immobile does not mean immobile in absolute coordinates. the major issue here is that using fix npt and immobilizing atoms at the same time is a conceptually problematic choice in the first place, since it does not give the expected behavior in the default case, but also applying the deformation to only the mobile atoms is problematic. it may indeed be preferable to have the “immobile” part not exactly immobile but kept rigid through using fix rigid instead.

axel.

In this respect, I realized that mineral structure has slightly moved in contrary to the setforce.

The best way to check your system, is too output the force per atom in your dump. Then you can check if your atom are set to the force you want.

Is it a bug?

I assume not.

it is documented behavior. what has to be considered in this case as well, is that by default fix npt will apply the cell deformation due to pressure to all atoms (dilate option). in this case, immobile does not mean immobile in absolute coordinates. the major issue here is that using fix npt and immobilizing atoms at the same time is a conceptually problematic choice in the first place, since it does not give the expected behavior in the default case, but also applying the deformation to only the mobile atoms is problematic. it may indeed be preferable to have the “immobile” part not exactly immobile but kept rigid through using fix rigid instead.

Exactly! In such situation, I always considered the “immobile” atoms to be immobile in the box referential, and was never bother by changes in their absolute coordinate.
When I first worked with IMD (after working with Lammps), I was shocked by the possibility to apply an NPT ensemble in non-periodic directions. Actually, it was making sense in one situation only, related to the current thread, that is some kind of force-controlled boundary conditions. In lammps, the same conditions can be achieve, by just using PBC and some vacuum.

Julien