understanding TIP4P bonded neighbours

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I’m just starting out with LAMMPS, but eventually I’ll need to make a potential somewhat similar to lj/cut/tip4p/long. For now, I’m just trying to understand exactly what pair_lj_cut_tip4p_long.cpp does.

I would be very grateful if someone could answer some (or all!) of my questions.

Before I get to my questions, I will briefly summarize my understanding of what’s going on just to make sure I’m even at the point of knowing what I’m asking.

For a parallel calculation, atoms are assigned to a particular processor. The neighbours of this atom may, or may not, be assigned to the same processor – if not, they are referred to as ghost atoms. Ghost atom positions are updated after they move, and the identities of ghost atoms are updated whenever neighbour lists are updated. Bonded neighbors are usually not included in the neighbour list of a given atom, although presumably bonded atoms will usually be assigned to the same processor (I think not necessarily?). In the TIP4P model, information about bonded hydrogen atoms is required to calculate the position of the negative charge site, which is located off of the oxygen atom. So, when looping over atom pairs, if one of the atoms is identified as oxygen, the index numbers of the bonded hydrogen atoms are obtained using:
hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);

Since there is a “+1” or “+2” in the above lines, atoms in the input data file must be ordered with oxygen first, then the two H atoms. If you wanted to, for whatever reason, have the hydrogen atoms listed first, and the oxygen atom listed last, you would have to change the above code to -1 (or -2) - is that correct?

Does atom->map(atom->tag[j] + 1) get the correct index number for the hydrogen atom whether the hydrogen atom is a ghost atom or not? My understanding is that the hydrogen atom does not need to be a ghost atom or an owned atom for the mapping to be correct - am I right?

I think the line hneigh[j][2] = 1; just indicates that the index numbers associated with the hydrogen atoms bonded to the jth oxygen atom are already known, and don’t need to be gotten with atom->map(atom->tag[j] + 1) again (until the positions change)?

The positions of owned and ghost atoms are already passed through the neighbour list (I’m sure I didn’t say that correctly, hopefully my meaning is clear) such that the minimum distance between the atoms is obtained when calculating, for example, delx = x[j][0] - xtmp;. Since the bonded hydrogen atoms obtained through the mapping described above might not be either ghost or owned atoms,

domain->minimum_image(delx1,dely1,delz1);
is used when calcuating the position of the charge site. Have I understood this correctly (I think minimum images are already handled properly by the neighbour lists for ghost and owned atoms)?

I have more questions, but this email is already too long. I have been printing things out to try to get a handle on everything, but I’m just not sure I’ve understood properly. I would really appreciate any help.

Cheers,

Jeff

2

Hi,

I'm just starting out with LAMMPS, but eventually I'll need to make a
potential somewhat similar to lj/cut/tip4p/long. For now, I'm just trying
to understand exactly what pair_lj_cut_tip4p_long.cpp does.

I would be very grateful if someone could answer some (or all!) of my
questions.
Before I get to my questions, I will briefly summarize my understanding of
what's going on just to make sure I'm even at the point of knowing what I'm
asking.

For a parallel calculation, atoms are assigned to a particular processor.
The neighbours of this atom may, or may not, be assigned to the same
processor -- if not, they are referred to as ghost atoms. Ghost atom

almost. you have local atoms and ghost atoms. independent of neighbor
lists. local atoms are the atoms that are *inside* a given subdomain
at the time of the neighborlist rebuild, ghost atoms are copies(!) of
atoms that are within the distance of the communication cutoff of the
subdomain (usually the largest cutoff plus neighborlist skin). when
neighborlists are build they are first filled with local atoms and
then with ghost atoms, as far as needed according to the cutoff and
neighborlist skin and exclusions.

positions are updated after they move, and the identities of ghost atoms are
updated whenever neighbour lists are updated. Bonded neighbors are usually
not included in the neighbour list of a given atom, although presumably
bonded atoms will usually be assigned to the same processor (I think not

this depends on the newton_bond flag. for newton bond yes, bonded
forces are computed once and then forces on the ghost atoms are then
communicated back to their originals.

necessarily?). In the TIP4P model, information about bonded hydrogen atoms
is required to calculate the position of the negative charge site, which is
located off of the oxygen atom. So, when looping over atom pairs, if one of
the atoms is identified as oxygen, the index numbers of the bonded hydrogen
atoms are obtained using:
   hneigh[j][0] = jH1 = atom->map(atom->tag[j] + 1);
   hneigh[j][1] = jH2 = atom->map(atom->tag[j] + 2);
Since there is a "+1" or "+2" in the above lines, atoms in the input data

file must be ordered with oxygen first, then the two H atoms. If you wanted
to, for whatever reason, have the hydrogen atoms listed first, and the
oxygen atom listed last, you would have to change the above code to -1 (or
-2) - is that correct?

yes.

Does atom->map(atom->tag[j] + 1) get the correct index number for the
hydrogen atom whether the hydrogen atom is a ghost atom or not? My

atom->tag[] returns the global index (what you assign in the data
file) of an atom, local or ghost.
atom->map() maps the global index to the local index or returns -1, if
the atom is neither local nor ghost.

understanding is that the hydrogen atom does not need to be a ghost atom or
an owned atom for the mapping to be correct - am I right?

no. see above.

I think the line hneigh[j][2] = 1; just indicates that the index numbers
associated with the hydrogen atoms bonded to the jth oxygen atom are already
known, and don't need to be gotten with atom->map(atom->tag[j] + 1) again
(until the positions change)?

yes. this pair style has seen some significant tuning and optimization
and this is one of the tricks used.

The positions of owned and ghost atoms are already passed through the
neighbour list (I'm sure I didn't say that correctly, hopefully my meaning
is clear) such that the minimum distance between the atoms is obtained when
calculating, for example, delx = x[j][0] - xtmp;. Since the bonded hydrogen
atoms obtained through the mapping described above might not be either ghost
or owned atoms,
  domain->minimum_image(delx1,dely1,delz1);
is used when calcuating the position of the charge site. Have I understood
this correctly (I think minimum images are already handled properly by the
neighbour lists for ghost and owned atoms)?

with few processors, a small system and a long cutoff, there can be
several copies of an atom with the same global index (aka "tag"). one
cannot know, if atom->map() returns the correct one, hence the minimum
image distance is enforced

HTH,
    axel.