import changes to various fixes by aidan to include virial contributions

doc/src/fix_cmap.txt

@@ -97,6 +97,11 @@ The "fix_modify"_fix_modify.html {energy} option is supported by this
 fix to add the potential "energy" of the CMAP interactions system's
 potential energy as part of "thermodynamic output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial yes}
+
 This fix computes a global scalar which can be accessed by various
 "output commands"_Section_howto.html#howto_15.  The scalar is the
 potential energy discussed above.  The scalar value calculated by this

doc/src/fix_efield.txt

@@ -124,6 +124,11 @@ can include the forces added by this fix in a consistent manner.
 I.e. there is a decrease in potential energy when atoms move in the
 direction of the added force due to the electric field.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix adding its forces. Default is the outermost level.

doc/src/fix_external.txt

@@ -131,6 +131,11 @@ forces added by this fix in a consistent manner.  I.e. there is a
 decrease in potential energy when atoms move in the direction of the
 added force.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial yes}
+
 This fix computes a global scalar which can be accessed by various
 "output commands"_Section_howto.html#howto_15.  The scalar is the
 potential energy discussed above.  The scalar stored by this fix

doc/src/fix_modify.txt

@@ -14,10 +14,11 @@ fix_modify fix-ID keyword value ... :pre
 
 fix-ID = ID of the fix to modify :ulb,l
 one or more keyword/value pairs may be appended :l
-keyword = {temp} or {press} or {energy} or {respa} or {dynamic/dof} :l
+keyword = {temp} or {press} or {energy} or {virial} or {respa} or {dynamic/dof} :l
   {temp} value = compute ID that calculates a temperature
   {press} value = compute ID that calculates a pressure
   {energy} value = {yes} or {no}
+  {virial} value = {yes} or {no}
   {respa} value = {1} to {max respa level} or {0} (for outermost level)
   {dynamic/dof} value = {yes} or {no}
     yes/no = do or do not recompute the number of degrees of freedom (DOF) contributing to the temperature :pre
@@ -52,11 +53,10 @@ define their own compute by default, as described in their
 documentation.  Thus this option allows the user to override the
 default method for computing P.
 
-For fixes that calculate a contribution to the potential energy of the
+The {energy} keyword can be used with fixes that support it.
-system, the {energy} keyword will include that contribution in
+{energy yes} adds a contribution to the potential energy of the
-thermodynamic output of potential energy.  This is because the {energy
+system. The fix's global and per-atom
-yes} setting must be specified to include the fix's global or per-atom
+energy is included in the calculation performed by the "compute
-energy in the calculation performed by the "compute
 pe"_compute_pe.html or "compute pe/atom"_compute_pe_atom.html
 commands.  See the "thermo_style"_thermo_style.html command for info
 on how potential energy is output.  For fixes that tally a global
@@ -69,6 +69,25 @@ are using it when performing an "energy minimization"_minimize.html
 and if you want the energy and forces it produces to be part of the
 optimization criteria.
 
+The {virial} keyword can be used with fixes that support it.
+{virial yes} adds a contribution to the virial of the
+system. The fix's global and per-atom
+virial is included in the calculation performed by the "compute
+pressure"_compute_pressure.html or
+"compute stress/atom"_compute_stress_atom.html
+commands.  See the "thermo_style"_thermo_style.html command for info
+on how pressure is output. 
+
+NOTE: You must specify the {virial yes} setting for a fix if you
+are doing "box relaxation"_fix_box_relax.html and
+if you want virial contribution of the fix to be part of the
+relaxation criteria, although this seems unlikely.
+
+NOTE: This option is only supported by fixes that explicitly say
+so. For some of these (e.g. the
+"fix shake"_fix_shake.html command) the default setting is
+{virial yes}, for others it is {virial no}.
+
 For fixes that set or modify forces, it may be possible to select at
 which "r-RESPA"_run_style.html level the fix operates via the {respa}
 keyword. The RESPA level at which the fix is active can be selected.
@@ -111,4 +130,4 @@ pressure"_compute_pressure.html, "thermo_style"_thermo_style.html
 [Default:]
 
 The option defaults are temp = ID defined by fix, press = ID defined
-by fix, energy = no, respa = 0.
+by fix, energy = no, virial = different for each fix style, respa = 0.

doc/src/fix_rigid.txt

