update several MANYBODY potentials for pair_modify nofdotr

This commit is contained in:
Axel Kohlmeyer
2021-07-11 08:20:20 -04:00
parent 9898941169
commit a4748b4c28
4 changed files with 31 additions and 38 deletions

View File

@ -227,7 +227,7 @@ void PairGW::compute(int eflag, int vflag)
f[k][1] += fk[1];
f[k][2] += fk[2];
if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2);
} // kk
} // jj
} // ii

View File

@ -67,7 +67,7 @@ PairLCBOP::~PairLCBOP()
{
memory->destroy(SR_numneigh);
memory->sfree(SR_firstneigh);
delete [] ipage;
delete[] ipage;
memory->destroy(N);
memory->destroy(M);
@ -167,7 +167,7 @@ void PairLCBOP::init_style()
if (oneatom != neighbor->oneatom) create = 1;
if (create) {
delete [] ipage;
delete[] ipage;
pgsize = neighbor->pgsize;
oneatom = neighbor->oneatom;
@ -386,7 +386,7 @@ void PairLCBOP::FSR(int eflag, int /*vflag*/)
del[0] = delx;
del[1] = dely;
del[2] = delz;
Bij = bondorder(i,j,del,rijmag,VA,f,vflag_atom);
Bij = bondorder(i,j,del,rijmag,VA,f);
dVAdi = Bij*dVA;
// F = (dVRdi+dVAdi)*(-grad rijmag)
@ -402,8 +402,7 @@ void PairLCBOP::FSR(int eflag, int /*vflag*/)
double evdwl=0.0;
if (eflag) evdwl = VR - Bij*VA;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}
@ -493,8 +492,7 @@ void PairLCBOP::FLR(int eflag, int /*vflag*/)
double evdwl=0.0;
if (eflag) evdwl = V;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);
}
}
}
@ -503,7 +501,7 @@ void PairLCBOP::FLR(int eflag, int /*vflag*/)
forces for Nij and Mij
------------------------------------------------------------------------- */
void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) {
void PairLCBOP::FNij( int i, int j, double factor, double **f) {
int atomi = i;
int atomj = j;
int *SR_neighs = SR_firstneigh[i];
@ -532,7 +530,7 @@ void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) {
f[atomk][1] -= rik[1]*fpair;
f[atomk][2] -= rik[2]*fpair;
if (vflag_atom) v_tally2(atomi,atomk,fpair,rik);
if (vflag_either) v_tally2(atomi,atomk,fpair,rik);
}
}
}
@ -540,7 +538,7 @@ void PairLCBOP::FNij( int i, int j, double factor, double **f, int vflag_atom) {
/* ---------------------------------------------------------------------- */
void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) {
void PairLCBOP::FMij( int i, int j, double factor, double **f) {
int atomi = i;
int atomj = j;
int *SR_neighs = SR_firstneigh[i];
@ -573,12 +571,12 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) {
f[atomk][0] -= rik[0]*fpair;
f[atomk][1] -= rik[1]*fpair;
f[atomk][2] -= rik[2]*fpair;
if (vflag_atom) v_tally2(atomi,atomk,fpair,rik);
if (vflag_either) v_tally2(atomi,atomk,fpair,rik);
}
if (dF > TOL) {
double factor2 = factor*f_c_ik*dF;
FNij( atomk, atomi, factor2, f, vflag_atom );
FNij(atomk, atomi, factor2, f);
}
}
}
@ -588,17 +586,15 @@ void PairLCBOP::FMij( int i, int j, double factor, double **f, int vflag_atom) {
Bij function
------------------------------------------------------------------------- */
double PairLCBOP::bondorder(int i, int j, double rij[3],
double rijmag, double VA,
double **f, int vflag_atom)
double PairLCBOP::bondorder(int i, int j, double rij[3],double rijmag, double VA,double **f)
{
double bij, bji;
/* bij & bji */{
double rji[3];
rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
bij = b(i,j,rij,rijmag,VA,f,vflag_atom);
bji = b(j,i,rji,rijmag,VA,f,vflag_atom);
bij = b(i,j,rij,rijmag,VA,f);
bji = b(j,i,rji,rijmag,VA,f);
}
double Fij_conj;
@ -663,31 +659,30 @@ double PairLCBOP::bondorder(int i, int j, double rij[3],
}
double dF_dNij, dF_dNji, dF_dNconj;
Fij_conj = F_conj( Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj );
Fij_conj = F_conj(Nij, Nji, Nconj, &dF_dNij, &dF_dNji, &dF_dNconj);
/*forces for Nij*/
if (3-Nij > TOL) {
double factor = -VA*0.5*( dF_dNij + dF_dNconj*( dNconj_dNij + dNconj_dNel*dNij_el_dNij ) );
FNij( i, j, factor, f, vflag_atom );
double factor = -VA*0.5*(dF_dNij + dF_dNconj*(dNconj_dNij + dNconj_dNel*dNij_el_dNij));
FNij(i, j, factor, f);
}
/*forces for Nji*/
if (3-Nji > TOL) {
double factor = -VA*0.5*( dF_dNji + dF_dNconj*( dNconj_dNji + dNconj_dNel*dNji_el_dNji ) );
FNij( j, i, factor, f, vflag_atom );
double factor = -VA*0.5*(dF_dNji + dF_dNconj*(dNconj_dNji + dNconj_dNel*dNji_el_dNji));
FNij(j, i, factor, f);
}
/*forces for Mij*/
if (3-Mij > TOL) {
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNij_el_dMij );
FMij( i, j, factor, f, vflag_atom );
double factor = -VA*0.5*(dF_dNconj*dNconj_dNel*dNij_el_dMij);
FMij(i, j, factor, f);
}
if (3-Mji > TOL) {
double factor = -VA*0.5*( dF_dNconj*dNconj_dNel*dNji_el_dMji );
FMij( j, i, factor, f, vflag_atom );
double factor = -VA*0.5*(dF_dNconj*dNconj_dNel*dNji_el_dMji);
FMij(j,i,factor,f);
}
}
double Bij = 0.5*( bij + bji + Fij_conj );
double Bij = 0.5*(bij + bji + Fij_conj);
return Bij;
}
@ -695,9 +690,7 @@ double PairLCBOP::bondorder(int i, int j, double rij[3],
bij function
------------------------------------------------------------------------- */
double PairLCBOP::b(int i, int j, double rij[3],
double rijmag, double VA,
double **f, int vflag_atom) {
double PairLCBOP::b(int i, int j, double rij[3], double rijmag, double VA, double **f) {
int *SR_neighs = SR_firstneigh[i];
double **x = atom->x;
int atomi = i;
@ -818,7 +811,7 @@ double PairLCBOP::b(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) {
double rji[3], rki[3];
rji[0] = -rij[0]; rji[1] = -rij[1]; rji[2] = -rij[2];
rki[0] = -rik[0]; rki[1] = -rik[1]; rki[2] = -rik[2];

View File

@ -76,10 +76,10 @@ class PairLCBOP : public Pair {
void FSR(int, int);
void FLR(int, int);
void FNij(int, int, double, double **, int);
void FMij(int, int, double, double **, int);
double bondorder(int, int, double *, double, double, double **, int);
double b(int, int, double *, double, double, double **, int);
void FNij(int, int, double, double **);
void FMij(int, int, double, double **);
double bondorder(int, int, double *, double, double, double **);
double b(int, int, double *, double, double, double **);
double gSpline(double, double *);
double hSpline(double, double *);

View File

@ -410,7 +410,7 @@ void PairExTeP::compute(int eflag, int vflag)
f[k][1] += fk[1];
f[k][2] += fk[2];
if (vflag_atom) v_tally3(i,j,k,fj,fk,delr1,delr2);
if (vflag_either) v_tally3(i,j,k,fj,fk,delr1,delr2);
}
}
}