2006-09-28 03:51:33 +08:00
|
|
|
/* ----------------------------------------------------------------------
|
|
|
|
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
2007-01-30 08:22:05 +08:00
|
|
|
http://lammps.sandia.gov, Sandia National Laboratories
|
|
|
|
Steve Plimpton, sjplimp@sandia.gov
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
|
|
|
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
2012-06-07 06:47:51 +08:00
|
|
|
certain rights in this software. This software is distributed under
|
2006-09-28 03:51:33 +08:00
|
|
|
the GNU General Public License.
|
|
|
|
|
|
|
|
See the README file in the top-level LAMMPS directory.
|
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
#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"
|
2006-11-14 06:17:52 +08:00
|
|
|
#include "lattice.h"
|
2011-01-12 03:33:08 +08:00
|
|
|
#include "input.h"
|
|
|
|
#include "variable.h"
|
2006-09-28 03:51:33 +08:00
|
|
|
#include "force.h"
|
2007-01-30 08:22:05 +08:00
|
|
|
#include "modify.h"
|
2013-09-05 01:03:05 +08:00
|
|
|
#include "fix.h"
|
2007-01-30 08:22:05 +08:00
|
|
|
#include "compute.h"
|
|
|
|
#include "compute_temp.h"
|
2006-09-28 03:51:33 +08:00
|
|
|
#include "random_park.h"
|
|
|
|
#include "group.h"
|
|
|
|
#include "comm.h"
|
|
|
|
#include "memory.h"
|
|
|
|
#include "error.h"
|
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
using namespace LAMMPS_NS;
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
enum{CREATE,SET,SCALE,RAMP,ZERO};
|
|
|
|
enum{ALL,LOCAL,GEOM};
|
2011-01-12 03:33:08 +08:00
|
|
|
enum{NONE,CONSTANT,EQUAL,ATOM};
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
#define WARMUP 100
|
|
|
|
#define SMALL 0.001
|
|
|
|
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
Velocity::Velocity(LAMMPS *lmp) : Pointers(lmp) {}
|
|
|
|
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
void Velocity::command(int narg, char **arg)
|
|
|
|
{
|
2011-09-24 02:06:55 +08:00
|
|
|
if (narg < 2) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2012-06-07 06:47:51 +08:00
|
|
|
if (domain->box_exist == 0)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Velocity command before simulation box is defined");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (atom->natoms == 0)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Velocity command with no atoms existing");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2009-08-17 22:52:16 +08:00
|
|
|
// atom masses must all be set
|
|
|
|
|
|
|
|
atom->check_mass();
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// identify group
|
|
|
|
|
|
|
|
igroup = group->find(arg[0]);
|
2011-09-24 02:06:55 +08:00
|
|
|
if (igroup == -1) error->all(FLERR,"Could not find velocity group ID");
|
2006-09-28 03:51:33 +08:00
|
|
|
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;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// set defaults
|
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
temperature = NULL;
|
2006-09-28 03:51:33 +08:00
|
|
|
dist_flag = 0;
|
|
|
|
sum_flag = 0;
|
|
|
|
momentum_flag = 1;
|
|
|
|
rotation_flag = 0;
|
2015-02-07 04:39:07 +08:00
|
|
|
bias_flag = 0;
|
2006-09-28 03:51:33 +08:00
|
|
|
loop_flag = ALL;
|
|
|
|
scale_flag = 1;
|
2013-09-05 01:03:05 +08:00
|
|
|
rfix = -1;
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// 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]);
|
|
|
|
|
2015-03-04 23:24:51 +08:00
|
|
|
// special cases where full init and border communication must be done first
|
|
|
|
// for ZERO if fix rigid/small is used
|
2015-03-06 23:35:35 +08:00
|
|
|
// for CREATE/SET if compute temp/cs is used
|
2015-03-04 23:24:51 +08:00
|
|
|
// b/c methods invoked in the compute/fix perform forward/reverse comm
|
|
|
|
|
|
|
|
int initcomm = 0;
|
|
|
|
if (style == ZERO && rfix >= 0 &&
|
|
|
|
strcmp(modify->fix[rfix]->style,"rigid/small") == 0) initcomm = 1;
|
|
|
|
if ((style == CREATE || style == SET) && temperature &&
|
2015-03-06 23:35:35 +08:00
|
|
|
strcmp(temperature->style,"temp/cs") == 0) initcomm = 1;
|
2015-03-04 23:24:51 +08:00
|
|
|
|
|
|
|
if (initcomm) {
|
|
|
|
lmp->init();
|
|
|
|
if (domain->triclinic) domain->x2lamda(atom->nlocal);
|
|
|
|
domain->pbc();
|
|
|
|
domain->reset_box();
|
|
|
|
comm->setup();
|
|
|
|
comm->exchange();
|
|
|
|
comm->borders();
|
|
|
|
if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
|
|
|
|
}
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
// initialize velocities based on style
|
2009-08-18 23:27:33 +08:00
|
|
|
// create() invoked differently, so can be called externally
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2009-08-17 22:52:16 +08:00
|
|
|
if (style == CREATE) {
|
2013-06-29 03:00:58 +08:00
|
|
|
double t_desired = force->numeric(FLERR,arg[2]);
|
|
|
|
int seed = force->inumeric(FLERR,arg[3]);
|
2009-08-17 22:52:16 +08:00
|
|
|
create(t_desired,seed);
|
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
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]);
|
|
|
|
}
|
|
|
|
|
2009-08-17 22:52:16 +08:00
|
|
|
/* ----------------------------------------------------------------------
|
|
|
|
initialization of defaults before calling velocity methods externaly
|
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
2011-12-17 05:15:49 +08:00
|
|
|
void Velocity::init_external(const char *extgroup)
|
2009-08-17 22:52:16 +08:00
|
|
|
{
|
|
|
|
igroup = group->find(extgroup);
|
2011-09-24 02:06:55 +08:00
|
|
|
if (igroup == -1) error->all(FLERR,"Could not find velocity group ID");
|
2009-08-17 22:52:16 +08:00
|
|
|
groupbit = group->bitmask[igroup];
|
|
|
|
|
|
|
|
temperature = NULL;
|
|
|
|
dist_flag = 0;
|
|
|
|
sum_flag = 0;
|
|
|
|
momentum_flag = 1;
|
|
|
|
rotation_flag = 0;
|
|
|
|
loop_flag = ALL;
|
|
|
|
scale_flag = 1;
|
|
|
|
}
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
|
2009-08-17 22:52:16 +08:00
|
|
|
void Velocity::create(double t_desired, int seed)
|
2006-09-28 03:51:33 +08:00
|
|
|
{
|
|
|
|
int i;
|
2015-02-07 04:39:07 +08:00
|
|
|
double **vhold;
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2011-09-24 02:06:55 +08:00
|
|
|
if (seed <= 0) error->all(FLERR,"Illegal velocity create command");
|
2007-10-10 07:40:36 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
// if sum_flag set, store a copy of current velocities
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
if (sum_flag) {
|
|
|
|
double **v = atom->v;
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
memory->create(vhold,nlocal,3,"velocity:vhold");
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
|
|
vhold[i][0] = v[i][0];
|
|
|
|
vhold[i][1] = v[i][1];
|
|
|
|
vhold[i][2] = v[i][2];
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// if temperature = NULL or bias_flag set,
|
|
|
|
// create a new ComputeTemp with the velocity group
|
|
|
|
|
|
|
|
int tcreate_flag = 0;
|
|
|
|
Compute *temperature_nobias = NULL;
|
|
|
|
|
|
|
|
if (temperature == NULL || bias_flag) {
|
2006-09-28 03:51:33 +08:00
|
|
|
char **arg = new char*[3];
|
2007-10-04 00:22:30 +08:00
|
|
|
arg[0] = (char *) "velocity_temp";
|
2006-09-28 03:51:33 +08:00
|
|
|
arg[1] = group->names[igroup];
|
2007-10-04 00:22:30 +08:00
|
|
|
arg[2] = (char *) "temp";
|
2015-02-07 04:39:07 +08:00
|
|
|
if (temperature == NULL) {
|
|
|
|
temperature = new ComputeTemp(lmp,3,arg);
|
|
|
|
tcreate_flag = 1;
|
|
|
|
} else temperature_nobias = new ComputeTemp(lmp,3,arg);
|
2006-09-28 03:51:33 +08:00
|
|
|
delete [] arg;
|
2007-01-30 08:22:05 +08:00
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
// initialize temperature computation(s)
|
2006-09-28 03:51:33 +08:00
|
|
|
// warn if groups don't match
|
|
|
|
|
|
|
|
if (igroup != temperature->igroup && comm->me == 0)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->warning(FLERR,"Mismatch between velocity and compute groups");
|
2006-09-28 03:51:33 +08:00
|
|
|
temperature->init();
|
2013-02-13 07:43:26 +08:00
|
|
|
temperature->setup();
|
2015-02-07 04:39:07 +08:00
|
|
|
if (temperature_nobias) {
|
|
|
|
temperature_nobias->init();
|
|
|
|
temperature_nobias->setup();
|
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
// if bias_flag set, remove bias velocity from all atoms
|
|
|
|
// for some temperature computes, must first calculate temp to do that
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
if (bias_flag) {
|
|
|
|
temperature->compute_scalar();
|
|
|
|
temperature->remove_bias_all();
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
// 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)
|
2008-01-19 01:18:46 +08:00
|
|
|
// via random->reset()
|
2006-09-28 03:51:33 +08:00
|
|
|
// will always produce same V, independent of P
|
|
|
|
// adjust by factor for atom mass
|
2015-02-07 04:39:07 +08:00
|
|
|
// set xdim,ydim,zdim = 1/0 for whether to create velocity in those dims
|
|
|
|
// zdim = 0 for 2d
|
|
|
|
// any dims can be 0 if bias temperature compute turns them off
|
|
|
|
// currently only temp/partial does
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
double **v = atom->v;
|
2009-07-14 21:50:43 +08:00
|
|
|
double *rmass = atom->rmass;
|
|
|
|
double *mass = atom->mass;
|
2006-09-28 03:51:33 +08:00
|
|
|
int *type = atom->type;
|
|
|
|
int *mask = atom->mask;
|
2015-02-07 04:39:07 +08:00
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
int dim = domain->dimension;
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
int m;
|
|
|
|
double vx,vy,vz,factor;
|
|
|
|
RanPark *random;
|
2012-06-07 06:47:51 +08:00
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
if (loop_flag == ALL) {
|
|
|
|
|
|
|
|
// create an atom map if one doesn't exist already
|
|
|
|
|
|
|
|
int mapflag = 0;
|
|
|
|
if (atom->map_style == 0) {
|
|
|
|
mapflag = 1;
|
2008-11-07 23:47:56 +08:00
|
|
|
atom->nghost = 0;
|
2006-09-28 03:51:33 +08:00
|
|
|
atom->map_init();
|
|
|
|
atom->map_set();
|
|
|
|
}
|
|
|
|
|
2008-04-09 06:55:40 +08:00
|
|
|
// error check
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2011-01-07 07:23:37 +08:00
|
|
|
if (atom->natoms > MAXSMALLINT)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Too big a problem to use velocity create loop all");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (atom->tag_enable == 0)
|
2012-06-07 06:47:51 +08:00
|
|
|
error->all(FLERR,
|
|
|
|
"Cannot use velocity create loop all unless atoms have IDs");
|
2008-04-09 06:55:40 +08:00
|
|
|
if (atom->tag_consecutive() == 0)
|
2012-06-07 06:47:51 +08:00
|
|
|
error->all(FLERR,
|
|
|
|
"Atom IDs must be consecutive for velocity create loop all");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// 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
|
|
|
|
|
2008-04-09 06:55:40 +08:00
|
|
|
random = new RanPark(lmp,seed);
|
|
|
|
int natoms = static_cast<int> (atom->natoms);
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
for (i = 1; i <= natoms; i++) {
|
|
|
|
if (dist_flag == 0) {
|
2014-11-18 23:41:05 +08:00
|
|
|
vx = random->uniform() - 0.5;
|
|
|
|
vy = random->uniform() - 0.5;
|
|
|
|
vz = random->uniform() - 0.5;
|
2006-09-28 03:51:33 +08:00
|
|
|
} else {
|
2012-06-07 06:47:51 +08:00
|
|
|
vx = random->gaussian();
|
|
|
|
vy = random->gaussian();
|
|
|
|
vz = random->gaussian();
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
m = atom->map(i);
|
|
|
|
if (m >= 0 && m < nlocal) {
|
2012-06-07 06:47:51 +08:00
|
|
|
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;
|
2015-02-07 05:59:28 +08:00
|
|
|
if (dim == 3) v[m][2] = vz * factor;
|
2012-06-07 06:47:51 +08:00
|
|
|
else v[m][2] = 0.0;
|
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// delete temporary atom map
|
|
|
|
|
|
|
|
if (mapflag) {
|
|
|
|
atom->map_delete();
|
|
|
|
atom->map_style = 0;
|
|
|
|
}
|
|
|
|
|
|
|
|
} else if (loop_flag == LOCAL) {
|
2007-01-30 08:22:05 +08:00
|
|
|
random = new RanPark(lmp,seed + comm->me);
|
2006-09-28 03:51:33 +08:00
|
|
|
for (i = 0; i < WARMUP; i++) random->uniform();
|
|
|
|
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
|
|
if (mask[i] & groupbit) {
|
2012-06-07 06:47:51 +08:00
|
|
|
if (dist_flag == 0) {
|
2014-11-18 23:41:05 +08:00
|
|
|
vx = random->uniform() - 0.5;
|
|
|
|
vy = random->uniform() - 0.5;
|
|
|
|
vz = random->uniform() - 0.5;
|
2012-06-07 06:47:51 +08:00
|
|
|
} 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;
|
2015-02-07 04:39:07 +08:00
|
|
|
if (dim == 3) v[i][2] = vz * factor;
|
2012-06-07 06:47:51 +08:00
|
|
|
else v[i][2] = 0.0;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
} else if (loop_flag == GEOM) {
|
2007-06-20 21:17:59 +08:00
|
|
|
random = new RanPark(lmp,1);
|
2006-09-28 03:51:33 +08:00
|
|
|
double **x = atom->x;
|
2012-06-07 06:47:51 +08:00
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
|
|
if (mask[i] & groupbit) {
|
2012-06-07 06:47:51 +08:00
|
|
|
random->reset(seed,x[i]);
|
|
|
|
if (dist_flag == 0) {
|
2014-11-18 23:41:05 +08:00
|
|
|
vx = random->uniform() - 0.5;
|
|
|
|
vy = random->uniform() - 0.5;
|
|
|
|
vz = random->uniform() - 0.5;
|
2012-06-07 06:47:51 +08:00
|
|
|
} 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;
|
2015-02-07 04:39:07 +08:00
|
|
|
if (dim == 3) v[i][2] = vz * factor;
|
2012-06-07 06:47:51 +08:00
|
|
|
else v[i][2] = 0.0;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// apply momentum and rotation zeroing
|
|
|
|
|
|
|
|
if (momentum_flag) zero_momentum();
|
|
|
|
if (rotation_flag) zero_rotation();
|
|
|
|
|
|
|
|
// scale temp to desired value
|
2015-02-07 04:39:07 +08:00
|
|
|
// if bias flag is set, bias velocities have already been removed:
|
|
|
|
// no-bias compute calculates temp only for new thermal velocities
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
double t;
|
|
|
|
if (bias_flag == 0) t = temperature->compute_scalar();
|
|
|
|
else t = temperature_nobias->compute_scalar();
|
2006-09-28 03:51:33 +08:00
|
|
|
rescale(t,t_desired);
|
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
// if bias_flag set, restore bias velocity to all atoms
|
2015-03-05 23:57:12 +08:00
|
|
|
// reapply needed for temperature computes where velocity
|
|
|
|
// creation has messed up the bias that was already removed:
|
|
|
|
// compute temp/partial needs to reset v dims to 0.0
|
|
|
|
// compute temp/cs needs to reset v to COM velocity of each C/S pair
|
2015-02-07 04:39:07 +08:00
|
|
|
|
|
|
|
if (bias_flag) {
|
|
|
|
temperature->reapply_bias_all();
|
|
|
|
temperature->restore_bias_all();
|
|
|
|
}
|
|
|
|
|
2012-06-07 06:47:51 +08:00
|
|
|
// if sum_flag set, add back in previous velocities
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
if (sum_flag) {
|
|
|
|
for (i = 0; i < nlocal; i++) {
|
|
|
|
if (mask[i] & groupbit) {
|
2012-06-07 06:47:51 +08:00
|
|
|
v[i][0] += vhold[i][0];
|
|
|
|
v[i][1] += vhold[i][1];
|
|
|
|
v[i][2] += vhold[i][2];
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
}
|
2015-02-07 04:39:07 +08:00
|
|
|
memory->destroy(vhold);
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
// free local memory
|
2015-02-07 04:39:07 +08:00
|
|
|
// if temperature compute was created, delete it
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
delete random;
|
2015-02-07 04:39:07 +08:00
|
|
|
if (tcreate_flag) delete temperature;
|
|
|
|
if (temperature_nobias) delete temperature_nobias;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/* ---------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::set(int narg, char **arg)
|
|
|
|
{
|
2011-01-12 03:33:08 +08:00
|
|
|
int xstyle,ystyle,zstyle,varflag;
|
2006-09-28 03:51:33 +08:00
|
|
|
double vx,vy,vz;
|
2011-01-12 03:33:08 +08:00
|
|
|
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;
|
2013-06-29 03:00:58 +08:00
|
|
|
else vx = force->numeric(FLERR,arg[0]);
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
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;
|
2013-06-29 03:00:58 +08:00
|
|
|
else vy = force->numeric(FLERR,arg[1]);
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
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;
|
2013-06-29 03:00:58 +08:00
|
|
|
else vz = force->numeric(FLERR,arg[2]);
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
// set and apply scale factors
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2011-01-12 03:33:08 +08:00
|
|
|
xscale = yscale = zscale = 1.0;
|
|
|
|
|
|
|
|
if (xstyle && !xstr) {
|
|
|
|
if (scale_flag) xscale = domain->lattice->xlattice;
|
|
|
|
vx *= xscale;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
if (ystyle && !ystr) {
|
|
|
|
if (scale_flag) yscale = domain->lattice->ylattice;
|
|
|
|
vy *= yscale;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
if (zstyle && !zstr) {
|
|
|
|
if (scale_flag) zscale = domain->lattice->zlattice;
|
|
|
|
vz *= zscale;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
2011-01-12 03:33:08 +08:00
|
|
|
// check variables
|
|
|
|
|
|
|
|
if (xstr) {
|
|
|
|
xvar = input->variable->find(xstr);
|
2012-06-07 06:47:51 +08:00
|
|
|
if (xvar < 0)
|
|
|
|
error->all(FLERR,"Variable name for velocity set does not exist");
|
2011-01-12 03:33:08 +08:00
|
|
|
if (input->variable->equalstyle(xvar)) xstyle = EQUAL;
|
|
|
|
else if (input->variable->atomstyle(xvar)) xstyle = ATOM;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Variable for velocity set is invalid style");
|
2011-01-12 03:33:08 +08:00
|
|
|
}
|
|
|
|
if (ystr) {
|
|
|
|
yvar = input->variable->find(ystr);
|
2012-06-07 06:47:51 +08:00
|
|
|
if (yvar < 0)
|
|
|
|
error->all(FLERR,"Variable name for velocity set does not exist");
|
2011-01-12 03:33:08 +08:00
|
|
|
if (input->variable->equalstyle(yvar)) ystyle = EQUAL;
|
|
|
|
else if (input->variable->atomstyle(yvar)) ystyle = ATOM;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Variable for velocity set is invalid style");
|
2011-01-12 03:33:08 +08:00
|
|
|
}
|
|
|
|
if (zstr) {
|
|
|
|
zvar = input->variable->find(zstr);
|
2012-06-07 06:47:51 +08:00
|
|
|
if (zvar < 0)
|
|
|
|
error->all(FLERR,"Variable name for velocity set does not exist");
|
2011-01-12 03:33:08 +08:00
|
|
|
if (input->variable->equalstyle(zvar)) zstyle = EQUAL;
|
|
|
|
else if (input->variable->atomstyle(zvar)) zstyle = ATOM;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Variable for velocity set is invalid style");
|
2011-01-12 03:33:08 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
|
|
|
|
varflag = ATOM;
|
|
|
|
else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
|
|
|
|
varflag = EQUAL;
|
|
|
|
else varflag = CONSTANT;
|
|
|
|
|
2011-02-11 23:37:55 +08:00
|
|
|
// error check for 2d models
|
|
|
|
|
|
|
|
if (domain->dimension == 2) {
|
|
|
|
if (zstyle == CONSTANT && vz != 0.0)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Cannot set non-zero z velocity for 2d simulation");
|
2011-02-11 23:37:55 +08:00
|
|
|
if (zstyle == EQUAL || zstyle == ATOM)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Cannot set variable z velocity for 2d simulation");
|
2011-02-11 23:37:55 +08:00
|
|
|
}
|
|
|
|
|
2011-01-12 03:33:08 +08:00
|
|
|
// allocate vfield array if necessary
|
|
|
|
|
|
|
|
double **vfield = NULL;
|
2011-03-26 05:13:51 +08:00
|
|
|
if (varflag == ATOM) memory->create(vfield,atom->nlocal,3,"velocity:vfield");
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
// set velocities via constants
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
double **v = atom->v;
|
|
|
|
int *mask = atom->mask;
|
|
|
|
int nlocal = atom->nlocal;
|
|
|
|
|
2011-01-12 03:33:08 +08:00
|
|
|
if (varflag == CONSTANT) {
|
|
|
|
for (int i = 0; i < nlocal; i++) {
|
|
|
|
if (mask[i] & groupbit) {
|
2012-06-07 06:47:51 +08:00
|
|
|
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;
|
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
// set velocities via variables
|
|
|
|
|
|
|
|
} else {
|
|
|
|
if (xstyle == EQUAL) vx = input->variable->compute_equal(xvar);
|
2014-05-02 05:10:56 +08:00
|
|
|
else if (xstyle == ATOM) {
|
|
|
|
if (vfield) input->variable->compute_atom(xvar,igroup,&vfield[0][0],3,0);
|
|
|
|
else input->variable->compute_atom(xvar,igroup,NULL,3,0);
|
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
if (ystyle == EQUAL) vy = input->variable->compute_equal(yvar);
|
2014-05-02 05:10:56 +08:00
|
|
|
else if (ystyle == ATOM) {
|
|
|
|
if (vfield) input->variable->compute_atom(yvar,igroup,&vfield[0][1],3,0);
|
|
|
|
else input->variable->compute_atom(yvar,igroup,NULL,3,0);
|
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
if (zstyle == EQUAL) vz = input->variable->compute_equal(zvar);
|
2014-05-02 05:10:56 +08:00
|
|
|
else if (zstyle == ATOM) {
|
|
|
|
if (vfield) input->variable->compute_atom(zvar,igroup,&vfield[0][2],3,0);
|
|
|
|
else input->variable->compute_atom(zvar,igroup,NULL,3,0);
|
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
for (int i = 0; i < nlocal; i++)
|
|
|
|
if (mask[i] & groupbit) {
|
2012-06-07 06:47:51 +08:00
|
|
|
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;
|
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
2011-01-12 03:33:08 +08:00
|
|
|
|
|
|
|
// clean up
|
|
|
|
|
|
|
|
delete [] xstr;
|
|
|
|
delete [] ystr;
|
|
|
|
delete [] zstr;
|
2011-03-26 05:13:51 +08:00
|
|
|
memory->destroy(vfield);
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/* ----------------------------------------------------------------------
|
2012-06-07 06:47:51 +08:00
|
|
|
rescale velocities of a group after computing its temperature
|
2006-09-28 03:51:33 +08:00
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::scale(int narg, char **arg)
|
|
|
|
{
|
2013-06-29 03:00:58 +08:00
|
|
|
double t_desired = force->numeric(FLERR,arg[0]);
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
// if temperature = NULL, create a new ComputeTemp with the velocity group
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
int tflag = 0;
|
|
|
|
if (temperature == NULL) {
|
2006-09-28 03:51:33 +08:00
|
|
|
char **arg = new char*[3];
|
2007-10-04 00:22:30 +08:00
|
|
|
arg[0] = (char *) "velocity_temp";
|
2006-09-28 03:51:33 +08:00
|
|
|
arg[1] = group->names[igroup];
|
2007-10-04 00:22:30 +08:00
|
|
|
arg[2] = (char *) "temp";
|
2007-01-30 08:22:05 +08:00
|
|
|
temperature = new ComputeTemp(lmp,3,arg);
|
|
|
|
tflag = 1;
|
2006-09-28 03:51:33 +08:00
|
|
|
delete [] arg;
|
2007-01-30 08:22:05 +08:00
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// initialize temperature computation
|
|
|
|
// warn if groups don't match
|
|
|
|
|
|
|
|
if (igroup != temperature->igroup && comm->me == 0)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->warning(FLERR,"Mismatch between velocity and compute groups");
|
2006-09-28 03:51:33 +08:00
|
|
|
temperature->init();
|
2013-02-13 07:43:26 +08:00
|
|
|
temperature->setup();
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// scale temp to desired value
|
2015-02-07 04:39:07 +08:00
|
|
|
// if bias flag is set:
|
|
|
|
// temperature calculation will be done accounting for bias
|
|
|
|
// remove/restore bias velocities before/after rescale
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2015-02-07 04:39:07 +08:00
|
|
|
if (bias_flag == 0) {
|
|
|
|
double t = temperature->compute_scalar();
|
|
|
|
rescale(t,t_desired);
|
|
|
|
} else {
|
|
|
|
double t = temperature->compute_scalar();
|
|
|
|
temperature->remove_bias_all();
|
|
|
|
rescale(t,t_desired);
|
|
|
|
temperature->restore_bias_all();
|
|
|
|
}
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// if temperature was created, delete it
|
|
|
|
|
2007-01-30 08:22:05 +08:00
|
|
|
if (tflag) delete temperature;
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/* ----------------------------------------------------------------------
|
2012-06-07 06:47:51 +08:00
|
|
|
apply a ramped set of velocities
|
2006-09-28 03:51:33 +08:00
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::ramp(int narg, char **arg)
|
|
|
|
{
|
2011-01-12 03:33:08 +08:00
|
|
|
// 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
|
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
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;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2012-06-07 06:47:51 +08:00
|
|
|
if (v_dim == 2 && domain->dimension == 2)
|
2011-09-24 02:06:55 +08:00
|
|
|
error->all(FLERR,"Velocity ramp in z for a 2d problem");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
double v_lo,v_hi;
|
|
|
|
if (v_dim == 0) {
|
2013-06-29 03:00:58 +08:00
|
|
|
v_lo = xscale*force->numeric(FLERR,arg[1]);
|
|
|
|
v_hi = xscale*force->numeric(FLERR,arg[2]);
|
2006-09-28 03:51:33 +08:00
|
|
|
} else if (v_dim == 1) {
|
2013-06-29 03:00:58 +08:00
|
|
|
v_lo = yscale*force->numeric(FLERR,arg[1]);
|
|
|
|
v_hi = yscale*force->numeric(FLERR,arg[2]);
|
2010-04-29 21:42:52 +08:00
|
|
|
} else if (v_dim == 2) {
|
2013-06-29 03:00:58 +08:00
|
|
|
v_lo = zscale*force->numeric(FLERR,arg[1]);
|
|
|
|
v_hi = zscale*force->numeric(FLERR,arg[2]);
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
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;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
double coord_lo,coord_hi;
|
|
|
|
if (coord_dim == 0) {
|
2013-06-29 03:00:58 +08:00
|
|
|
coord_lo = xscale*force->numeric(FLERR,arg[4]);
|
|
|
|
coord_hi = xscale*force->numeric(FLERR,arg[5]);
|
2006-09-28 03:51:33 +08:00
|
|
|
} else if (coord_dim == 1) {
|
2013-06-29 03:00:58 +08:00
|
|
|
coord_lo = yscale*force->numeric(FLERR,arg[4]);
|
|
|
|
coord_hi = yscale*force->numeric(FLERR,arg[5]);
|
2006-09-28 03:51:33 +08:00
|
|
|
} else if (coord_dim == 2) {
|
2013-06-29 03:00:58 +08:00
|
|
|
coord_lo = zscale*force->numeric(FLERR,arg[4]);
|
|
|
|
coord_hi = zscale*force->numeric(FLERR,arg[5]);
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
// vramp = ramped velocity component for v_dim
|
|
|
|
// add or set based on sum_flag
|
2012-06-07 06:47:51 +08:00
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
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
|
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::zero(int narg, char **arg)
|
|
|
|
{
|
2013-09-05 01:03:05 +08:00
|
|
|
if (strcmp(arg[0],"linear") == 0) {
|
|
|
|
if (rfix < 0) zero_momentum();
|
2015-03-04 23:24:51 +08:00
|
|
|
else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
|
|
|
|
modify->fix[rfix]->setup_pre_neighbor();
|
|
|
|
modify->fix[rfix]->zero_momentum();
|
|
|
|
} else if (strstr(modify->fix[rfix]->style,"rigid")) {
|
|
|
|
modify->fix[rfix]->zero_momentum();
|
|
|
|
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
|
2013-09-05 01:03:05 +08:00
|
|
|
|
|
|
|
} else if (strcmp(arg[0],"angular") == 0) {
|
|
|
|
if (rfix < 0) zero_rotation();
|
2015-03-04 23:24:51 +08:00
|
|
|
else if (strcmp(modify->fix[rfix]->style,"rigid/small") == 0) {
|
|
|
|
modify->fix[rfix]->setup_pre_neighbor();
|
|
|
|
modify->fix[rfix]->zero_rotation();
|
|
|
|
} else if (strstr(modify->fix[rfix]->style,"rigid")) {
|
|
|
|
modify->fix[rfix]->zero_rotation();
|
|
|
|
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
|
2013-09-05 01:03:05 +08:00
|
|
|
|
|
|
|
} else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
/* ----------------------------------------------------------------------
|
2012-06-07 06:47:51 +08:00
|
|
|
rescale velocities of group atoms to t_new from t_old
|
2015-02-07 04:39:07 +08:00
|
|
|
no bias applied here, since done in create() and scale()
|
2006-09-28 03:51:33 +08:00
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::rescale(double t_old, double t_new)
|
|
|
|
{
|
2011-09-24 02:06:55 +08:00
|
|
|
if (t_old == 0.0) error->all(FLERR,"Attempting to rescale a 0.0 temperature");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
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()
|
|
|
|
{
|
2013-09-05 01:03:05 +08:00
|
|
|
// cannot have no atoms in group
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2011-01-05 05:08:27 +08:00
|
|
|
if (group->count(igroup) == 0)
|
2013-09-05 01:03:05 +08:00
|
|
|
error->all(FLERR,"Cannot zero momentum of no atoms");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// 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()
|
|
|
|
{
|
2007-01-30 08:22:05 +08:00
|
|
|
int i;
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2013-09-05 01:03:05 +08:00
|
|
|
// cannot have no atoms in group
|
2006-09-28 03:51:33 +08:00
|
|
|
|
2011-01-05 05:08:27 +08:00
|
|
|
if (group->count(igroup) == 0)
|
2013-09-05 01:03:05 +08:00
|
|
|
error->all(FLERR,"Cannot zero momentum of no atoms");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// 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);
|
2007-01-30 08:22:05 +08:00
|
|
|
group->omega(angmom,inertia,omega);
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
// 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;
|
2014-01-15 00:17:20 +08:00
|
|
|
imageint *image = atom->image;
|
2006-09-28 03:51:33 +08:00
|
|
|
int nlocal = atom->nlocal;
|
2012-06-07 06:47:51 +08:00
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
double dx,dy,dz;
|
2012-10-26 00:24:37 +08:00
|
|
|
double unwrap[3];
|
2012-06-07 06:47:51 +08:00
|
|
|
|
2006-09-28 03:51:33 +08:00
|
|
|
for (i = 0; i < nlocal; i++)
|
|
|
|
if (mask[i] & groupbit) {
|
2012-10-26 00:24:37 +08:00
|
|
|
domain->unmap(x[i],image[i],unwrap);
|
2013-09-05 01:03:05 +08:00
|
|
|
dx = unwrap[0] - xcm[0];
|
2012-10-26 00:24:37 +08:00
|
|
|
dy = unwrap[1] - xcm[1];
|
|
|
|
dz = unwrap[2] - xcm[2];
|
2006-09-28 03:51:33 +08:00
|
|
|
v[i][0] -= omega[1]*dz - omega[2]*dy;
|
2007-08-15 22:24:56 +08:00
|
|
|
v[i][1] -= omega[2]*dx - omega[0]*dz;
|
2006-09-28 03:51:33 +08:00
|
|
|
v[i][2] -= omega[0]*dy - omega[1]*dx;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
/* ----------------------------------------------------------------------
|
2012-06-07 06:47:51 +08:00
|
|
|
parse optional parameters at end of velocity input line
|
2006-09-28 03:51:33 +08:00
|
|
|
------------------------------------------------------------------------- */
|
|
|
|
|
|
|
|
void Velocity::options(int narg, char **arg)
|
|
|
|
{
|
2011-09-24 02:06:55 +08:00
|
|
|
if (narg < 0) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
|
|
|
|
int iarg = 0;
|
|
|
|
while (iarg < narg) {
|
|
|
|
if (strcmp(arg[iarg],"dist") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (strcmp(arg[iarg+1],"uniform") == 0) dist_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"gaussian") == 0) dist_flag = 1;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"sum") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (strcmp(arg[iarg+1],"no") == 0) sum_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"yes") == 0) sum_flag = 1;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"mom") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (strcmp(arg[iarg+1],"no") == 0) momentum_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"yes") == 0) momentum_flag = 1;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"rot") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (strcmp(arg[iarg+1],"no") == 0) rotation_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"yes") == 0) rotation_flag = 1;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
|
|
|
} else if (strcmp(arg[iarg],"temp") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2007-01-30 08:22:05 +08:00
|
|
|
int icompute;
|
|
|
|
for (icompute = 0; icompute < modify->ncompute; icompute++)
|
2012-06-07 06:47:51 +08:00
|
|
|
if (strcmp(arg[iarg+1],modify->compute[icompute]->id) == 0) break;
|
|
|
|
if (icompute == modify->ncompute)
|
|
|
|
error->all(FLERR,"Could not find velocity temperature ID");
|
2007-01-30 08:22:05 +08:00
|
|
|
temperature = modify->compute[icompute];
|
|
|
|
if (temperature->tempflag == 0)
|
2012-06-07 06:47:51 +08:00
|
|
|
error->all(FLERR,
|
|
|
|
"Velocity temperature ID does not compute temperature");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
2015-02-07 04:39:07 +08:00
|
|
|
} else if (strcmp(arg[iarg],"bias") == 0) {
|
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
|
|
|
if (strcmp(arg[iarg+1],"no") == 0) bias_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"yes") == 0) bias_flag = 1;
|
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
|
|
|
iarg += 2;
|
2006-09-28 03:51:33 +08:00
|
|
|
} else if (strcmp(arg[iarg],"loop") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
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;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
2013-09-05 01:03:05 +08:00
|
|
|
} 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;
|
2006-09-28 03:51:33 +08:00
|
|
|
} else if (strcmp(arg[iarg],"units") == 0) {
|
2011-09-24 02:06:55 +08:00
|
|
|
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
if (strcmp(arg[iarg+1],"box") == 0) scale_flag = 0;
|
|
|
|
else if (strcmp(arg[iarg+1],"lattice") == 0) scale_flag = 1;
|
2011-09-24 02:06:55 +08:00
|
|
|
else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
iarg += 2;
|
2011-09-24 02:06:55 +08:00
|
|
|
} else error->all(FLERR,"Illegal velocity command");
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|
2015-02-07 04:39:07 +08:00
|
|
|
|
|
|
|
// error check
|
|
|
|
|
|
|
|
if (bias_flag && temperature == NULL)
|
|
|
|
error->all(FLERR,"Cannot use velocity bias command without temp keyword");
|
|
|
|
if (bias_flag && temperature->tempbias == 0)
|
|
|
|
error->all(FLERR,"Velocity temperature ID does calculate a velocity bias");
|
2006-09-28 03:51:33 +08:00
|
|
|
}
|