ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

make airebo compatible with pair_modify nofdotr

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-07-11 08:20:59 -04:00
parent a4748b4c28
commit 50da38722a
8 changed files with 228 additions and 175 deletions
Download Patch File Download Diff File Expand all files Collapse all files

98
src/MANYBODY/pair_airebo.cpp
View File

@ -85,10 +85,10 @@ PairAIREBO::~PairAIREBO()
{
memory->destroy(REBO_numneigh);
memory->sfree(REBO_firstneigh);
delete [] ipage;
delete[] ipage;
memory->destroy(nC);
memory->destroy(nH);
delete [] pvector;
delete[] pvector;
if (allocated) {
memory->destroy(setflag);
@ -100,7 +100,7 @@ PairAIREBO::~PairAIREBO()
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(lj4);
delete [] map;
delete[] map;
}
}
@ -112,9 +112,9 @@ void PairAIREBO::compute(int eflag, int vflag)
pvector[0] = pvector[1] = pvector[2] = 0.0;
REBO_neigh();
FREBO(eflag,vflag);
if (ljflag) FLJ(eflag,vflag);
if (torflag) TORSION(eflag,vflag);
FREBO(eflag);
if (ljflag) FLJ(eflag);
if (torflag) TORSION(eflag);
if (vflag_fdotr) virial_fdotr_compute();
}
@ -252,7 +252,7 @@ void PairAIREBO::init_style()
if (oneatom != neighbor->oneatom) create = 1;
if (create) {
delete [] ipage;
delete[] ipage;
pgsize = neighbor->pgsize;
oneatom = neighbor->oneatom;
@ -421,7 +421,7 @@ void PairAIREBO::REBO_neigh()
REBO forces and energy
------------------------------------------------------------------------- */
void PairAIREBO::FREBO(int eflag, int /*vflag*/)
void PairAIREBO::FREBO(int eflag)
{
int i,j,k,m,ii,inum,itype,jtype;
tagint itag,jtag;
@ -497,7 +497,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,vflag_atom);
bij = bondorder(i,j,del,rij,VA,f);
dVAdi = bij*dVA;
fpair = -(dVRdi+dVAdi) / rij;
@ -520,7 +520,7 @@ void PairAIREBO::FREBO(int eflag, int /*vflag*/)
find 3- and 4-step paths between atoms I,J via REBO neighbor lists
------------------------------------------------------------------------- */
void PairAIREBO::FLJ(int eflag, int /*vflag*/)
void PairAIREBO::FLJ(int eflag)
{
int i,j,k,m,ii,jj,kk,mm,inum,jnum,itype,jtype,ktype,mtype;
int atomi,atomj,atomk,atomm;
@ -786,8 +786,7 @@ 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,vflag_atom);
Stb = bondorderLJ(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,f);
} else Stb = 0.0;
fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) +
@ -815,7 +814,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
f[atomj][1] -= delij[1]*fpair;
f[atomj][2] -= delij[2]*fpair;
if (vflag_atom) v_tally2(atomi,atomj,fpair,delij);
if (vflag_either) v_tally2(atomi,atomj,fpair,delij);
} else if (npath == 3) {
fpair1 = dC*dwikS*wkjS / rikS;
@ -837,7 +836,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
f[atomk][1] -= fi[1] + fj[1];
f[atomk][2] -= fi[2] + fj[2];
if (vflag_atom)
if (vflag_either)
v_tally3(atomi,atomj,atomk,fi,fj,delikS,deljkS);
} else if (npath == 4) {
@ -873,7 +872,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
f[atomm][1] += fm[1];
f[atomm][2] += fm[2];
if (vflag_atom) {
if (vflag_either) {
delimS[0] = delikS[0] + delkmS[0];
delimS[1] = delikS[1] + delkmS[1];
delimS[2] = delikS[2] + delkmS[2];
@ -889,7 +888,7 @@ void PairAIREBO::FLJ(int eflag, int /*vflag*/)
torsional forces and energy
------------------------------------------------------------------------- */
void PairAIREBO::TORSION(int eflag, int /*vflag*/)
void PairAIREBO::TORSION(int eflag)
{
int i,j,k,l,ii,inum;
tagint itag,jtag;
@ -1249,9 +1248,7 @@ void PairAIREBO::TORSION(int eflag, int /*vflag*/)
Bij function
------------------------------------------------------------------------- */
double PairAIREBO::bondorder(int i, int j, double rij[3],
double rijmag, double VA,
double **f, int vflag_atom)
double PairAIREBO::bondorder(int i, int j, double rij[3], double rijmag, double VA, double **f)
{
int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
int itype,jtype,ktype,ltype,ntype;
@ -1435,7 +1432,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
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) {
if (vflag_either) {
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);
@ -1577,7 +1574,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
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) {
if (vflag_either) {
rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj);
}
@ -1613,7 +1610,7 @@ 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 (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -1627,7 +1624,7 @@ 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 (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -1649,7 +1646,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
}
}
}
@ -1680,7 +1677,7 @@ 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 (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -1694,7 +1691,7 @@ 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 (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -1716,7 +1713,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
}
}
}
@ -1928,7 +1925,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
f[atom4][2] += f4[2];
if (vflag_atom) {
if (vflag_either) {
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);
@ -1964,7 +1961,7 @@ 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 (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -1978,7 +1975,7 @@ 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 (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -2000,7 +1997,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
}
}
}
@ -2031,7 +2028,7 @@ 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 (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -2045,7 +2042,7 @@ 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 (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -2067,7 +2064,7 @@ double PairAIREBO::bondorder(int i, int j, double rij[3],
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
}
}
}
@ -2113,8 +2110,7 @@ but of the vector r_ij.
*/
double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rijmag_mod,
double VA, double rij[3], double rijmag,
double **f, int vflag_atom)
double VA, double rij[3], double rijmag, double **f)
{
int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
int itype,jtype,ktype,ltype,ntype;
@ -2438,7 +2434,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
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) {
if (vflag_either) {
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);
@ -2542,7 +2538,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
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) {
if (vflag_either) {
rlj[0] = -rjl[0]; rlj[1] = -rjl[1]; rlj[2] = -rjl[2];
v_tally3(atomi,atomj,atoml,fi,fl,rij,rlj);
}
@ -2577,7 +2573,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -2591,7 +2587,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -2613,7 +2609,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
}
}
}
@ -2644,7 +2640,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -2658,7 +2654,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -2680,7 +2676,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
}
}
}
@ -2885,7 +2881,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
f[atom4][2] += f4[2];
if (vflag_atom) {
if (vflag_either) {
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);
@ -2919,7 +2915,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -2933,7 +2929,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2(atomi,atomk,-tmp2,rik);
if (vflag_either) v_tally2(atomi,atomk,-tmp2,rik);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -2955,7 +2951,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2(atomk,atomn,-tmp2,rkn);
if (vflag_either) v_tally2(atomk,atomn,-tmp2,rkn);
}
}
}
@ -2986,7 +2982,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -3000,7 +2996,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2(atomj,atoml,-tmp2,rjl);
if (vflag_either) v_tally2(atomj,atoml,-tmp2,rjl);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -3022,7 +3018,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double /* rij_mod */[3], double rij
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2(atoml,atomn,-tmp2,rln);
if (vflag_either) v_tally2(atoml,atomn,-tmp2,rln);
}
}
}

