Interatomic potential - Indices for the local atoms

LAMMPS LAMMPS Mailing List Mirror
1

Hi,

I am developing an interatomic potential to perform MD simulations at Lammps. I have confusion regarding how energy and forces are updated in the respective pair style’s compute function and especially on multiprocessor systems. I am going through pair_lj_cut.cpp file and will post my questions referencing that potential.

At line 85 of this IP, compute function defines these variables, inum and ilist. From neigh_list.h, they correspond to number of atoms and list of central atoms at each processor.

At line 92, computation starts with this for loop and energy/ forces are updated according to i and j indices below.

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];

Now for different processors, variable ‘i’ remains same for a particular value of ‘ii’. it can be quickly checked by this conditional

if(ii==0){

std::cout<<"i at this processor is "<<i<<std::endl;

}

If the value of i remains same at all the processors, shouldn’t xtmp, ytmp and ztmp have same values across all the processors? Subsequently when we are updating force matrix, by calling f[i], shouldn’t all the processors update only a small set of atoms? I am sure thats not the case, Can you please explain me why its not happening that way?

I have another question regarding tag. When we say itag = tag[i], does it return the id of the atom? Is it a good idea to update forces through tags?

Regards,
Gurjot

2

Hi,

[…]

If the value of i remains same at all the processors, shouldn’t xtmp, ytmp and ztmp have same values across all the processors?

not at all. you are discarding the fact that LAMMPS is a parallel program that uses domain decomposition for parallelization to a large number of processors and for being able to handle very large systems efficiently and with good strong and weak parallel scaling. each processor has “owned” atoms and “ghost” atoms (copies from neighboring domains) and all arrays that represent properties of those atoms refer to the local index. i = 0 to i < atom->natoms addresses “local” atoms and i = atom->nlocal to i < atom->nlocal +atom->nghost to “ghost” atoms. for “ghost” atoms not all arrays are fully populated/allocated.

Subsequently when we are updating force matrix, by calling f[i], shouldn’t all the processors update only a small set of atoms?

no. each processor updates only the properties of its local atoms. for parallelization, then you have two steps of communication: 1) a so-called “reverse” communication, where properties computed and stored on ghost atoms (e.g. when a pair of atoms or bond/angle/dihedral/improper straddles a subdomain boundary and newton’s third law is applied, so each pair’s forces is computed only on one processor) is collected on their corresponding local atoms and 2) a so-called “forward” communication, where the updated positions is sent from the local atoms to their corresponding ghost atoms in neighboring domains. with this layout, communication is needed only to the 6 neighboring subdomains (assuming a “brick” style communication), that can be done efficiently in parallel. if each processor would hold the data to all atoms, you would need collective communications, where the overhead will grow with the number of CPUs used and thus parallelization will only work well for small systems and a small number of CPUs (plus data access with be less cache friendly).

I am sure thats not the case, Can you please explain me why its not happening that way?

see above. for more details, please consult the LAMMPS paper from 1995 and the various presentations about the LAMMPS code structure held at various LAMMPS-related workshops and conferences that are posted on the LAMMPS homepage.

I have another question regarding tag. When we say itag = tag[i], does it return the id of the atom? Is it a good idea to update forces through tags?

no, it would be a very bad idea. all per-atom arrays are indexed by their local index. there is a reverse mapping of tags to local indices, but due to the domain decomposition and distribution of atoms across processors, not all atom coordinates are accessible on all CPUs and thus the mapping function may return an index of -1, indicating that the atom with the given global index was not found.

again, there is plenty of information about this posted on lammps.sandia.gov and in the LAMMPS manual.

axel.

3

Hi Axel,

Thanks a lot for your help. I went through the paper you mentioned and couple of presentations given by Steve and other members. I also followed your advice at this thread https://lammps.sandia.gov/threads/msg13175.html. Unfortunately, I am unable to update the forces correctly with the methods outlined in the references.

I have three body potential with atom indices (i,j,k) and I followed the code structure as written in pair_tersoff.cpp. The force model that I am considering while developing IPs is from this equation

F_i = \sum_{j,k} dE(i,j,k)/dr_i + \sum_{j,k} dE(j,i,k)/dr_i + \sum_{j,k} dE(j,k,i)/dr_i

As you can see, I need to update the force vector for atom ‘i’ where atom ‘j’ or atom ‘k’ can be central atoms. Earlier, I implemented this by using tags which gave me forces corresponding to the model that we have developed. But when I switched to i,j,k obtained in the similar manner as in tersoff file, I am getting totally different forces from the model. Can you please help me in debugging this issue?

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
//localEnergy[i]=0.0;
itag = tag[i];
//i = itag-1; //Change -1
localEnergy[i]=0.0;
x_i = x[i][0]; //CENTRAL ATOM
y_i = x[i][1];
z_i = x[i][2];
itype = type[i];
jlist = firstneigh[i]; //Array of indexes of neighbours of central atom
jnum = numneigh[i]; //intger storing total number of neighbours for i atom
int numshort = 0;
//if(ii==0){
// std::cout<<"i j “<<i<<” "<<jlist[0]<<std::endl;
//}

