Fix Shake takes forever to run

LAMMPS LAMMPS General Discussion
1

I am running molecular simulations using EPM2 model of CO2. This is the force field data:


     units        real
    atom_style   full
    bond_style   harmonic
    angle_style  harmonic
    pair_style   lj/cut/coul/long  10.0
    kspace_style pppm  1.0e-5
    pair_modify  mix arithmetic

     bond_coeff   1        450   1.149  

    angle_coeff  1        147   180
    # explanation:  http://lammps.sandia.gov/doc/angle_harmonic.html

    pair_coeff   1 1     0.055898  2.757
    pair_coeff   2 2     0.159984  3.033
    group co2 type 1 2
    fix co2RigidShake co2 shake  1e-3 100 0 b 1 a 1

However, my simulation stops at the following step (After minimization) and does not go forward even though I can see that it is still running on the cluster:

Total # of neighbors = 274000
Ave neighs/atom = 91.3333
Ave special neighs/atom = 2
Neighbor list builds = 0
Dangerous builds = 0
3000 atoms in group co2
Finding SHAKE clusters ...
  0 = # of size 2 clusters
  0 = # of size 3 clusters
  0 = # of size 4 clusters
  1000 = # of frozen angles

I assume something must be wrong with the fix shake setup. Any ideas?

2

That assumption can be tested with a short run, substituting either fix rigid or a fully-flexible simulation. Many things that you (or I!) assume to be true, turn out not to be.

1 Like
3

How about checking out the documentation?
https://docs.lammps.org/fix_shake.html#restrictions

1 Like
4

Thanks. Using fix rigid/npt solved the problem, but now pressure oscillates a lot such that the system alternates between liquid and gas continuously. I also used fix rigid/npt/small but observed the same outcomes. I reverted back to fix npt and set bond and angle coefficients to zero but received the error below:

ERROR on proc 187: Bond atoms 547 548 missing on proc 187 at step 222 (src/ntopo_bond_all.cpp:63)
Last command: run ${nNptEql}
5

Is there a question here?

How do you determine that?

6

I use the time averaged density. My question is why this pressure oscilation hapoens with fix rigid and why setting bond and angle coefficients to zero does not work?

7

That is no good parameter. You could just have a dense gas.

Pressure oscillates, especially for small systems (1000 molecules is still fairly small, try replicate 2 2 2 or even bigger). How much the pressure oscillates depends on multiple parameters and settings. First order of business would be to determine if the system conserves energy properly with the given settings without thermostat and barostat.
With a gas, there is very little value to run a barostat, anyway. There have been plenty of discussions about pressure fluctuation and simulation with rigid bodies. So check out the archives.

Does not work how? What did you expect should happen?? and why?
Your choice seems particularly strange to me and it looks to me as if the behavior you observe should be expected when you have bonds defined between atoms without a force pulling them together. But then again, you were only sharing very few bits of information and parts only as vague descriptions, so a proper assessment of the situation is not possible.

1 Like