Reading and using bond information

LAMMPS LAMMPS Development
1

Hi,

I’m looking at LAMMPS source code and trying to figure out how LAMMPS reads in the bond information. Specifically, does LAMMPS read in bond information on different processors and store it there (i.e. each processor only has part of the bond information)? In atom.cpp, data_bonds function seems to build bond_type and bond_atom arrays based on the local particles on each processor as it indexes through m=map(atom1) that is the local ID, and the code only does this build once during read_data. But in ntopo_bond_all.cpp, my understanding is that bondlist is built using the new particle distributions on the processors (i.e. new nlocals) with the old bond_atom and bond_type that are generated when reading in from the data file. Are bond_atom and bond_type updated elsewhere after they are generated with read_data?

Thanks.

2

This is a complex question that requires multiple chunks of an answer.

  1. Bonded atom information is stored with atoms: this applies to bonds, angles, dihedrals, and impropers. How, depends on the Force::newton_bond setting: is it 1, then the information is stored with one of the atoms; is it 0, then the information is stored with all of the participating atoms.
  2. Bonded atom lists are generated by the neighbor list code from the information stored with atoms, after the “local” atoms have been redistributed into sub-domains and all “ghost” atoms recreated. Thus as the atoms migrate between sub-domains, so does the bonded interaction information. This causes the requirement, that there are sufficient number of ghost atoms (up to the communication cutoff, which is usually the same as the neighbor list cutoff) so that for all bonded interactions that straddle subdomain boundaries, all participating atoms can be found.
  3. The read_data command reads blocks of lines from the data file on MPI rank 0 and then broadcasts them to all participating MPI processes. Then the selected AtomVec class (from the atom style command) processes them and only keeps information that is relevant to the atoms at the individual MPI rank. See: 4.14. Utility functions — LAMMPS documentation
  4. The selected AtomVec class determines which per-atom information attached to atoms is communicated in the various communication patterns (border, forward, reverse, exchange).

If you need additional details, please check out: 4. Information for Developers — LAMMPS documentation and the LAMMMPS paper.

3

Hi Axel,

Thanks for your reply. The information here is very helpful, especially Point 3. I understand bonded atom info is stored with atom (presumably stored in bond_atom and bond_type), but I’m still confused about the implementation. Let me be more specific about my understanding of the general steps in saving the bonded info. You can point out which part I misunderstand.

Step 1:
In read_data.cpp, utils::read_lines_from_file reads in all bonded information on rank 0 from the data file and broadcast it to all MPI ranks (as you said).

Step 2:
atom->data_bonds is called to save bond info to all particles owned by each MPI rank. I’m looking at the following section in atom.cpp:

if ((m = map(atom1)) >= 0) {
  if (count) count[m]++; 
  else {
    bond_type[m][num_bond[m]] = itype;
    bond_atom[m][num_bond[m]] = atom2;
    num_bond[m]++; 
    avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset); 
  } 
}  
if (newton_bond == 0) {
  if ((m = map(atom2)) >= 0) {
    if (count) count[m]++; 
    else {
      bond_type[m][num_bond[m]] = itype;
      bond_atom[m][num_bond[m]] = atom1;
      num_bond[m]++;
      avec->data_bonds_post(m, num_bond[m], atom1, atom2, id_offset);
    }
  }
}

This section looks like it stores bonded info based on newton_bond. So say I’m on rank 1 that has nlocal_1 particles at the time of data reading, and m is the local ID that is mapped from global ID atom1 and loops over all particles on rank 1; i.e. bond_type and bond_atom are built for those nlocal_1 particles on rank 1 (and are saved on rank 1?).

Somewhere between Step 2 and 3, the time advancing loop starts.

Step 3:
At the time of data reading, the neighbor list is generated, and the topology is also built, through neigh_bond->build(). neigh_bond is a NTopo class, and I’m looking at ntopo_bond_all.cpp for build() method:

for (i = 0; i < nlocal; i++)
    for (m = 0; m < num_bond[i]; m++) {
        atom1 = atom->map(bond_atom[i][m]);
    ......
        atom1 = domain->closest_image(i, atom1);
        if (newton_bond || i < atom1) {
          if (nbondlist == maxbond) {
            maxbond += DELTA;
            memory->grow(bondlist, maxbond, 3, "neigh_topo:bondlist");
          }
          bondlist[nbondlist][0] = i;
          bondlist[nbondlist][1] = atom1;
          bondlist[nbondlist][2] = bond_type[i][m];
          nbondlist++;
      }

It builds bondlist and nbondlist using atom info and neighbor info on each rank. Say I’m still on rank 1, and at the time of data reading, nlocal in Step 3 is equal to nlocal_1 from Step 2. So all good.

Step 4:
bondlist and nbondlist are used in various Bond::compute() methods to calculate forces and update attributes in the atom class.

And somewhere here the time advancing loop ends.

Here is my question for Step 3. At a later time when particles move to new ranks and new neighbor lists are built, nlocal may represent different particles. Still using rank 1 as an example. Now rank 1 own nlocal_1_later particles, which in an extreme case can be completely different from the nlocal_1 particles in Step 2. As a result, Step 3 loops through nlocal_1_later particles. However, bond_atom is not updated, so the first index of bond_atom should represent nlocal_1 particles from Step 2; i.e. i in atom1 = atom->map(bond_atom[i][m]) should only loop through nlocal_1 particles, not nlocal_1_later particles at a later time. So what I see is a mismatched bond info in bond_atom and bond_type. What am I missing here? Thanks.

4

The bond information stays with local atoms. If an atom migrates from one subdomain to another, it takes all its per-atom information with it. Please see the various “pack” and “unpack” methods in the AtomVec classes corresponding to the atom styles. Lets take atom style “bond”. Have a look at the file src/MOLECULE/atom_vec_bond.cpp in its constructor it sets up the list of properties to be communicated: for the exchange step those are "molecule, num_bond, bond_type, bond_atom, nspecial, and special. Comm::exchange() is called during the reneighboring.

Thus, at any time, when LAMMPS decides to (re-)build the neighbor lists. It first deletes the bond list and then reconstructs it from the information stored with the local atoms.

5

Thank you for your help. I think atom_vec_bond.cpp and comm.cpp are what I need.

6

Need to do what?

7

I’m developing an inhouse code for a hybrid (particles+fluids) model. I just want to understand how LAMMPS deals with particles.