@@ -703,6 +703,11 @@ NVT, NPT, NPH rigid styles to add the energy change induced by the
 thermostatting to the system's potential energy as part of
 "thermodynamic output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial yes}
+
 The "fix_modify"_fix_modify.html {temp} and {press} options are
 supported by the 4 NPT and NPH rigid styles to change the computes
 used to calculate the instantaneous pressure tensor.  Note that the 2

doc/src/fix_shake.txt

@@ -186,6 +186,11 @@ to 1 and recompiling LAMMPS.
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial yes}
+
 No information about these fixes is written to "binary restart
 files"_restart.html.  None of the "fix_modify"_fix_modify.html options
 are relevant to these fixes.  No global or per-atom quantities are

doc/src/fix_wall.txt

@@ -252,6 +252,11 @@ fix to add the energy of interaction between atoms and each wall to
 the system's potential energy as part of "thermodynamic
 output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost level.

doc/src/fix_wall_region.txt

@@ -15,7 +15,7 @@ fix ID group-ID wall/region region-ID style epsilon sigma cutoff :pre
 ID, group-ID are documented in "fix"_fix.html command
 wall/region = style name of this fix command
 region-ID = region whose boundary will act as wall
-style = {lj93} or {lj126} or {colloid} or {harmonic}
+style = {lj93} or {lj126} or {lj1043} or {colloid} or {harmonic}
 epsilon = strength factor for wall-particle interaction (energy or energy/distance^2 units)
 sigma = size factor for wall-particle interaction (distance units)
 cutoff = distance from wall at which wall-particle interaction is cut off (distance units) :ul
@@ -112,6 +112,10 @@ For style {lj126}, the energy E is given by the 12/6 potential:
 
 :c,image(Eqs/pair_lj.jpg)
 
+For style {wall/lj1043}, the energy E is given by the 10/4/3 potential:
+
+:c,image(Eqs/fix_wall_lj1043.jpg)
+
 For style {colloid}, the energy E is given by an integrated form of
 the "pair_style colloid"_pair_colloid.html potential:
 
@@ -128,49 +132,8 @@ surface no longer interact.  The energy of the wall potential is
 shifted so that the wall-particle interaction energy is 0.0 at the
 cutoff distance.
 
-For the {lj93} and {lj126} styles, {epsilon} and {sigma} are the usual
+For a full description of these wall styles, see fix_style
-Lennard-Jones parameters, which determine the strength and size of the
+"wall"_fix_wall.html
-particle as it interacts with the wall.  Epsilon has energy units.
-Note that this {epsilon} and {sigma} may be different than any
-{epsilon} or {sigma} values defined for a pair style that computes
-particle-particle interactions.
-
-The {lj93} interaction is derived by integrating over a 3d
-half-lattice of Lennard-Jones 12/6 particles.  The {lj126} interaction
-is effectively a harder, more repulsive wall interaction.
-
-For the {colloid} style, {epsilon} is effectively a Hamaker constant
-with energy units for the colloid-wall interaction, {R} is the radius
-of the colloid particle, {D} is the distance from the surface of the
-colloid particle to the wall (r-R), and {sigma} is the size of a
-constituent LJ particle inside the colloid particle.  Note that the
-cutoff distance Rc in this case is the distance from the colloid
-particle center to the wall.
-
-The {colloid} interaction is derived by integrating over constituent
-LJ particles of size {sigma} within the colloid particle and a 3d
-half-lattice of Lennard-Jones 12/6 particles of size {sigma} in the
-wall.
-
-For the {wall/harmonic} style, {epsilon} is effectively the spring
-constant K, and has units (energy/distance^2).  The input parameter
-{sigma} is ignored.  The minimum energy position of the harmonic
-spring is at the {cutoff}.  This is a repulsive-only spring since the
-interaction is truncated at the {cutoff}
-
-NOTE: For all of the styles, you must insure that r is always > 0 for
-all particles in the group, or LAMMPS will generate an error.  This
-means you cannot start your simulation with particles on the region
-surface (r = 0) or with particles on the wrong side of the region
-surface (r < 0).  For the {wall/lj93} and {wall/lj126} styles, the
-energy of the wall/particle interaction (and hence the force on the
-particle) blows up as r -> 0.  The {wall/colloid} style is even more
-restrictive, since the energy blows up as D = r-R -> 0.  This means
-the finite-size particles of radius R must be a distance larger than R
-from the region surface.  The {harmonic} style is a softer potential
-and does not blow up as r -> 0, but you must use a large enough
-{epsilon} that particles always reamin on the correct side of the
-region surface (r > 0).
 
 [Restart, fix_modify, output, run start/stop, minimize info:]
 
