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.