//two-body loop
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj]; //j-th atom
jtag = tag[j];
//j = jtag-1; //CHANGE - 2
j &= NEIGHMASK; //NO-Clue about this line
if(ii==0){
if(jj==0){
std::cout<<"Inside i j “<<i<<” "<<j<<std::endl;
std::cout<<"itag jtag “<<itag<<” "<<jtag<<std::endl;
}
}
delxij = x_i - x[j][0]; //dr_ij
delyij = y_i - x[j][1]; //dr_ij
delzij = z_i - x[j][2]; //dr_ij
//domain->minimum_image(delxij,delyij,delzij);
rsq = delxijdelxij + delyijdelyij + delzij*delzij;
jtype = type[j];
/if (itag > jtag) {
if ((itag+jtag) 2 == 0) continue; } else if (itag < jtag) { if ((itag+jtag) 2 == 1) continue;
} else {
if (x[j][2] < z_i) continue;
if (x[j][2] == z_i && x[j][1] < y_i) continue;
if (x[j][2] == z_i && x[j][1] == y_i && x[j][0] < x_i) continue;
}/

if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
nIndices = nIndices +1;
neighshort[numshort++] = j;
if (numshort >= maxshort) {
maxshort +=maxshort/2;
memory->grow(neighshort,maxshort,“pair:neighshort”);
}
twobodyDEBUG(paramO, r, delxij,delyij,delzij, energyA,forceAX,forceAY,forceAZ);
}
}
/*if(i==0){

std::cout<<"Len neigshort "<<numshort<<std::endl;
}/
localEnergy[i] = energyATwoBWts;
dummy = forceAXTwoBWts;
f[i][0] += dummy[0,0];
dummy = forceAYTwoBWts;
f[i][1] += dummy[0,0];
dummy = forceAZ*TwoBWts;
f[i][2] += dummy[0,0];
jnumm1 = numshort- 1;
//f_I = arma::zeros(3,thD);
for (jj = 0; jj < jnumm1; jj++) {
//f_J = arma::zeros(3,thD);
j = neighshort[jj];
jtag = tag[j];
//if(j!=jtag){
// std::cout<<"j != jtag “<<std::endl;
// std::cout<<j<<” "<<jtag<<std::endl;
//}

jtype = type[j];
x_j = x[j][0];
y_j = x[j][1];
z_j = x[j][2];
delrij[0] = x_i-x_j;
delrij[1] = y_i-y_j;
delrij[2] = z_i-z_j;
rijSq = delrij[0]*delrij[0] + delrij[1]*delrij[1] + delrij[2]*delrij[2];
for (kk = jj+1; kk < numshort; kk++) {
if(jj==kk) continue;
k = neighshort[kk];
k &= NEIGHMASK;
ktag = tag[k];
ktype = type[k];
x_k = x[k][0];
y_k = x[k][1];
z_k = x[k][2];
delrik[0] = x_i-x_k;
delrik[1] = y_i-y_k;
delrik[2] = z_i-z_k;
//domain->minimum_image(delrik[0],delrik[1],delrik[2]);
delrjk[0] = x_j-x_k;
delrjk[1] = y_j-y_k;
delrjk[2] = z_j-z_k;
//domain->minimum_image(delrjk[0],delrjk[1],delrjk[2]);
rikSq = delrik[0]*delrik[0] + delrik[1]*delrik[1] + delrik[2]*delrik[2];
rjkSq = delrjk[0]*delrjk[0] + delrjk[1]*delrjk[1] + delrjk[2]delrjk[2];
if (rjkSq<=4.0cutsq[itype][jtype]){
rij = sqrt(rijSq);
rjk = sqrt(rjkSq);
rik = sqrt(rikSq);
if (rik>cutoff){
std::cout<<"rik greater than cutoff "<<std::endl;
}
loopTwoStart = clock();
threebodyDEBUG(paramO,rij,rik,rjk,delrij,delrik,delrjk,thEnergyA,f_I,f_J,f_K,insideTHduration);
for(int tt=0;tt<3;tt++){
dummy = arma::dot(f_I.row(tt),ThBWts);
f[i][tt]+=dummy(0,0);
dummy = arma::dot(f_J.row(tt),ThBWts);
//if(j<inum){
f[j][tt]+=dummy(0,0);
//}
//else{
// f[j+nlocal][tt]+=dummy(0,0);
//}
dummy = arma::dot(f_K.row(tt),ThBWts);
f[k][tt]+=dummy(0,0);
}

}
//std::cout<<“DONE WITH ONE K LOOP”<<std::endl;
//f_J = f_J/repmat(ThBrangeLE,3,1);
//for(int tt =0; tt<3;tt++){
// dummy = arma::dot(f_J.row(tt),ThBWts);
// f[j][tt]+=dummy(0,0);
//}
}
//f_I = f_I/repmat(ThBrangeLE,3,1);
//for(int tt =0; tt<3;tt++){
// dummy = arma::dot(f_I.row(tt),ThBWts);
// f[i][tt]+=dummy(0,0);
//}
}
dummy = thEnergyA*ThBWts;
localEnergy[i] = localEnergy[i] + dummy[0,0];

}

4

There is not much that can be debugged with the information provided, and I have no time (or get paid) to do your work. There is plenty of information available and lots of people have figured it out, otherwise there would not be as many contributed pair styles in LAMMPS as there are. LAMMPS is very much a community project and the majority of code in the LAMMPS distribution was written by contributors, not the core LAMMPS developers. Our main job is to maintain and improve or expand the infrastructure and core code and help integrating contributions.

Looking up atoms in the position or force arrays by their tag is incorrect. There is no fixed correlation between the tag and the index in those arrays managed by the AtomVec classes (even if it may look so at the beginning of a serial run). Tags may not be present, tags can be in an arbitrary order when read from a data file, or after a neighbor list update, there may be “holes” in the list of tags so that the largest tag is larger than the total number of atoms, and on a specific MPI rank there may be multiple atoms (as ghosts) present with the same tag. The neighbor lists do index atoms by their local index and nothing else, and also they will always refer to the closest image, which leads to LAMMPS not being subject to minimum image convention requirements. Designing a pair style code around the concept of knowing the data of all atoms is a mistake.

Please also note, that pair styles like tersoff request a so-called “full” neighbor list, where each pair is listed twice (i.e. if an atom A has an atom B as neighbor in its neighbor list, then atom A will also appear neighbor in the list of atom B). With the default “half” neighbor list, this is not the case and the code will then take advantage of newton’s third law and compute the force between those two atoms only once and apply it not just to the first atom, but - with the opposite sign - directly to the second atom as well. for the two-body term in tersoff, the full neighbor list is reduced to a half neighbor list with a chunk of code that you have commented out in the quoted code below. Apart from the incorrect indexing, this is another potential source of error.

Axel.

5

Thanks for your email. It has been very helpful. I have two follow-up questions, answers to which I wasn’t able to find. I hope the answers to them should resolve the implementation issues on my end.

When I am getting j and k indices like below, what values should I expect. For example, I ran the simulation for 200-atoms on a single processor and I got j and k as 393 and 406. Shouldn’t they be less than 200?

Secondly, what is the correct way to update forces on atoms in the neighbour list for ‘i’ i.e. atom ‘j and atom’k’? Is it through f[j] and f[k]?

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
jlist = firstneigh[i]; //Array of indexes of neighbours of central atom
jnum = numneigh[i]; //intger storing total number of neighbours for i atom

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

j = jlist[jj];
jtag = tag[j];
j &= NEIGHMASK

for (kk = 0; kk < jnum; kk++) {
if(jj==kk) continue;
k = jlist[kk];
ktag = tag[k];
k &= NEIGHMASK;
std::cout<<j<<" "<<k<<std::endl;

6

Thanks for your email. It has been very helpful. I have two follow-up questions, answers to which I wasn’t able to find. I hope the answers to them should resolve the implementation issues on my end.

When I am getting j and k indices like below, what values should I expect. For example, I ran the simulation for 200-atoms on a single processor and I got j and k as 393 and 406. Shouldn’t they be less than 200?

no. you have still not understood the concept of “ghost atoms” versus “local atoms” and how per-atom data is managed in general. you must keep in mind, that LAMMPS even treats the single processor case as a parallel case, and thus in the case of just a single MPI rank, the 6 neighboring subdomains (top, bottom, left, right, front, back) are just the same MPI rank, but you communicate “ghost” atoms just as with different MPI ranks. as i pointed out before, minimum image conventions, that are typical for “simple” MD codes, especially with replicated instead of distributed data, do not apply to LAMMPS (except for a couple of unusual and optional features like pair style list).
please go back and re-read my previous answers. everything is explained there or points to where it is explained.

i very much dislike having to explain the same thing to the same person multiple times, so expect that i will ignore further questions unless they address a new issue not yet covered.

Secondly, what is the correct way to update forces on atoms in the neighbour list for ‘i’ i.e. atom ‘j and atom’k’? Is it through f[j] and f[k]?

this should be obvious as soon as you understand the previous point.

axel.