diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp index bb750a2ada..fdc3029609 100644 --- a/src/compute_reduce.cpp +++ b/src/compute_reduce.cpp @@ -195,7 +195,8 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : "calculate a per-atom array"); if (argindex[i] && argindex[i] > modify->compute[icompute]->size_peratom_cols) - error->all(FLERR,"Compute reduce compute array is accessed out-of-range"); + error->all(FLERR, + "Compute reduce compute array is accessed out-of-range"); } else if (modify->compute[icompute]->local_flag) { flavor[i] = LOCAL; if (argindex[i] == 0 && @@ -207,8 +208,10 @@ ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) : "calculate a local array"); if (argindex[i] && argindex[i] > modify->compute[icompute]->size_local_cols) - error->all(FLERR,"Compute reduce compute array is accessed out-of-range"); - } else error->all(FLERR,"Compute reduce compute calculates global values"); + error->all(FLERR, + "Compute reduce compute array is accessed out-of-range"); + } else error->all(FLERR, + "Compute reduce compute calculates global values"); } else if (which[i] == FIX) { int ifix = modify->find_fix(ids[i]); diff --git a/src/respa.cpp b/src/respa.cpp index 19a6f7d800..339c965d1a 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -250,14 +250,6 @@ void Respa::init() if (modify->nfix == 0 && comm->me == 0) error->warning(FLERR,"No fixes defined, atoms won't move"); - // warn about incorrect pressures when using rRESPA with fix SHAKE - - int shakeflag = 0; - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style,"shake") == 0) shakeflag = 1; - if (shakeflag && comm->me == 0) - error->warning(FLERR,"Fix shake with rRESPA computes invalid pressures"); - // create fix needed for storing atom-based respa level forces // will delete it at end of run diff --git a/src/velocity.cpp b/src/velocity.cpp index a5c6a8673e..701a4f31eb 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -26,6 +26,9 @@ #include "variable.h" #include "force.h" #include "modify.h" +#include "fix.h" +#include "fix_rigid.h" +#include "fix_rigid_small.h" #include "compute.h" #include "compute_temp.h" #include "random_park.h" @@ -86,6 +89,7 @@ void Velocity::command(int narg, char **arg) rotation_flag = 0; loop_flag = ALL; scale_flag = 1; + rfix = -1; // read options from end of input line // change defaults as options specify @@ -617,13 +621,41 @@ void Velocity::ramp(int narg, char **arg) /* ---------------------------------------------------------------------- zero linear or angular momentum of a group + if using rigid/small requires init of entire system since + its methods perform forward/reverse comm, + comm::init needs neighbor::init needs pair::init needs kspace::init, etc + also requires setup_pre_neighbor call to setup bodies ------------------------------------------------------------------------- */ void Velocity::zero(int narg, char **arg) { - if (strcmp(arg[0],"linear") == 0) zero_momentum(); - else if (strcmp(arg[0],"angular") == 0) zero_rotation(); - else error->all(FLERR,"Illegal velocity command"); + if (strcmp(arg[0],"linear") == 0) { + if (rfix < 0) zero_momentum(); + else { + if (strcmp(modify->fix[rfix]->style,"rigid") == 0) + ((FixRigid *) modify->fix[rfix])->zero_momentum(igroup); + else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) { + lmp->init(); + ((FixRigidSmall *) modify->fix[rfix])->setup_pre_neighbor(); + ((FixRigidSmall *) modify->fix[rfix])->zero_momentum(igroup); + } + else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); + } + + } else if (strcmp(arg[0],"angular") == 0) { + if (rfix < 0) zero_rotation(); + else { + if (strcmp(modify->fix[rfix]->style,"rigid") == 0) + ((FixRigid *) modify->fix[rfix])->zero_rotation(igroup); + else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) { + lmp->init(); + ((FixRigidSmall *) modify->fix[rfix])->setup_pre_neighbor(); + ((FixRigidSmall *) modify->fix[rfix])->zero_rotation(igroup); + } + else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); + } + + } else error->all(FLERR,"Illegal velocity command"); } /* ---------------------------------------------------------------------- @@ -654,10 +686,10 @@ void Velocity::rescale(double t_old, double t_new) void Velocity::zero_momentum() { - // cannot have 0 atoms in group + // cannot have no atoms in group if (group->count(igroup) == 0) - error->all(FLERR,"Cannot zero momentum of 0 atoms"); + error->all(FLERR,"Cannot zero momentum of no atoms"); // compute velocity of center-of-mass of group @@ -687,10 +719,10 @@ void Velocity::zero_rotation() { int i; - // cannot have 0 atoms in group + // cannot have no atoms in group if (group->count(igroup) == 0) - error->all(FLERR,"Cannot zero momentum of 0 atoms"); + error->all(FLERR,"Cannot zero momentum of no atoms"); // compute omega (angular velocity) of group around center-of-mass @@ -717,7 +749,7 @@ void Velocity::zero_rotation() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { domain->unmap(x[i],image[i],unwrap); - dx = unwrap[0]- xcm[0]; + dx = unwrap[0] - xcm[0]; dy = unwrap[1] - xcm[1]; dz = unwrap[2] - xcm[2]; v[i][0] -= omega[1]*dz - omega[2]*dy; @@ -779,6 +811,11 @@ void Velocity::options(int narg, char **arg) else if (strcmp(arg[iarg+1],"geom") == 0) loop_flag = GEOM; else error->all(FLERR,"Illegal velocity command"); iarg += 2; + } else if (strcmp(arg[iarg],"rigid") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); + rfix = modify->find_fix(arg[iarg+1]); + if (rfix < 0) error->all(FLERR,"Fix ID for velocity does not exist"); + iarg += 2; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); if (strcmp(arg[iarg+1],"box") == 0) scale_flag = 0; diff --git a/src/velocity.h b/src/velocity.h index 076dc1d343..d5e175bbc9 100644 --- a/src/velocity.h +++ b/src/velocity.h @@ -35,7 +35,7 @@ class Velocity : protected Pointers { private: int igroup,groupbit; int style; - int dist_flag,sum_flag,momentum_flag,rotation_flag,loop_flag,scale_flag; + int dist_flag,sum_flag,momentum_flag,rotation_flag,loop_flag,scale_flag,rfix; double xscale,yscale,zscale; class Compute *temperature;