forked from lijiext/lammps
strip CR characters from source files
This commit is contained in:
Axel Kohlmeyer
parent
36f417ac14
commit
12ffa868b3
|
|
@ -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).
|
||||
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -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 <cstring>
|
||||
#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 <cstring>
|
||||
#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];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -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 <cmath>
|
||||
#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 <cmath>
|
||||
#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;
|
||||
}
|
||||
}
|
||||
|
|
@ -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.
|
||||
|
||||
*/
|
||||
|
|
|
|||
Loading…
Reference in New Issue