diff --git a/src/pair.cpp b/src/pair.cpp index 18d561bdb5..f76bb54c36 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -301,7 +301,7 @@ void Pair::init_style() specific pair style can override this function ------------------------------------------------------------------------- */ -void Pair::init_list(int /*which*/, NeighList *ptr) +void Pair::init_list(int which, NeighList *ptr) { list = ptr; } @@ -1740,6 +1740,18 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits, /* ---------------------------------------------------------------------- */ +void Pair::pairTensor(double fforce, double dfac, double delr[3], double phiTensor[6]) { + int m = 0; + for (int k=0; k<3; k++) { + phiTensor[m] = fforce; + for (int l=k; l<3; l++) { + if (l>k) phiTensor[m] = 0; + phiTensor[m++] += delr[k]*delr[l] * dfac; + } + } +} +/* ---------------------------------------------------------------------- */ + double Pair::memory_usage() { double bytes = comm->nthreads*maxeatom * sizeof(double); diff --git a/src/pair.h b/src/pair.h index c1a9e6eef8..a9506480be 100644 --- a/src/pair.h +++ b/src/pair.h @@ -145,6 +145,16 @@ class Pair : protected Pointers { return 0.0; } + void pairTensor(double fforce, double dfac, double delr[3], double phiTensor[6]); + + virtual double single2(int, int, int, int, + double, double[3], double, double, + double& fforce, double d2u[6]) { + fforce = 0.0; + for (int i=0; i<6; i++) d2u[i] = 0; + return 0.0; + } + virtual void settings(int, char **) = 0; virtual void coeff(int, char **) = 0; @@ -207,10 +217,6 @@ class Pair : protected Pointers { typedef union {int i; float f;} union_int_float_t; - // Accessor for the user-intel package to determine virial calc for hybrid - - inline int fdotr_is_set() const { return vflag_fdotr; } - protected: int vflag_fdotr; int maxeatom,maxvatom; diff --git a/src/pair_lj_smooth_linear.cpp b/src/pair_lj_smooth_linear.cpp index 17c789bcee..4cf8d6d3db 100644 --- a/src/pair_lj_smooth_linear.cpp +++ b/src/pair_lj_smooth_linear.cpp @@ -326,9 +326,9 @@ void PairLJSmoothLinear::read_restart_settings(FILE *fp) /* ---------------------------------------------------------------------- */ -double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype, +double PairLJSmoothLinear::single(int i, int j, int itype, int jtype, double rsq, - double /*factor_coul*/, double factor_lj, + double factor_coul, double factor_lj, double &fforce) { double r2inv,r6inv,forcelj,philj,r,rinv; @@ -347,3 +347,26 @@ double PairLJSmoothLinear::single(int /*i*/, int /*j*/, int itype, int jtype, return factor_lj*philj; } + +double PairLJSmoothLinear::single2(int i, int j, int itype, int jtype, double rsq, + double delr[3], double factor_coul, double factor_lj, + double &fforce, double d2u[6]) +{ + double r2inv,r6inv,forcelj,philj,r,rinv; + + r2inv = 1.0/rsq; + r6inv = r2inv*r2inv*r2inv; + rinv = sqrt(r2inv); + r = sqrt(rsq); + forcelj = r6inv*(lj1[itype][jtype]*r6inv-lj2[itype][jtype]); + forcelj = rinv*forcelj - dljcut[itype][jtype]; + fforce = factor_lj*forcelj*rinv; + + philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]); + philj = philj - ljcut[itype][jtype] + + (r-cut[itype][jtype])*dljcut[itype][jtype]; + + double d2r = factor_lj * r6inv * (13.0*lj1[itype][jtype]*r6inv - 7.0*lj2[itype][jtype])/rsq; + pairTensor(fforce, -(fforce + d2r) / rsq, delr, d2u); + return factor_lj*philj; +} diff --git a/src/pair_lj_smooth_linear.h b/src/pair_lj_smooth_linear.h index c18c442a18..9cee9f2951 100644 --- a/src/pair_lj_smooth_linear.h +++ b/src/pair_lj_smooth_linear.h @@ -38,6 +38,7 @@ class PairLJSmoothLinear : public Pair { void write_restart_settings(FILE *); void read_restart_settings(FILE *); double single(int, int, int, int, double, double, double, double &); + double single2(int, int, int, int, double, double[3], double, double, double&, double[6]); protected: double cut_global;