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