@@ -182,6 +145,11 @@ fix to add the energy of interaction between atoms and the wall to the
 system's potential energy as part of "thermodynamic
 output"_thermo_style.html.
 
+The "fix_modify"_fix_modify.html {virial} option is supported by this
+fix to add the contribution due to the interaction between
+atoms and each wall to the system's virial as part of "thermodynamic
+output"_thermo_style.html. The default is {virial no}
+
 The "fix_modify"_fix_modify.html {respa} option is supported by this
 fix. This allows to set at which level of the "r-RESPA"_run_style.html
 integrator the fix is adding its forces. Default is the outermost level.

src/COLLOID/fix_wall_colloid.cpp

@@ -81,6 +81,7 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
   double r3,rinv3,r2inv3,r4inv3;
   double rad,rad2,rad4,rad8,diam,new_coeff2;
   double eoffset;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -151,6 +152,12 @@ void FixWallColloid::wall_particle(int m, int which, double coord)
       ewall[0] -= eoffset;
 
       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");

src/KOKKOS/fix_wall_lj93_kokkos.cpp

@@ -31,6 +31,7 @@ FixWallLJ93Kokkos<DeviceType>::FixWallLJ93Kokkos(LAMMPS *lmp, int narg, char **a
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
   datamask_read = EMPTY_MASK;
   datamask_modify = EMPTY_MASK;
+  virial_flag = 0
 }
 
 /* ----------------------------------------------------------------------

src/MISC/fix_efield.cpp

@@ -55,6 +55,7 @@ FixEfield::FixEfield(LAMMPS *lmp, int narg, char **arg) :
   extscalar = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   qe2f = force->qe2f;
   xstr = ystr = zstr = NULL;
@@ -257,6 +258,11 @@ void FixEfield::post_force(int vflag)
   imageint *image = atom->image;
   int nlocal = atom->nlocal;
 
+  // energy and virial setup
+
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // reallocate efield array if necessary
 
   if (varflag == ATOM && atom->nmax > maxatom) {
@@ -281,6 +287,7 @@ void FixEfield::post_force(int vflag)
 
   double **x = atom->x;
   double fx,fy,fz;
+  double v[6];
 
   // constant efield
 
@@ -306,6 +313,15 @@ void FixEfield::post_force(int vflag)
           fsum[1] += fx;
           fsum[2] += fy;
           fsum[3] += fz;
+          if (evflag) {
+            v[0] = fx*unwrap[0];
+            v[1] = fy*unwrap[1];
+            v[2] = fz*unwrap[2];
+            v[3] = fx*unwrap[1];
+            v[4] = fx*unwrap[2];
+            v[5] = fy*unwrap[2];
+            v_tally(i, v);
+          }
         }
     }
 

src/MOLECULE/fix_cmap.cpp

@@ -76,6 +76,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
   restart_peratom = 1;
   peatom_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   peratom_freq = 1;
   scalar_flag = 1;
   global_freq = 1;

src/RIGID/fix_rigid.cpp

@@ -74,6 +74,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) :
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;

src/RIGID/fix_rigid_small.cpp

@@ -78,6 +78,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) :
   time_integrate = 1;
   rigid_flag = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
   enforce2d_flag = 1;

src/RIGID/fix_shake.cpp

@@ -60,6 +60,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
   MPI_Comm_size(world,&nprocs);
 
   virial_flag = 1;
+  thermo_virial = 1;
   create_attribute = 1;
   dof_flag = 1;
 

src/fix_external.cpp

@@ -37,6 +37,7 @@ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) :
   scalar_flag = 1;
   global_freq = 1;
   virial_flag = 1;
+  thermo_virial = 1;
   extscalar = 1;
 
   if (strcmp(arg[3],"pf/callback") == 0) {

src/fix_wall.cpp

@@ -45,6 +45,7 @@ FixWall::FixWall(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   // parse args
 
@@ -298,7 +299,12 @@ void FixWall::pre_force(int vflag)
 
 void FixWall::post_force(int vflag)
 {
+
+  // energy and virial setup
+
   eflag = 0;
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
   for (int m = 0; m <= nwall; m++) ewall[m] = 0.0;
 
   // coord = current position of wall

src/fix_wall_harmonic.cpp

@@ -34,6 +34,7 @@ FixWallHarmonic::FixWallHarmonic(LAMMPS *lmp, int narg, char **arg) :
 void FixWallHarmonic::wall_particle(int m, int which, double coord)
 {
   double delta,dr,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -60,6 +61,12 @@ void FixWallHarmonic::wall_particle(int m, int which, double coord)
       f[i][dim] -= fwall;
       ewall[0] += epsilon[m]*dr*dr;
       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");

src/fix_wall_lj1043.cpp

@@ -52,6 +52,7 @@ void FixWallLJ1043::precompute(int m)
 void FixWallLJ1043::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r4inv,r10inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -79,5 +80,11 @@ void FixWallLJ1043::wall_particle(int m, int which, double coord)
       ewall[0] += coeff1[m]*r10inv - coeff2[m]*r4inv -
 	coeff3[m]*pow(delta+coeff4[m],-3.0) - offset[m];
       ewall[m+1] += fwall;
+
+      if (evflag) {
+        if (side < 0) vn = -fwall*delta;
+        else vn = fwall*delta;
+        v_tally(dim, i, vn);
+      }
     }
 }

src/fix_wall_lj126.cpp

@@ -48,6 +48,7 @@ void FixWallLJ126::precompute(int m)
 void FixWallLJ126::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r6inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -76,6 +77,12 @@ void FixWallLJ126::wall_particle(int m, int which, double coord)
       f[i][dim] -= fwall;
       ewall[0] += r6inv*(coeff3[m]*r6inv - coeff4[m]) - offset[m];
       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");

src/fix_wall_lj93.cpp

@@ -49,6 +49,7 @@ void FixWallLJ93::precompute(int m)
 void FixWallLJ93::wall_particle(int m, int which, double coord)
 {
   double delta,rinv,r2inv,r4inv,r10inv,fwall;
+  double vn;
 
   double **x = atom->x;
   double **f = atom->f;
@@ -79,6 +80,12 @@ void FixWallLJ93::wall_particle(int m, int which, double coord)
       ewall[0] += coeff3[m]*r4inv*r4inv*rinv -
         coeff4[m]*r2inv*rinv - offset[m];
       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");

src/fix_wall_region.cpp

@@ -25,11 +25,13 @@
 #include "output.h"
 #include "respa.h"
 #include "error.h"
+#include "math_const.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
+using namespace MathConst;
 
-enum{LJ93,LJ126,COLLOID,HARMONIC};
+enum{LJ93,LJ126,LJ1043,COLLOID,HARMONIC};
 
 /* ---------------------------------------------------------------------- */
 
