diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 9b84a6c0f8..eb60f5c684 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -71,28 +71,29 @@ PairLJCutCoulLongTIP4P::~PairLJCutCoulLongTIP4P() void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype,itable; + int n,vlist[6]; + int iH1,iH2,jH1,jH2; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; double fraction,table; double delx1,dely1,delz1,delx2,dely2,delz2,delx3,dely3,delz3; double r,r2inv,r6inv,forcecoul,forcelj,cforce,negforce; double factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; - int iH1,iH2,jH1,jH2; - double xiM[3],xjM[3]; + double xiM[3],xjM[3],fO[3],fH[3],v[6]; double *x1,*x2; - double fO[3],fH[3]; int *ilist,*jlist,*numneigh,**firstneigh; float rsq; int *int_rsq = (int *) &rsq; + // if vflag_global = 2, reset vflag as if vflag_global = 1 + // necessary since TIP4P cannot compute virial as F dot r + // due to find_M() finding bonded H atoms which are not near O atom + + if (vflag % 4 == 2) vflag = 1 + vflag/4 * 4; + evdwl = ecoul = 0.0; if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style lj/cut/coul/long/tip4p does not yet support peratom energy/virial"); + else evflag = 0; double **f = atom->f; double **x = atom->x; @@ -105,23 +106,6 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; - // grow and zero temporary force array if necessary - - if (vflag && atom->nmax > nmax) { - memory->destroy_2d_double_array(ftmp); - nmax = atom->nmax; - ftmp = memory->create_2d_double_array(nmax,3,"pair:ftmp"); - } - - if (vflag) { - for (i = 0; i < 6; i++) virialtmp[i] = 0.0; - for (i = 0; i < nall; i++) { - ftmp[i][0] = 0.0; - ftmp[i][1] = 0.0; - ftmp[i][2] = 0.0; - } - } - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -179,8 +163,11 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - offset[itype][jtype]; evdwl *= factor_lj; - } - } else evdwl = 0.0; + } else evdwl = 0.0; + + if (evflag) ev_tally(i,j,nlocal,newton_pair, + evdwl,0.0,forcelj,delx,dely,delz); + } // adjust rsq for off-site O charge(s) @@ -227,27 +214,25 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) cforce = forcecoul * r2inv; // if i,j are not O atoms, force is applied directly - // if i or j are O atoms, force is on fictitious atoms - // spread force to all 3 atoms in water molecule + // if i or j are O atoms, force is on fictitious atom + // and is thus spread to all 3 atoms in water molecule // formulas due to Feenstra et al, J Comp Chem, 20, 786 (1999) + n = 0; + if (itype != typeO) { - if (vflag == 0) { - f[i][0] += delx * cforce; - f[i][1] += dely * cforce; - f[i][2] += delz * cforce; + f[i][0] += delx * cforce; + f[i][1] += dely * cforce; + f[i][2] += delz * cforce; - } else { - ftmp[i][0] += delx * cforce; - ftmp[i][1] += dely * cforce; - ftmp[i][2] += delz * cforce; - - virialtmp[0] += 0.5 * delx * delx * cforce; - virialtmp[1] += 0.5 * dely * dely * cforce; - virialtmp[2] += 0.5 * delz * delz * cforce; - virialtmp[3] += 0.5 * dely * delx * cforce; - virialtmp[4] += 0.5 * delz * delx * cforce; - virialtmp[5] += 0.5 * delz * dely * cforce; + if (vflag) { + v[0] = 0.5 * delx * delx * cforce; + v[1] = 0.5 * dely * dely * cforce; + v[2] = 0.5 * delz * delz * cforce; + v[3] = 0.5 * delx * dely * cforce; + v[4] = 0.5 * delx * delz * cforce; + v[5] = 0.5 * dely * delz * cforce; + vlist[n++] = i; } } else { @@ -259,32 +244,19 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) fH[1] = alpha * (dely*cforce); fH[2] = alpha * (delz*cforce); - if (vflag == 0) { - f[i][0] += fO[0]; - f[i][1] += fO[1]; - f[i][2] += fO[2]; + f[i][0] += fO[0]; + f[i][1] += fO[1]; + f[i][2] += fO[2]; - f[iH1][0] += fH[0]; - f[iH1][1] += fH[1]; - f[iH1][2] += fH[2]; + f[iH1][0] += fH[0]; + f[iH1][1] += fH[1]; + f[iH1][2] += fH[2]; - f[iH2][0] += fH[0]; - f[iH2][1] += fH[1]; - f[iH2][2] += fH[2]; - - } else { - ftmp[i][0] += fO[0]; - ftmp[i][1] += fO[1]; - ftmp[i][2] += fO[2]; - - ftmp[iH1][0] += fH[0]; - ftmp[iH1][1] += fH[1]; - ftmp[iH1][2] += fH[2]; - - ftmp[iH2][0] += fH[0]; - ftmp[iH2][1] += fH[1]; - ftmp[iH2][2] += fH[2]; + f[iH2][0] += fH[0]; + f[iH2][1] += fH[1]; + f[iH2][2] += fH[2]; + if (vflag) { delx1 = x[i][0] - x2[0]; dely1 = x[i][1] - x2[1]; delz1 = x[i][2] - x2[2]; @@ -300,32 +272,33 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) delz3 = x[iH2][2] - x2[2]; domain->minimum_image(delx3,dely3,delz3); - virialtmp[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); - virialtmp[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); - virialtmp[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); - virialtmp[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); - virialtmp[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); - virialtmp[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + v[0] = 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); + v[1] = 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); + v[2] = 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); + v[3] = 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); + v[4] = 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); + v[5] = 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + + vlist[n++] = i; + vlist[n++] = iH1; + vlist[n++] = iH2; } } if (jtype != typeO) { - if (vflag == 0) { - f[j][0] -= delx * cforce; - f[j][1] -= dely * cforce; - f[j][2] -= delz * cforce; + f[j][0] -= delx * cforce; + f[j][1] -= dely * cforce; + f[j][2] -= delz * cforce; - } else { - ftmp[j][0] -= delx * cforce; - ftmp[j][1] -= dely * cforce; - ftmp[j][2] -= delz * cforce; + if (vflag) { + v[0] += 0.5 * delx * delx * cforce; + v[1] += 0.5 * dely * dely * cforce; + v[2] += 0.5 * delz * delz * cforce; + v[3] += 0.5 * delx * dely * cforce; + v[4] += 0.5 * delx * delz * cforce; + v[5] += 0.5 * dely * delz * cforce; - virialtmp[0] += 0.5 * delx * delx * cforce; - virialtmp[1] += 0.5 * dely * dely * cforce; - virialtmp[2] += 0.5 * delz * delz * cforce; - virialtmp[3] += 0.5 * dely * delx * cforce; - virialtmp[4] += 0.5 * delz * delx * cforce; - virialtmp[5] += 0.5 * delz * dely * cforce; + vlist[n++] = j; } } else { @@ -339,32 +312,19 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) fH[1] = alpha * (dely*negforce); fH[2] = alpha * (delz*negforce); - if (vflag != 2) { - f[j][0] += fO[0]; - f[j][1] += fO[1]; - f[j][2] += fO[2]; + f[j][0] += fO[0]; + f[j][1] += fO[1]; + f[j][2] += fO[2]; - f[jH1][0] += fH[0]; - f[jH1][1] += fH[1]; - f[jH1][2] += fH[2]; + f[jH1][0] += fH[0]; + f[jH1][1] += fH[1]; + f[jH1][2] += fH[2]; - f[jH2][0] += fH[0]; - f[jH2][1] += fH[1]; - f[jH2][2] += fH[2]; - - } else { - ftmp[j][0] += fO[0]; - ftmp[j][1] += fO[1]; - ftmp[j][2] += fO[2]; - - ftmp[jH1][0] += fH[0]; - ftmp[jH1][1] += fH[1]; - ftmp[jH1][2] += fH[2]; - - ftmp[jH2][0] += fH[0]; - ftmp[jH2][1] += fH[1]; - ftmp[jH2][2] += fH[2]; + f[jH2][0] += fH[0]; + f[jH2][1] += fH[1]; + f[jH2][2] += fH[2]; + if (vflag) { delx1 = x[j][0] - x1[0]; dely1 = x[j][1] - x1[1]; delz1 = x[j][2] - x1[2]; @@ -380,12 +340,16 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) delz3 = x[jH2][2] - x1[2]; domain->minimum_image(delx3,dely3,delz3); - virialtmp[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); - virialtmp[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); - virialtmp[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); - virialtmp[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); - virialtmp[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); - virialtmp[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + v[0] += 0.5 * (delx1 * fO[0] + (delx2 + delx3) * fH[0]); + v[1] += 0.5 * (dely1 * fO[1] + (dely2 + dely3) * fH[1]); + v[2] += 0.5 * (delz1 * fO[2] + (delz2 + delz3) * fH[2]); + v[3] += 0.5 * (dely1 * fO[0] + (dely2 + dely3) * fH[0]); + v[4] += 0.5 * (delz1 * fO[0] + (delz2 + delz3) * fH[0]); + v[5] += 0.5 * (delz1 * fO[1] + (delz2 + delz3) * fH[1]); + + vlist[n++] = j; + vlist[n++] = jH1; + vlist[n++] = jH2; } } @@ -397,24 +361,13 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) ecoul = qtmp*q[j] * table; } if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; - } - } else ecoul = 0.0; + } else ecoul = 0.0; - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,ecoul,fpair,delx,dely,delz); + if (evflag) ev_tally_list(n,vlist,ecoul,v); + } } } } - - if (vflag_fdotr) { - virial_compute(); - for (int i = 0; i < 6; i++) virial[i] += virialtmp[i]; - for (int i = 0; i < nall; i++) { - f[i][0] += ftmp[i][0]; - f[i][1] += ftmp[i][1]; - f[i][2] += ftmp[i][2]; - } - } } /* ---------------------------------------------------------------------- @@ -459,68 +412,14 @@ void PairLJCutCoulLongTIP4P::init_style() error->all("Pair style lj/cut/coul/long/tip4p requires newton pair on"); if (!atom->q_flag) error->all("Pair style lj/cut/coul/long/tip4p requires atom attribute q"); - - // request regular or rRESPA neighbor lists - - int irequest; - - if (update->whichflag == 0 && strcmp(update->integrate_style,"respa") == 0) { - int respa = 0; - if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; - if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; - - if (respa == 0) irequest = neighbor->request(this); - else if (respa == 1) { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } else { - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 1; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respainner = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 2; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respamiddle = 1; - irequest = neighbor->request(this); - neighbor->requests[irequest]->id = 3; - neighbor->requests[irequest]->half = 0; - neighbor->requests[irequest]->respaouter = 1; - } - - } else irequest = neighbor->request(this); - - cut_coulsq = cut_coul * cut_coul; - - // set & error check interior rRESPA cutoffs - - if (strcmp(update->integrate_style,"respa") == 0) { - if (((Respa *) update->integrate)->level_inner >= 0) { - cut_respa = ((Respa *) update->integrate)->cutoff; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (MIN(cut_lj[i][j],cut_coul) < cut_respa[3]) - error->all("Pair cutoff < Respa interior cutoff"); - } - } else cut_respa = NULL; - - // insure use of correct KSpace long-range solver, set g_ewald - - if (force->kspace == NULL) + if (strcmp(force->kspace_style,"pppm/tip4p") != 0) error->all("Pair style is incompatible with KSpace style"); - if (strcmp(force->kspace_style,"pppm/tip4p") == 0) - g_ewald = force->kspace->g_ewald; - else error->all("Pair style is incompatible with KSpace style"); + if (force->bond == NULL) + error->all("Must use a bond style with TIP4P potential"); + if (force->angle == NULL) + error->all("Must use an angle style with TIP4P potential"); - // setup force tables - - if (ncoultablebits) init_tables(); + PairLJCutCoulLong::init_style(); // set alpha parameter diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 7cfc9662c4..3bb10935bd 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -90,11 +90,6 @@ void PairAIREBO::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style airebo does not yet support peratom energy/virial"); - REBO_neigh(); FREBO(eflag,vflag); if (ljflag) FLJ(eflag,vflag); @@ -507,7 +502,7 @@ void PairAIREBO::FREBO(int eflag, int vflag) del[0] = delx; del[1] = dely; del[2] = delz; - bij = bondorder(i,j,del,rij,VA,f); + bij = bondorder(i,j,del,rij,VA,f,vflag_atom); dVAdi = bij*dVA; fpair = -(dVRdi+dVAdi) / rij; @@ -537,15 +532,16 @@ void PairAIREBO::FLJ(int eflag, int vflag) int testpath,npath,done; double evdwl,fpair; double rsq,best,wik,wkm,cij,rij,dwij,dwik,dwkj,dwkm,dwmj; - double delij[3],rijsq,delik[3],rik,delkj[3]; + double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; - double delkm[3],rkm,delmj[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; + double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; - double delikS[3],delkjS[3],delkmS[3],delmjS[3]; + double delikS[3],deljkS[3],delkmS[3],deljmS[3],delimS[3]; double rikS,rkjS,rkmS,rmjS,wikS,dwikS; double wkjS,dwkjS,wkmS,dwkmS,wmjS,dwmjS; double fpair1,fpair2,fpair3; + double fi[3],fj[3],fk[3],fm[3]; // I-J interaction from full neighbor list // skip 1/2 of interactions since only consider each pair once @@ -629,10 +625,10 @@ void PairAIREBO::FLJ(int eflag, int vflag) } else wik = 0.0; if (wik > best) { - delkj[0] = x[k][0] - x[j][0]; - delkj[1] = x[k][1] - x[j][1]; - delkj[2] = x[k][2] - x[j][2]; - rsq = delkj[0]*delkj[0] + delkj[1]*delkj[1] + delkj[2]*delkj[2]; + deljk[0] = x[j][0] - x[k][0]; + deljk[1] = x[j][1] - x[k][1]; + deljk[2] = x[j][2] - x[k][2]; + rsq = deljk[0]*deljk[0] + deljk[1]*deljk[1] + deljk[2]*deljk[2]; if (rsq < rcmaxsq[ktype][jtype]) { rkj = sqrt(rsq); wkj = Sp(rkj,rcmin[ktype][jtype],rcmax[ktype][jtype],dwkj); @@ -646,9 +642,9 @@ void PairAIREBO::FLJ(int eflag, int vflag) rikS = rik; wikS = wik; dwikS = dwik; - delkjS[0] = delkj[0]; - delkjS[1] = delkj[1]; - delkjS[2] = delkj[2]; + deljkS[0] = deljk[0]; + deljkS[1] = deljk[1]; + deljkS[2] = deljk[2]; rkjS = rkj; wkjS = wkj; dwkjS = dwkj; @@ -679,11 +675,11 @@ void PairAIREBO::FLJ(int eflag, int vflag) } else wkm = 0.0; if (wik*wkm > best) { - delmj[0] = x[m][0] - x[j][0]; - delmj[1] = x[m][1] - x[j][1]; - delmj[2] = x[m][2] - x[j][2]; - rsq = delmj[0]*delmj[0] + delmj[1]*delmj[1] + - delmj[2]*delmj[2]; + deljm[0] = x[j][0] - x[m][0]; + deljm[1] = x[j][1] - x[m][1]; + deljm[2] = x[j][2] - x[m][2]; + rsq = deljm[0]*deljm[0] + deljm[1]*deljm[1] + + deljm[2]*deljm[2]; if (rsq < rcmaxsq[mtype][jtype]) { rmj = sqrt(rsq); wmj = Sp(rmj,rcmin[mtype][jtype],rcmax[mtype][jtype],dwmj); @@ -704,9 +700,9 @@ void PairAIREBO::FLJ(int eflag, int vflag) rkmS = rkm; wkmS = wkm; dwkmS = dwkm; - delmjS[0] = delmj[0]; - delmjS[1] = delmj[1]; - delmjS[2] = delmj[2]; + deljmS[0] = deljm[0]; + deljmS[1] = deljm[1]; + deljmS[2] = deljm[2]; rmjS = rmj; wmjS = wmj; dwmjS = dwmj; @@ -739,7 +735,8 @@ void PairAIREBO::FLJ(int eflag, int vflag) delscale[0] = scale * delij[0]; delscale[1] = scale * delij[1]; delscale[2] = scale * delij[2]; - Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,f); + Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA, + delij,rij,f,vflag_atom); } else Stb = 0.0; fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) + @@ -766,50 +763,71 @@ void PairAIREBO::FLJ(int eflag, int vflag) f[atomj][0] -= delij[0]*fpair; f[atomj][1] -= delij[1]*fpair; f[atomj][2] -= delij[2]*fpair; - if (vflag_atom) ev_tally(atomi,atomj,nlocal,newton_pair, - 0.0,0.0,fpair,delij[0],delij[1],delij[2]); + + if (vflag_atom) v_tally2(atomi,atomj,fpair,delij); + } else if (npath == 3) { fpair1 = dC*dwikS*wkjS / rikS; - f[atomi][0] += delikS[0]*fpair1; - f[atomi][1] += delikS[1]*fpair1; - f[atomi][2] += delikS[2]*fpair1; - f[atomk][0] -= delikS[0]*fpair1; - f[atomk][1] -= delikS[1]*fpair1; - f[atomk][2] -= delikS[2]*fpair1; + fi[0] = delikS[0]*fpair1; + fi[1] = delikS[1]*fpair1; + fi[2] = delikS[2]*fpair1; fpair2 = dC*wikS*dwkjS / rkjS; - f[atomk][0] += delkjS[0]*fpair2; - f[atomk][1] += delkjS[1]*fpair2; - f[atomk][2] += delkjS[2]*fpair2; - f[atomj][0] -= delkjS[0]*fpair2; - f[atomj][1] -= delkjS[1]*fpair2; - f[atomj][2] -= delkjS[2]*fpair2; + fj[0] = deljkS[0]*fpair2; + fj[1] = deljkS[1]*fpair2; + fj[2] = deljkS[2]*fpair2; + + f[atomi][0] += fi[0]; + f[atomi][1] += fi[1]; + f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; + f[atomj][1] += fj[1]; + f[atomj][2] += fj[2]; + f[atomk][0] -= fi[0] + fj[0]; + f[atomk][1] -= fi[1] + fj[1]; + f[atomk][2] -= fi[2] + fj[2]; + if (vflag_atom) - v_tally3(atomi,atomj,atomk,fpair1,fpair2,delikS,delkjS); + v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS); + } else { fpair1 = dC*dwikS*wkmS*wmjS / rikS; - f[atomi][0] += delikS[0]*fpair1; - f[atomi][1] += delikS[1]*fpair1; - f[atomi][2] += delikS[2]*fpair1; - f[atomk][0] -= delikS[0]*fpair1; - f[atomk][1] -= delikS[1]*fpair1; - f[atomk][2] -= delikS[2]*fpair1; + fi[0] = delikS[0]*fpair1; + fi[1] = delikS[1]*fpair1; + fi[2] = delikS[2]*fpair1; + fpair2 = dC*wikS*dwkmS*wmjS / rkmS; - f[atomk][0] += delkmS[0]*fpair2; - f[atomk][1] += delkmS[1]*fpair2; - f[atomk][2] += delkmS[2]*fpair2; - f[atomm][0] -= delkmS[0]*fpair2; - f[atomm][1] -= delkmS[1]*fpair2; - f[atomm][2] -= delkmS[2]*fpair2; + fk[0] = delkmS[0]*fpair2 - fi[0]; + fk[1] = delkmS[1]*fpair2 - fi[1]; + fk[2] = delkmS[2]*fpair2 - fi[2]; + fpair3 = dC*wikS*wkmS*dwmjS / rmjS; - f[atomm][0] += delmjS[0]*fpair3; - f[atomm][1] += delmjS[1]*fpair3; - f[atomm][2] += delmjS[2]*fpair3; - f[atomj][0] -= delmjS[0]*fpair3; - f[atomj][1] -= delmjS[1]*fpair3; - f[atomj][2] -= delmjS[2]*fpair3; - if (vflag_atom) - v_tally4(atomi,atomj,atomk,atomm, - fpair1,fpair2,fpair3,delikS,delkmS,delmjS); + fj[0] = deljmS[0]*fpair3; + fj[1] = deljmS[1]*fpair3; + fj[2] = deljmS[2]*fpair3; + + fm[0] = -delkmS[0]*fpair2 - fj[0]; + fm[1] = -delkmS[1]*fpair2 - fj[1]; + fm[2] = -delkmS[2]*fpair2 - fj[2]; + + f[atomi][0] += fi[0]; + f[atomi][1] += fi[1]; + f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; + f[atomj][1] += fj[1]; + f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; + f[atomk][1] += fk[1]; + f[atomk][2] += fk[2]; + f[atomm][0] += fm[0]; + f[atomm][1] += fm[1]; + f[atomm][2] += fm[2]; + + if (vflag_atom) { + delimS[0] = delikS[0] + delkmS[0]; + delimS[1] = delikS[1] + delkmS[1]; + delimS[2] = delikS[2] + delkmS[2]; + v_tally4(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS); + } } } } @@ -839,13 +857,11 @@ void PairAIREBO::TORSION(int eflag, int vflag) double dxidij,dxidik,dxidjk,dxjdji,dxjdjl,dxjdil; double ddndij,ddndik,ddndjk,ddndjl,ddndil,dcwddn,dcwdn,dvpdcw,Ftmp[3]; double del32[3],rsq,r32,del23[3],del21[3],r21; - double deljk[3],del34[3],delil[3],r23,r34; + double deljk[3],del34[3],delil[3],delkl[3],r23,r34; double fi[3],fj[3],fk[3],fl[3]; int itype,jtype,ktype,ltype,kk,ll,jj; int *REBO_neighs_i,*REBO_neighs_j; - evdwl = 0.0; - double **x = atom->x; double **f = atom->f; int *type = atom->type; @@ -1146,7 +1162,12 @@ void PairAIREBO::TORSION(int eflag, int vflag) f[k][0] += fk[0]; f[k][1] += fk[1]; f[k][2] += fk[2]; f[l][0] += fl[0]; f[l][1] += fl[1]; f[l][2] += fl[2]; - //if (evflag) ev_tally4(i,j,k,l,fi,fj,fk,fl); + if (evflag) { + delkl[0] = delil[0] - del21[0]; + delkl[1] = delil[1] - del21[1]; + delkl[2] = delil[2] - del21[2]; + ev_tally4(i,j,k,l,evdwl,fi,fj,fk,delil,del34,delkl); + } } } } @@ -1212,11 +1233,12 @@ double PairAIREBO::Sp2(double Xij, double Xmin, double Xmax, double &dX) ------------------------------------------------------------------------- */ double PairAIREBO::bondorder(int i, int j, double rij[3], - double rijmag, double VA, double **f) + double rijmag, double VA, + double **f, int vflag_atom) { int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int itype,jtype,ktype,ltype,ntype; - double rik[3], rjl[3], rkn[3],rknmag,dNki,dwjl,bij; + double rik[3],rjl[3],rkn[3],rji[3],rki[3],rlj[3],rknmag,dNki,dwjl,bij; double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl; double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3; double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ,Nki,Nlj,dS; @@ -1225,7 +1247,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], double dN2[2],dN3[3]; double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; double Tij; - double r32[3],r32mag,cos321; + double r32[3],r32mag,cos321,r43[3],r13[3]; double dNlj; double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; @@ -1236,7 +1258,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; - double F23[3],F12[3],F34[3],F31[3],F24[3]; + double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; + double f1[3],f2[3],f3[3],f4[4]; double dcut321,PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; @@ -1338,58 +1361,68 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - f[atomj][0] -= tmp2*dcosjikdrj[0]; - f[atomj][1] -= tmp2*dcosjikdrj[1]; - f[atomj][2] -= tmp2*dcosjikdrj[2]; - f[atomi][0] -= tmp2*dcosjikdri[0]; - f[atomi][1] -= tmp2*dcosjikdri[1]; - f[atomi][2] -= tmp2*dcosjikdri[2]; - f[atomk][0] -= tmp2*dcosjikdrk[0]; - f[atomk][1] -= tmp2*dcosjikdrk[1]; - f[atomk][2] -= tmp2*dcosjikdrk[2]; + 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]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); - f[atomj][0] -= tmp2*(-rij[0]/rijmag); - f[atomj][1] -= tmp2*(-rij[1]/rijmag); - f[atomj][2] -= tmp2*(-rij[2]/rijmag); - f[atomi][0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); - f[atomi][1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); - f[atomi][2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); - f[atomk][0] -= tmp2*(rik[0]/rikmag); - f[atomk][1] -= tmp2*(rik[1]/rikmag); - f[atomk][2] -= tmp2*(rik[2]/rikmag); + 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)); + fk[0] -= tmp2*(rik[0]/rikmag); + fk[1] -= tmp2*(rik[1]/rikmag); + fk[2] -= tmp2*(rik[2]/rikmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwik*g*exp(lamdajik))/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ktype]*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_atom) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); + } } } @@ -1470,58 +1503,67 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - f[atomi][0] -= tmp2*dcosijldri[0]; - f[atomi][1] -= tmp2*dcosijldri[1]; - f[atomi][2] -= tmp2*dcosijldri[2]; - f[atomj][0] -= tmp2*dcosijldrj[0]; - f[atomj][1] -= tmp2*dcosijldrj[1]; - f[atomj][2] -= tmp2*dcosijldrj[2]; - f[atoml][0] -= tmp2*dcosijldrl[0]; - f[atoml][1] -= tmp2*dcosijldrl[1]; - f[atoml][2] -= tmp2*dcosijldrl[2]; + 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]; tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); - f[atomi][0] -= tmp2*(rij[0]/rijmag); - f[atomi][1] -= tmp2*(rij[1]/rijmag); - f[atomi][2] -= tmp2*(rij[2]/rijmag); - f[atomj][0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); - f[atomj][1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); - f[atomj][2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); - f[atoml][0] -= tmp2*(rjl[0]/rjlmag); - f[atoml][1] -= tmp2*(rjl[1]/rjlmag); - f[atoml][2] -= tmp2*(rjl[2]/rjlmag); + 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)); + fl[0] -= tmp2*(rjl[0]/rjlmag); + fl[1] -= tmp2*(rjl[1]/rjlmag); + fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ltype]*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_atom) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); + } } } @@ -1554,6 +1596,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1562,6 +1606,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -1581,6 +1627,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -1611,6 +1659,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1619,6 +1669,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -1638,6 +1690,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -1809,38 +1863,53 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], F24[1] = (fcilpc*ril[1]); F24[2] = (fcilpc*ril[2]); - f[atom1][0] +=-F12[0]-F31[0]; - f[atom1][1] +=-F12[1]-F31[1]; - f[atom1][2] +=-F12[2]-F31[2]; - f[atom2][0] += F23[0]+F12[0]+F24[0]; - f[atom2][1] += F23[1]+F12[1]+F24[1]; - f[atom2][2] += F23[2]+F12[2]+F24[2]; - f[atom3][0] += -F23[0]+F34[0]+F31[0]; - f[atom3][1] += -F23[1]+F34[1]+F31[1]; - f[atom3][2] += -F23[2]+F34[2]+F31[2]; - f[atom4][0] += -F34[0]-F24[0]; - f[atom4][1] += -F34[1]-F24[1]; - f[atom4][2] += -F34[2]-F24[2]; + f1[0] = -F12[0]-F31[0]; + f1[1] = -F12[1]-F31[1]; + f1[2] = -F12[2]-F31[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; + f4[0] = -F34[0]-F24[0]; + f4[1] = -F34[1]-F24[1]; + f4[2] = -F34[2]-F24[2]; // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*dw21*w34/r21mag; - f[atom2][0] -= tmp2*r21[0]; - f[atom2][1] -= tmp2*r21[1]; - f[atom2][2] -= tmp2*r21[2]; - f[atom1][0] += tmp2*r21[0]; - f[atom1][1] += tmp2*r21[1]; - f[atom1][2] += tmp2*r21[2]; + f2[0] -= tmp2*r21[0]; + f2[1] -= tmp2*r21[1]; + f2[2] -= tmp2*r21[2]; + f1[0] += tmp2*r21[0]; + f1[1] += tmp2*r21[1]; + f1[2] += tmp2*r21[2]; tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*w21*dw34/r34mag; - f[atom3][0] -= tmp2*r34[0]; - f[atom3][1] -= tmp2*r34[1]; - f[atom3][2] -= tmp2*r34[2]; - f[atom4][0] += tmp2*r34[0]; - f[atom4][1] += tmp2*r34[1]; - f[atom4][2] += tmp2*r34[2]; + f3[0] -= tmp2*r34[0]; + f3[1] -= tmp2*r34[1]; + f3[2] -= tmp2*r34[2]; + f4[0] += tmp2*r34[0]; + f4[1] += tmp2*r34[1]; + f4[2] += tmp2*r34[2]; + + f[atom1][0] += f1[0]; f[atom1][1] += f1[1]; + f[atom1][2] += f1[2]; + f[atom2][0] += f2[0]; f[atom2][1] += f2[1]; + f[atom2][2] += f2[2]; + f[atom3][0] += f3[0]; f[atom3][1] += f3[1]; + f[atom3][2] += f3[2]; + f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; + f[atom4][2] += f4[2]; + + if (vflag_atom) { + r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; + r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; + v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); + } } } } @@ -1872,6 +1941,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -1880,6 +1951,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -1899,6 +1972,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -1929,6 +2004,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -1937,6 +2014,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -1956,6 +2035,8 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -1973,7 +2054,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, double VA, double rij0[3], double rij0mag, - double **f) + double **f, int vflag_atom) { int k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; int atomi,atomj,itype,jtype,ktype,ltype,ntype; @@ -1987,11 +2068,9 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, double dcosijldrj[3],dcosijldrl[3],dcosjikdrj[3],dwjl; double Tij,crosskij[3],crosskijmag; double crossijl[3],crossijlmag,omkijl; - double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; double bij,tmp3pij,tmp3pji,Stb,dStb; double r32[3],r32mag,cos321; - double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; double w21,dw21,r34[3],r34mag,cos234,w34,dw34; @@ -2006,6 +2085,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; double F12[3],F23[3],F34[3],F31[3],F24[3]; + double fi[3],fj[3],fk[3],fl[3],f1[3],f2[3],f3[3],f4[4]; + double rji[3],rki[3],rlj[3],r13[3],r43[3]; double **x = atom->x; int *type = atom->type; @@ -2251,58 +2332,68 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - f[atomj][0] -= tmp2*dcosjikdrj[0]; - f[atomj][1] -= tmp2*dcosjikdrj[1]; - f[atomj][2] -= tmp2*dcosjikdrj[2]; - f[atomi][0] -= tmp2*dcosjikdri[0]; - f[atomi][1] -= tmp2*dcosjikdri[1]; - f[atomi][2] -= tmp2*dcosjikdri[2]; - f[atomk][0] -= tmp2*dcosjikdrk[0]; - f[atomk][1] -= tmp2*dcosjikdrk[1]; - f[atomk][2] -= tmp2*dcosjikdrk[2]; + 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]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); - f[atomj][0] -= tmp2*(-rij[0]/rijmag); - f[atomj][1] -= tmp2*(-rij[1]/rijmag); - f[atomj][2] -= tmp2*(-rij[2]/rijmag); - f[atomi][0] -= tmp2*((-rik[0]/rikmag)+(rij[0]/rijmag)); - f[atomi][1] -= tmp2*((-rik[1]/rikmag)+(rij[1]/rijmag)); - f[atomi][2] -= tmp2*((-rik[2]/rikmag)+(rij[2]/rijmag)); - f[atomk][0] -= tmp2*(rik[0]/rikmag); - f[atomk][1] -= tmp2*(rik[1]/rikmag); - f[atomk][2] -= tmp2*(rik[2]/rikmag); + 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)); + fk[0] -= tmp2*(rik[0]/rikmag); + fk[1] -= tmp2*(rik[1]/rikmag); + fk[2] -= tmp2*(rik[2]/rikmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwik*g*exp(lamdajik))/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ktype]*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; // dgdN forces tmp2 = VA*.5*(tmp*tmp3*dwik)/rikmag; - f[atomi][0] -= tmp2*rik[0]; - f[atomi][1] -= tmp2*rik[1]; - f[atomi][2] -= tmp2*rik[2]; - f[atomk][0] += tmp2*rik[0]; - f[atomk][1] += tmp2*rik[1]; - f[atomk][2] += tmp2*rik[2]; + fi[0] -= tmp2*rik[0]; + fi[1] -= tmp2*rik[1]; + fi[2] -= tmp2*rik[2]; + fk[0] += tmp2*rik[0]; + fk[1] += tmp2*rik[1]; + fk[2] += tmp2*rik[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atomk][0] += fk[0]; f[atomk][1] += fk[1]; f[atomk][2] += fk[2]; + + if (vflag_atom) { + rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2]; + rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2]; + v_tally3(atomi,atomj,atomk,fj,fk,rji,rki); + } } } @@ -2350,57 +2441,66 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - f[atomi][0] -= tmp2*dcosijldri[0]; - f[atomi][1] -= tmp2*dcosijldri[1]; - f[atomi][2] -= tmp2*dcosijldri[2]; - f[atomj][0] -= tmp2*dcosijldrj[0]; - f[atomj][1] -= tmp2*dcosijldrj[1]; - f[atomj][2] -= tmp2*dcosijldrj[2]; - f[atoml][0] -= tmp2*dcosijldrl[0]; - f[atoml][1] -= tmp2*dcosijldrl[1]; - f[atoml][2] -= tmp2*dcosijldrl[2]; + 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]; tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); - f[atomi][0] -= tmp2*(rij[0]/rijmag); - f[atomi][1] -= tmp2*(rij[1]/rijmag); - f[atomi][2] -= tmp2*(rij[2]/rijmag); - f[atomj][0] -= tmp2*((-rjl[0]/rjlmag)-(rij[0]/rijmag)); - f[atomj][1] -= tmp2*((-rjl[1]/rjlmag)-(rij[1]/rijmag)); - f[atomj][2] -= tmp2*((-rjl[2]/rjlmag)-(rij[2]/rijmag)); - f[atoml][0] -= tmp2*(rjl[0]/rjlmag); - f[atoml][1] -= tmp2*(rjl[1]/rjlmag); - f[atoml][2] -= tmp2*(rjl[2]/rjlmag); + 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)); + fl[0] -= tmp2*(rjl[0]/rjlmag); + fl[1] -= tmp2*(rjl[1]/rjlmag); + fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // PIJ forces tmp2 = VA*.5*(tmp*dN2[ltype]*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; // dgdN forces tmp2=VA*.5*(tmp*tmp3*dwjl)/rjlmag; - f[atomj][0] -= tmp2*rjl[0]; - f[atomj][1] -= tmp2*rjl[1]; - f[atomj][2] -= tmp2*rjl[2]; - f[atoml][0] += tmp2*rjl[0]; - f[atoml][1] += tmp2*rjl[1]; - f[atoml][2] += tmp2*rjl[2]; + fj[0] -= tmp2*rjl[0]; + fj[1] -= tmp2*rjl[1]; + fj[2] -= tmp2*rjl[2]; + fl[0] += tmp2*rjl[0]; + fl[1] += tmp2*rjl[1]; + fl[2] += tmp2*rjl[2]; + + f[atomi][0] += fi[0]; f[atomi][1] += fi[1]; f[atomi][2] += fi[2]; + f[atomj][0] += fj[0]; f[atomj][1] += fj[1]; f[atomj][2] += fj[2]; + f[atoml][0] += fl[0]; f[atoml][1] += fl[1]; f[atoml][2] += fl[2]; + + if (vflag_atom) { + rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2]; + v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj); + } } } @@ -2432,6 +2532,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2440,6 +2542,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) > TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -2459,6 +2563,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2489,6 +2595,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2497,6 +2605,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -2516,6 +2626,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } @@ -2686,38 +2798,53 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, F24[1] = (fcilpc*ril[1]); F24[2] = (fcilpc*ril[2]); - f[atom1][0] +=-F12[0]-F31[0]; - f[atom1][1] +=-F12[1]-F31[1]; - f[atom1][2] +=-F12[2]-F31[2]; - f[atom2][0] += F23[0]+F12[0]+F24[0]; - f[atom2][1] += F23[1]+F12[1]+F24[1]; - f[atom2][2] += F23[2]+F12[2]+F24[2]; - f[atom3][0] += -F23[0]+F34[0]+F31[0]; - f[atom3][1] += -F23[1]+F34[1]+F31[1]; - f[atom3][2] += -F23[2]+F34[2]+F31[2]; - f[atom4][0] += -F34[0]-F24[0]; - f[atom4][1] += -F34[1]-F24[1]; - f[atom4][2] += -F34[2]-F24[2]; + f1[0] = -F12[0]-F31[0]; + f1[1] = -F12[1]-F31[1]; + f1[2] = -F12[2]-F31[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; + f4[0] = -F34[0]-F24[0]; + f4[1] = -F34[1]-F24[1]; + f4[2] = -F34[2]-F24[2]; // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*dw21*w34/r21mag; - f[atom2][0] -= tmp2*r21[0]; - f[atom2][1] -= tmp2*r21[1]; - f[atom2][2] -= tmp2*r21[2]; - f[atom1][0] += tmp2*r21[0]; - f[atom1][1] += tmp2*r21[1]; - f[atom1][2] += tmp2*r21[2]; + f2[0] -= tmp2*r21[0]; + f2[1] -= tmp2*r21[1]; + f2[2] -= tmp2*r21[2]; + f1[0] += tmp2*r21[0]; + f1[1] += tmp2*r21[1]; + f1[2] += tmp2*r21[2]; tmp2 = VA*Tij*((1.0-(om1234*om1234))) * (1.0-tspjik)*(1.0-tspijl)*w21*dw34/r34mag; - f[atom3][0] -= tmp2*r34[0]; - f[atom3][1] -= tmp2*r34[1]; - f[atom3][2] -= tmp2*r34[2]; - f[atom4][0] += tmp2*r34[0]; - f[atom4][1] += tmp2*r34[1]; - f[atom4][2] += tmp2*r34[2]; + f3[0] -= tmp2*r34[0]; + f3[1] -= tmp2*r34[1]; + f3[2] -= tmp2*r34[2]; + f4[0] += tmp2*r34[0]; + f4[1] += tmp2*r34[1]; + f4[2] += tmp2*r34[2]; + + f[atom1][0] += f1[0]; f[atom1][1] += f1[1]; + f[atom1][2] += f1[2]; + f[atom2][0] += f2[0]; f[atom2][1] += f2[1]; + f[atom2][2] += f2[2]; + f[atom3][0] += f3[0]; f[atom3][1] += f3[1]; + f[atom3][2] += f3[2]; + f[atom4][0] += f4[0]; f[atom4][1] += f4[1]; + f[atom4][2] += f4[2]; + + if (vflag_atom) { + r13[0] = -rjk[0]; r13[1] = -rjk[1]; r13[2] = -rjk[2]; + r43[0] = -r34[0]; r43[1] = -r34[1]; r43[2] = -r34[2]; + v_tally4(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43); + } } } } @@ -2747,6 +2874,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + tmp2 = VA*dN3[2]*(2.0*NconjtmpI*dwik*SpN)*Etmp/rikmag; f[atomi][0] -= tmp2*rik[0]; f[atomi][1] -= tmp2*rik[1]; @@ -2755,6 +2884,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomk][1] += tmp2*rik[1]; f[atomk][2] += tmp2*rik[2]; + if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik); + if (fabs(dNki) >TOL) { REBO_neighs_k = REBO_firstneigh[atomk]; for (n = 0; n < REBO_numneigh[atomk]; n++) { @@ -2774,6 +2905,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomn][0] += tmp2*rkn[0]; f[atomn][1] += tmp2*rkn[1]; f[atomn][2] += tmp2*rkn[2]; + + if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn); } } } @@ -2804,6 +2937,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + tmp2 = VA*dN3[2]*(2.0*NconjtmpJ*dwjl*SpN)*Etmp/rjlmag; f[atomj][0] -= tmp2*rjl[0]; f[atomj][1] -= tmp2*rjl[1]; @@ -2812,6 +2947,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atoml][1] += tmp2*rjl[1]; f[atoml][2] += tmp2*rjl[2]; + if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl); + if (fabs(dNlj) > TOL) { REBO_neighs_l = REBO_firstneigh[atoml]; for (n = 0; n < REBO_numneigh[atoml]; n++) { @@ -2831,6 +2968,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f[atomn][0] += tmp2*rln[0]; f[atomn][1] += tmp2*rln[1]; f[atomn][2] += tmp2*rln[2]; + + if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln); } } } diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 50d205fcaf..b56906ab92 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -82,9 +82,9 @@ class PairAIREBO : public Pair { void FLJ(int, int); void TORSION(int, int); - double bondorder(int, int, double *, double, double, double **); + double bondorder(int, int, double *, double, double, double **, int); double bondorderLJ(int, int, double *, double, double, - double *, double, double **); + double *, double, double **, int); double Sp(double, double, double, double &); double Sp2(double, double, double, double &); diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index d5b27f68ef..e92b5497de 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -226,7 +226,7 @@ void PairTersoff::compute(int eflag, int vflag) f[k][1] += fk[1]; f[k][2] += fk[2]; - if (evflag) ev_tally3(i,j,k,0.0,0.0,fj,fk,delr1,delr2); + if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2); } } } diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index b67e80737d..99d921b4c0 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -54,7 +54,6 @@ PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp) rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL; gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL; arho1 = arho2 = arho3 = arho3b = t_ave = NULL; - strssa = NULL; maxneigh = 0; scrfcn = dscrfcn = fcpair = NULL; @@ -96,8 +95,6 @@ PairMEAM::~PairMEAM() memory->destroy_2d_double_array(arho3b); memory->destroy_2d_double_array(t_ave); - memory->destroy_3d_double_array(strssa); - memory->sfree(scrfcn); memory->sfree(dscrfcn); memory->sfree(fcpair); @@ -127,11 +124,6 @@ void PairMEAM::compute(int eflag, int vflag) if (eflag || vflag) ev_setup(eflag,vflag); else evflag = vflag_fdotr = 0; - // error check (for now) - - if (eflag_atom || vflag_atom) - error->all("Pair style meam does not yet support peratom energy/virial"); - int newton_pair = force->newton_pair; // grow local arrays if necessary @@ -153,7 +145,6 @@ void PairMEAM::compute(int eflag, int vflag) memory->destroy_2d_double_array(arho3); memory->destroy_2d_double_array(arho3b); memory->destroy_2d_double_array(t_ave); - memory->destroy_3d_double_array(strssa); nmax = atom->nmax; @@ -173,7 +164,6 @@ void PairMEAM::compute(int eflag, int vflag) arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3"); arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b"); t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave"); - strssa = memory->create_3d_double_array(nmax,3,3,"pair:strssa"); } // neighbor list info @@ -238,7 +228,7 @@ void PairMEAM::compute(int eflag, int vflag) for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_dens_init_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], + meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], @@ -252,10 +242,10 @@ void PairMEAM::compute(int eflag, int vflag) offset += numneigh_half[i]; } - reverse_flag = 0; comm->reverse_comm_pair(this); - meam_dens_final_(&nlocal,&nmax,&eflag,&evdwl,&ntype,type,fmap, + meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &eng_vdwl,eatom,&ntype,type,fmap, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); @@ -269,16 +259,24 @@ void PairMEAM::compute(int eflag, int vflag) offset = 0; + // vptr is first value in vatom if it will be used by meam_force() + // else vatom may not exist, so pass dummy ptr + + double *vptr; + if (vflag_atom) vptr = &vatom[0][0]; + else vptr = &cutmax; + for (ii = 0; ii < inum_half; ii++) { i = ilist_half[ii]; ifort = i+1; - meam_force_(&ifort,&nmax,&eflag,&evdwl,&ntype,type,fmap,&x[0][0], + meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom, + &vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0], &numneigh_half[i],firstneigh_half[i], &numneigh_full[i],firstneigh_full[i], &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], - &t_ave[0][0],&f[0][0],&strssa[0][0][0],&errorflag); + &t_ave[0][0],&f[0][0],vptr,&errorflag); if (errorflag) { char str[128]; sprintf(str,"MEAM library error %d",errorflag); @@ -287,17 +285,11 @@ void PairMEAM::compute(int eflag, int vflag) offset += numneigh_half[i]; } - reverse_flag = 1; - comm->reverse_comm_pair(this); - // change neighbor list indices back to C indexing neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half); neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full); - // just sum global energy (for now) - - if (evflag) ev_tally(0,0,nlocal,newton_pair,evdwl,0.0,0.0,0.0,0.0,0.0); if (vflag_fdotr) virial_compute(); } @@ -811,45 +803,28 @@ int PairMEAM::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; - if (reverse_flag == 0) { - for (i = first; i < last; i++) { - buf[m++] = rho0[i]; - buf[m++] = arho2b[i]; - buf[m++] = arho1[i][0]; - buf[m++] = arho1[i][1]; - buf[m++] = arho1[i][2]; - buf[m++] = arho2[i][0]; - buf[m++] = arho2[i][1]; - buf[m++] = arho2[i][2]; - buf[m++] = arho2[i][3]; - buf[m++] = arho2[i][4]; - buf[m++] = arho2[i][5]; - for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; - buf[m++] = arho3b[i][0]; - buf[m++] = arho3b[i][1]; - buf[m++] = arho3b[i][2]; - buf[m++] = t_ave[i][0]; - buf[m++] = t_ave[i][1]; - buf[m++] = t_ave[i][2]; - } - size = 27; - - } else { - for (i = first; i < last; i++) { - buf[m++] = strssa[i][0][0]; - buf[m++] = strssa[i][0][1]; - buf[m++] = strssa[i][0][2]; - buf[m++] = strssa[i][1][0]; - buf[m++] = strssa[i][1][1]; - buf[m++] = strssa[i][1][2]; - buf[m++] = strssa[i][2][0]; - buf[m++] = strssa[i][2][1]; - buf[m++] = strssa[i][2][2]; - } - size = 9; + for (i = first; i < last; i++) { + buf[m++] = rho0[i]; + buf[m++] = arho2b[i]; + buf[m++] = arho1[i][0]; + buf[m++] = arho1[i][1]; + buf[m++] = arho1[i][2]; + buf[m++] = arho2[i][0]; + buf[m++] = arho2[i][1]; + buf[m++] = arho2[i][2]; + buf[m++] = arho2[i][3]; + buf[m++] = arho2[i][4]; + buf[m++] = arho2[i][5]; + for (k = 0; k < 10; k++) buf[m++] = arho3[i][k]; + buf[m++] = arho3b[i][0]; + buf[m++] = arho3b[i][1]; + buf[m++] = arho3b[i][2]; + buf[m++] = t_ave[i][0]; + buf[m++] = t_ave[i][1]; + buf[m++] = t_ave[i][2]; } - return size; + return m; } /* ---------------------------------------------------------------------- */ @@ -859,42 +834,26 @@ void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf) int i,j,k,m; m = 0; - if (reverse_flag == 0) { - for (i = 0; i < n; i++) { - j = list[i]; - rho0[j] += buf[m++]; - arho2b[j] += buf[m++]; - arho1[j][0] += buf[m++]; - arho1[j][1] += buf[m++]; - arho1[j][2] += buf[m++]; - arho2[j][0] += buf[m++]; - arho2[j][1] += buf[m++]; - arho2[j][2] += buf[m++]; - arho2[j][3] += buf[m++]; - arho2[j][4] += buf[m++]; - arho2[j][5] += buf[m++]; - for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; - arho3b[j][0] += buf[m++]; - arho3b[j][1] += buf[m++]; - arho3b[j][2] += buf[m++]; - t_ave[j][0] += buf[m++]; - t_ave[j][1] += buf[m++]; - t_ave[j][2] += buf[m++]; - } - - } else { - for (i = 0; i < n; i++) { - j = list[i]; - strssa[j][0][0] += buf[m++]; - strssa[j][0][1] += buf[m++]; - strssa[j][0][2] += buf[m++]; - strssa[j][1][0] += buf[m++]; - strssa[j][1][1] += buf[m++]; - strssa[j][1][2] += buf[m++]; - strssa[j][2][0] += buf[m++]; - strssa[j][2][1] += buf[m++]; - strssa[j][2][2] += buf[m++]; - } + for (i = 0; i < n; i++) { + j = list[i]; + rho0[j] += buf[m++]; + arho2b[j] += buf[m++]; + arho1[j][0] += buf[m++]; + arho1[j][1] += buf[m++]; + arho1[j][2] += buf[m++]; + arho2[j][0] += buf[m++]; + arho2[j][1] += buf[m++]; + arho2[j][2] += buf[m++]; + arho2[j][3] += buf[m++]; + arho2[j][4] += buf[m++]; + arho2[j][5] += buf[m++]; + for (k = 0; k < 10; k++) arho3[j][k] += buf[m++]; + arho3b[j][0] += buf[m++]; + arho3b[j][1] += buf[m++]; + arho3b[j][2] += buf[m++]; + t_ave[j][0] += buf[m++]; + t_ave[j][1] += buf[m++]; + t_ave[j][2] += buf[m++]; } } @@ -906,7 +865,6 @@ double PairMEAM::memory_usage() { double bytes = 11 * nmax * sizeof(double); bytes += (3 + 6 + 10 + 3 + 3) * nmax * sizeof(double); - bytes += 3*3 * nmax * sizeof(double); bytes += 3 * maxneigh * sizeof(double); return bytes; } diff --git a/src/MEAM/pair_meam.h b/src/MEAM/pair_meam.h index 0d8f552bd5..329e3d293e 100644 --- a/src/MEAM/pair_meam.h +++ b/src/MEAM/pair_meam.h @@ -22,20 +22,22 @@ extern "C" { void meam_setup_param_(int *, double *, int *, int *, int *); void meam_setup_done_(double *); - void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, + void meam_dens_init_(int *, int *, int *, int *, int *, double *, int *, int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, int *); - void meam_dens_final_(int *, int *, int *, double *, int *, int *, int *, + void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, + int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, int *); - void meam_force_(int *, int *, int *, double *, int *, int *, int *, + void meam_force_(int *, int *, int *, int *, int *, int *, + double *, double *, int *, int *, int *, double *, int *, int *, int *, int *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, double *, @@ -71,7 +73,6 @@ class PairMEAM : public Pair { int nelements; // # of unique elements char **elements; // names of unique elements double *mass; // mass of each element - int reverse_flag; // which pass of reverse comm is being done int *map; // mapping from atom types to elements int *fmap; // Fortran version of map array for MEAM lib @@ -83,7 +84,6 @@ class PairMEAM : public Pair { double *rho,*rho0,*rho1,*rho2,*rho3,*frhop; double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b; double **arho1,**arho2,**arho3,**arho3b,**t_ave; - double ***strssa; void allocate(); void read_files(char *, char *); diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index a9665098fd..25d2575fc1 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -57,6 +57,7 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) : rigid_flag = 1; virial_flag = 1; + MPI_Comm_rank(world,&me); // perform initial allocation of atom-based arrays @@ -254,10 +255,6 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) : fprintf(logfile,"%d clusters, %d bodies, %d joints, %d atoms\n", ncluster,nbody,njoint,nsum); } - - // zero fix_poems virial in case pressure uses it before 1st fix_poems call - - for (int n = 0; n < 6; n++) virial[n] = 0.0; } /* ---------------------------------------------------------------------- @@ -346,14 +343,6 @@ void FixPOEMS::init() error->all("POEMS fix must come before NPT/NPH fix"); } - // compute poems contribution to virial every step if fix NPT,NPH exists - - pressure_flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1; - if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1; - } - // timestep info dtv = update->dt; @@ -597,9 +586,9 @@ void FixPOEMS::init() make setup call to POEMS ------------------------------------------------------------------------- */ -void FixPOEMS::setup() +void FixPOEMS::setup(int vflag) { - int i,j,ibody; + int i,j,n,ibody; // vcm = velocity of center-of-mass of each rigid body // angmom = angular momentum of each rigid body @@ -654,29 +643,31 @@ void FixPOEMS::setup() angmom[ibody][2] = all[ibody][5]; } - // set omega from angmom + // virial setup before call to set_v + + if (vflag) v_setup(vflag); + else evflag = 0; + + // set velocities from angmom & omega for (ibody = 0; ibody < nbody; ibody++) omega_from_mq(angmom[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody],inertia[ibody],omega[ibody]); + set_v(); - // reset velocities from omega // guestimate virial as 2x the set_v contribution - // since post_force doesn't compute virial - for (int n = 0; n < 6; n++) virial[n] = 0.0; - set_v(1); - for (int n = 0; n < 6; n++) virial[n] *= 2.0; - - double ke = 0.0; - for (i = 0; i < nlocal; i++) - if (natom2body[i]) - ke += mass[atom->type[i]] * - (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + if (vflag_global) + for (n = 0; n < 6; n++) virial[n] *= 2.0; + if (vflag_atom) { + for (i = 0; i < nlocal; i++) + for (n = 0; n < 6; n++) + vatom[i][n] *= 2.0; + } // use post_force() to compute initial fcm & torque - post_force(1); + post_force(vflag); // setup for POEMS @@ -691,17 +682,20 @@ void FixPOEMS::setup() set x,v of body atoms accordingly /* ---------------------------------------------------------------------- */ -void FixPOEMS::initial_integrate() +void FixPOEMS::initial_integrate(int vflag) { // perform POEMS integration poems->LobattoOne(xcm,vcm,omega,torque,fcm,ex_space,ey_space,ez_space); - // set coords and velocities of atoms in rigid bodies + // virial setup before call to set_xv - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_xv(vflag); + if (vflag) v_setup(vflag); + else evflag = 0; + + // set coords and velocities of atoms in rigid bodies + + set_xv(); } /* ---------------------------------------------------------------------- @@ -712,6 +706,8 @@ void FixPOEMS::initial_integrate() void FixPOEMS::post_force(int vflag) { int i,j,ibody; + int xbox,ybox,zbox; + double dx,dy,dz; int *image = atom->image; double **x = atom->x; @@ -722,8 +718,6 @@ void FixPOEMS::post_force(int vflag) double yprd = domain->yprd; double zprd = domain->zprd; - int xbox,ybox,zbox; - double dx,dy,dz; for (ibody = 0; ibody < nbody; ibody++) for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; @@ -772,15 +766,14 @@ void FixPOEMS::final_integrate() poems->LobattoTwo(vcm,omega,torque,fcm); // set velocities of atoms in rigid bodies + // virial is already setup from initial_integrate - int vflag = 0; - if (pressure_flag || output->next_thermo == update->ntimestep) vflag = 1; - set_v(vflag); + set_v(); } /* ---------------------------------------------------------------------- */ -void FixPOEMS::initial_integrate_respa(int ilevel, int flag) +void FixPOEMS::initial_integrate_respa(int vflag, int ilevel, int flag) { if (flag) return; // only used by NPT,NPH @@ -788,7 +781,7 @@ void FixPOEMS::initial_integrate_respa(int ilevel, int flag) dtf = 0.5 * step_respa[ilevel] * force->ftm2v; dthalf = 0.5 * step_respa[ilevel]; - if (ilevel == 0) initial_integrate(); + if (ilevel == 0) initial_integrate(vflag); else final_integrate(); } @@ -1344,30 +1337,26 @@ void FixPOEMS::omega_from_mq(double *m, double *ex, double *ey, double *ez, v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixPOEMS::set_xv(int vflag) +void FixPOEMS::set_xv() { + int ibody; + int xbox,ybox,zbox; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; + int *image = atom->image; double **x = atom->x; double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; - int ibody; - int xbox,ybox,zbox; - - double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2; - double *mass = atom->mass; - double **f = atom->f; - int *type = atom->type; - - // zero out fix_poems virial - - if (vflag) for (int n = 0; n < 6; n++) virial[n] = 0.0; - - // set x and v + // set x and v of each atom // only set joint atoms for 1st rigid body they belong to for (int i = 0; i < nlocal; i++) { @@ -1378,18 +1367,21 @@ void FixPOEMS::set_xv(int vflag) ybox = (image[i] >> 10 & 1023) - 512; zbox = (image[i] >> 20) - 512; - // save old positions and velocities for virial contribution + // save old positions and velocities for virial - if (vflag) { + if (evflag) { x0 = x[i][0] + xbox*xprd; x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } + // x = displacement from center-of-mass, based on body orientation + // v = vcm + omega around center-of-mass + x[i][0] = ex_space[ibody][0]*displace[i][0] + ey_space[ibody][0]*displace[i][1] + ez_space[ibody][0]*displace[i][2]; @@ -1407,24 +1399,33 @@ void FixPOEMS::set_xv(int vflag) v[i][2] = omega[ibody][0]*x[i][1] - omega[ibody][1]*x[i][0] + vcm[ibody][2]; + // add center of mass to displacement + // map back into periodic box via xbox,ybox,zbox + x[i][0] += xcm[ibody][0] - xbox*xprd; x[i][1] += xcm[ibody][1] - ybox*yprd; x[i][2] += xcm[ibody][2] - zbox*zprd; - // compute body constraint forces for virial + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c final_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; - virial[0] += 0.5*fc0*x0; - virial[1] += 0.5*fc1*x1; - virial[2] += 0.5*fc2*x2; - virial[3] += 0.5*fc1*x0; - virial[4] += 0.5*fc2*x0; - virial[5] += 0.5*fc2*x1; + vr[0] = 0.5*fc0*x0; + vr[1] = 0.5*fc1*x1; + vr[2] = 0.5*fc2*x2; + vr[3] = 0.5*fc1*x0; + vr[4] = 0.5*fc2*x0; + vr[5] = 0.5*fc2*x1; + + v_tally(1,&i,1.0,vr); } } } @@ -1434,26 +1435,27 @@ void FixPOEMS::set_xv(int vflag) v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ -void FixPOEMS::set_v(int vflag) +void FixPOEMS::set_v() { - double **v = atom->v; - int nlocal = atom->nlocal; - int ibody; + int xbox,ybox,zbox; double dx,dy,dz; + double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; + double vr[6]; - double vold0,vold1,vold2,fc0,fc1,fc2,massone,x0,x1,x2; double *mass = atom->mass; double **f = atom->f; double **x = atom->x; + double **v = atom->v; int *type = atom->type; int *image = atom->image; + int nlocal = atom->nlocal; + double xprd = domain->xprd; double yprd = domain->yprd; double zprd = domain->zprd; - int xbox,ybox,zbox; - // set v + // set v of each atom // only set joint atoms for 1st rigid body they belong to for (int i = 0; i < nlocal; i++) { @@ -1472,24 +1474,27 @@ void FixPOEMS::set_v(int vflag) // save old velocities for virial - if (vflag) { - vold0 = v[i][0]; - vold1 = v[i][1]; - vold2 = v[i][2]; + if (evflag) { + v0 = v[i][0]; + v1 = v[i][1]; + v2 = v[i][2]; } v[i][0] = omega[ibody][1]*dz - omega[ibody][2]*dy + vcm[ibody][0]; v[i][1] = omega[ibody][2]*dx - omega[ibody][0]*dz + vcm[ibody][1]; v[i][2] = omega[ibody][0]*dy - omega[ibody][1]*dx + vcm[ibody][2]; - // compute body constraint forces for virial - // use unwrapped atom positions + // virial = unwrapped coords dotted into body constraint force + // body constraint force = implied force due to v change minus f external + // assume f does not include forces internal to body + // 1/2 factor b/c initial_integrate contributes other half + // assume per-atom contribution is due to constraint force on that atom - if (vflag) { + if (evflag) { massone = mass[type[i]]; - fc0 = massone*(v[i][0] - vold0)/dtf - f[i][0]; - fc1 = massone*(v[i][1] - vold1)/dtf - f[i][1]; - fc2 = massone*(v[i][2] - vold2)/dtf - f[i][2]; + fc0 = massone*(v[i][0] - v0)/dtf - f[i][0]; + fc1 = massone*(v[i][1] - v1)/dtf - f[i][1]; + fc2 = massone*(v[i][2] - v2)/dtf - f[i][2]; xbox = (image[i] & 1023) - 512; ybox = (image[i] >> 10 & 1023) - 512; @@ -1499,12 +1504,14 @@ void FixPOEMS::set_v(int vflag) x1 = x[i][1] + ybox*yprd; x2 = x[i][2] + zbox*zprd; - virial[0] += 0.5*fc0*x0; - virial[1] += 0.5*fc1*x1; - virial[2] += 0.5*fc2*x2; - virial[3] += 0.5*fc1*x0; - virial[4] += 0.5*fc2*x0; - virial[5] += 0.5*fc2*x1; + vr[0] = 0.5*fc0*x0; + vr[1] = 0.5*fc1*x1; + vr[2] = 0.5*fc2*x2; + vr[3] = 0.5*fc1*x0; + vr[4] = 0.5*fc2*x0; + vr[5] = 0.5*fc2*x1; + + v_tally(1,&i,1.0,vr); } } } diff --git a/src/POEMS/fix_poems.h b/src/POEMS/fix_poems.h index 16534d0bef..c4c52da4f5 100644 --- a/src/POEMS/fix_poems.h +++ b/src/POEMS/fix_poems.h @@ -24,11 +24,11 @@ class FixPOEMS : public Fix { ~FixPOEMS(); int setmask(); void init(); - void setup(); - void initial_integrate(); + void setup(int); + void initial_integrate(int); void post_force(int); void final_integrate(); - void initial_integrate_respa(int, int); + void initial_integrate_respa(int, int, int); void post_force_respa(int, int, int); void final_integrate_respa(int); @@ -48,7 +48,6 @@ class FixPOEMS : public Fix { double dtv,dtf,dthalf; double *step_respa; int nlevels_respa; - int pressure_flag; double total_ke; // atom assignment to rigid bodies @@ -100,8 +99,8 @@ class FixPOEMS : public Fix { void rotate(double **, int, int, int, int, double, double); void omega_from_mq(double *, double *, double *, double *, double *, double *); - void set_v(int); - void set_xv(int); + void set_v(); + void set_xv(); }; } diff --git a/src/compute_stress_atom.cpp b/src/compute_stress_atom.cpp index 0e472ab90f..3f1703704b 100644 --- a/src/compute_stress_atom.cpp +++ b/src/compute_stress_atom.cpp @@ -117,7 +117,7 @@ void ComputeStressAtom::compute_peratom() stress[i][j] = 0.0; // add in per-atom contributions from each force - + if (pairflag && force->pair) { double **vatom = force->pair->vatom; for (i = 0; i < npair; i++) diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp index 49ad4a0a72..d443530b8d 100644 --- a/src/fix_rigid.cpp +++ b/src/fix_rigid.cpp @@ -296,14 +296,6 @@ void FixRigid::init() error->all("Rigid fix must come before NPT/NPH fix"); } - // compute rigid contribution to virial every step if fix NPT,NPH exists - - pressure_flag = 0; - for (int i = 0; i < modify->nfix; i++) { - if (strcmp(modify->fix[i]->style,"npt") == 0) pressure_flag = 1; - if (strcmp(modify->fix[i]->style,"nph") == 0) pressure_flag = 1; - } - // timestep info dtv = update->dt; @@ -720,7 +712,7 @@ void FixRigid::initial_integrate(int vflag) if (vflag) v_setup(vflag); else evflag = 0; - // set coords and velocities if atoms in rigid bodies + // set coords and velocities of atoms in rigid bodies // from quarternion and omega set_xv(); @@ -1240,6 +1232,9 @@ void FixRigid::set_xv() int *image = atom->image; double **x = atom->x; double **v = atom->v; + double **f = atom->f; + double *mass = atom->mass; + int *type = atom->type; int nlocal = atom->nlocal; double xprd = domain->xprd; @@ -1251,10 +1246,8 @@ void FixRigid::set_xv() xz = domain->xz; yz = domain->yz; } - - double *mass = atom->mass; - double **f = atom->f; - int *type = atom->type; + + // set x and v of each atom for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; @@ -1340,7 +1333,7 @@ void FixRigid::set_xv() } /* ---------------------------------------------------------------------- - set space-frame velocity of each atom in rigid body + set space-frame velocity of each atom in a rigid body v = Vcm + (W cross (x - Xcm)) ------------------------------------------------------------------------- */ @@ -1348,7 +1341,7 @@ void FixRigid::set_v() { int ibody; int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; + double dx,dy,dz; double x0,x1,x2,v0,v1,v2,fc0,fc1,fc2,massone; double xy,xz,yz; double vr[6]; @@ -1370,6 +1363,8 @@ void FixRigid::set_v() yz = domain->yz; } + // set v of each atom + for (int i = 0; i < nlocal; i++) { if (body[i] < 0) continue; ibody = body[i]; diff --git a/src/fix_rigid.h b/src/fix_rigid.h index 46d26a4fb5..e4dff3f654 100644 --- a/src/fix_rigid.h +++ b/src/fix_rigid.h @@ -44,7 +44,6 @@ class FixRigid : public Fix { private: double dtv,dtf,dtq; double *step_respa; - int pressure_flag; int triclinic; int nbody; // # of rigid bodies diff --git a/src/lammps.cpp b/src/lammps.cpp index ad03cefc87..6860fe14b7 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -44,6 +44,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) memory = new Memory(this); error = new Error(this); universe = new Universe(this,communicator); + output = NULL; // parse input switches diff --git a/src/pair.cpp b/src/pair.cpp index 733a886562..495fdae821 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -485,13 +485,12 @@ void Pair::ev_tally_xyz(int i, int j, int nlocal, int newton_pair, /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators - called by SW and Tersoff potentials, newton_pair is always on - virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drij*fj + drik*fk - could just pass ri,rj,rk since coords are already unwrapped by PBC + called by SW potential, newton_pair is always on + virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, - double *fj, double *fk, double *drij, double *drik) + double *fj, double *fk, double *drji, double *drki) { double epairthird,v[6]; @@ -509,12 +508,12 @@ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, } if (vflag_atom) { - v[0] = THIRD * (drij[0]*fj[0] + drik[0]*fk[0]); - v[1] = THIRD * (drij[1]*fj[1] + drik[1]*fk[1]); - v[2] = THIRD * (drij[2]*fj[2] + drik[2]*fk[2]); - v[3] = THIRD * (drij[0]*fj[1] + drik[0]*fk[1]); - v[4] = THIRD * (drij[0]*fj[2] + drik[0]*fk[2]); - v[5] = THIRD * (drij[1]*fj[2] + drik[1]*fk[2]); + v[0] = THIRD * (drji[0]*fj[0] + drki[0]*fk[0]); + v[1] = THIRD * (drji[1]*fj[1] + drki[1]*fk[1]); + v[2] = THIRD * (drji[2]*fj[2] + drki[2]*fk[2]); + v[3] = THIRD * (drji[0]*fj[1] + drki[0]*fk[1]); + v[4] = THIRD * (drji[0]*fj[2] + drki[0]*fk[2]); + v[5] = THIRD * (drji[1]*fj[2] + drki[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; @@ -525,23 +524,134 @@ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, } } +/* ---------------------------------------------------------------------- + tally eng_vdwl and virial into global and per-atom accumulators + called by AIREBO potential, newton_pair is always on + ------------------------------------------------------------------------- */ + +void Pair::ev_tally4(int i, int j, int k, int m, double evdwl, + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm) +{ + double epairfourth,v[6]; + + if (eflag_either) { + if (eflag_global) eng_vdwl += evdwl; + if (eflag_atom) { + epairfourth = 0.25 * evdwl; + eatom[i] += epairfourth; + eatom[j] += epairfourth; + eatom[k] += epairfourth; + eatom[m] += epairfourth; + } + } + + if (vflag_atom) { + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); + + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; + vatom[k][0] += v[0]; vatom[k][1] += v[1]; vatom[k][2] += v[2]; + vatom[k][3] += v[3]; vatom[k][4] += v[4]; vatom[k][5] += v[5]; + vatom[m][0] += v[0]; vatom[m][1] += v[1]; vatom[m][2] += v[2]; + vatom[m][3] += v[3]; vatom[m][4] += v[4]; vatom[m][5] += v[5]; + } +} + +/* ---------------------------------------------------------------------- + tally ecoul and virial into each of n atoms in list + called by TIP4P potential, newton_pair is always on + changes v values by dividing by n + ------------------------------------------------------------------------- */ + +void Pair::ev_tally_list(int n, int *list, double ecoul, double *v) +{ + int i,j; + + if (eflag_either) { + if (eflag_global) eng_coul += ecoul; + if (eflag_atom) { + double epairatom = ecoul/n; + for (i = 0; i < n; i++) eatom[list[i]] += epairatom; + } + } + + if (vflag_either) { + if (vflag_global) { + virial[0] += v[0]; + virial[1] += v[1]; + virial[2] += v[2]; + virial[3] += v[3]; + virial[4] += v[4]; + virial[5] += v[5]; + } + + if (vflag_atom) { + v[0] /= n; + v[1] /= n; + v[2] /= n; + v[3] /= n; + v[4] /= n; + v[5] /= n; + for (i = 0; i < n; i++) { + j = list[i]; + vatom[j][0] += v[0]; + vatom[j][1] += v[1]; + vatom[j][2] += v[2]; + vatom[j][3] += v[3]; + vatom[j][4] += v[4]; + vatom[j][5] += v[5]; + } + } + } +} + /* ---------------------------------------------------------------------- tally virial into per-atom accumulators - called by airebo potential, newton_pair is always on - could just pass ri,rj,rk since coords are already unwrapped by PBC + called by AIREBO potential, newton_pair is always on + fpair is magnitude of force on atom I ------------------------------------------------------------------------- */ -void Pair::v_tally3(int i, int j, int k, - double f1, double f2, double *drik, double *drkj) +void Pair::v_tally2(int i, int j, double fpair, double *drij) { double v[6]; - v[0] = 0.0; - v[1] = 0.0; - v[2] = 0.0; - v[3] = 0.0; - v[4] = 0.0; - v[5] = 0.0; + v[0] = 0.5 * drij[0]*drij[0]*fpair; + v[1] = 0.5 * drij[1]*drij[1]*fpair; + v[2] = 0.5 * drij[2]*drij[2]*fpair; + v[3] = 0.5 * drij[0]*drij[1]*fpair; + v[4] = 0.5 * drij[0]*drij[2]*fpair; + v[5] = 0.5 * drij[1]*drij[2]*fpair; + + vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; + vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; + vatom[j][0] += v[0]; vatom[j][1] += v[1]; vatom[j][2] += v[2]; + vatom[j][3] += v[3]; vatom[j][4] += v[4]; vatom[j][5] += v[5]; +} + +/* ---------------------------------------------------------------------- + tally virial into per-atom accumulators + called by AIREBO and Tersoff potential, newton_pair is always on +------------------------------------------------------------------------- */ + +void Pair::v_tally3(int i, int j, int k, + double *fi, double *fj, double *drik, double *drjk) +{ + double v[6]; + + v[0] = THIRD * (drik[0]*fi[0] + drjk[0]*fj[0]); + v[1] = THIRD * (drik[1]*fi[1] + drjk[1]*fj[1]); + v[2] = THIRD * (drik[2]*fi[2] + drjk[2]*fj[2]); + v[3] = THIRD * (drik[0]*fi[1] + drjk[0]*fj[1]); + v[4] = THIRD * (drik[0]*fi[2] + drjk[0]*fj[2]); + v[5] = THIRD * (drik[1]*fi[2] + drjk[1]*fj[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; @@ -553,22 +663,21 @@ void Pair::v_tally3(int i, int j, int k, /* ---------------------------------------------------------------------- tally virial into per-atom accumulators - called by airebo potential, newton_pair is always on - could just pass ri,rj,rk,rm since coords are already unwrapped by PBC + called by AIREBO potential, newton_pair is always on ------------------------------------------------------------------------- */ void Pair::v_tally4(int i, int j, int k, int m, - double f1, double f2, double f3, - double *drik, double *drkm, double *drmj) + double *fi, double *fj, double *fk, + double *drim, double *drjm, double *drkm) { double v[6]; - v[0] = 0.0; - v[1] = 0.0; - v[2] = 0.0; - v[3] = 0.0; - v[4] = 0.0; - v[5] = 0.0; + v[0] = 0.25 * (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]); + v[1] = 0.25 * (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]); + v[2] = 0.25 * (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]); + v[3] = 0.25 * (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]); + v[4] = 0.25 * (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]); + v[5] = 0.25 * (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]); vatom[i][0] += v[0]; vatom[i][1] += v[1]; vatom[i][2] += v[2]; vatom[i][3] += v[3]; vatom[i][4] += v[4]; vatom[i][5] += v[5]; diff --git a/src/pair.h b/src/pair.h index 7a7c47efca..b50d685c9d 100644 --- a/src/pair.h +++ b/src/pair.h @@ -119,8 +119,12 @@ class Pair : protected Pointers { double, double, double, double, double, double); void ev_tally3(int, int, int, double, double, double *, double *, double *, double *); - void v_tally3(int, int, int, double, double, double *, double *); - void v_tally4(int, int, int, int, double, double, double, + void ev_tally4(int, int, int, int, double, + double *, double *, double *, double *, double *, double *); + void ev_tally_list(int, int *, double, double *); + void v_tally2(int, int, double, double *); + void v_tally3(int, int, int, double *, double *, double *, double *); + void v_tally4(int, int, int, int, double *, double *, double *, double *, double *, double *); void virial_compute(); }; diff --git a/src/style_user_ackland.h b/src/style_user_ackland.h index 6e7483a9f7..e69de29bb2 100644 --- a/src/style_user_ackland.h +++ b/src/style_user_ackland.h @@ -1,20 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef ComputeInclude -#include "compute_ackland_atom.h" -#endif - -#ifdef ComputeClass -ComputeStyle(ackland/atom,ComputeAcklandAtom) -#endif diff --git a/src/style_user_ewaldn.h b/src/style_user_ewaldn.h index 3eafa50744..e69de29bb2 100644 --- a/src/style_user_ewaldn.h +++ b/src/style_user_ewaldn.h @@ -1,30 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef KSpaceInclude -#include "ewald_n.h" -#endif - -#ifdef KSpaceClass -KSpaceStyle(ewald/n,EwaldN) -#endif - -#ifdef PairInclude -#include "pair_buck_coul.h" -#include "pair_lj_coul.h" -#endif - -#ifdef PairClass -PairStyle(buck/coul,PairBuckCoul) -PairStyle(lj/coul,PairLJCoul) -#endif diff --git a/src/thermo.cpp b/src/thermo.cpp index aeb684821a..f16c63696a 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -256,7 +256,7 @@ void Thermo::init() int ivariable; for (i = 0; i < nvariable; i++) { ivariable = input->variable->find(id_variable[i]); - if (ivariable < 0) error->all("Could not find thermo variable name"); + if (ivariable < 0) error->all("Could not find thermo custom variable name"); variables[i] = ivariable; } diff --git a/src/variable.cpp b/src/variable.cpp index 8298fecf66..f8a0b217ee 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -1094,7 +1094,7 @@ double Variable::eval_tree(Tree *tree, int i) return exp(eval_tree(tree->left,i)); if (tree->type == LN) { double arg = eval_tree(tree->left,i); - if (arg <= 0.0) error->all("Log of negative/zero in variable formula"); + if (arg <= 0.0) error->all("Log of zero/negative in variable formula"); return log(arg); }