What's wrong with the blowing-up pressure after restart?

Dear Everyone,

I have a system with water and silver. The system was heated in NVT (step 1)and then equilibrated in NVT (step 2). However when I restart the simulation at beginning of step 2, pressure blows up and atoms are lost shortly. Here are what I found about this problem:

  1. When I remove the silver atoms, there is no problem. So the problem lies in water and silver interaction. But there is no problem during step 1 where water and silver interaction also exists.

  2. I tried with “read_restart ** ** remap” option, it doesn’t help, as it applies to atoms out of domain immediately after restart.

  3. I tried the “fix nve/limit”, it doesn’t help. Pressure is still high and atoms or bonds still go lost.

  4. I am using “boundary p p p”.

Attached are my inputs. Waterbox is created with VMD-SOLVATE.
Is there anything wrong with my setting or any bug in the restart binary file?

Thanks a lot for help.

in.test (2.56 KB)

rest.test (2.23 KB)

Dear Everyone,

I have a system with water and silver. The system was heated in NVT (step
1)and then equilibrated in NVT (step 2). However when I restart the
simulation at beginning of step 2, pressure blows up and atoms are lost
shortly. Here are what I found about this problem:
1. When I remove the silver atoms, there is no problem. So the problem
lies in water and silver interaction. But there is no problem during step 1
where water and silver interaction also exists.

2. I tried with "read_restart ** ** remap" option, it doesn't help, as it
applies to atoms out of domain immediately after restart.

3. I tried the "fix nve/limit", it doesn't help. Pressure is still high
and atoms or bonds still go lost.

4. I am using "boundary p p p".

Attached are my inputs. Waterbox is created with VMD-SOLVATE.
Is there anything wrong with my setting or any bug in the restart binary
file?

​it is much more likely, that the issue is with your general system ​setup.

​have you tried running the water system by itself?
have you tried continuing your simulation without restarting?
are you sure you want to simulate a flexible water model? at a timestep of
1fs??

axel.

I have tried with water box only. No problem.
I tried running without restart (all in one script file). No problem.

I tried on Mac, it behaved the same.
I tried with single core and multi-cores, same results.

I am using Ubuntu 14.04 and mac Sierra.

Does it sound like a problem with restart binary file?

I have tried with water box only. No problem.
I tried running without restart (all in one script file). No problem.

I tried on Mac, it behaved the same.
I tried with single core and multi-cores, same results.

I am using Ubuntu 14.04 and mac Sierra.

Does it sound like a problem with restart binary file?

​please produce a complete and minimal test system (as few atoms as
possible and with as few lines of input as possible) that can clearly
reproduce what you are observing. your descriptions are too vague and there
are far too many doubts and needles complication in the inputs you
produced. i don't believe that you can run a stable simulation of flexible
water at room temperature with a 1fs timestep. i've done quite a few
simulations of rigid and flexible water in my time...

axel.

Hi Axel,

I have cleaned up the input files as attached.
I was using time step at 1fs following other papers. The water model is FPC. What would you suggest as a start?

My test results are:

  1. Pure water: no problem.
  2. water and silver without restart: no problem.
  3. water and silver with restart: pressure and energy blow up

System:
Ubuntu 14.04 in VM Player
LAMMPS: 2016 Nov and 2017 Mar.

Thanks for help.
Best,

in.test (1.86 KB)

rest.test (891 Bytes)

watercube.test (49.6 KB)

Hi Axel,

I have cleaned up the input files as attached.
I was using time step at 1fs following other papers. The water model is
FPC. What would you suggest as a start?

​the issue you are having is not caused by restarting directly, but rather
is in your create_box command. you are not reserving ​any storage for
exclusions/specials, which seems to be corrupting your topology data, but
in a way, in your special case, that only seems to become visible when you
are writing a restart or a data file.

as for the water, you need to check carefully. your input is at extremely
low temperature. there you may get away with 1fs. i'd be very surprised if
that will still work well without shake for room temperature.

​axel​

Axel,

Thanks you for valuable advice.
I am trying the special bonds and fix shake. But there is an error when I use “fix shake”:

ERROR: Did not find fix shake partner info (…/fix_shake.cpp:884)
Last command: fix fpcShake water shake 0.0001 20 100 t 1 2 b 1

I used “fix shake” AFTER “delete_atoms bond yes”.
For molecular systems like water, “compress yes” option doesn’t work.
I am not sure what else I should do.

Thanks!

Axel,

Thanks you for valuable advice.
I am trying the special bonds and fix shake. But there is an error when I
use "fix shake":

​did you read up on the water potential? is it a rigid model? does it
*require* SHAKE?

ERROR: Did not find fix shake partner info (../fix_shake.cpp:884)
Last command: fix fpcShake water shake 0.0001 20 100 t 1 2 b 1

I used "fix shake" AFTER "delete_atoms bond yes".
For molecular systems like water, "compress yes" option doesn't work.
I am not sure what else I should do.

start by not doing multiple changes at the same time.​ you are just making
your like needlessly hard.

BTW: you can also simplify your life, by creating the box through using
read_data on your water box data file (and just add one atom type for the
silver), as the first read_data command can automatically reserve type and
related info storage when it asks to create the box.

axel.

I double-checked my water model description. I shouldn’t use SHAKE as it is flexible model.
Now I have solved the problem with adding:
special_bonds lj/coul 1.0 1.0 0.0

I am still not quite sure if this is correct as few, if any, papers talked about this weighting factor.
Is there any general rule for this weighting factor, or is it just try and error?

Thanks a lot for help.

I double-checked my water model description. I shouldn't use SHAKE as it
is flexible model.
Now I have solved the problem with adding:
special_bonds lj/coul 1.0 1.0 0.0

​no. you just replaced one problem with a different one: you invalidated
your water force field by turning off exclusions.​ the proper correction
would be to augment your create_box command to reserve suitable storage for
flagging special atoms, as i have already told you.

you have now repeatedly violated the most important rule of using computer
software: do not use any commands or flags unless you know exactly what
they are doing and why.

I am still not quite sure if this is correct as few, if any, papers talked
about this weighting factor.
Is there any general rule for this weighting factor, or is it just try and
error?

​the fact, that you are asking this question makes me very sad, as it shows
how very little you know about what you are doing. of course it is
documented and part of your force field.​ for most water models you need to
have a lj/coul 0.0 0.0 0.0 settings (well, the last number is for 1-4
pairs, but since water molecules have no 1-4 pairs, it doesn't make a
difference what you set the 1-4 scaling to).

i am getting tired of having to spoon feed you all kinds of details about
MD simulations, that your adviser and your colleagues should have taught
you _before_ you even consider starting simulations with LAMMPS.

axel.

I am so sorry for bothering too much.
I did believe special_bonds should be “0.0 0.0 0.0” because the intramolecular potential consists of only the harmonic bond stretching and bond angle vibration terms.
But when I followed your advice to change “create_box” command, the blowing-up problem persists.
Then I guess I may play with setting 1-2 and 1-3 pair weighting factor to 1, “A value of 1.0 means include the full interaction”.
I’ll look again to find out more details. Thank you very much for help.
Best regards,