diff --git a/src/SHOCK/Install.sh b/src/SHOCK/Install.sh
index 84f40ddf7..1931770ee 100644
--- a/src/SHOCK/Install.sh
+++ b/src/SHOCK/Install.sh
@@ -1,19 +1,27 @@
# Install/unInstall package files in LAMMPS
if (test $1 = 1) then
+ cp fix_append_atoms.cpp ..
cp fix_msst.cpp ..
cp fix_nphug.cpp ..
+ cp fix_wall_piston.cpp ..
+ cp fix_append_atoms.h ..
cp fix_msst.h ..
cp fix_nphug.h ..
+ cp fix_wall_piston.h ..
elif (test $1 = 0) then
+ rm -f ../fix_append_atoms.cpp
rm -f ../fix_msst.cpp
rm -f ../fix_nphug.cpp
+ rm -f ../fix_wall_piston.cpp
+ rm -f ../fix_append_atoms.h
rm -f ../fix_msst.h
rm -f ../fix_nphug.h
+ rm -f ../fix_wall_piston.h
fi
diff --git a/src/SHOCK/fix_append_atoms.cpp b/src/SHOCK/fix_append_atoms.cpp
new file mode 100644
index 000000000..0bb38b48a
--- /dev/null
+++ b/src/SHOCK/fix_append_atoms.cpp
@@ -0,0 +1,483 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_append_atoms.h"
+#include "atom.h"
+#include "atom_vec.h"
+#include "comm.h"
+#include "modify.h"
+#include "domain.h"
+#include "lattice.h"
+#include "update.h"
+#include "random_mars.h"
+#include "error.h"
+#include "force.h"
+
+using namespace LAMMPS_NS;
+
+#define BIG 1.0e30
+#define EPSILON 1.0e-6
+
+/* ---------------------------------------------------------------------- */
+
+FixAppendAtoms::FixAppendAtoms(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ force_reneighbor = 1;
+ next_reneighbor = -1;
+ box_change = 1;
+ time_depend = 1;
+
+ if (narg < 4) error->all("Illegal fix append_atoms command");
+
+ scaleflag = 1;
+ spatflag=0;
+ xloflag = xhiflag = yloflag = yhiflag = zloflag = zhiflag = 0;
+
+ tempflag = 0;
+
+ ranflag = 0;
+ ranx = 0.0;
+ rany = 0.0;
+ ranz = 0.0;
+
+ randomx = NULL;
+ randomt = NULL;
+
+ int iarg = 0;
+ iarg = 3;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"xlo") == 0) {
+ error->all("Only zhi currently implemented for append_atoms");
+ xloflag = 1;
+ iarg++;
+ if (domain->boundary[0][0] != 3) error->all("Must shrink-wrap with minimum the append boundary");
+ } else if (strcmp(arg[iarg],"xhi") == 0) {
+ error->all("Only zhi currently implemented for append_atom");
+ xhiflag = 1;
+ iarg++;
+ if (domain->boundary[0][1] != 3) error->all("Must shrink-wrap with minimum th append boundary");
+ } else if (strcmp(arg[iarg],"ylo") == 0) {
+ error->all("Only zhi currently implemented for append_atom");
+ yloflag = 1;
+ iarg++;
+ if (domain->boundary[1][0] != 3) error->all("Must shrink-wrap with minimum th append boundary");
+ } else if (strcmp(arg[iarg],"yhi") == 0) {
+ error->all("Only zhi currently implemented for append_atom");
+ yhiflag = 1;
+ iarg++;
+ if (domain->boundary[1][1] != 3) error->all("Must shrink-wrap with minimum th append boundary");
+ } else if (strcmp(arg[iarg],"zlo") == 0) {
+ error->all("Only zhi currently implemented for append_atom");
+ zloflag = 1;
+ iarg++;
+ if (domain->boundary[2][0] != 3) error->all("Must shrink-wrap with minimum th append boundary");
+ } else if (strcmp(arg[iarg],"zhi") == 0) {
+ zhiflag = 1;
+ iarg++;
+ if (domain->boundary[2][1] != 3) error->all("Must shrink-wrap with minimum th append boundary");
+ } else if (strcmp(arg[iarg],"freq") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix append_atoms command");
+ freq = atof(arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"spatial") == 0) {
+ if (iarg+3 > narg) error->all("Illegal fix append_atoms command");
+ if (strcmp(arg[iarg+1],"f_") == 0) error->all("Bad fix ID in fix append_atoms command");
+ spatflag = 1;
+ int n = strlen(arg[iarg+1]);
+ spatlead = atof(arg[iarg+2]);
+ char *suffix = new char[n];
+ strcpy(suffix,&arg[iarg+1][2]);
+ n = strlen(suffix) + 1;
+ spatialid = new char[n];
+ strcpy(spatialid,suffix);
+ delete [] suffix;
+ iarg += 3;
+ // NEED TO CHECK TO MAKE SURE FIX IS AN AVE/SPATIAL
+ } else if (strcmp(arg[iarg],"size") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix append_atoms command");
+ size = atof(arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"units") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix append_atoms command");
+ if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+ else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+ else error->all("Illegal fix append_atoms command");
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"random") == 0) {
+ if (iarg+5 > narg) error->all("Illegal fix append_atoms command");
+ ranflag = 1;
+ ranx = atof(arg[iarg+1]);
+ rany = atof(arg[iarg+2]);
+ ranz = atof(arg[iarg+3]);
+ xseed = atoi(arg[iarg+4]);
+ if (xseed <= 0) error->all("Illegal fix append_atoms command");
+ randomx = new RanMars(lmp,xseed + comm->me);
+ iarg += 5;
+ } else if (strcmp(arg[iarg],"temp") == 0) {
+ if (iarg+5 > narg) error->all("Illegal fix append_atoms command");
+ tempflag = 1;
+ t_target = atof(arg[iarg+1]);
+ t_period = atof(arg[iarg+2]);
+ tseed = atoi(arg[iarg+3]);
+ t_extent = atof(arg[iarg+4]);
+ if (t_target <= 0) error->all("Illegal fix append_atoms command");
+ if (t_period <= 0) error->all("Illegal fix append_atoms command");
+ if (t_extent <= 0) error->all("Illegal fix append_atoms command");
+ if (tseed <= 0) error->all("Illegal fix append_atoms command");
+ randomt = new RanMars(lmp,tseed + comm->me);
+ gfactor1 = new double[atom->ntypes+1];
+ gfactor2 = new double[atom->ntypes+1];
+ iarg += 5;
+ } else error->all("Illegal fix append_atoms command");
+ }
+
+ if ((xloflag || xhiflag) && domain->xperiodic)
+ error->all("Cannot use append_atoms in periodic dimension");
+ if ((yloflag || yhiflag) && domain->yperiodic)
+ error->all("Cannot use append_atoms in periodic dimension");
+ if ((zloflag || zhiflag) && domain->zperiodic)
+ error->all("Cannot use append_atoms in periodic dimension");
+
+ if (domain->triclinic == 1) error->all("Cannot append atoms to a triclinic box");
+
+ // setup scaling
+
+ if (scaleflag && domain->lattice == NULL)
+ error->all("Use of fix append_atoms 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;
+
+ if (xloflag || xhiflag) size *= xscale;
+ if (yloflag || yhiflag) size *= yscale;
+ if (zloflag || zhiflag) size *= zscale;
+
+ if (ranflag) {
+ ranx *= xscale;
+ rany *= yscale;
+ ranz *= zscale;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAppendAtoms::~FixAppendAtoms()
+{
+ if (ranflag) delete randomx;
+ if (tempflag) {
+ delete randomt;
+ delete [] gfactor1;
+ delete [] gfactor2;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAppendAtoms::setmask()
+{
+ int mask = 0;
+ mask |= PRE_EXCHANGE;
+ mask |= INITIAL_INTEGRATE;
+ mask |= POST_FORCE;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAppendAtoms::initial_integrate(int vflag)
+{
+ if (update->ntimestep % freq == 0) {
+ next_reneighbor = update->ntimestep;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAppendAtoms::setup(int vflag)
+{
+ /*** CALL TO CREATE GROUP? SEE POST_FORCE ***/
+ post_force(vflag);
+}
+
+
+/* ---------------------------------------------------------------------- */
+
+int FixAppendAtoms::get_spatial()
+{
+ if (update->ntimestep % freq == 0) {
+ int ifix = modify->find_fix(spatialid);
+ if (ifix < 0)
+ error->all("Fix ID for fix ave/spatial does not exist");
+ Fix *fix = modify->fix[ifix];
+
+ int failed = 0;
+ int count = 0;
+ while (failed < 2) {
+ double tmp = fix->compute_vector(2*count);
+ if (tmp == 0.0) failed++;
+ else failed = 0;
+ count++;
+ }
+ double pos[count-2];
+ double val[count-2];
+ for (int loop=0; loop < count-2; loop++) {
+ pos[loop] = fix->compute_vector(2*loop);
+ val[loop] = fix->compute_vector(2*loop+1);
+ }
+
+// Always ignor the first and last
+// SHOULD HAVE A MEMORY OF MIN AND MAX ENERGY
+// CAPTURE BINSIZE FROM SPATIAL
+
+ double binsize = 2.0;
+ double min_energy=0.0;
+ double max_energy=0.0;
+ int header = size / binsize;
+ advance = 0;
+
+ for (int loop=1; loop <= header; loop++) {
+ max_energy += val[loop];
+ }
+ for (int loop=count-2-2*header; loop <=count-3-header; loop++) {
+ min_energy += val[loop];
+ }
+ max_energy /= header;
+ min_energy /= header;
+
+ double shockfront_min = 0.0;
+ double shockfront_max = 0.0;
+ double shockfront_loc = 0.0;
+ int front_found1 = 0;
+ for (int loop=count-3-header; loop > header; loop--) {
+ if (front_found1 == 1) continue;
+ if (val[loop] > min_energy + 0.1*(max_energy - min_energy)) {
+ shockfront_max = pos[loop];
+ front_found1=1;
+ }
+ }
+ int front_found2 = 0;
+ for (int loop=header+1; loop <=count-3-header; loop++) {
+ if (val[loop] > min_energy + 0.6*(max_energy - min_energy)) {
+ shockfront_min = pos[loop];
+ front_found2=1;
+ }
+ }
+ if (front_found1 + front_found2 == 0) shockfront_loc = 0.0;
+ else if (front_found1 + front_found2 == 1) shockfront_loc = shockfront_max + shockfront_min;
+ else if (front_found1 == 1 && front_found2 == 1 && shockfront_max-shockfront_min > spatlead/2.0) shockfront_loc = shockfront_max;
+ else shockfront_loc = (shockfront_max + shockfront_min) / 2.0;
+ if (comm->me == 0) printf("SHOCK: %g %g %g %g %g\n", shockfront_loc, shockfront_min, shockfront_max, domain->boxlo[2], domain->boxhi[2]);
+
+ if (domain->boxhi[2] - shockfront_loc < spatlead) advance = 1;
+ }
+
+ advance_sum = 0;
+ MPI_Allreduce(&advance,&advance_sum,1,MPI_INT,MPI_SUM,world);
+
+ if (advance_sum > 0) return 1;
+ else return 0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAppendAtoms::post_force(int vflag)
+{
+ double **f = atom->f;
+ double **v = atom->v;
+ double **x = atom->x;
+ int *type = atom->type;
+ int nlocal = atom->nlocal;
+
+ double gamma1,gamma2;
+ double tsqrt = sqrt(t_target);
+
+ if (atom->mass) {
+ if (tempflag) {
+ for (int i = 1; i <= atom->ntypes; i++) {
+ gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
+ gfactor2[i] = sqrt(atom->mass[i]) *
+ sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
+ force->ftm2v;
+ }
+ }
+ for (int i = 0; i < nlocal; i++) {
+ // SET TEMP AHEAD OF SHOCK
+ if (tempflag && x[i][2] >= domain->boxhi[2] - t_extent ) {
+ gamma1 = gfactor1[type[i]];
+ gamma2 = gfactor2[type[i]] * tsqrt;
+ f[i][0] += gamma1*v[i][0] + gamma2*(randomt->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + gamma2*(randomt->uniform()-0.5);
+ f[i][2] += gamma1*v[i][2] + gamma2*(randomt->uniform()-0.5);
+ }
+ // FREEZE ATOMS AT BOUNDARY
+ if (x[i][2] >= domain->boxhi[2] - size) {
+ f[i][0] = 0.0;
+ f[i][1] = 0.0;
+ f[i][2] = 0.0;
+ v[i][0] = 0.0;
+ v[i][1] = 0.0;
+ v[i][2] = 0.0;
+ }
+ }
+ } else {
+ double *rmass = atom->rmass;
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ for (int i = 0; i < nlocal; i++) {
+ // SET TEMP AHEAD OF SHOCK
+ if (tempflag && x[i][2] >= domain->boxhi[2] - t_extent ) {
+ gamma1 = -rmass[i] / t_period / ftm2v;
+ gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ gamma2 *= tsqrt;
+ f[i][0] += gamma1*v[i][0] + gamma2*(randomt->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + gamma2*(randomt->uniform()-0.5);
+ f[i][2] += gamma1*v[i][2] + gamma2*(randomt->uniform()-0.5);
+ }
+ // FREEZE ATOMS AT BOUNDARY
+ if (x[i][2] >= domain->boxhi[2] - size) {
+ f[i][0] = 0.0;
+ f[i][1] = 0.0;
+ f[i][2] = 0.0;
+ v[i][0] = 0.0;
+ v[i][1] = 0.0;
+ v[i][2] = 0.0;
+ }
+ }
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAppendAtoms::pre_exchange()
+{
+ int ntimestep = update->ntimestep;
+ int addnode = 0;
+
+ if (ntimestep % freq == 0) {
+ if (spatflag==1) if (get_spatial()==0) return;
+ if (comm->myloc[2] == comm->procgrid[2]-1) {
+ if (domain->lattice) {
+ nbasis = domain->lattice->nbasis;
+ basistype = new int[nbasis];
+ for (int i = 0; i < nbasis; i++) basistype[i] = 1;
+ } else error->all("must define lattice to append_atoms");
+
+ double bboxlo[3],bboxhi[3];
+
+ bboxlo[0] = domain->sublo[0]; bboxhi[0] = domain->subhi[0];
+ bboxlo[1] = domain->sublo[1]; bboxhi[1] = domain->subhi[1];
+ bboxlo[2] = domain->subhi[2]; bboxhi[2] = domain->subhi[2]+size;
+
+ double xmin,ymin,zmin,xmax,ymax,zmax;
+ xmin = ymin = zmin = BIG;
+ xmax = ymax = zmax = -BIG;
+
+ domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxlo[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxlo[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxlo[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxlo[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxlo[0],bboxlo[1],bboxhi[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxhi[0],bboxlo[1],bboxhi[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxlo[0],bboxhi[1],bboxhi[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+ domain->lattice->bbox(1,bboxhi[0],bboxhi[1],bboxhi[2],
+ xmin,ymin,zmin,xmax,ymax,zmax);
+
+ int ilo,ihi,jlo,jhi,klo,khi;
+ ilo = static_cast<int> (xmin);
+ jlo = static_cast<int> (ymin);
+ klo = static_cast<int> (zmin);
+ ihi = static_cast<int> (xmax);
+ jhi = static_cast<int> (ymax);
+ khi = static_cast<int> (zmax);
+
+ if (xmin < 0.0) ilo--;
+ if (ymin < 0.0) jlo--;
+ if (zmin < 0.0) klo--;
+
+ double **basis = domain->lattice->basis;
+ double x[3];
+ double *sublo = domain->sublo;
+ double *subhi = domain->subhi;
+ double *mass = atom->mass;
+
+ int i,j,k,m;
+ for (k = klo; k <= khi; k++) {
+ for (j = jlo; j <= jhi; j++) {
+ for (i = ilo; i <= ihi; i++) {
+ for (m = 0; m < nbasis; m++) {
+ x[0] = i + basis[m][0];
+ x[1] = j + basis[m][1];
+ x[2] = k + basis[m][2];
+
+ int flag = 0;
+ // convert from lattice coords to box coords
+ domain->lattice->lattice2box(x[0],x[1],x[2]);
+
+ if (x[0] >= sublo[0] && x[0] < subhi[0] &&
+ x[1] >= sublo[1] && x[1] < subhi[1] &&
+ x[2] >= subhi[2] && x[2] < subhi[2]+size) flag = 1;
+ else if (domain->dimension == 2 && x[1] >= domain->boxhi[1] &&
+ comm->myloc[1] == comm->procgrid[1]-1 &&
+ x[0] >= sublo[0] && x[0] < subhi[0]) flag = 1;
+
+ if (flag) {
+ if (ranflag) {
+ x[0] += ranx * 2.0*(randomx->uniform()-0.5);
+ x[1] += rany * 2.0*(randomx->uniform()-0.5);
+ x[2] += ranz * 2.0*(randomx->uniform()-0.5);
+ }
+ addnode++;
+ atom->avec->create_atom(basistype[m],x);
+ }
+ }
+ }
+ }
+ }
+ }
+ int addtotal = 0;
+ MPI_Barrier(world);
+ MPI_Allreduce(&addnode,&addtotal,1,MPI_INT,MPI_SUM,world);
+
+ if (addtotal) {
+ domain->reset_box();
+ if (atom->tag_enable) {
+ atom->tag_extend();
+ atom->natoms += addtotal;
+ if (atom->map_style) {
+ atom->nghost = 0;
+ atom->map_init();
+ atom->map_set();
+ }
+ }
+ }
+ }
+}
diff --git a/src/SHOCK/fix_append_atoms.h b/src/SHOCK/fix_append_atoms.h
new file mode 100644
index 000000000..3d101b921
--- /dev/null
+++ b/src/SHOCK/fix_append_atoms.h
@@ -0,0 +1,56 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(append_atoms,FixAppendAtoms)
+
+#else
+
+#ifndef FIX_APPEND_ATOMS_H
+#define FIX_APPEND_ATOMS_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixAppendAtoms : public Fix {
+ public:
+ FixAppendAtoms(class LAMMPS *, int, char **);
+ ~FixAppendAtoms();
+ int setmask();
+ void setup(int);
+ void pre_exchange();
+ void initial_integrate(int);
+ void post_force(int);
+
+ private:
+ int get_spatial();
+ int spatflag, xloflag, xhiflag, yloflag, yhiflag, zloflag, zhiflag;
+ int ranflag, tempflag, xseed, tseed;
+ double ranx, rany, ranz, t_target, t_period, t_extent;
+ class RanMars *randomx;
+ class RanMars *randomt;
+ int scaleflag, freq;
+ int *basistype, nbasis;
+ int advance, advance_sum;
+ double size, spatlead;
+ char *spatialid;
+ double tfactor;
+ double *gfactor1,*gfactor2;
+};
+
+}
+
+#endif
+#endif
diff --git a/src/SHOCK/fix_wall_piston.cpp b/src/SHOCK/fix_wall_piston.cpp
new file mode 100644
index 000000000..9dd6325bb
--- /dev/null
+++ b/src/SHOCK/fix_wall_piston.cpp
@@ -0,0 +1,325 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_wall_piston.h"
+#include "atom.h"
+#include "modify.h"
+#include "domain.h"
+#include "lattice.h"
+#include "update.h"
+#include "error.h"
+#include "random_mars.h"
+#include "force.h"
+#include "comm.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixWallPiston::FixWallPiston(LAMMPS *lmp, int narg, char **arg) :
+ Fix(lmp, narg, arg)
+{
+ force_reneighbor = 1;
+ next_reneighbor = -1;
+ box_change = 1;
+ time_depend = 1;
+
+ if (narg < 4) error->all("Illegal fix wall/piston command");
+
+ randomt = NULL;
+ tempflag = 0;
+ scaleflag = 1;
+ roughflag = 0;
+ rampflag = 0;
+ rampNL1flag = 0;
+ rampNL2flag = 0;
+ rampNL3flag = 0;
+ rampNL4flag = 0;
+ rampNL5flag = 0;
+ x0 = y0 = z0 = vx = vy = vz = 0.0;
+ xloflag = xhiflag = yloflag = yhiflag = zloflag = zhiflag = 0;
+ int iarg = 0;
+ iarg = 3;
+ while (iarg < narg) {
+ if (strcmp(arg[iarg],"xlo") == 0) { error->all("Fix wall/piston command only available at zlo");
+ } else if (strcmp(arg[iarg],"ylo") == 0) { error->all("Fix wall/piston command only available at zlo");
+ } else if (strcmp(arg[iarg],"zlo") == 0) {
+ zloflag = 1;
+ iarg++;
+ if (domain->boundary[2][0] != 2) error->all("Must shrink-wrap piston boundary");
+ } else if (strcmp(arg[iarg],"xhi") == 0) { error->all("Fix wall/piston command only available at zlo");
+ } else if (strcmp(arg[iarg],"yhi") == 0) { error->all("Fix wall/piston command only available at zlo");
+ } else if (strcmp(arg[iarg],"zhi") == 0) { error->all("Fix wall/piston command only available at zlo");
+ } else if (strcmp(arg[iarg],"vel") == 0) {
+ if (iarg+4 > narg) error->all("Illegal fix wall/piston command");
+ vx = atof(arg[iarg+1]);
+ vy = atof(arg[iarg+2]);
+ vz = atof(arg[iarg+3]);
+ iarg += 4;
+ } else if (strcmp(arg[iarg],"pos") == 0) {
+ if (iarg+4 > narg) error->all("Illegal fix wall/piston command");
+ x0 = atof(arg[iarg+1]);
+ y0 = atof(arg[iarg+2]);
+ z0 = atof(arg[iarg+3]);
+ iarg += 4;
+ } else if (strcmp(arg[iarg],"temp") == 0) {
+ if (iarg+5 > narg) error->all("Illegal fix wall/pistons command");
+ tempflag = 1;
+ t_target = atof(arg[iarg+1]);
+ t_period = atof(arg[iarg+2]);
+ tseed = atoi(arg[iarg+3]);
+ t_extent = atof(arg[iarg+4]);
+ if (t_target <= 0) error->all("Illegal fix wall/piston command");
+ if (t_period <= 0) error->all("Illegal fix wall/piston command");
+ if (t_extent <= 0) error->all("Illegal fix wall/piston command");
+ if (tseed <= 0) error->all("Illegal fix wall/pistons command");
+ randomt = new RanMars(lmp,tseed + comm->me);
+ gfactor1 = new double[atom->ntypes+1];
+ gfactor2 = new double[atom->ntypes+1];
+ iarg += 5;
+ } else if (strcmp(arg[iarg],"rough") == 0) {
+ roughflag = 1;
+ roughdist = atof(arg[iarg+1]);
+ iarg += 2;
+ } else if (strcmp(arg[iarg],"ramp") == 0) {
+ rampflag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"rampNL1") == 0) {
+ rampNL1flag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"rampNL2") == 0) {
+ rampNL2flag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"rampNL3") == 0) {
+ rampNL3flag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"rampNL4") == 0) {
+ rampNL4flag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"rampNL5") == 0) {
+ rampNL5flag = 1;
+ iarg++;
+ } else if (strcmp(arg[iarg],"units") == 0) {
+ if (iarg+2 > narg) error->all("Illegal fix wall/piston command");
+ if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+ else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+ else error->all("Illegal fix wall/piston command");
+ iarg += 2;
+ } else error->all("Illegal fix wall/piston command");
+ }
+
+ if (vx < 0.0 || vy < 0.0 || vz < 0.0) error->all("Illegal fix wall/piston velocity");
+ 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");
+
+ // setup scaling
+
+ if (scaleflag && domain->lattice == NULL)
+ error->all("Use of fix wall/piston 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;
+
+ vx *= xscale;
+ vy *= yscale;
+ vz *= zscale;
+ x0 *= xscale;
+ y0 *= yscale;
+ z0 *= zscale;
+ roughdist *= zscale;
+
+ if (rampflag || rampNL1flag || rampNL2flag || rampNL3flag || rampNL4flag || rampNL5flag) {
+ maxvx = vx;
+ maxvy = vy;
+ maxvz = vz;
+ }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixWallPiston::setmask()
+{
+ int mask = 0;
+ mask |= POST_INTEGRATE;
+ mask |= POST_INTEGRATE_RESPA;
+ mask |= INITIAL_INTEGRATE;
+ return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallPiston::initial_integrate(int vflag)
+{
+ next_reneighbor = update->ntimestep;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixWallPiston::post_integrate()
+{
+ double xlo, xhi, ylo, yhi, zlo, zhi;
+
+ double **x = atom->x;
+ double **v = atom->v;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
+ double t = (update->ntimestep - update->beginstep) * update->dt;
+ double tott = (update->endstep - update->beginstep) * update->dt;
+ double tt = t * t;
+ double ttt = tt * t;
+ double tttt = tt * tt;
+ double t0p5 = sqrt(t/tott);
+ double t1p5 = t0p5*t0p5*t0p5;
+ double t2p5 = t1p5*t0p5*t0p5;
+
+ if (rampflag) {
+ paccelx = maxvx / tott;
+ paccely = maxvy / tott;
+ paccelz = maxvz / tott;
+
+ if (zloflag) { zlo = z0 + 0.5 * paccelz * tt; vz = paccelz * t; }
+ }
+ else if (rampNL1flag) {
+ paccelz = maxvz / tott;
+ angfreq = 6.283185 / (0.5 * tott);
+
+ if (zloflag) {
+ zlo = z0 + paccelz * (0.5*tt + 1.0/(angfreq*angfreq) - 1.0/(angfreq*angfreq)*cos(angfreq*t));
+ vz = paccelz * (t + 1.0/angfreq*sin(angfreq*t));
+ }
+ else { error->all("NL ramp in wall/piston only implemented in zlo for now"); }
+ }
+ else if (rampNL2flag) {
+ paccelz = maxvz / tott;
+ angfreq = 18.84956 / tott;
+
+ if (zloflag) {
+ zlo = z0 + paccelz * (0.5*tt + 4.0/(3.0*angfreq*angfreq)*(1.0-cos(angfreq*t)) + 1.0/(6.0*angfreq*angfreq)*(1.0-cos(2.0*angfreq*t)));
+ vz = paccelz * (t + 4.0/(3.0*angfreq)*sin(angfreq*t) + 1.0/(3.0*angfreq)*sin(2.0*angfreq*t));
+ }
+ else { error->all("NL ramp in wall/piston only implemented in zlo for now"); }
+ }
+ else if (rampNL3flag) {
+ paccelz = maxvz / tott;
+
+ if (zloflag) {
+ zlo = z0 + paccelz*tott*tott/2.5 * (t2p5 );
+ vz = paccelz * tott * (t1p5 );
+ }
+ else { error->all("NL ramp in wall/piston only implemented in zlo for now"); }
+ }
+ else if (rampNL4flag) {
+ paccelz = maxvz / tott;
+
+ if (zloflag) {
+ zlo = z0 + paccelz/tott/3.0 * (ttt);
+ vz = paccelz / tott * (tt);
+ }
+ else { error->all("NL ramp in wall/piston only implemented in zlo for now"); }
+ }
+ else if (rampNL5flag) {
+ paccelz = maxvz / tott;
+
+ if (zloflag) {
+ zlo = z0 + paccelz/tott/tott/4.0 * (tttt);
+ vz = paccelz / tott / tott * (ttt);
+ }
+ else { error->all("NL ramp in wall/piston only implemented in zlo for now"); }
+ }
+ else {
+ if (zloflag) { zlo = z0 + vz * t; }
+ }
+
+ if (update->ntimestep % 1000 == 0)
+ if (comm->me == 0) {
+ if (screen)
+ fprintf(screen,"SHOCK: step %d t %g zpos %g vz %g az %g zlo %g\n",
+ update->ntimestep, t, zlo, vz, paccelz, domain->boxlo[2]);
+ if (logfile)
+ fprintf(logfile,"SHOCK: step %d t %g zpos %g vz %g az %g zlo %g\n",
+ update->ntimestep, t, zlo, vz, paccelz, domain->boxlo[2]);
+ }
+
+ // VIRIAL PRESSURE CONTRIBUTION?
+
+ for (int i = 0; i < nlocal; i++) {
+ if (mask[i] & groupbit) {
+ roughoff = 0.0;
+ if (roughflag) {
+ roughoff += roughdist*fabs((x[i][0] - domain->boxlo[0])/(domain->boxhi[0]-domain->boxlo[0])-0.5);
+ roughoff += roughdist*fabs((x[i][1] - domain->boxlo[1])/(domain->boxhi[1]-domain->boxlo[1])-0.5);
+ }
+ if (zloflag && x[i][2] < zlo - roughoff) {
+ x[i][2] = 2.0 * (zlo - roughoff) - x[i][2];
+ v[i][2] = 2.0 * vz - v[i][2];
+ }
+ }
+ }
+ double **f = atom->f;
+ int *type = atom->type;
+
+ double gamma1,gamma2;
+ double tsqrt = sqrt(t_target);
+
+ if (atom->mass) {
+ if (tempflag) {
+ for (int i = 1; i <= atom->ntypes; i++) {
+ gfactor1[i] = -atom->mass[i] / t_period / force->ftm2v;
+ gfactor2[i] = sqrt(atom->mass[i]) *
+ sqrt(24.0*force->boltz/t_period/update->dt/force->mvv2e) /
+ force->ftm2v;
+ }
+ }
+ for (int i = 0; i < nlocal; i++) {
+ // SET TEMP AHEAD OF PISTON
+ if (tempflag && x[i][2] <= domain->boxlo[2] + t_extent ) {
+ gamma1 = gfactor1[type[i]];
+ gamma2 = gfactor2[type[i]] * tsqrt;
+ f[i][0] += gamma1*v[i][0] + gamma2*(randomt->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + gamma2*(randomt->uniform()-0.5);
+ f[i][2] += gamma1*(v[i][2]-vz) + gamma2*(randomt->uniform()-0.5);
+ }
+ }
+ } else {
+ double *rmass = atom->rmass;
+ double boltz = force->boltz;
+ double dt = update->dt;
+ double mvv2e = force->mvv2e;
+ double ftm2v = force->ftm2v;
+
+ for (int i = 0; i < nlocal; i++) {
+ // SET TEMP AHEAD OF PISTON
+ if (tempflag && x[i][2] <= domain->boxlo[2] + t_extent ) {
+ gamma1 = -rmass[i] / t_period / ftm2v;
+ gamma2 = sqrt(rmass[i]) * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
+ gamma2 *= tsqrt;
+ f[i][0] += gamma1*v[i][0] + gamma2*(randomt->uniform()-0.5);
+ f[i][1] += gamma1*v[i][1] + gamma2*(randomt->uniform()-0.5);
+ f[i][2] += gamma1*v[i][2] + gamma2*(randomt->uniform()-0.5);
+ }
+ }
+ }
+
+}
diff --git a/src/SHOCK/fix_wall_piston.h b/src/SHOCK/fix_wall_piston.h
new file mode 100644
index 000000000..46d4cdee6
--- /dev/null
+++ b/src/SHOCK/fix_wall_piston.h
@@ -0,0 +1,46 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+#ifdef FIX_CLASS
+
+FixStyle(wall/piston,FixWallPiston)
+
+#else
+
+#ifndef LMP_FIX_WALL_PISTON_H
+#define LMP_FIX_WALL_PISTON_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixWallPiston : public Fix {
+ public:
+ FixWallPiston(class LAMMPS *, int, char **);
+ int setmask();
+ void post_integrate();
+ void initial_integrate(int);
+
+ private:
+ int xloflag,xhiflag,yloflag,yhiflag,zloflag,zhiflag;
+ int scaleflag, roughflag, rampflag, rampNL1flag, rampNL2flag, rampNL3flag, rampNL4flag, rampNL5flag;
+ double roughdist,roughoff,x0,y0,z0,vx,vy,vz,maxvx,maxvy,maxvz,paccelx,paccely,paccelz, angfreq;
+ int tempflag, tseed;
+ double t_target, t_period, t_extent;
+ class RanMars *randomt;
+ double *gfactor1,*gfactor2;
+};
+
+}
+
+#endif
+#endif