@@ -47,6 +49,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
   extvector = 1;
   respa_level_support = 1;
   ilevel_respa = 0;
+  virial_flag = 1;
 
   // parse args
 
@@ -59,6 +62,7 @@ FixWallRegion::FixWallRegion(LAMMPS *lmp, int narg, char **arg) :
 
   if (strcmp(arg[4],"lj93") == 0) style = LJ93;
   else if (strcmp(arg[4],"lj126") == 0) style = LJ126;
+  else if (strcmp(arg[4],"lj1043") == 0) style = LJ1043;
   else if (strcmp(arg[4],"colloid") == 0) style = COLLOID;
   else if (strcmp(arg[4],"harmonic") == 0) style = HARMONIC;
   else error->all(FLERR,"Illegal fix wall/region command");
@@ -143,6 +147,20 @@ void FixWallRegion::init()
     double r2inv = 1.0/(cutoff*cutoff);
     double r6inv = r2inv*r2inv*r2inv;
     offset = r6inv*(coeff3*r6inv - coeff4);
+  } else if (style == LJ1043) {
+    coeff1 = MY_2PI * 2.0/5.0 * epsilon * pow(sigma,10.0);
+    coeff2 = MY_2PI * epsilon * pow(sigma,4.0);
+    coeff3 = MY_2PI * pow(2.0,1/2.0) / 3 * epsilon * pow(sigma,3.0);
+    coeff4 = 0.61 / pow(2.0,1/2.0) * sigma;
+    coeff5 = coeff1 * 10.0;
+    coeff6 = coeff2 * 4.0;
+    coeff7 = coeff3 * 3.0;
+
+    double rinv = 1.0/cutoff;
+    double r2inv = rinv*rinv;
+    double r4inv = r2inv*r2inv;
+    offset = coeff1*r4inv*r4inv*r2inv - coeff2*r4inv -
+      coeff3*pow(cutoff+coeff4,-3.0);
   } else if (style == COLLOID) {
     coeff1 = -4.0/315.0 * epsilon * pow(sigma,6.0);
     coeff2 = -2.0/3.0 * epsilon;
@@ -186,6 +204,7 @@ void FixWallRegion::post_force(int vflag)
 {
   int i,m,n;
   double rinv,fx,fy,fz,tooclose;
+  double delx, dely, delz, v[6];
 
   double **x = atom->x;
   double **f = atom->f;
@@ -198,13 +217,18 @@ void FixWallRegion::post_force(int vflag)
 
   int onflag = 0;
 
+  // energy and virial setup
+
+  eflag = 0;
+  if (vflag) v_setup(vflag);
+  else evflag = 0;
+
   // region->match() insures particle is in region or on surface, else error
   // if returned contact dist r = 0, is on surface, also an error
   // in COLLOID case, r <= radius is an error
   // initilize ewall after region->prematch(),
   //   so a dynamic region can access last timestep values
 
-  eflag = 0;
   ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0;
 
   for (i = 0; i < nlocal; i++)
@@ -226,19 +250,32 @@ void FixWallRegion::post_force(int vflag)
 
         if (style == LJ93) lj93(region->contact[m].r);
         else if (style == LJ126) lj126(region->contact[m].r);
+        else if (style == LJ1043) lj1043(region->contact[m].r);
         else if (style == COLLOID) colloid(region->contact[m].r,radius[i]);
         else harmonic(region->contact[m].r);
 
-        ewall[0] += eng;
+        delx = region->contact[m].delx;
-        fx = fwall * region->contact[m].delx * rinv;
+        dely = region->contact[m].dely;
-        fy = fwall * region->contact[m].dely * rinv;
+        delz = region->contact[m].delz;
-        fz = fwall * region->contact[m].delz * rinv;
+        fx = fwall * delx * rinv;
+        fy = fwall * dely * rinv;
+        fz = fwall * delz * rinv;
         f[i][0] += fx;
         f[i][1] += fy;
         f[i][2] += fz;
         ewall[1] -= fx;
         ewall[2] -= fy;
         ewall[3] -= fz;
+        ewall[0] += eng;
+        if (evflag) {
+          v[0] = fx*delx;
+          v[1] = fy*dely;
+          v[2] = fz*delz;
+          v[3] = fx*dely;
+          v[4] = fx*delz;
+          v[5] = fy*delz;
+          v_tally(i, v);
+        }
       }
     }
 
