From 544dde674bfbd8dd7004cd345f084819cad8bc1a Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 12 Oct 2010 17:58:40 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5023
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/COLLOID/fix_wall_colloid.cpp | 10 +-
src/COLLOID/fix_wall_colloid.h | 2 +-
src/fix_gravity.cpp | 4 +-
src/fix_wall.cpp | 220 ++++++-------
src/fix_wall.h | 12 +-
src/fix_wall_harmonic.cpp | 10 +-
src/fix_wall_harmonic.h | 2 +-
src/fix_wall_lj126.cpp | 10 +-
src/fix_wall_lj126.h | 2 +-
src/fix_wall_lj93.cpp | 10 +-
src/fix_wall_lj93.h | 2 +-
src/fix_wall_reflect.cpp | 282 ++++++----------
src/fix_wall_reflect.h | 9 +-
src/variable.cpp | 533 +++++++++++++++++++++----------
src/variable.h | 3 +
15 files changed, 608 insertions(+), 503 deletions(-)
diff --git a/src/COLLOID/fix_wall_colloid.cpp b/src/COLLOID/fix_wall_colloid.cpp
index 8fcef726c6..4859378bcf 100644
--- a/src/COLLOID/fix_wall_colloid.cpp
+++ b/src/COLLOID/fix_wall_colloid.cpp
@@ -84,11 +84,13 @@ void FixWallColloid::precompute(int m)
}
/* ----------------------------------------------------------------------
- interaction of all particles in group with all 6 walls (if defined)
+ interaction of all particles in group with a wall
+ m = index of wall coeffs
+ which = xlo,xhi,ylo,yhi,zlo,zhi
error if any finite-size particle is touching or penetrating wall
------------------------------------------------------------------------- */
-void FixWallColloid::wall_particle(int m, double coord)
+void FixWallColloid::wall_particle(int m, int which, double coord)
{
double delta,delta2,rinv,r2inv,r4inv,r8inv,fwall;
double r2,rinv2,r2inv2,r4inv2,r6inv2;
@@ -102,8 +104,8 @@ void FixWallColloid::wall_particle(int m, double coord)
int *type = atom->type;
int nlocal = atom->nlocal;
- int dim = m/2;
- int side = m % 2;
+ int dim = which / 2;
+ int side = which % 2;
if (side == 0) side = -1;
int onflag = 0;
diff --git a/src/COLLOID/fix_wall_colloid.h b/src/COLLOID/fix_wall_colloid.h
index 1d09975e18..4b1734cf78 100644
--- a/src/COLLOID/fix_wall_colloid.h
+++ b/src/COLLOID/fix_wall_colloid.h
@@ -29,7 +29,7 @@ class FixWallColloid : public FixWall {
FixWallColloid(class LAMMPS *, int, char **);
void init();
void precompute(int);
- void wall_particle(int, double);
+ void wall_particle(int, int, double);
private:
double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp
index 8070ccc51a..7d7c362875 100644
--- a/src/fix_gravity.cpp
+++ b/src/fix_gravity.cpp
@@ -62,8 +62,8 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
zdir = atof(arg[7]);
} else error->all("Illegal fix gravity command");
- double PI = 2.0 * asin(1.0);
- degree2rad = 2.0*PI / 360.0;
+ double PI = 4.0*atan(1.0);
+ degree2rad = PI/180.0;
if (style == CHUTE || style == SPHERICAL || style == GRADIENT) {
if (domain->dimension == 3) {
diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp
index fee8d27fa2..78632c99d8 100644
--- a/src/fix_wall.cpp
+++ b/src/fix_wall.cpp
@@ -16,16 +16,18 @@
#include "string.h"
#include "fix_wall.h"
#include "atom.h"
+#include "input.h"
+#include "variable.h"
#include "domain.h"
#include "lattice.h"
#include "update.h"
-#include "output.h"
#include "respa.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{XLO,XHI,YLO,YHI,ZLO,ZHI};
+enum{NONE,EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
@@ -34,17 +36,13 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
{
scalar_flag = 1;
vector_flag = 1;
- size_vector = 6;
global_freq = 1;
extscalar = 1;
extvector = 1;
- time_depend = 1;
// parse args
- for (int m = 0; m < 6; m++) wallflag[m] = 0;
- velflag = 0;
- wigflag = 0;
+ nwall = 0;
int scaleflag = 1;
int iarg = 3;
@@ -53,39 +51,42 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
(strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
(strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
if (iarg+5 > narg) error->all("Illegal fix wall command");
- int m;
- if (strcmp(arg[iarg],"xlo") == 0) m = XLO;
- else if (strcmp(arg[iarg],"xhi") == 0) m = XHI;
- else if (strcmp(arg[iarg],"ylo") == 0) m = YLO;
- else if (strcmp(arg[iarg],"yhi") == 0) m = YHI;
- else if (strcmp(arg[iarg],"zlo") == 0) m = ZLO;
- else if (strcmp(arg[iarg],"zhi") == 0) m = ZHI;
- wallflag[m] = 1;
- coord0[m] = atof(arg[iarg+1]);
- epsilon[m] = atof(arg[iarg+2]);
- sigma[m] = atof(arg[iarg+3]);
- cutoff[m] = atof(arg[iarg+4]);
+
+ int newwall;
+ if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
+ else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
+ else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
+ else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
+ else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
+ else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
+
+ for (int m = 0; m < nwall; m++)
+ if (newwall == wallwhich[m])
+ error->all("Wall defined twice in fix wall command");
+
+ wallwhich[nwall] = newwall;
+ if (strcmp(arg[iarg+1],"EDGE") == 0) {
+ wallstyle[nwall] = EDGE;
+ int dim = wallwhich[nwall] / 2;
+ int side = wallwhich[nwall] % 2;
+ if (side == 0) coord0[nwall] = domain->boxlo[dim];
+ else coord0[nwall] = domain->boxhi[dim];
+ } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
+ wallstyle[nwall] = VARIABLE;
+ int n = strlen(&arg[iarg+1][2]) + 1;
+ varstr[nwall] = new char[n];
+ strcpy(varstr[nwall],&arg[iarg+1][2]);
+ } else {
+ wallstyle[nwall] = CONSTANT;
+ coord0[nwall] = atof(arg[iarg+1]);
+ }
+
+ epsilon[nwall] = atof(arg[iarg+2]);
+ sigma[nwall] = atof(arg[iarg+3]);
+ cutoff[nwall] = atof(arg[iarg+4]);
+ nwall++;
iarg += 5;
- } else if (strcmp(arg[iarg],"vel") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall command");
- velflag = 1;
- double vtmp = atof(arg[iarg+1]);
- for (int m = 0; m < 6; m++) vel[m] = vtmp;
- iarg += 2;
- } else if (strcmp(arg[iarg],"wiggle/sin") == 0) {
- if (iarg+3 > narg) error->all("Illegal fix wall command");
- wigflag = 1;
- double atmp = atof(arg[iarg+1]);
- period = atof(arg[iarg+2]);
- for (int m = 0; m < 6; m++) amplitude[m] = atmp;
- iarg += 3;
- } else if (strcmp(arg[iarg],"wiggle/cos") == 0) {
- if (iarg+3 > narg) error->all("Illegal fix wall command");
- wigflag = 2;
- double atmp = atof(arg[iarg+1]);
- period = atof(arg[iarg+2]);
- for (int m = 0; m < 6; m++) amplitude[m] = atmp;
- iarg += 3;
+
} else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all("Illegal fix wall command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
@@ -95,66 +96,70 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix wall command");
}
+ size_vector = nwall;
+
// error check
- int flag = 0;
- for (int m = 0; m < 6; m++) if (wallflag[m]) flag = 1;
- if (!flag) error->all("Illegal fix wall command");
-
- for (int m = 0; m < 6; m++)
- if (wallflag[m] && cutoff[m] <= 0.0)
+ if (nwall == 0) error->all("Illegal fix wall command");
+ for (int m = 0; m < nwall; m++)
+ if (cutoff[m] <= 0.0)
error->all("Fix wall cutoff <= 0.0");
- if (velflag && wigflag)
- error->all("Cannot set both vel and wiggle in fix wall command");
-
- if ((wallflag[XLO] || wallflag[XHI]) && domain->xperiodic)
- error->all("Cannot use fix wall in periodic dimension");
- if ((wallflag[YLO] || wallflag[YHI]) && domain->yperiodic)
- error->all("Cannot use fix wall in periodic dimension");
- if ((wallflag[ZLO] || wallflag[ZHI]) && domain->zperiodic)
- error->all("Cannot use fix wall in periodic dimension");
-
- if ((wallflag[ZLO] || wallflag[ZHI]) && domain->dimension == 2)
- error->all("Cannot use fix wall zlo/zhi for a 2d simulation");
-
- // setup scaling
-
- if (scaleflag && domain->lattice == NULL)
- error->all("Use of fix wall with undefined lattice");
-
- double xscale,yscale,zscale;
- if (scaleflag) {
- xscale = domain->lattice->xlattice;
- yscale = domain->lattice->ylattice;
- zscale = domain->lattice->zlattice;
- }
- else xscale = yscale = zscale = 1.0;
-
- // apply scaling factors to coord0, vel, amplitude
-
- double scale;
- for (int m = 0; m < 6; m++) {
- if (wallflag[m] == 0) continue;
- if (m < 2) scale = xscale;
- else if (m < 4) scale = yscale;
- else scale = zscale;
- coord0[m] *= scale;
- if (velflag) vel[m] *= scale;
- if (wigflag) amplitude[m] *= scale;
+ for (int m = 0; m < nwall; m++) {
+ if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
+ error->all("Cannot use fix wall in periodic dimension");
+ if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
+ error->all("Cannot use fix wall in periodic dimension");
+ if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
+ error->all("Cannot use fix wall in periodic dimension");
}
- // setup oscillations
+ for (int m = 0; m < nwall; m++)
+ if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
+ error->all("Cannot use fix wall zlo/zhi for a 2d simulation");
+
+ // scale coord for CONSTANT walls
- if (wigflag) {
- double PI = 4.0 * atan(1.0);
- omega = 2.0*PI / period;
+ int flag = 0;
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == CONSTANT) flag = 1;
+
+ if (flag) {
+ if (scaleflag && domain->lattice == NULL)
+ error->all("Use of fix wall with undefined lattice");
+
+ double xscale,yscale,zscale;
+ if (scaleflag) {
+ xscale = domain->lattice->xlattice;
+ yscale = domain->lattice->ylattice;
+ zscale = domain->lattice->zlattice;
+ }
+ else xscale = yscale = zscale = 1.0;
+
+ double scale;
+ for (int m = 0; m < nwall; m++) {
+ if (wallwhich[m] < YLO) scale = xscale;
+ else if (wallwhich[m] < ZLO) scale = yscale;
+ else scale = zscale;
+ if (wallstyle[m] == CONSTANT) coord0[m] *= scale;
+ }
}
+ // set time_depend if any wall positions are variable
+
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == VARIABLE) time_depend = 1;
+
eflag = 0;
- for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+ for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
+}
- time_origin = update->ntimestep;
+/* ---------------------------------------------------------------------- */
+
+FixWall::~FixWall()
+{
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == VARIABLE) delete [] varstr[m];
}
/* ---------------------------------------------------------------------- */
@@ -175,10 +180,18 @@ void FixWall::init()
{
dt = update->dt;
+ for (int m = 0; m < nwall; m++) {
+ if (wallstyle[m] != VARIABLE) continue;
+ varindex[m] = input->variable->find(varstr[m]);
+ if (varindex[m] < 0)
+ error->all("Variable name for fix wall does not exist");
+ if (!input->variable->equalstyle(varindex[m]))
+ error->all("Variable for fix wall is invalid style");
+ }
+
// setup coefficients
- for (int m = 0; m < 6; m++)
- if (wallflag[m]) precompute(m);
+ for (int m = 0; m < nwall; m++) precompute(m);
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
@@ -209,29 +222,18 @@ void FixWall::min_setup(int vflag)
void FixWall::post_force(int vflag)
{
eflag = 0;
- for (int m = 0; m < 7; m++) ewall[m] = 0.0;
+ for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
- // coord0 = initial position of wall
// coord = current position of wall
+ // evaluate variable if necessary
- double delta = (update->ntimestep - time_origin) * dt;
double coord;
+ for (int m = 0; m < nwall; m++) {
+ if (wallstyle[m] == VARIABLE)
+ coord = input->variable->compute_equal(varindex[m]);
+ else coord = coord0[m];
- for (int m = 0; m < 6; m++) {
- if (wallflag[m] == 0) continue;
-
- if (velflag) {
- if (m % 2 == 0) coord = coord0[m] + delta*vel[m];
- else coord = coord0[m] - delta*vel[m];
- } else if (wigflag == 1) {
- if (m % 2 == 0) coord = coord0[m] + amplitude[m]*sin(omega*delta);
- else coord = coord0[m] - amplitude[m]*sin(omega*delta);
- } else if (wigflag == 2) {
- if (m % 2 == 0) coord = coord0[m] + amplitude[m]*(1.0-cos(omega*delta));
- else coord = coord0[m] - amplitude[m]*(1.0-cos(omega*delta));
- } else coord = coord0[m];
-
- wall_particle(m,coord);
+ wall_particle(m,wallwhich[m],coord);
}
}
@@ -258,7 +260,7 @@ double FixWall::compute_scalar()
// only sum across procs one time
if (eflag == 0) {
- MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+ MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
eflag = 1;
}
return ewall_all[0];
@@ -273,7 +275,7 @@ double FixWall::compute_vector(int n)
// only sum across procs one time
if (eflag == 0) {
- MPI_Allreduce(ewall,ewall_all,7,MPI_DOUBLE,MPI_SUM,world);
+ MPI_Allreduce(ewall,ewall_all,nwall+1,MPI_DOUBLE,MPI_SUM,world);
eflag = 1;
}
return ewall_all[n+1];
diff --git a/src/fix_wall.h b/src/fix_wall.h
index dc7c26963a..04abd315f2 100644
--- a/src/fix_wall.h
+++ b/src/fix_wall.h
@@ -21,7 +21,7 @@ namespace LAMMPS_NS {
class FixWall : public Fix {
public:
FixWall(class LAMMPS *, int, char **);
- virtual ~FixWall() {}
+ virtual ~FixWall();
int setmask();
virtual void init();
void setup(int);
@@ -33,14 +33,14 @@ class FixWall : public Fix {
double compute_vector(int);
virtual void precompute(int) = 0;
- virtual void wall_particle(int, double) = 0;
+ virtual void wall_particle(int, int, double) = 0;
protected:
- int wallflag[6];
+ int nwall;
+ int wallwhich[6],wallstyle[6];
double coord0[6],epsilon[6],sigma[6],cutoff[6];
- int velflag,wigflag;
- double vel[6],amplitude[6];
- double period,omega;
+ char *varstr[6];
+ int varindex[6];
int eflag;
double ewall[7],ewall_all[7];
int nlevels_respa;
diff --git a/src/fix_wall_harmonic.cpp b/src/fix_wall_harmonic.cpp
index d8b18aebd5..4d0660a0f4 100644
--- a/src/fix_wall_harmonic.cpp
+++ b/src/fix_wall_harmonic.cpp
@@ -24,11 +24,13 @@ FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) :
FixWall(lmp, narg, arg) {}
/* ----------------------------------------------------------------------
- interaction of all particles in group with all 6 walls (if defined)
+ interaction of all particles in group with a wall
+ m = index of wall coeffs
+ which = xlo,xhi,ylo,yhi,zlo,zhi
error if any particle is on or behind wall
------------------------------------------------------------------------- */
-void FixWallHarmonic::wall_particle(int m, double coord)
+void FixWallHarmonic::wall_particle(int m, int which, double coord)
{
double delta,dr,fwall;
@@ -37,8 +39,8 @@ void FixWallHarmonic::wall_particle(int m, double coord)
int *mask = atom->mask;
int nlocal = atom->nlocal;
- int dim = m/2;
- int side = m % 2;
+ int dim = which / 2;
+ int side = which % 2;
if (side == 0) side = -1;
int onflag = 0;
diff --git a/src/fix_wall_harmonic.h b/src/fix_wall_harmonic.h
index 391a2880e3..a5aa575e01 100644
--- a/src/fix_wall_harmonic.h
+++ b/src/fix_wall_harmonic.h
@@ -28,7 +28,7 @@ class FixWallHarmonic : public FixWall {
public:
FixWallHarmonic(class LAMMPS *, int, char **);
void precompute(int) {}
- void wall_particle(int, double);
+ void wall_particle(int, int, double);
private:
double offset[6];
diff --git a/src/fix_wall_lj126.cpp b/src/fix_wall_lj126.cpp
index 1379676a49..132cea9c0b 100644
--- a/src/fix_wall_lj126.cpp
+++ b/src/fix_wall_lj126.cpp
@@ -38,11 +38,13 @@ void FixWallLJ126::precompute(int m)
}
/* ----------------------------------------------------------------------
- interaction of all particles in group with all 6 walls (if defined)
+ interaction of all particles in group with a wall
+ m = index of wall coeffs
+ which = xlo,xhi,ylo,yhi,zlo,zhi
error if any particle is on or behind wall
------------------------------------------------------------------------- */
-void FixWallLJ126::wall_particle(int m, double coord)
+void FixWallLJ126::wall_particle(int m, int which, double coord)
{
double delta,rinv,r2inv,r6inv,fwall;
@@ -51,8 +53,8 @@ void FixWallLJ126::wall_particle(int m, double coord)
int *mask = atom->mask;
int nlocal = atom->nlocal;
- int dim = m/2;
- int side = m % 2;
+ int dim = which / 2;
+ int side = which % 2;
if (side == 0) side = -1;
int onflag = 0;
diff --git a/src/fix_wall_lj126.h b/src/fix_wall_lj126.h
index 1a4f17db47..eba0ffa3e5 100644
--- a/src/fix_wall_lj126.h
+++ b/src/fix_wall_lj126.h
@@ -28,7 +28,7 @@ class FixWallLJ126 : public FixWall {
public:
FixWallLJ126(class LAMMPS *, int, char **);
void precompute(int);
- void wall_particle(int, double);
+ void wall_particle(int, int, double);
private:
double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
diff --git a/src/fix_wall_lj93.cpp b/src/fix_wall_lj93.cpp
index 99f5f90903..e106ec5e6a 100644
--- a/src/fix_wall_lj93.cpp
+++ b/src/fix_wall_lj93.cpp
@@ -39,11 +39,13 @@ void FixWallLJ93::precompute(int m)
}
/* ----------------------------------------------------------------------
- interaction of all particles in group with all 6 walls (if defined)
+ interaction of all particles in group with a wall
+ m = index of wall coeffs
+ which = xlo,xhi,ylo,yhi,zlo,zhi
error if any particle is on or behind wall
------------------------------------------------------------------------- */
-void FixWallLJ93::wall_particle(int m, double coord)
+void FixWallLJ93::wall_particle(int m, int which, double coord)
{
double delta,rinv,r2inv,r4inv,r10inv,fwall;
@@ -52,8 +54,8 @@ void FixWallLJ93::wall_particle(int m, double coord)
int *mask = atom->mask;
int nlocal = atom->nlocal;
- int dim = m/2;
- int side = m % 2;
+ int dim = which / 2;
+ int side = which % 2;
if (side == 0) side = -1;
int onflag = 0;
diff --git a/src/fix_wall_lj93.h b/src/fix_wall_lj93.h
index 3a6cd77957..78df357c64 100644
--- a/src/fix_wall_lj93.h
+++ b/src/fix_wall_lj93.h
@@ -28,7 +28,7 @@ class FixWallLJ93 : public FixWall {
public:
FixWallLJ93(class LAMMPS *, int, char **);
void precompute(int);
- void wall_particle(int, double);
+ void wall_particle(int, int, double);
private:
double coeff1[6],coeff2[6],coeff3[6],coeff4[6],offset[6];
diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp
index ba9aea7b4a..56bb32e270 100644
--- a/src/fix_wall_reflect.cpp
+++ b/src/fix_wall_reflect.cpp
@@ -25,6 +25,7 @@
using namespace LAMMPS_NS;
+enum{XLO,XHI,YLO,YHI,ZLO,ZHI};
enum{NONE,EDGE,CONSTANT,VARIABLE};
/* ---------------------------------------------------------------------- */
@@ -34,101 +35,48 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
{
if (narg < 4) error->all("Illegal fix wall/reflect command");
- xloflag = xhiflag = yloflag = yhiflag = zloflag = zhiflag = NONE;
- xlostr = xhistr = ylostr = yhistr = zlostr = zhistr = NULL;
+ // parse args
+
+ nwall = 0;
int scaleflag = 1;
int iarg = 3;
while (iarg < narg) {
- if (strcmp(arg[iarg],"xlo") == 0) {
+ if ((strcmp(arg[iarg],"xlo") == 0) || (strcmp(arg[iarg],"xhi") == 0) ||
+ (strcmp(arg[iarg],"ylo") == 0) || (strcmp(arg[iarg],"yhi") == 0) ||
+ (strcmp(arg[iarg],"zlo") == 0) || (strcmp(arg[iarg],"zhi") == 0)) {
if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
+
+ int newwall;
+ if (strcmp(arg[iarg],"xlo") == 0) newwall = XLO;
+ else if (strcmp(arg[iarg],"xhi") == 0) newwall = XHI;
+ else if (strcmp(arg[iarg],"ylo") == 0) newwall = YLO;
+ else if (strcmp(arg[iarg],"yhi") == 0) newwall = YHI;
+ else if (strcmp(arg[iarg],"zlo") == 0) newwall = ZLO;
+ else if (strcmp(arg[iarg],"zhi") == 0) newwall = ZHI;
+
+ for (int m = 0; m < nwall; m++)
+ if (newwall == wallwhich[m])
+ error->all("Wall defined twice in fix wall/reflect command");
+
+ wallwhich[nwall] = newwall;
if (strcmp(arg[iarg+1],"EDGE") == 0) {
- xloflag = EDGE;
- xlo = domain->boxlo[0];
+ wallstyle[nwall] = EDGE;
+ int dim = wallwhich[nwall] / 2;
+ int side = wallwhich[nwall] % 2;
+ if (side == 0) coord0[nwall] = domain->boxlo[dim];
+ else coord0[nwall] = domain->boxhi[dim];
} else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- xloflag = VARIABLE;
+ wallstyle[nwall] = VARIABLE;
int n = strlen(&arg[iarg+1][2]) + 1;
- xlostr = new char[n];
- strcpy(xlostr,&arg[iarg+1][2]);
+ varstr[nwall] = new char[n];
+ strcpy(varstr[nwall],&arg[iarg+1][2]);
} else {
- xloflag = CONSTANT;
- xlo = atof(arg[iarg+1]);
- }
- iarg += 2;
- } else if (strcmp(arg[iarg],"xhi") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
- if (strcmp(arg[iarg+1],"EDGE") == 0) {
- xhiflag = EDGE;
- xhi = domain->boxhi[0];
- } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- xhiflag = VARIABLE;
- int n = strlen(&arg[iarg+1][2]) + 1;
- xhistr = new char[n];
- strcpy(xhistr,&arg[iarg+1][2]);
- } else {
- xhiflag = CONSTANT;
- xhi = atof(arg[iarg+1]);
- }
- iarg += 2;
- } else if (strcmp(arg[iarg],"ylo") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
- if (strcmp(arg[iarg+1],"EDGE") == 0) {
- yloflag = EDGE;
- ylo = domain->boxlo[1];
- } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- yloflag = VARIABLE;
- int n = strlen(&arg[iarg+1][2]) + 1;
- ylostr = new char[n];
- strcpy(ylostr,&arg[iarg+1][2]);
- } else {
- yloflag = CONSTANT;
- ylo = atof(arg[iarg+1]);
- }
- iarg += 2;
- } else if (strcmp(arg[iarg],"yhi") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
- if (strcmp(arg[iarg+1],"EDGE") == 0) {
- yhiflag = EDGE;
- yhi = domain->boxhi[1];
- } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- yhiflag = VARIABLE;
- int n = strlen(&arg[iarg+1][2]) + 1;
- yhistr = new char[n];
- strcpy(yhistr,&arg[iarg+1][2]);
- } else {
- yhiflag = CONSTANT;
- yhi = atof(arg[iarg+1]);
- }
- iarg += 2;
- } else if (strcmp(arg[iarg],"zlo") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
- if (strcmp(arg[iarg+1],"EDGE") == 0) {
- zloflag = EDGE;
- zlo = domain->boxlo[2];
- } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- zloflag = VARIABLE;
- int n = strlen(&arg[iarg+1][2]) + 1;
- zlostr = new char[n];
- strcpy(zlostr,&arg[iarg+1][2]);
- } else {
- zloflag = CONSTANT;
- zlo = atof(arg[iarg+1]);
- }
- iarg += 2;
- } else if (strcmp(arg[iarg],"zhi") == 0) {
- if (iarg+2 > narg) error->all("Illegal fix wall/reflect command");
- if (strcmp(arg[iarg+1],"EDGE") == 0) {
- zhiflag = EDGE;
- zhi = domain->boxhi[2];
- } else if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
- zhiflag = VARIABLE;
- int n = strlen(&arg[iarg+1][2]) + 1;
- zhistr = new char[n];
- strcpy(zhistr,&arg[iarg+1][2]);
- } else {
- zhiflag = CONSTANT;
- zhi = atof(arg[iarg+1]);
+ wallstyle[nwall] = CONSTANT;
+ coord0[nwall] = atof(arg[iarg+1]);
}
+
+ nwall++;
iarg += 2;
} else if (strcmp(arg[iarg],"units") == 0) {
@@ -140,15 +88,32 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
} else error->all("Illegal fix wall/reflect command");
}
- // setup scaling if any wall positions were set as CONSTANT
- // apply scaling factors to wall positions
+ // error check
- if (xloflag == CONSTANT || xhiflag == CONSTANT ||
- yloflag == CONSTANT || yhiflag == CONSTANT ||
- zloflag == CONSTANT || zhiflag == CONSTANT) {
+ if (nwall == 0) error->all("Illegal fix wall command");
+ for (int m = 0; m < nwall; m++) {
+ if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic)
+ error->all("Cannot use fix wall/reflect in periodic dimension");
+ if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic)
+ error->all("Cannot use fix wall/reflect in periodic dimension");
+ if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic)
+ error->all("Cannot use fix wall/reflect in periodic dimension");
+ }
+
+ for (int m = 0; m < nwall; m++)
+ if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2)
+ error->all("Cannot use fix wall/reflect zlo/zhi for a 2d simulation");
+
+ // scale coord for CONSTANT walls
+
+ int flag = 0;
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == CONSTANT) flag = 1;
+
+ if (flag) {
if (scaleflag && domain->lattice == NULL)
- error->all("Use of fix wall/reflect with undefined lattice");
+ error->all("Use of fix wall with undefined lattice");
double xscale,yscale,zscale;
if (scaleflag) {
@@ -158,37 +123,27 @@ FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) :
}
else xscale = yscale = zscale = 1.0;
- if (xloflag == CONSTANT) xlo *= xscale;
- if (xhiflag == CONSTANT) xhi *= xscale;
- if (yloflag == CONSTANT) ylo *= xscale;
- if (yhiflag == CONSTANT) yhi *= xscale;
- if (zloflag == CONSTANT) zlo *= xscale;
- if (zhiflag == CONSTANT) zhi *= xscale;
+ double scale;
+ for (int m = 0; m < nwall; m++) {
+ if (wallwhich[m] < YLO) scale = xscale;
+ else if (wallwhich[m] < ZLO) scale = yscale;
+ else scale = zscale;
+ if (wallstyle[m] == CONSTANT) coord0[m] *= scale;
+ }
}
- // error checks
+ // set time_depend if any wall positions are variable
- if (!xloflag && !xhiflag && !yloflag && !yhiflag && !zloflag && !zhiflag)
- error->all("Illegal fix wall/reflect command");
-
- if ((xloflag || xhiflag) && domain->xperiodic)
- error->all("Cannot use wall in periodic dimension");
- if ((yloflag || yhiflag) && domain->yperiodic)
- error->all("Cannot use wall in periodic dimension");
- if ((zloflag || zhiflag) && domain->zperiodic)
- error->all("Cannot use wall in periodic dimension");
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == VARIABLE) time_depend = 1;
}
/* ---------------------------------------------------------------------- */
FixWallReflect::~FixWallReflect()
{
- delete [] xlostr;
- delete [] xhistr;
- delete [] ylostr;
- delete [] yhistr;
- delete [] zlostr;
- delete [] zhistr;
+ for (int m = 0; m < nwall; m++)
+ if (wallstyle[m] == VARIABLE) delete [] varstr[m];
}
/* ---------------------------------------------------------------------- */
@@ -205,46 +160,12 @@ int FixWallReflect::setmask()
void FixWallReflect::init()
{
- if (xlostr) {
- xlovar = input->variable->find(xlostr);
- if (xlovar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(xlovar))
- error->all("Variable for fix wall/reflect is invalid style");
- }
- if (xhistr) {
- xhivar = input->variable->find(xhistr);
- if (xhivar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(xhivar))
- error->all("Variable for fix wall/reflect is invalid style");
- }
- if (ylostr) {
- ylovar = input->variable->find(ylostr);
- if (ylovar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(ylovar))
- error->all("Variable for fix wall/reflect is invalid style");
- }
- if (yhistr) {
- yhivar = input->variable->find(yhistr);
- if (yhivar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(yhivar))
- error->all("Variable for fix wall/reflect is invalid style");
- }
- if (zlostr) {
- zlovar = input->variable->find(zlostr);
- if (zlovar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(zlovar))
- error->all("Variable for fix wall/reflect is invalid style");
- }
- if (zhistr) {
- zhivar = input->variable->find(zhistr);
- if (zhivar < 0)
- error->all("Variable name for fix wall/relect does not exist");
- if (!input->variable->equalstyle(zhivar))
+ for (int m = 0; m < nwall; m++) {
+ if (wallstyle[m] != VARIABLE) continue;
+ varindex[m] = input->variable->find(varstr[m]);
+ if (varindex[m] < 0)
+ error->all("Variable name for fix wall/reflect does not exist");
+ if (!input->variable->equalstyle(varindex[m]))
error->all("Variable for fix wall/reflect is invalid style");
}
@@ -261,44 +182,35 @@ void FixWallReflect::init()
void FixWallReflect::post_integrate()
{
- if (xloflag == VARIABLE) xlo = input->variable->compute_equal(xlovar);
- if (xhiflag == VARIABLE) xhi = input->variable->compute_equal(xhivar);
- if (yloflag == VARIABLE) ylo = input->variable->compute_equal(ylovar);
- if (yhiflag == VARIABLE) yhi = input->variable->compute_equal(yhivar);
- if (zloflag == VARIABLE) zlo = input->variable->compute_equal(zlovar);
- if (zhiflag == VARIABLE) zhi = input->variable->compute_equal(zhivar);
+ int i,dim,side;
+ double coord;
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
- for (int i = 0; i < nlocal; i++) {
- if (mask[i] & groupbit) {
- if (xloflag && x[i][0] < xlo) {
- x[i][0] = xlo + (xlo - x[i][0]);
- v[i][0] = -v[i][0];
+ for (int m = 0; m < nwall; m++) {
+ if (wallstyle[m] == VARIABLE)
+ coord = input->variable->compute_equal(varindex[m]);
+ else coord = coord0[m];
+
+ dim = wallwhich[m] / 2;
+ side = wallwhich[m] % 2;
+
+ for (i = 0; i < nlocal; i++)
+ if (mask[i] & groupbit) {
+ if (side == 0) {
+ if (x[i][dim] < coord) {
+ x[i][dim] = coord + (coord - x[i][dim]);
+ v[i][dim] = -v[i][dim];
+ }
+ } else {
+ if (x[i][0] > coord) {
+ x[i][dim] = coord - (x[i][dim] - coord);
+ v[i][dim] = -v[i][dim];
+ }
+ }
}
- if (xhiflag && x[i][0] > xhi) {
- x[i][0] = xhi - (x[i][0] - xhi);
- v[i][0] = -v[i][0];
- }
- if (yloflag && x[i][1] < ylo) {
- x[i][1] = ylo + (ylo - x[i][1]);
- v[i][1] = -v[i][1];
- }
- if (yhiflag && x[i][1] > yhi) {
- x[i][1] = yhi - (x[i][1] - yhi);
- v[i][1] = -v[i][1];
- }
- if (zloflag && x[i][2] < zlo) {
- x[i][2] = zlo + (zlo - x[i][2]);
- v[i][2] = -v[i][2];
- }
- if (zhiflag && x[i][2] > zhi) {
- x[i][2] = zhi - (x[i][2] - zhi);
- v[i][2] = -v[i][2];
- }
- }
}
}
diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h
index f0f7ef41fb..b0bda7b465 100644
--- a/src/fix_wall_reflect.h
+++ b/src/fix_wall_reflect.h
@@ -33,10 +33,11 @@ class FixWallReflect : public Fix {
void post_integrate();
private:
- int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag;
- char *xlostr,*xhistr,*ylostr,*yhistr,*zlostr,*zhistr;
- int xlovar,xhivar,ylovar,yhivar,zlovar,zhivar;
- double xlo,xhi,ylo,yhi,zlo,zhi;
+ int nwall;
+ int wallwhich[6],wallstyle[6];
+ double coord0[6];
+ char *varstr[6];
+ int varindex[6];
};
}
diff --git a/src/variable.cpp b/src/variable.cpp
index a943000735..5e1aba8270 100644
--- a/src/variable.cpp
+++ b/src/variable.cpp
@@ -43,11 +43,18 @@ using namespace LAMMPS_NS;
enum{INDEX,LOOP,WORLD,UNIVERSE,ULOOP,STRING,EQUAL,ATOM};
enum{ARG,OP};
+
+// customize by adding a math function
+
enum{DONE,ADD,SUBTRACT,MULTIPLY,DIVIDE,CARAT,UNARY,
EQ,NE,LT,LE,GT,GE,AND,OR,
SQRT,EXP,LN,LOG,SIN,COS,TAN,ASIN,ACOS,ATAN,ATAN2,
- RANDOM,NORMAL,CEIL,FLOOR,ROUND,
+ RANDOM,NORMAL,CEIL,FLOOR,ROUND,RAMP,STAGGER,LOGFREQ,
+ VLINEAR,SWIGGLE,CWIGGLE,
VALUE,ATOMARRAY,TYPEARRAY,INTARRAY};
+
+// customize by adding a special function
+
enum{SUM,XMIN,XMAX,AVE,TRAP};
#define INVOKED_SCALAR 1
@@ -83,6 +90,8 @@ Variable::Variable(LAMMPS *lmp) : Pointers(lmp)
precedence[MULTIPLY] = precedence[DIVIDE] = 6;
precedence[CARAT] = 7;
precedence[UNARY] = 8;
+
+ PI = 4.0*atan(1.0);
}
/* ---------------------------------------------------------------------- */
@@ -595,6 +604,7 @@ void Variable::copy(int narg, char **from, char **to)
sqrt(x),exp(x),ln(x),log(x),
sin(x),cos(x),tan(x),asin(x),atan2(y,x),...
group function = count(group), mass(group), xcm(group,x), ...
+ special function = sum(x),min(x), ...
atom value = x[i], y[i], vx[i], ...
atom vector = x, y, vx, ...
compute = c_ID, c_ID[i], c_ID[i][j]
@@ -686,10 +696,10 @@ double Variable::evaluate(char *str, Tree **tree)
// ----------------
// letter: c_ID, c_ID[], c_ID[][], f_ID, f_ID[], f_ID[][],
- // v_name, v_name[], exp(), xcm(,), x, x[], vol
+ // v_name, v_name[], exp(), xcm(,), x, x[], PI, vol
// ----------------
- } else if (islower(onechar)) {
+ } else if (isalpha(onechar)) {
if (expect == OP) error->all("Invalid syntax in variable formula");
expect = OP;
@@ -1134,7 +1144,8 @@ double Variable::evaluate(char *str, Tree **tree)
delete [] id;
// ----------------
- // math/group/special function or atom value/vector or thermo keyword
+ // math/group/special function or atom value/vector or
+ // constant or thermo keyword
// ----------------
} else {
@@ -1183,6 +1194,20 @@ double Variable::evaluate(char *str, Tree **tree)
atom_vector(word,tree,treestack,ntreestack);
+ // ----------------
+ // constant
+ // ----------------
+
+ } else if (is_constant(word)) {
+ value1 = constant(word);
+ if (tree) {
+ Tree *newtree = new Tree();
+ newtree->type = VALUE;
+ newtree->value = value1;
+ newtree->left = newtree->middle = newtree->right = NULL;
+ treestack[ntreestack++] = newtree;
+ } else argstack[nargstack++] = value1;
+
// ----------------
// thermo keyword
// ----------------
@@ -1350,7 +1375,9 @@ double Variable::evaluate(char *str, Tree **tree)
process an evaulation tree
customize by adding a math function:
sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
- atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round()
+ atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
+ ramp(x,y),stagger(x,y),logfreq(x,y,z),
+ vlinear(x,y),swiggle(x,y,z),cwiggle(x,y,z)
---------------------------------------------------------------------- */
double Variable::eval_tree(Tree *tree, int i)
@@ -1368,12 +1395,12 @@ double Variable::eval_tree(Tree *tree, int i)
return eval_tree(tree->left,i) * eval_tree(tree->right,i);
if (tree->type == DIVIDE) {
double denom = eval_tree(tree->right,i);
- if (denom == 0.0) error->all("Divide by 0 in variable formula");
+ if (denom == 0.0) error->one("Divide by 0 in variable formula");
return eval_tree(tree->left,i) / denom;
}
if (tree->type == CARAT) {
double exponent = eval_tree(tree->right,i);
- if (exponent == 0.0) error->all("Power by 0 in variable formula");
+ if (exponent == 0.0) error->one("Power by 0 in variable formula");
return pow(eval_tree(tree->left,i),exponent);
}
if (tree->type == UNARY)
@@ -1416,7 +1443,7 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == SQRT) {
double arg = eval_tree(tree->left,i);
- if (arg < 0.0) error->all("Sqrt of negative value in variable formula");
+ if (arg < 0.0) error->one("Sqrt of negative value in variable formula");
return sqrt(arg);
}
if (tree->type == EXP)
@@ -1424,13 +1451,13 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == LN) {
double arg = eval_tree(tree->left,i);
if (arg <= 0.0)
- error->all("Log of zero/negative value in variable formula");
+ error->one("Log of zero/negative value in variable formula");
return log(arg);
}
if (tree->type == LOG) {
double arg = eval_tree(tree->left,i);
if (arg <= 0.0)
- error->all("Log of zero/negative value in variable formula");
+ error->one("Log of zero/negative value in variable formula");
return log10(arg);
}
@@ -1444,13 +1471,13 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == ASIN) {
double arg = eval_tree(tree->left,i);
if (arg < -1.0 || arg > 1.0)
- error->all("Arcsin of invalid value in variable formula");
+ error->one("Arcsin of invalid value in variable formula");
return asin(arg);
}
if (tree->type == ACOS) {
double arg = eval_tree(tree->left,i);
if (arg < -1.0 || arg > 1.0)
- error->all("Arccos of invalid value in variable formula");
+ error->one("Arccos of invalid value in variable formula");
return acos(arg);
}
if (tree->type == ATAN)
@@ -1463,6 +1490,7 @@ double Variable::eval_tree(Tree *tree, int i)
double upper = eval_tree(tree->middle,i);
if (randomatom == NULL) {
int seed = static_cast<int> (eval_tree(tree->right,i));
+ if (seed <= 0) error->one("Invalid math function in variable formula");
randomatom = new RanMars(lmp,seed+me);
}
return randomatom->uniform()*(upper-lower)+lower;
@@ -1470,11 +1498,13 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == NORMAL) {
double mu = eval_tree(tree->left,i);
double sigma = eval_tree(tree->middle,i);
+ if (sigma < 0.0) error->one("Invalid math function in variable formula");
if (randomatom == NULL) {
int seed = static_cast<int> (eval_tree(tree->right,i));
+ if (seed <= 0) error->one("Invalid math function in variable formula");
randomatom = new RanMars(lmp,seed+me);
}
- return randomatom->gaussian()*sigma + mu;
+ return mu + sigma*randomatom->gaussian();
}
if (tree->type == CEIL)
@@ -1484,6 +1514,76 @@ double Variable::eval_tree(Tree *tree, int i)
if (tree->type == ROUND)
return MYROUND(eval_tree(tree->left,i));
+ if (tree->type == RAMP) {
+ double arg1 = eval_tree(tree->left,i);
+ double arg2 = eval_tree(tree->right,i);
+ double delta = update->ntimestep - update->beginstep;
+ delta /= update->endstep - update->beginstep;
+ double arg = arg1 + delta*(arg2-arg1);
+ return arg;
+ }
+
+ if (tree->type == STAGGER) {
+ int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
+ int ivalue2 = static_cast<int> (eval_tree(tree->right,i));
+ if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
+ error->one("Invalid math function in variable formula");
+ int lower = update->ntimestep/ivalue1 * ivalue1;
+ int delta = update->ntimestep - lower;
+ double arg;
+ if (delta < ivalue2) arg = lower+ivalue2;
+ else arg = lower+ivalue1;
+ return arg;
+ }
+
+ if (tree->type == LOGFREQ) {
+ int ivalue1 = static_cast<int> (eval_tree(tree->left,i));
+ int ivalue2 = static_cast<int> (eval_tree(tree->middle,i));
+ int ivalue3 = static_cast<int> (eval_tree(tree->right,i));
+ if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+ error->one("Invalid math function in variable formula");
+ double arg;
+ if (update->ntimestep < ivalue1) arg = ivalue1;
+ else {
+ int lower = ivalue1;
+ while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
+ int multiple = update->ntimestep/lower;
+ if (multiple < ivalue2) arg = (multiple+1)*lower;
+ else arg = lower*ivalue3;
+ }
+ return arg;
+ }
+
+ if (tree->type == VLINEAR) {
+ double arg1 = eval_tree(tree->left,i);
+ double arg2 = eval_tree(tree->right,i);
+ double delta = update->ntimestep - update->beginstep;
+ double arg = arg1 + arg2*delta*update->dt;
+ return arg;
+ }
+
+ if (tree->type == SWIGGLE) {
+ double arg1 = eval_tree(tree->left,i);
+ double arg2 = eval_tree(tree->middle,i);
+ double arg3 = eval_tree(tree->right,i);
+ if (arg3 == 0.0) error->one("Invalid math function in variable formula");
+ double delta = update->ntimestep - update->beginstep;
+ double omega = 2.0*PI/arg3;
+ double arg = arg1 + arg2*sin(omega*delta*update->dt);
+ return arg;
+ }
+
+ if (tree->type == CWIGGLE) {
+ double arg1 = eval_tree(tree->left,i);
+ double arg2 = eval_tree(tree->middle,i);
+ double arg3 = eval_tree(tree->right,i);
+ if (arg3 == 0.0) error->one("Invalid math function in variable formula");
+ double delta = update->ntimestep - update->beginstep;
+ double omega = 2.0*PI/arg3;
+ double arg = arg1 + arg2*(1.0-cos(omega*delta*update->dt));
+ return arg;
+ }
+
return 0.0;
}
@@ -1563,9 +1663,11 @@ int Variable::int_between_brackets(char *&ptr)
word = math function
contents = str between parentheses with one,two,three args
return 0 if not a match, 1 if successfully processed
- customize by adding a math function in 2 places:
+ customize by adding a math function:
sqrt(),exp(),ln(),log(),sin(),cos(),tan(),asin(),acos(),atan(),
- atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round()
+ atan2(y,x),random(x,y,z),normal(x,y,z),ceil(),floor(),round(),
+ ramp(x,y),stagger(x,y),logfreq(x,y,z),
+ vlinear(x,y),swiggle(x,y,z),cwiggle(x,y,z)
------------------------------------------------------------------------- */
int Variable::math_function(char *word, char *contents, Tree **tree,
@@ -1581,7 +1683,10 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
strcmp(word,"acos") && strcmp(word,"atan") &&
strcmp(word,"atan2") && strcmp(word,"random") &&
strcmp(word,"normal") && strcmp(word,"ceil") &&
- strcmp(word,"floor") && strcmp(word,"round"))
+ strcmp(word,"floor") && strcmp(word,"round") &&
+ strcmp(word,"ramp") && strcmp(word,"stagger") &&
+ strcmp(word,"logfreq") && strcmp(word,"vlinear") &&
+ strcmp(word,"swiggle") && strcmp(word,"cwiggle"))
return 0;
// parse contents for arg1,arg2,arg3 separated by commas
@@ -1728,6 +1833,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
else {
if (randomequal == NULL) {
int seed = static_cast<int> (value3);
+ if (seed <= 0) error->all("Invalid math function in variable formula");
randomequal = new RanMars(lmp,seed);
}
argstack[nargstack++] = randomequal->uniform()*(value2-value1) + value1;
@@ -1736,11 +1842,14 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
if (narg != 3) error->all("Invalid math function in variable formula");
if (tree) newtree->type = NORMAL;
else {
+ if (value2 < 0.0)
+ error->all("Invalid math function in variable formula");
if (randomequal == NULL) {
int seed = static_cast<int> (value3);
+ if (seed <= 0) error->all("Invalid math function in variable formula");
randomequal = new RanMars(lmp,seed);
}
- argstack[nargstack++] = randomequal->gaussian()*value2 + value1;
+ argstack[nargstack++] = value1 + value2*randomequal->gaussian();
}
} else if (strcmp(word,"ceil") == 0) {
@@ -1757,6 +1866,94 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
if (narg != 1) error->all("Invalid math function in variable formula");
if (tree) newtree->type = ROUND;
else argstack[nargstack++] = MYROUND(value1);
+
+ } else if (strcmp(word,"ramp") == 0) {
+ if (narg != 2) error->all("Invalid math function in variable formula");
+ if (update->whichflag == 0)
+ error->all("Cannot use ramp in variable formula between runs");
+ if (tree) newtree->type = RAMP;
+ else {
+ double delta = update->ntimestep - update->beginstep;
+ delta /= update->endstep - update->beginstep;
+ double value = value1 + delta*(value2-value1);
+ argstack[nargstack++] = value;
+ }
+
+ } else if (strcmp(word,"stagger") == 0) {
+ if (narg != 2) error->all("Invalid math function in variable formula");
+ if (tree) newtree->type = STAGGER;
+ else {
+ int ivalue1 = static_cast<int> (value1);
+ int ivalue2 = static_cast<int> (value2);
+ if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
+ error->all("Invalid math function in variable formula");
+ int lower = update->ntimestep/ivalue1 * ivalue1;
+ int delta = update->ntimestep - lower;
+ double value;
+ if (delta < ivalue2) value = lower+ivalue2;
+ else value = lower+ivalue1;
+ argstack[nargstack++] = value;
+ }
+
+ } else if (strcmp(word,"logfreq") == 0) {
+ if (narg != 3) error->all("Invalid math function in variable formula");
+ if (tree) newtree->type = LOGFREQ;
+ else {
+ int ivalue1 = static_cast<int> (value1);
+ int ivalue2 = static_cast<int> (value2);
+ int ivalue3 = static_cast<int> (value3);
+ if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
+ error->all("Invalid math function in variable formula");
+ double value;
+ if (update->ntimestep < ivalue1) value = ivalue1;
+ else {
+ int lower = ivalue1;
+ while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
+ int multiple = update->ntimestep/lower;
+ if (multiple < ivalue2) value = (multiple+1)*lower;
+ else value = lower*ivalue3;
+ }
+ argstack[nargstack++] = value;
+ }
+
+ } else if (strcmp(word,"vlinear") == 0) {
+ if (narg != 2) error->all("Invalid math function in variable formula");
+ if (update->whichflag == 0)
+ error->all("Cannot use vlinear in variable formula between runs");
+ if (tree) newtree->type = VLINEAR;
+ else {
+ double delta = update->ntimestep - update->beginstep;
+ double value = value1 + value2*delta*update->dt;
+ argstack[nargstack++] = value;
+ }
+
+ } else if (strcmp(word,"swiggle") == 0) {
+ if (narg != 3) error->all("Invalid math function in variable formula");
+ if (update->whichflag == 0)
+ error->all("Cannot use swiggle in variable formula between runs");
+ if (tree) newtree->type = CWIGGLE;
+ else {
+ if (value3 == 0.0)
+ error->all("Invalid math function in variable formula");
+ double delta = update->ntimestep - update->beginstep;
+ double omega = 2.0*PI/value3;
+ double value = value1 + value2*sin(omega*delta*update->dt);
+ argstack[nargstack++] = value;
+ }
+
+ } else if (strcmp(word,"cwiggle") == 0) {
+ if (narg != 3) error->all("Invalid math function in variable formula");
+ if (update->whichflag == 0)
+ error->all("Cannot use cwiggle in variable formula between runs");
+ if (tree) newtree->type = CWIGGLE;
+ else {
+ if (value3 == 0.0)
+ error->all("Invalid math function in variable formula");
+ double delta = update->ntimestep - update->beginstep;
+ double omega = 2.0*PI/value3;
+ double value = value1 + value2*(1.0-cos(omega*delta*update->dt));
+ argstack[nargstack++] = value;
+ }
}
delete [] arg1;
@@ -1772,7 +1969,7 @@ int Variable::math_function(char *word, char *contents, Tree **tree,
word = group function
contents = str between parentheses with one,two,three args
return 0 if not a match, 1 if successfully processed
- customize by adding a group function in 2 places with optional region arg:
+ customize by adding a group function with optional region arg:
count(group),mass(group),charge(group),
xcm(group,dim),vcm(group,dim),fcm(group,dim),
bound(group,xmin),gyration(group),ke(group),angmom(group),
@@ -2017,9 +2214,8 @@ int Variable::region_function(char *id)
word = special function
contents = str between parentheses with one,two,three args
return 0 if not a match, 1 if successfully processed
- customize by adding a special function in 2 places
- ramp(x,y),stagger(x,y),logfreq(x,y,z),sum(x),min(x),max(x),
- ave(x),trap(x)
+ customize by adding a special function:
+ sum(x),min(x),max(x),ave(x),trap(x)
------------------------------------------------------------------------- */
int Variable::special_function(char *word, char *contents, Tree **tree,
@@ -2028,9 +2224,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
{
// word not a match to any special function
- if (strcmp(word,"ramp") && strcmp(word,"stagger") &&
- strcmp(word,"logfreq") && strcmp(word,"sum") &&
- strcmp(word,"min") && strcmp(word,"max") &&
+ if (strcmp(word,"sum") && strcmp(word,"min") && strcmp(word,"max") &&
strcmp(word,"ave") && strcmp(word,"trap"))
return 0;
@@ -2066,163 +2260,122 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
// match word to special function
- double value;
+ int method;
+ if (strcmp(word,"sum") == 0) method = SUM;
+ else if (strcmp(word,"min") == 0) method = XMIN;
+ else if (strcmp(word,"max") == 0) method = XMAX;
+ else if (strcmp(word,"ave") == 0) method = AVE;
+ else if (strcmp(word,"trap") == 0) method = TRAP;
- if (strcmp(word,"ramp") == 0) {
- if (narg != 2) error->all("Invalid special function in variable formula");
- if (update->whichflag == 0)
- error->all("Cannot use ramp in variable formula between runs");
- double value1 = numeric(arg1);
- double value2 = numeric(arg2);
- double delta = update->ntimestep - update->beginstep;
- delta /= update->endstep - update->beginstep;
- value = value1 + delta*(value2-value1);
+ if (narg != 1) error->all("Invalid special function in variable formula");
+
+ Compute *compute = NULL;
+ Fix *fix = NULL;
+ int index,nvec,nstride;
- } else if (strcmp(word,"stagger") == 0) {
- if (narg != 2) error->all("Invalid special function in variable formula");
- int ivalue1 = inumeric(arg1);
- int ivalue2 = inumeric(arg2);
- if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue1 <= ivalue2)
- error->all("Invalid special function in variable formula");
- int lower = update->ntimestep/ivalue1 * ivalue1;
- int delta = update->ntimestep - lower;
- if (delta < ivalue2) value = lower+ivalue2;
- else value = lower+ivalue1;
+ if (strstr(arg1,"c_") == arg1) {
+ ptr1 = strchr(arg1,'[');
+ if (ptr1) {
+ ptr2 = ptr1;
+ index = int_between_brackets(ptr2);
+ *ptr1 = '\0';
+ } else index = 0;
- } else if (strcmp(word,"logfreq") == 0) {
- if (narg != 3) error->all("Invalid special function in variable formula");
- int ivalue1 = inumeric(arg1);
- int ivalue2 = inumeric(arg2);
- int ivalue3 = inumeric(arg3);
- if (ivalue1 <= 0 || ivalue2 <= 0 || ivalue3 <= 0 || ivalue2 >= ivalue3)
- error->all("Invalid special function in variable formula");
- if (update->ntimestep < ivalue1) value = ivalue1;
- else {
- int lower = ivalue1;
- while (update->ntimestep >= ivalue3*lower) lower *= ivalue3;
- int multiple = update->ntimestep/lower;
- if (multiple < ivalue2) value = (multiple+1)*lower;
- else value = lower*ivalue3;
- }
-
- } else {
- if (narg != 1) error->all("Invalid special function in variable formula");
-
- int method;
- if (strcmp(word,"sum") == 0) method = SUM;
- else if (strcmp(word,"min") == 0) method = XMIN;
- else if (strcmp(word,"max") == 0) method = XMAX;
- else if (strcmp(word,"ave") == 0) method = AVE;
- else if (strcmp(word,"trap") == 0) method = TRAP;
-
- Compute *compute = NULL;
- Fix *fix = NULL;
- int index,nvec,nstride;
-
- if (strstr(arg1,"c_") == arg1) {
- ptr1 = strchr(arg1,'[');
- if (ptr1) {
- ptr2 = ptr1;
- index = int_between_brackets(ptr2);
- *ptr1 = '\0';
- } else index = 0;
-
- int icompute = modify->find_compute(&arg1[2]);
- if (icompute < 0) error->all("Invalid compute ID in variable formula");
- compute = modify->compute[icompute];
- if (index == 0 && compute->vector_flag) {
- if (update->whichflag == 0) {
- if (compute->invoked_vector != update->ntimestep)
- error->all("Compute used in variable between runs is not current");
- } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
- compute->compute_vector();
- compute->invoked_flag |= INVOKED_VECTOR;
- }
- nvec = compute->size_vector;
- nstride = 1;
- } else if (index && compute->array_flag) {
- if (index > compute->size_array_cols)
- error->all("Variable formula compute array "
- "is accessed out-of-range");
- if (update->whichflag == 0) {
- if (compute->invoked_array != update->ntimestep)
- error->all("Compute used in variable between runs is not current");
- } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
- compute->compute_array();
- compute->invoked_flag |= INVOKED_ARRAY;
- }
- nvec = compute->size_array_rows;
- nstride = compute->size_array_cols;
- } else error->all("Mismatched compute in variable formula");
-
- } else if (strstr(arg1,"f_") == arg1) {
- ptr1 = strchr(arg1,'[');
- if (ptr1) {
- ptr2 = ptr1;
- index = int_between_brackets(ptr2);
- *ptr1 = '\0';
- } else index = 0;
-
- int ifix = modify->find_fix(&arg1[2]);
- if (ifix < 0) error->all("Invalid fix ID in variable formula");
- fix = modify->fix[ifix];
- if (index == 0 && fix->vector_flag) {
- if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
- error->all("Fix in variable not computed at compatible time");
- nvec = fix->size_vector;
- nstride = 1;
- } else if (index && fix->array_flag) {
- if (index > fix->size_array_cols)
- error->all("Variable formula fix array is accessed out-of-range");
- if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
- error->all("Fix in variable not computed at compatible time");
- nvec = fix->size_array_rows;
- nstride = fix->size_array_cols;
- } else error->all("Mismatched fix in variable formula");
-
- } else error->all("Invalid special function in variable formula");
+ int icompute = modify->find_compute(&arg1[2]);
+ if (icompute < 0) error->all("Invalid compute ID in variable formula");
+ compute = modify->compute[icompute];
+ if (index == 0 && compute->vector_flag) {
+ if (update->whichflag == 0) {
+ if (compute->invoked_vector != update->ntimestep)
+ error->all("Compute used in variable between runs is not current");
+ } else if (!(compute->invoked_flag & INVOKED_VECTOR)) {
+ compute->compute_vector();
+ compute->invoked_flag |= INVOKED_VECTOR;
+ }
+ nvec = compute->size_vector;
+ nstride = 1;
+ } else if (index && compute->array_flag) {
+ if (index > compute->size_array_cols)
+ error->all("Variable formula compute array "
+ "is accessed out-of-range");
+ if (update->whichflag == 0) {
+ if (compute->invoked_array != update->ntimestep)
+ error->all("Compute used in variable between runs is not current");
+ } else if (!(compute->invoked_flag & INVOKED_ARRAY)) {
+ compute->compute_array();
+ compute->invoked_flag |= INVOKED_ARRAY;
+ }
+ nvec = compute->size_array_rows;
+ nstride = compute->size_array_cols;
+ } else error->all("Mismatched compute in variable formula");
- value = 0.0;
- if (method == XMIN) value = BIG;
- if (method == XMAX) value = -BIG;
-
- if (compute) {
- double *vec;
- if (index) vec = &compute->array[0][index];
- else vec = compute->vector;
-
- int j = 0;
- for (int i = 0; i < nvec; i++) {
- if (method == SUM) value += vec[j];
- else if (method == XMIN) value = MIN(value,vec[j]);
- else if (method == XMAX) value = MAX(value,vec[j]);
- else if (method == AVE) value += vec[j];
- else if (method == TRAP) {
- if (i > 0 && i < nvec-1) value += vec[j];
- else value += 0.5*vec[j];
- }
- j += nstride;
+ } else if (strstr(arg1,"f_") == arg1) {
+ ptr1 = strchr(arg1,'[');
+ if (ptr1) {
+ ptr2 = ptr1;
+ index = int_between_brackets(ptr2);
+ *ptr1 = '\0';
+ } else index = 0;
+
+ int ifix = modify->find_fix(&arg1[2]);
+ if (ifix < 0) error->all("Invalid fix ID in variable formula");
+ fix = modify->fix[ifix];
+ if (index == 0 && fix->vector_flag) {
+ if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+ error->all("Fix in variable not computed at compatible time");
+ nvec = fix->size_vector;
+ nstride = 1;
+ } else if (index && fix->array_flag) {
+ if (index > fix->size_array_cols)
+ error->all("Variable formula fix array is accessed out-of-range");
+ if (update->whichflag > 0 && update->ntimestep % fix->global_freq)
+ error->all("Fix in variable not computed at compatible time");
+ nvec = fix->size_array_rows;
+ nstride = fix->size_array_cols;
+ } else error->all("Mismatched fix in variable formula");
+
+ } else error->all("Invalid special function in variable formula");
+
+ double value = 0.0;
+ if (method == XMIN) value = BIG;
+ if (method == XMAX) value = -BIG;
+
+ if (compute) {
+ double *vec;
+ if (index) vec = &compute->array[0][index];
+ else vec = compute->vector;
+
+ int j = 0;
+ for (int i = 0; i < nvec; i++) {
+ if (method == SUM) value += vec[j];
+ else if (method == XMIN) value = MIN(value,vec[j]);
+ else if (method == XMAX) value = MAX(value,vec[j]);
+ else if (method == AVE) value += vec[j];
+ else if (method == TRAP) {
+ if (i > 0 && i < nvec-1) value += vec[j];
+ else value += 0.5*vec[j];
}
+ j += nstride;
}
-
- if (fix) {
- double one;
- for (int i = 0; i < nvec; i++) {
- if (index) one = fix->compute_array(i,index-1);
- else one = fix->compute_vector(i);
- if (method == SUM) value += one;
- else if (method == XMIN) value = MIN(value,one);
- else if (method == XMAX) value = MAX(value,one);
- else if (method == AVE) value += one;
- else if (method == TRAP) {
- if (i > 1 && i < nvec) value += one;
- else value += 0.5*one;
- }
- }
- }
-
- if (method == AVE) value /= nvec;
}
+
+ if (fix) {
+ double one;
+ for (int i = 0; i < nvec; i++) {
+ if (index) one = fix->compute_array(i,index-1);
+ else one = fix->compute_vector(i);
+ if (method == SUM) value += one;
+ else if (method == XMIN) value = MIN(value,one);
+ else if (method == XMAX) value = MAX(value,one);
+ else if (method == AVE) value += one;
+ else if (method == TRAP) {
+ if (i > 1 && i < nvec) value += one;
+ else value += 0.5*one;
+ }
+ }
+ }
+
+ if (method == AVE) value /= nvec;
delete [] arg1;
delete [] arg2;
@@ -2247,7 +2400,8 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
flag = 1 -> vector is a per-atom compute or fix quantity with nstride
id = positive global ID of atom, converted to local index
push result onto tree or arg stack
- customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+ customize by adding an atom vector:
+ mass,type,x,y,z,vx,vy,vz,fx,fy,fz
------------------------------------------------------------------------- */
void Variable::peratom2global(int flag, char *word,
@@ -2300,7 +2454,8 @@ void Variable::peratom2global(int flag, char *word,
/* ----------------------------------------------------------------------
check if word matches an atom vector
return 1 if yes, else 0
- customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+ customize by adding an atom vector:
+ mass,type,x,y,z,vx,vy,vz,fx,fy,fz
------------------------------------------------------------------------- */
int Variable::is_atom_vector(char *word)
@@ -2323,7 +2478,8 @@ int Variable::is_atom_vector(char *word)
process an atom vector in formula
push result onto tree
word = atom vector
- customize by adding an atom vector: mass,type,x,y,z,vx,vy,vz,fx,fy,fz
+ customize by adding an atom vector:
+ mass,type,x,y,z,vx,vy,vz,fx,fy,fz
------------------------------------------------------------------------- */
void Variable::atom_vector(char *word, Tree **tree,
@@ -2362,6 +2518,29 @@ void Variable::atom_vector(char *word, Tree **tree,
else if (strcmp(word,"fz") == 0) newtree->array = &atom->f[0][2];
}
+/* ----------------------------------------------------------------------
+ check if word matches a constant
+ return 1 if yes, else 0
+ customize by adding a constant: PI
+------------------------------------------------------------------------- */
+
+int Variable::is_constant(char *word)
+{
+ if (strcmp(word,"PI") == 0) return 1;
+ return 0;
+}
+
+/* ----------------------------------------------------------------------
+ process a constant in formula
+ customize by adding a constant: PI
+------------------------------------------------------------------------- */
+
+double Variable::constant(char *word)
+{
+ if (strcmp(word,"PI") == 0) return PI;
+ return 0.0;
+}
+
/* ----------------------------------------------------------------------
read a floating point value from a string
generate an error if not a legitimate floating point value
diff --git a/src/variable.h b/src/variable.h
index 4d075f0818..f5cec503c5 100644
--- a/src/variable.h
+++ b/src/variable.h
@@ -43,6 +43,7 @@ class Variable : protected Pointers {
int *which; // next available value for each variable
int *pad; // 1 = pad loop/uloop variables with 0s, 0 = no pad
char ***data; // str value of each variable's values
+ double PI;
class RanMars *randomequal; // random number generator for equal-style vars
class RanMars *randomatom; // random number generator for atom-style vars
@@ -74,6 +75,8 @@ class Variable : protected Pointers {
Tree **, Tree **, int &, double *, int &);
int is_atom_vector(char *);
void atom_vector(char *, Tree **, Tree **, int &);
+ int is_constant(char *);
+ double constant(char *);
double numeric(char *);
int inumeric(char *);
void print_tree(Tree *, int);