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 = delxij*delxij + delyij*delyij + 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 = forceAX*TwoBWts;*

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];

}