Segmentation fault with TIP5P model

Dear LAMMPS users and developers,

I have been trying to create a TIP5P water model in LAMMPS(15 May 2015) using "fix RIGID all rigid/nve/small molecule.", since SHAKE and RATTLE are saying "ERROR: Shake cluster of more than 4 atoms".

It seems like to me that the temperature calculation for different bins in the simulation is non-trivial because the number of degree-of-freedoms are not same for different type of atoms (O, H, and lone pairs). Therefore, it may not be very possible to define it using "compute temp/chunk with adof x" for chunk calculations. ("compute temp" for group of bins also estimated very inconsistent temperature results.)

As a result, I decided on trying this using ke/rigid and erotate/rigid commands where I can calculate temperature specifying all degree-of-freedoms manually. Everything works fine except the case with "compute rke all erotate/rigid RIGID", which leads to segmentation fault. I looked at the previous mail list comments about this and I read that this is not a structural LAMMPS error. So is there anyone with some advice if it is a bug or how to overcome (I can provide an inout file)? I will also appreciate any LAMMPS trick about the temperature profile of rigid TIP5P model.

Thank you to all.

Dear LAMMPS users and developers,

I have been trying to create a TIP5P water model in LAMMPS(15 May 2015)
using "fix RIGID all rigid/nve/small molecule.", since SHAKE and RATTLE
are saying "ERROR: Shake cluster of more than 4 atoms".

this is a documented feature of the SHAKE and RATTLE implementation in LAMMPS.
this restriction increases parallel efficiency significantly.

It seems like to me that the temperature calculation for different bins
in the simulation is non-trivial because the number of
degree-of-freedoms are not same for different type of atoms (O, H, and
lone pairs).

that reasoning makes no sense.

Therefore, it may not be very possible to define it using
"compute temp/chunk with adof x" for chunk calculations. ("compute temp"
for group of bins also estimated very inconsistent temperature results.)

i disagree. the documentation even explains how to determine "x". what
value of "x" have you been using for "adof x" in compute "temp/chunk"?

what do you mean by "inconsistent" temperature results?

As a result, I decided on trying this using ke/rigid and erotate/rigid
commands where I can calculate temperature specifying all
degree-of-freedoms manually. Everything works fine except the case with
"compute rke all erotate/rigid RIGID", which leads to segmentation
fault. I looked at the previous mail list comments about this and I read
that this is not a structural LAMMPS error. So is there anyone with some
advice if it is a bug or how to overcome (I can provide an inout file)?

the first advice for this kind of issue is always the same: if you
don't use the latest patchlevel of LAMMPS, update to the very latest
developer version and see, if the issue persists.
i cannot reproduce it, which suggests, that the problem has been resolved.

I will also appreciate any LAMMPS trick about the temperature profile of
rigid TIP5P model.

there is no special "trick" needed. just use fix ave/chunk and compute
chunk/atom bin/1d

axel.

Dear LAMMPS users and developers,

I have been trying to create a TIP5P water model in LAMMPS(15 May 2015)
using "fix RIGID all rigid/nve/small molecule.", since SHAKE and RATTLE
are saying "ERROR: Shake cluster of more than 4 atoms".

this is a documented feature of the SHAKE and RATTLE implementation in
LAMMPS.
this restriction increases parallel efficiency significantly.

It seems like to me that the temperature calculation for different bins
in the simulation is non-trivial because the number of
degree-of-freedoms are not same for different type of atoms (O, H, and
lone pairs).

that reasoning makes no sense.

Therefore, it may not be very possible to define it using
"compute temp/chunk with adof x" for chunk calculations. ("compute temp"
for group of bins also estimated very inconsistent temperature results.)

i disagree. the documentation even explains how to determine "x". what
value of "x" have you been using for "adof x" in compute "temp/chunk"?

what do you mean by "inconsistent" temperature results?

######################################################################################################

compute cc1 all chunk/atom bin/1d x lower 0.2 units reduced
compute myChunk all temp/chunk cc1 temp
fix 1 all ave/time 50 10 500 c_myChunk file tmp.out mode vector

######################################################################################################

When I use chunks to calculate temperature with the commands above, this is
the resulting temperature profile:

# Time-averaged data for fix 1
# TimeStep Number-of-rows
# Row c_myChunk[1]
500 5
1 117.931
2 126.961
3 123.127
4 116.555
5 118.099

But in the same time, this is the resulting last 10 timesteps of
"thermo_style custom step temp "

     400 299.82212
     410 302.68525
     420 301.02635
     430 309.38812
     440 301.86262
     450 302.93716
     460 313.96966
     470 298.85252
     480 294.06505
     490 308.0981
     500 311.66066

These are obtained by "rigid/nve/small" after equilibration with
"rigid/npt/small" targeted at 300 K and 1 atm. So as you see, temperatures
of the chunks are not similar with temperature results of thermo_style temp
(which is the equivalent of compute temp with group all, as far as I know).
This was what I meant by inconsistent temperatures (sorry about bad
explanation in previous). I think this is happening because I did not

indeed, this output is consistent with the documentation and that is
what matters.

specify adof in temp/chunk, but adof is different for oxygens and
hydrogens+lone pairs for TIP5P model, so how can I treat them differently?

this is nonsense. why should the degrees of freedom by different?
please explain.

from the perspective of statistical mechanics a particle is a
particle. it doesn't make a difference what it stands for.

This is why I thought that "compute temp/chunk" may not be useful in this
particular problem. I will be very pleased, if you can tell me why this does
not make sense.

be pleased then.

axel.

I’m not following all the details here, but if

you are using compute temp with a fix rigid (for a rigid

TIP5P model), then the temperature compute will

remove DOFs for each molecule to account for

the rigid constraints.

If you are also using

compute temp/chunk with spatial bins (or whatever)

as the chunk, the the doc page says this:

Note that currently the global and per-chunk temperatures calculated by this compute only include translational degrees of freedom for each atom. No rotational degrees of freedom are included for finite-size particles. Also no degrees of freedom are subtracted for any velocity bias or constraints that are applied, such as compute temp/partial, or fix shake or fix rigid. This is because those degrees of freedom (e.g. a constrained bond) could apply to sets of atoms that are both included and excluded from a specific chunk, and hence the concept is somewhat ill-defined. In some cases, you can use the adof and cdof keywords to adjust the calculated degress of freedom appropriately, as explained below.

So you would get more DOF with compute temp/chunk

and hence a lower temperature.

Steve

Thank you so much for this detailed explanation Dr. Plimpton. An alternative solution to that ill-defined concept is discussed in the paper below, fyi. Maybe I can try to modify compute_temp_chunk files accordingly.

Zhang, M., Lussetti, E., de Souza, L. E., & Müller-Plathe, F. (2005). Thermal conductivities of molecular liquids by reverse nonequilibrium molecular dynamics. The Journal of Physical Chemistry B, 109(31), 15060-15067.

Dr. Kohlmeyer, you can be sure that I am pleased, since you made me think about the definition of DOF for atoms with constraints again. Please let me read and learn more about this, so that I will try explaining, as you told so.