From 12ffa868b3c4dacec847970203417d5a81b8a94b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 24 Sep 2019 12:08:16 -0400 Subject: [PATCH] strip CR characters from source files --- doc/src/fix_wall_stochastic.txt | 214 ++++++------ src/USER-MISC/fix_wall_stochastic.h | 168 +++++----- src/fix_wall_reflect.cpp | 498 ++++++++++++++-------------- src/fix_wall_reflect.h | 170 +++++----- src/random_mars.cpp | 336 +++++++++---------- src/random_mars.h | 100 +++--- 6 files changed, 743 insertions(+), 743 deletions(-) diff --git a/doc/src/fix_wall_stochastic.txt b/doc/src/fix_wall_stochastic.txt index 876a329124..0e54eac847 100644 --- a/doc/src/fix_wall_stochastic.txt +++ b/doc/src/fix_wall_stochastic.txt @@ -1,107 +1,107 @@ -"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c - -:link(lws,http://lammps.sandia.gov) link(ld,Manual.html) -:link(lc,Section_commands.html#comm) - -:line - -fix wall/stochastic command :h3 - -[Syntax:] - -fix ID group-ID wall/stochastic reflectionstyle face pos temp velx vely velz accomx accomy accomz ... :pre - -ID, group-ID are documented in "fix"_fix.html command :ulb,l - -wall/stochastic = style name of this fix command :l - -one or more face/arg pairs may be appended :l - -reflectionstyle = diffusive or maxwell or cercignanilampis - -face = {xlo} or {xhi} or {ylo} or {yhi} or {zlo} or {zhi} :l - - {xlo},{ylo},{zlo} pos = EDGE or constant - - EDGE = current lo edge of simulation box - - constant = number like 0.0 or -30.0 (distance units) - - temp = wall temperature - - velx, vely, velz = wall velocity along x, y, z directions - - accomx, accomy, accomz = accommodation coefficients along x, y, z directions (only for maxwell and cercignanilampis reflection style) - - {xhi},{yhi},{zhi} pos = EDGE or constant or variable - - EDGE = current hi edge of simulation box - - constant = number like 50.0 or 100.3 (distance units) - - temp = wall temperature - - velx, vely, velz = wall velocity along x, y, z directions - - accomx, accomy, accomz = accommodation coefficients along x, y, z directions (only for maxwell and cercignanilampis style) - - :pre - -zero or more keyword/value pairs may be appended :l keyword = {units} :l {units} -value = {lattice} or {box} - - {lattice} = the wall position is defined in lattice units - - {box} = the wall position is defined in simulation box units :pre - -:ule - -[Examples:] - -fix zwalls all wall/stochastic diffusive zlo EDGE 300 0.1 0.1 0 zhi EDGE 200 0.1 0.1 0 - -fix ywalls all wall/stochastic maxwell ylo 5.0 300 0.1 0.0 0.0 0.8 yhi 10.0 300 0.1 0.0 0.0 0.8 - -fix xwalls all wall/stochastic cercignanilampis xlo 0.0 300 0.0 0.1 0.9 0.8 0.7 xhi EDGE 300 0.0 0.1 0 0.9 0.8 0.7 units box :pre - -[Description:] - -Different from the command fix wall/reflect which is equivalent to the mirror reflection, the post collision velocity of the atoms in the stochastic wall is random. The randomness can come from many sources: thermal motion of solid atoms, the surface roughness, etc... Three stochastic wall models are implemented. - -Diffusive wall (reflectionstyle = diffusive): all the gas atoms are reflected diffusively. Their velocity distribution corresponds to the equilibrium distribution at the wall temperature. - -Maxwell model (reflectionstyle = maxwell): the reflection is partially diffusive and partially specular. The accommodation coefficient is the portion of the diffusive reflection "Maxwell"_#Maxwell. - -Cercignani Lampis model: 3 accommodations coefficients are used, two velocity accommodation coefficients and one normal kinetic energy accommodation "CL"_#CL, "To"_#To. The code will recognize automatically the normal direction of the wall and the other two tangential directions. For example, if accomx = 0.5, accomy = 0.2, accomz = 0.3 and the wall is normal to direction y, accomx and accomz are the tangential velocity accommodation coefficients and accomy is the normal kinetic energy accommodation coefficients. - -The {units} keyword determines the distance units used to define a wall position - -A {box} value selects standard distance units as defined by the -"units"_units.html command, e.g. Angstroms for units = real or metal. - -A {lattice} value means the distance units are in lattice spacings. The -"lattice"_lattice.html command must have been previously used to define the -lattice spacings. - -:line - -[Restrictions:] - -Due to the similarities, this fix has the same limitation as the the -wall/reflect. Any dimension (xyz) that has a wall must be non-periodic. It -shoudld not be used with rigid bodies such as those defined by a "fix rigid" -command. Furthermore, the wall velocity must lie on the same plane as the wall. Package USER-MISC is required to use this fix. - -[Related commands:] - -"fix wall/reflect"_fix_wall_reflect.html, - -[Default:] none - -:line -:link(Maxwell) [(Maxwell)] J. C. Maxwell. Philos. Trans. Royal Soc. London. 157: 49–88 (1867). - -:link(CL) [(CL)] C. Cercignani and M. Lampis. Transp. Theory Stat. Phys. 1, 2, 101 (1971). - -:link(To) [(To)] Q.D. To, V.H. Vu, G. Lauriat, and C. Leonard. J. Math. Phys. 56, 103101 (2015). - +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix wall/stochastic command :h3 + +[Syntax:] + +fix ID group-ID wall/stochastic reflectionstyle face pos temp velx vely velz accomx accomy accomz ... :pre + +ID, group-ID are documented in "fix"_fix.html command :ulb,l + +wall/stochastic = style name of this fix command :l + +one or more face/arg pairs may be appended :l + +reflectionstyle = diffusive or maxwell or cercignanilampis + +face = {xlo} or {xhi} or {ylo} or {yhi} or {zlo} or {zhi} :l + + {xlo},{ylo},{zlo} pos = EDGE or constant + + EDGE = current lo edge of simulation box + + constant = number like 0.0 or -30.0 (distance units) + + temp = wall temperature + + velx, vely, velz = wall velocity along x, y, z directions + + accomx, accomy, accomz = accommodation coefficients along x, y, z directions (only for maxwell and cercignanilampis reflection style) + + {xhi},{yhi},{zhi} pos = EDGE or constant or variable + + EDGE = current hi edge of simulation box + + constant = number like 50.0 or 100.3 (distance units) + + temp = wall temperature + + velx, vely, velz = wall velocity along x, y, z directions + + accomx, accomy, accomz = accommodation coefficients along x, y, z directions (only for maxwell and cercignanilampis style) + + :pre + +zero or more keyword/value pairs may be appended :l keyword = {units} :l {units} +value = {lattice} or {box} + + {lattice} = the wall position is defined in lattice units + + {box} = the wall position is defined in simulation box units :pre + +:ule + +[Examples:] + +fix zwalls all wall/stochastic diffusive zlo EDGE 300 0.1 0.1 0 zhi EDGE 200 0.1 0.1 0 + +fix ywalls all wall/stochastic maxwell ylo 5.0 300 0.1 0.0 0.0 0.8 yhi 10.0 300 0.1 0.0 0.0 0.8 + +fix xwalls all wall/stochastic cercignanilampis xlo 0.0 300 0.0 0.1 0.9 0.8 0.7 xhi EDGE 300 0.0 0.1 0 0.9 0.8 0.7 units box :pre + +[Description:] + +Different from the command fix wall/reflect which is equivalent to the mirror reflection, the post collision velocity of the atoms in the stochastic wall is random. The randomness can come from many sources: thermal motion of solid atoms, the surface roughness, etc... Three stochastic wall models are implemented. + +Diffusive wall (reflectionstyle = diffusive): all the gas atoms are reflected diffusively. Their velocity distribution corresponds to the equilibrium distribution at the wall temperature. + +Maxwell model (reflectionstyle = maxwell): the reflection is partially diffusive and partially specular. The accommodation coefficient is the portion of the diffusive reflection "Maxwell"_#Maxwell. + +Cercignani Lampis model: 3 accommodations coefficients are used, two velocity accommodation coefficients and one normal kinetic energy accommodation "CL"_#CL, "To"_#To. The code will recognize automatically the normal direction of the wall and the other two tangential directions. For example, if accomx = 0.5, accomy = 0.2, accomz = 0.3 and the wall is normal to direction y, accomx and accomz are the tangential velocity accommodation coefficients and accomy is the normal kinetic energy accommodation coefficients. + +The {units} keyword determines the distance units used to define a wall position + +A {box} value selects standard distance units as defined by the +"units"_units.html command, e.g. Angstroms for units = real or metal. + +A {lattice} value means the distance units are in lattice spacings. The +"lattice"_lattice.html command must have been previously used to define the +lattice spacings. + +:line + +[Restrictions:] + +Due to the similarities, this fix has the same limitation as the the +wall/reflect. Any dimension (xyz) that has a wall must be non-periodic. It +shoudld not be used with rigid bodies such as those defined by a "fix rigid" +command. Furthermore, the wall velocity must lie on the same plane as the wall. Package USER-MISC is required to use this fix. + +[Related commands:] + +"fix wall/reflect"_fix_wall_reflect.html, + +[Default:] none + +:line +:link(Maxwell) [(Maxwell)] J. C. Maxwell. Philos. Trans. Royal Soc. London. 157: 49–88 (1867). + +:link(CL) [(CL)] C. Cercignani and M. Lampis. Transp. Theory Stat. Phys. 1, 2, 101 (1971). + +:link(To) [(To)] Q.D. To, V.H. Vu, G. Lauriat, and C. Leonard. J. Math. Phys. 56, 103101 (2015). + diff --git a/src/USER-MISC/fix_wall_stochastic.h b/src/USER-MISC/fix_wall_stochastic.h index 9eabc8fdca..6fc56f0b6e 100644 --- a/src/USER-MISC/fix_wall_stochastic.h +++ b/src/USER-MISC/fix_wall_stochastic.h @@ -1,84 +1,84 @@ -/* -*- 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/stochastic,FixWallStochastic) - -#else - -#ifndef LMP_FIX_WALL_STOCHASTIC_H -#define LMP_FIX_WALL_STOCHASTIC_H - - -#include "random_mars.h" -#include "fix_wall_reflect.h" - -namespace LAMMPS_NS { - -class FixWallStochastic : public FixWallReflect { - public: - FixWallStochastic(class LAMMPS *, int, char **); - void wall_particle(int m,int which, double coord); - - - - - - 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/stochastic,FixWallStochastic) + +#else + +#ifndef LMP_FIX_WALL_STOCHASTIC_H +#define LMP_FIX_WALL_STOCHASTIC_H + + +#include "random_mars.h" +#include "fix_wall_reflect.h" + +namespace LAMMPS_NS { + +class FixWallStochastic : public FixWallReflect { + public: + FixWallStochastic(class LAMMPS *, int, char **); + void wall_particle(int m,int which, double coord); + + + + + + 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. + +*/ diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index cffc97413a..b6a2f214c1 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -1,250 +1,250 @@ - /* ---------------------------------------------------------------------- - 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") == 0){ // 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]; - } - } - } + /* ---------------------------------------------------------------------- + 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") == 0){ // 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 diff --git a/src/fix_wall_reflect.h b/src/fix_wall_reflect.h index cb549e218b..992eabbd09 100644 --- a/src/fix_wall_reflect.h +++ b/src/fix_wall_reflect.h @@ -1,85 +1,85 @@ -/* -*- 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,FixWallReflect) - -#else - -#ifndef LMP_FIX_WALL_REFLECT_H -#define LMP_FIX_WALL_REFLECT_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixWallReflect : public Fix { - public: - 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; - int wallwhich[6],wallstyle[6]; - double coord0[6]; - char *varstr[6]; - int varindex[6]; - int varflag; - double xscale,yscale,zscale; -}; - -} - -#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/reflect command - -Self-explanatory. - -E: Cannot use fix wall/reflect in periodic dimension - -Self-explanatory. - -E: Cannot use fix wall/reflect zlo/zhi for a 2d simulation - -Self-explanatory. - -E: Variable name for fix wall/reflect does not exist - -Self-explanatory. - -E: Variable for fix wall/reflect 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,FixWallReflect) + +#else + +#ifndef LMP_FIX_WALL_REFLECT_H +#define LMP_FIX_WALL_REFLECT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixWallReflect : public Fix { + public: + 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; + int wallwhich[6],wallstyle[6]; + double coord0[6]; + char *varstr[6]; + int varindex[6]; + int varflag; + double xscale,yscale,zscale; +}; + +} + +#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/reflect command + +Self-explanatory. + +E: Cannot use fix wall/reflect in periodic dimension + +Self-explanatory. + +E: Cannot use fix wall/reflect zlo/zhi for a 2d simulation + +Self-explanatory. + +E: Variable name for fix wall/reflect does not exist + +Self-explanatory. + +E: Variable for fix wall/reflect 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/random_mars.cpp b/src/random_mars.cpp index 633d57e062..0ee62db9d6 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -1,169 +1,169 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -// Marsaglia random number generator -// see RANMAR in F James, Comp Phys Comm, 60, 329 (1990) - -#include "random_mars.h" -#include -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp), - u(NULL) -{ - int ij,kl,i,j,k,l,ii,jj,m; - double s,t; - - if (seed <= 0 || seed > 900000000) - error->one(FLERR,"Invalid seed for Marsaglia random # generator"); - - save = 0; - u = new double[97+1]; - - ij = (seed-1)/30082; - kl = (seed-1) - 30082*ij; - i = (ij/177) % 177 + 2; - j = ij %177 + 2; - k = (kl/169) % 178 + 1; - l = kl % 169; - for (ii = 1; ii <= 97; ii++) { - s = 0.0; - t = 0.5; - for (jj = 1; jj <= 24; jj++) { - m = ((i*j) % 179)*k % 179; - i = j; - j = k; - k = m; - l = (53*l+1) % 169; - if ((l*m) % 64 >= 32) s = s + t; - t = 0.5*t; - } - u[ii] = s; - } - c = 362436.0 / 16777216.0; - cd = 7654321.0 / 16777216.0; - cm = 16777213.0 / 16777216.0; - i97 = 97; - j97 = 33; - uniform(); -} - -/* ---------------------------------------------------------------------- */ - -RanMars::~RanMars() -{ - delete [] u; -} - -/* ---------------------------------------------------------------------- - uniform RN -------------------------------------------------------------------------- */ - -double RanMars::uniform() -{ - double uni = u[i97] - u[j97]; - if (uni < 0.0) uni += 1.0; - u[i97] = uni; - i97--; - if (i97 == 0) i97 = 97; - j97--; - if (j97 == 0) j97 = 97; - c -= cd; - if (c < 0.0) c += cm; - uni -= c; - if (uni < 0.0) uni += 1.0; - return uni; -} - - -/* ---------------------------------------------------------------------- - gaussian RN -------------------------------------------------------------------------- */ - -double RanMars::gaussian() -{ - double first,v1,v2,rsq,fac; - - if (!save) { - do { - v1 = 2.0*uniform()-1.0; - v2 = 2.0*uniform()-1.0; - rsq = v1*v1 + v2*v2; - } while ((rsq >= 1.0) || (rsq == 0.0)); - fac = sqrt(-2.0*log(rsq)/rsq); - second = v1*fac; - first = v2*fac; - save = 1; - } else { - first = second; - save = 0; - } - return first; -} - -/* ---------------------------------------------------------------------- - Gaussian RN -------------------------------------------------------------------------- */ - -double RanMars::gaussian(double mu, double sigma) -{ - double v1; - v1 = mu+sigma*gaussian(); - return v1; -} -/* ---------------------------------------------------------------------- */ - - -/* ---------------------------------------------------------------------- - Rayleigh RN -------------------------------------------------------------------------- */ - -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; - } -} - - -/* ---------------------------------------------------------------------- - Bessel exponential RN -------------------------------------------------------------------------- */ - -double RanMars::besselexp(double theta, double alpha, double cp) -{ - double first,v1,v2; - - - if ((theta < 0) || (alpha < 0) || (alpha >1)) - error->all(FLERR,"Invalid Bessel exponential distribution parameters"); - else { - v1 = uniform(); - v2 = uniform(); - if (cp < 0) - first = sqrt((1-alpha)*cp*cp-2*alpha*theta*log(v1)+2*sqrt(-2*theta*(1-alpha)*alpha*log(v1))*cos(2*M_PI*v2)*cp); - else { - first = - sqrt((1-alpha)*cp*cp-2*alpha*theta*log(v1)-2*sqrt(-2*theta*(1-alpha)*alpha*log(v1))*cos(2*M_PI*v2)*cp) ;} - return first; - } +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +// Marsaglia random number generator +// see RANMAR in F James, Comp Phys Comm, 60, 329 (1990) + +#include "random_mars.h" +#include +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +RanMars::RanMars(LAMMPS *lmp, int seed) : Pointers(lmp), + u(NULL) +{ + int ij,kl,i,j,k,l,ii,jj,m; + double s,t; + + if (seed <= 0 || seed > 900000000) + error->one(FLERR,"Invalid seed for Marsaglia random # generator"); + + save = 0; + u = new double[97+1]; + + ij = (seed-1)/30082; + kl = (seed-1) - 30082*ij; + i = (ij/177) % 177 + 2; + j = ij %177 + 2; + k = (kl/169) % 178 + 1; + l = kl % 169; + for (ii = 1; ii <= 97; ii++) { + s = 0.0; + t = 0.5; + for (jj = 1; jj <= 24; jj++) { + m = ((i*j) % 179)*k % 179; + i = j; + j = k; + k = m; + l = (53*l+1) % 169; + if ((l*m) % 64 >= 32) s = s + t; + t = 0.5*t; + } + u[ii] = s; + } + c = 362436.0 / 16777216.0; + cd = 7654321.0 / 16777216.0; + cm = 16777213.0 / 16777216.0; + i97 = 97; + j97 = 33; + uniform(); +} + +/* ---------------------------------------------------------------------- */ + +RanMars::~RanMars() +{ + delete [] u; +} + +/* ---------------------------------------------------------------------- + uniform RN +------------------------------------------------------------------------- */ + +double RanMars::uniform() +{ + double uni = u[i97] - u[j97]; + if (uni < 0.0) uni += 1.0; + u[i97] = uni; + i97--; + if (i97 == 0) i97 = 97; + j97--; + if (j97 == 0) j97 = 97; + c -= cd; + if (c < 0.0) c += cm; + uni -= c; + if (uni < 0.0) uni += 1.0; + return uni; +} + + +/* ---------------------------------------------------------------------- + gaussian RN +------------------------------------------------------------------------- */ + +double RanMars::gaussian() +{ + double first,v1,v2,rsq,fac; + + if (!save) { + do { + v1 = 2.0*uniform()-1.0; + v2 = 2.0*uniform()-1.0; + rsq = v1*v1 + v2*v2; + } while ((rsq >= 1.0) || (rsq == 0.0)); + fac = sqrt(-2.0*log(rsq)/rsq); + second = v1*fac; + first = v2*fac; + save = 1; + } else { + first = second; + save = 0; + } + return first; +} + +/* ---------------------------------------------------------------------- + Gaussian RN +------------------------------------------------------------------------- */ + +double RanMars::gaussian(double mu, double sigma) +{ + double v1; + v1 = mu+sigma*gaussian(); + return v1; +} +/* ---------------------------------------------------------------------- */ + + +/* ---------------------------------------------------------------------- + Rayleigh RN +------------------------------------------------------------------------- */ + +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; + } +} + + +/* ---------------------------------------------------------------------- + Bessel exponential RN +------------------------------------------------------------------------- */ + +double RanMars::besselexp(double theta, double alpha, double cp) +{ + double first,v1,v2; + + + if ((theta < 0) || (alpha < 0) || (alpha >1)) + error->all(FLERR,"Invalid Bessel exponential distribution parameters"); + else { + v1 = uniform(); + v2 = uniform(); + if (cp < 0) + first = sqrt((1-alpha)*cp*cp-2*alpha*theta*log(v1)+2*sqrt(-2*theta*(1-alpha)*alpha*log(v1))*cos(2*M_PI*v2)*cp); + else { + first = - sqrt((1-alpha)*cp*cp-2*alpha*theta*log(v1)-2*sqrt(-2*theta*(1-alpha)*alpha*log(v1))*cos(2*M_PI*v2)*cp) ;} + return first; + } } \ No newline at end of file diff --git a/src/random_mars.h b/src/random_mars.h index 5ca7cc0ee6..af653e201b 100644 --- a/src/random_mars.h +++ b/src/random_mars.h @@ -1,50 +1,50 @@ -/* -*- 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. -------------------------------------------------------------------------- */ - -#ifndef LMP_RANMARS_H -#define LMP_RANMARS_H - -#include "pointers.h" - -namespace LAMMPS_NS { - -class RanMars : protected Pointers { - public: - RanMars(class LAMMPS *, int); - ~RanMars(); - double uniform(); - double gaussian(); - double gaussian(double mu, double sigma); - double rayleigh(double sigma); - double besselexp(double theta, double alpha, double cp); - - private: - int save; - double second; - double *u; - int i97,j97; - double c,cd,cm; -}; - -} - -#endif - -/* ERROR/WARNING messages: - -E: Invalid seed for Marsaglia random # generator - -The initial seed for this random number generator must be a positive -integer less than or equal to 900 million. - -*/ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifndef LMP_RANMARS_H +#define LMP_RANMARS_H + +#include "pointers.h" + +namespace LAMMPS_NS { + +class RanMars : protected Pointers { + public: + RanMars(class LAMMPS *, int); + ~RanMars(); + double uniform(); + double gaussian(); + double gaussian(double mu, double sigma); + double rayleigh(double sigma); + double besselexp(double theta, double alpha, double cp); + + private: + int save; + double second; + double *u; + int i97,j97; + double c,cd,cm; +}; + +} + +#endif + +/* ERROR/WARNING messages: + +E: Invalid seed for Marsaglia random # generator + +The initial seed for this random number generator must be a positive +integer less than or equal to 900 million. + +*/