Question occures during the use of fix npt

Dear Lammps users/developers,

I want to simulate a particle containing formic acid and water, and the NPT ensemble is used. The formic acid molecule is fully rigid. I use the method recommended in the manual. Therefore, two kinds of “in” files are listed:

fix a2 formic rigid/npt molecule temp 180 180 100 iso 1.0 1.0 1000 dilate all

fix shk water shake 0.0001 50 100 b 5 a 5 t 6 7

fix b2 water nvt temp 180 180 100.0

fix shk water shake 0.0001 50 100 b 5 a 5 t 6 7

fix b2 water npt temp 180 180 100.0 iso 1.0 1.0 1000 dilate all

fix a2 formic rigid/nvt molecule temp 180 180 100

I get two different results, of which the first crashes after hundreds of thousands of steps, and the second method can operate normally. The change of system’s pressure have been tracked and stated as follow.

DDD1C82E@...5901....jpgDDD1992E@...5901....jpg

a. Method ① b. Method ②

When I tack the volume of the system, I have found that the volume of method ① decreases monotonically from a=8nm to a=3.8nm. Then the simulation crashes. In addition, when the Pdamp is changed from 1000 to 2000, the result remains the same. However, the volume of method ② increases from a=8nm to about a=22nm, and the error is not occurring.

I am confused by two question: ⑴ why results of these two methods is different, I have learnt from the manual that a global pressure or stress tensor is computed for all atoms, regardless of what atoms are in the fix group.⑵ why the simulation crashes when method ① is used. I want to know whether the reason is the wrong use of fix npt command in the “in” file or the wrong physical model.

I would appreciate any help and suggestion.

Thanks!

Zhang Chao

  1. Barostats can only be counted on to work in simple and/or commonly encountered use cases. For the rest, it is often easier to make do without a barostat. Barostats are only strictly necessary when N is small, where NVT will not sample the physically-important range of density fluctuations.

  2. I think the problem in the first script is that the rigid/npt fix is also acting on the water molecules, and this is interfering with the action of fix shake.

  3. In the second script, fix npt and fix shake are acting on the same group and so play nice together.

  4. You may still have a problem with the formic acid molecule being dilated by fix npt. I f that is unacceptable, you can fix it using dilate water keyword, but that may cause other problems.

  5. You can probably use fix rigid/nve instead of rigid/nvt, as the rigid molecule will be thermostated by the water molecules.

DDD1992E@...5901....jpg

DDD1C82E@...5901....jpg