[lammps-users] Unexpected wall/region harmonic behavior

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am having problems using the wall/region harmonic fix, and closer
inspection reveals that the wall-particle interaction energies are
much larger than expected.

I may be misunderstanding how the wall potential is defined. Say I
have a spherical wall with an r_cutoff of 20 A. Based on the potential
function in the manual ...

E = \epsilon * (r - r_cutoff)^2, for r < r_cutoff

... I expected that the wall–particle energy when an atom reaches
within 20 A of the wall would be zero (since r - r_cutoff = 0), and
increase as the particle gets closer to the wall (due to decreasing r,
and therefore increasing (r - r_cutoff)^2).

However, what I find is that the wall–particle energy (as output by
f_wall-ID) when an atom is at the cutoff distance (i.e. just entering
the wall potential), is not zero, but actually equals ...

E_cutoff = \epsilon * r_cutoff^2

... causing the atom to bounce off violently (often blowing up the system).

Can anyone confirm this? Is this a bug, or — much more likely, I'm
sure — am I missing something?

Thanks!

Note: The E_cutoff values reported here are seen only in the
21-Dec-2010 Windows build of LAMMPS. The 15-Mar-2011 version on a
Linux machine returns "nan" for the energy, but a similar trajectory
to the Windows run. As expected after the 6-Feb-2011 bug fix, the wall
forces reported by the later version are different, but still high for
a particle just entering the wall potential — if this is a bug, I
suspect that the 6-Feb-2011 fix for the forces did not fix it.

2

I'll take a look at this next week.

Steve

3

The only bug I see is that the harmonic() routine
in fix_wall_region.cpp is adding an offset to the energy,
which isn't needed (unlike the other cases). So the
offset variable is unintialized, which could lead to
funny thermo output energy values. But it doesn't affect dynamics.

Two tests: delete the offset term from the harmonic()
routine, leaving eng = epsilon*r*r; and see if that
helps.

Also, you can try the script in.test below, which
works fine for me. Themo output and viz of the system
are both as expected.

Steve

4

Thanks for the quick response, Dr. Plimpton. I tried both your proposed tests.

The script you provided in your reply wouldn't run without first modifying the z-boundary to be periodic, and, when it did run, it prematurely ended with a "Particle on or inside surface of region used in fix wall/region" error. It was not clear to me whether the particles behaved as expected prior to the error.

For the second test, I just ran my own test system using the 22-Mar-11 build of LAMMPS, in which the offset was removed from the harmonic term. In the harmonic() routine, you now have:

  fwall = 2.0*epsilon*r;
  eng = epsilon*r*r;

I am not very adept at debugging C, so it isn't obvious to me how "r" was defined here. If I understand the manual correctly, one must ensure that eng and fwall = 0 when the atom is at a distance r_cutoff from the wall and increase with decreasing distance from the wall -- this is only satisfied if:

   r = (r_cutoff) - (distance from wall)

This is not what I see in my test system, however: at a distance of r_cutoff from the wall, both fwall and eng are non-zero. Specifically, at the point where an atom is r_cutoff away from the wall, the following energies and x, y, z forces are reported:

  3123.6239 -33.566305 -13.287303 247.32414

Note that I was using epsilon = 5 and r_cutoff = 25, and so the reported energy is given by

  eng = epsilon*r_cutoff*r_cutoff

And the magnitude of the force vector is given by

  fwall = 2.0*epsilon*r_cutoff

i.e., in calculating r for the formulae in the harmonic() routine, the distance from the wall is not being subtracted from r_cutoff.

This aberrant force certainly appears to affect the trajectory as well: the particle almost immediate jets off the wall (within 1fs), which is quite impressive for such a heavy particle.

I have attached the config and data files for my test system. Please let me know what, if anything, I'm missing here.

Thanks again.

t4.data (307 Bytes)

t4.in (487 Bytes)

5

you're right - there was a bug. The code was not
offsetting the origin of the harmonic spring in fix wall/region
from the cutoff distance. Try these 3 lines in fix_wall_region::harmonic()

  double dr = cutoff - r;
  fwall = 2.0*epsilon*dr;
  eng = epsilon*dr*dr;

A modified input test script is also below.
I'll post a patch in a moment for this.

Thanks,
Steve

# LJ example with fix wall/region

units lj
atom_style atomic
boundary s s p

lattice sc 0.5

# region for atom creation

region r1 sphere 0 0 0 5
create_box 1 r1
create_atoms 1 region r1

mass 1 1.0

velocity all create 1.44 87287 loop geom

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

neighbor 0.3 bin
neigh_modify delay 0 every 20 check no

fix 1 all nve

# region for bounding wall

region r1a sphere 0 0 0 6
fix 2 all wall/region r1a harmonic 20.0 0.0 1.0

dump 1 all atom 100 tmp.dump
thermo 100
thermo_style custom step temp epair press f_2

run 30000

6

Thanks, Dr. Plimpton. Looks to me like the modifications you described fix the problem.

7

P.S.: I haven't used the other potentials in the wall/region fix (lj12-6, lj9-3, etc), but, glancing over the source file, it looks like you may have to add the cutoff offset to those as well.

No, I think those measure "r" from the wall, not from the spring minimum
which is the cutoff distance.

Steve