diff --git a/src/compute_reduce.cpp b/src/compute_reduce.cpp
index bb750a2ad..fdc302960 100644
--- a/src/compute_reduce.cpp
+++ b/src/compute_reduce.cpp
@@ -1,653 +1,656 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "compute_reduce.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "modify.h"
#include "fix.h"
#include "force.h"
#include "comm.h"
#include "group.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{SUM,MINN,MAXX,AVE};
enum{X,V,F,COMPUTE,FIX,VARIABLE};
enum{PERATOM,LOCAL};
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
#define INVOKED_PERATOM 8
#define INVOKED_LOCAL 16
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
ComputeReduce::ComputeReduce(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
int iarg;
if (strcmp(style,"reduce") == 0) {
if (narg < 5) error->all(FLERR,"Illegal compute reduce command");
idregion = NULL;
iarg = 3;
} else if (strcmp(style,"reduce/region") == 0) {
if (narg < 6) error->all(FLERR,"Illegal compute reduce/region command");
iregion = domain->find_region(arg[3]);
if (iregion == -1)
error->all(FLERR,"Region ID for compute reduce/region does not exist");
int n = strlen(arg[3]) + 1;
idregion = new char[n];
strcpy(idregion,arg[3]);
iarg = 4;
}
if (strcmp(arg[iarg],"sum") == 0) mode = SUM;
else if (strcmp(arg[iarg],"min") == 0) mode = MINN;
else if (strcmp(arg[iarg],"max") == 0) mode = MAXX;
else if (strcmp(arg[iarg],"ave") == 0) mode = AVE;
else error->all(FLERR,"Illegal compute reduce command");
iarg++;
MPI_Comm_rank(world,&me);
// parse remaining values until one isn't recognized
which = new int[narg-4];
argindex = new int[narg-4];
flavor = new int[narg-4];
ids = new char*[narg-4];
value2index = new int[narg-4];
nvalues = 0;
while (iarg < narg) {
ids[nvalues] = NULL;
if (strcmp(arg[iarg],"x") == 0) {
which[nvalues] = X;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"y") == 0) {
which[nvalues] = X;
argindex[nvalues++] = 1;
} else if (strcmp(arg[iarg],"z") == 0) {
which[nvalues] = X;
argindex[nvalues++] = 2;
} else if (strcmp(arg[iarg],"vx") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"vy") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 1;
} else if (strcmp(arg[iarg],"vz") == 0) {
which[nvalues] = V;
argindex[nvalues++] = 2;
} else if (strcmp(arg[iarg],"fx") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 0;
} else if (strcmp(arg[iarg],"fy") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 1;
} else if (strcmp(arg[iarg],"fz") == 0) {
which[nvalues] = F;
argindex[nvalues++] = 2;
} else if (strncmp(arg[iarg],"c_",2) == 0 ||
strncmp(arg[iarg],"f_",2) == 0 ||
strncmp(arg[iarg],"v_",2) == 0) {
if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
else if (arg[iarg][0] == 'v') which[nvalues] = VARIABLE;
int n = strlen(arg[iarg]);
char *suffix = new char[n];
strcpy(suffix,&arg[iarg][2]);
char *ptr = strchr(suffix,'[');
if (ptr) {
if (suffix[strlen(suffix)-1] != ']')
error->all(FLERR,"Illegal compute reduce command");
argindex[nvalues] = atoi(ptr+1);
*ptr = '\0';
} else argindex[nvalues] = 0;
n = strlen(suffix) + 1;
ids[nvalues] = new char[n];
strcpy(ids[nvalues],suffix);
nvalues++;
delete [] suffix;
} else break;
iarg++;
}
// optional args
replace = new int[nvalues];
for (int i = 0; i < nvalues; i++) replace[i] = -1;
while (iarg < narg) {
if (strcmp(arg[iarg],"replace") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal compute reduce command");
if (mode != MINN && mode != MAXX)
error->all(FLERR,"Compute reduce replace requires min or max mode");
int col1 = atoi(arg[iarg+1]) - 1;
int col2 = atoi(arg[iarg+2]) - 1;
if (col1 < 0 || col1 >= nvalues || col2 < 0 || col2 >= nvalues)
error->all(FLERR,"Illegal compute reduce command");
if (col1 == col2) error->all(FLERR,"Illegal compute reduce command");
if (replace[col1] >= 0 || replace[col2] >= 0)
error->all(FLERR,"Invalid replace values in compute reduce");
replace[col1] = col2;
iarg += 3;
} else error->all(FLERR,"Illegal compute reduce command");
}
// delete replace if not set
int flag = 0;
for (int i = 0; i < nvalues; i++)
if (replace[i] >= 0) flag = 1;
if (!flag) {
delete [] replace;
replace = NULL;
}
// setup and error check
for (int i = 0; i < nvalues; i++) {
if (which[i] == X || which[i] == V || which[i] == F)
flavor[i] = PERATOM;
else if (which[i] == COMPUTE) {
int icompute = modify->find_compute(ids[i]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute reduce does not exist");
if (modify->compute[icompute]->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 &&
modify->compute[icompute]->size_peratom_cols != 0)
error->all(FLERR,"Compute reduce compute does not "
"calculate a per-atom vector");
if (argindex[i] && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute reduce compute does not "
"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 &&
modify->compute[icompute]->size_local_cols != 0)
error->all(FLERR,"Compute reduce compute does not "
"calculate a local vector");
if (argindex[i] && modify->compute[icompute]->size_local_cols == 0)
error->all(FLERR,"Compute reduce compute does not "
"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]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute reduce does not exist");
if (modify->fix[ifix]->peratom_flag) {
flavor[i] = PERATOM;
if (argindex[i] == 0 &&
modify->fix[ifix]->size_peratom_cols != 0)
error->all(FLERR,"Compute reduce fix does not "
"calculate a per-atom vector");
if (argindex[i] && modify->fix[ifix]->size_peratom_cols == 0)
error->all(FLERR,"Compute reduce fix does not "
"calculate a per-atom array");
if (argindex[i] &&
argindex[i] > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,"Compute reduce fix array is accessed out-of-range");
} else if (modify->fix[ifix]->local_flag) {
flavor[i] = LOCAL;
if (argindex[i] == 0 &&
modify->fix[ifix]->size_local_cols != 0)
error->all(FLERR,"Compute reduce fix does not "
"calculate a local vector");
if (argindex[i] && modify->fix[ifix]->size_local_cols == 0)
error->all(FLERR,"Compute reduce fix does not "
"calculate a local array");
if (argindex[i] &&
argindex[i] > modify->fix[ifix]->size_local_cols)
error->all(FLERR,"Compute reduce fix array is accessed out-of-range");
} else error->all(FLERR,"Compute reduce fix calculates global values");
} else if (which[i] == VARIABLE) {
int ivariable = input->variable->find(ids[i]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute reduce does not exist");
if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Compute reduce variable is not atom-style variable");
flavor[i] = PERATOM;
}
}
// this compute produces either a scalar or vector
if (nvalues == 1) {
scalar_flag = 1;
if (mode == SUM) extscalar = 1;
else extscalar = 0;
vector = onevec = NULL;
indices = owner = NULL;
} else {
vector_flag = 1;
size_vector = nvalues;
if (mode == SUM) extvector = 1;
else extvector = 0;
vector = new double[size_vector];
onevec = new double[size_vector];
indices = new int[size_vector];
owner = new int[size_vector];
}
maxatom = 0;
varatom = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeReduce::~ComputeReduce()
{
delete [] which;
delete [] argindex;
delete [] flavor;
for (int m = 0; m < nvalues; m++) delete [] ids[m];
delete [] ids;
delete [] value2index;
delete [] replace;
delete [] idregion;
delete [] vector;
delete [] onevec;
delete [] indices;
delete [] owner;
memory->destroy(varatom);
}
/* ---------------------------------------------------------------------- */
void ComputeReduce::init()
{
// set indices and check validity of all computes,fixes,variables
for (int m = 0; m < nvalues; m++) {
if (which[m] == COMPUTE) {
int icompute = modify->find_compute(ids[m]);
if (icompute < 0)
error->all(FLERR,"Compute ID for compute reduce does not exist");
value2index[m] = icompute;
} else if (which[m] == FIX) {
int ifix = modify->find_fix(ids[m]);
if (ifix < 0)
error->all(FLERR,"Fix ID for compute reduce does not exist");
value2index[m] = ifix;
} else if (which[m] == VARIABLE) {
int ivariable = input->variable->find(ids[m]);
if (ivariable < 0)
error->all(FLERR,"Variable name for compute reduce does not exist");
value2index[m] = ivariable;
} else value2index[m] = -1;
}
// set index and check validity of region
if (idregion) {
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR,"Region ID for compute reduce/region does not exist");
}
}
/* ---------------------------------------------------------------------- */
double ComputeReduce::compute_scalar()
{
invoked_scalar = update->ntimestep;
double one = compute_one(0,-1);
if (mode == SUM) {
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
} else if (mode == MINN) {
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MIN,world);
} else if (mode == MAXX) {
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_MAX,world);
} else if (mode == AVE) {
MPI_Allreduce(&one,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
bigint n = count(0);
if (n) scalar /= n;
}
return scalar;
}
/* ---------------------------------------------------------------------- */
void ComputeReduce::compute_vector()
{
invoked_vector = update->ntimestep;
for (int m = 0; m < nvalues; m++)
if (!replace || replace[m] < 0) {
onevec[m] = compute_one(m,-1);
indices[m] = index;
}
if (mode == SUM) {
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
} else if (mode == MINN) {
if (!replace) {
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MIN,world);
} else {
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MINLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
if (me == owner[replace[m]])
vector[m] = compute_one(m,indices[replace[m]]);
MPI_Bcast(&vector[m],1,MPI_DOUBLE,owner[replace[m]],world);
}
}
} else if (mode == MAXX) {
if (!replace) {
for (int m = 0; m < nvalues; m++)
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_MAX,world);
} else {
for (int m = 0; m < nvalues; m++)
if (replace[m] < 0) {
pairme.value = onevec[m];
pairme.proc = me;
MPI_Allreduce(&pairme,&pairall,1,MPI_DOUBLE_INT,MPI_MAXLOC,world);
vector[m] = pairall.value;
owner[m] = pairall.proc;
}
for (int m = 0; m < nvalues; m++)
if (replace[m] >= 0) {
if (me == owner[replace[m]])
vector[m] = compute_one(m,indices[replace[m]]);
MPI_Bcast(&vector[m],1,MPI_DOUBLE,owner[replace[m]],world);
}
}
} else if (mode == AVE) {
for (int m = 0; m < nvalues; m++) {
MPI_Allreduce(&onevec[m],&vector[m],1,MPI_DOUBLE,MPI_SUM,world);
bigint n = count(m);
if (n) vector[m] /= n;
}
}
}
/* ----------------------------------------------------------------------
calculate reduced value for one input M and return it
if flag = -1:
sum/min/max/ave all values in vector
for per-atom quantities, limit to atoms in group
if mode = MIN or MAX, also set index to which vector value wins
if flag >= 0: simply return vector[flag]
------------------------------------------------------------------------- */
double ComputeReduce::compute_one(int m, int flag)
{
int i;
// invoke the appropriate attribute,compute,fix,variable
// for flag = -1, compute scalar quantity by scanning over atom properties
// only include atoms in group for atom properties and per-atom quantities
index = -1;
int vidx = value2index[m];
int aidx = argindex[m];
int *mask = atom->mask;
int nlocal = atom->nlocal;
double one;
if (mode == SUM) one = 0.0;
else if (mode == MINN) one = BIG;
else if (mode == MAXX) one = -BIG;
else if (mode == AVE) one = 0.0;
if (which[m] == X) {
double **x = atom->x;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,x[i][aidx],i);
} else one = x[flag][aidx];
} else if (which[m] == V) {
double **v = atom->v;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,v[i][aidx],i);
} else one = v[flag][aidx];
} else if (which[m] == F) {
double **f = atom->f;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,f[i][aidx],i);
} else one = f[flag][aidx];
// invoke compute if not previously invoked
} else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[vidx];
if (flavor[m] == PERATOM) {
if (!(compute->invoked_flag & INVOKED_PERATOM)) {
compute->compute_peratom();
compute->invoked_flag |= INVOKED_PERATOM;
}
if (aidx == 0) {
double *comp_vec = compute->vector_atom;
int n = nlocal;
if (flag < 0) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one,comp_vec[i],i);
} else one = comp_vec[flag];
} else {
double **carray_atom = compute->array_atom;
int n = nlocal;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one,carray_atom[i][aidxm1],i);
} else one = carray_atom[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (!(compute->invoked_flag & INVOKED_LOCAL)) {
compute->compute_local();
compute->invoked_flag |= INVOKED_LOCAL;
}
if (aidx == 0) {
double *comp_vec = compute->vector_local;
int n = compute->size_local_rows;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,comp_vec[i],i);
else one = comp_vec[flag];
} else {
double **carray_local = compute->array_local;
int n = compute->size_local_rows;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,carray_local[i][aidxm1],i);
else one = carray_local[flag][aidxm1];
}
}
// access fix fields, check if fix frequency is a match
} else if (which[m] == FIX) {
if (update->ntimestep % modify->fix[vidx]->peratom_freq)
error->all(FLERR,"Fix used in compute reduce not computed at compatible time");
Fix *fix = modify->fix[vidx];
if (flavor[m] == PERATOM) {
if (aidx == 0) {
double *fix_vector = fix->vector_atom;
int n = nlocal;
if (flag < 0) {
for (i = 0; i < n; i++)
if (mask[i] & groupbit) combine(one,fix_vector[i],i);
} else one = fix_vector[flag];
} else {
double **fix_array = fix->array_atom;
int aidxm1 = aidx - 1;
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,fix_array[i][aidxm1],i);
} else one = fix_array[flag][aidxm1];
}
} else if (flavor[m] == LOCAL) {
if (aidx == 0) {
double *fix_vector = fix->vector_local;
int n = fix->size_local_rows;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,fix_vector[i],i);
else one = fix_vector[flag];
} else {
double **fix_array = fix->array_local;
int n = fix->size_local_rows;
int aidxm1 = aidx - 1;
if (flag < 0)
for (i = 0; i < n; i++)
combine(one,fix_array[i][aidxm1],i);
else one = fix_array[flag][aidxm1];
}
}
// evaluate atom-style variable
} else if (which[m] == VARIABLE) {
if (nlocal > maxatom) {
maxatom = atom->nmax;
memory->destroy(varatom);
memory->create(varatom,maxatom,"reduce:varatom");
}
input->variable->compute_atom(vidx,igroup,varatom,1,0);
if (flag < 0) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) combine(one,varatom[i],i);
} else one = varatom[flag];
}
return one;
}
/* ---------------------------------------------------------------------- */
bigint ComputeReduce::count(int m)
{
int vidx = value2index[m];
if (which[m] == X || which[m] == V || which[m] == F)
return group->count(igroup);
else if (which[m] == COMPUTE) {
Compute *compute = modify->compute[vidx];
if (flavor[m] == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
bigint ncount = compute->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount,&ncountall,1,MPI_LMP_BIGINT,MPI_SUM,world);
return ncountall;
}
} else if (which[m] == FIX) {
Fix *fix = modify->fix[vidx];
if (flavor[m] == PERATOM) {
return group->count(igroup);
} else if (flavor[m] == LOCAL) {
bigint ncount = fix->size_local_rows;
bigint ncountall;
MPI_Allreduce(&ncount,&ncountall,1,MPI_LMP_BIGINT,MPI_SUM,world);
return ncountall;
}
} else if (which[m] == VARIABLE)
return group->count(igroup);
bigint dummy = 0;
return dummy;
}
/* ----------------------------------------------------------------------
combine two values according to reduction mode
for MIN/MAX, also update index with winner
------------------------------------------------------------------------- */
void ComputeReduce::combine(double &one, double two, int i)
{
if (mode == SUM || mode == AVE) one += two;
else if (mode == MINN) {
if (two < one) {
one = two;
index = i;
}
} else if (mode == MAXX) {
if (two > one) {
one = two;
index = i;
}
}
}
/* ----------------------------------------------------------------------
memory usage of varatom
------------------------------------------------------------------------- */
double ComputeReduce::memory_usage()
{
double bytes = maxatom * sizeof(double);
return bytes;
}
diff --git a/src/respa.cpp b/src/respa.cpp
index 19a6f7d80..339c965d1 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -1,705 +1,697 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Mark Stevens (SNL), Paul Crozier (SNL)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "respa.h"
#include "neighbor.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "output.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "fix_respa.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
Respa::Respa(LAMMPS *lmp, int narg, char **arg) : Integrate(lmp, narg, arg)
{
if (narg < 1) error->all(FLERR,"Illegal run_style respa command");
nlevels = force->inumeric(FLERR,arg[0]);
if (nlevels < 1) error->all(FLERR,"Respa levels must be >= 1");
if (narg < nlevels) error->all(FLERR,"Illegal run_style respa command");
loop = new int[nlevels];
for (int iarg = 1; iarg < nlevels; iarg++) {
loop[iarg-1] = force->inumeric(FLERR,arg[iarg]);
if (loop[iarg-1] <= 0) error->all(FLERR,"Illegal run_style respa command");
}
loop[nlevels-1] = 1;
// set level at which each force is computed
// argument settings override defaults
level_bond = level_angle = level_dihedral = level_improper = -1;
level_pair = level_kspace = -1;
level_inner = level_middle = level_outer = -1;
int iarg = nlevels;
while (iarg < narg) {
if (strcmp(arg[iarg],"bond") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_bond = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"angle") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_angle = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_dihedral = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"improper") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_improper = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"pair") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_pair = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"inner") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
level_inner = force->inumeric(FLERR,arg[iarg+1]) - 1;
cutoff[0] = force->numeric(FLERR,arg[iarg+2]);
cutoff[1] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"middle") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal run_style respa command");
level_middle = force->inumeric(FLERR,arg[iarg+1]) - 1;
cutoff[2] = force->numeric(FLERR,arg[iarg+2]);
cutoff[3] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"outer") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_outer = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else if (strcmp(arg[iarg],"kspace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal run_style respa command");
level_kspace = force->inumeric(FLERR,arg[iarg+1]) - 1;
iarg += 2;
} else error->all(FLERR,"Illegal run_style respa command");
}
// cannot specify both pair and inner/middle/outer
if (level_pair >= 0 &&
(level_inner >= 0 || level_middle >= 0 || level_outer >= 0))
error->all(FLERR,"Cannot set both respa pair and inner/middle/outer");
// if either inner and outer is specified, then both must be
if ((level_inner >= 0 && level_outer == -1) ||
(level_outer >= 0 && level_inner == -1))
error->all(FLERR,"Must set both respa inner and outer");
// middle cannot be set without inner/outer
if (level_middle >= 0 && level_inner == -1)
error->all(FLERR,"Cannot set respa middle without inner/outer");
// set defaults if user did not specify level
// bond to innermost level
// angle same as bond, dihedral same as angle, improper same as dihedral
// pair to outermost level if no inner/middle/outer
// inner/middle/outer have no defaults
// kspace same as pair or outer
if (level_bond == -1) level_bond = 0;
if (level_angle == -1) level_angle = level_bond;
if (level_dihedral == -1) level_dihedral = level_angle;
if (level_improper == -1) level_improper = level_dihedral;
if (level_pair == -1 && level_inner == -1) level_pair = nlevels-1;
if (level_kspace == -1 && level_pair >= 0) level_kspace = level_pair;
if (level_kspace == -1 && level_pair == -1) level_kspace = level_outer;
// print respa levels
if (comm->me == 0) {
if (screen) {
fprintf(screen,"Respa levels:\n");
for (int i = 0; i < nlevels; i++) {
fprintf(screen," %d =",i+1);
if (level_bond == i) fprintf(screen," bond");
if (level_angle == i) fprintf(screen," angle");
if (level_dihedral == i) fprintf(screen," dihedral");
if (level_improper == i) fprintf(screen," improper");
if (level_pair == i) fprintf(screen," pair");
if (level_inner == i) fprintf(screen," pair-inner");
if (level_middle == i) fprintf(screen," pair-middle");
if (level_outer == i) fprintf(screen," pair-outer");
if (level_kspace == i) fprintf(screen," kspace");
fprintf(screen,"\n");
}
}
if (logfile) {
fprintf(logfile,"Respa levels:\n");
for (int i = 0; i < nlevels; i++) {
fprintf(logfile," %d =",i+1);
if (level_bond == i) fprintf(logfile," bond");
if (level_angle == i) fprintf(logfile," angle");
if (level_dihedral == i) fprintf(logfile," dihedral");
if (level_improper == i) fprintf(logfile," improper");
if (level_pair == i) fprintf(logfile," pair");
if (level_inner == i) fprintf(logfile," pair-inner");
if (level_middle == i) fprintf(logfile," pair-middle");
if (level_outer == i) fprintf(logfile," pair-outer");
if (level_kspace == i) fprintf(logfile," kspace");
fprintf(logfile,"\n");
}
}
}
// check that levels are in correct order
if (level_angle < level_bond || level_dihedral < level_angle ||
level_improper < level_dihedral)
error->all(FLERR,"Invalid order of forces within respa levels");
if (level_pair >= 0) {
if (level_pair < level_improper || level_kspace < level_pair)
error->all(FLERR,"Invalid order of forces within respa levels");
}
if (level_pair == -1 && level_middle == -1) {
if (level_inner < level_improper || level_outer < level_inner ||
level_kspace < level_outer)
error->all(FLERR,"Invalid order of forces within respa levels");
}
if (level_pair == -1 && level_middle >= 0) {
if (level_inner < level_improper || level_middle < level_inner ||
level_outer < level_inner || level_kspace < level_outer)
error->all(FLERR,"Invalid order of forces within respa levels");
}
// warn if any levels are devoid of forces
int flag = 0;
for (int i = 0; i < nlevels; i++)
if (level_bond != i && level_angle != i && level_dihedral != i &&
level_improper != i && level_pair != i && level_inner != i &&
level_middle != i && level_outer != i && level_kspace != i) flag = 1;
if (flag && comm->me == 0)
error->warning(FLERR,"One or more respa levels compute no forces");
// check cutoff consistency if inner/middle/outer are enabled
if (level_inner >= 0 && cutoff[1] < cutoff[0])
error->all(FLERR,"Respa inner cutoffs are invalid");
if (level_middle >= 0 && (cutoff[3] < cutoff[2] || cutoff[2] < cutoff[1]))
error->all(FLERR,"Respa middle cutoffs are invalid");
// set outer pair of cutoffs to inner pair if middle is not enabled
if (level_inner >= 0 && level_middle < 0) {
cutoff[2] = cutoff[0];
cutoff[3] = cutoff[1];
}
// allocate other needed arrays
newton = new int[nlevels];
step = new double[nlevels];
}
/* ---------------------------------------------------------------------- */
Respa::~Respa()
{
delete [] loop;
delete [] newton;
delete [] step;
}
/* ----------------------------------------------------------------------
initialization before run
------------------------------------------------------------------------- */
void Respa::init()
{
Integrate::init();
// warn if no fixes
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
char **fixarg = new char*[4];
fixarg[0] = (char *) "RESPA";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "RESPA";
fixarg[3] = new char[8];
sprintf(fixarg[3],"%d",nlevels);
modify->add_fix(4,fixarg);
delete [] fixarg[3];
delete [] fixarg;
fix_respa = (FixRespa *) modify->fix[modify->nfix-1];
// insure respa inner/middle/outer is using Pair class that supports it
if (level_inner >= 0)
if (force->pair && force->pair->respa_enable == 0)
error->all(FLERR,"Pair style does not support rRESPA inner/middle/outer");
// virial_style = 1 (explicit) since never computed implicitly like Verlet
virial_style = 1;
// setup lists of computes for global and per-atom PE and pressure
ev_setup();
// detect if fix omp is present and will clear force arrays
int ifix = modify->find_fix("package_omp");
if (ifix >= 0) external_force_clear = 1;
// set flags for what arrays to clear in force_clear()
// need to clear additionals arrays if they exist
torqueflag = 0;
if (atom->torque_flag) torqueflag = 1;
erforceflag = 0;
if (atom->erforce_flag) erforceflag = 1;
e_flag = 0;
if (atom->e_flag) e_flag = 1;
rho_flag = 0;
if (atom->rho_flag) rho_flag = 1;
// step[] = timestep for each level
step[nlevels-1] = update->dt;
for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
step[ilevel] = step[ilevel+1]/loop[ilevel];
// set newton flag for each level
for (int ilevel = 0; ilevel < nlevels; ilevel++) {
newton[ilevel] = 0;
if (force->newton_bond) {
if (level_bond == ilevel || level_angle == ilevel ||
level_dihedral == ilevel || level_improper == ilevel)
newton[ilevel] = 1;
}
if (force->newton_pair) {
if (level_pair == ilevel || level_inner == ilevel ||
level_middle == ilevel || level_outer == ilevel)
newton[ilevel] = 1;
}
}
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
}
/* ----------------------------------------------------------------------
setup before run
------------------------------------------------------------------------- */
void Respa::setup()
{
if (comm->me == 0 && screen) fprintf(screen,"Setting up run ...\n");
update->setupflag = 1;
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
atom->setup();
modify->setup_pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
if (atom->sortfreq > 0) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
neighbor->ncalls = 0;
// compute all forces
ev_set(update->ntimestep);
for (int ilevel = 0; ilevel < nlevels; ilevel++) {
force_clear(newton[ilevel]);
modify->setup_pre_force_respa(vflag,ilevel);
if (level_pair == ilevel && pair_compute_flag)
force->pair->compute(eflag,vflag);
if (level_inner == ilevel && pair_compute_flag)
force->pair->compute_inner();
if (level_middle == ilevel && pair_compute_flag)
force->pair->compute_middle();
if (level_outer == ilevel && pair_compute_flag)
force->pair->compute_outer(eflag,vflag);
if (level_bond == ilevel && force->bond)
force->bond->compute(eflag,vflag);
if (level_angle == ilevel && force->angle)
force->angle->compute(eflag,vflag);
if (level_dihedral == ilevel && force->dihedral)
force->dihedral->compute(eflag,vflag);
if (level_improper == ilevel && force->improper)
force->improper->compute(eflag,vflag);
if (level_kspace == ilevel && force->kspace) {
force->kspace->setup();
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
modify->setup(vflag);
sum_flevel_f();
output->setup();
update->setupflag = 0;
}
/* ----------------------------------------------------------------------
setup without output
flag = 0 = just force calculation
flag = 1 = reneighbor and force calculation
------------------------------------------------------------------------- */
void Respa::setup_minimal(int flag)
{
update->setupflag = 1;
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
if (flag) {
modify->setup_pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
domain->image_check();
domain->box_too_small_check();
modify->setup_pre_neighbor();
neighbor->build();
neighbor->ncalls = 0;
}
// compute all forces
ev_set(update->ntimestep);
for (int ilevel = 0; ilevel < nlevels; ilevel++) {
force_clear(newton[ilevel]);
modify->setup_pre_force_respa(vflag,ilevel);
if (level_pair == ilevel && pair_compute_flag)
force->pair->compute(eflag,vflag);
if (level_inner == ilevel && pair_compute_flag)
force->pair->compute_inner();
if (level_middle == ilevel && pair_compute_flag)
force->pair->compute_middle();
if (level_outer == ilevel && pair_compute_flag)
force->pair->compute_outer(eflag,vflag);
if (level_bond == ilevel && force->bond)
force->bond->compute(eflag,vflag);
if (level_angle == ilevel && force->angle)
force->angle->compute(eflag,vflag);
if (level_dihedral == ilevel && force->dihedral)
force->dihedral->compute(eflag,vflag);
if (level_improper == ilevel && force->improper)
force->improper->compute(eflag,vflag);
if (level_kspace == ilevel && force->kspace) {
force->kspace->setup();
if (kspace_compute_flag) force->kspace->compute(eflag,vflag);
}
if (newton[ilevel]) comm->reverse_comm();
copy_f_flevel(ilevel);
}
modify->setup(vflag);
sum_flevel_f();
update->setupflag = 0;
}
/* ----------------------------------------------------------------------
run for N steps
------------------------------------------------------------------------- */
void Respa::run(int n)
{
bigint ntimestep;
for (int i = 0; i < n; i++) {
ntimestep = ++update->ntimestep;
ev_set(ntimestep);
recurse(nlevels-1);
if (modify->n_end_of_step) modify->end_of_step();
if (ntimestep == output->next) {
timer->stamp();
sum_flevel_f();
output->write(update->ntimestep);
timer->stamp(TIME_OUTPUT);
}
}
}
/* ----------------------------------------------------------------------
delete rRESPA fix at end of run, so its atom arrays won't persist
------------------------------------------------------------------------- */
void Respa::cleanup()
{
modify->post_run();
modify->delete_fix("RESPA");
domain->box_too_small_check();
update->update_time();
}
/* ---------------------------------------------------------------------- */
void Respa::reset_dt()
{
step[nlevels-1] = update->dt;
for (int ilevel = nlevels-2; ilevel >= 0; ilevel--)
step[ilevel] = step[ilevel+1]/loop[ilevel];
}
/* ---------------------------------------------------------------------- */
void Respa::recurse(int ilevel)
{
copy_flevel_f(ilevel);
for (int iloop = 0; iloop < loop[ilevel]; iloop++) {
modify->initial_integrate_respa(vflag,ilevel,iloop);
if (modify->n_post_integrate_respa)
modify->post_integrate_respa(ilevel,iloop);
if (ilevel) recurse(ilevel-1);
// at outermost level, check on rebuilding neighbor list
// at innermost level, communicate
// at middle levels, do nothing
if (ilevel == nlevels-1) {
int nflag = neighbor->decide();
if (nflag) {
if (modify->n_pre_exchange) modify->pre_exchange();
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
if (atom->sortfreq > 0 &&
update->ntimestep >= atom->nextsort) atom->sort();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
if (modify->n_pre_neighbor) modify->pre_neighbor();
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
} else if (ilevel == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(TIME_COMM);
}
} else if (ilevel == 0) {
timer->stamp();
comm->forward_comm();
timer->stamp(TIME_COMM);
}
// force computations
// important that ordering is same as Verlet
// so that any order dependencies are the same
// when potentials are invoked at same level
force_clear(newton[ilevel]);
if (modify->n_pre_force_respa)
modify->pre_force_respa(vflag,ilevel,iloop);
timer->stamp();
if (level_pair == ilevel && pair_compute_flag) {
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (level_inner == ilevel && pair_compute_flag) {
force->pair->compute_inner();
timer->stamp(TIME_PAIR);
}
if (level_middle == ilevel && pair_compute_flag) {
force->pair->compute_middle();
timer->stamp(TIME_PAIR);
}
if (level_outer == ilevel && pair_compute_flag) {
force->pair->compute_outer(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (level_bond == ilevel && force->bond) {
force->bond->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_angle == ilevel && force->angle) {
force->angle->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_dihedral == ilevel && force->dihedral) {
force->dihedral->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_improper == ilevel && force->improper) {
force->improper->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (level_kspace == ilevel && kspace_compute_flag) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
if (newton[ilevel]) {
comm->reverse_comm();
timer->stamp(TIME_COMM);
}
if (modify->n_post_force_respa)
modify->post_force_respa(vflag,ilevel,iloop);
modify->final_integrate_respa(ilevel,iloop);
}
copy_f_flevel(ilevel);
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
------------------------------------------------------------------------- */
void Respa::force_clear(int newtonflag)
{
if (external_force_clear) return;
// clear global force array
// nall includes ghosts only if newton flag is set
int nall;
if (newtonflag) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
size_t nbytes = sizeof(double) * nall;
if (nbytes > 0 ) {
memset(&(atom->f[0][0]),0,3*nbytes);
if (torqueflag) memset(&(atom->torque[0][0]),0,3*nbytes);
if (erforceflag) memset(&(atom->erforce[0]), 0, nbytes);
if (e_flag) memset(&(atom->de[0]), 0, nbytes);
if (rho_flag) memset(&(atom->drho[0]), 0, nbytes);
}
}
/* ----------------------------------------------------------------------
copy force components from atom->f to FixRespa->f_level
------------------------------------------------------------------------- */
void Respa::copy_f_flevel(int ilevel)
{
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int i = 0; i < n; i++) {
f_level[i][ilevel][0] = f[i][0];
f_level[i][ilevel][1] = f[i][1];
f_level[i][ilevel][2] = f[i][2];
}
}
/* ----------------------------------------------------------------------
copy force components from FixRespa->f_level to atom->f
------------------------------------------------------------------------- */
void Respa::copy_flevel_f(int ilevel)
{
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int i = 0; i < n; i++) {
f[i][0] = f_level[i][ilevel][0];
f[i][1] = f_level[i][ilevel][1];
f[i][2] = f_level[i][ilevel][2];
}
}
/* ----------------------------------------------------------------------
sum all force components from FixRespa->f_level to create full atom->f
------------------------------------------------------------------------- */
void Respa::sum_flevel_f()
{
copy_flevel_f(0);
double ***f_level = fix_respa->f_level;
double **f = atom->f;
int n = atom->nlocal;
for (int ilevel = 1; ilevel < nlevels; ilevel++) {
for (int i = 0; i < n; i++) {
f[i][0] += f_level[i][ilevel][0];
f[i][1] += f_level[i][ilevel][1];
f[i][2] += f_level[i][ilevel][2];
}
}
}
diff --git a/src/velocity.cpp b/src/velocity.cpp
index a5c6a8673..701a4f31e 100644
--- a/src/velocity.cpp
+++ b/src/velocity.cpp
@@ -1,790 +1,827 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "velocity.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "lattice.h"
#include "input.h"
#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"
#include "group.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{CREATE,SET,SCALE,RAMP,ZERO};
enum{ALL,LOCAL,GEOM};
enum{NONE,CONSTANT,EQUAL,ATOM};
#define WARMUP 100
#define SMALL 0.001
/* ---------------------------------------------------------------------- */
Velocity::Velocity(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
void Velocity::command(int narg, char **arg)
{
if (narg < 2) error->all(FLERR,"Illegal velocity command");
if (domain->box_exist == 0)
error->all(FLERR,"Velocity command before simulation box is defined");
if (atom->natoms == 0)
error->all(FLERR,"Velocity command with no atoms existing");
// atom masses must all be set
atom->check_mass();
// identify group
igroup = group->find(arg[0]);
if (igroup == -1) error->all(FLERR,"Could not find velocity group ID");
groupbit = group->bitmask[igroup];
// identify style
if (strcmp(arg[1],"create") == 0) style = CREATE;
else if (strcmp(arg[1],"set") == 0) style = SET;
else if (strcmp(arg[1],"scale") == 0) style = SCALE;
else if (strcmp(arg[1],"ramp") == 0) style = RAMP;
else if (strcmp(arg[1],"zero") == 0) style = ZERO;
else error->all(FLERR,"Illegal velocity command");
// set defaults
temperature = NULL;
dist_flag = 0;
sum_flag = 0;
momentum_flag = 1;
rotation_flag = 0;
loop_flag = ALL;
scale_flag = 1;
+ rfix = -1;
// read options from end of input line
// change defaults as options specify
if (style == CREATE) options(narg-4,&arg[4]);
else if (style == SET) options(narg-5,&arg[5]);
else if (style == SCALE) options(narg-3,&arg[3]);
else if (style == RAMP) options(narg-8,&arg[8]);
else if (style == ZERO) options(narg-3,&arg[3]);
// initialize velocities based on style
// create() invoked differently, so can be called externally
if (style == CREATE) {
double t_desired = force->numeric(FLERR,arg[2]);
int seed = force->inumeric(FLERR,arg[3]);
create(t_desired,seed);
}
else if (style == SET) set(narg-2,&arg[2]);
else if (style == SCALE) scale(narg-2,&arg[2]);
else if (style == RAMP) ramp(narg-2,&arg[2]);
else if (style == ZERO) zero(narg-2,&arg[2]);
}
/* ----------------------------------------------------------------------
initialization of defaults before calling velocity methods externaly
------------------------------------------------------------------------- */
void Velocity::init_external(const char *extgroup)
{
igroup = group->find(extgroup);
if (igroup == -1) error->all(FLERR,"Could not find velocity group ID");
groupbit = group->bitmask[igroup];
temperature = NULL;
dist_flag = 0;
sum_flag = 0;
momentum_flag = 1;
rotation_flag = 0;
loop_flag = ALL;
scale_flag = 1;
}
/* ---------------------------------------------------------------------- */
void Velocity::create(double t_desired, int seed)
{
int i;
if (seed <= 0) error->all(FLERR,"Illegal velocity create command");
// if temperature = NULL, create a new ComputeTemp with the velocity group
int tflag = 0;
if (temperature == NULL) {
char **arg = new char*[3];
arg[0] = (char *) "velocity_temp";
arg[1] = group->names[igroup];
arg[2] = (char *) "temp";
temperature = new ComputeTemp(lmp,3,arg);
tflag = 1;
delete [] arg;
}
// initialize temperature computation
// warn if groups don't match
if (igroup != temperature->igroup && comm->me == 0)
error->warning(FLERR,"Mismatch between velocity and compute groups");
temperature->init();
temperature->setup();
// store a copy of current velocities
double **v = atom->v;
int nlocal = atom->nlocal;
double **vhold;
memory->create(vhold,nlocal,3,"velocity:vnew");
for (i = 0; i < nlocal; i++) {
vhold[i][0] = v[i][0];
vhold[i][1] = v[i][1];
vhold[i][2] = v[i][2];
}
// create new velocities, in uniform or gaussian distribution
// loop option determines looping style, ALL is default
// ALL = loop over all natoms, only set those I own via atom->map
// cannot do this if atom IDs do not span 1-Natoms (some were deleted)
// will produce same V, independent of P, if atoms were read-in
// will NOT produce same V, independent of P, if used create_atoms
// LOCAL = only loop over my atoms, adjust RNG to be proc-specific
// will never produce same V, independent of P
// GEOM = only loop over my atoms
// choose RNG for each atom based on its xyz coord (geometry)
// via random->reset()
// will always produce same V, independent of P
// adjust by factor for atom mass
// for 2d, set Vz to 0.0
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int dimension = domain->dimension;
int m;
double vx,vy,vz,factor;
RanPark *random;
if (loop_flag == ALL) {
// create an atom map if one doesn't exist already
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_style = 1;
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// error check
if (atom->natoms > MAXSMALLINT)
error->all(FLERR,"Too big a problem to use velocity create loop all");
if (atom->tag_enable == 0)
error->all(FLERR,
"Cannot use velocity create loop all unless atoms have IDs");
if (atom->tag_consecutive() == 0)
error->all(FLERR,
"Atom IDs must be consecutive for velocity create loop all");
// loop over all atoms in system
// generate RNGs for all atoms, only assign to ones I own
// use either per-type mass or per-atom rmass
random = new RanPark(lmp,seed);
int natoms = static_cast<int> (atom->natoms);
for (i = 1; i <= natoms; i++) {
if (dist_flag == 0) {
vx = random->uniform();
vy = random->uniform();
vz = random->uniform();
} else {
vx = random->gaussian();
vy = random->gaussian();
vz = random->gaussian();
}
m = atom->map(i);
if (m >= 0 && m < nlocal) {
if (mask[m] & groupbit) {
if (rmass) factor = 1.0/sqrt(rmass[m]);
else factor = 1.0/sqrt(mass[type[m]]);
v[m][0] = vx * factor;
v[m][1] = vy * factor;
if (dimension == 3) v[m][2] = vz * factor;
else v[m][2] = 0.0;
}
}
}
// delete temporary atom map
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
} else if (loop_flag == LOCAL) {
random = new RanPark(lmp,seed + comm->me);
for (i = 0; i < WARMUP; i++) random->uniform();
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (dist_flag == 0) {
vx = random->uniform();
vy = random->uniform();
vz = random->uniform();
} else {
vx = random->gaussian();
vy = random->gaussian();
vz = random->gaussian();
}
if (rmass) factor = 1.0/sqrt(rmass[i]);
else factor = 1.0/sqrt(mass[type[i]]);
v[i][0] = vx * factor;
v[i][1] = vy * factor;
if (dimension == 3) v[i][2] = vz * factor;
else v[i][2] = 0.0;
}
}
} else if (loop_flag == GEOM) {
random = new RanPark(lmp,1);
double **x = atom->x;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
random->reset(seed,x[i]);
if (dist_flag == 0) {
vx = random->uniform();
vy = random->uniform();
vz = random->uniform();
} else {
vx = random->gaussian();
vy = random->gaussian();
vz = random->gaussian();
}
if (rmass) factor = 1.0/sqrt(rmass[i]);
else factor = 1.0/sqrt(mass[type[i]]);
v[i][0] = vx * factor;
v[i][1] = vy * factor;
if (dimension == 3) v[i][2] = vz * factor;
else v[i][2] = 0.0;
}
}
}
// apply momentum and rotation zeroing
if (momentum_flag) zero_momentum();
if (rotation_flag) zero_rotation();
// scale temp to desired value
double t = temperature->compute_scalar();
rescale(t,t_desired);
// if sum_flag set, add back in previous velocities
if (sum_flag) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] += vhold[i][0];
v[i][1] += vhold[i][1];
v[i][2] += vhold[i][2];
}
}
}
// free local memory
// if temperature was created, delete it
memory->destroy(vhold);
delete random;
if (tflag) delete temperature;
}
/* ---------------------------------------------------------------------- */
void Velocity::set(int narg, char **arg)
{
int xstyle,ystyle,zstyle,varflag;
double vx,vy,vz;
char *xstr,*ystr,*zstr;
int xvar,yvar,zvar;
// parse 3 args
xstyle = ystyle = zstyle = CONSTANT;
xstr = ystr = zstr = NULL;
if (strstr(arg[0],"v_") == arg[0]) {
int n = strlen(&arg[0][2]) + 1;
xstr = new char[n];
strcpy(xstr,&arg[0][2]);
} else if (strcmp(arg[0],"NULL") == 0) xstyle = NONE;
else vx = force->numeric(FLERR,arg[0]);
if (strstr(arg[1],"v_") == arg[1]) {
int n = strlen(&arg[1][2]) + 1;
ystr = new char[n];
strcpy(ystr,&arg[1][2]);
} else if (strcmp(arg[1],"NULL") == 0) ystyle = NONE;
else vy = force->numeric(FLERR,arg[1]);
if (strstr(arg[2],"v_") == arg[2]) {
int n = strlen(&arg[2][2]) + 1;
zstr = new char[n];
strcpy(zstr,&arg[2][2]);
} else if (strcmp(arg[2],"NULL") == 0) zstyle = NONE;
else vz = force->numeric(FLERR,arg[2]);
// set and apply scale factors
xscale = yscale = zscale = 1.0;
if (xstyle && !xstr) {
if (scale_flag) xscale = domain->lattice->xlattice;
vx *= xscale;
}
if (ystyle && !ystr) {
if (scale_flag) yscale = domain->lattice->ylattice;
vy *= yscale;
}
if (zstyle && !zstr) {
if (scale_flag) zscale = domain->lattice->zlattice;
vz *= zscale;
}
// check variables
if (xstr) {
xvar = input->variable->find(xstr);
if (xvar < 0)
error->all(FLERR,"Variable name for velocity set does not exist");
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
else error->all(FLERR,"Variable for velocity set is invalid style");
}
if (ystr) {
yvar = input->variable->find(ystr);
if (yvar < 0)
error->all(FLERR,"Variable name for velocity set does not exist");
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
else error->all(FLERR,"Variable for velocity set is invalid style");
}
if (zstr) {
zvar = input->variable->find(zstr);
if (zvar < 0)
error->all(FLERR,"Variable name for velocity set does not exist");
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
else error->all(FLERR,"Variable for velocity set is invalid style");
}
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
varflag = ATOM;
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
varflag = EQUAL;
else varflag = CONSTANT;
// error check for 2d models
if (domain->dimension == 2) {
if (zstyle == CONSTANT && vz != 0.0)
error->all(FLERR,"Cannot set non-zero z velocity for 2d simulation");
if (zstyle == EQUAL || zstyle == ATOM)
error->all(FLERR,"Cannot set variable z velocity for 2d simulation");
}
// allocate vfield array if necessary
double **vfield = NULL;
if (varflag == ATOM) memory->create(vfield,atom->nlocal,3,"velocity:vfield");
// set velocities via constants
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (varflag == CONSTANT) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (sum_flag == 0) {
if (xstyle) v[i][0] = vx;
if (ystyle) v[i][1] = vy;
if (zstyle) v[i][2] = vz;
} else {
if (xstyle) v[i][0] += vx;
if (ystyle) v[i][1] += vy;
if (zstyle) v[i][2] += vz;
}
}
}
// set velocities via variables
} else {
if (xstyle == EQUAL) vx = input->variable->compute_equal(xvar);
else if (xstyle == ATOM && vfield)
input->variable->compute_atom(xvar,igroup,&vfield[0][0],3,0);
if (ystyle == EQUAL) vy = input->variable->compute_equal(yvar);
else if (ystyle == ATOM && vfield)
input->variable->compute_atom(yvar,igroup,&vfield[0][1],3,0);
if (zstyle == EQUAL) vz = input->variable->compute_equal(zvar);
else if (zstyle == ATOM && vfield)
input->variable->compute_atom(zvar,igroup,&vfield[0][2],3,0);
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (sum_flag == 0) {
if (xstyle == ATOM) v[i][0] = vfield[i][0];
else if (xstyle) v[i][0] = vx;
if (ystyle == ATOM) v[i][1] = vfield[i][1];
else if (ystyle) v[i][1] = vy;
if (zstyle == ATOM) v[i][2] = vfield[i][2];
else if (zstyle) v[i][2] = vz;
} else {
if (xstyle == ATOM) v[i][0] += vfield[i][0];
else if (xstyle) v[i][0] += vx;
if (ystyle == ATOM) v[i][1] += vfield[i][1];
else if (ystyle) v[i][1] += vy;
if (zstyle == ATOM) v[i][2] += vfield[i][2];
else if (zstyle) v[i][2] += vz;
}
}
}
// clean up
delete [] xstr;
delete [] ystr;
delete [] zstr;
memory->destroy(vfield);
}
/* ----------------------------------------------------------------------
rescale velocities of a group after computing its temperature
------------------------------------------------------------------------- */
void Velocity::scale(int narg, char **arg)
{
double t_desired = force->numeric(FLERR,arg[0]);
// if temperature = NULL, create a new ComputeTemp with the velocity group
int tflag = 0;
if (temperature == NULL) {
char **arg = new char*[3];
arg[0] = (char *) "velocity_temp";
arg[1] = group->names[igroup];
arg[2] = (char *) "temp";
temperature = new ComputeTemp(lmp,3,arg);
tflag = 1;
delete [] arg;
}
// initialize temperature computation
// warn if groups don't match
if (igroup != temperature->igroup && comm->me == 0)
error->warning(FLERR,"Mismatch between velocity and compute groups");
temperature->init();
temperature->setup();
// scale temp to desired value
double t = temperature->compute_scalar();
rescale(t,t_desired);
// if temperature was created, delete it
if (tflag) delete temperature;
}
/* ----------------------------------------------------------------------
apply a ramped set of velocities
------------------------------------------------------------------------- */
void Velocity::ramp(int narg, char **arg)
{
// set scale factors
if (scale_flag) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
// parse args
int v_dim;
if (strcmp(arg[0],"vx") == 0) v_dim = 0;
else if (strcmp(arg[0],"vy") == 0) v_dim = 1;
else if (strcmp(arg[0],"vz") == 0) v_dim = 2;
else error->all(FLERR,"Illegal velocity command");
if (v_dim == 2 && domain->dimension == 2)
error->all(FLERR,"Velocity ramp in z for a 2d problem");
double v_lo,v_hi;
if (v_dim == 0) {
v_lo = xscale*force->numeric(FLERR,arg[1]);
v_hi = xscale*force->numeric(FLERR,arg[2]);
} else if (v_dim == 1) {
v_lo = yscale*force->numeric(FLERR,arg[1]);
v_hi = yscale*force->numeric(FLERR,arg[2]);
} else if (v_dim == 2) {
v_lo = zscale*force->numeric(FLERR,arg[1]);
v_hi = zscale*force->numeric(FLERR,arg[2]);
}
int coord_dim;
if (strcmp(arg[3],"x") == 0) coord_dim = 0;
else if (strcmp(arg[3],"y") == 0) coord_dim = 1;
else if (strcmp(arg[3],"z") == 0) coord_dim = 2;
else error->all(FLERR,"Illegal velocity command");
double coord_lo,coord_hi;
if (coord_dim == 0) {
coord_lo = xscale*force->numeric(FLERR,arg[4]);
coord_hi = xscale*force->numeric(FLERR,arg[5]);
} else if (coord_dim == 1) {
coord_lo = yscale*force->numeric(FLERR,arg[4]);
coord_hi = yscale*force->numeric(FLERR,arg[5]);
} else if (coord_dim == 2) {
coord_lo = zscale*force->numeric(FLERR,arg[4]);
coord_hi = zscale*force->numeric(FLERR,arg[5]);
}
// vramp = ramped velocity component for v_dim
// add or set based on sum_flag
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double fraction,vramp;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
fraction = (x[i][coord_dim] - coord_lo) / (coord_hi - coord_lo);
fraction = MAX(fraction,0.0);
fraction = MIN(fraction,1.0);
vramp = v_lo + fraction*(v_hi - v_lo);
if (sum_flag) v[i][v_dim] += vramp;
else v[i][v_dim] = vramp;
}
}
/* ----------------------------------------------------------------------
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");
}
/* ----------------------------------------------------------------------
rescale velocities of group atoms to t_new from t_old
------------------------------------------------------------------------- */
void Velocity::rescale(double t_old, double t_new)
{
if (t_old == 0.0) error->all(FLERR,"Attempting to rescale a 0.0 temperature");
double factor = sqrt(t_new/t_old);
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
v[i][0] *= factor;
v[i][1] *= factor;
v[i][2] *= factor;
}
}
/* ----------------------------------------------------------------------
zero the linear momentum of a group of atoms by adjusting v by -Vcm
------------------------------------------------------------------------- */
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
double masstotal = group->mass(igroup);
double vcm[3];
group->vcm(igroup,masstotal,vcm);
// adjust velocities by vcm to zero linear momentum
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
v[i][0] -= vcm[0];
v[i][1] -= vcm[1];
v[i][2] -= vcm[2];
}
}
/* ----------------------------------------------------------------------
zero the angular momentum of a group of atoms by adjusting v by -(w x r)
------------------------------------------------------------------------- */
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
double xcm[3],angmom[3],inertia[3][3],omega[3];
double masstotal = group->mass(igroup);
group->xcm(igroup,masstotal,xcm);
group->angmom(igroup,xcm,angmom);
group->inertia(igroup,xcm,inertia);
group->omega(angmom,inertia,omega);
// adjust velocities to zero omega
// vnew_i = v_i - w x r_i
// must use unwrapped coords to compute r_i correctly
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
tagint *image = atom->image;
int nlocal = atom->nlocal;
double dx,dy,dz;
double unwrap[3];
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;
v[i][1] -= omega[2]*dx - omega[0]*dz;
v[i][2] -= omega[0]*dy - omega[1]*dx;
}
}
/* ----------------------------------------------------------------------
parse optional parameters at end of velocity input line
------------------------------------------------------------------------- */
void Velocity::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR,"Illegal velocity command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"dist") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
if (strcmp(arg[iarg+1],"uniform") == 0) dist_flag = 0;
else if (strcmp(arg[iarg+1],"gaussian") == 0) dist_flag = 1;
else error->all(FLERR,"Illegal velocity command");
iarg += 2;
} else if (strcmp(arg[iarg],"sum") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
if (strcmp(arg[iarg+1],"no") == 0) sum_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) sum_flag = 1;
else error->all(FLERR,"Illegal velocity command");
iarg += 2;
} else if (strcmp(arg[iarg],"mom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
if (strcmp(arg[iarg+1],"no") == 0) momentum_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) momentum_flag = 1;
else error->all(FLERR,"Illegal velocity command");
iarg += 2;
} else if (strcmp(arg[iarg],"rot") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
if (strcmp(arg[iarg+1],"no") == 0) rotation_flag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) rotation_flag = 1;
else error->all(FLERR,"Illegal velocity command");
iarg += 2;
} else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
int icompute;
for (icompute = 0; icompute < modify->ncompute; icompute++)
if (strcmp(arg[iarg+1],modify->compute[icompute]->id) == 0) break;
if (icompute == modify->ncompute)
error->all(FLERR,"Could not find velocity temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,
"Velocity temperature ID does not compute temperature");
iarg += 2;
} else if (strcmp(arg[iarg],"loop") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
if (strcmp(arg[iarg+1],"all") == 0) loop_flag = ALL;
else if (strcmp(arg[iarg+1],"local") == 0) loop_flag = LOCAL;
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;
else if (strcmp(arg[iarg+1],"lattice") == 0) scale_flag = 1;
else error->all(FLERR,"Illegal velocity command");
iarg += 2;
} else error->all(FLERR,"Illegal velocity command");
}
}
diff --git a/src/velocity.h b/src/velocity.h
index 076dc1d34..d5e175bbc 100644
--- a/src/velocity.h
+++ b/src/velocity.h
@@ -1,140 +1,140 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMMAND_CLASS
CommandStyle(velocity,Velocity)
#else
#ifndef LMP_VELOCITY_H
#define LMP_VELOCITY_H
#include "pointers.h"
namespace LAMMPS_NS {
class Velocity : protected Pointers {
public:
Velocity(class LAMMPS *);
void command(int, char **);
void init_external(const char *);
void options(int, char **);
void create(double, int);
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;
void set(int, char **);
void scale(int, char **);
void ramp(int, char **);
void zero(int, char **);
void rescale(double, double);
void zero_momentum();
void zero_rotation();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Velocity command before simulation box is defined
The velocity command cannot be used before a read_data, read_restart,
or create_box command.
E: Velocity command with no atoms existing
A velocity command has been used, but no atoms yet exist.
E: Could not find velocity group ID
A group ID used in the velocity command does not exist.
W: Mismatch between velocity and compute groups
The temperature computation used by the velocity command will not be
on the same group of atoms that velocities are being set for.
E: Too big a problem to use velocity create loop all
The system size must fit in a 32-bit integer to use this option.
E: Cannot use velocity create loop all unless atoms have IDs
Atoms in the simulation to do not have IDs, so this style
of velocity creation cannot be performed.
E: Atom IDs must be consecutive for velocity create loop all
Self-explanatory.
E: Variable name for velocity set does not exist
Self-explanatory.
E: Variable for velocity set is invalid style
Only atom-style variables can be used.
E: Cannot set non-zero z velocity for 2d simulation
Self-explanatory.
E: Cannot set variable z velocity for 2d simulation
Self-explanatory.
E: Velocity ramp in z for a 2d problem
Self-explanatory.
E: Attempting to rescale a 0.0 temperature
Cannot rescale a temperature that is already 0.0.
E: Cannot zero momentum of 0 atoms
The collection of atoms for which momentum is being computed has no
atoms.
E: Could not find velocity temperature ID
The compute ID needed by the velocity command to compute temperature
does not exist.
E: Velocity temperature ID does not compute temperature
The compute ID given to the velocity command must compute
temperature.
U: Use of velocity with undefined lattice
If units = lattice (the default) for the velocity set or velocity ramp
command, then a lattice must first be defined via the lattice command.
*/