modifying group/group command

Hi All,

I would like to calculate this quantity, F_ij*V_i, F and V are force and velocity, i and j run over the atoms in group 1 and group 2.
Using the compute group/group command in lammps we can calculate F_ij. I have used the compute_group_group.cpp to create a new compute compute_mygroup_mygroup.cpp. My new compute is working correctly for a few thousand steps and then giving wrong values.

The changes I have made in the compute_group_group.cpp file are
After this code

xtmp = x[i][0];

ytmp = x[i][1];

ztmp = x[i][2];

I have added

vx = v[i][0];

vy = v[i][1];

vz = v[i][2];

And then modified this part of the code

if (ij_flag) {

one[1] += delx*fpair;

one[2] += dely*fpair;

one[3] += delz*fpair;

}

if (ji_flag) {

one[1] -= delx*fpair;

one[2] -= dely*fpair;

one[3] -= delz*fpair;

}

to

if (ij_flag) {

one[1] += delxfpairvx;

one[2] += delyfpairvy;

one[3] += delzfpairvz;

}

if (ji_flag) {

one[1] -= delxfpair(-vx);

one[2] -= delyfpair(-vy);

one[3] -= delzfpair(-vz);

}

I have tested it using a simple system with 4 atoms, 2 atoms belongs to group-1 and the other two belongs to group-2.

lattice sc 1.0
region box block 0 21 0 21 0 21 units box
create_box 2 box

create_atoms 1 random 2 1234567 box
create_atoms 2 random 2 7654321 box
mass * 1.0

pair_style lj/cut 2.5
pair_coeff * * 0.0 0.0
pair_coeff 1 2 1.0 1.0 (only atom type 1-2 interact).

group g1 type 1
group g2 type 2

compute cmp1 g1 mygroup/mygroup g2
variable ans1 equal c_cmp1[1]

variable ans2 equal “vx[1]*fx[1]+vx[2]*fx[2]”

thermo_style custom step v_ans1 v_ans2
thermo 100

fix fxnvt all nvt temp 1.0 1.0 1.0

run 20000

Since there are 4 atoms possibly interacting in ans1 (via your

compute) and only 2 atoms interacting in ans2, why

should they be the same?

You could perform a dump of velocities and forces

(see dump local for pairwise forces) and calculate

the answer yourself as a check.

Steve

Hi Steve and All,

I have set
pair_coeff * * 0.0 0.0

pair_coeff 1 2 1.0 1.0

So for ans1- also the interactions are same as ans2 and the result should be identical.

Thanks.

Hi Steve and All,

I have set
pair_coeff * * 0.0 0.0
pair_coeff 1 2 1.0 1.0

So for ans1- also the interactions are same as ans2 and the result should be
identical.

your implementation is wrong. while you can assume F_ij = -F_ji, you
cannot assume that v_i = - v_j but your implementation does.

axel.

Thanks Axel,

By printing velocities and forces of atoms from compute_group_group.cpp I figured out what I’m doing wrong. As you have pointed out the first part of the code (ij_flag) is ok but the velocities in the second part (ji_flag) are wrong.

if (ij_flag) {

one[1] += delxfpairvx;

one[2] += delyfpairvy;

one[3] += delzfpairvz;

}

if (ji_flag) {

one[1] -= delxfpair(-vx);

one[2] -= delyfpair(-vy);

one[3] -= delzfpair(-vz);

}

Now my question is simple, how do I get the actual atom id j so that I can access its velocity in the ji_flag part of the code.

In the below code, i is the actual atom id, but j is some sort of neighbour index but not atom id. When I printed j value, sometimes it is more than the number of atoms (4) in my system.

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];

jlist = firstneigh[i];
jnum = numneigh[i];

for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;

In this part of the code I want to get the actual atom id j so that I can access its velocity.

Many thanks.

Thanks Axel,

By printing velocities and forces of atoms from compute_group_group.cpp I
figured out what I'm doing wrong. As you have pointed out the first part of
the code (ij_flag) is ok but the velocities in the second part (ji_flag) are
wrong.

if (ij_flag) {

one[1] += delx*fpair*vx;

one[2] += dely*fpair*vy;

one[3] += delz*fpair*vz;

}

if (ji_flag) {

one[1] -= delx*fpair*(-vx);

one[2] -= dely*fpair*(-vy);

one[3] -= delz*fpair*(-vz);

}

Now my question is simple, how do I get the actual atom id j so that I can
access its velocity in the ji_flag part of the code.

In the below code, i is the actual atom id, but j is some sort of neighbour
index but not atom id. When I printed j value, sometimes it is more than the

that is nonsense. you are misunderstanding and misreading the LAMMPS
source code at a very fundamental level. both, i and j are "local"
indices into the per-atom storage of the subdomain. there are no
global arrays in LAMMPS indexed by atom id (which is in atom->tag[i]
and atom->tag[j], BTW). it is only by chance and due to the fact that
your test system is extremely small, that "i" index and tag may be
identical. while "i" indexed atom *must* be "owned" atoms, i.e.
between 0 and atom->nlocal, "j" indexed atoms may be "ghost" atoms, if
those are the closest periodic image. whether this is the case, is
determined during neigbor list build.

