TIP5P with OPLSAA ethanol mixure, error:Non-numerical pressure-simulation unstable

LAMMPS LAMMPS Mailing List Mirror
1

Dear lammps user

I am using lammps (Stable version (13 Mar 2018) ) to simulate viscosity of binary mixture of water (TIP5P, fix/rigid/small molecule) and ethanol (OPLSAA force field). The box are established by Moltemplate.sh. The lammps come cross the error on the first 100 step : Non-numerical pressure-simulation unstable. By the way, when i use the same simulation settings to test the pure TIP5P water model, Lammps can run successfully and no errors. The main script are listed as follows, I appreciate if any one can give me some advice.

timestep 0.2

group H2O type 73 74 75
fix rigid H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000

group ethanol type 80 85 96 97 99
fix npt ethanol npt temp 298 298 100 iso 1 1 1000

best regards
cheng

2

Dear lammps user

I am using lammps (Stable version (13 Mar 2018) ) to simulate viscosity of binary mixture of water (TIP5P, fix/rigid/small molecule) and ethanol (OPLSAA force field). The box are established by Moltemplate.sh. The lammps come cross the error on the first 100 step : Non-numerical pressure-simulation unstable. By the way, when i use the same simulation settings to test the pure TIP5P water model, Lammps can run successfully and no errors. The main script are listed as follows, I appreciate if any one can give me some advice.

please study your simulation output carefully. you *must not* have two
fixes that modify the box dimensions. LAMMPS should report this issue
in the output. it is discussed in the manual, too.

axel.

3

