Angles and ghost particles

LAMMPS LAMMPS Mailing List Mirror
1

Dear lammps users,

I’m computing total volume of triangulated mesh in my own angle (the
code of the routine is at the end of the message). Mesh is stored in
lammps using angles. For instance, I have cube made of 12 triangles.
Cube is 4x4 with the center in origin. Simulation domain is cube 6x6.
Cube vertices are frozen particles.
When I compute cube's volume on one processor, it works properly. But
when on 2 - I have weird triangles:
on the second processor I have a triangle constructed from one “real”
vertex and 2 “ghosts” with the following coordinates:
-2 -2 -4
2 - 2 -2
2 -2 -4

This triangle is not even similar to the original one which (to my mind) is
-2 -2 2
2 -2 -2
2 -2 2

Thus it’s area is different from what I would expect. Do you have an
idea about handling such situations? Do I need to get
domain->minimum_image from all triangular's points participating in
calculations(is this the only way)?

Code of the routine:

  double computeTotalVolume(double** x, int** anglelist, int
nanglelist, Domain*& domain, MPI_Comm comm)
  {
    double volume = 0.0;

    for (int n = 0; n < nanglelist; n++) {
      int i1 = anglelist[n][0];
      int i2 = anglelist[n][1];
      int i3 = anglelist[n][2];

      // 1st bond - 21
      double x21[3];
      MathExtra::sub3(x[i2], x[i1], x21);

      // 3nd bond - 31
      double x31[3];
      MathExtra::sub3(x[i3], x[i1], x31);

      // normal to the triangle
      double normal[3];
      MathExtra::cross3(x21, x31, normal);

      double massCenter[3];
      computeMassCenterTriangle(x, i1, i2, i3, massCenter);
      domain->minimum_image(massCenter);

      double mcOnNorm = MathExtra::dot3(massCenter, normal);

      assert(mcOnNorm > 0.0); //if it is < 0 then mesh is not oriented properly
      volume += mcOnNorm / 6.;
    }

    double totalVolume;
    MPI_Allreduce(&volume, &totalVolume, 1, MPI_DOUBLE, MPI_SUM, comm);

    assert(totalVolume > 0.); //check that mesh is positively oriented
    return totalVolume;
  }

the mesh:

Atoms

1 1 1 -2.0 2.0 -2.0
2 1 1 2.0 2.0 -2.0
3 1 1 2.0 -2.0 -2.0
4 1 1 -2.0 -2.0 -2.0
5 1 1 -2.0 2.0 2.0
6 1 1 2.0 2.0 2.0
7 1 1 2.0 -2.0 2.0
8 1 1 -2.0 -2.0 2.0

Angles

1 1 1 2 3
2 1 3 4 1
3 1 5 7 6
4 1 7 5 8
5 1 1 6 2.
6 1 6 1 5
7 1 2 7 3
8 1 7 2 6
9 1 8 3 7
10 1 3 8 4
11 1 4 8 5
12 1 4 5 1

2

Unless the coords of your frozen particles need to
be unwrapped across periodic boundaires, I don't
know why you would need minimum images.

LAMMPS will add periodic offsets to ghost particles
that are across PBCs. So that each processor
owns ghost particles that are "near" the particles
it owns.

Steve

3

Because the original angle:
-2 -2 2
2 -2 -2
2 -2 2
Is not similar to the angle which belongs to one of the processes:
-2 -2 -4 (ghost)
2 - 2 -2 (real)
2 -2 -4 (ghost)
That's why the area of this triangle is different from the area of the
original one. The center of mass is also different(and cannot be
mapped to the original domain by adding vector (0,0, 6)).
I'm not sure what I shall use - minimal_image or something else
(something like Domain::unmap, for instance). I just don't know
another way to compute real volume of the particular tetrahedron(where
the triangle is basement).

4

it is difficult to discuss items like this without
being able to reproduce it. there is a lot of
(crucial) information missing, would you mind
posting this as well?

thanks,
axel.

5

I have created a small fix which prints computed volume. The input
script is in.cube2, data file is cube.atom.
If it is started in one process, it will print "Volume: 64"
But in case if it is run in np > 1, there will be assertion fail
because volume of tetrahedron will be less than 0.0.
I attached also fix_freeze because I modified it a little bit.

fix_printVolume.cpp (2.66 KB)

fix_printVolume.h (444 Bytes)

cube.atom (860 Bytes)

in.cube (1.01 KB)

fix_freeze.cpp (3.68 KB)

fix_freeze.h (1.07 KB)

6

kirill,

please find attached a version of your fix that
corrects for the inconsistencies that you had
in your implementations. in particular.

  • the middle atom is the one with the two bonds
    not the first.

  • you need to apply minimum image conventions
    not on positions, but on position differences

  • you should use those distances to compute
    the minimum image corrected center of mass.

this code gives consistent results for any number
of processors uses.

HTH,
axel.

fix_printVolume.cpp (2.71 KB)

7

Thank you for the detailed reply.
There is a mistake in the code you sent me. And it is not clear how to fix it.
If I apply minimum_image to the differences of positions I will get
wrong area of a triangle.

For instance,
x12 = (4, 0 , 0) -> (-2, 0, 0)
x32 = (0, 4, 0) -> (0, -2, 0)
Where '->' means 'apply minimum image'

Thus the normal=[a32, a12] is (0, 0, -4) and it's norm is 4, while the
real normal (without applying minimum_image) would be (0, 0, -16).
So because of this error in area computation, this code returns 16 for
the volume of 4x4x4 cube instead of 64.

8

Thank you for the detailed reply.
There is a mistake in the code you sent me. And it is not clear how to fix it.

Use a larger box. Your periodic images are closer than your original coordinates. Under consideration of PBC with minimum image conventions your volume is that small.

Axel

9

i should add a few more comments, before you run into the next
problem, or people reading this thread get confused.

the test system has "bonds" that are 4 angstrom long,
but the closest image is only 2 angstrom away. since
the system was indicated to be periodic, volume of 16
cubic angstrom is correct.

doubling the size of the simulation box rights this and
the result is the expected one of 64 cubic angstrom.

furthermore, since the bond length is much larger
than the cutoff plus neighbor list skin, one needs to
add the command:

communicate cutoff 4.0

so that atoms belonging to the same
(tri)angle are not "lost" in a parallel run.

cheers,
    axel.