diff --git a/src/compute_pair.cpp b/src/compute_pair.cpp index 951df069ee..9617e016aa 100644 --- a/src/compute_pair.cpp +++ b/src/compute_pair.cpp @@ -21,34 +21,42 @@ using namespace LAMMPS_NS; +enum{EPAIR,EVDWL,ECOUL}; + /* ---------------------------------------------------------------------- */ ComputePair::ComputePair(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg != 4) error->all("Illegal compute pair command"); + if (narg < 4 || narg > 5) error->all("Illegal compute pair command"); if (igroup) error->all("Compute pair must use group all"); + scalar_flag = 1; + extscalar = 1; + peflag = 1; + timeflag = 1; + int n = strlen(arg[3]) + 1; pstyle = new char[n]; strcpy(pstyle,arg[3]); + if (narg == 5) { + if (strcmp(arg[4],"epair") == 0) evalue = EPAIR; + if (strcmp(arg[4],"evdwl") == 0) evalue = EVDWL; + if (strcmp(arg[4],"ecoul") == 0) evalue = ECOUL; + } else evalue = EPAIR; + pair = force->pair_match(pstyle,1); if (!pair) error->all("Unrecognized pair style in compute pair command"); npair = pair->nextra; - if (!npair) - error->all("Pair style in compute pair command stores no values"); - // settings - - vector_flag = 1; - size_vector = npair; - extvector = 1; - peflag = 1; - timeflag = 1; - - one = new double[npair]; - vector = new double[npair]; + if (npair) { + vector_flag = 1; + size_vector = npair; + extvector = 1; + one = new double[npair]; + vector = new double[npair]; + } else one = vector = NULL; } /* ---------------------------------------------------------------------- */ @@ -72,6 +80,23 @@ void ComputePair::init() /* ---------------------------------------------------------------------- */ +double ComputePair::compute_scalar() +{ + invoked_scalar = update->ntimestep; + if (update->eflag_global != invoked_scalar) + error->all("Energy was not tallied on needed timestep"); + + double eng; + if (evalue == EPAIR) eng = pair->eng_vdwl + pair->eng_coul; + else if (evalue == EVDWL) eng = pair->eng_vdwl; + else if (evalue == ECOUL) eng = pair->eng_coul; + + MPI_Allreduce(&eng,&scalar,1,MPI_DOUBLE,MPI_SUM,world); + return scalar; +} + +/* ---------------------------------------------------------------------- */ + void ComputePair::compute_vector() { invoked_vector = update->ntimestep; diff --git a/src/compute_pair.h b/src/compute_pair.h index 3d7d5dccfa..cf0e955ad9 100644 --- a/src/compute_pair.h +++ b/src/compute_pair.h @@ -29,10 +29,11 @@ class ComputePair : public Compute { ComputePair(class LAMMPS *, int, char **); ~ComputePair(); void init(); + double compute_scalar(); void compute_vector(); private: - int npair; + int evalue,npair; char *pstyle; class Pair *pair; double *one; diff --git a/src/pair.cpp b/src/pair.cpp index 5091d1f62f..d48a004997 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -426,6 +426,7 @@ void Pair::ev_tally(int i, int j, int nlocal, int newton_pair, void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, double delx, double dely, double delz) { + double v[6]; if (eflag_either) { if (eflag_global) { @@ -436,23 +437,30 @@ void Pair::ev_tally_full(int i, double evdwl, double ecoul, double fpair, } if (vflag_either) { - if (vflag_global) { - virial[0] += 0.5*delx*delx*fpair; - virial[1] += 0.5*dely*dely*fpair; - virial[2] += 0.5*delz*delz*fpair; - virial[3] += 0.5*delx*dely*fpair; - virial[4] += 0.5*delx*delz*fpair; - virial[5] += 0.5*dely*delz*fpair; - } - } + v[0] = 0.5*delx*delx*fpair; + v[1] = 0.5*dely*dely*fpair; + v[2] = 0.5*delz*delz*fpair; + v[3] = 0.5*delx*dely*fpair; + v[4] = 0.5*delx*delz*fpair; + v[5] = 0.5*dely*delz*fpair; - if (vflag_atom) { - vatom[i][0] += 0.5*delx*delx*fpair; - vatom[i][1] += 0.5*dely*dely*fpair; - vatom[i][2] += 0.5*delz*delz*fpair; - vatom[i][3] += 0.5*delx*dely*fpair; - vatom[i][4] += 0.5*delx*delz*fpair; - vatom[i][5] += 0.5*dely*delz*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) { + 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]; + } } } @@ -576,36 +584,36 @@ void Pair::ev_tally_xyz_full(int i, double evdwl, double ecoul, } if (vflag_either) { - v[0] = delx*fx; - v[1] = dely*fy; - v[2] = delz*fz; - v[3] = delx*fy; - v[4] = delx*fz; - v[5] = dely*fz; - } + v[0] = 0.5*delx*fx; + v[1] = 0.5*dely*fy; + v[2] = 0.5*delz*fz; + v[3] = 0.5*delx*fy; + v[4] = 0.5*delx*fz; + v[5] = 0.5*dely*fz; - if (vflag_global) { - virial[0] += 0.5*v[0]; - virial[1] += 0.5*v[1]; - virial[2] += 0.5*v[2]; - virial[3] += 0.5*v[3]; - virial[4] += 0.5*v[4]; - virial[5] += 0.5*v[5]; - } - - if (vflag_atom) { - vatom[i][0] += 0.5*v[0]; - vatom[i][1] += 0.5*v[1]; - vatom[i][2] += 0.5*v[2]; - vatom[i][3] += 0.5*v[3]; - vatom[i][4] += 0.5*v[4]; - vatom[i][5] += 0.5*v[5]; + 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) { + 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]; + } } } /* ---------------------------------------------------------------------- tally eng_vdwl and virial into global and per-atom accumulators - called by SW potential, newton_pair is always on + called by SW and hbond potentials, newton_pair is always on virial = riFi + rjFj + rkFk = (rj-ri) Fj + (rk-ri) Fk = drji*fj + drki*fk ------------------------------------------------------------------------- */ @@ -627,20 +635,36 @@ void Pair::ev_tally3(int i, int j, int k, double evdwl, double ecoul, } } - if (vflag_atom) { - 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]); + if (vflag_either) { + v[0] = drji[0]*fj[0] + drki[0]*fk[0]; + v[1] = drji[1]*fj[1] + drki[1]*fk[1]; + v[2] = drji[2]*fj[2] + drki[2]*fk[2]; + v[3] = drji[0]*fj[1] + drki[0]*fk[1]; + v[4] = drji[0]*fj[2] + drki[0]*fk[2]; + v[5] = drji[1]*fj[2] + drki[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]; + } - 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]; + if (vflag_atom) { + vatom[i][0] += THIRD*v[0]; vatom[i][1] += THIRD*v[1]; + vatom[i][2] += THIRD*v[2]; vatom[i][3] += THIRD*v[3]; + vatom[i][4] += THIRD*v[4]; vatom[i][5] += THIRD*v[5]; + + vatom[j][0] += THIRD*v[0]; vatom[j][1] += THIRD*v[1]; + vatom[j][2] += THIRD*v[2]; vatom[j][3] += THIRD*v[3]; + vatom[j][4] += THIRD*v[4]; vatom[j][5] += THIRD*v[5]; + + vatom[k][0] += THIRD*v[0]; vatom[k][1] += THIRD*v[1]; + vatom[k][2] += THIRD*v[2]; vatom[k][3] += THIRD*v[3]; + vatom[k][4] += THIRD*v[4]; vatom[k][5] += THIRD*v[5]; + } } }