thanks for your reply, and after searched the previous mailing list (https://lammps.sandia.gov/threads/msg33197.html ) and manual. i found my problem is about " performing a hybrid simulation in NPT with water atoms in rigid body, and ethanol not". So I tried the advice in the above mailing list and modified my input file to “fix npt ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol, fix nvt H2O rigid/nvt/small molecule temp 298 298 100”. and the error is also, high temperature peaks appears even after long equilibration runs and the system failed at last with bond atom missing. So i am confused about how to perform the hybrid simulation with water atom in rigid body and ethanol in flexibel in Lammps

best regards
Yalishanda

------------------ Original ------------------

4

thanks for your reply, and after searched the previous mailing list
(https://lammps.sandia.gov/threads/msg33197.html ) and manual. i found my
problem is about " performing a hybrid simulation in NPT with water atoms in
rigid body, and ethanol not". So I tried the advice in the above mailing
list and modified my input file to "fix npt ethanol npt temp 298 298 100
iso 1 1 1000 dilate ethanol, fix nvt H2O rigid/nvt/small molecule temp 298
298 100". and the error is also, high temperature peaks appears even after
long equilibration runs and the system failed at last with bond atom
missing. So i am confused about how to perform the hybrid simulation with
water atom in rigid body and ethanol in flexibel in Lammps

best regards
Yalishanda

Dear Yalishanda

    These kinds of hybrid simulations are complicated to run. I
honestly don't know how to run these hybrid simulations in LAMMPS
correctly. The best help I can provide is to send you a list of
alternative methods you can try. If someone else has an opinion,
please post a reply. If you figure out which of these methods works,
please post again to the list.

As you know, Axel is referring to these two lines from your input file:

fix rigid H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000
fix npt ethanol npt temp 298 298 100 iso 1 1 1000

Both of these lines will update the simulation box size in order to
preserve the pressure. You must change your input file so that only
one of these commands controls the box size. If my understanding is
correct, then it does not matter whether you control the pressure
using the fix which moves the rigid H2O molecules, or using the fix
which integrates other flexible molecules in your system. (From the
documentation: "Regardless of what atoms are in the fix group (the
only atoms which are time integrated), a global pressure or stress
tensor is computed for all atoms.")

Here are some alternatives to your approach. I'm not sure which one
of these approaches will produce the correct physical behavior. You
will have to try each, and compare the behavior.

There are multiple reasons a simulation can crash. Try the simplest
simulation first and verify it does not crash. "Method 1" and "Method
2" are constant-volume simulations (not constant pressure). Please
try them and verify that your simulation does not crash before adding
pressure regulation (methods 3-6):

Method 1: NVE simulation. This should run without crashing.

velocity all create 298 4928459
fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000

Method 2: NVT simulation. This should run without crashing.

fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

I personally prefer to use "fix langevin" to regulate the temperature
(instead the Nose-Hoover based thermostats you are using). Langevin
dynamics has a much simpler implementation than Nose-Hoover. "fix
langevin" simply applies a random force to the particles similar to
collisions they would experience when immersed in a heat bath. Both
"fix nvt" or "fix npt" use the Nose-Hoover algorithm. I confess I
have forgotten how Nose-Hoover thermostats work, but I worry that
strange behavior could occur when you are applying a Nose-Hoover
thermostat to PART of your system. (Not the whole system. If you
apply Langevin to part of your system, nothing too strange should
happen.)

https://lammps.sandia.gov/doc/fix_langevin.html

If methods 1 and 2 are working, then continue:

Method 3: Apply the pressure regulation to the H2O molecules. Apply a
random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/nph/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

Method 4: Apply the pressure regulation to the ethanol molecules.
Apply a random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nve iso 1 1 1000 dilate mobile
fix fxLangevin all langevin 298 298 100.0 48279

Method 5: Regulate the pressure of the system using the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nvt temp 298 298 100

https://lammps.sandia.gov/doc/fix_rigid.html
https://lammps.sandia.gov/doc/fix_nve.html

Method 6: Regulate the pressure of the system using the ethanol molecules

You tried this already?

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

https://lammps.sandia.gov/doc/fix_rigid.html

  If you tried this and it did not work, there are other reasons a
simulation can crash. In that case, probably "Method 1" (described
earlier) should also crash. So try method 1 first.

   HOPEFULLY one of these methods work (method 3-6).
   If not, then try method 7 and 8 below.

Is it correct to apply a Nose-Hoover thermostat to both the H2O
molecules and the ethanol molecules? I honestly don't remember. The
next 2 methods apply temperature regulation to EITHER H2O OR ethanol.
send your
Method 7: Regulate the temperature and pressure of the system using
the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nve

compute tempH20 H20 temp
fix_modify fxH20 temp tempH20

Method 8: Regulate the temperature and pressure of the system using
the ethanol molecules

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

I also suggest you create a version TIP5P which is ALMOST rigid.
Use strong, stiff bonds and angles to preserve the shape of the TIP5P
molecule (at least 2000.0), and reduce the timestep to 0.05fs or less
(to prevent numerical instability resulting from these stiff bonds).
You will not use this in your final simulations, because it is
inefficient to use such a small timestep. However you can compare the
simulations you run using methods 3-6 with a simulation using no rigid
molecules, to see which of these methods reproduces the correct
potential energy and volume.

I am sorry that I cannot tell you which method to use.

Good luck.

Andrew

P.S.
    TIP5P is a water model I would love to include with moltemplate.
If you get this working, would you be willing to submit "tip5p.lt"
file so I can put it on the moltemplate web page and git download
page? (I can add your name to the list of contributors if you like.)

5

I corrected some mistakes below

thanks for your reply, and after searched the previous mailing list
(https://lammps.sandia.gov/threads/msg33197.html ) and manual. i found my
problem is about " performing a hybrid simulation in NPT with water atoms in
rigid body, and ethanol not". So I tried the advice in the above mailing
list and modified my input file to "fix npt ethanol npt temp 298 298 100
iso 1 1 1000 dilate ethanol, fix nvt H2O rigid/nvt/small molecule temp 298
298 100". and the error is also, high temperature peaks appears even after
long equilibration runs and the system failed at last with bond atom
missing. So i am confused about how to perform the hybrid simulation with
water atom in rigid body and ethanol in flexibel in Lammps

best regards
Yalishanda

Dear Yalishanda

    These kinds of hybrid simulations are complicated to run. I
honestly don't know how to run these hybrid simulations in LAMMPS
correctly. The best help I can provide is to send you a list of
alternative methods you can try. If someone else has an opinion,
please post a reply. If you figure out which of these methods works,
please post again to the list.

As you know, Axel is referring to these two lines from your input file:

fix rigid H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000
fix npt ethanol npt temp 298 298 100 iso 1 1 1000

Both of these lines will update the simulation box size in order to
preserve the pressure. You must change your input file so that only
one of these commands controls the box size. If my understanding is
correct, then it does not matter whether you control the pressure
using the fix which moves the rigid H2O molecules, or using the fix
which integrates other flexible molecules in your system. (From the
documentation: "Regardless of what atoms are in the fix group (the
only atoms which are time integrated), a global pressure or stress
tensor is computed for all atoms.")

Here are some alternatives to your approach. I'm not sure which one
of these approaches will produce the correct physical behavior. You
will have to try each, and compare the behavior.

There are multiple reasons a simulation can crash. Try the simplest
simulation first and verify it does not crash. "Method 1" and "Method
2" are constant-volume simulations (not constant pressure). Please
try them and verify that your simulation does not crash before adding
pressure regulation (methods 3-6):

Method 1: NVE simulation. This should run without crashing.

velocity all create 298 4928459
fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000

Method 2: NVT simulation. This should run without crashing.

fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

I personally prefer to use "fix langevin" to regulate the temperature
(instead the Nose-Hoover based thermostats you are using). Langevin
dynamics has a much simpler implementation than Nose-Hoover. "fix
langevin" simply applies a random force to the particles similar to
collisions they would experience when immersed in a heat bath. Both
"fix nvt" or "fix npt" use the Nose-Hoover algorithm. I confess I
have forgotten how Nose-Hoover thermostats work, but I worry that
strange behavior could occur when you are applying a Nose-Hoover
thermostat to PART of your system. (Not the whole system. If you
apply Langevin to part of your system, nothing too strange should
happen.)

https://lammps.sandia.gov/doc/fix_langevin.html

If methods 1 and 2 are working, then continue:

Method 3: Apply the pressure regulation to the H2O molecules. Apply a
random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/nph/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

Method 4: Apply the pressure regulation to the ethanol molecules.
Apply a random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nve iso 1 1 1000 dilate mobile

I'm sorry. That was a mistake. That should be:

fix fxEthanol ethanol nve iso 1 1 1000 dilate ethanol

fix fxLangevin all langevin 298 298 100.0 48279

Method 5: Regulate the pressure of the system using the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nvt temp 298 298 100

https://lammps.sandia.gov/doc/fix_rigid.html
https://lammps.sandia.gov/doc/fix_nve.html

Method 6: Regulate the pressure of the system using the ethanol molecules

You tried this already?

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

https://lammps.sandia.gov/doc/fix_rigid.html

  If you tried this and it did not work, there are other reasons a
simulation can crash. In that case, probably "Method 1" (described
earlier) should also crash. So try method 1 first.

   HOPEFULLY one of these methods work (method 3-6).
   If not, then try method 7 and 8 below.

Is it correct to apply a Nose-Hoover thermostat to both the H2O
molecules and the ethanol molecules? I honestly don't remember. The
next 2 methods apply temperature regulation to EITHER H2O OR ethanol.

Method 7: Regulate the temperature and pressure of the system using
the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nve

compute tempH20 H20 temp
fix_modify fxH20 temp tempH20

Another mistake. That should be:

compute tempH2O H2O temp
fix_modify fxH2O temp tempH2O

These last two commands might be unnecessary.
Try the simulation without these commands. If the potential energy or
volume of the system is not correct, then add these two commands.

Method 8: Regulate the temperature and pressure of the system using
the ethanol molecules

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

The next two commands might also be helpful.

compute tempEthanol ethanol temp
fix_modify fxEthanol temp tempEthanol

Again, if you solve your problem, please reply and tell us which
method worked the best.

6

Many more mistakes. I was so sloppy!
Here are more corrections...

I corrected some mistakes below

thanks for your reply, and after searched the previous mailing list
(https://lammps.sandia.gov/threads/msg33197.html ) and manual. i found my
problem is about " performing a hybrid simulation in NPT with water atoms in
rigid body, and ethanol not". So I tried the advice in the above mailing
list and modified my input file to "fix npt ethanol npt temp 298 298 100
iso 1 1 1000 dilate ethanol, fix nvt H2O rigid/nvt/small molecule temp 298
298 100". and the error is also, high temperature peaks appears even after
long equilibration runs and the system failed at last with bond atom
missing. So i am confused about how to perform the hybrid simulation with
water atom in rigid body and ethanol in flexibel in Lammps

best regards
Yalishanda

Dear Yalishanda

    These kinds of hybrid simulations are complicated to run. I
honestly don't know how to run these hybrid simulations in LAMMPS
correctly. The best help I can provide is to send you a list of
alternative methods you can try. If someone else has an opinion,
please post a reply. If you figure out which of these methods works,
please post agai>>>>n to the list.

As you know, Axel is referring to these two lines from your input file:

fix rigid H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000
fix npt ethanol npt temp 298 298 100 iso 1 1 1000

Both of these lines will update the simulation box size in order to
preserve the pressure. You must change your input file so that only
one of these commands controls the box size. If my understanding is
correct, then it does not matter whether you control the pressure
using the fix which moves the rigid H2O molecules, or using the fix
which integrates other flexible molecules in your system. (From the
documentation: "Regardless of what atoms are in the fix group (the
only atoms which are time integrated), a global pressure or stress
tensor is computed for all atoms.")

Here are some alternatives to your approach. I'm not sure which one
of these approaches will produce the correct physical behavior. You
will have to try each, and compare the behavior.

There are multiple reasons a simulation can crash. Try the simplest
simulation first and verify it does not crash. "Method 1" and "Method
2" are constant-volume simulations (not constant pressure). Please
try them and verify that your simulation does not crash before adding
pressure regulation (methods 3-6):

Method 1: NVE simulation. This should run without crashing.

velocity all create 298 4928459
fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000

This should be:

velocity all create 298 4928459
fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nve

Method 2: NVT simulation. This should run without crashing.

fix fxH2O H2O rigid/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

Another mistake. This should be:

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nve
fix fxLanvegin all langevin 298 298 100.0 48279

I personally prefer to use "fix langevin" to regulate the temperature
(instead the Nose-Hoover based thermostats you are using). Langevin
dynamics has a much simpler implementation than Nose-Hoover. "fix
langevin" simply applies a random force to the particles similar to
collisions they would experience when immersed in a heat bath. Both
"fix nvt" or "fix npt" use the Nose-Hoover algorithm. I confess I
have forgotten how Nose-Hoover thermostats work, but I worry that
strange behavior could occur when you are applying a Nose-Hoover
thermostat to PART of your system. (Not the whole system. If you
apply Langevin to part of your system, nothing too strange should
happen.)

https://lammps.sa>>ndia.gov/doc/fix_langevin.html

If methods 1 and 2 are working, then continue:

Method 3: Apply the pressure regulation to the H2O molecules. Apply a
random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/nph/small molecule iso 1 1 1000
fix fxEthanol ethanol nve iso 1 1 1000
fix fxLanvegin all langevin 298 298 100.0 48279

This should be:

fix fxH2O H2O rigid/nph/small molecule iso 1 1 1000 #"dilate ethanol"?
fix fxEthanol ethanol nve
fix fxLanvegin all langevin 298 298 100.0 48279

I am not sure if you need to add "dilate ethanol" to the end of "fix
fxH2O H2O rigid/nph/small..."
(I suspect not. Try it without "dilate ethanol" first.)

Method 4: Apply the pressure regulation to the ethanol molecules.
Apply a random force to both H2O and ethanol (using fix langevin).

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nve iso 1 1 1000 dilate mobile

I'm sorry. That was a mistake. That should be:

fix fxEthanol ethanol nve iso 1 1 1000 dilate ethanol

fix fxLangevin all langevin 298 298 100.0 48279

Method 4 should be:

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol nph iso 1 1 1000 dilate ethanol
fix fxLangevin all langevin 298 298 100.0 48279

I am not sure if you need to add "dilate ethanol" to the end of "fix
fxEthanol ethanol nph ..."
(I suspect you do. Try it with "dilate ethanol" first.)

Method 5: Regulate the pressure of the system using the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nvt temp 298 298 100

This should be:

fix fxH2O H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000
#"dilate ethanol"?
fix fxEthanol ethanol nvt temp 298 298 100

https://lammps.sandia.gov/doc/fix_rigid.html
https://lammps.sandia.gov/doc/fix_nve.html

Method 6: Regulate the pressure of the system using the ethanol molecules

You tried this already?

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

This looks correct. I don't know why it did not work for you.

https://lammps.sandia.gov/doc/fix_rigid.html

  If you tried this and it did not work, there are other reasons a
simulation can crash. In that case, probably "Method 1" (described
earlier) should also crash. So try method 1 first.

   HOPEFULLY one of these methods work (method 3-6).
   If not, then try method 7 and 8 below.

Is it correct to apply a Nose-Hoover thermostat to both the H2O
molecules and the ethanol molecules? I honestly don't remember. The
next 2 methods apply temperature regulation to EITHER H2O OR ethanol.

Method 7: Regulate the temperature and pressure of the system using
the H2O molecules

fix fxH2O H2O rigid/nvt/small molecule temp 298 298 100 iso 1 1 1000
fix fxEthanol ethanol nve

This should be:

fix fxH2O H2O rigid/npt/small molecule temp 298 298 100 iso 1 1 1000
#"dilate ethanol"?
fix fxEthanol ethanol nve

compute tempH20 H20 temp
fix_modify fxH20 temp tempH20

Another mistake. That should be:

compute tempH2O H2O temp
fix_modify fxH2O temp tempH2O

These last two commands might be unnecessary.
Try the simulation without these commands. If the potential energy or
volume of the system is not correct, then add these two commands.

Method 8: Regulate the temperature and pressure of the system using
the ethanol molecules

fix fxH2O H2O rigid/small molecule
fix fxEthanol ethanol npt temp 298 298 100 iso 1 1 1000 dilate ethanol

This looks correct.

I am very sorry. I was in a hurry when I typed my original message.
There may be more mistakes. But hopefully you understand my intentions.

I am not sure if my emails are helpful at all.
Good luck. Please let us know what method works.

Andrew

7

My preference for these kind of systems is to use the variable cell (barostat) integrator for the rigid particles and - unless there is a good reason not to - let it dilate group all. Partial dilation can easily create overlaps, and using variable cell with the rigid fix does the dilation on the rigid bodies, not on individual atoms (which is likely to be ignored by the rigid fix anyway, as it moves and maintains the rigid bodies across the run.