From 937cf0b996e5a31ebaeecbe412f03a4abaa508e0 Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 31 May 2017 12:20:12 +0200 Subject: [PATCH] Bugfix: Kronecker term ignored in spline forces. The code ignored the kronecker(ktype, 0) or kronecker(ltype, 0) terms in the contributing terms to NconjtmpI and NconjtmpJ. The issue was present both in ::bondorder and ::bondorderLJ and led to energy conservation issues. It has been fixed by checking for the atom type before entering the offending calculations and adding clarifying comments. --- src/MANYBODY/pair_airebo.cpp | 32 ++++++++++++++++++++++++++++++++ src/USER-OMP/pair_airebo_omp.cpp | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index cc7efbcaa6..d83f5a39a8 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -1615,6 +1615,10 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1678,6 +1682,10 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1960,6 +1968,10 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2023,6 +2035,10 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2560,6 +2576,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2623,6 +2643,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2895,6 +2919,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2958,6 +2986,10 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index 95f9d8b401..206e8e86e6 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -1387,6 +1387,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1450,6 +1454,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1732,6 +1740,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1795,6 +1807,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2332,6 +2348,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2395,6 +2415,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2667,6 +2691,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr); + // due to kronecker(ktype, 0) term in contribution + // to NconjtmpI and later Nijconj + if (ktype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2730,6 +2758,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr); + // due to kronecker(ltype, 0) term in contribution + // to NconjtmpJ and later Nijconj + if (ltype != 0) continue; + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1];