Commit 0e16dc3e authored by Markus Hoehnerbach's avatar Markus Hoehnerbach
Browse files

In PairAIREBO::bondorderLJ: p^sigma pi account for d/drij derivatives. 

The bonderorderLJ function operates on a facticious distance |rij|, i.e. everything gets calculated "as if" atoms i and j were a given distance alpha apart.
Mathematically, bondorderLJ is a function of rij (a vector), that is (in terms of the real distance Rij) rij = alpha * Rij/|Rij|.
When we calculate the forces in bondorderLJ, we have to make sure to chain in this derivative whenever we calculate derivatives w.r.t. rij.
The right correction, as it turns our, is Fij = alpha / |Rij| * (Identity(3,3) - Rij * Rij^T / |Rij|^2) * fij.
This commit only fixes this for the p_ij^sigma pi terms, which were modified to separate out the d/drij derivative in the cosine calculation.
Now, derivatives are taken w.r.t. the connecting edges instead of the edge points.
parent 1b3f6e25

src/MANYBODY/pair_airebo.cpp

+88 −73
Original line number Diff line number Diff line
@@ -2130,6 +2130,12 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, 
  Stb = 0.0;

  dStb = 0.0;



  double realrij[3];

  realrij[0] = x[atomi][0] - x[atomj][0];

  realrij[1] = x[atomi][1] - x[atomj][1];

  realrij[2] = x[atomi][2] - x[atomj][2];

  double realrijmag = sqrt(realrij[0] * realrij[0] + realrij[1] * realrij[1] + realrij[2] * realrij[2]);



  REBO_neighs = REBO_firstneigh[i];

  for (k = 0; k < REBO_numneigh[i]; k++) {

    atomk = REBO_neighs[k];
@@ -2319,50 +2325,55 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, 
        lamdajik = 4.0*kronecker(itype,1) *

          ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag));

        wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik);

        cosjik = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) /

          (rijmag*rikmag);

        double rjk[3];

        rjk[0] = rik[0] - rij[0];

        rjk[1] = rik[1] - rij[1];

        rjk[2] = rik[2] - rij[2];

        double rjkmag = sqrt(rjk[0] * rjk[0] + rjk[1] * rjk[1] + rjk[2] * rjk[2]);

        double rijrik = 2 * rijmag * rikmag;

        double rr = rijmag * rijmag - rikmag * rikmag;

        cosjik = (rr + rikmag * rikmag - rjkmag * rjkmag) / rijrik;

        cosjik = MIN(cosjik,1.0);

        cosjik = MAX(cosjik,-1.0);



        dcosjikdri[0] = ((rij[0]+rik[0])/(rijmag*rikmag)) -

          (cosjik*((rij[0]/(rijmag*rijmag))+(rik[0]/(rikmag*rikmag))));

        dcosjikdri[1] = ((rij[1]+rik[1])/(rijmag*rikmag)) -

          (cosjik*((rij[1]/(rijmag*rijmag))+(rik[1]/(rikmag*rikmag))));

        dcosjikdri[2] = ((rij[2]+rik[2])/(rijmag*rikmag)) -

          (cosjik*((rij[2]/(rijmag*rijmag))+(rik[2]/(rikmag*rikmag))));

        dcosjikdrk[0] = (-rij[0]/(rijmag*rikmag)) +

          (cosjik*(rik[0]/(rikmag*rikmag)));

        dcosjikdrk[1] = (-rij[1]/(rijmag*rikmag)) +

          (cosjik*(rik[1]/(rikmag*rikmag)));

        dcosjikdrk[2] = (-rij[2]/(rijmag*rikmag)) +

          (cosjik*(rik[2]/(rikmag*rikmag)));

        dcosjikdrj[0] = (-rik[0]/(rijmag*rikmag)) +

          (cosjik*(rij[0]/(rijmag*rijmag)));

        dcosjikdrj[1] = (-rik[1]/(rijmag*rikmag)) +

          (cosjik*(rij[1]/(rijmag*rijmag)));

        dcosjikdrj[2] = (-rik[2]/(rijmag*rikmag)) +

          (cosjik*(rij[2]/(rijmag*rijmag)));

        double dctdjk = -2 / rijrik;

        double dctdik = (-rr + rjkmag * rjkmag) / (rijrik * rikmag * rikmag);

        // evaluate splines g and derivatives dg



        g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN);



        tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik));

        fj[0] = -tmp2*dcosjikdrj[0];

        fj[1] = -tmp2*dcosjikdrj[1];

        fj[2] = -tmp2*dcosjikdrj[2];

        fi[0] = -tmp2*dcosjikdri[0];

        fi[1] = -tmp2*dcosjikdri[1];

        fi[2] = -tmp2*dcosjikdri[2];

        fk[0] = -tmp2*dcosjikdrk[0];

        fk[1] = -tmp2*dcosjikdrk[1];

        fk[2] = -tmp2*dcosjikdrk[2];



        fi[0] =  tmp2 * dctdik * rik[0];

        fi[1] =  tmp2 * dctdik * rik[1];

        fi[2] =  tmp2 * dctdik * rik[2];

        fk[0] = -tmp2 * dctdik * rik[0];

        fk[1] = -tmp2 * dctdik * rik[1];

        fk[2] = -tmp2 * dctdik * rik[2];

        fj[0] = 0;

        fj[1] = 0;

        fj[2] = 0;

        double fjk[3];

        fjk[0] = tmp2 * dctdjk * rjk[0];

        fjk[1] = tmp2 * dctdjk * rjk[1];

        fjk[2] = tmp2 * dctdjk * rjk[2];

        fi[0] += fjk[0];

        fi[1] += fjk[1];

        fi[2] += fjk[2];

        fk[0] -= fjk[0];

        fk[1] -= fjk[1];

        fk[2] -= fjk[2];

        double rijmbr = rcmin[itype][jtype] / realrijmag;

        fj[0] += rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));

        fj[1] += rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));

        fj[2] += rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));

        fi[0] -= rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));

        fi[1] -= rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));

        fi[2] -= rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag));



        tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1));

        fj[0] -= tmp2*(-rij[0]/rijmag);

        fj[1] -= tmp2*(-rij[1]/rijmag);

        fj[2] -= tmp2*(-rij[2]/rijmag);

        fi[0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag));

        fi[1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag));

        fi[2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag));

        fi[0] += tmp2*(rik[0]/rikmag);

        fi[1] += tmp2*(rik[1]/rikmag);

        fi[2] += tmp2*(rik[2]/rikmag);

        fk[0] -= tmp2*(rik[0]/rikmag);

        fk[1] -= tmp2*(rik[1]/rikmag);

        fk[2] -= tmp2*(rik[2]/rikmag);
