# EDIP forces bug

Hello all,

[This bug is confined completely to pair_edip.cpp]

I was recently writing a version of the EDIP potential for KIM (see pair_style kim documentation page for more information/links about KIM) and was comparing my forces to the forces of the EDIP potential currently implemented in LAMMPS (i.e. pair_style edip). Although my energy matches that of the LAMMPS EDIP, I have ascertained that my forces differ from those of the LAMMPS EDIP in that the signs of the coordination forces are reversed.

I have checked via numerical differentiation of the energy via Ridder’s method that my forces are correct and those of the LAMMPS EDIP are incorrect. Moreover, after working things out again analytically, I still believe that the sign of the coordination forces is incorrect in the LAMMPS EDIP.

For two atoms i and j, let vec_Rij be the vector defined by vec_Rij := Ri-Rj and let Rij be its magnitude. Then, as we loop through all of the atoms and their neighbors, we want to sum up components of f_i as f_i += f_ij where f_ij = - (dE /dRij * dRij/dvec_Rij) while adjusting f_j by f_j -= f_ij in the usual way. In this case, since we have defined vec_Rij := Ri-Rj, we have dRij/dvec_Rij = (vec_Rij)/Rij, whereas if we had used vec_Rij := Rj-Ri, we would get dRij/dvec_Rij = -(vec_Rij)/Rij.

Examining lines 469-479 in the most recent (as of July 25, 2012) version of pair_edip.cpp where the coordination forces are calculated, we see the partial coordination force resulting on atom i from atom j calculated as

f_ij[0] = forceModCoord_ij * dr_ij[0];
f_ij[1] = forceModCoord_ij * dr_ij[1];
f_ij[2] = forceModCoord_ij * dr_ij[2];

f[j][0] -= f_ij[0];
f[j][1] -= f_ij[1];
f[j][2] -= f_ij[2];

f[i][0] += f_ij[0];
f[i][1] += f_ij[1];
f[i][2] += f_ij[2];

In the above code snippet, dr_ij has essentially been defined as dr_ij = Ri-Rj (the factor of 1/Rij has already been accounted for inside of forceModCoord_ij). Thus, there should still be a minus sign on the front of lines 469-471 since it was not accounted for in calculating forceModCoord_ij (you’ll have to check this yourself by working out all of the derivatives analytically).

I suspect this problem went unnoticed because it only becomes evident when there are a substantial number of atoms which are only partial neighbors of one another. This may not necessarily occur when computing e.g. the formation energy of a common defect or a bulk energy.

For comparison purposes, you may want to take a look at Prof. Martin Z. Bazant’s own EDIP code at http://web.mit.edu/bazant/www/EDIP/patch/forces-edip.1.1.c This is what I based my KIM EDIP code on.

I apologize if any of my explanation is unclear. I will promptly respond to feedback on this matter.

Regards,
Daniel S. Karls
University of Minnesota

Sorry, part of what I wrote is nonsense. Please disregard the following:

“For two atoms i and j, let vec_Rij be the vector defined by vec_Rij := Ri-Rj and let Rij be its magnitude. Then, as we loop through all of the atoms and their neighbors, we want to sum up components of f_i as f_i += f_ij where f_ij = - (dE /dRij * dRij/dvec_Rij) while adjusting f_j by f_j -= f_ij in the usual way. In this case, since we have defined vec_Rij := Ri-Rj, we have dRij/dvec_Rij = (vec_Rij)/Rij, whereas if we had used vec_Rij := Rj-Ri, we would get dRij/dvec_Rij = -(vec_Rij)/Rij.”

What I meant to say is that we, of course, compute the partial forces on atom i through f_ij = -(dE/dRij * dRij/dRi). Now, whether you’re using vec_Rij = Ri-Rj or vec_Rij = Rj-Ri, you have dRij/dRi = [Ri_x - Rj_x , Ri_y - Rj_y, Ri_z - Rj_z]/Rij. So when we compute the partial coordination force on atom i, we want

f_ij = - (forceModCoord_ij * [Ri_x - Rj_x , Ri_y - Rj_y, Ri_z - Rj_z]/Rij)

Since, at the point in pair_edip.cpp where the coordination forces are calculated, dr_ij = Ri-Rj, we should still get a minus sign out front, which is missing in the current version of the source.

Sorry for the mistake.

Daniel S. Karls
University of Minnesota

Sounds like you might have found a bug. I'm CCing the
author of pair edip, Luca, to let him weigh in on this ...

Steve

Hi Daniel, I'm currently on vacation. I'll be back at work in August and I'll have a deeper look at the code and my analytical derivation. Many thanks for your work!