which reminds me, you will also have to use "comm_modify vel yes" to
make your modified compute work, since velocities for ghost atoms are
usually not communicated (as they are typically not needed) and thus
may be bogus.

i suggest you have a look at the source to pair style dpd, which uses
per-atom velocities (actually the velocity difference) in its pairwise
force computation.

number of atoms (4) in my system.

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];

jlist = firstneigh[i];
jnum = numneigh[i];

    for (jj = 0; jj < jnum; jj++) {
      j = jlist[jj];
      factor_lj = special_lj[sbmask(j)];
      j &= NEIGHMASK;

     In this part of the code I want to get the actual atom id j so that I
can access its velocity.

no, you do *not* want this. it will be wrong and you will be
"rewarded" with segmentation faults for larger systems.

axel.

Thanks Axel for pointing out

  1. comm_modify vel (that velocities are not communicated between procs by default),
  2. pair_dpd.cpp which uses velocities,
  3. both i and j are local.
    My earlier statements are indeed correct only by chance as I have only 4 atoms in my test system and I was using only one proc.

My new compute_mygroup_mygroup.cpp looks like this now

for (ii = 0; ii < inum; ii++) {

i = ilist[ii];
vxi = v[i][0];

for (jj = 0; jj < jnum; jj++) {

j = jlist[jj];
vxj = v[j][0];

if (ij_flag) {

one[1] += delxfpairvxi;
// printf(“velocity: %f\n”, vxi);
}

if (ji_flag) {

one[1] -= delxfpairvxj;
// printf(“velocity: %f\n", vxj);
}

and my lammps test script is below. This time I have 4 atoms in group-1 and 100 atoms in group-2. Again, atoms in group-1

don’t interact with the atoms in the same group.

atom_style full
lattice sc 1.0
region box block 0 10 0 10 0 10 units box
create_box 2 box

create_atoms 1 random 4 1234567 box
create_atoms 2 random 100 7654321 box
mass * 1.0
velocity all create 1.0 1234567

pair_style lj/cut 2.5
pair_coeff * * 0.0 0.0
pair_coeff 1 2 1.0 1.0

group g1 type 1
group g2 type 2

comm_modify vel yes

compute cmp1 g1 mygroup/mygroup g2
variable ans1 equal c_cmp1[1]
variable ans2 equal “vx[1]*fx[1]+vx[2]*fx[2]+vx[3]*fx[3]+vx[4]*fx[4]”

thermo_style custom step v_ans1 v_ans2
thermo 1000

fix fxnvt all nvt temp 1.0 1.0 1.0
run 1000000

From this I expect ans1 and ans2 to be equal, but they are slightly different as shown below

step ans1 ans2

0 -12.705741 -12.705741
1000 5.1573062 5.156231
2000 6.0659603 6.0654934
3000 -11.449794 -11.451709
4000 0.89251503 0.88322259
5000 -2.4597239 -2.4565729
–------------------------------------
996000 5.4019888 5.4017708
997000 0.85114559 0.85150923
998000 3.0144552 3.0176506
999000 3.1692967 3.1715986
1000000 -0.64502429 -0.64450581

To test what’s going on, I have uncommented the two printf statements from the above compute_mygroup_mygroup.cpp

and used only one atom in group-1.

The velocities at a particular time printed out from compute_mygroup_mygroup.cpp are

velocity: 0.313027

velocity: 0.309425
velocity: 0.313027
velocity: 0.313027
velocity: 0.309425

These supposed to be equal as they are velocities of the only one atom in group-1 printed as many times as it has neighbours at any given step. Not sure why/where such difference is coming from …

I’m I still missing something or doing wrong …

Many thanks.

If you have a periodic system, then some of the IJ interactions

may involve ghost atoms J. The velocities of ghost

atoms are not updated by default. You can enforce that

they will be by using comm_modify vel yes. Your compute

should require that setting be in place, else throw an error.

Some pair styles, like dpd, enforce this, as they also use

the velocities of ghost atoms.

Steve

Thanks Steve,

I already have comm_modify vel yes in my lammps script. I have rechecked the pair_dpd.cpp file to see if something related to comm_modify needs to be enforced in this file, but I didn’t find anything.

Do I need to make any additional changes related to comm_modify in my compute_mygroup_mygroup.cpp file.

Yes I have PBC in all directions.

Thanks again.

Here is the error check:

void PairDPD::init_style()
{
if (comm->ghost_velocity == 0)
error->all(FLERR,“Pair dpd requires ghost atoms store velocity”);

I don’t know why your compute isn’t working as you expect.

But you can easily debug it yourself. There are only a handful

of interactions for 4 atoms. You can print out what it computes

for each one and see why the total is not what you expect.

Again, I’ll point out that if you have 2 type1 and 2 type2 atoms

there are 4 possible type 1/2 interactions. I don’t see why

that should be equal to this variable:

variable ans2 equal
"vx[1]*fx[1]+vx[2]*fx[2]+vx[3]*fx[3]+vx[4]*fx[4]

Steve