【questions】about group.cpp

Dear All,

Recently, I am trying to modify fix_addforce.cpp to get rid of the limit imposed by group, basically, I want to be able to add different forces to each individual molecules based on their position and direction. So in the new addforce(name it fix_push), I will use coordinate information of some atoms belong to each molecule to decide the direct of my addforce.

the main part of code I modify is in post_force:

"

double unwrap[3];
double unwrap2[3];
double dx,dy,dz, dist;
int j;
for(int i = 0; i < nlocal; i++){ // i is the local index, not the id
if(mask[i] & groupbit){ // if the particle is in the right group
if(region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
if( i-mol_length*int(i/mol_length) == 0){ // this is the first particle of a molecule
domain->unmap(x[i],image[i],unwrap);
j = i+1;
domain->unmap(x[j],image[j],unwrap2);

// compute the versor
dx = unwrap[0]-unwrap2[0];
dy = unwrap[1]-unwrap2[1];
dz = unwrap[2]-unwrap2[2];
dist = sqrt(dxdx + dydy + dzdz);
dx /= dist;
dy /= dist;
dz /= dist;
// now we need to find the position of particles i and j, and bla bla
foriginal[0] -= xvalue
unwrap[0] + yvalueunwrap[1] + zvalueunwrap[2]; // energy change
foriginal[1] += fvaluedx;
foriginal[2] += fvalue
dy;
foriginal[3] += fvaluedz;
f[i][0] += fvalue
dx;
f[i][1] += fvaluedy;
f[i][2] += fvalue
dz;
}
}
}
"

It went on good for most of the case, but sometimes the molecule is not moving as I hoped. Sometimes, it is moving forward and move backward and is oscillating all the time. So I’m thinking is it because I didn’t consider what the difference of i of nlocal and id of atoms. I’m not very clear about how lammps communication and organized in multi-processors(parallel).

Will the id of atoms change during the process of running ? and is it ok that I manipulate as i show in the code? and anyone can tell me the use of foriginal here and how to understand “domain-> unmap(x[j],image[j],unwarp)”, I know the function of it, but I don’t know for example, what is image[i].

anyone who can tell me any of my questions are welcome to reply. I may ask too many questions at one time.

Dear All,

Recently, I am trying to modify fix_addforce.cpp to get rid of the limit
imposed by group, basically, I want to be able to add different forces to
each individual molecules based on their position and direction. So in the
new addforce(name it fix_push), I will use coordinate information of some
atoms belong to each molecule to decide the direct of my addforce.

[...]

It went on good for most of the case, but sometimes the molecule is not
moving as I hoped. Sometimes, it is moving forward and move backward and is
oscillating all the time. So I'm thinking is it because I didn't consider
what the difference of i of nlocal and id of atoms. I'm not very clear about
how lammps communication and organized in multi-processors(parallel).

when running in parallel, atoms are distributed across MPI ranks (i.e.
processors) according to their position into subdomains (-> domain
decomposition).
those atoms are referred to in LAMMPS as "local" atoms and have the
index i from 0 to atom->nlocal-1. in addition, LAMMPS maintains copies
of atoms from neighboring subdomains, which are referred to as "ghost"
atoms. those have the index i from atom->nlocal to atom->nlocal +
atom->nghost - 1.
the atom ID of individual atoms is stored in atom->tag[]. since atoms
move between subdomains, the identity of atoms with the same local
index between time steps is not preserved. you can find out the atom
ID of an atom via atom->tag[i] and you can get the local index for a
known atom ID via atom->map(id).

please note, that a Fix::post_force() method may only add forces to
"local" atoms, so you will likely need to do a so-called reverse
communication to send force data assigned to ghost atoms back to the
MPI rank that owns them.

Will the id of atoms change during the process of running ? and is it ok
that I manipulate as i show in the code? and anyone can tell me the use of
foriginal here and how to understand "domain-> unmap(x[j],image[j],unwarp)",
I know the function of it, but I don't know for example, what is image[i].

atom->image[] contains the image flags, i.e. the count how many times
an atom has passed through the periodic boundary in x, y, or z. this
is normally encode in 3x10bits of a single 32-bit integer (or 3x
22bits of a 64-bit integer when LAMMPS is compiled with
-DLAMMPS_BIGBIG)

axel.

Dear Axel,

Thank you so much for your reply! you solve me a lot of doubts. And one more questions, can I run for loop in fix_addforce like “for(int id; id <= ntotal; id ++)” ?

Besides, I just add a force not from other atoms, which means there is no opposite force on other atoms, in this case maybe I don’t need send force data assigned to ghost atoms? Because I am thinking , once the atom is in local, an additional force will add to the atoms.

Best Regards,
Pin

Dear Axel,

Thank you so much for your reply! you solve me a lot of doubts. And one more
questions, can I run for loop in fix_addforce like "for(int id; id <=
ntotal; id ++)" ?

hard to say. depends on what you do inside that loop. but remember,
you have only a subset of the total system on each MPI rank. you
really should read the original LAMMPS paper and some of steve's
workshop presentations about the parallelization and code structure in
LAMMPS. please also keep in mind, that atom ids are typically not
guaranteed to be consecutive and starting at 1.

Besides, I just add a force not from other atoms, which means there is no
opposite force on other atoms, in this case maybe I don't need send force
data assigned to ghost atoms? Because I am thinking , once the atom is in
local, an additional force will add to the atoms.

if you want to do the same thing to all atoms within the same
molecule, there is a good chance, that some of those atoms are ghost
atoms. without knowing more details of what and how you want to do
things, it is impossible to say for certain. keep in mind, that is
said you "may" need reverse communication.

it looks to me, that you need to think some more about what working on
*distributed* data in parallel means.

axel.

Just want to be sure you realize the fix addforce command
has an option to use an atom-style variable for the added
force, which can be a function of each atom’s position
and velocity.

Steve

Hi Steve,

I do realize there is a way to use atom-style variable, but I really don’t understand how to use it. And I think it cannot give a velocity depends on the atom and it’s near by atoms…(I hope to give a velocity depends on the molecule direction) if it is possible, could you please let me know?

Thanks&Regards,

Pin

If you want to add a force to each atom that depends

on attributes of its nearby atoms, or the molecule it

belongs to, then no, you can’t use an atom-style

variable to do that. You’d have to write code

to do that, as in code up a fix of your own.

Steve