Commit 937cf0b9 authored by Markus Hoehnerbach's avatar Markus Hoehnerbach
Browse files

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.
parent 286d4f27

src/MANYBODY/pair_airebo.cpp

+32 −0
Original line number Diff line number Diff line
@@ -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];

src/USER-OMP/pair_airebo_omp.cpp

+32 −0
Original line number Diff line number Diff line
@@ -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];

CRA Git | Maintained and supported by SUSTech CRA and CCSE