diff --git a/doc/utils/converters/lammpsdoc/txt2rst.py b/doc/utils/converters/lammpsdoc/txt2rst.py index 6949cb690e..9fac2f0635 100755 --- a/doc/utils/converters/lammpsdoc/txt2rst.py +++ b/doc/utils/converters/lammpsdoc/txt2rst.py @@ -104,7 +104,7 @@ class RSTMarkup(Markup): anchor_pos = href.find('#') - if anchor_pos >= 0: + if anchor_pos >= 0 and not href.startswith('http'): href = href[anchor_pos+1:] return ":ref:`%s <%s>`" % (content, href) diff --git a/doc/utils/converters/tests/test_txt2rst.py b/doc/utils/converters/tests/test_txt2rst.py index f74f8b70dd..9bdb5b2e48 100644 --- a/doc/utils/converters/tests/test_txt2rst.py +++ b/doc/utils/converters/tests/test_txt2rst.py @@ -424,6 +424,11 @@ class TestSpecialCommands(unittest.TestCase): "one \n\n" "a :ref:`link ` to above\n\n", s) + def test_external_anchor_link(self): + s = self.txt2rst.convert('some text "containing a\n' + 'link"_http://lammps.sandia.gov/movies.html#granregion with an anchor') + self.assertEqual('some text `containing a link `_ with an anchor\n\n', s) + def test_define_link_alias(self): s = self.txt2rst.convert("one :link(alias,value)\n" "\"test\"_alias\n") diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 050d4a2b96..ee5db06e30 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1557,5 +1557,5 @@ size_t const colvarmodule::cv_prec = 14; size_t const colvarmodule::cv_width = 21; size_t const colvarmodule::en_prec = 14; size_t const colvarmodule::en_width = 21; -std::string const colvarmodule::line_marker = +const char * const colvarmodule::line_marker = (const char *) "----------------------------------------------------------------------\n"; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 454fb65b6c..0b19d1f779 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -366,7 +366,7 @@ public: /// Number of characters to represent the collective variables energy static size_t const en_width; /// Line separator in the log output - static std::string const line_marker; + static const char * const line_marker; // proxy functions diff --git a/src/CLASS2/improper_class2.cpp b/src/CLASS2/improper_class2.cpp index 942f46aa22..6c599276ca 100644 --- a/src/CLASS2/improper_class2.cpp +++ b/src/CLASS2/improper_class2.cpp @@ -819,7 +819,8 @@ void ImproperClass2::angleangle(int eflag, int vflag) if (evflag) ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper, fabcd[0],fabcd[2],fabcd[3], - delxAB,delyAB,delzAB,delxBC,delyBC,delzBC,delxBD,delyBD,delzBD); + delxAB,delyAB,delzAB,delxBC,delyBC,delzBC, + delxBD-delxBC,delyBD-delyBC,delzBD-delzBC); } } diff --git a/src/GRANULAR/fix_wall_gran.cpp b/src/GRANULAR/fix_wall_gran.cpp index 781ec91447..8ba6a2402d 100644 --- a/src/GRANULAR/fix_wall_gran.cpp +++ b/src/GRANULAR/fix_wall_gran.cpp @@ -37,7 +37,7 @@ using namespace MathConst; // XYZ PLANE need to be 0,1,2 -enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; +enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY}; enum{NONE,CONSTANT,EQUAL}; @@ -78,21 +78,21 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : kn = force->numeric(FLERR,arg[4]); if (strcmp(arg[5],"NULL") == 0) kt = kn * 2.0/7.0; else kt = force->numeric(FLERR,arg[5]); - + gamman = force->numeric(FLERR,arg[6]); if (strcmp(arg[7],"NULL") == 0) gammat = 0.5 * gamman; else gammat = force->numeric(FLERR,arg[7]); - + xmu = force->numeric(FLERR,arg[8]); int dampflag = force->inumeric(FLERR,arg[9]); if (dampflag == 0) gammat = 0.0; - + if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 || xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1) error->all(FLERR,"Illegal fix wall/gran command"); - + // convert Kn and Kt from pressure units to force/distance^2 if Hertzian - + if (pairstyle == HERTZ_HISTORY) { kn /= force->nktv2p; kt /= force->nktv2p; @@ -119,7 +119,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : } // wallstyle args - + idregion = NULL; if (strcmp(arg[iarg],"xplane") == 0) { @@ -160,12 +160,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : strcpy(idregion,arg[iarg+1]); iarg += 2; } - + // optional args wiggle = 0; wshear = 0; - + while (iarg < narg) { if (strcmp(arg[iarg],"wiggle") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix wall/gran command"); @@ -188,7 +188,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : iarg += 3; } else error->all(FLERR,"Illegal fix wall/gran command"); } - + if (wallstyle == XPLANE && domain->xperiodic) error->all(FLERR,"Cannot use wall in periodic dimension"); if (wallstyle == YPLANE && domain->yperiodic) @@ -197,7 +197,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Cannot use wall in periodic dimension"); if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic)) error->all(FLERR,"Cannot use wall in periodic dimension"); - + if (wiggle && wshear) error->all(FLERR,"Cannot wiggle and shear fix wall/gran"); if (wiggle && wallstyle == ZCYLINDER && axis != 2) @@ -212,15 +212,15 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Cannot wiggle or shear with fix wall/gran/region"); // setup oscillations - + if (wiggle) omega = 2.0*MY_PI / period; - + // perform initial allocation of atom-based arrays // register with Atom class if (pairstyle == BONDED_HISTORY) sheardim = 7; else sheardim = 3; - + shearone = NULL; grow_arrays(atom->nmax); atom->add_callback(0); @@ -348,7 +348,7 @@ void FixWallGran::post_force(int vflag) } vwall[axis] = amplitude*omega*sin(arg); } else if (wshear) vwall[axis] = vshear; - + // loop over all my atoms // rsq = distance from wall // dx,dy,dz = signed distance from wall @@ -358,7 +358,7 @@ void FixWallGran::post_force(int vflag) // compute force and torque on atom if close enough to wall // via wall potential matched to pair potential // set shear if pair potential stores history - + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -368,14 +368,14 @@ void FixWallGran::post_force(int vflag) double *rmass = atom->rmass; int *mask = atom->mask; int nlocal = atom->nlocal; - + rwall = 0.0; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - + dx = dy = dz = 0.0; - + if (wallstyle == XPLANE) { del1 = x[i][0] - wlo; del2 = whi - x[i][0]; @@ -401,7 +401,7 @@ void FixWallGran::post_force(int vflag) dx = -delr/delxy * x[i][0]; dy = -delr/delxy * x[i][1]; // rwall = -2r_c if inside cylinder, 2r_c outside - rwall = 2*(1-2*(delxy < cylradius))*cylradius; + rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius; if (wshear && axis != 2) { vwall[0] += vshear * x[i][1]/delxy; vwall[1] += -vshear * x[i][0]/delxy; @@ -409,12 +409,12 @@ void FixWallGran::post_force(int vflag) } } } - + rsq = dx*dx + dy*dy + dz*dz; - + if (rsq > radius[i]*radius[i]) { if (history) - for (j = 0; j < sheardim; j++) + for (j = 0; j < sheardim; j++) shearone[i][j] = 0.0; } else { @@ -467,66 +467,66 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz, rsqinv = 1.0/rsq; // relative translational velocity - + vr1 = v[0] - vwall[0]; vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - + // normal component - + vnnr = vr1*dx + vr2*dy + vr3*dz; vn1 = dx*vnnr * rsqinv; vn2 = dy*vnnr * rsqinv; vn3 = dz*vnnr * rsqinv; - + // tangential component - + vt1 = vr1 - vn1; vt2 = vr2 - vn2; vt3 = vr3 - vn3; - + // relative rotational velocity - + wr1 = radius*omega[0] * rinv; wr2 = radius*omega[1] * rinv; wr3 = radius*omega[2] * rinv; - + // normal forces = Hookian contact + normal velocity damping - + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - + // relative velocities - + vtr1 = vt1 - (dz*wr2-dy*wr3); vtr2 = vt2 - (dx*wr3-dz*wr1); vtr3 = vt3 - (dy*wr1-dx*wr2); vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; vrel = sqrt(vrel); - + // force normalization - + fn = xmu * fabs(ccel*r); fs = meff*gammat*vrel; if (vrel != 0.0) ft = MIN(fn,fs) / vrel; else ft = 0.0; - + // tangential force due to tangential velocity damping - + fs1 = -ft*vtr1; fs2 = -ft*vtr2; fs3 = -ft*vtr3; - + // forces & torques - + fx = dx*ccel + fs1; fy = dy*ccel + fs2; fz = dz*ccel + fs3; - + f[0] += fx; f[1] += fy; f[2] += fz; - + tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -546,60 +546,60 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel; double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; double shrmag,rsht,rinv,rsqinv; - + r = sqrt(rsq); rinv = 1.0/r; rsqinv = 1.0/rsq; - + // relative translational velocity - + vr1 = v[0] - vwall[0]; vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - + // normal component - + vnnr = vr1*dx + vr2*dy + vr3*dz; vn1 = dx*vnnr * rsqinv; vn2 = dy*vnnr * rsqinv; vn3 = dz*vnnr * rsqinv; - + // tangential component - + vt1 = vr1 - vn1; vt2 = vr2 - vn2; vt3 = vr3 - vn3; - + // relative rotational velocity - + wr1 = radius*omega[0] * rinv; wr2 = radius*omega[1] * rinv; wr3 = radius*omega[2] * rinv; - + // normal forces = Hookian contact + normal velocity damping - + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; - + // relative velocities - + vtr1 = vt1 - (dz*wr2-dy*wr3); vtr2 = vt2 - (dx*wr3-dz*wr1); vtr3 = vt3 - (dy*wr1-dx*wr2); vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; vrel = sqrt(vrel); - + // shear history effects - + if (shearupdate) { shear[0] += vtr1*dt; shear[1] += vtr2*dt; shear[2] += vtr3*dt; } shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - + // rotate shear displacements - + rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; rsht = rsht*rsqinv; if (shearupdate) { @@ -607,18 +607,18 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, shear[1] -= rsht*dy; shear[2] -= rsht*dz; } - + // tangential forces = shear + tangential velocity damping - + fs1 = - (kt*shear[0] + meff*gammat*vtr1); fs2 = - (kt*shear[1] + meff*gammat*vtr2); fs3 = - (kt*shear[2] + meff*gammat*vtr3); - + // rescale frictional displacements and forces if needed - + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fn = xmu * fabs(ccel*r); - + if (fs > fn) { if (shrmag != 0.0) { shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - @@ -632,17 +632,17 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz, fs3 *= fn/fs; } else fs1 = fs2 = fs3 = 0.0; } - + // forces & torques - + fx = dx*ccel + fs1; fy = dy*ccel + fs2; fz = dz*ccel + fs3; - + f[0] += fx; f[1] += fy; f[2] += fz; - + tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -662,56 +662,56 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, double wr1,wr2,wr3,damp,ccel,vtr1,vtr2,vtr3,vrel; double fn,fs,fs1,fs2,fs3,fx,fy,fz,tor1,tor2,tor3; double shrmag,rsht,polyhertz,rinv,rsqinv; - + r = sqrt(rsq); rinv = 1.0/r; rsqinv = 1.0/rsq; - + // relative translational velocity - + vr1 = v[0] - vwall[0]; vr2 = v[1] - vwall[1]; vr3 = v[2] - vwall[2]; - + // normal component - + vnnr = vr1*dx + vr2*dy + vr3*dz; vn1 = dx*vnnr / rsq; vn2 = dy*vnnr / rsq; vn3 = dz*vnnr / rsq; - + // tangential component - + vt1 = vr1 - vn1; vt2 = vr2 - vn2; vt3 = vr3 - vn3; - + // relative rotational velocity - + wr1 = radius*omega[0] * rinv; wr2 = radius*omega[1] * rinv; wr3 = radius*omega[2] * rinv; - + // normal forces = Hertzian contact + normal velocity damping // rwall = 0 is flat wall case // rwall positive or negative is curved wall - // will break (as it should) if rwall is negative and + // will break (as it should) if rwall is negative and // its absolute value < radius of particle - + damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radius-r)*rinv - damp; if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius); else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius)); ccel *= polyhertz; - + // relative velocities - + vtr1 = vt1 - (dz*wr2-dy*wr3); vtr2 = vt2 - (dx*wr3-dz*wr1); vtr3 = vt3 - (dy*wr1-dx*wr2); vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; vrel = sqrt(vrel); - + // shear history effects if (shearupdate) { @@ -720,9 +720,9 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, shear[2] += vtr3*dt; } shrmag = sqrt(shear[0]*shear[0] + shear[1]*shear[1] + shear[2]*shear[2]); - + // rotate shear displacements - + rsht = shear[0]*dx + shear[1]*dy + shear[2]*dz; rsht = rsht*rsqinv; if (shearupdate) { @@ -730,18 +730,18 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, shear[1] -= rsht*dy; shear[2] -= rsht*dz; } - + // tangential forces = shear + tangential velocity damping - + fs1 = -polyhertz * (kt*shear[0] + meff*gammat*vtr1); fs2 = -polyhertz * (kt*shear[1] + meff*gammat*vtr2); fs3 = -polyhertz * (kt*shear[2] + meff*gammat*vtr3); - + // rescale frictional displacements and forces if needed - + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); fn = xmu * fabs(ccel*r); - + if (fs > fn) { if (shrmag != 0.0) { shear[0] = (fn/fs) * (shear[0] + meff*gammat*vtr1/kt) - @@ -757,7 +757,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, } // forces & torques - + fx = dx*ccel + fs1; fy = dy*ccel + fs2; fz = dz*ccel + fs3; @@ -765,7 +765,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz, f[0] += fx; f[1] += fy; f[2] += fz; - + tor1 = rinv * (dy*fs3 - dz*fs2); tor2 = rinv * (dz*fs1 - dx*fs3); tor3 = rinv * (dx*fs2 - dy*fs1); @@ -864,8 +864,8 @@ void FixWallGran::bonded_history(double rsq, double dx, double dy, double dz, // use Tsuji et al form - polyhertz = 1.2728- 4.2783*0.9 + 11.087*0.9*0.9 - 22.348*0.9*0.9*0.9 + - 27.467*0.9*0.9*0.9*0.9 - 18.022*0.9*0.9*0.9*0.9*0.9 + + polyhertz = 1.2728- 4.2783*0.9 + 11.087*0.9*0.9 - 22.348*0.9*0.9*0.9 + + 27.467*0.9*0.9*0.9*0.9 - 18.022*0.9*0.9*0.9*0.9*0.9 + 4.8218*0.9*0.9*0.9*0.9*0.9*0.9; gammatsuji = 0.2*sqrt(meff*kn); @@ -969,7 +969,7 @@ void FixWallGran::bonded_history(double rsq, double dx, double dy, double dz, shear[6] += magtwist*dt; ktwist = 0.5*kt*(a0*aovera0)*(a0*aovera0); - magtortwist = -ktwist*shear[6] - + magtortwist = -ktwist*shear[6] - 0.5*polyhertz*gammatsuji*(a0*aovera0)*(a0*aovera0)*magtwist; twistcrit=TWOTHIRDS*a0*aovera0*Fcrit; @@ -1113,11 +1113,11 @@ void FixWallGran::unpack_restart(int nlocal, int nth) double **extra = atom->extra; // skip to Nth set of extra values - + int m = 0; for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); m++; - + for (int i = 0; i < sheardim; i++) shearone[nlocal][i] = extra[nlocal][m++]; } diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp index 30a62ffb18..a2f176e732 100644 --- a/src/MOLECULE/improper_umbrella.cpp +++ b/src/MOLECULE/improper_umbrella.cpp @@ -231,9 +231,25 @@ void ImproperUmbrella::compute(int eflag, int vflag) f[i4][2] += f4[2]*a; } - if (evflag) + if (evflag) { + + // get correct 4-body geometry for virial tally + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + ev_tally(i1,i2,i3,i4,nlocal,newton_bond,eimproper,f1,f3,f4, vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z); + } } } diff --git a/src/Makefile b/src/Makefile index 92c0029dc0..12466e01fa 100755 --- a/src/Makefile +++ b/src/Makefile @@ -55,8 +55,8 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \ user-quip user-reaxc user-smd user-smtbq user-sph user-tally \ user-vtk -PACKLIB = compress gpu kim kokkos meam poems python reax voronoi \ - user-atc user-awpmd user-colvars user-h5md user-molfile \ +PACKLIB = compress gpu kim kokkos meam mpiio poems python reax voronoi \ + user-atc user-awpmd user-colvars user-h5md user-lb user-molfile \ user-qmmm user-quip user-vtk PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index 5e63d07c3b..dc97db1942 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -193,7 +193,7 @@ double colvarproxy_lammps::compute() previous_step = _lmp->update->ntimestep; if (cvm::debug()) { - cvm::log(cvm::line_marker+ + cvm::log(std::string(cvm::line_marker)+ "colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+ "Updating internal data.\n"); } diff --git a/src/USER-MISC/fix_flow_gauss.cpp b/src/USER-MISC/fix_flow_gauss.cpp index 129a4aa0e2..ad4c78f87f 100644 --- a/src/USER-MISC/fix_flow_gauss.cpp +++ b/src/USER-MISC/fix_flow_gauss.cpp @@ -27,6 +27,7 @@ #include "domain.h" #include "error.h" #include "citeme.h" +#include "respa.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -63,6 +64,8 @@ FixFlowGauss::FixFlowGauss(LAMMPS *lmp, int narg, char **arg) : extvector = 1; size_vector = 3; global_freq = 1; //data available every timestep + respa_level_support = 1; + //default respa level=outermost level is set in init() dimension = domain->dimension; @@ -109,9 +112,23 @@ int FixFlowGauss::setmask() int mask = 0; mask |= POST_FORCE; mask |= THERMO_ENERGY; + mask |= POST_FORCE_RESPA; return mask; } +/* ---------------------------------------------------------------------- */ + +void FixFlowGauss::init() +{ + //if respa level specified by fix_modify, then override default (outermost) + //if specified level too high, set to max level + if (strstr(update->integrate_style,"respa")) { + ilevel_respa = ((Respa *) update->integrate)->nlevels-1; + if (respa_level >= 0) + ilevel_respa = MIN(respa_level,ilevel_respa); + } +} + /* ---------------------------------------------------------------------- setup is called after the initial evaluation of forces before a run, so we must remove the total force here too @@ -127,7 +144,13 @@ void FixFlowGauss::setup(int vflag) if (mTot <= 0.0) error->all(FLERR,"Invalid group mass in fix flow/gauss"); - post_force(vflag); + if (strstr(update->integrate_style,"respa")) { + ((Respa *) update->integrate)->copy_flevel_f(ilevel_respa); + post_force_respa(vflag,ilevel_respa,0); + ((Respa *) update->integrate)->copy_f_flevel(ilevel_respa); + } + else + post_force(vflag); } /* ---------------------------------------------------------------------- @@ -148,7 +171,6 @@ void FixFlowGauss::post_force(int vflag) int ii,jj; //find the total force on all atoms - //initialize to zero double f_thisProc[3]; for (ii=0; ii<3; ii++) @@ -201,6 +223,11 @@ void FixFlowGauss::post_force(int vflag) } +void FixFlowGauss::post_force_respa(int vflag, int ilevel, int iloop) +{ + if (ilevel == ilevel_respa) post_force(vflag); +} + /* ---------------------------------------------------------------------- negative of work done by this fix This is only computed if requested, either with fix_modify energy yes, or with the energy keyword. Otherwise returns 0. diff --git a/src/USER-MISC/fix_flow_gauss.h b/src/USER-MISC/fix_flow_gauss.h index 5337679542..75bdf8c1f9 100644 --- a/src/USER-MISC/fix_flow_gauss.h +++ b/src/USER-MISC/fix_flow_gauss.h @@ -30,10 +30,12 @@ FixStyle(flow/gauss,FixFlowGauss) public: FixFlowGauss(class LAMMPS *, int, char **); int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); double compute_scalar(); double compute_vector(int n); - void post_force(int); - void setup(int); protected: int dimension; @@ -44,6 +46,7 @@ FixStyle(flow/gauss,FixFlowGauss) double pe_tot; //total added energy double dt; //timestep bool workflag; //if calculate work done by fix + int ilevel_respa; }; diff --git a/src/USER-OMP/fix_wall_gran_omp.cpp b/src/USER-OMP/fix_wall_gran_omp.cpp index f80d8040bb..55ae8593a7 100644 --- a/src/USER-OMP/fix_wall_gran_omp.cpp +++ b/src/USER-OMP/fix_wall_gran_omp.cpp @@ -25,8 +25,8 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER}; // XYZ PLANE need to be 0,1,2 -enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY}; +enum{XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION}; // XYZ PLANE need to be 0,1,2 +enum{HOOKE,HOOKE_HISTORY,HERTZ_HISTORY,BONDED_HISTORY}; #define BIG 1.0e20 @@ -107,6 +107,7 @@ void FixWallGranOMP::post_force(int vflag) if (mask[i] & groupbit) { double dx,dy,dz,del1,del2,delxy,delr,rsq; + double rwall = 0.0; dx = dy = dz = 0.0; @@ -128,13 +129,16 @@ void FixWallGranOMP::post_force(int vflag) } else if (wallstyle == ZCYLINDER) { delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]); delr = cylradius - delxy; - if (delr > radius[i]) dz = cylradius; - else { + if (delr > radius[i]) { + dz = cylradius; + rwall = 0.0; + } else { dx = -delr/delxy * x[i][0]; dy = -delr/delxy * x[i][1]; + rwall = (delxy < cylradius) ? -2*cylradius : 2*cylradius; if (wshear && axis != 2) { - vwall[0] = vshear * x[i][1]/delxy; - vwall[1] = -vshear * x[i][0]/delxy; + vwall[0] += vshear * x[i][1]/delxy; + vwall[1] += -vshear * x[i][0]/delxy; vwall[2] = 0.0; } } @@ -143,11 +147,10 @@ void FixWallGranOMP::post_force(int vflag) rsq = dx*dx + dy*dy + dz*dz; if (rsq > radius[i]*radius[i]) { - if (pairstyle != HOOKE) { - shear[i][0] = 0.0; - shear[i][1] = 0.0; - shear[i][2] = 0.0; - } + if (history) + for (int j = 0; j < sheardim; j++) + shearone[i][j] = 0.0; + } else { // meff = effective mass of sphere @@ -156,17 +159,20 @@ void FixWallGranOMP::post_force(int vflag) double meff = rmass[i]; if (fix_rigid && mass_rigid[i] > 0.0) meff = mass_rigid[i]; - // inovke sphere/wall interaction + // invoke sphere/wall interaction if (pairstyle == HOOKE) - hooke(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],meff); + hooke(rsq,dx,dy,dz,vwall,v[i],f[i], + omega[i],torque[i],radius[i],meff); else if (pairstyle == HOOKE_HISTORY) - hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],meff,shear[i]); + hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i], + omega[i],torque[i],radius[i],meff,shearone[i]); else if (pairstyle == HERTZ_HISTORY) - hertz_history(rsq,dx,dy,dz,vwall,v[i],f[i],omega[i],torque[i], - radius[i],meff,shear[i]); + hertz_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], + omega[i],torque[i],radius[i],meff,shearone[i]); + else if (pairstyle == BONDED_HISTORY) + bonded_history(rsq,dx,dy,dz,vwall,rwall,v[i],f[i], + omega[i],torque[i],radius[i],meff,shearone[i]); } } } diff --git a/src/USER-OMP/improper_class2_omp.cpp b/src/USER-OMP/improper_class2_omp.cpp index 5e32795371..3beeb41a26 100644 --- a/src/USER-OMP/improper_class2_omp.cpp +++ b/src/USER-OMP/improper_class2_omp.cpp @@ -696,7 +696,8 @@ void ImproperClass2OMP::angleangle_thr(int nfrom, int nto, ThrData * const thr) if (EVFLAG) ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper, - fabcd[0],fabcd[2],fabcd[3],delxAB,delyAB,delzAB, - delxBC,delyBC,delzBC,delxBD,delyBD,delzBD,thr); + fabcd[0],fabcd[2],fabcd[3], + delxAB,delyAB,delzAB,delxBC,delyBC,delzBC, + delxBD-delxBC,delyBD-delyBC,delzBD-delzBC,thr); } } diff --git a/src/USER-OMP/improper_umbrella_omp.cpp b/src/USER-OMP/improper_umbrella_omp.cpp index 9d28e25053..9ea94c5705 100644 --- a/src/USER-OMP/improper_umbrella_omp.cpp +++ b/src/USER-OMP/improper_umbrella_omp.cpp @@ -254,8 +254,24 @@ void ImproperUmbrellaOMP::eval(int nfrom, int nto, ThrData * const thr) f[i4].z += f4[2]*a; } - if (EVFLAG) + if (EVFLAG) { + + // get correct 4-body geometry for virial tally + + vb1x = x[i1].x - x[i2].x; + vb1y = x[i1].y - x[i2].y; + vb1z = x[i1].z - x[i2].z; + + vb2x = x[i3].x - x[i2].x; + vb2y = x[i3].y - x[i2].y; + vb2z = x[i3].z - x[i2].z; + + vb3x = x[i4].x - x[i3].x; + vb3y = x[i4].y - x[i3].y; + vb3z = x[i4].z - x[i3].z; + ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,eimproper,f1,f3,f4, vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,thr); + } } }