Atom->map(global_id) always returning -1 in custom compute

LAMMPS LAMMPS Development
1

Dear LAMMPS users and developers,

I am working on a custom compute function where I need to calculate the displacement of an atom relative to the center of its neighboring atoms. For instance, I want to compute the displacement vector of a carbon atom in a methane (CH₄) molecule, shifting away from the center of the CH₄ tetrahedral structure. The neighbor list is provided by an external file, which records the global IDs of both the central and neighboring atoms.

However, I’ve encountered a problem: when I use atom->map(global_id) to retrieve the local index of an atom from its global ID, it consistently returns -1. This occurs even though I am running LAMMPS with only one processor (also, I have set OMP_NUM_THREADS=1).

Here is the core part of the code:

void ComputeCustomDisp::compute_peratom() 
{
    // check number of atoms
    if (atom->nmax > nmax) {
        memory->destroy(array_atom);
        nmax = atom->nmax;
        memory->create(array_atom, nmax, size_peratom_cols, "custom_disp:array_atom");
    }

    // reset the array_atom to 0
    for (int i = 0; i < atom->natoms; i++) {
        array_atom[i][0] = 0.0;
        array_atom[i][1] = 0.0;
        array_atom[i][2] = 0.0;
    }

    // assign the coord
    double **x = atom->x;
    
    // loop over all central atoms
    size_t i,j;
    for(i=0; i < central_id.size(); i++) {
        int central_global_id = central_id[i];
        int central_local_id = atom->map(central_global_id);
        printf("central_global_id: %d\n", central_global_id);
        printf("central_local_id: %d\n", central_local_id);

        if (central_local_id < 0) continue; // skip if the central atom isn't in the processor

        // loop over all neighbors of the central atom
        double xn=0.0, yn=0.0, zn=0.0;
        int neighbor_count = 0;
        for(j = 0; j < neighbor_id[i].size(); j++) {
            int neighbor_global_id = neighbor_id[i][j];
            int neighbor_local_id = atom->map(neighbor_global_id);
            printf("neighbor_global_id: %d\n", neighbor_global_id);
            printf("neighbor_local_id: %d\n", neighbor_local_id);

            if (neighbor_local_id < 0) continue; // skip if the neighbor atom isn't in the processor

            xn += x[neighbor_local_id][0];
            yn += x[neighbor_local_id][1];
            zn += x[neighbor_local_id][2];
            neighbor_count++;
        }

        xn /= neighbor_count;
        yn /= neighbor_count;
        zn /= neighbor_count;

        double dx = x[central_local_id][0] - xn;
        double dy = x[central_local_id][1] - yn;
        double dz = x[central_local_id][2] - zn;

        array_atom[central_local_id][0] = dx;
        array_atom[central_local_id][1] = dy;
        array_atom[central_local_id][2] = dz;
    }
}

The issue is that for every central atom, the atom->map() function always returns -1, even though I am running on a single processor. Below is the output from my code:

central_global_id: 5
central_local_id: -1
central_global_id: 6
central_local_id: -1
central_global_id: 7
central_local_id: -1
central_global_id: 8
central_local_id: -1
central_global_id: 9
central_local_id: -1
central_global_id: 10
central_local_id: -1
central_global_id: 11
central_local_id: -1
central_global_id: 12
central_local_id: -1

I’m new to both C++ and LAMMPS development, so I would greatly appreciate any advice or suggestions to resolve this issue. Specifically, why is atom->map() returning -1 even when running with a single processor?

Thank you in advance for your help!

Best regards,
denan

2

This is wrong. The loop should be over atom->nlocal or atom->nlocal+atom->nghost.

Does your system actually have a map? By default a map is only created, if it is needed.
It seems atom->map_style is Atom::MAP_NONE in your case.

3

Thank you for your suggestion, Professor!

Everything is working perfectly now after I added the atom->map_init() call.

4

Hi axel!

After setting atom->map_init(), my simulations started running smoothly. However, I’m facing a new issue:

In my custom MD simulations, I will run multiple loops until a certain condition is met. Here is key part of my input:

variable        j loop ${nloop}
label           loop_mark
 
... # some variable define

thermo_style    custom step temp etotal pe ke vol lx ly lz
run             20

... # some variable define 

if "some condition" then "jump  SELF break"
next            j
jump            SELF loop_mark

label           break

I’ve noticed that each loop takes longer to complete over time.

In my original implementation, I called atom->init() in void ComputeCustomDisp::init(), and when looping over all atoms in void ComputeCustomDisp::compute_peratom() function, I used atom->map(global index) to get local index on-the-fly:

void ComputeCustomDisp::init()
{
    // read the neighbor list file, initialize atom->map()
    read_file();
    atom->map_init();
}

void ComputeCustomDisp::compute_peratom() 
{   

    // check number of atoms
    if (atom->nmax > nmax) {
        memory->destroy(array_atom);
        nmax = atom->nmax;
        memory->create(array_atom, nmax, size_peratom_cols, "custom_disp:array_atom");
    }

    // reset the array_atom to 0
    for (int i = 0; i < atom->nlocal + atom->nghost; i++) {
        array_atom[i][0] = 0.0;
        array_atom[i][1] = 0.0;
        array_atom[i][2] = 0.0;
    }

    // assign the coord
    double **x = atom->x;
    
    // loop over all central atoms
    size_t i,j;
    for(i=0; i < central_id.size(); i++) {
        int central_global_id = central_id[i];
        int central_local_id = atom->map(central_global_id);

        if (central_local_id < 0 || central_local_id > atom->nlocal) continue; // skip if the central atom isn't in the processor

        // loop over all neighbors of the central atom
        double dx=0, dy=0, dz=0;
        double tmpx, tmpy, tmpz;
        int neighbor_count = 0;
        for(j = 0; j < neighbor_id[i].size(); j++) {
            int neighbor_global_id = neighbor_id[i][j];
            int neighbor_local_id = atom->map(neighbor_global_id);

            if (neighbor_local_id < 0) {
                error->warning(FLERR, "Cannot find the local index of the neighbor atom");
                continue;
            }

            tmpx = x[central_local_id][0] - x[neighbor_local_id][0];
            tmpy = x[central_local_id][1] - x[neighbor_local_id][1];
            tmpz = x[central_local_id][2] - x[neighbor_local_id][2];
            domain->minimum_image(tmpx, tmpy, tmpz);
            dx += tmpx;
            dy += tmpy;
            dz += tmpz;
            neighbor_count++;
        }
        dx /= neighbor_count;
        dy /= neighbor_count;
        dz /= neighbor_count;

        array_atom[central_local_id][0] = dx;
        array_atom[central_local_id][1] = dy;
        array_atom[central_local_id][2] = dz;
    }
}

To try to fix this, I stored local indexes in a vector during the void ComputeCustomDisp::init() and used them directly in void ComputeCustomDisp::compute_peratom(), which seemed to solve the problem.

According to LAMMPS’s domain decomposition scheme, the atoms in a certain region are held by a processor. Are atoms at the boundaries of a domain reassigned to different processors as they move across the boundaries? or are they attached to a fixed processor throughout the simulation?

Could you offer any insights or suggestions on why the original implementation will slow the simulations? If the atoms held by a processor are changing dynamically, my new solution will fail, is there any more effective to handle it?

Best regards,
Denan

5

Please create a new topic for new issues with a corresponding subject line.