diff --git a/doc/src/fix_press_langevin.rst b/doc/src/fix_press_langevin.rst
index 8438d72192..02437bd731 100644
--- a/doc/src/fix_press_langevin.rst
+++ b/doc/src/fix_press_langevin.rst
@@ -54,7 +54,7 @@ the Langevin equation such as:
 
    f_P = & \frac{N k_B T_{target}}{V} + \frac{1}{V d}\sum_{i=1}^{N} \vec r_i \cdot \vec f_i - P_{target} \\
    Q\ddot{L} + \alpha{}\dot{L} = & f_P + \beta(t)\\
-   L^{n+1} = & L^{n} + bdt\dot{L}^{n} \frac{bdt^{2}}{2Q} \\
+   L^{n+1} = & L^{n} + bdt\dot{L}^{n} + \frac{bdt^{2}}{2Q} f^{n}_{P} + \frac{bdt}{2Q} \beta^{n+1}  \\
    \dot{L}^{n+1} = & \alpha\dot{L}^{n} + \frac{dt}{2Q}\left(a f^{n}_{P} + f^{n+1}_{P}\right) + \frac{b}{Q}\beta^{n+1} \\
    a = & \frac{1-\frac{\alpha{}dt}{2Q}}{1+\frac{\alpha{}dt}{2Q}} \\
    b = & \frac{1}{1+\frac{\alpha{}dt}{2Q}} \\
diff --git a/doc/src/region.rst b/doc/src/region.rst
index 6d4bbaca4b..dbe14360ff 100644
--- a/doc/src/region.rst
+++ b/doc/src/region.rst
@@ -40,7 +40,7 @@ Syntax
        *plane* args = px py pz nx ny nz
          px,py,pz = point on the plane (distance units)
          nx,ny,nz = direction normal to plane (distance units)
-           px,py,pz can be a variable (see below)
+           px,py,pz,nx,ny,nz can be a variable (see below)
        *prism* args = xlo xhi ylo yhi zlo zhi xy xz yz
          xlo,xhi,ylo,yhi,zlo,zhi = bounds of untilted prism (distance units)
          xy = distance to tilt y in x direction (distance units)
@@ -211,9 +211,11 @@ equal-style :doc:`variable <variable>`.  Likewise, for style *sphere*
 and *ellipsoid* the x-, y-, and z- coordinates of the center of the
 sphere/ellipsoid can be specified as an equal-style variable.  And for
 style *cylinder* the two center positions c1 and c2 for the location of
-the cylinder axes can be specified as a equal-style variable.  For style
-*cone* and *prism* all properties can be defined via equal-style variables.  For
-style *plane* the point can be defined via equal-style variables.
+the cylinder axes can be specified as a equal-style variable.  For styles
+*block*, *cone*, *prism*, and *plane* all properties can be defined via
+equal-style variables.  For style *plane*, the components of the direction
+vector normal to plane should be either all constants or all defined by
+equal-style variables.
 
 If the value is a variable, it should be specified as v_name, where
 name is the variable name.  In this case, the variable will be
@@ -226,6 +228,21 @@ keywords for the simulation box parameters and timestep and elapsed
 time.  Thus it is easy to specify a time-dependent radius or have
 a time dependent position of the sphere or cylinder region.
 
+.. note::
+
+   Whenever a region property, such as a coordinate or an upper/lower
+   bound, is defined via an equal-style variable, the variable should
+   not cause any of the region boundaries to move
+   too far within a single timestep. Otherwise, bad dynamics will occur.
+   "Too far" means a small fraction of the approximate distance of
+   closest approach between two particles, which for the case of Lennard-Jones
+   particles is the distance of the energy minimum while for granular
+   particles it is their diameter. An example is a rapidly varying direction
+   vector in region plane since a small change in the normal to plane will
+   shift the region surface far away from the region point by a large displacement.
+   Similarly, bad dynamics can also occur for fast changing variables employed
+   in the move/rotate options.
+
 See the :doc:`Howto tricilinc <Howto_triclinic>` page for a
 geometric description of triclinic boxes, as defined by LAMMPS, and
 how to transform these parameters to and from other commonly used
diff --git a/src/region_plane.cpp b/src/region_plane.cpp
index 11d96b8266..597f586d8a 100644
--- a/src/region_plane.cpp
+++ b/src/region_plane.cpp
@@ -25,9 +25,11 @@ using namespace LAMMPS_NS;
 /* ---------------------------------------------------------------------- */
 
 RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
