Accessing periodic image of atoms in my own fix

Just checked it and it’s the same problem.
Apprently, for the 0th timestep the unmapped coordinates of the edge/corner-atoms are incorrect (fix postforce)

I’ll try to delay my fix by 1 timestep.
If you have some hints, feel free to response.

best regards,
frank.

Thanks for the advice.
In fact I’ve just taken a look into the fix_spring_self.cpp file and it seems like domain->unmap(…) is just what I’m looking for. Im gonna test it right now.

Thank you and
best regards,
frank.

I’ve written a fix for pressure-absorbing boundary conditions for fcc-lattices which I’m using for my shockwave simulations.

It uses a monochromatic wave-solution approach for integrating the equations of motion for the atoms which represent the boundary group.

One input parameter is the distance between the current atomic positions of the boundary atoms and the equilibrium position of these atoms at the very first timestep.

Actually it works surprisingly well.

However, if the simulations runs for a while, it might happen that some of my the atoms travel across the periodic boundaries and a simple distance-calculation according to “pos_new- pos_old” is obiously wrong. I have to take into account the periodic images.

so this is what I’m doing (pseudocode):

loop “i” over all atoms:
if atom is part of the boundary-group:
imageint *imgflag = atom->image;

flag_x=static_cast((( imgflag[i] & IMGMASK)-IMGMAX);

flag_y=static_cast(( imgflag[i] >> IMGBITS & IMGMASK) - IMGMAX);
flag_z=static_cast(( imgflag[i] >> IMG2BITS) - IMGMAX);

//get the “true” positions

x_pos=x[i][0]+flag_x*Ly;

y_pos=x[i][1]+flag_yLy;
z_pos=x[i][2]+flag_y
Ly;

this chunk seems to have a cut and paste error.
why don’t you save yourself some grief and use one of the domain->unmap() methods instead?

​axel.

Just checked it and it's the same problem.
Apprently, for the 0th timestep the unmapped coordinates of the
edge/corner-atoms are incorrect (fix postforce)

I'll try to delay my fix by 1 timestep.
If you have some hints, feel free to response.

​at this point, i would need a simple/small test case, so i can review what
is happening with a number of debugging tools and determine whether this is
a general issue in LAMMPS or a particular problem of your input deck or
code.

i assume you are working off an up-to-date git or subversion checkout of
LAMMPS, right?

axel.

Dear lammps users,

Sorry for reopening this old thread. We are facing a similar problem. (16Feb 2016 version of lammps)

We are simulating biological cell suspensions which are meshed in a triangular fashion (PFA).

The cell dynamics depend on its volume, hence an additional potential function is required. We added a custom class angle_custom.cpp for the same.

The volume of the cell is calculated using the signed tetrahedron value and incorporated in angle_custom.cpp file,

the volume of the tetrahedron with vertices (0,0,0), (a1,a2,a3), (b1,b2,b3) and (c1,c2,c3) is

[a1b2c3 +a2b3c1 + a3b1c2 -a1b3c2 -a2b1c3 - a3b2c1]/6.

Code segment for calculating volume (part of the added custom class):

>>>

double **x = atom->x;
double **f = atom->f;
int **anglelist = neighbor->anglelist;
int nanglelist = neighbor->nanglelist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
imageint *image = atom->image;

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


double ui1[3], ui2[3], ui3[3];

// Setting a flag to print the selected atoms and timesteps

int z = 0;
if ((atom->tag[i1] == 171) && (atom->tag[i2] == 165) && (atom->tag[i3] == 69) && (update->ntimestep >= 3820)) {
char str3[512];
z = 1;
}

domain->unmap(x[i1],image[i1],ui1, &z);
domain->unmap(x[i2],image[i2],ui2, &z);
domain->unmap(x[i3],image[i3],ui3, &z);

vol += (ui1[0]*ui2[1]*ui3[2] + ui1[1]*ui2[2]*ui3[0] + ui1[2]*ui2[0]*ui3[1] - ui1[0]*ui2[2]*ui3[1] - ui1[1]*ui2[0]*ui3[2] - ui1[2]*ui2[1]*ui3[0])/6;

}

We used unmap function from domain.cpp to obtain unwrapped coordinates, but simulation is running fine until atoms start crossing the boundary.

