I’m writing a pair style based off lj/cut/coul/long and discovered that some ghost atoms aren’t being included in the neighbor lists despite being well within the cutoff.
The setup is: There are 4 atoms fairly close to each other in a line on the x axis. The x direction of the box is also fairly small so ghost atoms should appear on either side of the line of atoms and be included in the neighbor lists. The box is more than twice the size of the cutoff in the y and z directions so it is effectively a one dimensional simulation for simplicity.
I noticed I was getting weird errors in how the forces were being calculated compared to another program so I decided to investigate. I added this line to the j loop of the compute function of the pair style:
if (i==0) fprintf(screen,“x: %f y: %f z: %f q: %f\n”,x[j][0],x[j][1],x[j][2],q[j]);
So, it would tell me the location (and charge) of every neighbor of the very first atom in the simulation. It turned out that the first three neighbors were the other three local atoms. The next four were the images of the four local atoms in the positive x direction. The last one was an image of itself two box lengths away (right on the edge of the cutoff) in the positive x direction. So, a total of eight neighbors. The other atoms all had eight neighbors too, and all only from the positive x direction.
There were no atoms in the negative x direction included in the neighbor list even though there should have been. I exported the location of every ghost atom to check if they were there and the ghost atoms did exist in the negative x direction. They just weren’t included on any neighbor list.
Here is the input file:
"
units real
boundary p p p
atom_style full
timestep 1
lattice sc 1
region box block 0 8 0 40 0 40
create_box 2 box
create_atoms 1 single 0 0 0
create_atoms 2 single 2.4 0 0
create_atoms 1 single 4 0 0
create_atoms 2 single 6.4 0 0
pair_style lj/cut 16
mass 1 1
mass 2 1
pair_coeff * * 1 1
neighbor 0.3 bin
neigh_modify every 1
run 1
"
Run it and it should say that “Ave neighs/atom = 8” when the real number should be around 12ish. Tested on lammps versions 11jan12 and 6jun12.
(I edited down the input file to remove unnecessary commands, if you want I can post the actual input script used but they’re functionally identical and the results/neighbor lists are exactly the same)
Also, here is another version of the input script with only one atom, again, the results and neighbor lists are the same, only ghost atoms in the positive x direction are included in the neighbor list, not the negative x direction, and the neighs/atom is still 8.
"
units real
boundary p p p
atom_style full
timestep 1
lattice sc 1
region box block 0 2 0 40 0 40
create_box 1 box
create_atoms 1 single 0 0 0
pair_style lj/cut 16
mass 1 1
pair_coeff * * 1 1
neighbor 0.3 bin
neigh_modify every 1
run 1
"
-Adam Hogan
Student, University of South Florida