10
src/MANYBODY/pair_airebo.h
View File

@ -96,12 +96,12 @@ class PairAIREBO : public Pair {
double Tf[5][5][10], Tdfdx[5][5][10], Tdfdy[5][5][10], Tdfdz[5][5][10];
void REBO_neigh();
void FREBO(int, int);
void FLJ(int, int);
void TORSION(int, int);
void FREBO(int);
void FLJ(int);
void TORSION(int);
double bondorder(int, int, double *, double, double, double **, int);
double bondorderLJ(int, int, double *, double, double, double *, double, double **, int);
double bondorder(int, int, double *, double, double, double **);
double bondorderLJ(int, int, double *, double, double, double *, double, double **);
double gSpline(double, double, int, double *, double *);
double PijSpline(double, double, int, int, double *);

87
src/OPENMP/pair_airebo_omp.cpp
View File

@ -70,9 +70,9 @@ void PairAIREBOOMP::compute(int eflag, int vflag)
thr->timer(Timer::START);
ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
FREBO_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv0,thr);
if (ljflag) FLJ_thr(ifrom,ito,evflag,eflag,vflag_atom,&pv1,thr);
if (torflag) TORSION_thr(ifrom,ito,evflag,eflag,&pv2,thr);
FREBO_thr(ifrom,ito,eflag,&pv0,thr);
if (ljflag) FLJ_thr(ifrom,ito,eflag,&pv1,thr);
if (torflag) TORSION_thr(ifrom,ito,eflag,&pv2,thr);
thr->timer(Timer::PAIR);
reduce_thr(this, eflag, vflag, thr);
@ -184,8 +184,7 @@ void PairAIREBOOMP::REBO_neigh_thr()
REBO forces and energy
------------------------------------------------------------------------- */
void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
int vflag_atom, double *pv0, ThrData * const thr)
void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData * const thr)
{
int i,j,k,m,ii,itype,jtype;
tagint itag,jtag;
@ -259,7 +258,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
del[0] = delx;
del[1] = dely;
del[2] = delz;
bij = bondorder_thr(i,j,del,rij,VA,vflag_atom,thr);
bij = bondorder_thr(i,j,del,rij,VA,thr);
dVAdi = bij*dVA;
fpair = -(dVRdi+dVAdi) / rij;
@ -282,8 +281,7 @@ void PairAIREBOOMP::FREBO_thr(int ifrom, int ito, int evflag, int eflag,
find 3- and 4-step paths between atoms I,J via REBO neighbor lists
------------------------------------------------------------------------- */
void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
int vflag_atom, double *pv1, ThrData * const thr)
void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData * const thr)
{
int i,j,k,m,ii,jj,kk,mm,jnum,itype,jtype,ktype,mtype;
tagint itag,jtag;
@ -547,8 +545,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
delscale[0] = scale * delij[0];
delscale[1] = scale * delij[1];
delscale[2] = scale * delij[2];
Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,
delij,rij,vflag_atom,thr);
Stb = bondorderLJ_thr(i,j,delscale,rcmin[itype][jtype],VA,delij,rij,thr);
} else Stb = 0.0;
fpair = -(dStr * (Stb*cij*VLJ - cij*VLJ) +
@ -576,7 +573,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
f[atomj][1] -= delij[1]*fpair;
f[atomj][2] -= delij[2]*fpair;
if (vflag_atom) v_tally2_thr(atomi,atomj,fpair,delij,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomj,fpair,delij,thr);
} else if (npath == 3) {
fpair1 = dC*dwikS*wkjS / rikS;
@ -634,11 +631,11 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
f[atomm][1] += fm[1];
f[atomm][2] += fm[2];
if (vflag_atom) {
if (vflag_either) {
delimS[0] = delikS[0] + delkmS[0];
delimS[1] = delikS[1] + delkmS[1];
delimS[2] = delikS[2] + delkmS[2];
v_tally4_thr(atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr);
v_tally4_thr(this,atomi,atomj,atomk,atomm,fi,fj,fk,delimS,deljmS,delkmS,thr);
}
}
}
@ -650,8 +647,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
torsional forces and energy
------------------------------------------------------------------------- */
void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag,
double *pv2, ThrData * const thr)
void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData * const thr)
{
int i,j,k,l,ii;
tagint itag,jtag;
@ -1011,7 +1007,7 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, int evflag, int eflag,
------------------------------------------------------------------------- */
double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
double VA, int vflag_atom, ThrData * const thr)
double VA, ThrData * const thr)
{
int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
int itype,jtype,ktype,ltype,ntype;
@ -1378,7 +1374,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -1392,7 +1388,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -1414,7 +1410,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
}
}
}
@ -1445,7 +1441,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -1459,7 +1455,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -1481,7 +1477,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
}
}
}
@ -1693,10 +1689,10 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
f[atom4][2] += f4[2];
if (vflag_atom) {
if (vflag_either) {
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_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
}
}
}
@ -1729,7 +1725,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -1743,7 +1739,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -1765,7 +1761,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
}
}
}
@ -1796,7 +1792,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -1810,7 +1806,7 @@ double PairAIREBOOMP::bondorder_thr(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_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -1832,7 +1828,7 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag,
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
}
}
}
@ -1870,8 +1866,7 @@ there probably also need to be performed here.
*/
double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], double rijmag_mod,
double VA, double rij[3], double rijmag,
int vflag_atom, ThrData * const thr)
double VA, double rij[3], double rijmag, ThrData * const thr)
{
int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4;
int itype,jtype,ktype,ltype,ntype;
@ -2335,7 +2330,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -2349,7 +2344,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -2371,7 +2366,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
}
}
}
@ -2402,7 +2397,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -2416,7 +2411,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -2438,7 +2433,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
}
}
}
@ -2643,10 +2638,10 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atom4][0] += f4[0]; f[atom4][1] += f4[1];
f[atom4][2] += f4[2];
if (vflag_atom) {
if (vflag_either) {
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_thr(atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
v_tally4_thr(this,atom1,atom2,atom3,atom4,f1,f2,f4,r13,r23,r43,thr);
}
}
}
@ -2677,7 +2672,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
// due to kronecker(ktype, 0) term in contribution
// to NconjtmpI and later Nijconj
@ -2691,7 +2686,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomk][1] += tmp2*rik[1];
f[atomk][2] += tmp2*rik[2];
if (vflag_atom) v_tally2_thr(atomi,atomk,-tmp2,rik,thr);
if (vflag_either) v_tally2_thr(this,atomi,atomk,-tmp2,rik,thr);
if (fabs(dNki) > TOL) {
REBO_neighs_k = REBO_firstneigh[atomk];
@ -2713,7 +2708,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomn][1] += tmp2*rkn[1];
f[atomn][2] += tmp2*rkn[2];
if (vflag_atom) v_tally2_thr(atomk,atomn,-tmp2,rkn,thr);
if (vflag_either) v_tally2_thr(this,atomk,atomn,-tmp2,rkn,thr);
}
}
}
@ -2744,7 +2739,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
// due to kronecker(ltype, 0) term in contribution
// to NconjtmpJ and later Nijconj
@ -2758,7 +2753,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atoml][1] += tmp2*rjl[1];
f[atoml][2] += tmp2*rjl[2];
if (vflag_atom) v_tally2_thr(atomj,atoml,-tmp2,rjl,thr);
if (vflag_either) v_tally2_thr(this,atomj,atoml,-tmp2,rjl,thr);
if (fabs(dNlj) > TOL) {
REBO_neighs_l = REBO_firstneigh[atoml];
@ -2780,7 +2775,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double /* rij_mod */[3], dou
f[atomn][1] += tmp2*rln[1];
f[atomn][2] += tmp2*rln[2];
if (vflag_atom) v_tally2_thr(atoml,atomn,-tmp2,rln,thr);
if (vflag_either) v_tally2_thr(this,atoml,atomn,-tmp2,rln,thr);
}
}
}