As soon as the value returned by *image (>>>imageint *image = atom->image) goes to “0”, the following values in domain.cpp changes as follows resulting in an error in the unwrapped coordinates. xbox=-512, ybox=-512, zbox=-512

Modified unmap function in domain.cpp for reference:
>>> void Domain::unmap(double *x, imageint image, double *y, int *z)
{

if (*z == 1) {

char str3[512];
sprintf(str3,“inside unmap initial: “” %g” " %g" " %g" " %i",
x[0],x[1],x[2],image);
error->warning(FLERR,str3,0);
}

int xbox = (image & IMGMASK) - IMGMAX;
int ybox = (image >> IMGBITS & IMGMASK) - IMGMAX;
int zbox = (image >> IMG2BITS) - IMGMAX;

if (triclinic == 0) {
y[0] = x[0] + xbox*xprd;
y[1] = x[1] + ybox*yprd;
y[2] = x[2] + zbox*zprd;
} else {
y[0] = x[0] + h[0]*xbox + h[5]*ybox + h[4]*zbox;
y[1] = x[1] + h[1]*ybox + h[3]*zbox;
y[2] = x[2] + h[2]*zbox;
}

if (*z == 1) {

char str4[512];
sprintf(str4,“inside unmap after: “” %g” " %g" " %g" " %g" " %g" " %g" " %i",
x[0],y[0],x[1],y[1],x[2],y[2],image);
error->warning(FLERR,str4,0);
}
}

Output of the same for reference:

WARNING: inside unmap initial: 6.03009 0.293295 0.0675477 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 6.03009 6.03009 0.293295 0.293295 0.0675477 0.0675477 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: 6.03893 -0.332659 -0.0224239 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 6.03893 6.03893 -0.332659 -0.332659 -0.0224239 -0.0224239 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: 5.95966 0.00464437 -0.546244 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 5.95966 5.95966 0.00464437 0.00464437 -0.546244 -0.546244 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: 6.03122 0.294202 0.0678906 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 6.03122 6.03122 0.294202 0.294202 0.0678906 0.0678906 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: 6.03935 -0.333786 -0.0221702 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 6.03935 6.03935 -0.333786 -0.333786 -0.0221702 -0.0221702 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: 5.9602 0.00525877 -0.547087 537395712 (…/domain.cpp:1443)
WARNING: inside unmap after: 5.9602 5.9602 0.00525877 0.00525877 -0.547087 -0.547087 537395712 (…/domain.cpp:1466)
WARNING: inside unmap initial: -5.96886 0.293208 0.0671074 537395713 (…/domain.cpp:1443)
WARNING: inside unmap after: -5.96886 6.03114 0.293208 0.293208 0.0671074 0.0671074 537395713 (…/domain.cpp:1466)
WARNING: inside unmap initial: -5.959 -0.333037 -0.0222073 537395713 (…/domain.cpp:1443)
WARNING: inside unmap after: -5.959 6.041 -0.333037 -0.333037 -0.0222073 -0.0222073 537395713 (…/domain.cpp:1466)
WARNING: inside unmap initial: -6.03818 0.00535828 -0.546978 0 (…/domain.cpp:1443)
WARNING: inside unmap after: -6.03818 -6150.04 0.00535828 -6143.99 -0.546978 -3584.55 0 (…/domain.cpp:1466)

The box size is:
-6.0000000000000000e+00 6.0000000000000000e+00 xlo xhi
-6.0000000000000000e+00 6.0000000000000000e+00 ylo yhi
-3.5000000000000000e+00 3.5000000000000000e+00 zlo zhi

Kindly update if this is a bug in the version of lammps. If so suggest the possible fix for the same.

Thanks and Regards,

Abhijith A

PhD Scholar,

Applied Mechanics Department,

Indian Institute of Technology, Madras.

cell_structure.png

This is more likely a bug in your code. It is impossible to say from just the segment you provide, but the identity of atoms pointed to by those indices i1 i2 i3 can change any time when the neighbor list is rebuild (you would have to do the reverse lookup by using global indices and then atom->map()), and you also need to make certain, that your communication cutoff is sufficient to hold the correct ghost atoms (and not be pointed to a different periodic image with atom->map() or to -1 in case the atom is not present at all).

on top of that, none of the developers will have a serious look at things unless you can not only prove that this is not an issue with your code, but also that the issue is not already resolved with a newer version of LAMMPS.

axel.