diff --git a/doc/src/fix_wall.txt b/doc/src/fix_wall.txt index 162e32c68c..b8e380fb16 100644 --- a/doc/src/fix_wall.txt +++ b/doc/src/fix_wall.txt @@ -12,22 +12,25 @@ fix wall/lj126 command :h3 fix wall/lj1043 command :h3 fix wall/colloid command :h3 fix wall/harmonic command :h3 +fix wall/morse command :h3 [Syntax:] fix ID group-ID style face args ... keyword value ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l -style = {wall/lj93} or {wall/lj126} or {wall/lj1043} or {wall/colloid} or {wall/harmonic} :l +style = {wall/lj93} or {wall/lj126} or {wall/lj1043} or {wall/colloid} or {wall/harmonic} or {wall/morse} :l one or more face/arg pairs may be appended :l face = {xlo} or {xhi} or {ylo} or {yhi} or {zlo} or {zhi} :l - args = coord epsilon sigma cutoff + args = coord epsilon \[alpha\] sigma cutoff coord = position of wall = EDGE or constant or variable EDGE = current lo or hi edge of simulation box constant = number like 0.0 or -30.0 (distance units) variable = "equal-style variable"_variable.html like v_x or v_wiggle epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units) epsilon can be a variable (see below) + alpha = width factor for wall-particle interaction (1/distance units) + [only] for {wall/morse}. alpha can be a variable (see below) sigma = size factor for wall-particle interaction (distance units) sigma can be a variable (see below) cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :pre @@ -48,6 +51,7 @@ keyword = {units} or {fld} :l fix wallhi all wall/lj93 xlo -1.0 1.0 1.0 2.5 units box fix wallhi all wall/lj93 xhi EDGE 1.0 1.0 2.5 +fix wallhi all wall/morse xhi EDGE 1.0 1.0 1.0 2.5 units box fix wallhi all wall/lj126 v_wiggle 23.2 1.0 1.0 2.5 fix zwalls all wall/colloid zlo 0.0 1.0 1.0 0.858 zhi 40.0 1.0 1.0 0.858 :pre @@ -80,6 +84,10 @@ potential: :c,image(Eqs/fix_wall_harmonic.jpg) +For style {wall/morse}, the energy E is given by a Morse potential: + +:c,image(Eqs/pair_morse.jpg) + In all cases, {r} is the distance from the particle to the wall at position {coord}, and Rc is the {cutoff} distance at which the particle and wall no longer interact. The energy of the wall @@ -147,7 +155,15 @@ constant K, and has units (energy/distance^2). The input parameter spring is at the {cutoff}. This is a repulsive-only spring since the interaction is truncated at the {cutoff} -For any wall, the {epsilon} and/or {sigma} parameter can be specified +For the {wall/morse} style, one additional parameter {alpha} is required. +Thus the parameters are in this order: {epsilon} as the depth of the +Morse potential (D_0), {alpha} as the width parameter of the Morse +potential, and {sigma} the location of the minimum (r_0) +the wall. {D_0} has energy units, {alpha} inverse distance units, and +{r_0} distance units. + +For any wall, the {epsilon} and/or {sigma} and/or {alpha} parameter can +be specified as an "equal-style variable"_variable.html, in which case it should be specified as v_name, where name is the variable name. As with a variable wall position, the variable is evaluated each timestep and diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index 8ceafb4bfc..a30b14c231 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -23,6 +23,7 @@ #include "respa.h" #include "error.h" #include "force.h" +#include "utils.h" using namespace LAMMPS_NS; using namespace FixConst; @@ -99,6 +100,19 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) : estyle[nwall] = CONSTANT; } + if (utils::strmatch(style,"^wall/morse")) { + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { + int n = strlen(&arg[iarg+3][2]) + 1; + astr[nwall] = new char[n]; + strcpy(astr[nwall],&arg[iarg+3][2]); + astyle[nwall] = VARIABLE; + } else { + alpha[nwall] = force->numeric(FLERR,arg[iarg+3]); + astyle[nwall] = CONSTANT; + } + ++iarg; + } + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) { int n = strlen(&arg[iarg+3][2]) + 1; sstr[nwall] = new char[n]; diff --git a/src/fix_wall.h b/src/fix_wall.h index 82634c349f..b2bdc2e27c 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -45,12 +45,12 @@ class FixWall : public Fix { virtual void wall_particle(int, int, double) = 0; protected: - double epsilon[6],sigma[6],cutoff[6]; + double epsilon[6],sigma[6],alpha[6],cutoff[6]; double ewall[7],ewall_all[7]; double xscale,yscale,zscale; - int estyle[6],sstyle[6],wstyle[6]; + int estyle[6],sstyle[6],astyle[6],wstyle[6]; int eindex[6],sindex[6]; - char *estr[6],*sstr[6]; + char *estr[6],*sstr[6],*astr[6]; int varflag; // 1 if any wall position,epsilon,sigma is a var int eflag; // per-wall flag for energy summation int ilevel_respa; diff --git a/src/fix_wall_morse.cpp b/src/fix_wall_morse.cpp new file mode 100644 index 0000000000..f854341e11 --- /dev/null +++ b/src/fix_wall_morse.cpp @@ -0,0 +1,87 @@ +/* ---------------------------------------------------------------------- + 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_morse.h" +#include +#include "atom.h" +#include "error.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixWallMorse::FixWallMorse(LAMMPS *lmp, int narg, char **arg) : + FixWall(lmp, narg, arg) +{ + dynamic_group_allow = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixWallMorse::precompute(int m) +{ + coeff1[m] = 2.0 * epsilon[m] * alpha[m]; + const double alpha_dr = -alpha[m] * (cutoff[m] - sigma[m]); + offset[m] = epsilon[m] * (exp(2.0*alpha_dr) - 2.0*exp(alpha_dr)); +} + +/* ---------------------------------------------------------------------- + interaction of all particles in group with a wall + m = index of wall coeffs + which = xlo,xhi,ylo,yhi,zlo,zhi + error if any particle is on or behind wall +------------------------------------------------------------------------- */ + +void FixWallMorse::wall_particle(int m, int which, double coord) +{ + double delta,fwall; + double vn; + + double **x = atom->x; + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + int dim = which / 2; + int side = which % 2; + if (side == 0) side = -1; + + int onflag = 0; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (side < 0) delta = x[i][dim] - coord; + else delta = coord - x[i][dim]; + if (delta >= cutoff[m]) continue; + if (delta <= 0.0) { + onflag = 1; + continue; + } + double dr = delta - sigma[m]; + double dexp = exp(-alpha[m] * dr); + fwall = side * coeff1[m] * (dexp*dexp - dexp) / delta; + ewall[0] += epsilon[m] * (dexp*dexp - 2.0*dexp) - offset[m]; + f[i][dim] -= fwall; + ewall[m+1] += fwall; + + if (evflag) { + if (side < 0) vn = -fwall*delta; + else vn = fwall*delta; + v_tally(dim, i, vn); + } + } + } + + if (onflag) error->one(FLERR,"Particle on or inside fix wall surface"); +} diff --git a/src/fix_wall_morse.h b/src/fix_wall_morse.h new file mode 100644 index 0000000000..a7359e9cde --- /dev/null +++ b/src/fix_wall_morse.h @@ -0,0 +1,49 @@ +/* -*- 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/morse,FixWallMorse) + +#else + +#ifndef LMP_FIX_WALL_MORSE_H +#define LMP_FIX_WALL_MORSE_H + +#include "fix_wall.h" + +namespace LAMMPS_NS { + +class FixWallMorse : public FixWall { + public: + FixWallMorse(class LAMMPS *, int, char **); + void precompute(int); + void wall_particle(int, int, double); + + private: + double coeff1[6],offset[6]; +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Particle on or inside fix wall surface + +Particles must be "exterior" to the wall in order for energy/force to +be calculated. + +*/