@@ -2427,51 +2438,55 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, 
        lamdaijl = 4.0*kronecker(jtype,1) *

          ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag));

        wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl);

        cosijl = (-1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2]))) /

          (rijmag*rjlmag);

        double ril[3];

        ril[0] = rij[0] + rjl[0];

        ril[1] = rij[1] + rjl[1];

        ril[2] = rij[2] + rjl[2];

        double rilmag = sqrt(ril[0] * ril[0] + ril[1] * ril[1] + ril[2] * ril[2]);

        double rijrjl = 2 * rijmag * rjlmag;

        double rr = rijmag * rijmag - rjlmag * rjlmag;

        cosijl = (rr + rjlmag * rjlmag - rilmag * rilmag) / rijrjl;

        cosijl = MIN(cosijl,1.0);

        cosijl = MAX(cosijl,-1.0);



        dcosijldri[0] = (-rjl[0]/(rijmag*rjlmag)) -

          (cosijl*rij[0]/(rijmag*rijmag));

        dcosijldri[1] = (-rjl[1]/(rijmag*rjlmag)) -

          (cosijl*rij[1]/(rijmag*rijmag));

        dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) -

          (cosijl*rij[2]/(rijmag*rijmag));

        dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) +

          (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag))));

        dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) +

          (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag))));

        dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) +

          (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag))));

        dcosijldrl[0] = (rij[0]/(rijmag*rjlmag)) +

          (cosijl*rjl[0]/(rjlmag*rjlmag));

        dcosijldrl[1] = (rij[1]/(rijmag*rjlmag)) +

          (cosijl*rjl[1]/(rjlmag*rjlmag));

        dcosijldrl[2] = (rij[2]/(rijmag*rjlmag)) +

          (cosijl*rjl[2]/(rjlmag*rjlmag));

        double dctdil = -2 / rijrjl;

        double dctdjl = (-rr + rilmag * rilmag) / (rijrjl * rjlmag * rjlmag);



        // evaluate splines g and derivatives dg



        g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN);

        tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl));

        fi[0] = -tmp2*dcosijldri[0];

        fi[1] = -tmp2*dcosijldri[1];

        fi[2] = -tmp2*dcosijldri[2];

        fj[0] = -tmp2*dcosijldrj[0];

        fj[1] = -tmp2*dcosijldrj[1];

        fj[2] = -tmp2*dcosijldrj[2];

        fl[0] = -tmp2*dcosijldrl[0];

        fl[1] = -tmp2*dcosijldrl[1];

        fl[2] = -tmp2*dcosijldrl[2];

        fj[0] =  tmp2 * dctdjl * rjl[0];

        fj[1] =  tmp2 * dctdjl * rjl[1];

        fj[2] =  tmp2 * dctdjl * rjl[2];

        fl[0] = -tmp2 * dctdjl * rjl[0];

        fl[1] = -tmp2 * dctdjl * rjl[1];

        fl[2] = -tmp2 * dctdjl * rjl[2];

        fi[0] = 0;

        fi[1] = 0;

        fi[2] = 0;

        double fil[3];

        fil[0] = tmp2 * dctdil * ril[0];

        fil[1] = tmp2 * dctdil * ril[1];

        fil[2] = tmp2 * dctdil * ril[2];

        fj[0] += fil[0];

        fj[1] += fil[1];

        fj[2] += fil[2];

        fl[0] -= fil[0];

        fl[1] -= fil[1];

        fl[2] -= fil[2];

        double rijmbr = rcmin[itype][jtype] / realrijmag;

        fi[0] += rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag));

        fi[1] += rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag));

        fi[2] += rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag));

        fj[0] -= rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag));

        fj[1] -= rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag));

        fj[2] -= rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag));



 

        tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1));

        fi[0] -= tmp2*(rij[0]/rijmag);

        fi[1] -= tmp2*(rij[1]/rijmag);

        fi[2] -= tmp2*(rij[2]/rijmag);

        fj[0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag));

        fj[1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag));

        fj[2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag));

        fj[0] += tmp2*(rjl[0]/rjlmag);

        fj[1] += tmp2*(rjl[1]/rjlmag);

        fj[2] += tmp2*(rjl[2]/rjlmag);

        fl[0] -= tmp2*(rjl[0]/rjlmag);

        fl[1] -= tmp2*(rjl[1]/rjlmag);

        fl[2] -= tmp2*(rjl[2]/rjlmag);

CRA Git | Maintained and supported by SUSTech CRA and CCSE