[lammps-users] bug with free boundary condition 'fm' for LJ solid

Dear colleagues,

I have the following problem. The model is a Lennard-Jones solid film on a
frozen solid substrate. Below there is the correspondiing Lammps script. This
script runs normally with the 17Jan05 Lammps version, but with 12Arp06 version
it crashes with the 'lost atoms' error after the first several hundred steps.

I should mention that initial velocity distributions generated by the
'velocity' command differ for two Lammps versions. Hence the atoms
trajectories differ as well. However if I arbitrarily change the seed number
in the 'velocity' command the results remain always the same: the 12Apr06
version crashes, the 17Jan05 version runs normally..

I made the data-file from the restart-file saved after the first 100 steps in
order to compare the runs from this identical initial condition. The result is
again the same: the 12Apr06 version crashes, the 17Jan05 version runs normally.

Presumably this problem may be connected with the free ('fm') boundary conditions.

There is no such problem while using eam potential for the same task (the same
system geometry).

Could you give me some advice where the bug can be? Can it be the Lammps bug?

Thank you.

Vladimir

units lj
boundary p p fm

atom_style atomic
lattice fcc 1.03383

region box block 0.0 5.0 0.0 5.0 0.0 14.0
region freezedreg block 0.0 5.0 0.0 5.0 0.0 2.01
region heatedreg block 0.0 5.0 0.0 5.0 2.02 10.0

create_box 1 box

create_atoms 1 box
mass 1 1.0
group freezed region freezedreg
group heated region heatedreg

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

velocity all create 0.95 552233 mom yes rot yes dist gaussian
velocity freezed create 0.0 2745

thermo 50

fix 1 all nve
fix 3 freezed setforce 0.0 0.0 0.0

dump 1 all atom 100 dump_atom

timestep 0.005
run 20000

Vladimir,

Trust the more recent version of LAMMPS. You are
indeed losing atoms. The problem is that you're fixing
the lower z bound, but doing nothing to ensure that no
atoms can move below that z bound. You need to put in
a wall or something to prevent atoms from moving
beyond the lower z bound. Alternatively, you could use
shrink-wrapped (the "s" option) boundary conditions,
which seems to work fine.

There was a problem with earlier versions of LAMMPS
that caused this atom loss to not be detected. That's
why you should trust the more recent version. See:

http://www.cs.sandia.gov/~sjplimp/lammps/bug.10Nov05.html

"22 Mar 2006

Added some logic to the routine that communicates
atoms to new processors, to insure that atoms that
cross non-periodic boundaries are deleted and become
"lost" as they should be. Previously, this would
happen in parallel, but sometimes not in serial (on a
single processor), which could be confusing. Now the
behavior is the same on one or many procs.

This will be part of the next LAMMPS release (by end
of March 06).

Thanks to Jaime Sanchez (U Kentucky) for identifying
the inconsistency."

I hope this helps!

Paul

Paul,

I would like to mention that (as you have recommended) from the beginning I
made the layer of atoms near the lower z bound freezed.

Please could you clarify two further questions connected with this problem:

1. Is the following correct? If atoms in some group has zero velocities then
it is equivivalent i) to use 'fix group setforce 0.0 0.0 0.0' and 'fix all
nve', or ii) simply exclude this atoms from integration 'fix (all-group) nve'.

2. I have found that the model under consideration (see the script below)
works fine in the Lammps 12Apr06 version if I use 'fs' boundary conditions
instead of the 'fm' condition, i.e. 'fs' does not cause 'lost atoms error' but
'fm' does cause such an error under the same other conditions. It looks
strange for me because as far as I understand the only difference of the 'm'
b.c. is the minimum threshold to the maximum z-extent that is determined by
the initial box size.

Am I missing something?

Vladimir

units lj
#boundary p p fs #works fine
boundary p p fm #cause 'lost atoms' error

atom_style atomic
lattice fcc 1.03383

region box block 0.0 5.0 0.0 5.0 0.0 14.0
region freezedreg block 0.0 5.0 0.0 5.0 0.0 2.01

create_box 1 box

create_atoms 1 box
mass 1 1.0
group freezed region freezedreg
group moving subtract all freezed

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

velocity moving create 0.95 552233 mom yes rot yes dist gaussian

thermo 50

fix 1 moving nve

dump 1 all atom 100 dump_atom

timestep 0.005
run 200000

Paul Crozier <[email protected]> said:

Vladimir,

1. Is the following correct? If atoms in some group
has zero velocities then
it is equivivalent i) to use 'fix group setforce 0.0
0.0 0.0' and 'fix all
nve', or ii) simply exclude this atoms from
integration 'fix (all-group) nve'.

Yes, these two approaches should yield identical
results.

2. I have found that the model under consideration
(see the script below)
works fine in the Lammps 12Apr06 version if I use
'fs' boundary conditions
instead of the 'fm' condition, i.e. 'fs' does not
cause 'lost atoms error' but
'fm' does cause such an error under the same other
conditions. It looks
strange for me because as far as I understand the
only difference of the 'm'
b.c. is the minimum threshold to the maximum
z-extent that is determined by
the initial box size.

Looks like you found a bug. LAMMPS slightly expands
shrink-wrapped boundaries to insure atoms are
enclosed, but was apparently doing this only for "s"
and not for "m". In domain.cpp, changing a few "==" to
">=" makes it do the right thing. So instead of:

  if (nonperiodic > 0) {
    if (xperiodic == 0) {
      if (boundary[0][0] == 2) boxxlo -= SMALL;
      if (boundary[0][1] == 2) boxxhi += SMALL;
    }
    if (yperiodic == 0) {
      if (boundary[1][0] == 2) boxylo -= SMALL;
      if (boundary[1][1] == 2) boxyhi += SMALL;
    }
    if (zperiodic == 0) {
      if (boundary[2][0] == 2) boxzlo -= SMALL;
      if (boundary[2][1] == 2) boxzhi += SMALL;
    }
  }

We have:

  if (nonperiodic > 0) {
    if (xperiodic == 0) {
      if (boundary[0][0] >= 2) boxxlo -= SMALL;
      if (boundary[0][1] >= 2) boxxhi += SMALL;
    }
    if (yperiodic == 0) {
      if (boundary[1][0] >= 2) boxylo -= SMALL;
      if (boundary[1][1] >= 2) boxyhi += SMALL;
    }
    if (zperiodic == 0) {
      if (boundary[2][0] >= 2) boxzlo -= SMALL;
      if (boundary[2][1] >= 2) boxzhi += SMALL;
    }
  }

I'll attach the modified domain.cpp file.

Paul

domain.cpp (18.4 KB)