diff --git a/src/PRD/prd.cpp b/src/PRD/prd.cpp index 81ec477514..27d8d327db 100644 --- a/src/PRD/prd.cpp +++ b/src/PRD/prd.cpp @@ -67,7 +67,7 @@ void PRD::command(int narg, char **arg) error->all("PRD command before simulation box is defined"); if (universe->nworlds != universe->nprocs && atom->map_style == 0) - error->all("Cannot use PRD with multi-proc replicas " + error->all("Cannot use PRD with multi-processor replicas " "unless atom map exists"); if (universe->nworlds == 1 && comm->me == 0) error->warning("Running PRD with only one replica"); diff --git a/src/fix_heat.cpp b/src/fix_heat.cpp index 5e8f75e6c5..19f7aea004 100644 --- a/src/fix_heat.cpp +++ b/src/fix_heat.cpp @@ -21,6 +21,7 @@ #include "fix_heat.h" #include "atom.h" #include "domain.h" +#include "region.h" #include "group.h" #include "force.h" #include "update.h" @@ -43,6 +44,20 @@ FixHeat::FixHeat(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) heat_input = atof(arg[4]); + // optional args + + iregion = -1; + + int iarg = 6; + while (iarg < narg) { + if (strcmp(arg[iarg],"region") == 0) { + if (iarg+2 > narg) error->all("Illegal fix heat command"); + iregion = domain->find_region(arg[iarg+1]); + if (iregion == -1) error->all("Fix heat region ID does not exist"); + iarg += 2; + } else error->all("Illegal fix heat command"); + } + scale = 1.0; } @@ -69,15 +84,23 @@ void FixHeat::init() void FixHeat::end_of_step() { - double **v = atom->v; - int *mask = atom->mask; - int nlocal = atom->nlocal; - + double heat,ke; double vsub[3],vcm[3]; + Region *region = NULL; + if (iregion >= 0) region = domain->regions[iregion]; + + if (region < 0) { + heat = heat_input*nevery*update->dt*force->ftm2v; + ke = group->ke(igroup)*force->ftm2v; + group->vcm(igroup,masstotal,vcm); + } else { + masstotal = group->mass(igroup,iregion); + if (masstotal == 0.0) error->all("Fix heat group has no atoms"); + heat = heat_input*nevery*update->dt*force->ftm2v; + ke = group->ke(igroup,iregion)*force->ftm2v; + group->vcm(igroup,masstotal,vcm,iregion); + } - double heat = heat_input*nevery*update->dt*force->ftm2v; - double ke = group->ke(igroup)*force->ftm2v; - group->vcm(igroup,masstotal,vcm); double vcmsq = vcm[0]*vcm[0] + vcm[1]*vcm[1] + vcm[2]*vcm[2]; double escale = (ke + heat - 0.5*vcmsq*masstotal)/(ke - 0.5*vcmsq*masstotal); if (escale < 0.0) error->all("Fix heat kinetic energy went negative"); @@ -86,13 +109,27 @@ void FixHeat::end_of_step() vsub[0] = (scale-1.0) * vcm[0]; vsub[1] = (scale-1.0) * vcm[1]; vsub[2] = (scale-1.0) * vcm[2]; - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - v[i][0] = scale*v[i][0] - vsub[0]; - v[i][1] = scale*v[i][1] - vsub[1]; - v[i][2] = scale*v[i][2] - vsub[2]; - } + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + if (iregion < 0) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + v[i][0] = scale*v[i][0] - vsub[0]; + v[i][1] = scale*v[i][1] - vsub[1]; + v[i][2] = scale*v[i][2] - vsub[2]; + } + } else { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { + v[i][0] = scale*v[i][0] - vsub[0]; + v[i][1] = scale*v[i][1] - vsub[1]; + v[i][2] = scale*v[i][2] - vsub[2]; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_heat.h b/src/fix_heat.h index 971cbf67eb..cb59bf0f00 100644 --- a/src/fix_heat.h +++ b/src/fix_heat.h @@ -33,6 +33,7 @@ class FixHeat : public Fix { double compute_scalar(); private: + int iregion; double heat_input; double masstotal; double scale; diff --git a/src/variable.cpp b/src/variable.cpp index 90a44f1a48..9ec17ab2ba 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -203,6 +203,7 @@ void Variable::set(int narg, char **arg) num[nvar] != num[jvar]) error->all("All universe/uloop variables must have same # of values"); + /* if (me == 0) { if (universe->uscreen) fprintf(universe->uscreen, @@ -213,6 +214,7 @@ void Variable::set(int narg, char **arg) "Initial ${%s} setting: value %d on partition %d\n", arg[0],index[nvar]+1,universe->iworld); } + */ // ATOM // remove pre-existing var if also style ATOM (allows it to be reset)