diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 3db5c808e7..0751e66a19 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -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);