@@ -319,6 +356,23 @@ void FixWallRegion::lj126(double r)
   eng = r6inv*(coeff3*r6inv - coeff4) - offset;
 }
 
+/* ----------------------------------------------------------------------
+   LJ 10/4/3 interaction for particle with wall
+   compute eng and fwall = magnitude of wall force
+------------------------------------------------------------------------- */
+
+void FixWallRegion::lj1043(double r)
+{
+  double rinv = 1.0/r;
+  double r2inv = rinv*rinv;
+  double r4inv = r2inv*r2inv;
+  double r10inv = r4inv*r4inv*r2inv;
+  fwall = coeff5*r10inv*rinv - coeff6*r4inv*rinv -
+    coeff7*pow(r+coeff4,-4.0);
+  eng = coeff1*r10inv - coeff2*r4inv -
+    coeff3*pow(r+coeff4,-3.0) - offset;
+}
+
 /* ----------------------------------------------------------------------
    colloid interaction for finite-size particle of rad with wall
    compute eng and fwall = magnitude of wall force

src/fix_wall_region.h

@@ -47,10 +47,12 @@ class FixWallRegion : public Fix {
   char *idregion;
 
   double coeff1,coeff2,coeff3,coeff4,offset;
+  double coeff5,coeff6,coeff7;
   double eng,fwall;
 
   void lj93(double);
   void lj126(double);
+  void lj1043(double);
   void colloid(double, double);
   void harmonic(double);
 };