From aa70f8cc6b52c94b15e0288e5f032ead27594781 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Mon, 13 Jan 2020 14:26:22 -0700 Subject: [PATCH] small edits for reformatting --- src/USER-MISC/fix_wall_reflect_stochastic.cpp | 540 +++++++++--------- src/USER-MISC/fix_wall_reflect_stochastic.h | 165 +++--- src/fix_wall_reflect.cpp | 493 ++++++++-------- src/fix_wall_reflect.h | 8 +- src/random_mars.cpp | 46 +- 5 files changed, 623 insertions(+), 629 deletions(-) diff --git a/src/USER-MISC/fix_wall_reflect_stochastic.cpp b/src/USER-MISC/fix_wall_reflect_stochastic.cpp index dda2a9296f..346db84464 100644 --- a/src/USER-MISC/fix_wall_reflect_stochastic.cpp +++ b/src/USER-MISC/fix_wall_reflect_stochastic.cpp @@ -1,266 +1,274 @@ -/* ---------------------------------------------------------------------- - 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 "fix_wall_reflect_stochastic.h" -#include -#include -#include "atom.h" -#include "comm.h" -#include "update.h" -#include "modify.h" -#include "domain.h" -#include "lattice.h" -#include "input.h" -#include "variable.h" -#include "error.h" -#include "force.h" -#include "random_mars.h" - - -using namespace LAMMPS_NS; -using namespace FixConst; - -enum{NONE=0,DIFFUSIVE=1,MAXWELL=2,CERCIGNANILAMPIS=3}; - - -/* ---------------------------------------------------------------------- */ - -FixWallReflectStochastic::FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : - FixWallReflect(lmp, narg, arg), random(NULL) -{ - int arginc,dir; - - if (narg < 8) error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); - - if (domain->triclinic != 0) error->all(FLERR, "This fix wall/reflect/stochastic is not incompatible with triclinic simulation box"); - - dynamic_group_allow = 1; - - // parse args - - nwall = 0; - int scaleflag = 1; - reflectionstyle = NONE; - - - if (strcmp(arg[3],"diffusive") == 0) { - reflectionstyle = DIFFUSIVE; - arginc = 6; - } else if (strcmp(arg[3],"maxwell") == 0) { - reflectionstyle = MAXWELL; - arginc = 7; - } else if (strcmp(arg[3],"cercignanilampis") == 0) { - reflectionstyle = CERCIGNANILAMPIS; - arginc = 9; - } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); - - - seedfix = force->inumeric(FLERR,arg[4]); - if (seedfix <=0) error->all(FLERR,"Random seed must be a postive number"); - random = new RanMars(lmp,seedfix + comm->me); - - int iarg = 5; - while (iarg < narg) { - 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(FLERR,"Illegal fix wall/reflect/stochastic 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 < 6); m++) - if (newwall == wallwhich[m]) - error->all(FLERR,"Wall defined twice in fix wall/reflect/stochastic 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 { - wallstyle[nwall] = CONSTANT; - coord0[nwall] = force->numeric(FLERR,arg[iarg+1]); - } - - - walltemp[nwall]= force->numeric(FLERR,arg[iarg+2]); - - for (dir = 0; dir < 3; dir++) { - wallvel[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+3]); - int dim = wallwhich[nwall] / 2; - if ((wallvel[nwall][dir] !=0) & (dir == dim)) error->all(FLERR,"The wall velocity must be tangential"); - - if (reflectionstyle == CERCIGNANILAMPIS) { - wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+6]); - } else if (reflectionstyle == MAXWELL) { - wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+6]);// One accommodation coefficient for all directions - } - } - - nwall++; - iarg += arginc; - - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect/stochastic command"); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); - iarg += 2; - } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); - } - - // error check - - if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); - - for (int m = 0; m < nwall; m++) { - if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) - error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension"); - if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) - error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension"); - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) - error->all(FLERR,"Cannot use fix wall/reflect/stochastic in periodic dimension"); - } - - for (int m = 0; m < nwall; m++) - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) - error->all(FLERR, - "Cannot use fix wall/reflect/stochastic zlo/zhi for a 2d simulation"); - - // scale factors for CONSTANT walls, VARIABLE wall is not allowed - - int flag = 0; - for (int m = 0; m < nwall; m++) - if (wallstyle[m] != EDGE) flag = 1; - - if (flag) { - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; - - for (int m = 0; m < nwall; m++) { - if (wallstyle[m] != CONSTANT) continue; - if (wallwhich[m] < YLO) coord0[m] *= xscale; - else if (wallwhich[m] < ZLO) coord0[m] *= yscale; - else coord0[m] *= zscale; - } - } - - } - - -FixWallReflectStochastic::~FixWallReflectStochastic() -{ - -delete random ; - -} - -/* ---------------------------------------------------------------------- */ - -void FixWallReflectStochastic::wall_particle(int m, int which, double coord) -{ - int i, dir, dim, side, sign; - double vsave,factor,timecol,difftest,theta; - - double *rmass; - double *mass = atom->mass; - - // coord = current position of wall - // evaluate variable if necessary, wrap with clear/add - - double **x = atom->x; - double **v = atom->v; - int *mask = atom->mask; - int *type = atom->type; - int nlocal = atom->nlocal; - - dim = which / 2; - side = which % 2; - - for (i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - - sign = 0; - - if ((side == 0) & (x[i][dim] < coord)){ - sign = 1; - } else if ((side == 1) & (x[i][dim] > coord)){ - sign = -1; - } - - - if (sign != 0) { - if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e);// theta = kT/m - else theta = force->boltz*walltemp[m]/(mass[type[i]]*force->mvv2e); - factor = sqrt(theta); - - timecol= (x[i][dim] - coord)/v[i][dim]; // time travelling after collision - - if (reflectionstyle == MAXWELL) {difftest = random->uniform(); }// for Maxwell model only - - for (dir = 0; dir < 3; dir++) { - - x[i][dir] -= v[i][dir]*timecol; // pushing back atoms to the wall position and assign new reflected velocity - - // Diffusive reflection - if (reflectionstyle == DIFFUSIVE){ - - if (dir != dim) { - v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); - } else { v[i][dir] = sign*random->rayleigh(factor) ; } - - } - - // Maxwell reflection - if (reflectionstyle == MAXWELL){ - - if (difftest < wallaccom[m][dir]) { - if (dir != dim) { - v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); - } else { v[i][dir] = sign*random->rayleigh(factor) ; } - } else { - if (dir == dim) {v[i][dir] = -v[i][dir]; - } - } - } - - - // Cercignani Lampis reflection - if (reflectionstyle == CERCIGNANILAMPIS){ - if (dir != dim) { - v[i][dir] = wallvel[m][dir] + random->gaussian((1-wallaccom[m][dir])*(v[i][dir]-wallvel[m][dir]),sqrt((2-wallaccom[m][dir])*wallaccom[m][dir]*theta)) ; - } else { v[i][dir] = random->besselexp(theta, wallaccom[m][dir], v[i][dir]) ; } - - } - - // update new position due to the new velocity - if (dir != dim) { - x[i][dir] += v[i][dir]*timecol; } - else {x[i][dir] = coord + v[i][dir]*timecol;} - }// for loop - }// if sign - }// if mask - } -} \ No newline at end of file +/* ---------------------------------------------------------------------- + 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 "fix_wall_reflect_stochastic.h" +#include +#include +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "lattice.h" +#include "input.h" +#include "variable.h" +#include "random_mars.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +enum{NONE=0,DIFFUSIVE=1,MAXWELL=2,CERCIGNANILAMPIS=3}; + +/* ---------------------------------------------------------------------- */ + +FixWallReflectStochastic:: +FixWallReflectStochastic(LAMMPS *lmp, int narg, char **arg) : + FixWallReflect(lmp, narg, arg), random(NULL) +{ + if (narg < 8) error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + + if (domain->triclinic != 0) + error->all(FLERR, "Fix wall/reflect/stochastic cannot be used with " + "triclinic simulation box"); + + dynamic_group_allow = 1; + + // parse args + + int arginc; + + nwall = 0; + int scaleflag = 1; + reflectionstyle = NONE; + + if (strcmp(arg[3],"diffusive") == 0) { + reflectionstyle = DIFFUSIVE; + arginc = 6; + } else if (strcmp(arg[3],"maxwell") == 0) { + reflectionstyle = MAXWELL; + arginc = 7; + } else if (strcmp(arg[3],"cercignanilampis") == 0) { + reflectionstyle = CERCIGNANILAMPIS; + arginc = 9; + } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + + + seedfix = force->inumeric(FLERR,arg[4]); + if (seedfix <= 0) error->all(FLERR,"Random seed must be a postive number"); + + int iarg = 5; + while (iarg < narg) { + 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(FLERR,"Illegal fix wall/reflect/stochastic 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 < 6); m++) + if (newwall == wallwhich[m]) + error->all(FLERR, + "Fix wall/reflect/stochastic command wall defined twice"); + + 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 { + wallstyle[nwall] = CONSTANT; + coord0[nwall] = force->numeric(FLERR,arg[iarg+1]); + } + + walltemp[nwall]= force->numeric(FLERR,arg[iarg+2]); + + for (int dir = 0; dir < 3; dir++) { + wallvel[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+3]); + int dim = wallwhich[nwall] / 2; + if ((wallvel[nwall][dir] !=0) & (dir == dim)) + error->all(FLERR,"The wall velocity must be tangential"); + + if (reflectionstyle == CERCIGNANILAMPIS) { + wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+dir+6]); + } else if (reflectionstyle == MAXWELL) { + // one accommodation coefficient for all directions + wallaccom[nwall][dir]= force->numeric(FLERR,arg[iarg+6]); + } + } + + nwall++; + iarg += arginc; + + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) + error->all(FLERR,"Illegal wall/reflect/stochastic command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + iarg += 2; + } else error->all(FLERR,"Illegal fix wall/reflect/stochastic command"); + } + + // error check + + if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); + + for (int m = 0; m < nwall; m++) { + if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) + error->all(FLERR,"Cannot use fix wall/reflect/stochastic " + "in periodic dimension"); + if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) + error->all(FLERR,"Cannot use fix wall/reflect/stochastic " + "in periodic dimension"); + if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) + error->all(FLERR,"Cannot use fix wall/reflect/stochastic " + "in periodic dimension"); + } + + for (int m = 0; m < nwall; m++) + if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->dimension == 2) + error->all(FLERR, + "Cannot use fix wall/reflect/stochastic zlo/zhi " + "for a 2d simulation"); + + // scale factors for CONSTANT walls, VARIABLE wall is not allowed + + int flag = 0; + for (int m = 0; m < nwall; m++) + if (wallstyle[m] != EDGE) flag = 1; + + if (flag) { + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + for (int m = 0; m < nwall; m++) { + if (wallstyle[m] != CONSTANT) continue; + if (wallwhich[m] < YLO) coord0[m] *= xscale; + else if (wallwhich[m] < ZLO) coord0[m] *= yscale; + else coord0[m] *= zscale; + } + } + + // random number generator + + random = new RanMars(lmp,seedfix + comm->me); +} + +/* ---------------------------------------------------------------------- */ + +FixWallReflectStochastic::~FixWallReflectStochastic() +{ + delete random; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallReflectStochastic::wall_particle(int m, int which, double coord) +{ + int i, dir, dim, side, sign; + double vsave,factor,timecol,difftest,theta; + + double *rmass; + double *mass = atom->mass; + + // coord = current position of wall + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + dim = which / 2; + side = which % 2; + + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + + sign = 0; + if ((side == 0) & (x[i][dim] < coord)) sign = 1; + else if ((side == 1) & (x[i][dim] > coord)) sign = -1; + if (sign == 0) continue; + + // theta = kT/m + + if (rmass) theta = force->boltz*walltemp[m]/(rmass[i]*force->mvv2e); + else theta = force->boltz*walltemp[m]/(mass[type[i]]*force->mvv2e); + factor = sqrt(theta); + + // time travelling after collision + + timecol = (x[i][dim]-coord) / v[i][dim]; + + // only needed for Maxwell model + + if (reflectionstyle == MAXWELL) difftest = random->uniform(); + + for (dir = 0; dir < 3; dir++) { + + // pushing back atoms to wall position, assign new reflected velocity + + x[i][dir] -= v[i][dir]*timecol; + + // diffusive reflection + + if (reflectionstyle == DIFFUSIVE) { + if (dir != dim) + v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); + else v[i][dir] = sign*random->rayleigh(factor); + + // Maxwell reflection + + } else if (reflectionstyle == MAXWELL) { + if (difftest < wallaccom[m][dir]) { + if (dir != dim) + v[i][dir] = wallvel[m][dir] + random->gaussian(0,factor); + else v[i][dir] = sign*random->rayleigh(factor); + } else { + if (dir == dim) v[i][dir] = -v[i][dir]; + } + + // Cercignani Lampis reflection + + } else if (reflectionstyle == CERCIGNANILAMPIS) { + if (dir != dim) + v[i][dir] = wallvel[m][dir] + + random->gaussian((1-wallaccom[m][dir]) * + (v[i][dir] - wallvel[m][dir]), + sqrt((2-wallaccom[m][dir]) * + wallaccom[m][dir]*theta)); + else v[i][dir] = random->besselexp(theta,wallaccom[m][dir],v[i][dir]); + } + + // update new position due to the new velocity + + if (dir != dim) x[i][dir] += v[i][dir]*timecol; + else x[i][dir] = coord + v[i][dir]*timecol; + } + } +} diff --git a/src/USER-MISC/fix_wall_reflect_stochastic.h b/src/USER-MISC/fix_wall_reflect_stochastic.h index c6957bb533..a253daea7a 100644 --- a/src/USER-MISC/fix_wall_reflect_stochastic.h +++ b/src/USER-MISC/fix_wall_reflect_stochastic.h @@ -1,85 +1,80 @@ -/* -*- c++ -*- ---------------------------------------------------------- - 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/reflect/stochastic,FixWallReflectStochastic) - -#else - -#ifndef LMP_FIX_WALL_REFLECT_STOCHASTIC_H -#define LMP_FIX_WALL_REFLECT_STOCHASTIC_H - - -#include "random_mars.h" -#include "fix_wall_reflect.h" - -namespace LAMMPS_NS { - -class FixWallReflectStochastic : public FixWallReflect { - public: - FixWallReflectStochastic(class LAMMPS *, int, char **); - void wall_particle(int m,int which, double coord); - virtual ~FixWallReflectStochastic(); - - - - - - protected: - class RanMars *random; - int seedfix; - double walltemp[6],wallvel[6][3],wallaccom[6][3]; - int reflectionstyle; - - -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Wall defined twice in fix wall/stochastic command - -Self-explanatory. - -E: Cannot use fix wall/stochastic in periodic dimension - -Self-explanatory. - -E: Cannot use fix wall/stochastic zlo/zhi for a 2d simulation - -Self-explanatory. - -E: Variable name for fix wall/stochastic does not exist - -Self-explanatory. - -E: Variable for fix wall/stochastic is invalid style - -Only equal-style variables can be used. - -W: Should not allow rigid bodies to bounce off relecting walls - -LAMMPS allows this, but their dynamics are not computed correctly. - -*/ +/* -*- c++ -*- ---------------------------------------------------------- + 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/reflect/stochastic,FixWallReflectStochastic) + +#else + +#ifndef LMP_FIX_WALL_REFLECT_STOCHASTIC_H +#define LMP_FIX_WALL_REFLECT_STOCHASTIC_H + +#include "random_mars.h" +#include "fix_wall_reflect.h" + +namespace LAMMPS_NS { + +class FixWallReflectStochastic : public FixWallReflect { + public: + FixWallReflectStochastic(class LAMMPS *, int, char **); + virtual ~FixWallReflectStochastic(); + + private: + int seedfix; + double walltemp[6],wallvel[6][3],wallaccom[6][3]; + int reflectionstyle; + + class RanMars *random; + + void wall_particle(int m,int which, double coord); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Wall defined twice in fix wall/stochastic command + +Self-explanatory. + +E: Cannot use fix wall/stochastic in periodic dimension + +Self-explanatory. + +E: Cannot use fix wall/stochastic zlo/zhi for a 2d simulation + +Self-explanatory. + +E: Variable name for fix wall/stochastic does not exist + +Self-explanatory. + +E: Variable for fix wall/stochastic is invalid style + +Only equal-style variables can be used. + +W: Should not allow rigid bodies to bounce off relecting walls + +LAMMPS allows this, but their dynamics are not computed correctly. + +*/ diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index 48e67fe33f..988179d06f 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -1,250 +1,243 @@ - /* ---------------------------------------------------------------------- - 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 "fix_wall_reflect.h" -#include -#include "atom.h" -#include "comm.h" -#include "update.h" -#include "modify.h" -#include "domain.h" -#include "lattice.h" -#include "input.h" -#include "variable.h" -#include "error.h" -#include "force.h" - - -using namespace LAMMPS_NS; -using namespace FixConst; - - - -/* ---------------------------------------------------------------------- */ - -FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - nwall(0) -{ - - - if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command"); - - if (strcmp(arg[2],"wall/reflect/stochastic") == 0) return;// child class can be stochastic - - - dynamic_group_allow = 1; - - // parse args - - nwall = 0; - int scaleflag = 1; - - int iarg = 3; - while (iarg < narg) { - 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(FLERR,"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 < 6); m++) - if (newwall == wallwhich[m]) - error->all(FLERR,"Wall defined twice in fix wall/reflect 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] = force->numeric(FLERR,arg[iarg+1]); - } - - nwall++; - iarg += 2; - - - } else if (strcmp(arg[iarg],"units") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command"); - if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; - else error->all(FLERR,"Illegal fix wall/reflect command"); - iarg += 2; - } else error->all(FLERR,"Illegal fix wall/reflect command"); - } - - // error check - - if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); - - for (int m = 0; m < nwall; m++) { - if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) - error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); - if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) - error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); - if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) - error->all(FLERR,"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(FLERR, - "Cannot use fix wall/reflect zlo/zhi for a 2d simulation"); - - // scale factors for CONSTANT and VARIABLE walls - - int flag = 0; - for (int m = 0; m < nwall; m++) - if (wallstyle[m] != EDGE) flag = 1; - - if (flag) { - if (scaleflag) { - xscale = domain->lattice->xlattice; - yscale = domain->lattice->ylattice; - zscale = domain->lattice->zlattice; - } - else xscale = yscale = zscale = 1.0; - - for (int m = 0; m < nwall; m++) { - if (wallstyle[m] != CONSTANT) continue; - if (wallwhich[m] < YLO) coord0[m] *= xscale; - else if (wallwhich[m] < ZLO) coord0[m] *= yscale; - else coord0[m] *= zscale; - } - } - - // set varflag if any wall positions are variable - - varflag = 0; - for (int m = 0; m < nwall; m++) - if (wallstyle[m] == VARIABLE) varflag = 1; -} - - - -/* ---------------------------------------------------------------------- */ - -FixWallReflect::~FixWallReflect() -{ - if (copymode) return; - - for (int m = 0; m < nwall; m++) - if (wallstyle[m] == VARIABLE) delete [] varstr[m]; -} - -/* ---------------------------------------------------------------------- */ - -int FixWallReflect::setmask() -{ - int mask = 0; - mask |= POST_INTEGRATE; - mask |= POST_INTEGRATE_RESPA; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixWallReflect::init() -{ - 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(FLERR,"Variable name for fix wall/reflect does not exist"); - if (!input->variable->equalstyle(varindex[m])) - error->all(FLERR,"Variable for fix wall/reflect is invalid style"); - } - - int nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - - if (nrigid && comm->me == 0) - error->warning(FLERR,"Should not allow rigid bodies to bounce off " - "relecting walls"); -} - -/* ---------------------------------------------------------------------- */ - -void FixWallReflect::post_integrate() -{ - //int m; - double coord; - - // coord = current position of wall - // evaluate variable if necessary, wrap with clear/add - - - - if (varflag) modify->clearstep_compute(); - - for (int m = 0; m < nwall; m++) { - if (wallstyle[m] == VARIABLE) { - coord = input->variable->compute_equal(varindex[m]); - if (wallwhich[m] < YLO) coord *= xscale; - else if (wallwhich[m] < ZLO) coord *= yscale; - else coord *= zscale; - } else coord = coord0[m]; - - - wall_particle(m,wallwhich[m],coord); - - - if (varflag) modify->addstep_compute(update->ntimestep + 1); -} -} - -void FixWallReflect::wall_particle(int m, int which, double coord) -{ - int i, dim, side,sign; - - double **x = atom->x; - double **v = atom->v; - int *mask = atom->mask; - int nlocal = atom->nlocal; - - dim = which / 2; - side = which % 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][dim] > coord) { - x[i][dim] = coord - (x[i][dim] - coord); - v[i][dim] = -v[i][dim]; - } - } - } -} \ No newline at end of file + /* ---------------------------------------------------------------------- + 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 "fix_wall_reflect.h" +#include +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "modify.h" +#include "domain.h" +#include "force.h" +#include "lattice.h" +#include "input.h" +#include "variable.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixWallReflect::FixWallReflect(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), + nwall(0) +{ + if (narg < 4) error->all(FLERR,"Illegal fix wall/reflect command"); + + // let child class process all args + + if (strcmp(arg[2],"wall/reflect/stochastic") == 0) return; + + dynamic_group_allow = 1; + + // parse args + + nwall = 0; + int scaleflag = 1; + + int iarg = 3; + while (iarg < narg) { + 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(FLERR,"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 < 6); m++) + if (newwall == wallwhich[m]) + error->all(FLERR,"Wall defined twice in fix wall/reflect 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] = force->numeric(FLERR,arg[iarg+1]); + } + + nwall++; + iarg += 2; + + } else if (strcmp(arg[iarg],"units") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal wall/reflect command"); + if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0; + else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1; + else error->all(FLERR,"Illegal fix wall/reflect command"); + iarg += 2; + + } else error->all(FLERR,"Illegal fix wall/reflect command"); + } + + // error check + + if (nwall == 0) error->all(FLERR,"Illegal fix wall command"); + + for (int m = 0; m < nwall; m++) { + if ((wallwhich[m] == XLO || wallwhich[m] == XHI) && domain->xperiodic) + error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); + if ((wallwhich[m] == YLO || wallwhich[m] == YHI) && domain->yperiodic) + error->all(FLERR,"Cannot use fix wall/reflect in periodic dimension"); + if ((wallwhich[m] == ZLO || wallwhich[m] == ZHI) && domain->zperiodic) + error->all(FLERR,"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(FLERR, + "Cannot use fix wall/reflect zlo/zhi for a 2d simulation"); + + // scale factors for CONSTANT and VARIABLE walls + + int flag = 0; + for (int m = 0; m < nwall; m++) + if (wallstyle[m] != EDGE) flag = 1; + + if (flag) { + if (scaleflag) { + xscale = domain->lattice->xlattice; + yscale = domain->lattice->ylattice; + zscale = domain->lattice->zlattice; + } + else xscale = yscale = zscale = 1.0; + + for (int m = 0; m < nwall; m++) { + if (wallstyle[m] != CONSTANT) continue; + if (wallwhich[m] < YLO) coord0[m] *= xscale; + else if (wallwhich[m] < ZLO) coord0[m] *= yscale; + else coord0[m] *= zscale; + } + } + + // set varflag if any wall positions are variable + + varflag = 0; + for (int m = 0; m < nwall; m++) + if (wallstyle[m] == VARIABLE) varflag = 1; +} + +/* ---------------------------------------------------------------------- */ + +FixWallReflect::~FixWallReflect() +{ + if (copymode) return; + + for (int m = 0; m < nwall; m++) + if (wallstyle[m] == VARIABLE) delete [] varstr[m]; +} + +/* ---------------------------------------------------------------------- */ + +int FixWallReflect::setmask() +{ + int mask = 0; + mask |= POST_INTEGRATE; + mask |= POST_INTEGRATE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallReflect::init() +{ + 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(FLERR,"Variable name for fix wall/reflect does not exist"); + if (!input->variable->equalstyle(varindex[m])) + error->all(FLERR,"Variable for fix wall/reflect is invalid style"); + } + + int nrigid = 0; + for (int i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) nrigid++; + + if (nrigid && comm->me == 0) + error->warning(FLERR,"Should not allow rigid bodies to bounce off " + "relecting walls"); +} + +/* ---------------------------------------------------------------------- */ + +void FixWallReflect::post_integrate() +{ + double coord; + + // coord = current position of wall + // evaluate variable if necessary, wrap with clear/add + + if (varflag) modify->clearstep_compute(); + + for (int m = 0; m < nwall; m++) { + if (wallstyle[m] == VARIABLE) { + coord = input->variable->compute_equal(varindex[m]); + if (wallwhich[m] < YLO) coord *= xscale; + else if (wallwhich[m] < ZLO) coord *= yscale; + else coord *= zscale; + } else coord = coord0[m]; + + wall_particle(m,wallwhich[m],coord); + } + + if (varflag) modify->addstep_compute(update->ntimestep + 1); +} + +/* ---------------------------------------------------------------------- + this method may be overwritten by a child class +------------------------------------------------------------------------- */ + +void FixWallReflect::wall_particle(int m, int which, double coord) +{ + int i,dim,side; + + double **x = atom->x; + double **v = atom->v; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + dim = which / 2; + side = which % 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][dim] > coord) { + x[i][dim] = coord - (x[i][dim] - coord); + v[i][dim] = -v[i][dim]; + } + } + } + } +} diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h index 992eabbd09..0917cafd7b 100644 --- a/src/fix_wall_reflect.h +++ b/src/fix_wall_reflect.h @@ -26,14 +26,14 @@ namespace LAMMPS_NS { class FixWallReflect : public Fix { public: + enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5}; + enum{NONE=0,EDGE,CONSTANT,VARIABLE}; + FixWallReflect(class LAMMPS *, int, char **); virtual ~FixWallReflect(); - virtual void wall_particle(int m, int which, double coord); int setmask(); void init(); void post_integrate(); - enum{XLO=0,XHI=1,YLO=2,YHI=3,ZLO=4,ZHI=5}; - enum{NONE=0,EDGE,CONSTANT,VARIABLE}; protected: int nwall; @@ -43,6 +43,8 @@ class FixWallReflect : public Fix { int varindex[6]; int varflag; double xscale,yscale,zscale; + + virtual void wall_particle(int m, int which, double coord); }; } diff --git a/src/random_mars.cpp b/src/random_mars.cpp index c458913077..227a790a62 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -90,7 +90,6 @@ double RanMars::uniform() return uni; } - /* ---------------------------------------------------------------------- gaussian RN ------------------------------------------------------------------------- */ @@ -127,7 +126,6 @@ double RanMars::gaussian(double mu, double sigma) return v1; } - /* ---------------------------------------------------------------------- Rayleigh RN ------------------------------------------------------------------------- */ @@ -136,15 +134,12 @@ double RanMars::rayleigh(double sigma) { double first,v1; - if (sigma <= 0) { - error->all(FLERR,"Invalid Rayleigh parameter"); - } else { - v1 = uniform(); - first = sigma*sqrt(-2.0*log(v1)); - return first; - } -} + if (sigma <= 0) error->all(FLERR,"Invalid Rayleigh parameter"); + v1 = uniform(); + first = sigma*sqrt(-2.0*log(v1)); + return first; +} /* ---------------------------------------------------------------------- Bessel exponential RN @@ -153,20 +148,21 @@ double RanMars::rayleigh(double sigma) double RanMars::besselexp(double theta, double alpha, double cp) { double first,v1,v2; - if ((theta < 0) || (alpha < 0) || (alpha >1)) { + + if (theta < 0.0 || alpha < 0.0 || alpha < 1.0) error->all(FLERR,"Invalid Bessel exponential distribution parameters"); - } else { - v1 = uniform(); - v2 = uniform(); - if (cp < 0) { - first = sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) - + 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) - * cos(2.0*MathConst::MY_PI*v2)*cp); - } else { - first = -sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) - - 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) - * cos(2.0*MathConst::MY_PI*v2)*cp); - } - return first; - } + + v1 = uniform(); + v2 = uniform(); + + if (cp < 0.0) + first = sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) + + 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) * + cos(2.0*MathConst::MY_PI*v2)*cp); + else + first = -sqrt((1.0-alpha)*cp*cp - 2.0*alpha*theta*log(v1) - + 2.0*sqrt(-2.0*theta*(1.0-alpha)*alpha*log(v1)) * + cos(2.0*MathConst::MY_PI*v2)*cp); + + return first; }