git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@4014 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp
2010-04-21 19:47:39 +00:00
parent 72b018027b
commit ff8c37dfe1
4 changed files with 55 additions and 15 deletions

View File

@ -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");

View File

@ -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];
}
}
}
/* ---------------------------------------------------------------------- */

View File

@ -33,6 +33,7 @@ class FixHeat : public Fix {
double compute_scalar();
private:
int iregion;
double heat_input;
double masstotal;
double scale;

View File

@ -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)