propagate DPD exclusion changes to INTEL and KOKKOS packages

This commit is contained in:
Axel Kohlmeyer
2023-01-02 11:33:08 -05:00
parent 8fb1193a0b
commit 37b3ba827f
9 changed files with 74 additions and 63 deletions

View File

@ -325,8 +325,11 @@ void PairDPDIntel::eval(const int offload, const int vflag,
const flt_t rinv = (flt_t)1.0/sqrt(rsq); const flt_t rinv = (flt_t)1.0/sqrt(rsq);
if (rinv > icut) { if (rinv > icut) {
flt_t factor_dpd; flt_t factor_dpd, factor_sqrt;
if (!ONETYPE) factor_dpd = special_lj[sbindex]; if (!ONETYPE) {
factor_dpd = special_lj[sbindex];
factor_sqrt = special_lj[sbindex];
}
flt_t delvx = vxtmp - v[j].x; flt_t delvx = vxtmp - v[j].x;
flt_t delvy = vytmp - v[j].y; flt_t delvy = vytmp - v[j].y;
@ -342,8 +345,11 @@ void PairDPDIntel::eval(const int offload, const int vflag,
gamma = parami[jtype].gamma; gamma = parami[jtype].gamma;
sigma = parami[jtype].sigma; sigma = parami[jtype].sigma;
} }
flt_t fpair = a0 - iwd * gamma * dot + sigma * randnum * dtinvsqrt; flt_t fpair = a0 - iwd * gamma * dot;
if (!ONETYPE) fpair *= factor_dpd; if (!ONETYPE) {
fpair *= factor_dpd;
fpair += factor_sqrt * sigma * randnum * dtinvsqrt;
} else fpair += sigma * randnum * dtinvsqrt;
fpair *= iwd; fpair *= iwd;
const flt_t fpx = fpair * delx; const flt_t fpx = fpair * delx;
@ -493,8 +499,7 @@ void PairDPDIntel::init_style()
fix->pair_init_check(); fix->pair_init_check();
#ifdef _LMP_INTEL_OFFLOAD #ifdef _LMP_INTEL_OFFLOAD
if (fix->offload_balance() != 0.0) if (fix->offload_balance() != 0.0)
error->all(FLERR, error->all(FLERR, "Offload for dpd/intel is not yet available. Set balance to 0.");
"Offload for dpd/intel is not yet available. Set balance to 0.");
#endif #endif
if (fix->precision() == FixIntel::PREC_MODE_MIXED) if (fix->precision() == FixIntel::PREC_MODE_MIXED)

View File

@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
template<class DeviceType> template<class DeviceType>
PairDPDExtKokkos<DeviceType>::PairDPDExtKokkos(class LAMMPS *lmp) : PairDPDExtKokkos<DeviceType>::PairDPDExtKokkos(class LAMMPS *_lmp) :
PairDPDExt(lmp) , PairDPDExt(_lmp) ,
#ifdef DPD_USE_RAN_MARS #ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp) rand_pool(0 /* unused */, lmp)
#else #else
@ -134,6 +134,10 @@ void PairDPDExtKokkos<DeviceType>::compute(int eflagin, int vflagin)
special_lj[1] = force->special_lj[1]; special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2]; special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3]; special_lj[3] = force->special_lj[3];
special_rf[0] = sqrt(force->special_lj[0]);
special_rf[1] = sqrt(force->special_lj[1]);
special_rf[2] = sqrt(force->special_lj[2]);
special_rf[3] = sqrt(force->special_lj[3]);
nlocal = atom->nlocal; nlocal = atom->nlocal;
dtinvsqrt = 1.0/sqrt(update->dt); dtinvsqrt = 1.0/sqrt(update->dt);
@ -232,7 +236,8 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
int i,j,jj,jnum,itype,jtype; int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair; double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz; double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd; double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz;
double prefactor_g,prefactor_s,factor_dpd,factor_sqrt;
double fx = 0,fy = 0,fz = 0; double fx = 0,fy = 0,fz = 0;
double evdwl = 0; double evdwl = 0;
i = d_ilist[ii]; i = d_ilist[ii];
@ -249,6 +254,7 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
double P[3][3]; double P[3][3];
j = d_neighbors(i,jj); j = d_neighbors(i,jj);
factor_dpd = special_lj[sbmask(j)]; factor_dpd = special_lj[sbmask(j)];
factor_sqrt = special_rf[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
delx = xtmp - x(j,0); delx = xtmp - x(j,0);
@ -291,33 +297,26 @@ void PairDPDExtKokkos<DeviceType>::operator() (TagDPDExtKokkos<NEIGHFLAG,EVFLAG>
// drag force - parallel // drag force - parallel
fpair -= params(itype,jtype).gamma*wdPar*wdPar*dot*rinv; fpair -= params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
fpair *= factor_dpd;
// random force - parallel // random force - parallel
fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt; fpair += factor_sqrt*params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
fpairx = fpair*rinv*delx; fpairx = fpair*rinv*delx;
fpairy = fpair*rinv*dely; fpairy = fpair*rinv*dely;
fpairz = fpair*rinv*delz; fpairz = fpair*rinv*delz;
// drag force - perpendicular // drag force - perpendicular
fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp* prefactor_g = factor_dpd*params(itype,jtype).gammaT*wdPerp*wdPerp;
(P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz); fpairx -= prefactor_g * (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp* fpairy -= prefactor_g * (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
(P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz); fpairz -= prefactor_g * (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
// random force - perpendicular // random force - perpendicular
fpairx += params(itype,jtype).sigmaT*wdPerp* prefactor_s = factor_sqrt*params(itype,jtype).sigmaT*wdPerp;
(P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt; fpairx += prefactor_s * (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
fpairy += params(itype,jtype).sigmaT*wdPerp* fpairy += prefactor_s * (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
(P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt; fpairz += prefactor_s * (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairz += params(itype,jtype).sigmaT*wdPerp*
(P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairx *= factor_dpd;
fpairy *= factor_dpd;
fpairz *= factor_dpd;
fx += fpairx; fx += fpairx;
fy += fpairy; fy += fpairy;

View File

@ -80,7 +80,7 @@ class PairDPDExtKokkos : public PairDPDExt {
const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz, const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz,
const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const; const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
private: private:
double special_lj[4]; double special_lj[4], special_rf[4];
int eflag,vflag; int eflag,vflag;
int neighflag,nlocal; int neighflag,nlocal;
double dtinvsqrt; double dtinvsqrt;

View File

@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
template<class DeviceType> template<class DeviceType>
PairDPDExtTstatKokkos<DeviceType>::PairDPDExtTstatKokkos(class LAMMPS *lmp) : PairDPDExtTstatKokkos<DeviceType>::PairDPDExtTstatKokkos(class LAMMPS *_lmp) :
PairDPDExtTstat(lmp) , PairDPDExtTstat(_lmp) ,
#ifdef DPD_USE_RAN_MARS #ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp) rand_pool(0 /* unused */, lmp)
#else #else
@ -149,6 +149,10 @@ void PairDPDExtTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
special_lj[1] = force->special_lj[1]; special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2]; special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3]; special_lj[3] = force->special_lj[3];
special_rf[0] = sqrt(force->special_lj[0]);
special_rf[1] = sqrt(force->special_lj[1]);
special_rf[2] = sqrt(force->special_lj[2]);
special_rf[3] = sqrt(force->special_lj[3]);
nlocal = atom->nlocal; nlocal = atom->nlocal;
dtinvsqrt = 1.0/sqrt(update->dt); dtinvsqrt = 1.0/sqrt(update->dt);
@ -233,11 +237,11 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f); auto v_f = ScatterViewHelper<NeedDup_v<NEIGHFLAG,DeviceType>,decltype(dup_f),decltype(ndup_f)>::get(dup_f,ndup_f);
auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>(); auto a_f = v_f.template access<AtomicDup_v<NEIGHFLAG,DeviceType>>();
int i,j,jj,jnum,itype,jtype; int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair; double xtmp,ytmp,ztmp,delx,dely,delz,fpairx,fpairy,fpairz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz; double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz,factor_dpd; double rsq,r,rinv,dot,wd,wdPar,wdPerp,randnum,randnumx,randnumy,randnumz;
double prefactor_g,prefactor_s,factor_dpd,factor_sqrt;
double fx = 0,fy = 0,fz = 0; double fx = 0,fy = 0,fz = 0;
i = d_ilist[ii]; i = d_ilist[ii];
@ -254,6 +258,7 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
double P[3][3]; double P[3][3];
j = d_neighbors(i,jj); j = d_neighbors(i,jj);
factor_dpd = special_lj[sbmask(j)]; factor_dpd = special_lj[sbmask(j)];
factor_sqrt = special_rf[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
delx = xtmp - x(j,0); delx = xtmp - x(j,0);
@ -292,34 +297,27 @@ void PairDPDExtTstatKokkos<DeviceType>::operator() (TagDPDExtTstatKokkos<NEIGHFL
randnumz = rand_gen.normal(); randnumz = rand_gen.normal();
// drag force - parallel // drag force - parallel
fpair = -params(itype,jtype).gamma*wdPar*wdPar*dot*rinv; fpair = params(itype,jtype).gamma*wdPar*wdPar*dot*rinv;
fpair *= factor_dpd;
// random force - parallel // random force - parallel
fpair += params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt; fpair += factor_sqrt*params(itype,jtype).sigma*wdPar*randnum*dtinvsqrt;
fpairx = fpair*rinv*delx; fpairx = fpair*rinv*delx;
fpairy = fpair*rinv*dely; fpairy = fpair*rinv*dely;
fpairz = fpair*rinv*delz; fpairz = fpair*rinv*delz;
// drag force - perpendicular // drag force - perpendicular
fpairx -= params(itype,jtype).gammaT*wdPerp*wdPerp* prefactor_g = factor_dpd*params(itype,jtype).gammaT*wdPerp*wdPerp;
(P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz); fpairx -= prefactor_g * (P[0][0]*delvx + P[0][1]*delvy + P[0][2]*delvz);
fpairy -= params(itype,jtype).gammaT*wdPerp*wdPerp* fpairy -= prefactor_g * (P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz);
(P[1][0]*delvx + P[1][1]*delvy + P[1][2]*delvz); fpairz -= prefactor_g * (P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
fpairz -= params(itype,jtype).gammaT*wdPerp*wdPerp*
(P[2][0]*delvx + P[2][1]*delvy + P[2][2]*delvz);
// random force - perpendicular // random force - perpendicular
fpairx += params(itype,jtype).sigmaT*wdPerp* prefactor_s = factor_sqrt*params(itype,jtype).sigmaT*wdPerp;
(P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt; fpairx += prefactor_s * (P[0][0]*randnumx + P[0][1]*randnumy + P[0][2]*randnumz)*dtinvsqrt;
fpairy += params(itype,jtype).sigmaT*wdPerp* fpairy += prefactor_s * (P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt;
(P[1][0]*randnumx + P[1][1]*randnumy + P[1][2]*randnumz)*dtinvsqrt; fpairz += prefactor_s * (P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairz += params(itype,jtype).sigmaT*wdPerp*
(P[2][0]*randnumx + P[2][1]*randnumy + P[2][2]*randnumz)*dtinvsqrt;
fpairx *= factor_dpd;
fpairy *= factor_dpd;
fpairz *= factor_dpd;
fx += fpairx; fx += fpairx;
fy += fpairy; fy += fpairy;

View File

@ -79,7 +79,7 @@ class PairDPDExtTstatKokkos : public PairDPDExtTstat {
const F_FLOAT &fx,const F_FLOAT &fy, const F_FLOAT &fz, const F_FLOAT &fx,const F_FLOAT &fy, const F_FLOAT &fz,
const F_FLOAT &delx,const F_FLOAT &dely, const F_FLOAT &delz) const; const F_FLOAT &delx,const F_FLOAT &dely, const F_FLOAT &delz) const;
private: private:
double special_lj[4]; double special_lj[4], special_rf[4];
int eflag,vflag; int eflag,vflag;
int neighflag,nlocal; int neighflag,nlocal;
double dtinvsqrt; double dtinvsqrt;

View File

@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
template<class DeviceType> template<class DeviceType>
PairDPDKokkos<DeviceType>::PairDPDKokkos(class LAMMPS *lmp) : PairDPDKokkos<DeviceType>::PairDPDKokkos(class LAMMPS *_lmp) :
PairDPD(lmp) , PairDPD(_lmp) ,
#ifdef DPD_USE_RAN_MARS #ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp) rand_pool(0 /* unused */, lmp)
#else #else
@ -134,6 +134,10 @@ void PairDPDKokkos<DeviceType>::compute(int eflagin, int vflagin)
special_lj[1] = force->special_lj[1]; special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2]; special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3]; special_lj[3] = force->special_lj[3];
special_rf[0] = sqrt(force->special_lj[0]);
special_rf[1] = sqrt(force->special_lj[1]);
special_rf[2] = sqrt(force->special_lj[2]);
special_rf[3] = sqrt(force->special_lj[3]);
nlocal = atom->nlocal; nlocal = atom->nlocal;
dtinvsqrt = 1.0/sqrt(update->dt); dtinvsqrt = 1.0/sqrt(update->dt);
@ -248,7 +252,7 @@ void PairDPDKokkos<DeviceType>::operator() (TagDPDKokkos<NEIGHFLAG,EVFLAG>, cons
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = d_neighbors(i,jj); j = d_neighbors(i,jj);
factor_dpd = special_lj[sbmask(j)]; factor_dpd = special_lj[sbmask(j)];
factor_sqrt = special_sqrt[sbmask(i)]; factor_sqrt = special_rf[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
delx = xtmp - x(j,0); delx = xtmp - x(j,0);

View File

@ -76,10 +76,10 @@ class PairDPDKokkos : public PairDPD {
template<int NEIGHFLAG> template<int NEIGHFLAG>
KOKKOS_INLINE_FUNCTION KOKKOS_INLINE_FUNCTION
void ev_tally(EV_FLOAT &ev, const int &i, const int &j, void ev_tally(EV_FLOAT &ev, const int &i, const int &j,
const F_FLOAT &epair, const F_FLOAT &fpair, const F_FLOAT &delx, const F_FLOAT &epair, const F_FLOAT &fpair,
const F_FLOAT &dely, const F_FLOAT &delz) const; const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const;
private: private:
double special_lj[4]; double special_lj[4], special_rf[4];
int eflag,vflag; int eflag,vflag;
int neighflag,nlocal; int neighflag,nlocal;
double dtinvsqrt; double dtinvsqrt;

View File

@ -41,8 +41,8 @@ using namespace LAMMPS_NS;
template<class DeviceType> template<class DeviceType>
PairDPDTstatKokkos<DeviceType>::PairDPDTstatKokkos(class LAMMPS *lmp) : PairDPDTstatKokkos<DeviceType>::PairDPDTstatKokkos(class LAMMPS *_lmp) :
PairDPDTstat(lmp) , PairDPDTstat(_lmp) ,
#ifdef DPD_USE_RAN_MARS #ifdef DPD_USE_RAN_MARS
rand_pool(0 /* unused */, lmp) rand_pool(0 /* unused */, lmp)
#else #else
@ -128,7 +128,6 @@ void PairDPDTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
memory->create(eatom,maxeatom,"pair:eatom"); memory->create(eatom,maxeatom,"pair:eatom");
memset(&eatom[0], 0, maxeatom * sizeof(double)); memset(&eatom[0], 0, maxeatom * sizeof(double));
} }
if (vflag_atom) { if (vflag_atom) {
memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->destroy_kokkos(k_vatom,vatom);
memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom");
@ -149,6 +148,10 @@ void PairDPDTstatKokkos<DeviceType>::compute(int eflagin, int vflagin)
special_lj[1] = force->special_lj[1]; special_lj[1] = force->special_lj[1];
special_lj[2] = force->special_lj[2]; special_lj[2] = force->special_lj[2];
special_lj[3] = force->special_lj[3]; special_lj[3] = force->special_lj[3];
special_rf[0] = sqrt(force->special_lj[0]);
special_rf[1] = sqrt(force->special_lj[1]);
special_rf[2] = sqrt(force->special_lj[2]);
special_rf[3] = sqrt(force->special_lj[3]);
nlocal = atom->nlocal; nlocal = atom->nlocal;
dtinvsqrt = 1.0/sqrt(update->dt); dtinvsqrt = 1.0/sqrt(update->dt);
@ -236,7 +239,7 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
int i,j,jj,jnum,itype,jtype; int i,j,jj,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double xtmp,ytmp,ztmp,delx,dely,delz,fpair;
double vxtmp,vytmp,vztmp,delvx,delvy,delvz; double vxtmp,vytmp,vztmp,delvx,delvy,delvz;
double rsq,r,rinv,dot,wd,randnum,factor_dpd; double rsq,r,rinv,dot,wd,randnum,factor_dpd,factor_sqrt;
double fx = 0,fy = 0,fz = 0; double fx = 0,fy = 0,fz = 0;
i = d_ilist[ii]; i = d_ilist[ii];
@ -252,6 +255,7 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
for (jj = 0; jj < jnum; jj++) { for (jj = 0; jj < jnum; jj++) {
j = d_neighbors(i,jj); j = d_neighbors(i,jj);
factor_dpd = special_lj[sbmask(j)]; factor_dpd = special_lj[sbmask(j)];
factor_sqrt = special_rf[sbmask(j)];
j &= NEIGHMASK; j &= NEIGHMASK;
delx = xtmp - x(j,0); delx = xtmp - x(j,0);
@ -274,10 +278,11 @@ void PairDPDTstatKokkos<DeviceType>::operator() (TagDPDTstatKokkos<NEIGHFLAG,VFL
// drag force - parallel // drag force - parallel
fpair = -params(itype,jtype).gamma*wd*wd*dot*rinv; fpair = -params(itype,jtype).gamma*wd*wd*dot*rinv;
fpair *= factor_dpd;
// random force - parallel // random force - parallel
fpair += params(itype,jtype).sigma*wd*randnum*dtinvsqrt; fpair += factor_sqrt*params(itype,jtype).sigma*wd*randnum*dtinvsqrt;
fpair *= factor_dpd*rinv; fpair *= rinv;
fx += fpair*delx; fx += fpair*delx;
fy += fpair*dely; fy += fpair*dely;

View File

@ -79,7 +79,7 @@ class PairDPDTstatKokkos : public PairDPDTstat {
const F_FLOAT &fpair, const F_FLOAT &delx, const F_FLOAT &fpair, const F_FLOAT &delx,
const F_FLOAT &dely, const F_FLOAT &delz) const; const F_FLOAT &dely, const F_FLOAT &delz) const;
private: private:
double special_lj[4]; double special_lj[4], special_rf[4];
int eflag,vflag; int eflag,vflag;
int neighflag,nlocal; int neighflag,nlocal;
double dtinvsqrt; double dtinvsqrt;