implement fix wall/morse based on fix wall/region and fix wall/harmonic

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

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

: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

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

#include "respa.h"

#include "error.h"

#include "force.h"

#include "utils.h"

using namespace LAMMPS_NS;

using namespace FixConst;

        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];

  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;

/* ----------------------------------------------------------------------

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

#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");

}

/* -*- 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.

*/