13
src/OPENMP/pair_airebo_omp.h
View File

@ -33,16 +33,13 @@ class PairAIREBOOMP : public PairAIREBO, public ThrOMP {
virtual double memory_usage();
protected:
double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, int vflag_atom,
ThrData *const thr);
double bondorder_thr(int i, int j, double rij[3], double rijmag, double VA, ThrData *const thr);
double bondorderLJ_thr(int i, int j, double rij[3], double rijmag, double VA, double rij0[3],
double rijmag0, int vflag_atom, ThrData *const thr);
double rijmag0, ThrData *const thr);
void FREBO_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv0,
ThrData *const thr);
void FLJ_thr(int ifrom, int ito, int evflag, int eflag, int vflag_atom, double *pv1,
ThrData *const thr);
void TORSION_thr(int ifrom, int ito, int evflag, int eflag, double *pv2, ThrData *const thr);
void FREBO_thr(int ifrom, int ito, int eflag, double *pv0, ThrData *const thr);
void FLJ_thr(int ifrom, int ito, int eflag, double *pv1, ThrData *const thr);
void TORSION_thr(int ifrom, int ito, int eflag, double *pv2, ThrData *const thr);
void REBO_neigh_thr();
};

50
src/OPENMP/thr_omp.cpp
View File

@ -1465,20 +1465,31 @@ void ThrOMP::ev_tally_thr(Improper * const imprp, const int i1, const int i2,
fpair is magnitude of force on atom I
------------------------------------------------------------------------- */
void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
void ThrOMP::v_tally2_thr(Pair *const pair, const int i, const int j, const double fpair,
const double * const drij, ThrData * const thr)
{
double v[6];
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;
v[0] = drij[0]*drij[0]*fpair;
v[1] = drij[1]*drij[1]*fpair;
v[2] = drij[2]*drij[2]*fpair;
v[3] = drij[0]*drij[1]*fpair;
v[4] = drij[0]*drij[2]*fpair;
v[5] = drij[1]*drij[2]*fpair;
if (pair->vflag_global) v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
v[0] *= 0.5;
v[1] *= 0.5;
v[2] *= 0.5;
v[3] *= 0.5;
v[4] *= 0.5;
v[5] *= 0.5;
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
}
}
/* ----------------------------------------------------------------------
@ -1486,7 +1497,7 @@ void ThrOMP::v_tally2_thr(const int i, const int j, const double fpair,
called by AIREBO and Tersoff potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k,
void ThrOMP::v_tally3_thr(Pair *const pair, const int i, const int j, const int k,
const double * const fi, const double * const fj,
const double * const drik, const double * const drjk,
ThrData * const thr)
@ -1521,7 +1532,7 @@ void ThrOMP::v_tally3_thr(Pair *pair, const int i, const int j, const int k,
called by AIREBO potential, newton_pair is always on
------------------------------------------------------------------------- */
void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m,
void ThrOMP::v_tally4_thr(Pair *const pair, const int i, const int j, const int k, const int m,
const double * const fi, const double * const fj,
const double * const fk, const double * const drim,
const double * const drjm, const double * const drkm,
@ -1529,17 +1540,26 @@ void ThrOMP::v_tally4_thr(const int i, const int j, const int k, const int m,
{
double v[6];
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]);
v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
if (pair->vflag_global) v_tally(thr->virial_pair,v);
if (pair->vflag_atom) {
v[0] *= 0.25;
v[1] *= 0.25;
v[2] *= 0.25;
v[3] *= 0.25;
v[4] *= 0.25;
v[5] *= 0.25;
v_tally(thr->vatom_pair[i],v);
v_tally(thr->vatom_pair[j],v);
v_tally(thr->vatom_pair[k],v);
v_tally(thr->vatom_pair[m],v);
}
}
/* ---------------------------------------------------------------------- */