-    Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr)
+    Region(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr),
+    nxstr(nullptr), nystr(nullptr), nzstr(nullptr)
 {
   xvar = yvar = zvar = 0.0;
+  nxvar = nyvar = nzvar = 0.0;
 
   options(narg - 8, &arg[8]);
 
@@ -61,15 +63,34 @@ RegPlane::RegPlane(LAMMPS *lmp, int narg, char **arg) :
     zstyle = CONSTANT;
   }
 
+  int nxstyle = (utils::strmatch(arg[5], "^v_")) ? 1 : 0;
+  int nystyle = (utils::strmatch(arg[6], "^v_")) ? 1 : 0;
+  int nzstyle = (utils::strmatch(arg[7], "^v_")) ? 1 : 0;
+
+  if (!nxstyle && !nystyle && !nzstyle) nstyle = CONSTANT;
+  if (nxstyle && nystyle && nzstyle) nstyle = VARIABLE;
+
+  if (nstyle == CONSTANT) {
+    normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
+    normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
+    normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
+  } else if (nstyle == VARIABLE) {
+    normal[0] = 0.0;
+    normal[1] = 0.0;
+    normal[2] = 0.0;
+    nxstr = utils::strdup(arg[5] + 2);
+    nystr = utils::strdup(arg[6] + 2);
+    nzstr = utils::strdup(arg[7] + 2);
+    varshape = 1;
+  } else {
+    error->all(FLERR, "The components of the normal vector should be either all variables or all constants");
+  }
+
   if (varshape) {
     variable_check();
     RegPlane::shape_update();
   }
 
-  normal[0] = xscale * utils::numeric(FLERR, arg[5], false, lmp);
-  normal[1] = yscale * utils::numeric(FLERR, arg[6], false, lmp);
-  normal[2] = zscale * utils::numeric(FLERR, arg[7], false, lmp);
-
   // enforce unit normal
 
   double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
@@ -95,6 +116,9 @@ RegPlane::~RegPlane()
   delete[] xstr;
   delete[] ystr;
   delete[] zstr;
+  delete[] nxstr;
+  delete[] nystr;
+  delete[] nzstr;
   delete[] contact;
 }
 
@@ -174,6 +198,20 @@ void RegPlane::shape_update()
   if (ystyle == VARIABLE) yp = yscale * input->variable->compute_equal(yvar);
 
   if (zstyle == VARIABLE) zp = zscale * input->variable->compute_equal(zvar);
+
+  if (nstyle == VARIABLE) {
+    normal[0] = xscale * input->variable->compute_equal(nxvar);
+    normal[1] = yscale * input->variable->compute_equal(nyvar);
+    normal[2] = zscale * input->variable->compute_equal(nzvar);
+
+    // enforce unit normal
+    double rsq = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
+
+    if (rsq == 0.0) error->all(FLERR, "Illegal region plane normal vector: {} {} {}", normal[0], normal[1], normal[2]);
+    normal[0] /= sqrt(rsq);
+    normal[1] /= sqrt(rsq);
+    normal[2] /= sqrt(rsq);
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -202,4 +240,21 @@ void RegPlane::variable_check()
     if (!input->variable->equalstyle(zvar))
       error->all(FLERR, "Variable {} for region plane is invalid style", zstr);
   }
+
+  if (nstyle == VARIABLE) {
+    nxvar = input->variable->find(nxstr);
+    if (nxvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nxstr);
+    if (!input->variable->equalstyle(nxvar))
+      error->all(FLERR, "Variable {} for region plane is invalid style", nxstr);
+
+    nyvar = input->variable->find(nystr);
+    if (nyvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nystr);
+    if (!input->variable->equalstyle(nyvar))
+      error->all(FLERR, "Variable {} for region plane is invalid style", nystr);
+
+    nzvar = input->variable->find(nzstr);
+    if (nzvar < 0) error->all(FLERR, "Variable {} for region plane does not exist", nzstr);
+    if (!input->variable->equalstyle(nzvar))
+      error->all(FLERR, "Variable {} for region plane is invalid style", nzstr);
+  }
 }
diff --git a/src/region_plane.h b/src/region_plane.h
index 0988fda812..118bd88793 100644
--- a/src/region_plane.h
+++ b/src/region_plane.h
@@ -45,6 +45,10 @@ class RegPlane : public Region {
   int zstyle, zvar;
   char *xstr, *ystr, *zstr;
 
+  int nstyle;
+  int nxvar, nyvar, nzvar;
+  char *nxstr, *nystr, *nzstr;
+
   void variable_check();
 };