ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added Matt's piston and append code

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6783 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
2011-08-24 16:53:00 +00:00
parent 029a4d2df2
commit 0cb40558c5
5 changed files with 918 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/SHOCK/Install.sh
View File

@ -2,18 +2,26 @@
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

483
src/SHOCK/fix_append_atoms.cpp Normal file
View File

@ -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();
}
}
}
}
}

56
src/SHOCK/fix_append_atoms.h Normal file
View File

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

325
src/SHOCK/fix_wall_piston.cpp Normal file
View File

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

46
src/SHOCK/fix_wall_piston.h Normal file
View File

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