9
src/OPENMP/thr_omp.h
View File

@ -135,6 +135,8 @@ class ThrOMP {
void ev_tally_xyz_full_thr(Pair *const, const int, const double, const double, const double,
const double, const double, const double, const double, const double,
ThrData *const);
void v_tally2_thr(Pair *const, const int, const int, const double, const double *const,
ThrData *const);
void ev_tally3_thr(Pair *const, const int, const int, const int, const double, const double,
const double *const, const double *const, const double *const,
const double *const, ThrData *const);
@ -143,6 +145,9 @@ class ThrOMP {
void ev_tally4_thr(Pair *const, const int, const int, const int, const int, const double,
const double *const, const double *const, const double *const,
const double *const, const double *const, const double *const, ThrData *const);
void v_tally4_thr(Pair *const, const int, const int, const int, const int, const double *const,
const double *const, const double *const, const double *const,
const double *const, const double *const, ThrData *const);
// Bond
void ev_tally_thr(Bond *const, const int, const int, const int, const int, const double,
@ -171,10 +176,6 @@ class ThrOMP {
ThrData *const);
// style independent versions
void v_tally2_thr(const int, const int, const double, const double *const, ThrData *const);
void v_tally4_thr(const int, const int, const int, const int, const double *const,
const double *const, const double *const, const double *const,
const double *const, const double *const, ThrData *const);
void ev_tally_list_thr(Pair *const, const int, const int *const, const double *const,
const double, const double, ThrData *const);
};

90
src/pair.cpp
View File

@ -1499,17 +1499,35 @@ void Pair::v_tally2(int i, int j, double fpair, double *drij)
{
double v[6];
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;
v[0] = drij[0]*drij[0]*fpair;
v[1] = drij[1]*drij[1]*fpair;
v[2] = drij[2]*drij[2]*fpair;
v[3] = drij[0]*drij[1]*fpair;
v[4] = drij[0]*drij[2]*fpair;
v[5] = drij[1]*drij[2]*fpair;
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] *= 0.5;
v[1] *= 0.5;
v[2] *= 0.5;
v[3] *= 0.5;
v[4] *= 0.5;
v[5] *= 0.5;
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];
}
}
/* ----------------------------------------------------------------------
@ -1522,23 +1540,29 @@ void Pair::v_tally3(int i, int j, int k,
{
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]);
v[0] = (drik[0]*fi[0] + drjk[0]*fj[0]);
v[1] = (drik[1]*fi[1] + drjk[1]*fj[1]);
v[2] = (drik[2]*fi[2] + drjk[2]*fj[2]);
v[3] = (drik[0]*fi[1] + drjk[0]*fj[1]);
v[4] = (drik[0]*fi[2] + drjk[0]*fj[2]);
v[5] = (drik[1]*fi[2] + drjk[1]*fj[2]);
if (vflag_global) {
virial[0] += 3.0*v[0];
virial[1] += 3.0*v[1];
virial[2] += 3.0*v[2];
virial[3] += 3.0*v[3];
virial[4] += 3.0*v[4];
virial[5] += 3.0*v[5];
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] *= THIRD;
v[1] *= THIRD;
v[2] *= THIRD;
v[3] *= THIRD;
v[4] *= THIRD;
v[5] *= THIRD;
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];
@ -1559,12 +1583,29 @@ void Pair::v_tally4(int i, int j, int k, int m,
{
double v[6];
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]);
v[0] = (drim[0]*fi[0] + drjm[0]*fj[0] + drkm[0]*fk[0]);
v[1] = (drim[1]*fi[1] + drjm[1]*fj[1] + drkm[1]*fk[1]);
v[2] = (drim[2]*fi[2] + drjm[2]*fj[2] + drkm[2]*fk[2]);
v[3] = (drim[0]*fi[1] + drjm[0]*fj[1] + drkm[0]*fk[1]);
v[4] = (drim[0]*fi[2] + drjm[0]*fj[2] + drkm[0]*fk[2]);
v[5] = (drim[1]*fi[2] + drjm[1]*fj[2] + drkm[1]*fk[2]);
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] *= 0.25;
v[1] *= 0.25;
v[2] *= 0.25;
v[3] *= 0.25;
v[4] *= 0.25;
v[5] *= 0.25;
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];
@ -1574,6 +1615,7 @@ void Pair::v_tally4(int i, int j, int k, int m,
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];
}
}
/* ----------------------------------------------------------------------

10
unittest/force-styles/test_pair_style.cpp
View File

@ -1089,12 +1089,14 @@ TEST(PairStyle, intel)
EXPECT_FP_LE_WITH_EPS((pair->eng_vdwl + pair->eng_coul), energy, epsilon);
if (print_stats) std::cerr << "run_energy stats:" << stats << std::endl;
// pair styles sw and tersoff require newton on, but that also requires fdotr for /intel
std::cerr << "pair style : " << test_config.pair_style << "\n";
if ((test_config.pair_style != "sw") && (test_config.pair_style != "tersoff")) {
// pair styles sw and tersoff and airebo INTEL package variants require newton on,
// but that also requires fdotr for /intel
if ((test_config.pair_style != "sw") && (test_config.pair_style != "tersoff")
&& (test_config.pair_style != "rebo") && (test_config.pair_style != "airebo")
&& (test_config.pair_style != "airebo/morse")) {
if (!verbose) ::testing::internal::CaptureStdout();
restart_lammps(lmp, test_config, true, false);
restart_lammps(lmp, test_config, true, true);
if (!verbose) ::testing::internal::GetCapturedStdout();
f = lmp->atom->f;
tag = lmp->atom->tag;